# Amplitude Prediction from Uplink to Downlink CSI against Receiver Distortion
in FDD Systems

Chaojin Qinga,∗, Zilong Wanga, Qing Yea, Wenhui Liua, and Linsi Hea

###### Abstract

In frequency division duplex (FDD) massive multiple-input multiple-output
(mMIMO) systems, the reciprocity mismatch caused by receiver distortion
seriously degrades the amplitude prediction performance of channel state
information (CSI). To tackle this issue, from the perspective of distortion
suppression and reciprocity calibration, a lightweight neural network-based
amplitude prediction method is proposed in this paper. Specifically, with the
receiver distortion at the base station (BS), conventional methods are
employed to extract the amplitude feature of uplink CSI. Then, learning along
the direction of the uplink wireless propagation channel, a dedicated and
lightweight distortion-learning network (Dist-LeaNet) is designed to restrain
the receiver distortion and calibrate the amplitude reciprocity between the
uplink and downlink CSI. Subsequently, by cascading, a single hidden layer-
based amplitude-prediction network (Amp-PreNet) is developed to accomplish
amplitude prediction of downlink CSI based on the strong amplitude
reciprocity. Simulation results show that, considering the receiver distortion
in FDD systems, the proposed scheme effectively improves the amplitude
prediction accuracy of downlink CSI while reducing the transmission and
processing delay.

###### keywords:

CSI feedback , massive MIMO , amplitude prediction , receiver distortion ,
lightweight network

††journal: Physical Communication

organization=School of Electrical Engineering and Electronic Information,
Xihua University,city=Chengdu, postcode=610039, country=China

##  1 Introduction

As one of the key techniques in the fifth generation (5G) communications, the
massive multiple-input multiple-output (mMIMO) has shown great prospects in
providing high spectrum and energy efficiency [1, 2, 3]. In frequency division
duplex (FDD) systems, the downlink channel state information (CSI) estimated
by user equipment (UE) usually needs to be fed back to the base station (BS)
[4]. However, due to the large number of antennas, the CSI feedback overhead
in mMIMO systems increases sharply, which results in large transmission and
processing delay, energy consumption, and transmission resource occupation,
etc [5]. Especially, in high-speed scenarios, the transmission and processing
delay may cause the downlink CSI obtained at the BS to be outdated [6].
Therefore, it is crucial to reduce the feedback overhead and processing delay
in FDD mMIMO systems.

Recently, some studies have shown that there is a strong amplitude correlation
(or reciprocity 111Note that, due to the small enough difference in frequency-
independent parameters between the uplink CSI and downlink CSI in FDD systems,
we assume the reciprocity in this paper.) between the uplink CSI and downlink
CSI in FDD systems [7, 8]. Hence, it becomes popular to use deep learning (DL)
to directly predict the downlink CSI from the uplink CSI to reduce/eliminate
the feedback process [9, 10, 11]. In [9], considering the position-to-channel
mapping is bijective, a sparse complex-valued neural network (SCNet) is
proposed to approximate the uplink-to-downlink mapping function. In [10],
according to the spatial correlation in time-varying scenarios, a
convolutional neural network (CNN)-based downlink channel prediction method is
investigated. Under the premise of considering the channel time invariance, an
attention-based deep learning network is proposed in [11] to directly predict
the downlink CSI from the uplink CSI.

In [7, 8, 9, 10, 11], the feedback and prediction methods are mainly based on
the reciprocity between the uplink and downlink CSIs. However, this
reciprocity is vulnerable in practical systems. This is because the overall
channel consists of not only the wireless propagation channel, but also the
radio frequency (RF) front-end, e.g., analog-to-digital converters (ADCs),
filters and low-noise amplifiers (LNAs), etc [12, 13]. Although uplink and
downlink wireless propagation channels may be reciprocal, the hardware
imperfection (HI) of RF front-end inevitably introduces nonlinear distortion
into the amplitude and phase of the transmitted and received signals [14].
This nonlinear distortion causes the amplitude mismatch and phase mismatch,
thereby resulting in the reciprocity mismatch between the uplink and downlink
CSIs.

Due to the significant impact of reciprocity mismatch on system performance,
reciprocity calibration is crucial for communication systems [12]. In both
regular MIMO and mMIMO systems, the conventional reciprocity calibration
method is based on the dedicated hardware circuits [15], which increases the
energy consumption and hardware cost that comes from RF chains required to
support a number of antennas [16]. Besides, the existing reciprocity
calibration is usually investigated in time division duplex (TDD) systems. In
contrast, there is limited literature addressing the issue of reciprocity
calibration in FDD systems when utilizing the reciprocity. That is, in
practical FDD mMIMO systems, the reciprocity between uplink and downlink
channels is destroyed by the nonlinear distortion due to the difference
between uplink and downlink hardware, making the reciprocity-based prediction
results inaccurate, and even resulting in the prediction process impossible to
achieve.

Therefore, it is vital to consider the nonlinear distortion before utilizing
reciprocity and the reciprocity calibration is also essential for the
reciprocity-based amplitude prediction in FDD systems. To suppress the impact
of distortion on reciprocity and calibrate the reciprocity, this paper
proposes an amplitude prediction scheme against receiver distortion. To the
best of our knowledge, the amplitude reciprocity-based CSI prediction by
considering distortion has not been investigated in FDD systems. The main
contributions of this paper are summarized as follows:

  * 1.

We propose a more practical scenario which considers the distortion before
utilizing the reciprocity in FDD systems. Specifically, we take the receiver
distortion at BS as an example to illustrate that reciprocity will be affected
by distortion, which is a valuable reference for both UE and BS.

  * 2.

We design a network architecture cascading distortion learning and amplitude
prediction to improve the practicality and accuracy of amplitude prediction.
Specifically, learning along the direction of the uplink wireless propagation
channel, a distortion-learning network (Dist-LeaNet) is designed to restrain
the receiver distortion and calibrate the amplitude reciprocity between the
uplink and downlink CSI. Subsequently, based on the channel amplitude
reciprocity, an amplitude-prediction network (Amp-PreNet) is developed to
predict the amplitude of downlink CSI directly at the BS, thus avoiding the
overhead and transmission delay caused by feedback.

  * 3.

We construct a lightweight learning and prediction network architecture to
reduce the processing delay and computational complexity for the BS receiver.
Due to that the nonlinear distortion varies slowly compared with the wireless
propagation channel, the features of distortion are easy to capture. Hence,
Dist-LeaNet is constructed with a lightweight network architecture. With the
assistance of Dist-LeaNet, Amp-PreNet is also constructed with lightweight
network architecture based on the strong amplitude correlation.

The rest of this paper is organized as follows. In Section 2, we introduce the
system model. Then, the amplitude prediction scheme against receiver
distortion is presented in Section 3 and followed by numerical results in
Section 4. Finally, Section 5 concludes our work.

Notation: Boldface upper case and lower case letters denote matrix and vector
respectively. ${\cal N}\left({\mu,{\sigma^{2}}}\right)$ stands for normal
distribution with mean $\mu$ and variance ${{\sigma^{2}}}$; ${\cal
U}\left({a,b}\right)$ stands for uniform distribution on the interval
$\left({a,b}\right)$; $|\cdot|$ denotes the operation of taking the modulus of
a complex value; ${\left(\cdot\right)^{T}}$ denotes transpose; $E[\cdot]$
represents the expectation operation; ${\left\|\cdot\right\|}$ is the
Euclidean norm.

##  2 System Model

Figure 1: System model.

The system model is given in Fig. 1, in which an FDD massive MIMO system that
consists of a BS with $N$ antennas and $U$ single-antenna users in speed $v$
is considered. At the BS, the received uplink signal from user-$u$, denoted as
${\widetilde{\bf{Y}}_{u}}\in\mathbb{C}^{N\times N}$, is given by

${\widetilde{\bf{Y}}_{u}}={{\bf{g}}_{u}}{{\bf{x}}_{u}}+{{\bf{N}}_{u}},$ (1)

where ${{\bf{g}}_{u}}\in{{\mathbb{C}}^{N\times 1}}$ denotes the uplink channel
(i.e., uplink CSI) from the user-$u$ to the BS in the angular domain,
${{\bf{x}}_{u}}\in\mathbb{C}^{1\times N}$ stands $N$-length uplink pilot and
data of user-$u$, and ${{\bf{N}}_{u}}\in{{\mathbb{C}}^{N\times N}}$ is the
circularly symmetric complex Gaussian (CSCG) noise with zero-mean and variance
${\sigma_{u}^{2}}$. With the antenna diversity [17], the uplink pilot and data
of each UE is received by $N$ BS antennas to form the $N\times N$ complex
signal ${\widetilde{\bf{Y}}_{u}}$. From [18], there exists a correlation
between the uplink and downlink channels due to the shared common physical
paths and similar spatial propagation characteristics. For example, the
downlink CSI is constructed by utilizing frequency-independent parameters
between the uplink and downlink channels in the angular domain [19].
Therefore, from [20, 21], the downlink CSI of user-$u$ (i.e.,
${\bf{h}}_{u}\in\mathbb{C}^{N\times 1}$) can be recovered from the uplink CSI
${{\bf{g}}_{u}}$. However, the inevitable uplink distortion (e.g., caused by
imperfect hardware) makes this processing difficult. By denoting the mapping
function of equivalent distortion at the BS as
${f_{{\text{R-dis}}}}\left(\cdot\right)$, the distorted signal of
${\widetilde{\bf{Y}}_{u}}$, denoted by ${{\bf{Y}}_{u}}\in\mathbb{C}^{N\times
N}$, is expressed as

${{\bf{Y}}_{u}}\triangleq{f_{{\text{R-dis}}}}\left({{\widetilde{\bf{Y}}_{u}}}\right).$
(2)

Then, the uplink channel ${\bf{g}}_{u}$ is estimated according to
${{\bf{Y}}_{u}}$ and the uplink pilot in ${{\bf{x}}_{u}}$. By denoting the
estimated ${\bf{g}}_{u}$ as $\widetilde{\bf{g}}_{u}$
($\widetilde{\bf{g}}_{u}\in\mathbb{C}^{N\times 1}$), our work aims to utilize
the strong amplitude correlation to predict the amplitude of downlink CSI
$\bf{h}_{u}$ from $\widetilde{\bf{g}}_{u}$ [7, 8] in this paper. First, the
amplitude feature of $\widetilde{\bf{g}}_{u}$ is extracted. Subsequently, we
build two dedicated networks, Dist-LeaNet and Amp-PreNet, to restrain the
distortion of receiver and enhance the prediction accuracy of downlink CSI
amplitude, respectively. The details are described in Section 3. However, the
phase information exhibits unique importance due to its frequency-dependent
nature, which is directly fed back to the BS [8].

##  3 Amplitude Prediction Scheme against Receiver Distortion

To effectively utilize the reciprocity in practical scenarios, we present the
proposed amplitude prediction scheme against receiver distortion in this
section. In Section 3.1, the model of uplink receiver distortion is presented.
With the uplink receiver distortion at the BS, we develop Dist-LeaNet and Amp-
PreNet to restrain the distortion and predict the downlink CSI amplitude,
respectively. Both Dist-LeaNet and Amp-PreNet are elaborated in Section 3.2.

###  3.1 Uplink Receiver Distortion

In the uplink communication, nonlinear distortion is inevitably encountered
[22], e.g., the distortion of user’s power amplifiers (PAs), the distortion of
BS’s LNAs and ADCs, etc [23]. We mainly take the receiver distortion of BS as
an example to represent the uplink distortion, which has reference value for
both UE and BS.

Specifically, the nonlinear distortion varies slowly compared with the
wireless propagation channel [24]. Therefore, to further represent the
distortion function ${f_{{\text{R-dis}}}}\left(\cdot\right)$, the receiver
distortion at the BS is denoted as
${{\bf{D}}_{{\text{R-BS}}}}\in{{\mathbb{C}}^{N\times N}}$, wherein its
diagonal elements represent the amplitude and phase distortion of each
hardware at different antennas, and the off-diagonal elements correspond to
the crosstalk and mutual coupling effect between different antennas [23]. The
proper hardware circuit design can ensure the nearly-zero crosstalk, and the
antenna mutual coupling effect is often ignored [25]. Therefore, the receiver
distortion matrix can be regarded to be diagonal, which is expressed as [26]

${{\bf{D}}_{{\text{R-BS}}}}={\text{diag}}\left({{r_{1{\text{,BS}}}},\cdots,{r_{n{\text{,BS}}}},\cdots,{r_{N{\text{,BS}}}}}\right),$
(3)

where
${r_{n,{\text{BS}}}}=\left|{{r_{n,{\text{BS}}}}}\right|{e^{j\phi_{n,{\text{BS}}}^{\text{r}}}}$
$\left({n=1,2,\cdots,N}\right)$. According to [26], the amplitudes of the
distortion obey log-normal distribution, and the phases of the distortion obey
uniform distribution, i.e.,

${\rm{ln}}\left|{{r_{n,{\rm{BS}}}}}\right|\sim{\cal
N}\left({0,\delta_{{\rm{r,BS}}}^{2}}\right),\phi_{n,{\rm{BS}}}^{\rm{r}}\sim{\cal
U}\left[{-{\theta_{{\rm{r,BS}}}},{\theta_{{\rm{r,BS}}}}}\right].$

###  3.2 Dist-LeaNet and Amp-PreNet

In order to restrain the receiver distortion and obtain amplitude feature of
the uplink wireless propagation channel, we construct the lightweight and
effective Dist-LeaNet, which is supposed to be considered when channel
reciprocity is involved. Then, a recovered uplink CSI amplitude feature
${\widehat{\bf{g}}_{u,{\text{amp}}}}$, is learned from Dist-LeaNet.
Subsequently, to predict the downlink CSI amplitude feature, we design the
lightweight Amp-PreNet, which utilizes the amplitude correlation of CSI in the
angular domain [8]. The corresponding network design, training and deployment
are as follows.

####  3.2.1 Network Design

According to [27], choosing the appropriate number of layers and hidden
neurons is still a challenge in the neural network (NN). That is, for a
specific network design, there is currently no established theoretical
guidance on the optimal number of layers or the number of neurons to be
included at each layer. Typically, complex hyper-parameter tuning is
necessary. Based on plenty of experimental results, we design the lightweight
Dist-LeaNet and Amp-PreNet, both of which are single hidden-layer NN.
Specifically, considering the trade-off between performance and complexity, we
train the network with different number of layers and neurons. After verifying
the performance of the trained network, we select a suitable lightweight
network architecture to reduce the computational complexity while improve the
prediction performance compared with [7]. The network architectures of Dist-
LeaNet and Amp-PreNet are summarized in Table 1, and the detailed descriptions
are given as follows.

Table 1: Architecture of Dist-LeaNet and Amp-PreNet.

Layer | Input | Hidden | Output
---|---|---|---
Dist-LeaNet | Amp-PreNet | Dist-LeaNet | Amp-PreNet | Dist-LeaNet | Amp-PreNet
Batch normalization | $\surd$ | $\times$ | $\times$ | $\times$ | $\times$ | $\times$
Neuron number | $N$ | $N$ | $2N$ | $2N$ | $N$ | $N$
Activation function | None | None | Linear | LReLU | Linear | Linear

In both Dist-LeaNet and Amp-PreNet, the neurons of the input layer, hidden
layer, and output layer are $N$, $2N$, and $N$, respectively. In Dist-LeaNet,
a batch normalization (BN) is employed for the input layer, which normalizes
the network input as zero mean and unit variance. For the hidden layer and
output layer of Dist-LeaNet, the linear activation is employed. Then, the
Dist-LeaNet is followed by Amp-PreNet with the cascaded mode, i.e., the output
of Dist-LeaNet is the input of Amp-PreNet. Without BN, the leaky rectified
unit (LReLU) [28] and linear activation are adopted for the hidden layer and
output layer of Amp-PreNet, respectively.

With the estimated uplink CSI $\widetilde{\bf{g}}_{u}$, we extract its
amplitude feature (denoted as ${\widetilde{\bf{g}}_{u,{\text{amp}}}}$)
according to

${\widetilde{\bf{g}}_{u,{\text{amp}}}}={\left[{\left|{{{\widetilde{g}}_{u,1}}}\right|,\left|{{{\widetilde{g}}_{u,2}}}\right|,\ldots,\left|{{{\widetilde{g}}_{u,N}}}\right|}\right]^{T}}.$
(4)

Due to the BS’s nonlinear distortion (e.g., the LNA and ADC in BS hardware),
${\widetilde{\bf{g}}_{u,{\text{amp}}}}$ cannot use to map the amplitude of
downlink CSI of user-$u$ (i.e., the amplitude of $\textbf{h}_{u}$). This
results in that the methods of CSI prediction and recovery in [7, 8, 9, 10,
11, 20], cannot be applied directly. Thus, we develop Dist-LeaNet to learn
along the direction of the amplitude of $\textbf{g}_{u}$ and restrain the
nonlinear distortion at the BS, which is expressed as

${\widehat{\bf{g}}_{u,{\text{amp}}}}={f_{{\text{Dist-
Lea}}}}\left({{{\widetilde{\bf{g}}}_{u,{\text{amp}}}},{\bf{\Theta}_{{\text{Dist-
Lea}}}}}\right),$ (5)

where ${f_{{\text{Dist-Lea}}}}\left(\cdot\right)$ and
${{\bf{\Theta}_{{\text{Dist-Lea}}}}}$ denote the mapping function of
distortion suppression and the training parameters of Dist-LeaNet,
respectively.

On the basis of the obtained ${\widehat{\bf{g}}_{u,{\text{amp}}}}$, the
amplitude of downlink CSI of user-$u$ (i.e., the amplitude of
$\textbf{h}_{u}$) can be mapped. Thus, based on the strong amplitude
correlation, we construct Amp-PreNet to predict the amplitude feature of
downlink CSI ${\widehat{\bf{h}}_{u,{\text{amp}}}}$, which can be expressed as

${\widehat{\bf{h}}_{u,{\text{amp}}}}={f_{{\text{Amp-
Pre}}}}\left({{{\widehat{\bf{g}}}_{u,{\text{amp}}}},{\bf{\Theta}_{{\text{Amp-
Pre}}}}}\right),$ (6)

where ${f_{{\text{Amp-Pre}}}}\left(\cdot\right)$ and
${{\bf{\Theta}_{{\text{Amp-Pre}}}}}$ denote the mapping function of amplitude
prediction and the training parameters of Amp-PreNet, respectively.

####  3.2.2 Training and Deployment

The training sets are acquired by simulation, and a significant amount of data
samples are collected to train Dist-LeaNet and Amp-PreNet. Specifically, these
data samples are generated as follows.

Amplitude correlated channels are generated by MATLAB 5G Toolbox, which is
subject to specifications of the Clustered-Delay-Line (CDL) channel model in
3GPP TR 38.901 [29]. Similar to the setting in [18], the frequency-independent
parameters (e.g., the angle of departure (AoD)) are fixed, while varying the
complex gain of each path between the uplink and downlink channels. The AoD of
downlink CSI is approximately the same as the angle of arrival (AoA) of the
uplink CSI in a short time slot [18], showing a relatively strong amplitude
correlation in the angle domain. The amplitude attenuation of clusters also
reflects the amplitude reciprocity [8], due to the similar geographical
environment in a short time slot. Thus, ${{\bf{g}}_{u}}$ and ${{\bf{h}}_{u}}$
are obtained by transforming the generated uplink and downlink channels to the
angular domain, respectively [18].

To train Dist-LeaNet and Amp-PreNet, we use the amplitude of ${{\bf{g}}_{u}}$
and ${{\bf{h}}_{u}}$ as network labels, respectively, with the joint training
method. The optimization goal of Dist-LeaNet is to minimize the mean squared
error (MSE) between ${\widehat{\bf{g}}_{u,{\text{amp}}}}$ and
${{\bf{g}}_{u,{\text{amp}}}}$, which is derived as

$\mathop{\min}\limits_{{\bf{\Theta}_{{\text{Dist-
Lea}}}}}E\left[{{{\left\|{{f_{{\text{Dist-
Lea}}}}\left({{{\widetilde{\bf{g}}}_{u,{\text{amp}}}},{\bf{\Theta}_{{\text{Dist-
Lea}}}}}\right)-{{\bf{g}}_{u,{\text{amp}}}}}\right\|}^{2}}}\right].$ (7)

Similarly, the Amp-PreNet minimizes the MSE of the downlink CSI amplitude,
i.e.,
$E\left[{{{\left\|{{{\widehat{\bf{h}}}_{u,{\text{amp}}}}-{{\bf{h}}_{u,{\text{amp}}}}}\right\|}^{2}}}\right]$,
which is further expressed by

$\mathop{\min}\limits_{{\bf{\Theta}_{{\text{Amp-
Pre}}}}}E\left[{{{\left\|{{f_{{\text{Amp-
Pre}}}}\left({{{\widehat{\bf{g}}}_{u,{\text{amp}}}},{\bf{\Theta}_{{\text{Amp-
Pre}}}}}\right)-{{\bf{h}}_{u,{\text{amp}}}}}\right\|}^{2}}}\right].$ (8)

We perform the joint training once for both Dist-LeaNet and Amp-PreNet, and
save the trained network parameters for testing.

By using the Dist-LeaNet, the high precision uplink CSI amplitude
${\widehat{\bf{g}}_{u,{\text{amp}}}}$ is obtained. Then,
${\widehat{\bf{g}}_{u,{\text{amp}}}}$ is used to predict the downlink CSI
amplitude ${\widehat{\bf{h}}_{u,{\text{amp}}}}$ in the Amp-PreNet. We consider
the distortion before utilizing the channel reciprocity to accomplish
amplitude prediction of downlink CSI. The proposed scheme demonstrates a
better prediction accuracy and reduces the impact of time delay effectively in
a practical scenario.

##  4 Experiment results

In this section, we provide numerical results of the proposed scheme.
Definitions and basic parameters involved in simulation are first given in
Section 4.1. Subsequently, to verify the effectiveness of the proposed scheme,
the normalized mean squared error (NMSE) of the predicted downlink CSI
amplitude is given in Section 4.2. Finally, computational complexity and
online running time comparison analysis are shown in Section 4.3.

###  4.1 Parameters Setting

Definitions involved in simulations are given as follows. The equivalent
signal-to-noise ratio (SNR) and NMSE are defined as similar to [7]. During the
experiments, $v=300$ km/h, $\delta_{{\text{r,BS}}}^{2}=1$, and
${\theta_{{\text{r,BS}}}}=\pi$ are considered, respectively. The probability
density functions (PDF) of the amplitude and phase of receiver distortion are
shown in Fig. 2. Following the setting in [8], we set the uplink frequency to
5.1 GHz and the downlink frequency to 5.3 GHz. Thus, according to
${f_{m}}={{v{f_{c}}}\mathord{\left/{\vphantom{{v{f_{c}}}c}}\right.\kern-1.2pt}c}$
[17] with ${f_{m}}$, ${f_{c}}$, and ${c}$ being the maximum doppler shift, the
carrier frequency, and the speed of light, respectively, the maximum doppler
shift for the uplink CSI and downlink CSI are 1418 Hz and 1473 Hz,
respectively. The complex-valued Zadoff-Chu (ZC) sequence [30] is employed as
the pilot for uplink channel estimation with least squares (LS) criterion in
the simulation. For Dist-LeaNet and Amp-PreNet, their training and testing
data-sets are generated according to (4). The sample numbers of training set,
validation set, and testing set are 30,000, 5,000, and 15,000, respectively.
In this paper, the NMSE performance of the proposed scheme is compared with
those of [7] and [8]. In addition, to verify the effectiveness of Dist-LeaNet,
the proposed scheme without Dist-LeaNet, denoted as “Proposed (without Dist-
LeaNet)”, is also simulated. It is worth noting that inspired by signal
detection [31, 32], the detection-based amplitude prediction is an interesting
topic. However, this is beyond the scope of this paper and prompts us to
conduct exploratory research in the future.

Figure 2: (a) The PDF of the amplitude distortion at the receiver, where
$\delta_{{\text{r,BS}}}^{2}=1$. (b) The PDF of the phase distortion at the
receiver, where ${\theta_{{\text{r,BS}}}}=\pi$.

###  4.2 NMSE Performance

To validate the effectiveness of amplitude prediction, NMSE curves of the
recovered downlink CSI amplitude are plotted in Fig. 3, where $N=128$ is
considered. From Fig. 3, it can be observed that the NMSE of “Proposed” is
smaller than those of “Ref [7]”, “Ref [8]”, and “Proposed (without Dist-
LeaNet)”, showing the effectiveness of the proposed scheme in recovering the
downlink CSI amplitude. Specifically, the NMSE of “Proposed” is smaller than
that of “Proposed (without Dist-LeaNet)”, which confirms that Dist-LeaNet
plays an essential role for the proposed scheme in distortion suppression and
reciprocity calibration. In addition, for each given SNR, the NMSE of
“Proposed (without Dist-LeaNet)” is lower than those of “Ref [7]” and “Ref
[8]”, which indicates the effectiveness of Amp-PreNet in CSI prediction.
Overall, the proposed scheme proves to be advantageous in improving the NMSE
performance in various SNR scenarios.

Figure 3: NMSE of downlink CSI amplitude versus SNR, where $N=128$.

To verify the NMSE performance against the impact of $N$, NMSE curves of “Ref
[7]”, “Ref [8]”, “Proposed (without Dist-LeaNet)”, and “Proposed” are plotted
in Fig. 4, where $N=64$, $N=128$, and $N=256$ are considered. For each given
$N$, the NMSE of downlink CSI amplitude of “Proposed” is smaller than those of
“Ref [7]”, “Ref [8]”, and “Proposed (without Dist-LeaNet)”. As the increase of
$N$ (i.e., the number of antennas increases), the NMSE increases due to the
more nonlinear distortion introduced on antennas at the same SNR. Thus,
compared with “Ref [7]”, “Ref [8]”, and “Proposed (without Dist-LeaNet)”, the
proposed Dist-LeaNet restrains the distortion and Amp-PreNet predicts the CSI
amplitude effectively against varying $N$.

Figure 4: NMSE of downlink CSI amplitude versus SNR under different $N$.

###  4.3 Computational Complexity and Online Running Time

In this subsection, the computational complexity and online running time of
“Ref [7]”, “Ref [8]”, and “Proposed” are presented and analyzed as follows.

####  4.3.1 Computational Complexity Analysis

The number of floating-point operations (FLOPs) is considered as the metric of
computational complexity, which can be used to describe the NN complexity
[18]. According to [18], the FLOPs of “Ref [7]”, “Ref [8]” and “Proposed” are
$20{N^{2}}-6N$,
${{{N^{2}}}\mathord{\left/{\vphantom{{{N^{2}}}2}}\right.\kern-1.2pt}2}+{{26841N}\mathord{\left/{\vphantom{{26841N}8}}\right.\kern-1.2pt}8}$,
and $16{N^{2}}-6N$, respectively. The comparison and case details of
computational complexity are given in Table 2 and Fig. 5 (a). For $N<217$, the
proposed scheme demonstrates the lowest computational complexity. When
$N\geqslant 217$, the FLOPs number of “Ref [8]” is lower than those of
“Proposed” and “Ref [7]”. Nevertheless, the proposed scheme improves the
prediction performance of downlink CSI amplitude greatly at the expense of a
tolerable computational complexity.

Figure 5: (a) The FLOPs number versus $N$. (b) Online running time comparison
of “Ref [7]”, “Ref [8]” and “Proposed” for $10^{5}$ experiments. Table 2:
Analysis of Computational Complexity.

Method | Complexity | Case1 $\left({N=64}\right)$ | Case2 $\left({N=128}\right)$ | Case3 $\left({N=256}\right)$
---|---|---|---|---
Ref [7] | $20{N^{2}}-6N$ | 81,536 | 326,912 | 1,309,180
Ref [8] | ${{{N^{2}}}\mathord{\left/{\vphantom{{{N^{2}}}2}}\right.\kern-1.2pt}2}+{{26841N}\mathord{\left/{\vphantom{{26841N}8}}\right.\kern-1.2pt}8}$ | 216,776 | 437,648 | 891,680
Proposed | $16{N^{2}}-6N$ | 65,152 | 261,376 | 1,047,040

####  4.3.2 Online Running Time

The comparison of online running time is given in Fig. 5 (b), where $N=64$,
$N=128$, and $N=256$ are considered. For a fair comparison, $10^{5}$
experiments of online running are conducted for “Ref [7]”, “Ref [8]” and
“Proposed” on the same computer. From Fig. 5 (b), for each given $N$, the
online running time of “Proposed” is shorter than those of “Ref [7]” and “Ref
[8]”. This reflects that the proposed scheme reduces the transmission and
processing delay effectively, due to the application of prediction method and
lightweight network architecture. Additionally, to demonstrate that the
proposed scheme can prevent the downlink CSI outdated, the online running time
is compared with the coherence time of downlink CSI. Considering that the
correlation coefficient of the channel at any two time points within the
coherent time is not less than 0.5 [33], the maximum Doppler frequency shift
${f_{m}}$ is used to measure the coherence time of the channel. According to
$M={9\mathord{\left/{\vphantom{9{\left({16\pi{f_{m}}}\right)}}}\right.\kern-1.2pt}{\left({16\pi{f_{m}}}\right)}}$
[17], the coherence time $M$ of the downlink CSI is 0.122 ms. However, when
$N=64$, $N=128$, and $N=256$ are considered, the online running times for each
experiment of the proposed scheme are 0.0069 ms, 0.0099 ms, and 0.0196 ms,
respectively, which are smaller than the coherence time of downlink CSI. This
indicates that although there is a certain delay, the proposed scheme can
effectively prevent the predicted downlink CSI amplitude from becoming
outdated.

##  5 Conclusion

This paper presents an amplitude prediction scheme from uplink to downlink CSI
against receiver distortion in FDD systems. By using a lightweight and
dedicated Dist-LeaNet, the amplitude feature of the uplink wireless
propagation channel is obtained after distortion suppression and reciprocity
calibration. Then, with the uplink CSI, the downlink CSI amplitude is
predicted by a lightweight Amp-PreNet. Experiments show that, compared with
methods that don’t consider the distortion in communication systems, the
proposed scheme is more practical and achieves a better prediction accuracy on
NMSE performance of the downlink CSI amplitude. This idea of considering and
handling distortion has reference significance for both UE and BS. In our
future work, we will conduct exploratory research on detection-based amplitude
prediction methods.

## Acknowledgements

This work is supported in part by the Sichuan Science and Technology Program
(Grant No. 2023YFG0316), the Industry-University Research Innovation Fund of
China University (Grant No. 2021ITA10016), the Key Scientific Research Fund of
Xihua University (Grant No. Z1320929), and the Special Funds of Industry
Development of Sichuan Province (Grant No. zyf-2018-056).

## References

  * Fu et al. [2019] S. Fu, Z. Su, Y. Jia, H. Zhou, Y. Jin, J. Ren, B. Wu, K. M. S. Huq,  Interference cooperation via distributed game in 5g networks,  IEEE Internet of Things Journal 6 (2019) 311–320. doi:10.1109/JIOT.2017.2743116.
  * Da Silva et al. [2023] H. T. P. Da Silva, H. S. Silva, M. S. Alencar, W. J. L. De Queiroz, U. S. Dias,  Energy and spectral efficiencies of cell-free millimeter-wave massive mimo systems under rain attenuation based on ray tracing simulations,  IEEE Access 11 (2023) 26979–26995. doi:10.1109/ACCESS.2023.3257847.
  * Vasudevan et al. [2019] K. Vasudevan, S. Singh, A. P. K. Reddy,  Coherent receiver for turbo coded single-user massive mimo-ofdm with retransmissions,  in: S. Mohammady (Ed.), Multiplexing, IntechOpen, Rijeka, 2019. doi:10.5772/intechopen.85893.
  * Qing et al. [2020] C. Qing, B. Cai, Q. Yang, J. Wang, C. Huang,  Elm-based superimposed csi feedback for fdd massive mimo system,  IEEE Access 8 (2020) 53408–53418. doi:10.1109/ACCESS.2020.2980969.
  * Zhang et al. [2023] B. Zhang, H. Li, X. Liang, X. Gu, L. Zhang,  Model transmission-based online updating approach for massive mimo csi feedback,  IEEE Communications Letters (2023) 1–1. doi:10.1109/LCOMM.2023.3265680.
  * Wang et al. [2022] X. Wang, Y. Shi, W. Xin, G. Yang, Z. Jiang,  A polynomial fourier transform based channel prediction scheme for high-mobility massive mimo scenarios,  in: 2022 IEEE/CIC International Conference on Communications in China (ICCC), 2022, pp. 488–493. doi:10.1109/ICCC55456.2022.9880832.
  * Qing et al. [2022] C. Qing, Q. Ye, W. Liu, J. Wang,  Fusion learning for 1-bit cs-based superimposed csi feedback with bi-directional channel reciprocity,  IEEE Communications Letters 26 (2022) 813–817. doi:10.1109/LCOMM.2022.3145099.
  * Liu et al. [2019] Z. Liu, L. Zhang, Z. Ding,  Exploiting bi-directional channel reciprocity in deep learning for low rate massive mimo csi feedback,  IEEE Wireless Communications Letters 8 (2019) 889–892. doi:10.1109/LWC.2019.2898662.
  * Yang et al. [2019] Y. Yang, F. Gao, G. Y. Li, M. Jian,  Deep learning-based downlink channel prediction for fdd massive mimo system,  IEEE Communications Letters 23 (2019) 1994–1998. doi:10.1109/LCOMM.2019.2934851.
  * Zhang et al. [2022] S. Zhang, Y. Zhang, Z. Zhang, J. Zhang, L. Xia, T. Jiang,  Deep learning-based time-varying channel prediction for mimo systems,  in: 2022 IEEE 95th Vehicular Technology Conference: (VTC2022-Spring), 2022, pp. 1–5. doi:10.1109/VTC2022-Spring54318.2022.9860388.
  * Jiang et al. [2023] Y. Jiang, W. Lin, W. Zhao, C. Wang,  Acsinet: Attention-based deep learning network for csi prediction in fdd mimo systems,  IEEE Wireless Communications Letters 12 (2023) 471–475. doi:10.1109/LWC.2022.3230950.
  * Nie et al. [2020] R. Nie, L. Chen, N. Zhao, Y. Chen, F. R. Yu, G. Wei,  Relaying systems with reciprocity mismatch: Impact analysis and calibration,  IEEE Transactions on Communications 68 (2020) 4035–4049. doi:10.1109/TCOMM.2020.2982632.
  * Rogalin et al. [2014] R. Rogalin, O. Y. Bursalioglu, H. Papadopoulos, G. Caire, A. F. Molisch, A. Michaloliakos, V. Balan, K. Psounis,  Scalable synchronization and reciprocity calibration for distributed multiuser mimo,  IEEE Transactions on Wireless Communications 13 (2014) 1815–1831. doi:10.1109/TWC.2014.030314.130474.
  * Aghdam and Eriksson [2019] S. R. Aghdam, T. Eriksson,  On the performance of distortion-aware linear receivers in uplink massive mimo systems,  in: 2019 16th International Symposium on Wireless Communication Systems (ISWCS), 2019, pp. 208–212. doi:10.1109/ISWCS.2019.8877331.
  * Song et al. [2022] G. Song, Z. Wang, Y. Du, L. Xiao, T. Jiang,  Enabling affordable implicit channel feedback for internet of things,  IEEE Internet of Things Journal 9 (2022) 16606–16618. doi:10.1109/JIOT.2022.3151344.
  * Janjua et al. [2023] M. B. Janjua, H. T. Abbas, K. A. Qaraqe, H. Arslan,  Beam selection for ambient backscatter communication in beamspace mmwave symbiotic radio,  IEEE Wireless Communications Letters 12 (2023) 560–564. doi:10.1109/LWC.2023.3235104.
  * Rappaport [1996] T. S. Rappaport, Wireless communications: Principles and practice, prentice hall PTR New Jersey, 1996\.
  * Guo et al. [2022] J. Guo, C.-K. Wen, S. Jin,  Canet: Uplink-aided downlink channel acquisition in fdd massive mimo using deep learning,  IEEE Transactions on Communications 70 (2022) 199–214. doi:10.1109/TCOMM.2021.3120294.
  * Han et al. [2019] Y. Han, T.-H. Hsu, C.-K. Wen, K.-K. Wong, S. Jin,  Efficient downlink channel reconstruction for fdd multi-antenna systems,  IEEE Transactions on Wireless Communications 18 (2019) 3161–3176. doi:10.1109/TWC.2019.2911497.
  * Arnold et al. [2019] M. Arnold, S. Dörner, S. Cammerer, S. Yan, J. Hoydis, S. t. Brink,  Enabling fdd massive mimo through deep learning-based channel prediction,  arXiv preprint arXiv:1901.03664 (2019).
  * Hou et al. [2021] W. Hou, J. Sun, G. Gui, T. Ohtsuki, A. M. Elbir, H. Gacanin, H. Sari,  Federated learning for dl-csi prediction in fdd massive mimo systems,  IEEE Wireless Communications Letters 10 (2021) 1810–1814. doi:10.1109/LWC.2021.3081695.
  * Qing et al. [2022] C. Qing, L. Dong, L. Wang, J. Wang, C. Huang,  Joint model and data-driven receiver design for data-dependent superimposed training scheme with imperfect hardware,  IEEE Transactions on Wireless Communications 21 (2022) 3779–3791. doi:10.1109/TWC.2021.3123948.
  * Wei et al. [2022] X. Wei, Y. Jiang, X. Wang, C. Shen,  Tx-rx reciprocity calibration for hybrid massive mimo systems,  IEEE Wireless Communications Letters 11 (2022) 431–435. doi:10.1109/LWC.2021.3132066.
  * Morte Palacios et al. [2020] J. Morte Palacios, O. Raeesi, A. Gokceoglu, M. Valkama,  Impact of channel non-reciprocity in cell-free massive mimo,  IEEE Wireless Communications Letters 9 (2020) 344–348. doi:10.1109/LWC.2019.2954513.
  * Jiang et al. [2015] X. Jiang, M. Čirkić, F. Kaltenberger, E. G. Larsson, L. Deneire, R. Knopp,  Mimo-tdd reciprocity under hardware imbalances: Experimental results,  in: 2015 IEEE International Conference on Communications (ICC), 2015, pp. 4949–4953. doi:10.1109/ICC.2015.7249107.
  * Wei et al. [2016] H. Wei, D. Wang, H. Zhu, J. Wang, S. Sun, X. You,  Mutual coupling calibration for multiuser massive mimo systems,  IEEE Transactions on Wireless Communications 15 (2016) 606–619. doi:10.1109/TWC.2015.2476467.
  * Zappone et al. [2019] A. Zappone, M. Di Renzo, M. Debbah,  Wireless networks design in the era of deep learning: Model-based, ai-based, or both?,  IEEE Transactions on Communications 67 (2019) 7331–7376. doi:10.1109/TCOMM.2019.2924010.
  * Guo et al. [2020] J. Guo, C.-K. Wen, S. Jin, G. Y. Li,  Convolutional neural network-based multiple-rate compressive sensing for massive mimo csi feedback: Design, simulation, and analysis,  IEEE Transactions on Wireless Communications 19 (2020) 2827–2840. doi:10.1109/TWC.2020.2968430.
  * 3GPP [2022] 3GPP,  Study on channel model for frequencies from 0.5 to 100 GHz (Tech. Rep. 3GPP TR 38.901 (V17.0.0), Mar. 2022).
  * Benzaghta and Rabie [2021] M. Benzaghta, K. M. Rabie,  Massive mimo systems for 5g: A systematic mapping study on antenna design challenges and channel estimation open issues,  IET Communications 15 (2021) 1677–1690.
  * Vasudevan [2012] K. Vasudevan,  Iterative detection of turbo-coded offset qpsk in the presence of frequency and clock offsets and awgn,  Signal, Image and Video Processing 6 (2012) 557–567.
  * Vasudevan et al. [2019] K. Vasudevan, K. Madhu, S. Singh,  Data detection in single user massive mimo using re-transmissions,  The Open Signal Processing Journal 6 (2019) 15–26.
  * Wang et al. [2018] S. Wang, K. Guan, D. He, G. Li, X. Lin, B. Ai, Z. Zhong,  Doppler shift and coherence time of 5g vehicular channels at 3.5 ghz,  in: 2018 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting, 2018, pp. 2005–2006. doi:10.1109/APUSNCURSINRSM.2018.8608340.


# Neutron-star seismology with realistic, finite-temperature nuclear matter

Fabian Gittins<EMAIL_ADDRESS>Mathematical Sciences and STAG
Research Centre, University of Southampton, Southampton SO17 1BJ, United
Kingdom  Nils Andersson<EMAIL_ADDRESS>Mathematical Sciences and
STAG Research Centre, University of Southampton, Southampton SO17 1BJ, United
Kingdom

###### Abstract

The oscillation spectrum of a neutron star is notably rich and intrinsically
dependent on the equation of state of nuclear matter. With recent advancements
in gravitational-wave and electromagnetic astronomy, we are nearing the
capability to perform neutron-star asteroseismology and probe the complex
physics of neutron stars. With this in mind, we explore the implementation of
three-parameter finite-temperature matter models in the computation of
neutron-star oscillations. We consider in detail the thermodynamics of nuclear
matter and show how this information enters the problem. Our realistic
treatment takes into account entropy and composition gradients that exist in
the nuclear matter, giving rise to buoyant g-mode oscillations. To illustrate
the implementation, we determine the oscillation spectrum of a low-temperature
neutron star. In addition to the expected compositional and thermal g-modes,
we find perturbations sourced by phase transitions in the equation of state.
We also examine the thermal evolution of the oscillation spectrum of a neutron
star with constant redshifted temperature.

##  I Introduction

Neutron stars are among the most complex objects in our Universe. A realistic
description of a neutron star inevitably requires physics across the extremes;
from Einstein’s general relativity—to describe the immense gravitational
fields—to quantum chromodynamics—in order to characterise the strong nuclear
interactions in the dense interiors. The situation is no less complicated by
the cataclysmic events that neutron stars are involved in, such as binary
mergers. Indeed, the celebrated event GW170817 [1] marked the first
gravitational-wave observation of a neutron-star coalescence and among the
main science goals with these detections is to constrain the properties of
ultra-dense nuclear matter [2, 3]. This effort will be further bolstered by
the construction of the next generation of ground-based instruments, Cosmic
Explorer [4] and the Einstein Telescope [5].

The full complexity of neutron stars is arguably best evidenced by their
remarkable oscillation spectra [6, 7]. Neutron stars possess oscillations that
are linked to the pressure (acoustic f\- and p-modes), thermal and
compositional stratification (buoyant g-modes), phase transitions (interfacial
i-modes), the presence of an elastic crust (shear s-modes) and so on. As a
guiding principle, each ingredient of physics that is added to a stellar model
gives rise to a new class of oscillation mode.

The vibrational modes of neutron stars may be relevant in a variety of
physical scenarios. There currently exists a promising connection between
theory and measured quasi-periodic oscillations in the X-ray tails of giant
magnetar flares [8, 9]. Along a similar vein, about 20 low-mass X-ray binaries
have exhibited thermonuclear burst oscillations at frequencies up to 3
different to the stellar rotation [10, 11], which have been conjectured to
arise from stellar pulsations [12, 13]. (However, the association may be
problematic as oscillation frequencies identical to the neutron star’s spin
cannot have a mode origin [14].) To complement these electromagnetic
detections, neutron-star modes may leave their mark on gravitational-wave
astronomy [15]. The gravitational-wave instability associated with the r-mode
[16, 17] has lead to recent searches [18] and oscillation modes will also be
involved with the tidal response in compact-binary inspirals [19], where the
predominant contribution comes from the f-mode and low-frequency oscillations
may be excited at a lower level. Conducting neutron-star asteroseismology with
these observational channels has the exciting potential to constrain the rich
physics associated with the modes, including the dense nuclear-matter equation
of state.

The goal of this work is to carefully develop the machinery required to
incorporate realistic, finite-temperature nuclear-matter models into
calculations of neutron-star oscillations. The bulk of the astrophysical
literature until now has focused on simple, barotropic models for the stellar
material. As we discuss here, this treatment of the matter will only apply for
equilibrium neutron stars. For dynamical contexts, we must consider departures
from the barotropic description.

There are notable exceptions to the trend of assuming one-parameter matter
models. This includes work on incorporating superfluidity [20, 21, 22, 23] and
crustal physics [7, 24] into the modelling. Additionally, there are
oscillation studies for supernovae [25] and proto-neutron stars [26], where
the matter is hot and out of equilibrium, which rely on data extracted from
general-relativistic hydrodynamical simulations. In our calculation, we
determine the oscillation modes of fluid neutron stars with realistic
stratification in the matter. But importantly the background is taken to be a
fully relativistic equilibrium. This makes the numerical computation more
tractable (and potentially more precise), as well as appropriate for isolated
neutron stars.

This paper is organised as follows. In Sec. II, we lay out the thermodynamical
framework that forms the basis of three-parameter nuclear-matter equations of
state—connecting realistic equations of state with the barotropic models that
are usually implemented. The mode calculation is detailed in Sec. III. In this
section, we explicitly show how the matter enters the mathematical formulation
of the problem. Adopting a finite-temperature nuclear-matter model in Sec. IV,
we compute the relevant thermodynamics directly from the equation of state and
present oscillation results for a low-temperature neutron star. We also show
the thermal evolution of the perturbations of a neutron star with uniform
redshifted temperature. Finally, we summarise and conclude in Sec. V.

In this article, we use a positive metric signature $(-,+,+,+)$ and adopt the
Einstein summation convention, where repeated indices denote a summation. We
use early Latin characters $a,b,c,\ldots$ for spacetime indices and later
symbols $j,k,l\ldots$ to represent spatial indices. The labels $(\ell,{\sf
m})$ represent the spherical-harmonic numbers. Unless specified otherwise, we
will assume natural units, where $G=c=k_{\text{B}}=1$.

##  II Thermodynamics and the equation of state

We begin with a discussion on the thermodynamics of the nuclear matter. This
is a natural place to start since the thermodynamics encode the microscopic
nuclear interactions and determine the material properties of neutron stars.
We want to examine the information we need in order to construct realistic,
finite-temperature stellar models and how this relates to the fundamental
thermodynamics.

Consider a small element of fluid in the neutron star from the point of view
of a local inertial reference frame that moves along with the fluid.111The
fluid element is small in the sense that its extensive variables are
infinitesimal relative to their totals across the neutron star. The
thermodynamical state of the fluid parcel is fully characterised by its
extensive variables—total energy $E$, entropy $S$, proper volume $V$ and
particle numbers $\\{N_{\text{x}}\\}$, where x labels the particle species—and
intensive variables—temperature $T$, pressure $p$ and particle chemical
potentials $\\{\mu_{\text{x}}\\}$. An infinitesimal change in the total energy
$dE$ of the fluid element is related to variations in the other extensive
variables $(dS,dV,\\{dN_{\text{x}}\\})$ through the first law of
thermodynamics. In the inertial frame, the first law can be expressed as [27]

$dE=T\,dS-p\,dV+\sum_{\text{x}}\mu_{\text{x}}\,dN_{\text{x}}.$ (1)

By virtue of the fluid forming an extensive system, we also have the
corresponding Euler relation

$E=TS-pV+\sum_{\text{x}}\mu_{\text{x}}N_{\text{x}}.$ (2)

In order to make progress in describing a neutron star, we must depart from
generality and make assumptions about the nuclear composition. For our
purposes, we will assume the fluid comprises neutrons, protons and electrons,
with respective particle numbers
$\\{N_{\text{n}},N_{\text{p}},N_{\text{e}}\\}$ and chemical potentials
$\\{\mu_{\text{n}},\mu_{\text{p}},\mu_{\text{e}}\\}$ in the fluid element. We
also assume that the fluid parcel is charge neutral
$N_{\text{p}}=N_{\text{e}}$ and contains the number of baryons
$N_{\text{b}}=N_{\text{n}}+N_{\text{p}}$. In the context of fluid dynamics, it
is convenient to express the thermodynamics in terms of densities.222In doing
so, we also reduce the number of independent variables by one. This is a
generic feature of an extensive system. Therefore, we write

$\displaystyle
d\varepsilon=n_{\text{b}}T\,ds+\frac{\varepsilon+p}{n_{\text{b}}}\,dn_{\text{b}}+n_{\text{b}}\mu_{\Delta}\,dY_{\text{e}},$
(3)
$\displaystyle\varepsilon=n_{\text{b}}Ts-p+n_{\text{b}}(\mu_{\text{n}}+\mu_{\Delta}Y_{\text{e}}),$
(4)

where $\varepsilon=E/V$ is the total energy density, $s=S/N_{\text{b}}$ is the
entropy per baryon, $n_{\text{b}}=N_{\text{b}}/V$ is the number density of
baryons, $Y_{\text{e}}=N_{\text{e}}/N_{\text{b}}$ is the electron fraction and
$\mu_{\Delta}=\mu_{\text{p}}+\mu_{\text{e}}-\mu_{\text{n}}$ quantifies
deviations from cold $\beta$-equilibrium. (It is worth noting that, at the
high temperatures and densities reached in neutron-star mergers, this relation
is no longer sufficient to describe $\beta$ equilibrium [28] but we ignore
this issue here.)

The first law (3) shows that the system possesses a fundamental potential
$\varepsilon=\varepsilon(s,n_{\text{b}},Y_{\text{e}})$ from which all the
thermodynamical information can be derived. This is illustrated in the
following. From Eq. (3), we see that the potentials $T$, $p$ and
$\mu_{\Delta}$ can be directly obtained as partial derivatives of
$\varepsilon$. These are

$T(s,n_{\text{b}},Y_{\text{e}})=\frac{1}{n_{\text{b}}}\left(\frac{\partial\varepsilon}{\partial
s}\right)_{n_{\text{b}},Y_{\text{e}}},\qquad
p(s,n_{\text{b}},Y_{\text{e}})=n_{\text{b}}^{2}\left[\frac{\partial(\varepsilon/n_{\text{b}})}{\partial
n_{\text{b}}}\right]_{s,Y_{\text{e}}},\qquad\mu_{\Delta}(s,n_{\text{b}},Y_{\text{e}})=\frac{1}{n_{\text{b}}}\left(\frac{\partial\varepsilon}{\partial
Y_{\text{e}}}\right)_{s,n_{\text{b}}}.$ (5)

Then, the final state variable $\mu_{\text{n}}$ is found using the Euler
relation (4), which can be expressed as

$\mu_{\text{n}}(s,n_{\text{b}},Y_{\text{e}})=-\frac{s}{n_{\text{b}}}\left(\frac{\partial\varepsilon}{\partial
s}\right)_{n_{\text{b}},Y_{\text{e}}}+\left(\frac{\partial\varepsilon}{\partial
n_{\text{b}}}\right)_{s,Y_{\text{e}}}-\frac{Y_{\text{e}}}{n_{\text{b}}}\left(\frac{\partial\varepsilon}{\partial
Y_{\text{e}}}\right)_{s,n_{\text{b}}}.$ (6)

This serves as a demonstration of how all the thermodynamical information sits
in $\varepsilon=\varepsilon(s,n_{\text{b}},Y_{\text{e}})$. For this reason,
$\varepsilon=\varepsilon(s,n_{\text{b}},Y_{\text{e}})$ serves as the equation
of state of the nuclear matter that connects the microphysics to the local
fluid behaviour.

In many neutron-star calculations, further assumptions about the nuclear
physics are made. One common assumption is that the matter maintains $\beta$
equilibrium, $\mu_{\Delta}=0$, which reduces the number of independent
variables by one and leads to a two-parameter equation of state
$\varepsilon=\varepsilon(s,n_{\text{b}})$. This equilibrium corresponds to
when electron-capture and neutron-decay reactions balance [28]. Chemical
equilibrium is a reasonable approximation when the dynamics are much slower
than the nuclear-reaction timescales, but will inevitably break down in the
opposite regime [29, 30]. A second common assumption is that the matter is
cold, $T=0$, such that we end up with

$d\varepsilon=\frac{\varepsilon+p}{n_{\text{b}}}\,dn_{\text{b}}=\mu_{\text{n}}\,dn_{\text{b}}$
(7)

and thus the fluid is barotropic—it depends on only one variable,
$\varepsilon=\varepsilon(n_{\text{b}})$. Under these assumptions, we can
freely invert $p=p(n_{\text{b}})$ to write the energy density as
$\varepsilon=\varepsilon[n_{\text{b}}(p)]=\varepsilon(p)$.333Another way to
retrieve this function without an inversion is to make the pressure $p$ a
natural variable of the system. This can be achieved with a Legendre transform
by introducing the enthalpy per baryon $h=(\varepsilon+p)/n_{\text{b}}$. Under
this transformation, the fundamental potential becomes
$h=h(s,p,Y_{\text{e}})$, where $p$ has replaced the baryon-number density
$n_{\text{b}}$ as an independent variable. Thence, the energy density is
determined from
$\varepsilon(s,p,Y_{\text{e}})=-p+n_{\text{b}}(s,p,Y_{\text{e}})h(s,p,Y_{\text{e}})$,
which reduces to $\varepsilon=\varepsilon(p)$ for cold, $\beta$-equilibrated
matter. The latter function is commonly encountered in the astrophysical
literature since it closes the structure equations [see Eqs. (12) below] and
is used (for example) to determine the (static) tidal Love numbers [31, 32,
33]. Indeed, recent constraints on the function $\varepsilon=\varepsilon(p)$
have been obtained from the inspiral of GW170817 [1, 2, 3], along with
complementary information obtained from X-ray data [34, 35, 36].

While this simplified view of the nuclear matter is perfectly appropriate for
describing cold, equilibrium neutron stars, it is important to note that a
relation $\varepsilon=\varepsilon(p)$ is in fact insufficient to describe the
full thermodynamical state of the star. (We illustrate and elaborate on this
point later.) Our objective should be to obtain the fundamental
thermodynamical description at the level of the equation of state.

More severely, the assumptions of $\mu_{\Delta}=T=0$ will fail for some
physical situations of interest—three of which immediately come to mind:
pulsating neutron stars (the focus of this article), binary mergers and proto-
neutron stars. In all these cases, the dynamical timescales may be
sufficiently fast such that $\beta$ equilibrium can no longer be assumed [37,
38]. Furthermore, the mergers and supernovae that we wish to examine with
gravitational-wave and electromagnetic instruments are hot, energetic events.
For these reasons, realistic calculations of neutron-star dynamics must go
beyond these assumptions.

##  III The pulsation problem

###  III.1 The equilibrium neutron star

With the thermodynamics accounted for, we are now in a position to connect the
equation of state with the stellar model. We treat the neutron star as a non-
rotating, perfect fluid. The background is spherically symmetric and its
geometry is characterised by the metric $g_{ab}$, given by the usual line
element

$ds^{2}=g_{ab}\,dx^{a}dx^{b}=-e^{\nu}\,dt^{2}+e^{\lambda}\,dr^{2}+r^{2}\,(d\theta^{2}+\sin^{2}\theta\,d\phi^{2}),$
(8)

with time coordinate $t$, spherical polar coordinates $(r,\theta,\phi)$ and
metric potentials $\nu(r)$ and $\lambda(r)$. By continuity with the exterior
Schwarzschild spacetime, the gravitational mass $m(r)$ enclosed in areal
radius $r$ is defined by

$e^{\lambda}=\frac{1}{1-2m/r}.$ (9)

In the proper reference frame, the perfect fluid is described by the
conservative stress-energy tensor

$T_{ab}=(\varepsilon+p)u_{a}u_{b}+pg_{ab},$ (10)

where the fluid four-velocity $u^{a}$ has components

$u^{t}=e^{-\nu/2},\qquad u^{j}=0.$ (11)

Finally, the Einstein field equations lead to the well-known relativistic
equations of stellar structure

$\frac{dm}{dr}=4\pi r^{2}\varepsilon,\qquad\frac{d\nu}{dr}=2\frac{m+4\pi
r^{3}p}{r(r-2m)},\qquad\frac{dp}{dr}=-\frac{\varepsilon+p}{2}\frac{d\nu}{dr}.$
(12)

Equations (12) can be integrated provided a relationship between the energy
density $\varepsilon$ and pressure $p$. From the solution to Eqs. (12), one
obtains the stellar radius $R$ and total gravitational mass $M=m(R)$, which
characterise the exterior spacetime of the spherical star.

At this point, we pick up on the issue of the equation of state (as introduced
in Sec. II) and how it enters the macroscopic description of the neutron star.
For this discussion, we will assume cold, $\beta$-equilibrated nuclear matter.
Recall that, in this regime, the fluid is barotropic and the thermodynamical
system can be fully characterised by the equation of state
$\varepsilon=\varepsilon(n_{\text{b}})$ [see Eq. (7)]. From the equation of
state, we can straightforwardly derive the thermodynamical potentials
$p=p(n_{\text{b}})$ and $\mu_{\text{n}}=\mu_{\text{n}}(n_{\text{b}})$.

Although $\varepsilon=\varepsilon(p)$ is the only relation required from the
thermodynamics to solve Eqs. (12), it is not sufficient to describe the
complete thermodynamical state of the neutron star. For example, suppose we
wanted to know the distribution of baryons $N_{\text{b}}(r)$ in the star. This
would be necessary to find the total number of baryons $N_{\text{b}}(R)$ and
the star’s baryonic mass $M_{\text{b}}=m_{\text{b}}N_{\text{b}}(R)$, where
$m_{\text{b}}$ is the mass of a baryon. To determine $N_{\text{b}}(r)$, we
must integrate the corresponding number density $n_{\text{b}}$ over the proper
volume of the star. This is equivalent to solving the differential equation

$\frac{dN_{\text{b}}}{dr}=\frac{4\pi r^{2}n_{\text{b}}}{\sqrt{1-2m/r}}.$ (13)

Here the issue is apparent: to solve Eq. (13), we require a thermodynamical
relation in addition to $\varepsilon=\varepsilon(p)$. A function
$n_{\text{b}}=n_{\text{b}}(p)$ [or equivalently
$n_{\text{b}}=n_{\text{b}}(\varepsilon)$] is necessary to close Eq. (13).
However, $n_{\text{b}}=n_{\text{b}}(p)$ is not derivable from
$\varepsilon=\varepsilon(p)$. The resolution to this is that both functions
would follow naturally from a fundamental thermodynamical potential, such as
$\varepsilon=\varepsilon(n_{\text{b}})$. Thus, while
$\varepsilon=\varepsilon(p)$ is perfectly sufficient to solve Eqs. (12), it
cannot describe the full thermodynamics of the star.

Before we move on, it is worth briefly commenting on the physics that cannot
be extracted in the barotropic regime, even with the fundamental potential
$\varepsilon=\varepsilon(n_{\text{b}})$. To illustrate this, we relax the
previous assumptions and assume a three-dimensional equation of state, where
we can freely compute the state variables. From this starting point, if we now
choose to consider zero-temperature, equilibrium matter, we may solve the
coupled, non-linear system of equations $T(s,n_{\text{b}},Y_{\text{e}})=0$ and
$\mu_{\Delta}(s,n_{\text{b}},Y_{\text{e}})=0$ to obtain $s=s(n_{\text{b}})$
and $Y_{\text{e}}=Y_{\text{e}}(n_{\text{b}})$. However, under the barotropic
assumption, we observe that the first law (7) no longer contains information
on the entropy and chemical abundances in the fluid. Hence, this information
will not be accessible to us if we only focus our attention on constraining
$\varepsilon=\varepsilon(n_{\text{b}})$.

###  III.2 The oscillations

Let us move on to the natural vibrations of the neutron star. To calculate the
oscillation modes, we will use the formulation of the perturbation equations
of Detweiler and Lindblom [39, 40], along with the augmentation introduced by
Krüger  _et al._ [7] to address numerical noise in the low-frequency spectrum.
For the exterior problem, we use the method of Andersson  _et al._ [41]. (We
refer the interested reader to these papers for the full system of equations.)
Much of this work rests on the framework of relativistic Lagrangian
perturbation theory for fluids [42, 43, 27].

The linear perturbations of a relativistic star can be described by the
Eulerian perturbation of the metric $h_{ab}$ and the Lagrangian displacement
vector of the fluid $\xi^{a}$. Since our focus is on the non-radial
oscillations of a spherically symmetric neutron star, we can decompose the
perturbations into the form $Y_{\ell}^{\sf m}e^{i\omega t}$, where $\omega$ is
the natural frequency of a mode and $Y_{\ell}^{\sf m}(\theta,\phi)$ is a
spherical harmonic with quantum numbers $(\ell,{\sf m})$. For the spacetime,
we adopt the Regge-Wheeler gauge, where high angular derivatives are chosen to
vanish [44]. In the polar sector,444It is well known that (with the exception
of the rapidly damped gravitational-wave w-modes [45]) the axial perturbations
of a non-rotating, fluid star are trivial solutions to the perturbation
problem (see, e.g., Refs. [46, 47]). Other axial oscillation modes exist in
instances where there is anisotropy in the neutron-star structure, such as
rotation [47] or the presence of an elastic crust [7]. For this calculation,
we can ignore them. the metric perturbation reads

$[h_{ab}]=-\begin{bmatrix}e^{\nu}r^{\ell}H_{0}&i\omega r^{\ell+1}H_{1}&0&0\\\
i\omega r^{\ell+1}H_{1}&e^{\lambda}r^{\ell}H_{2}&0&0\\\ 0&0&r^{\ell+2}K&0\\\
0&0&0&r^{\ell+2}\sin^{2}\theta K\end{bmatrix}Y_{\ell}^{\sf m}e^{i\omega t},$
(14)

which depends on the functions $H_{0}(r)$, $H_{1}(r)$, $H_{2}(r)$ and $K(r)$.
The displacement vector describes how the fluid elements in the equilibrium
configuration are displaced in the perturbation. We can exploit the gauge
freedom associated with the Lagrangian formulation to set $u_{a}\xi^{a}=0$
[48] and decompose the non-trivial components as

$\xi^{r}=e^{-\lambda/2}r^{\ell-1}WY_{\ell}^{\sf m}e^{i\omega
t},\qquad\xi^{\theta}=-r^{\ell-2}V\partial_{\theta}Y_{\ell}^{\sf m}e^{i\omega
t},\qquad\xi^{\phi}=-\frac{r^{\ell-2}}{\sin^{2}\theta}V\partial_{\phi}Y_{\ell}^{\sf
m}e^{i\omega t},$ (15)

with functions $W(r)$ and $V(r)$. The perturbation is described by the
spacetime $(H_{0},H_{1},H_{2},K)$ and fluid functions $(W,V)$. However, some
additional simplifications can be made.

Detweiler and Lindblom [40] reduced the pulsation equations (initially derived
by Thorne and Campolattaro [49]) for the interior of the star to a system of
four first-order differential equations for the functions $(H_{1},K,W,X)$,
where $X$ is an auxiliary variable proportional to the Lagrangian pressure
perturbation. These equations correspond to the two degrees of freedom in the
fluid and two for the gravitational radiation and are solved subject to
boundary conditions. Only two of the four linearly independent solutions
satisfy regularity at the stellar centre. At the surface, a physically
acceptable solution must have a vanishing Lagrangian perturbation of the
pressure, $X(R)=0$, which three linearly independent solutions can satisfy.
These boundary conditions are sufficient to determine $(H_{1},K,W,X)$ in the
stellar interior for any given frequency $\omega$ (up to an arbitrary
amplitude).

A generic perturbation of the star at frequency $\omega$ will have ingoing and
outgoing gravitational waves at infinity. Indeed, ingoing radiation with
frequency $\omega$ will cause the neutron star to oscillate at that frequency.
Therefore, a free oscillation mode of the star corresponds to a solution with
purely outgoing waves.

In the exterior, the functions associated with the fluid vanish and we only
have perturbations of the spacetime. The metric perturbation functions $H_{1}$
and $K$ can be combined into a single differential equation, known as the
Zerilli equation [50, 51]. This second-order equation represents the two
degrees of freedom in the gravitational field—ingoing and outgoing radiation.
An oscillation mode corresponds to when the amplitude of ingoing gravitational
waves vanishes. This defines an eigenvalue problem, with eigenfrequency
$\omega$.

Just as the background Einstein equations needed a relation connecting
$\varepsilon$ and $p$, the perturbations require a relationship between the
(Lagrangian) perturbations $\Delta\varepsilon$ and $\Delta p$. This
information comes from the thermodynamics of the nuclear matter. If we
continued our assumptions that the matter is cold and the nuclear reactions
are sufficiently fast to maintain $\beta$ equilibrium as the star oscillates,
then we can assert $\Delta p/p=\Gamma\Delta\varepsilon/(\varepsilon+p)$, where

$\Gamma=\frac{\varepsilon+p}{p}\frac{dp}{d\varepsilon}.$ (16)

(These were the assumptions adopted by Detweiler and Lindblom [39, 40].) The
index $\Gamma$ can be immediately inferred from the function
$\varepsilon=\varepsilon(p)$. However, as we have already alluded, we must be
more careful for realistic neutron stars.

It was first argued by Reisenegger and Goldreich [29] that the characteristic
timescale associated with the relevant weak interactions is actually much
longer than a typical oscillation period of a neutron star. This means that,
as a fluid element in the star is displaced from its equilibrium position
during a perturbation, the reactions do not act fast enough to restore $\beta$
equilibrium between the element and its new environment. Thus, bouyancy forces
associated with the differing composition make the fluid parcel oscillate.
These oscillations are (low-frequency) g-modes, sourced by composition
gradients.555The compositional g-modes join the spectrum alongside the well-
known class of g-modes that oscillate due to entropy gradients in the stellar
material. For a finite-temperature neutron star, we expect both effects to
contribute to the buoyancy for g-modes to exist. Therefore, for neutron stars,
it is instead more appropriate to regard the chemical composition of the fluid
as frozen during an oscillation, $\Delta Y_{\text{e}}=0$. (See Ref. [52] for a
recent study on the g-mode spectrum in the intermediate regime of finite
reaction times.)

Assuming that no heat is transferred, $\Delta s=0$, we have

$\frac{\Delta p}{p}=\Gamma_{1}\frac{\Delta
n_{\text{b}}}{n_{\text{b}}}=\Gamma_{1}\frac{\Delta\varepsilon}{\varepsilon+p},$
(17)

where

$\Gamma_{1}=\left(\frac{\partial\ln p}{\partial\ln
n_{\text{b}}}\right)_{s,Y_{\text{e}}}=\frac{\varepsilon+p}{p}\left(\frac{\partial
p}{\partial\varepsilon}\right)_{s,Y_{\text{e}}}$ (18)

is the adiabatic index, which holds information about composition and entropy
gradients in the nuclear matter. Motivated by Eq. (17), it is useful to
introduce the (relativistic) Schwarzschild discriminant

$A=\frac{1}{\varepsilon+p}\frac{d\varepsilon}{dr}-\frac{1}{\Gamma_{1}p}\frac{dp}{dr}=\left(\frac{1}{\Gamma}-\frac{1}{\Gamma_{1}}\right)\frac{d\ln
p}{dr},$ (19)

which determines the convective stability of the star [53]. If $A\leq 0$
everywhere in the fluid (equivalently, $\Gamma_{1}\geq\Gamma$), then the star
is convectively stable and supports g-mode oscillations. Otherwise, if $A=0$
(when $\Gamma_{1}=\Gamma$), the g-modes instead reside in the zero-frequency
subspace of perturbations [46, 47].

##  IV Implementing a realistic matter model

To illustrate the steps involved in incorporating a realistic description of
the nuclear matter, we adopt the finite-temperature APR model [54] as our
nuclear-matter equation of state (as implemented in the CompOSE library [55,
56, 57]). The current state of the art is based on three-parameter models with
temperature $T$, baryon-number density $n_{\text{b}}$ and fraction of
electrons $Y_{\text{e}}$ as the natural variables of the system [56].666The
equations of state in CompOSE in general use the charge fraction of strongly
interacting particles $Y_{\text{q}}$ as an independent variable. This is
because this quantity is well-defined for pure hadronic matter and extends
naturally to incorporate muons. For the composition we assume,
$Y_{\text{q}}=Y_{\text{e}}$. This is accomplished by using the (Helmholtz)
free-energy density $f=\varepsilon-n_{\text{b}}Ts$ as the fundamental
potential. With this potential, the first law (3) and Euler relation (4)
become

$\displaystyle
df=-n_{\text{b}}s\,dT+\frac{f+p}{n_{\text{b}}}\,dn_{\text{b}}+n_{\text{b}}\mu_{\Delta}\,dY_{\text{e}},$
(20) $\displaystyle
f=-p+n_{\text{b}}(\mu_{\text{n}}+\mu_{\Delta}Y_{\text{e}})$ (21)

and $f=f(T,n_{\text{b}},Y_{\text{e}})$ is the equation of state. In practice,
the relevant quantities that describe the matter’s thermodynamical state are
usually provided in tabulated form. We now describe how one calculates a
pulsating star from a three-parameter equation-of-state table.

We start from the table, which includes the pressure $p$, entropy per baryon
$s$, deviation from $\beta$ equilibrium $\mu_{\Delta}$ and total energy
density $\varepsilon$ on a grid of $(T,n_{\text{b}},Y_{\text{e}})$. We
interpolate this information to find the functions
$p=p(T,n_{\text{b}},Y_{\text{e}})$, $s=s(T,n_{\text{b}},Y_{\text{e}})$,
$\mu_{\Delta}=\mu_{\Delta}(T,n_{\text{b}},Y_{\text{e}})$ and
$\varepsilon=\varepsilon(T,n_{\text{b}},Y_{\text{e}})$.777To account for the
wide spread in magnitudes, we actually use logarithms to interpolate the
quantities $\ln p$, $\ln s$, $\mu_{\Delta}$ and $\ln\varepsilon$ as functions
of $(\ln T,\ln n_{\text{b}},Y_{\text{e}})$ using splines. For the APR model
that we work with, the grid is given as $0.00999\leq T/\unit{\mega}\leq 251$,
$$5.43\text{\times}{10}^{-13}$\leq n_{\text{b}}/\unit{\per\femto\cubed}\leq
6.32$ and $0.005\leq Y_{\text{e}}\leq 0.655$. In addition to the functions, we
also require their partial derivatives to determine the indices
$\Gamma=\Gamma(T,n_{\text{b}},Y_{\text{e}})$ and
$\Gamma_{1}=\Gamma_{1}(T,n_{\text{b}},Y_{\text{e}})$. We use Eq. (16) to write

$\Gamma=\frac{\varepsilon+p}{p}\frac{dp/dn_{\text{b}}}{d\varepsilon/dn_{\text{b}}},$
(22)

where the derivatives are given by

$\displaystyle\frac{dp}{dn_{\text{b}}}$ $\displaystyle=\left(\frac{\partial
p}{\partial
T}\right)_{n_{\text{b}},Y_{\text{e}}}\frac{dT}{dn_{\text{b}}}+\left(\frac{\partial
p}{\partial n_{\text{b}}}\right)_{T,Y_{\text{e}}}+\left(\frac{\partial
p}{\partial
Y_{\text{e}}}\right)_{T,n_{\text{b}}}\frac{dY_{\text{e}}}{dn_{\text{b}}},$
(23) $\displaystyle\frac{d\varepsilon}{dn_{\text{b}}}$
$\displaystyle=\left(\frac{\partial\varepsilon}{\partial
T}\right)_{n_{\text{b}},Y_{\text{e}}}\frac{dT}{dn_{\text{b}}}+\left(\frac{\partial\varepsilon}{\partial
n_{\text{b}}}\right)_{T,Y_{\text{e}}}+\left(\frac{\partial\varepsilon}{\partial
Y_{\text{e}}}\right)_{T,n_{\text{b}}}\frac{dY_{\text{e}}}{dn_{\text{b}}}.$
(24)

For Eq. (18), we use $\Delta s=\Delta Y_{\text{e}}=0$ to infer

$\Gamma_{1}=\frac{n_{\text{b}}}{p}\left[\left(\frac{\partial p}{\partial
n_{\text{b}}}\right)_{T,Y_{\text{e}}}-\left(\frac{\partial p}{\partial
T}\right)_{n_{\text{b}},Y_{\text{e}}}\left(\frac{\partial s}{\partial
n_{\text{b}}}\right)_{T,Y_{\text{e}}}\left(\frac{\partial s}{\partial
T}\right)_{n_{\text{b}},Y_{\text{e}}}^{-1}\right].$ (25)

We obtain the partial thermodynamic derivatives in $\Gamma$ and $\Gamma_{1}$
from the interpolation functions.888It is worth noting that second derivatives
of the fundamental potential are not independent. They are related via the
Maxwell relations. The Maxwell relation that simplifies Eq. (25) is $(\partial
p/\partial T)_{n_{\text{b}},Y_{\text{e}}}=-n_{\text{b}}^{2}(\partial
s/\partial n_{\text{b}})_{T,Y_{\text{e}}}$. For this reason, it is important
to use an interpolation routine that ensures continuity of the functions and
their derivatives in the grid. It should be noted that the partial derivatives
are not explicitly provided in the table. Thus, a numerical technique to
approximate the derivatives (like the one we use) is necessary.

Next, we assume the background star is in chemical equilibrium to determine
the composition. This is carried out by solving

$\mu_{\Delta}(T,n_{\text{b}},Y_{\text{e}})=0$ (26)

to obtain the function $Y_{\text{e}}=Y_{\text{e}}(T,n_{\text{b}})$. Although
Eq. (26) represents a simple root-search, it must be carried out at every
integration step in the background and perturbation equations. Thus, to speed
up the computations, we determine $Y_{\text{e}}$ at each $(T,n_{\text{b}})$ in
the original grid and pass it to the interpolation functions—e.g.,
$p=p[T,n_{\text{b}},Y_{\text{e}}(T,n_{\text{b}})]=p(T,n_{\text{b}})$—to form a
two-dimensional table for the matter model. (A similar approach is adopted in
some numerical-relativity simulations of neutron-star mergers.) For the APR
model, we find that, for certain values of $(T,n_{\text{b}})$, $\mu_{\Delta}$
does not change sign with varying $Y_{\text{e}}$. We therefore restricted the
two-dimensional grid to $0.01\lesssim T/\unit{\mega}\lesssim 230$ and
$$1.3\text{\times}{10}^{-7}$\lesssim
n_{\text{b}}/\unit{\per\femto\cubed}\lesssim 2.5$, where equilibrium can be
enforced.

Then we interpolate the quantities on the new, two-dimensional grid. To
calculate the equilibrium star, we invert $p=p(T,n_{\text{b}})$ to obtain
$n_{\text{b}}=n_{\text{b}}(T,p)$. Given the temperature $T$, this is
sufficient information to integrate the structure equations (12), as well as
Eq. (13). The perturbation equations need $\Gamma=\Gamma(T,n_{\text{b}})$ and
$\Gamma_{1}=\Gamma_{1}(T,n_{\text{b}})$. To determine $\Gamma$, we require

$\frac{dY_{\text{e}}}{dn_{\text{b}}}=\left(\frac{\partial
Y_{\text{e}}}{\partial
T}\right)_{n_{\text{b}}}\frac{dT}{dn_{\text{b}}}+\left(\frac{\partial
Y_{\text{e}}}{\partial n_{\text{b}}}\right)_{T},$ (27)

which is inferred from $Y_{\text{e}}=Y_{\text{e}}(T,n_{\text{b}})$. This
concludes the description of how a three-parameter matter model enters the
perturbation problem.

###  IV.1 A cold neutron star

We are now in a position to determine the oscillation modes of a neutron star
described by the three-parameter APR model [54]. First, we will explore the
$l=2$ oscillation spectrum of a neutron star of $M_{\text{b}}=1.4M_{\odot}$
with constant $T=\qty{0.02}{\mega}$. Since the temperature is uniform, we have
$dT/dn_{\text{b}}=0$. Due to the low-temperature cut-off in the table (which
is a typical feature in nuclear-matter models), we must assume temperatures
above that of mature neutron stars (below
$T\sim\qty{0.1}{\kilo}\sim\qty{e6}{}$). We did not consider a zero-temperature
model, since these equations of state generally do not contain information on
$s$, which is needed to compute $\Gamma_{1}$ [see Eq. (25)], and partial
derivatives with respect to $T$. Still, at the chosen temperature the thermal
pressure is fairly weak compared to the fluid contribution. We show in Fig. 1
the energy density and adiabatic indices of the equilibrium neutron star,
indicating the core-crust phase transition and the location of neutron drip.

Figure 1: The profiles of the energy density (left panel) and adiabatic
indices (right panel) in the cold neutron star described by the APR equation
of state. The two vertical, dashed lines in each panel correspond to the core-
crust phase transition at high density (right) and neutron drip at low density
(left). The core-crust transition produces a discontinuity in the density,
which results in sharp behaviour in the indices. Neutron drip sources a change
in the density derivative, corresponding to the second sharp feature (at lower
densities) in the background index $\Gamma$.

Our calculation is carried out in full general relativity and therefore
gravitational radiation is emitted by the modes. The waves carry away energy
and damp the oscillations, which result in the eigensolutions being complex.
Although we allow for complex numbers in our numerical determination of the
perturbations, we will only show the real parts in our results, since the
g-modes we focus on are weakly damped and thus have small imaginary
contributions to the frequency.

In Fig. 2, we show a visualisation of the oscillation spectrum for quadrupolar
perturbations. The figure presents the amplitude of ingoing radiation
$\tilde{A}_{\text{in}}$ (with a particular normalisation to flatten out the
spectrum) versus real-valued frequencies [41]. A mode solution corresponds to
zeros of $\tilde{A}_{\text{in}}$ and appear as sharp singularities in Fig. 2.

Figure 2: The oscillation spectrum for quadrupolar perturbations of the cold
neutron star, visualised with the ingoing gravitational-wave amplitude against
(real) oscillation frequency. The left panel focuses on a low-frequency
portion of the spectrum, while the right panel displays high frequencies. A
narrow singularity corresponds to a slowly damped mode solution. The low-
frequency singularities represent the g-modes of the star, restored by
stratification in the thermal and compositional gradients, and the core-crust
phase transition. Vertical, dotted lines indicate the real frequencies of the
oscillation modes summarised in Table 1. The two highest frequency g-modes are
not visible in the low-frequency range, but can be seen in the right panel.
The right panel shows the f-mode and the first few p-modes at successively
higher frequencies.

With the spectrum visualised, we are in a position to precisely determine the
eigensolutions of the oscillation modes. We start by searching for the
fundamental f-mode of the star. The f-mode is identified by its broad
singularity in the spectrum (see the right panel of Fig. 2). We find it with
frequency $\text{Re}(\omega)/(2\pi)=\qty{1.9123}{\kilo}$ and show its radial
eigenfunction $W$ in Fig. 3. As is characteristic of f-modes, the
eigenfunction gradually increases away from the centre and is largest at the
stellar surface.

Figure 3: The radial-displacement eigenfunction of the quadrupolar f-mode of
the cold neutron star. The eigenfunction has no nodes and reaches its maximum
at the stellar surface. The f-mode sources large (Eulerian) density
disturbances in the fluid and couples strongly with outgoing gravitational
radiation.

Next, we examine the low-frequency perturbations. We show the eigenfunctions
of the first four g-modes in Fig. 4. The first g-mode—the g1-mode—has
frequency $\text{Re}(\omega)/(2\pi)=\qty{508.5}{}$. The sharp feature in its
eigenfunction indicates its association with the energy-density discontinuity
at the core-crust phase transition (around
$\varepsilon/c^{2}\sim\qty{e14}{\per\cubic\centi}$), shown in the left panel
of Fig. 1. The frequency for this oscillation mode is higher than a typical
g-mode. We will return to this point in a moment.

Figure 4: The radial-displacement eigenfunction of the first four quadrupolar
g-modes of the cold neutron star. The vertical, dashed lines correspond to the
features indicated in Fig. 1. The g1-mode eigenfunction has a particularly
prominent spike at the energy-density discontinuity (the core-crust phase
transition). Some of the other g-modes have notable features at neutron drip
in the crust.

The g-mode perturbations arise when there is stratification in the fluid. As
we have discussed, this naturally occurs as the entropy and chemical
composition vary in the stellar profile. However, there is another form of
stratification that is expected to occur in neutron stars: phase transitions
that lead to discontinuities in the energy density. Such a phase transition
will exist at the interface between the fluid core and elastic crust, and may
arise in the high-density core if there is, for example, quark deconfinement.
These interfaces lead to the presence of discontinuity g-modes, also commonly
known as i-modes [58, 6, 59]. To clarify this identification, we have also
computed the mode spectrum when the stratification is artificially removed by
setting $\Gamma_{1}=\Gamma$. For the unstratified stellar model, the high-
frequency portion is largely unchanged; the f-mode moves to
$\text{Re}(\omega)/(2\pi)=\qty{1.9122}{\kilo}$. However, by removing the
composition gradients, the entire low-frequency spectrum is eliminated; the
only exception to this is a solution with
$\text{Re}(\omega)/(2\pi)=\qty{471.5}{}$. We see from its eigenfunction, shown
in Fig. 5, that it has two features that directly correspond to the energy-
density profile, including the sharp feature at the core-crust boundary. Thus,
the i-mode of the unstratified star and the g1-mode likely have the same
origin. However, their eigenfrequencies substantially differ. We will continue
referring to the oscillation of the stratified neutron star as the g1-mode.

Finn [58] derived a formula to estimate the typical frequencies of
discontinuity modes. He found that the mode eigenfrequency is proportional to
the density jump and the depth at which the interface occurs. We observe in
the left panel of Fig. 1 that the discontinuity that gives rise to the mode is
large. For comparison, we computed the corresponding mode using the DD2
nuclear-matter model [60] at $T=\qty{0.2}{\mega}$. With this equation of
state, the density discontinuity is smaller and the mode oscillates at a lower
frequency, $\text{Re}(\omega)/(2\pi)=\qty{285.5}{}$. Hence, the large
frequency of the g1-mode can be explained by the relatively large
discontinuity in the APR equation of state.

Figure 5: The radial-displacement eigenfunction of the quadrupolar i-mode of
the unstratified, cold neutron star. The vertical, dashed lines correspond to
the features indicated in Fig. 1. The eigenfunction is qualitatively similar
to that of the g1-mode encountered in the realistic, stratified neutron star
shown in Fig. 4.

Another (less prominent) feature is seen in the g2-mode eigenfunction at a
different point in the star. This mode has frequency
$\text{Re}(\omega)/(2\pi)=\qty{304.5}{}$ and the sign change in its
eigenfunction derivative corresponds to neutron drip in the crust at
$\varepsilon/c^{2}\sim\qty{4e11}{\per\cubic\centi}$. At neutron drip, Fig. 1
shows that the derivative of the energy density is discontinuous, which is
seen clearly in the profile of $\Gamma$. The eigenfunctions in Fig. 4 show
that other g-modes are also sensitive to this feature.

The left panel of Fig. 2 shows the spectrum of the rest of the g-modes. We see
that these oscillations are consistent with typical g-mode frequencies of
$\text{Re}(\omega)/(2\pi)\sim\qtyrange{10}{100}{}$ [29, 19, 7]. We summarise
the oscillations of the cold neutron star in Table 1.

Table 1: The (real) eigenfrequencies of the quadrupolar f-mode and first several g-modes of the cold neutron star described by the APR equation of state. |  $\text{Re}(\omega)/(2\pi)$ /
---|---
f-mode | $1912.3$
g1-mode | $508.5$
g2-mode | $304.5$
g3-mode | $126.1$
g4-mode | $116.6$
g5-mode | $112.6$
g6-mode | $89.1$
g7-mode | $72.1$
g8-mode | $71.7$
g9-mode | $61.8$

###  IV.2 Uniform redshifted temperature

Before we conclude, we shall briefly lay out the steps in implementing a
realistic thermal treatment. In general relativity, stars are not expected to
have constant temperature. Instead, isolated stars evolve towards uniform
redshifted temperatures, where $Te^{\nu/2}=\text{const}$ [61]. Since a typical
neutron star is expected to be formed from the remnant of a supernova, the
study of the thermal evolution and its impact on the neutron-star mode
spectrum has attracted some attention [26, 7].

In order to compute the stellar equilibrium, we now specify the redshifted
temperature and use the shooting method to satisfy the boundary condition on
the metric potential at the surface, $e^{\nu(R)}=1-2M/R$. The star’s
temperature profile results in a non-zero $dT/dn_{\text{b}}$. To determine
this quantity, we first note that

$\frac{dT}{dn_{\text{b}}}=-\frac{T}{2}\frac{d\nu}{dn_{\text{b}}}=\frac{T}{\varepsilon+p}\frac{dp}{dn_{\text{b}}}.$
(28)

Then, combining Eqs. (23) and (27), we obtain

$\frac{dT}{dn_{\text{b}}}=\left[\left(\frac{\partial p}{\partial
n_{\text{b}}}\right)_{T,Y_{\text{e}}}+\left(\frac{\partial p}{\partial
Y_{\text{e}}}\right)_{T,n_{\text{b}}}\left(\frac{\partial
Y_{\text{e}}}{\partial
n_{\text{b}}}\right)_{T}\right]\left[\frac{\varepsilon+p}{T}-\left(\frac{\partial
p}{\partial T}\right)_{n_{\text{b}},Y_{\text{e}}}-\left(\frac{\partial
p}{\partial Y_{\text{e}}}\right)_{T,n_{\text{b}}}\left(\frac{\partial
Y_{\text{e}}}{\partial T}\right)_{n_{\text{b}}}\right]^{-1}.$ (29)

This information is required to determine $\Gamma$.

We calculated the perturbations of a neutron star with
$Te^{\nu/2}=\qty{0.02}{\mega}$ and found very few differences as compared to
the cold constant-temperature model. The eigenfrequencies as reported in Table
1 are the same to the accuracy provided.

Although these temperatures are already much higher than that of a mature,
isolated neutron star, it is informative to consider how the oscillation modes
vary. To see this, we determined the frequencies of the f-mode and first three
g-modes as functions of the redshifted temperature in the range
$Te^{\nu/2}=\qtyrange{0.02}{15}{\mega}$. This sequence of neutron stars share
the same the baryon mass of $M_{\text{b}}=1.4M_{\odot}$. As the star is heated
up, it goes from an equilibrium radius of $R=\qty{11.6}{\kilo}$ up to
$R=\qty{16.6}{\kilo}$. The oscillation frequencies are displayed in Fig. 6.

Figure 6: The (real) quadrupolar mode eigenfrequencies of the neutron star
against redshifted temperature. The eigenfrequencies correspond to the f\- and
first three g-modes. The f-mode is relatively insensitive to the temperature,
even up to $Te^{\nu/2}\sim\qty{1}{\mega}$. In contrast, the g-modes change
drastically as the higher temperatures source larger entropy gradients.

As expected, the frequency of the f-mode—approximately proportional to
$\sqrt{M/R^{3}}$—decreases as the star expands. Figure 6 shows that the f-mode
frequency is relatively insensitive to the temperature up to
$Te^{\nu/2}\sim\qty{1}{\mega}$. The situation is quite different for the
g-modes.

Thermodynamically, temperature and entropy are conjugates of each other. Thus,
as the star is heated, the entropy gradients in the fluid become larger. This
provides a stronger restoring force for buoyant perturbations and causes the
g-modes to oscillate more rapidly. This is consistent with the behaviour seen
in Fig. 6. These results are in qualitative agreement with other studies [26,
7].

##  V Conclusions

We have explored how realistic, finite-temperature nuclear matter enters the
computation of neutron-star oscillation modes. In paying close attention to
the thermodynamics, we have illustrated the implicit assumptions with one-
parameter (barotropic) matter models; zero temperature $T=0$ and the particle
species maintaining chemical equilibrium $\mu_{\Delta}=0$. These assumptions
will not hold for dynamical situations, such as oscillating neutron stars,
proto-neutron stars and during compact-binary coalescences.

With a view to constructing realistic neutron-star models, we have focused our
attention on three-parameter equations of state, which represents the present
state of the art in nuclear astrophysics. With these equations of state, the
stellar material can support entropy and composition gradients, giving rise to
low-frequency g-mode oscillations.

As a demonstration, we adopted the APR equation of state with finite
temperature and examined the resultant oscillation spectrum for
$T=\qty{0.02}{\mega}$. APR nuclear matter involves a discontinuity in the
energy density that corresponds with the phase transition at the core-crust
boundary. We identified a discontinuity g-mode sourced by this interface in
the star. The likely explanation for its high frequency
$\text{Re}(\omega)/(2\pi)=\qty{508.5}{}$ is the large density jump in the APR
equation of state. We also calculated the first few g-modes that arise from
composition and entropy stratification in the fluid.

Realistic, mature neutron stars are expected to possess temperature profiles
that are uniformly redshifted. We detailed how the oscillation spectrum of
such a star may be computed and showed how the eigenfrequencies of the f-mode
and first three g-modes vary for different temperatures.

As we wrap up, it is worth making some final comments on the implementation of
the thermodynamics. Matter models are typically provided as tables, owing to
the numerical complexity of computing nuclear interactions at high densities.
Since neutron-star models require this information in functional form, we must
resort to numerical interpolation of the tabular data. However, as it stands,
there is no codified standard for how to conduct this interpolation and
different choices for the order and method, as well as the quantities to
interpolate, lead to quantitatively different results for the macroscopic
stellar model. It is worth noting that in the past there has been work on
developing interpolation methods that are thermodynamically consistent [62].
More effort in this direction, aimed at meeting the requirements of seismology
work, would be welcome.

Additionally, the oscillation modes depend on thermodynamic derivatives—such
as the first partial derivatives of $p=p(T,n_{\text{b}},Y_{\text{e}})$—that
are not often explicitly provided in equation-of-state tables. When this is
the case, numerical differentiation is the only option. We determined this
information by interpolating the quantities using splines and used the spline
coefficients to estimate the derivatives. As can be seen in Fig. 1, this gives
rise to singular features in the neutron-star profile, which leave their mark
on the oscillation spectrum. With the goal of conducting reliable neutron-star
asteroseismology in the era of precision astronomy, these are important issues
to address.

###### Acknowledgements.

NA acknowledges support from STFC via grant number ST/R00045X/1. The authors
are grateful to Ian Hawke for useful discussions on the numerical component of
this work and Christian Krüger for sharing code to compare against. The
software developed to support this article is available in a GitHub repository
[63] and is written in the Julia programming language [64, 65, 66, 67, 68].
The figures were generated using Matplotlib [69, 70].

## References

  * Abbott  _et al._ [2017] B. P. Abbott _et al._ , GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017).
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ , GW170817: Measurements of Neutron Star Radii and Equation of State, Phys. Rev. Lett. 121, 161101 (2018).
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ , Properties of the Binary Neutron Star Merger GW170817, Physical Review X 9, 011001 (2019).
  * Reitze  _et al._ [2019] D. Reitze _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, arXiv:1907.04833 [astro-ph.IM] (2019).
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The Einstein Telescope: a third-generation gravitational wave observatory, Classical Quant. Grav. 27, 194002 (2010).
  * McDermott  _et al._ [1988] P. N. McDermott, H. M. van Horn, and C. J. Hansen, Nonradial Oscillations of Neutron Stars, Astrophys. J. 325, 725 (1988).
  * Krüger  _et al._ [2015] C. J. Krüger, W. C. G. Ho, and N. Andersson, Seismology of adolescent neutron stars: Accounting for thermal effects and crust elasticity, Phys. Rev. D 92, 063009 (2015).
  * Watts and Strohmayer [2006] A. L. Watts and T. E. Strohmayer, Detection with RHESSI of High-Frequency X-Ray Oscillations in the Tail of the 2004 Hyperflare from SGR 1806-20, Astrophys. J. 637, L117 (2006).
  * Samuelsson and Andersson [2007] L. Samuelsson and N. Andersson, Neutron star asteroseismology. Axial crust oscillations in the Cowling approximation, Mon. Not. R. Astron. Soc. 374, 256 (2007).
  * Watts [2012] A. L. Watts, Thermonuclear Burst Oscillations, Annu. Rev. Astron. Astrophys. 50, 609 (2012).
  * Watts and Kini [2021] A. L. Watts and Y. Kini, Thermonuclear burst oscillations from accreting neutron stars: catalogue of confirmed sources (2021).
  * Bildsten  _et al._ [1996] L. Bildsten, G. Ushomirsky, and C. Cutler, Ocean g-Modes on Rotating Neutron Stars, Astrophys. J. 460, 827 (1996).
  * Heyl [2004] J. S. Heyl, r-Modes on Rapidly Rotating, Relativistic Stars. I. Do Type I Bursts Excite Modes in the Neutron Star Ocean?, Astrophys. J. 600, 939 (2004).
  * Gittins  _et al._ [2023] F. Gittins, T. Celora, A. Beri, and N. Andersson, Modelling Neutron-Star Ocean Dynamics, Universe 9, 226 (2023).
  * Andersson and Kokkotas [1998] N. Andersson and K. D. Kokkotas, Towards gravitational wave asteroseismology, Mon. Not. R. Astron. Soc. 299, 1059 (1998).
  * Andersson [1998] N. Andersson, A New Class of Unstable Modes of Rotating Relativistic Stars, Astrophys. J. 502, 708 (1998).
  * Friedman and Morsink [1998] J. L. Friedman and S. M. Morsink, Axial Instability of Rotating Relativistic Stars, Astrophys. J. 502, 714 (1998).
  * Abbott  _et al._ [2021] R. Abbott _et al._ , Constraints from LIGO O3 Data on Gravitational-wave Emission Due to R-modes in the Glitching Pulsar PSR J0537–6910, Astrophys. J. 922, 71 (2021).
  * Lai [1994] D. Lai, Resonant Oscillations and Tidal Heating in Coalescing Binary Neutron Stars, Mon. Not. R. Astron. Soc. 270, 611 (1994).
  * Andersson and Comer [2001] N. Andersson and G. L. Comer, On the dynamics of superfluid neutron star cores, Mon. Not. R. Astron. Soc. 328, 1129 (2001).
  * Andersson  _et al._ [2002] N. Andersson, G. L. Comer, and D. Langlois, Oscillations of general relativistic superfluid neutron stars, Phys. Rev. D 66, 104002 (2002).
  * Prix  _et al._ [2004] R. Prix, G. L. Comer, and N. Andersson, Inertial modes of non-stratified superfluid neutron stars, Mon. Not. R. Astron. Soc. 348, 625 (2004).
  * Lin  _et al._ [2008] L.-M. Lin, N. Andersson, and G. L. Comer, Oscillations of general relativistic multifluid/multilayer compact stars, Phys. Rev. D 78, 083008 (2008).
  * Andersson  _et al._ [2019] N. Andersson, B. Haskell, G. L. Comer, and L. Samuelsson, The dynamics of neutron star crusts: Lagrangian perturbation theory for a relativistic superfluid-elastic system, Classical Quant. Grav. 36, 105004 (2019).
  * Morozova  _et al._ [2018] V. Morozova, D. Radice, A. Burrows, and D. Vartanyan, The Gravitational Wave Signal from Core-collapse Supernovae, Astrophys. J. 861, 10 (2018).
  * Ferrari  _et al._ [2003] V. Ferrari, G. Miniutti, and J. A. Pons, Gravitational waves from newly born, hot neutron stars, Mon. Not. R. Astron. Soc. 342, 629 (2003).
  * Andersson and Comer [2021] N. Andersson and G. L. Comer, Relativistic fluid dynamics: physics for many different scales, Living Rev. Relativity 24, 3 (2021).
  * Alford and Harris [2018] M. G. Alford and S. P. Harris, $\beta$ equilibrium in neutron-star mergers, Phys. Rev. C 98, 065806 (2018).
  * Reisenegger and Goldreich [1992] A. Reisenegger and P. Goldreich, A New Class of g-Modes in Neutron Stars, Astrophys. J. 395, 240 (1992).
  * Andersson and Pnigouras [2019] N. Andersson and P. Pnigouras, The g-mode spectrum of reactive neutron star cores, Mon. Not. R. Astron. Soc. 489, 4043 (2019).
  * Hinderer [2008] T. Hinderer, Tidal Love Numbers of Neutron Stars, Astrophys. J. 677, 1216 (2008).
  * Binnington and Poisson [2009] T. Binnington and E. Poisson, Relativistic theory of tidal Love numbers, Phys. Rev. D 80, 084018 (2009).
  * Damour and Nagar [2009] T. Damour and A. Nagar, Relativistic tidal properties of neutron stars, Phys. Rev. D 80, 084035 (2009).
  * Raaijmakers  _et al._ [2019] G. Raaijmakers, T. E. Riley, A. L. Watts, S. K. Greif, S. M. Morsink, K. Hebeler, A. Schwenk, T. Hinderer, S. Nissanke, S. Guillot, Z. Arzoumanian, S. Bogdanov, D. Chakrabarty, K. C. Gendreau, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, and M. T. Wolff, A Nicer View of PSR J0030+0451: Implications for the Dense Matter Equation of State, Astrophys. J. 887, L22 (2019).
  * Raaijmakers  _et al._ [2020] G. Raaijmakers, S. K. Greif, T. E. Riley, T. Hinderer, K. Hebeler, A. Schwenk, A. L. Watts, S. Nissanke, S. Guillot, J. M. Lattimer, and R. M. Ludlam, Constraining the Dense Matter Equation of State with Joint Analysis of NICER and LIGO/Virgo Measurements, Astrophys. J. 893, L21 (2020).
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, Constraints on the Dense Matter Equation of State and Neutron Star Properties from NICER’s Mass-Radius Estimate of PSR J0740+6620 and Multimessenger Observations, Astrophys. J. 918, L29 (2021).
  * Hammond  _et al._ [2021] P. Hammond, I. Hawke, and N. Andersson, Thermal aspects of neutron star mergers, Phys. Rev. D 104, 103006 (2021).
  * Hammond  _et al._ [2023] P. Hammond, I. Hawke, and N. Andersson, Impact of nuclear reactions on gravitational waves from neutron star mergers, Phys. Rev. D 107, 043023 (2023).
  * Lindblom and Detweiler [1983] L. Lindblom and S. L. Detweiler, The quadrupole oscillations of neutron stars., Astrophys. J. Suppl. Ser. 53, 73 (1983).
  * Detweiler and Lindblom [1985] S. Detweiler and L. Lindblom, On the nonradial pulsations of general relativistic stellar models, Astrophys. J. 292, 12 (1985).
  * Andersson  _et al._ [1995] N. Andersson, K. D. Kokkotas, and B. F. Schutz, A new numerical approach to the oscillation modes of relativistic stars, Mon. Not. R. Astron. Soc. 274, 1039 (1995).
  * Friedman [1978] J. L. Friedman, Generic instability of rotating relativistic stars, Commun. Math. Phys. 62, 247 (1978).
  * Friedman and Stergioulas [2013] J. L. Friedman and N. Stergioulas,  _Rotating Relativistic Stars_ (Cambridge University Press, Cambridge, United Kingdom, 2013).
  * Regge and Wheeler [1957] T. Regge and J. A. Wheeler, Stability of a Schwarzschild Singularity, Phys. Rev. 108, 1063 (1957).
  * Kokkotas [1994] K. D. Kokkotas, Axial Modes for Relativistic Stars, Mon. Not. R. Astron. Soc. 268, 1015 (1994).
  * Lockitch and Friedman [1999] K. H. Lockitch and J. L. Friedman, Where are the r-Modes of Isentropic Stars?, Astrophys. J. 521, 764 (1999).
  * Lockitch  _et al._ [2000] K. H. Lockitch, N. Andersson, and J. L. Friedman, Rotational modes of relativistic stars: Analytic results, Phys. Rev. D 63, 024019 (2000).
  * Schutz and Sorkin [1977] B. F. Schutz and R. Sorkin, Variational aspects of relativistic field theories, with application to perfect fluids., Ann. Phys. 107, 1 (1977).
  * Thorne and Campolattaro [1967] K. S. Thorne and A. Campolattaro, Non-Radial Pulsation of General-Relativistic Stellar Models. I. Analytic Analysis for $L\geq 2$, Astrophys. J. 149, 591 (1967).
  * Zerilli [1970] F. J. Zerilli, Gravitational Field of a Particle Falling in a Schwarzschild Geometry Analyzed in Tensor Harmonics, Phys. Rev. D 2, 2141 (1970).
  * Fackerell [1971] E. D. Fackerell, Solutions of Zerilli’s Equation for Even-Parity Gravitational Perturbations, Astrophys. J. 166, 197 (1971).
  * Counsell  _et al._ [2024] A. R. Counsell, F. Gittins, and N. Andersson, The impact of nuclear reactions on the neutron-star g-mode spectrum, Mon. Not. R. Astron. Soc. 531, 1721 (2024).
  * Detweiler and Ipser [1973] S. L. Detweiler and J. R. Ipser, A Variational Principle and a Stability Criterion for the Non-Radial Modes of Pulsation of Stellar Models in General Relativity, Astrophys. J. 185, 685 (1973).
  * Schneider  _et al._ [2019] A. S. Schneider, C. Constantinou, B. Muccioli, and M. Prakash, Akmal-Pandharipande-Ravenhall equation of state for simulations of supernovae, neutron stars, and binary mergers, Phys. Rev. C 100, 025803 (2019).
  * Typel  _et al._ [2015] S. Typel, M. Oertel, and T. Klähn, CompOSE CompStar online supernova equations of state harmonising the concert of nuclear physics and astrophysics compose.obspm.fr, Physics of Particles and Nuclei 46, 633 (2015).
  * Oertel  _et al._ [2017] M. Oertel, M. Hempel, T. Klähn, and S. Typel, Equations of state for supernovae and compact stars, Rev. Mod. Phys. 89, 015007 (2017).
  * Typel  _et al._ [2022] S. Typel, M. Oertel, T. Klähn, D. Chatterjee, V. Dexheimer, C. Ishizuka, M. Mancini, J. Novak, H. Pais, C. Providencia, A. Raduta, M. Servillat, and L. Tolos, CompOSE Reference Manual, arXiv:2203.03209 [astro-ph.HE] (2022).
  * Finn [1987] L. S. Finn, g-modes in zero-temperature neutron stars, Mon. Not. R. Astron. Soc. 227, 265 (1987).
  * Tsang  _et al._ [2012] D. Tsang, J. S. Read, T. Hinderer, A. L. Piro, and R. Bondarescu, Resonant Shattering of Neutron Star Crusts, Phys. Rev. Lett. 108, 011102 (2012).
  * Hempel and Schaffner-Bielich [2010] M. Hempel and J. Schaffner-Bielich, A statistical model for a complete supernova equation of state, Nucl. Phys. A837, 210 (2010).
  * Thorne [1967] K. S. Thorne, The General Relativistic Theory of Stellar Structure and Dynamics, in  _High Energy Astrophysics, Volume 3_ , Vol. 3, edited by C. Dewitt, E. Schatzman, and P. Véron (1967) p. V.
  * Haensel and Proszynski [1982] P. Haensel and M. Proszynski, Pion condensation in cold dense matter and neutron stars, Astrophys. J. 258, 306 (1982).
  * Gittins [2024] F. Gittins, RealisticSeismology (2024).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Rev. 59, 65 (2017).
  * Rackauckas and Nie [2017] C. Rackauckas and Q. Nie, DifferentialEquations.jl–A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia, J. Open Res. Software 5, 15 (2017).
  * Kittisopikul  _et al._ [2022] M. Kittisopikul, T. E. Holy, and T. Aschan, Interpolations.jl (2022).
  * Pal  _et al._ [2024] A. Pal, F. Holtorf, A. Larsson, T. Loman, Utkarsh, F. Schäefer, Q. Qu, A. Edelman, and C. Rackauckas, NonlinearSolve.jl: High-Performance and Robust Solvers for Systems of Nonlinear Equations in Julia, arXiv:2403.16341 [math.NA] (2024).
  * Mogensen and Riseth [2018] P. K. Mogensen and A. N. Riseth, Optim: A mathematical optimization package for Julia, J. Open Source Software 3, 615 (2018).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng. 9, 90 (2007).
  * Rowley [2022] C. Rowley, PythonCall.jl (2022).