# Prandtl number effects on extreme mixing events in forced stratified
turbulence

Nicolaos Petropoulos1<EMAIL_ADDRESS>Miles M. P. Couchman2  Ali Mashayek3
Stephen M. de Bruyn Kops4  Colm-cille P. Caulfield5,1  1Department of Applied
Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3
0WA, UK 2Department of Mathematics and Statistics, York University, Toronto,
ON M3J 1P3, Canada 3Department of Earth Sciences, University of Cambridge,
Cambridge CB2 3EQ, UK 4Department of Mechanical and Industrial Engineering,
University of Massachusetts Amherst, Amherst, MA 01003, USA 5Institute for
Energy and Environmental Flows, University of Cambridge, Cambridge CB3 0EZ, UK

###### Abstract

Relatively strongly stratified turbulent flows tend to self-organise into a
‘layered anisotropic stratified turbulence’ (LAST) regime, characterised by
relatively deep and well-mixed density ‘layers’ separated by relatively thin
‘interfaces’ of enhanced density gradient. Understanding the associated mixing
dynamics is a central problem in geophysical fluid dynamics. It is challenging
to study ‘LAST’ mixing, as it is associated with Reynolds numbers
$Re:=UL/\nu\gg 1$ and Froude numbers $Fr:=(2\pi U)/(LN)\ll 1$, ($U$ and $L$
being characteristic velocity and length scales, $\nu$ being the kinematic
viscosity and $N$ the buoyancy frequency). Since a sufficiently large dynamic
range (largely) unaffected by stratification and viscosity is required, it is
also necessary for the buoyancy Reynolds number $Re_{b}:=\epsilon/(\nu
N^{2})\gg 1$ where $\epsilon$ is the (appropriately volume-averaged) turbulent
kinetic energy dissipation rate. This requirement is exacerbated for
oceanically relevant flows, as the Prandtl number
$Pr:=\nu/\kappa=\mathcal{O}(10)$ in thermally-stratified water (where $\kappa$
is the thermal diffusivity), thus leading (potentially) to even finer density
field structures. We report here on four forced fully resolved direct
numerical simulations of stratified turbulence at various Froude ($Fr=0.5,2$)
and Prandtl numbers ($Pr=1,7$) forced so that $Re_{b}=50$, with resolutions up
to $30240\times 30240\times 3780$. We find that, as $Pr$ increases, emergent
‘interfaces’ become finer and their contribution to bulk mixing
characteristics decreases at the expense of the small-scale density structures
populating the well-mixed ‘layers’. However, extreme mixing events (as
quantified by significantly elevated local destruction rates of buoyancy
variance $\chi_{0}$) are always preferentially found in the (statically
stable) interfaces, irrespective of the value of $Pr$.

##  1 Introduction

In density-stratified flows, small-scale turbulence is known to enhance the
rate at which density gradients are irreversibly smoothed by diffusive
processes (i.e. mixed) and quantifying this rate in geophysical, environmental
and industrial settings is of crucial importance. For instance, turbulent
mixing is known to have a leading order impact on global oceanic circulations
(Wunsch & Ferrari, 2004). Stratified turbulent flows are characterised by a
variety of length scales associated with the structure of the turbulent
velocity and density fields. Velocity fluctuations are smoothed by molecular
viscosity and dissipated into heat below the Kolmogorov scale
$L_{K}:=(\nu^{3}/\epsilon)^{1/4}$ (where $\nu$ is the molecular viscosity and
$\epsilon$ is the appropriately volume-averaged dissipation rate of turbulent
kinetic energy). Similarly, density fluctuations are dissipated by molecular
diffusion below the Batchelor scale
$L_{\rho}:=(\nu\kappa^{2}/\epsilon)^{1/4}$, $\kappa$ being the molecular
diffusivity. The (molecular) Prandtl number $Pr:=\nu/\kappa$, that effectively
quantifies the relative strength of viscous dissipation to diffusion, is the
ratio of the square of these two length scales. Unlike the atmosphere, in
which $Pr\simeq 0.7$ ($L_{K}\approx L_{\rho}$), in the ocean
$Pr\sim\mathcal{O}(10)$ when thermally stratified and the equivalent Schmidt
number $Sc:=\nu/D\sim\mathcal{O}(1000)$ ($D$ being the salt diffusivity) when
dominantly salt-stratified.

For turbulent flows with $Pr\gg 1$, the inevitable scale separation with
$L_{\rho}\ll L_{K}$ is such that fully resolved numerical simulations are
highly challenging. However, there is increasing evidence indicating that the
value of $Pr$ has (perhaps unsurprisingly) a leading order impact on
irreversible scalar-mixing properties of stratified turbulent flows. For
instance, the effects of $Pr$ variations on the properties of secondary
instabilities arising from the breakdown of Kelvin-Helmholtz billows were
reported by Mashayek & Peltier (2011) and Salehipour et al. (2015). Similarly,
Zhou et al. (2017) studied the influence of $Pr$ variations on fully developed
turbulence in stratified plane Couette flows using direct numerical
simulations (DNS) and reported significant effects on density and momentum
fluxes. Also using DNS, Legaspi & Waite (2020) analysed the effects of $Pr$
variations on homogeneous forced stratified turbulence and showed, at least
for the range of parameters they considered, that the $Pr=1$ simulations are
able to describe $Pr>1$ dynamics at large scales and also that kinetic energy
spectra (which importantly do not directly contain information from the
density field), are largely unaffected by variations in $Pr$. However, they
did observe $Pr$ effects at small scales below $L_{K}$ and in the density flux
spectra. Therefore, here we focus on the influence of $Pr$ variations on
small-scale structures of the density field and their associated effects on
the mixing properties of forced stratified turbulent flows.

Another important length scale associated with stratified turbulent flows is
the buoyancy length scale $L_{B}:=2\pi U/N_{0}$ where $U$ is a characteristic
velocity of the flow and $N_{0}$ is a characteristic (background) value of the
buoyancy frequency. Scaled using a characteristic horizontal length scale
$L_{h}$, this parameter then defines a horizontal Froude number
$Fr_{h}:=L_{B}/L_{h}=2\pi U/N_{0}L_{h}$ that quantifies the relative strength
of the stratification of the flow. For small horizontal Froude numbers
$Fr_{h}$ (i.e. relatively strongly stratified flows), the density field is
known to self-organise into relatively well-mixed ‘layers’ (whose size scales
as $L_{B}$; the buoyancy scale can therefore be thought of as the largest
energetically possible overturn) separated by relatively thin ‘interfaces’ of
enhanced density gradient (Billant & Chomaz, 2001; Waite, 2011). Such density
‘staircase’ structure has important implications for irreversible scalar
mixing in density-stratified turbulent flows. Couchman et al. (2023) showed
that (in the $Pr=1$ case) while static instabilities are most prevalent within
the well-mixed layers (that are hence characterised by relatively high values
of $\epsilon$), much of the scalar mixing, as described by the (appropriately
volume-averaged) buoyancy variance destruction rate $\chi$, is located in the
relatively strongly stratified interfaces, a phenomenon that is not apparent
when considering $\epsilon$ only. Hence, we focus here on the contribution of
the different structures described above to the overall mixing properties of
the flow, as described by $\chi$.

Fundamentally, we aim to understand how the density field’s small-scale
organisation shapes the bulk mixing properties of forced stratified turbulent
flows at different Prandtl numbers. To this end, we analyse fully resolved DNS
data of forced stratified turbulent flows at $Pr=1$ and $Pr=7$, each with
$Fr=0.5$ and $Fr=2$, for values of $Re$ such that the emergent $Re_{b}=50$.
The rest of this paper is organised as follows. We describe the DNS data in
section 2 and then present the methodology used to extract distinct classes of
density field structures in section 3. In section 4, we apply this methodology
to the DNS data to segment the density field into weakly and strongly
stratified interfaces, relatively small-scale ‘lamella-like’ structures and
larger scale density inversions and analyse the associated mixing properties
of the extracted density structures. Brief conclusions are drawn in section 5.

##  2 Summary of the DNS datasets

Name | $Pr$ | $Fr$ | $Re$ | $Re_{\lambda}$ | Domain size | Grid size | $L_{K}$ | $L_{T}$ | $L_{B}$ | $\epsilon$ | $\chi$
---|---|---|---|---|---|---|---|---|---|---|---
P1F200 | 1 | 2 | 1717 | 184 | $2\pi\times 2\pi\times\pi$ | 4096 $\times$ 4096 $\times$ 2048 | 0.0037 | 0.099 | 1.99 | 1.05 | 0.54
P1F050 | 1 | 0.5 | 24096 | 826 | $2\pi\times 2\pi\times\pi/4$ | 11264 $\times$ 11264 $\times$ 1408 | 0.00052 | 0.030 | 0.52 | 0.95 | 0.38
P7F200 | 7 | 2 | 1717 | 178 | $2\pi\times 2\pi\times\pi$ | 6144 $\times$ 6144 $\times$ 3072 | 0.0037 | 0.096 | 1.95 | 1.09 | 0.44
P7F050 | 7 | 0.5 | 24096 | 655 | $2\pi\times 2\pi\times\pi/4$ | 30240 $\times$ 30240 $\times$ 3780 | 0.00048 | 0.024 | 0.51 | 1.36 | 0.43
Table 1: Summary of the DNS data. $L_{K}$, $L_{T}$ and $L_{B}$ denote the
Kolmogorov, Taylor and buoyancy scales.

We consider statistically steady, forced, fully-resolved DNS of stratified
turbulence from the simulation campaign originally reported by Almalkie & de
Bruyn Kops (2012). The non-hydrostatic dimensionless Navier-Stokes equations
under the Boussinesq approximation,

$\begin{cases}\partial_{t}\boldsymbol{u}+\boldsymbol{u}\bcdot\bnabla\boldsymbol{u}=-\bnabla
p+\frac{1}{Re}\nabla^{2}\boldsymbol{u}-\left(\frac{2\pi}{Fr}\right)^{2}\rho\hat{\boldsymbol{z}}+\mathcal{F},\quad\bnabla\bcdot\boldsymbol{u}=0,\\\
\partial_{t}\rho+\boldsymbol{u}\bcdot\bnabla\rho=\frac{1}{PrRe}\nabla^{2}\rho,\end{cases}$
(1)

are numerically integrated using a pseudospectral code in a triply periodic
domain, where $\hat{\boldsymbol{z}}$ is the (upward) vertical unit vector. The
dimensionless parameters are: the above-defined Prandtl number $Pr$; the
Froude number $Fr:=2\pi U/(N_{0}L)$ (where $U$ and $L$ are characteristic
velocity and length scales associated with the forcing and $N_{0}$ is the
background buoyancy frequency); and the Reynolds number $Re:=UL/\nu$. The
forcing term $\mathcal{F}$ corresponds to the deterministic ‘Rf’ scheme
described in Rao & de Bruyn Kops (2011) and is designed to match a target low
wavenumber kinetic energy spectrum at steady-state. This ensures a buoyancy
Reynolds number $Re_{b}:=\epsilon/(\nu N_{0}^{2})\simeq 50$ for these
simulations, where $\epsilon$ is the volume-average of the point-wise local
dissipation rate of turbulent kinetic energy $\epsilon_{0}$, defined in terms
of the symmetric part of the strain-rate tensor $s_{ij}$ as:

$\epsilon_{0}:=\frac{2}{Re}s_{ij}s_{ij};\
s_{ij}:=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial
u_{j}}{\partial x_{i}}\right).$ (2)

The density field is a superposition of a background linear density field
$\overline{\rho}(z)$, characterised by a reference density $\rho_{0}$ and
reference density gradient $-N_{0}^{2}\rho_{0}/g$, and a perturbation
$\rho^{\prime}(\boldsymbol{x},t)$ so that, in dimensional form,

$\rho(\boldsymbol{x},t)=\overline{\rho}(z)+\rho^{\prime}(\boldsymbol{x},t):=\rho_{0}(1-N_{0}^{2}z/g)+\rho^{\prime}(\boldsymbol{x},t).$
(3)

We can use $\rho^{\prime}$ to compute the (dimensional) local destruction rate
of buoyancy variance:

$\chi_{0}:=\frac{g^{2}\kappa}{\rho_{0}^{2}N_{0}^{2}}\bnabla\rho^{\prime}\bcdot\bnabla\rho^{\prime}.$
(4)

This quantity, appropriately scaled in this way by $-\rho_{0}N_{0}^{2}/g$
(i.e. the density gradient against which the turbulence acts), corresponds to
the local destruction rate of available potential energy (density) and can
therefore be used as a proxy of local irreversible mixing (Howland et al.,
2021) which can still be calculated pointwise-locally in the flow domain.
Henceforth, we consider dimensionless quantities (the prefactor being
$1/(PrRe)$ in our system) and denote by $\chi$ the volume-average of local
$\chi_{0}$ across the entire computational domain.

Here, we consider two values of the Prandtl number ($Pr=\left\\{1,7\right\\}$)
as well as two values of the Froude number ($Fr=\left\\{0.5,2\right\\}$). The
simulation parameters are summarised in table 1. We choose the grid spacing
$\Delta\lesssim L_{\rho}$, and a vertical domain dimension of approximately
$1.5L_{B}$. We consider a single snapshot in time (at statistical steady-
state) of the various simulations.

For the relatively weakly stratified $Fr=2$ simulations, a patch of elevated
turbulence, described by relatively high values of the local dissipation rate
$\epsilon_{0}$, is generated by a relatively large-scale vertically aligned
vortex such as the one reported in (Couchman et al., 2023) for the P1F200
dataset. Outside the vortex, the density field self-organises into relatively
well-mixed layers separated by interfaces characterised by enhanced density
gradients (see figures 1 and 2). In the relatively strongly stratified
$Fr=0.5$ case, the vortex does not appear. Couchman et al. (2023) showed (at
least for the P1F200 simulation) that the interfaces account for relatively
large values of $\chi_{0}$ and low values of $\epsilon_{0}$ and hence
‘extreme’ values of the (local) flux coefficient
$\Gamma_{0}:=\chi_{0}/\epsilon_{0}$, whereas the well-mixed layers were
potentially more turbulent (relatively large values of $\epsilon_{0}$) but
characterised by relatively low values of $\chi_{0}$ (owing to relatively low
values of the density gradient) and hence relatively low values of
$\Gamma_{0}$. Here we aim to understand whether such a picture holds for flows
with different values of $Fr$ and more importantly of $Pr$, and understand how
different structures in the density field (defined more precisely in the next
section) influence the mixing within the flow, quantified by values of (local)
$\chi_{0}$.

Figure 1: Flow segmentation methodology. For the P1F200 simulation, the
density field, represented here by its vertical gradient as shown in panel
$(a)$, is vertically sorted into a statically stable field, whose vertical
gradient is denoted $N^{2}_{\ast}$ (up to a multiplicative constant), as shown
in panel $(b)$. The sorting algorithm highlights the stable interfaces of the
density field, i.e. the points of the density field unaltered by the sorting
procedure. The strongly stratified interfaces (SINT, in blue with
$N^{2}_{\ast}>1$) are then extracted in panel $(c)$. The entire segmented
field is shown in panel $(d)$, including the weakly stratified interfaces with
$N^{2}_{\ast}\leq 1$ (WINT, in purple), and the relatively well-mixed regions
between interfaces, further sub-divided into small-scale ‘lamella’ structures
(LAM, in orange) and larger-scale density inversions (INV, in red) using the
procedure described in section 3 and shown schematically in panel $(e)$, based
on the Taylor microscale $L_{T}$.

##  3 Density field segmentation methodology

Following Couchman et al. (2023), we are interested in the contribution to
mixing of the emergent stably stratified interfaces appearing in the density
field and how this contribution varies as a function of $Pr$ and $Fr$. As $Pr$
increases, we expect finer structures to arise in the intervening well-mixed
regions and we also are interested in understanding how these structures shape
the mixing properties of the flow. Therefore, prior to any analysis, the flow
is segmented into four different categories based on the local (vertical)
density gradient $\partial_{z}\rho$ as well as the neighbouring structure of
the flow. We apply the following methodology, as summarized in figure 1:

  * •

Stably stratified interfaces: We sort vertical density profiles by values of
the density $\rho$ in order to create statically stable density profiles
$\rho^{\ast}(\boldsymbol{x})$. We identify points of the dataset whose values
of the density remain unaltered by the sorting procedure, i.e. points that lie
in statically stable regions of the dataset. These points have relatively high
values of the background buoyancy gradient
$N_{\ast}^{2}:=-g/(\rho_{0}N_{0}^{2})\partial_{z}\rho^{\ast}$ (see figure
1(_b-c_)) and correspond to stably stratified interfaces lying between
relatively well-mixed regions, such as those reported by Couchman et al.
(2023) for the P1F200 dataset. We thus refer to such points as belonging to
the INT (‘interface’) cluster, which is further subdivided into relatively
strongly stratified (SINT) and relatively weakly stratified (WINT) interfaces,
for $N_{\ast}^{2}>1$ and $N_{\ast}^{2}\leq 1$ respectively.

  * •

Relatively well-mixed layers: We subdivide points that are altered by the
sorting procedure (i.e. are not within an interface INT) into two categories
depending on the local density gradient $\partial_{z}\rho$ and neighbouring
structure of the (unsorted) density field $\rho$, in particular relative to
$L_{T}$, the Taylor microscale of the flow. $L_{T}$ describes the scale at
which viscosity starts to affect the development of turbulent eddies
significantly, defined (dimensionally) as
$L_{T}:=\sqrt{10\nu{\mathcal{K}}/\epsilon}$, where ${\mathcal{K}}$ is the
volume averaged turbulent kinetic energy. As illustrated in figure 1(e), by
considering a point in the dataset at vertical coordinate $z_{c}$ and density
$\rho_{c}$ satisfying $\partial_{z}\rho(z_{c})>0$, we define $z_{+}$ as the
closest point moving upwards for which $\rho(z_{+})=\rho_{c}$. As equality is
difficult to ensure numerically, we define

$z_{+}:=\text{min}\\{z>z_{c}\;|\;\rho(z)\leq\rho(z_{c})\\}.$ (5)

Similarly, if $\partial_{z}\rho(z_{c})<0$, we define
$z_{-}:=\text{min}\\{z<z_{c}\;|\;\rho(z)\geq\rho(z_{c})\\}$. If $\lvert
z_{c}-z_{+}\rvert\leq L_{T}$ (or $\lvert z_{c}-z_{-}\rvert\leq L_{T}$
depending on the local value of the density gradient), we have identified a
relatively small scale structure in the density field, such as a blob or
stretched ‘lamella’ as discussed in detail by Villermaux (2019), which is
strongly affected by viscosity. Therefore we classify these structures as
belonging to the LAM (‘lamella’) cluster. Conversely, if $\lvert
z_{c}-z_{+}\rvert>L_{T}$ (or $\lvert z_{c}-z_{-}\rvert>L_{T}$), we have
identified a point belonging to a relatively large scale density inversion
largely unaffected by viscosity, and so we classify it as being in the INV
(‘inversion’) cluster. This segmentation methodology, summarised in figure 1,
is applied to the four datasets considered in this work, as shown in figure 2.

Figure 2: Normalized vertical density gradient
$\partial_{z}\rho^{\prime}/|\partial_{z}\overline{\rho}|$ (panels $(a,c,e,g)$)
and associated segmented fields (panels $(b,d,f,h)$), for simulation P1F200
($a,b$), P1F050 ($c,d$), P7F200 ($e,f$) and P7F050 ($g,h$). The strongly
stratified interfaces (SINT) are depicted in blue, the small-scale structures
(LAM) in orange, the larger-scale density inversions (INV) in red and the
weakly stably stratified regions (WINT) in purple.

##  4 Cluster properties

###  4.1 Contributions to $\chi$

We consider the relative importance of each cluster defined in section 3 in
terms of their contribution to the local and bulk mixing properties of the
flow, as described by local $\chi_{0}$ and volume-averaged $\chi$. We
particularly focus on the role of ‘extreme’ mixing events and therefore sort
the data by decreasing values of $\chi_{0}$, defining a sorted vector
$\boldsymbol{\chi}_{0}^{*}=(\chi_{0}^{0*},\cdots,\chi_{0}^{M*})$ where $M$ is
the number of points in the dataset and $\chi_{0}^{0*}$ and $\chi_{0}^{M*}$
are the greatest and lowest values of $\chi_{0}$ found in the dataset,
respectively. This vector can be used to construct the normalized cumulative
contribution to $\chi$ as

$\forall
n\in\\{1,\cdots,M\\},\;\chi^{c}(n):=\frac{1}{\chi}\sum_{i=1}^{n}\chi_{0}^{i*}.$
(6)

The cumulative contribution $\chi^{c}$ is plotted in figure 3 for each dataset
(solid black line), highlighting the fact that the dominant contribution to
$\chi$ is generated by a relatively small set of localised ‘extreme’ events,
regardless of $Fr$ and $Pr$. More specifically, for the four sets of
parameters considered in this work, approximately $80\%$ of the contribution
to $\chi$ is contained within the first $10\%$ of the points sorted by
decreasing values of $\chi_{0}$ (and hence $10\%$ of the box volume). Similar
statistics are observed in oceanographic data (Couchman et al., 2021).

The total contribution of each cluster to bulk $\chi$ is shown in figure
3(_e_). We also assign the data points sorted by values of $\chi_{0}$ into
$n=20$ equal volume bins, in order to calculate a relative contribution of
points in each cluster to the bin-volume-averaged
$\langle\chi\rangle_{\text{bin}}:=\sum_{\text{bin}}\chi_{0}/(V/n)$, where $V$
is the total number of points in the domain, as illustrated by colored shading
in figures 3(_a-d_). Focusing on the $Pr=1$ case, the contribution (figure
3(_e_)) from the strongly stratified interfaces (SINT cluster) to bulk $\chi$
is roughly $25-35\%$, whereas the contribution from small-scale structures
(LAM) reaches $\sim 60\%$. As $Fr$ decreases and stratification strengthens,
the total contribution from large-scale inversions (INV) shrinks from $\sim
10\%$ for $Fr=2$ to less than $3\%$ for $Fr=0.5$, possibly due to the
suppression of large-scale overturnings by relatively strong stratification.
We observe similar trends for the relative contributions in each bin. In
increasing the Prandtl number from $Pr=1$ to $7$, strong interfaces become
(perhaps unsurprisingly) finer (figure 2(_f,h_)) and their total contribution
shrinks to about $10\%$, whereas the total contribution from small-scale
structures in ‘lamella’ (LAM) increases to about $80-90\%$. The total
contribution from large-scale inversions does not exceed $10\%$. Again, we
observe similar trends in the relative contributions as $\chi_{0}$ decreases.

Figure 3: $(a$-$d)$ Normalized cumulative contribution to $\chi$ (black line;
see equation (6)) for each simulation. Data points are assigned to 20 equal
volume bins, sorted by decreasing $\chi_{0}$ and clustered using the method
presented in section 3. For each bin, we compute the relative contribution of
each cluster to $\langle\chi\rangle_{\text{bin}}$, as shown by the heights of
the colored regions. For each bin, we compute the fraction of points within
each cluster to the bin’s total number of points (red circled dots), along
with the (arithmetic) mean value of $\Gamma_{0}:=\chi_{0}/\epsilon_{0}$ for
each cluster (dashed lines). (_e_) Total contributions to $\chi$ from each
cluster for the four simulations.  Figure 4: Probability density functions
(PDF) for $\log_{10}(\chi_{0})$ (panel (_a_), (_e_)),
$\log_{10}(\epsilon_{0})$ (panel (_b_), (_f_)) and $\log_{10}(\Gamma_{0})$
(panel (_c_), (_g_)) for the different clusters defined in section 3 and for
the four simulations. The statistical mean of each field (without the
logarithm) is represented by a colour-coded circle ($Pr=1$) or triangle
($Pr=7$). The green dotted vertical line corresponds to the canonical value
$\Gamma=0.2$.

###  4.2 Statistics of $\Gamma$

For each bin in figure 3, we compute a mean value of the (local) flux
coefficient $\Gamma_{0}:=\chi_{0}/\epsilon_{0}$ associated with each cluster
(see dashed lines), although such mean values of ratios of local quantities
should always be treated with caution. The leftmost bins (corresponding to the
most ‘extreme’ values of $\chi_{0}$) have the largest values of $\Gamma_{0}$
(of order unity, significantly above the canonical value $\Gamma=0.2$)
suggesting strongly that the extreme mixing events (large $\chi_{0}$) are not
necessarily correlated with extreme turbulent events (large $\epsilon_{0}$).
Among the extreme events in $\chi_{0}$, those in the strongly stratified
interfaces (SINT) cluster correspond to the largest mean values of
$\Gamma_{0}$ (as suggested by Couchman et al. (2023) for the $Pr=1$, $Fr=2$
case), independently of variations in $Pr$ and $Fr$. This point is further
emphasized in figure 4 where we plot the probability density functions (PDFs)
of $\log_{10}(\chi_{0})$, $\log_{10}(\epsilon_{0})$ and
$\log_{10}(\Gamma_{0})$ for the different clusters and simulations. We also
present the statistical mean of $\chi_{0}$, $\epsilon_{0}$ and $\Gamma_{0}$.
For each simulation, the statistical mean of the flux coefficient $\Gamma_{0}$
is maximized for the strongly stratified interfaces (almost twice as large as
the statistical mean for the entire dataset), because of relatively low values
of $\epsilon_{0}$ and large values of $\chi_{0}$. As the Prandtl number
increases, the statistical mean of $\Gamma_{0}$ decreases (a result also
observed by Salehipour et al. (2015)) for both the weakly and strongly
stratified cases, for the following (different) reasons.

In the weakly stratified simulations with $Fr=2$ (figures 4(_a-c_)), density
effectively acts as a passive scalar (as can be seen from the weakly
stratified scaling of Riley et al. (1981), for instance) and hence the
statistics of $\epsilon_{0}$ do not depend on $Pr$, as can be seen in figure
4(_b_). However, as $Pr$ increases, the left tails of the $\chi_{0}$
distributions become more significant, as might be expected for the mixing of
an effectively passive scalar, as $\chi_{0}$ is multiplied by $(PrRe)^{-1}$
where $Re$ is fixed. We note that a widening of the tails of
$\partial_{z}\rho^{\prime}$ (not shown here, see for instance Riley et al.
(2023)) effectively rebuilds the right tail of $\chi_{0}$, explaining why the
PDF of $\chi_{0}$ is not just shifted toward lower values. A more in-depth
discussion of this phenomenon is provided by Bragg & de Bruyn Kops (2023). As
a result of the statistics of $\epsilon_{0}$ remaining roughly independent of
$Pr$ but those of $\chi_{0}$ decreasing with increasing $Pr$, the statistics
of $\Gamma_{0}$ decrease as $Pr$ increases for the weakly stratified ($Fr=2$)
case.

Conversely, in the strongly stratified simulations with $Fr=0.5$, buoyancy
acts ‘actively’ on the momentum field, as demonstrated for instance by the
strongly stratified scaling analysis of Billant & Chomaz (2001). Moreover, as
$Pr$ increases, the volume contribution of the small-scale structures (LAM
cluster) increases at the expense of the strongly stratified interfaces (SINT)
(see figure 3). LAM structures, which locally are only weakly affected by
stratification because they populate the relatively well-mixed regions of the
flow, are viscously affected (as their vertical extent is, by definition,
smaller than $L_{T}$) but still significantly disordered. As a result, their
local static instability inevitably encourages increased viscous dissipation,
and so their enhanced prevalence actually shifts the statistics of
$\epsilon_{0}$ towards higher values for the $Pr=7$ flows as compared to
$Pr=1$ (yellow curves, figure 4(_f_)). In this sense, enhanced stratification
at higher $Pr$ actually makes the flow ‘more turbulent’, reminiscent of the
prediction by Pearson & Linden (1983) of relatively long-lived ‘approximately
horizontal striations’ in high-$Pr$ decaying stratified turbulence. Also,
though this is a second order effect, for the strongly stratified simulations
the left tail of the $\chi_{0}$ distribution once again becomes somewhat more
significant as $Pr$ increases. We note that while the forcing scheme endeavors
to maintain a constant $Re_{b}$, bulk $\epsilon$ is slightly mismatched
between the two $Pr$ simulations at $Fr=0.5$ (see Table 1). Therefore, while
increased $Pr$ dramatically changes the relative prevalence of LAM versus SINT
structures, part of the observed increase in local $\epsilon_{0}$ statistics
at $Pr=7,Fr=0.5$ is due to the slight mismatch in targeted bulk $\epsilon$
rather than being a purely $Pr$ effect. Importantly, it is apparent in figure
4(_e_) that the total amount of irreversible mixing for the flow with $Pr=7$
actually slightly _increases_ compared to $Pr=1$ in the more strongly
stratified $Fr=0.5$ simulations, essentially because the $Pr=7$ flow is more
vigorous, despite the fact that $\Gamma_{0}$ appears to decrease with
increasing $Pr$. Therefore, consideration of $\Gamma_{0}$, or even
$\Gamma:=\chi/\epsilon$ (constructed from volume-averaged dissipation rates),
in isolation should be treated with caution.

##  5 Discussion

We have analysed the influence of the Prandtl number $Pr$ and Froude number
$Fr$ on density structures in forced turbulent stratified flows and their
contribution to irreversible scalar mixing properties. Using fully resolved
DNS data and a flow segmentation algorithm based on the local value of the
(vertical) density gradient and the local (vertical) structure of the density
field, we have extracted distinct regions of the turbulent density field –
interfaces (both relatively strong and relatively weak) separating well-mixed
density layers made up of both small-scale lamellar structures and larger
scale density inversions – and analysed their contribution to the bulk value
of the destruction rate of buoyancy variance $\chi$ and to the statistics of
the (locally evaluated) flux coefficient $\Gamma_{0}$.

As $Pr$ increases, the strongly stratified density interfaces become finer and
their contribution to average values of $\chi$ decreases at the expense of the
small-scale lamellar structures in the relatively well-mixed density layers.
However, similarly to the flow with $Pr=1$ (Couchman et al., 2023), these
structures are ‘quiescent’, yet mixing hotspots. The points in these
structures associated with the largest values of (local) $\chi_{0}$ are
associated with extreme values of $\Gamma_{0}$ and hence with relatively low
values of $\epsilon_{0}$. More generally, the flux coefficient associated with
strongly stratified structures is (in an averaged sense) essentially twice as
large as its value for the three other segmented regions, regardless of the
values of $Pr$ and $Fr$. All in all, strongly stratified interfaces are
therefore characterised by relatively large values of $\chi_{0}$ and
relatively weak values of $\epsilon_{0}$ compared to the other density
structures considered in this work and are therefore likely to be overlooked
if considering $\epsilon_{0}$ (or indeed the volume average $\epsilon$) alone
as a proxy for significant irreversible mixing.

As is becoming increasingly well-appreciated, a universal value of $\Gamma$ is
not able to capture the complex and inhomogeneous structures of the density
and density gradient turbulent fields. Moreover, our description of the
density fields as well as consideration of the considered mixing regime (i.e.
‘passive’ mixing at higher $Fr$ as compared to ‘active’ mixing at lower $Fr$)
offers a potential explanation for the empirical observation that the ‘mixing
efficiency’ $\Gamma/(1+\Gamma)$ decreases with $Pr$. On the one hand, mixing
in the relatively weakly stratified case can be described as ‘passive’: the
velocity field and hence $\epsilon_{0}$ are largely unaffected by $Pr$ but the
small-scale structure of the density field evolves as $Pr$ increases from $1$
to $7$ in a way that results in a decrease in $\chi_{0}$ and thus the flux
coefficient $\Gamma_{0}$. On the other hand, mixing in the relatively strongly
stratified case can be thought of as ‘active’, with the density field and its
small-scale lamellar structure having a leading order impact on the rate at
which the turbulent kinetic energy is dissipated ($\epsilon_{0}$), ultimately
leading to a decrease in $\Gamma_{0}$ as $Pr$ increases. This shows another
way in which a focus on the flux coefficient can be misleading. Since the
decrease in $\Gamma_{0}$ is actually principally related to an increase in
$\epsilon_{0}$, the total amount of mixing (quantified by $\chi$) actually
_increases_ in the higher $Pr$, strongly stratified flow considered here.
Another key result of our analysis is that the Prandtl number $Pr$ has a
leading order impact (at least compared to the Froude number $Fr$) on the
density field structures and their mixing properties, further emphasizing the
need to take this parameter into account when ‘measuring mixing’, a process
that is inherently a diffusive one (Villermaux, 2019; Caulfield, 2021) and
that affects the density field.

We here considered steady-state forced stratified turbulent flows, and so it
would now be interesting to apply our segmentation and associated mixing
contribution analysis to time-dependent decaying stratified turbulent flows.
For example, what is the fate of the different emerging structures? Can a
separate analysis of the temporal evolution of each structure help us better
understand the mixing history of a stratified turbulent flow as a whole and
unveil the fundamental rules of stratified mixing? These are questions left
for future work.

[Funding]This project received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie Grant
Agreement No. 956457 and used resources of the Oak Ridge Leadership Computing
Facility at the Oak Ridge National Laboratory, supported by the Office of
Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725.
S.deB.K. was supported under U.S. Office of Naval Research Grant number
N00014-19-1-2152. For the purpose of open access, the authors have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising from this submission.

## References

  * Almalkie & de Bruyn Kops (2012) Almalkie, S. & de Bruyn Kops, S. M. 2012 Kinetic energy dynamics in forced, homogeneous, and axisymmetric stably stratified turbulence. J. Turb. (13), N29.
  * Billant & Chomaz (2001) Billant, P. & Chomaz, J.-M. 2001 Self-similarity of strongly stratified inviscid flows. Phys. Fluids 13 (6), 1645–1651.
  * Bragg & de Bruyn Kops (2023) Bragg, A. D. & de Bruyn Kops, S. M. 2023 Understanding the effect of Prandtl number on momentum and scalar mixing rates in neutral and stably stratified flows using gradient field dynamics. arXiv:2308.00518 .
  * Caulfield (2021) Caulfield, C.P. 2021 Layering, instabilities, and mixing in turbulent stratified flows. Annu. Rev. Fluid Mech. 53 (1), 113–145.
  * Couchman et al. (2023) Couchman, M. M. P., de Bruyn Kops, S. M. & Caulfield, C. P. 2023 Mixing across stable density interfaces in forced stratified turbulence. J. Fluid Mech. 961, A20.
  * Couchman et al. (2021) Couchman, M. M. P., Wynne-Cattanach, B., Alford, M. H., Caulfield, C. P., Kerswell, R. R., MacKinnon, J. A. & Voet, G. 2021 Data-driven identification of turbulent oceanic mixing from observational microstructure data. Geophys. Res. Lett. 48 (23), e2021GL094978.
  * Howland et al. (2021) Howland, C. J., Taylor, J. R. & Caulfield, C. P. 2021 Quantifying mixing and available potential energy in vertically periodic simulations of stratified flows. J. Fluid Mech. 914, A12.
  * Legaspi & Waite (2020) Legaspi, J. D. & Waite, M. L. 2020 Prandtl number dependence of stratified turbulence. J. Fluid Mech. 903, A12.
  * Mashayek & Peltier (2011) Mashayek, A. & Peltier, W. R. 2011 Turbulence transition in stratified atmospheric and oceanic shear flows: Reynolds and Prandtl number controls upon the mechanism. Geophys. Res. Lett. 38 (16).
  * Pearson & Linden (1983) Pearson, H. J. & Linden, P. F. 1983 The final stage of decay of turbulence in stably stratified fluid. J. Fluid Mech. 134, 195–203.
  * Rao & de Bruyn Kops (2011) Rao, K. J. & de Bruyn Kops, S. M. 2011 A mathematical framework for forcing turbulence applied to horizontally homogeneous stratified flow. Phys. Fluids 23 (6).
  * Riley et al. (2023) Riley, J. J., Couchman, M. M. P. & de Bruyn Kops, S. M. 2023 The effect of Prandtl number on decaying stratified turbulence. J. Turb. 0 (0), 1–19.
  * Riley et al. (1981) Riley, J. J., Metcalfe, R. W. & Weissman, M. A. 1981 Direct numerical simulations of homogeneous turbulence in density-stratified fluids. In AIP Conference Proceedings, , vol. 76, pp. 79–112.
  * Salehipour et al. (2015) Salehipour, H., Peltier, W. R. & Mashayek, A. 2015 Turbulent diapycnal mixing in stratified shear flows: the influence of Prandtl number on mixing efficiency and transition at high Reynolds number. J. Fluid Mech. 773, 178–223.
  * Villermaux (2019) Villermaux, E. 2019 Mixing versus stirring. Annu. Rev. Fluid Mech. 51, 245–273.
  * Waite (2011) Waite, M. L. 2011 Stratified turbulence at the buoyancy scale. Phys. Fluids 23 (6).
  * Wunsch & Ferrari (2004) Wunsch, C. & Ferrari, R. 2004 Vertical mixing, energy, and the general circulation of the oceans. Annu. Rev. Fluid Mech. 36, 281–314.
  * Zhou et al. (2017) Zhou, Q., Taylor, J. R. & Caulfield, C. P. 2017 Self-similar mixing in stratified plane Couette flow for varying Prandtl number. J. Fluid Mech. 820, 86–120.


# Universal Relations for the Increase in the Mass and Radius of a Rotating
Neutron Star

Andreas Konstantinou Department of Physics
University of Alberta
Edmonton AB, T6G 2E1, Canada  Sharon M. Morsink Department of Physics
University of Alberta
Edmonton AB, T6G 2E1, Canada

(Accepted by ApJ, June 21, 2022)

###### Abstract

Rotation causes an increase in a neutron star’s mass and equatorial radius.
The mass and radius depend sensitively on the unknown equation of state (EOS)
of cold, dense matter. However, the increases in mass and radius due to
rotation are almost independent of the EOS. The EOS independence leads to the
idea of neutron star universality. In this paper, we compute sequences of
rotating neutron stars with constant central density. We use a collection of
randomly generated EOS to construct simple correction factors to the mass and
radius computed from the equations of hydrostatic equilibrium for non-rotating
neutron stars. The correction factors depend only on the non-rotating star’s
mass and radius and are almost independent of the EOS. This makes it
computationally inexpensive to include observations of rotating neutron stars
in EOS inference codes. We also construct a mapping from the measured mass and
radius of a rotating neutron star to a corresponding non-rotating star. The
mapping makes it possible to construct a zero-spin mass-radius curve if the
masses and radii of many neutron stars with different spins are measured.

We show that the changes in polar and equatorial radii are symmetric, in that
the polar radius shrinks at the same rate that the equatorial radius grows.
This symmetry is related to the observation that the equatorial compactness
(the ratio of mass to radius) is almost constant on one of the constant-
density sequences.

neutron stars, relativistic stars, pulsars, stellar rotation, compact objects,
millisecond pulsars

##  1 Introduction

The cores of neutron stars provide a unique combination of low temperature and
high density that can provide insights about the equation of state (EOS) of
supranuclear density matter (Watts et al., 2016). While the EOS of the matter
in a neutron star is still undetermined, each proposed EOS is mapped to a
unique mass-radius curve by the relativistic equations of stellar structure.
As a result, the measurement of the masses and radii of many neutron stars
can, in principle, be used to determine the neutron star’s EOS (Lindblom,
1992; Lindblom & Indik, 2012, 2014). This is one of the main motivations for
the Neutron Star Interior Composition Explorer (NICER) X-ray telescope’s
observations of rotation-powered pulsars (Gendreau et al., 2012).

The use of mass and radius measurements of neutron stars to constrain the cold
dense EOS is typically done through the framework of statistical inference
(see, e.g. Miller et al. (2020)). A set of $N$ neutron stars, each with known
spin are observed, and a likelihood function for each neutron star’s mass and
equatorial radius is found. These likelihood functions could come from pulse-
profile modelling of x-ray timing data, as is done with NICER observations
(Miller et al., 2019; Riley et al., 2019; Miller et al., 2021; Riley et al.,
2021), for example. Alternatively, spectroscopic observations of neutron stars
in quiescent low-mass x-ray binaries or thermonuclear bursts on the surfaces
of accreting neutron star could yield likelihood functions (Özel & Freire,
2016). Given a set of parameters describing an EOS, the relativistic stellar
structure equations predict values of mass and equatorial radius for each
neutron star’s particular value of spin. Each neutron star is assumed to be
described by the same set of EOS parameters, but each will have a different
density at its centre. This allows the calculation of the likelihood that a
particular set of EOS parameters correctly describes the set of $N$ neutron
star observations.

Although the procedure described in the previous paragraph is well-defined, if
many sets of EOS parameters are to be tested, the computation of rapidly
rotating neutron star models would make this time-consuming and
computationally expensive. In practice, the neutron stars observed by NICER
are rotating slowly, with spin frequencies near $200$ Hz. As a result, the
computation of rotating models is not necessary, since the masses and radii do
not change drastically at slow rotation rates. EOS inference for the NICER
neutron stars makes use of the spherically symmetric TOV equations (Miller et
al., 2019; Raaijmakers et al., 2019; Miller et al., 2021; Raaijmakers et al.,
2021) to map between EOS parameters and the mass-radius values. While the use
of the TOV equations does not introduce any bias in the NICER results, future
observations, either with better precision, or of more rapidly rotating
neutron stars could lead to biased results. The goal of this paper is to
evaluate the magnitude of the inaccuracies that result from the TOV
approximation in EOS inference, and to provide a simple, computationally
inexpensive method for including rotation in EOS inference codes.

The concept of universality is useful when considering the effects of rotation
on neutron stars. Black holes are truly universal: specifying only the mass
and dimensionless spin parameter leads to a complete description of the black
hole’s gravitational field and properties. Neutron stars, at first glance, do
not seem to have universal properties. Each proposed EOS predicts a different
curve of mass as a function of radius, $M(R)$. However, other properties, such
as the moment of inertia (Ravenhall & Pethick, 1994) and the gravitational
binding energy (Lattimer & Prakash, 2001), have an approximately universal
dependence on the compactness, $C=M/R$, independent of the EOS. When the
neutron star rotates at angular velocity, $\Omega$, many properties depend
approximately on only the dimensionless angular velocity, defined by
$\bar{\Omega}^{2}=\Omega^{2}GM/R_{e}^{3}$, where $R_{e}$ is the equatorial
radius, along with the compactness. For instance, approximately universal
relations depending on $\bar{\Omega}$ and $C_{e}=M/R_{e}$ for the oblate shape
of the neutron star (Morsink et al., 2007; Silva et al., 2021) and the
effective acceleration due to gravity over the neutron star’s surface (AlGendy
& Morsink, 2014) have been used to determine the equatorial radii of neutron
stars by NICER (Bogdanov et al., 2019). An approximately universal relation
between the maximum mass rotating neutron star and the corresponding maximum
mass for a non-rotating star (Lasota et al., 1996; Breu & Rezzolla, 2016) has
proven to be useful for placing limits on the properties of the remnants of
neutron star mergers such as GW170817 (Rezzolla et al., 2018). The oscillation
modes of rotating neutron stars also exhibit universal properties (Krüger &
Kokkotas, 2020), which could be used to constrain the EOS (Völkel & Krüger,
2022). An even stronger universality is the “I-Love-Q” relation between the
moments of inertia, Love numbers, and quadrupole moments of rotating neutron
stars (Yagi & Yunes, 2013). A more detailed review of neutron star
universality is provided by Yagi & Yunes (2017).

In this paper we explore the universality of the change in the mass and radius
of uniformly rotating neutron stars in order to introduce some approximate
relations for the changes that occur when a neutron star spins. These
increases are compared to non-rotating neutron stars with the same central
density. The increases in mass and radius are written as simple expressions,
independent of the EOS, depending only on the mass and radius of the non-
rotating reference star as well as the angular velocity of the rotating star.
This allows the construction of a set of simple correction factors for the
masses and radii of neutron stars calculated using the spherically symmetric
TOV equations. We also construct an inverse mapping from rotating stars to
non-rotating stars to allow the reconstruction of the zero-spin mass-radius
curve from observations of rotating neutron stars.

The structure of this paper is as follows. In Section 2 we introduce our
library of randomly generated EOS, and compute sequences of rapidly rotating
neutron stars connected by a constant value of central density. In Section 3
we explore the dependence of the Kepler frequency on the mass and radius of
the non-rotating member of a sequence, which we exploit to construct universal
formulae for the increases in mass and radius due to rotation. In this section
we also discuss a simple physical model for the universal behaviour seen and
its application to EOS inference using only the solutions to the TOV
equations. In Section 4 we construct a mapping from the measured properties of
a rotating star to the properties of non-rotating neutron stars as a method to
determine the zero-spin mass-radius curve.

##  2 Spin-up sequences

We use the computer program Rapidly Rotating Neutron Star, rns, (Stergioulas &
Friedman, 1995) to solve the relativistic stellar structure equations for an
axisymmetric, rotating neutron star given a zero-temperature EOS. The rns
program is an adaptation of the Green function method (Komatsu et al., 1989;
Cook et al., 1992) for solving the Einstein field equations for axisymmetric
stars, which is based on the Newtonian self-consistent field method (Hachisu,
1986).

###  2.1 Equations of state

We created two separate families of randomly generated EOS. The first family
is the piecewise polytrope (PP) (Read et al., 2009), constructed with the same
range of parameters given by Hebeler et al. (2013). The second family of EOS
is constructed using the speed-of-sound (CS) parametrization introduced by
Greif et al. (2019). For both families we randomly chose values for the
various parameters using flat priors for the ranges specified in the
respective papers. Initially we constructed separate universal relations for
each family of EOS in order to test if there is any bias introduced by the
choice of EOS family. However, the differences were minor, so we present
results in this paper based on all of the randomly generated EOS from both
families. We determined the optimal number of EOS to generate by first
choosing a small number of EOS, computing the full set of rotating neutron
stars, and then fitting the increase in radius caused by rotation. We then
iteratively generated a new EOS and then fit the radius increase again,
checking the convergence of the best-fit coefficients. We found that a set of
32 randomly generated EOS provided a fit to the mass that converged. The CS
method creates a more restricted range of radii, so we limited our choices to
13 EOS from this family and used 19 EOS from the PP family in order to find
empirical relations that are valid for a wider range of possible EOS.

The maximum masses (of non-rotating stars) generated from these 32 EOS range
from $1.97M_{\odot}$ to $2.88M_{\odot}$. The radius of spherical 1.4
$M_{\odot}$ stars computed with the PP EOS range from 10.9 - 14.7 km, while
the CS EOS stars range from 11.2 - 12.4 km.

Other methods for constructing phenomenological parametric EOS exist, such as
the spectral method (Lindblom & Indik, 2014). Legred et al. (2022) have shown
that parametric EOS can introduce correlations between the pressures at
different densities that add artificial restrictions to the EOS that can bias
inference results, and show that non-parametric methods, such as the Gaussian
process method (Landry & Essick, 2019) do not suffer from this type of bias.

It is beyond the scope of this paper to explore all other types of
phenomenological EOS. Instead we test the empirical relations derived in the
following section against a small set of tabulated EOS to look for extreme EOS
models where our results might fail. We included the following hadronic EOS:
the APR EOS (Akmal et al., 1998) includes 3-nucleon interactions and special
relativistic corrections; the two BBB EOS (Baldo et al., 1997) as hadronic
examples that have maximum masses smaller than $2M_{\odot}$; an example of an
EOS that includes hyperons (Table 5.8 of Glendenning (1997)), with a maximum
mass below $1.6M_{\odot}$; the ultra stiff EOS L (Pandharipande & Smith, 1975)
with a pion condensate. We also include three quark/hadron hybrid EOS (Alford
et al., 2005; Kojo et al., 2021), and a bare quark star (Alcock et al., 1986).

###  2.2 Constant central-density spin-up sequences for one EOS

In our application, we compute sequences of stars with the same EOS and
central-density. The first star of the sequence is non-rotating, with a mass
and radius denoted by $M_{*}$ and $R_{*}$ given by the solution of the TOV
equations. The subsequent stars in the same sequence have the same central
density, but increasing eccentricity. The values of mass and equatorial radius
for each subsequent star in the sequence, $M$ and $R_{e}$, are larger than the
corresponding values of $M_{*}$ and $R_{*}$ for the same central density.

As an example we provide a full set of constant central-density spin-up
sequences for one sample EOS, the piecewise polytrope PP0, in Figure 1. The
choice of EOS to illustrate the sequences is not important, since the general
features and shapes of the plots for each EOS are similar, even if the
particular values of mass and radius are different. Each stable stellar model
computed is represented as a point on these plots. Figure 1 (left pane) shows
each spin-up sequence as a vertical series of dots. The lowest mass star in
each sequence corresponds to a solution of the TOV equations. Solutions of the
TOV equations are shown with a bold black curve. sThe star with the largest
mass on this curve is the TOV limit star. Stars with larger mass and same the
central density in a sequence correspond to a star with a larger spin
frequency. The two coloured vertical lines show the minimum central densities
for the two densest polytropic regions for this particular PP.

Each spin-up sequence terminates at the highest mass shown, which is also the
highest frequency possible for a stable star in the sequence. The termination
of the sequence is either due to the Kepler limit (also known as the mass-shed
limit), shown with a bold dashed line, or the onset of the quasi-radial
instability, shown with a thin dashed line. In the case of the maximum mass
non-rotating star, there is only one member of this sequence, since a rotating
star with the same central density as the maximum mass non-rotating star is
unstable to quasi-radial modes. In our calculations we approximate the quasi-
radial instability curve as a straight line in mass versus the logarithm of
the central density as suggested by the results of Cook et al. (1994) (see
their Figure 1, for example). While this is not exact, our results don’t
depend on the exact determination of the quasi-radial instability onset.

This procedure allows us to compute a set of sequences of rotating neutron
stars that include both normal and supramassive neutron stars (Cook et al.,
1992). A solid curve dividing the models into normal and supramassive stars is
shown in Figure 1a, corresponding to stars with the same baryon mass as the
maximum mass non-rotating stellar model. Stars below this curve are normal,
while stars above the curve are supramassive.

The same constant-density spin-up sequences are shown in the right-hand pane
of Figure 1, where the mass and equatorial radius of each stable star are
shown, along with the same curves described for the left-hand pane of Figure
1. Each constant-density sequence is approximately a straight line with
positive slope, so that the star’s equatorial compactness is approximately
constant as it spins faster. We will return to this property later in Section
3.2. Additionally we add a gold curve, representing stars that are rotating at
95% of the Kepler frequency, and gold points that overlap the TOV curve,
representing an inverse map which will be explained in Section 4.

Figure 1: A set of constant central-density spin-up sequences for EOS PP0
showing mass versus central density (left) and mass versus equatorial radius
(right). Each blue point represents the mass, equatorial radius, and central
energy density of a rotating neutron star constructed with this EOS.

Note that the sequences shown in Figure 1 are different from the constant
baryon-mass sequences used in many applications (eg. Cook et al. (1994)) to
model the spin-down of an isolated star that doesn’t lose or accrete matter.
In particular, it may be useful to compare the plots shown in Figure 1 with
Figures 1 and 2 of Cook et al. (1994) (for a different EOS). In this paper, we
are investigating the size of spin corrections to the mass and radius obtained
from the TOV equations, so a comparison between models with the same central
energy density is more useful.

##  3 Universal relations for the change in mass and radius

We now investigate some approximately universal properties of neutron stars as
they spin. Our approach is to consider dimensionless quantities, such as
$M/M_{*}$, the mass of a rotating star divided by the mass of a non-rotating
star with the same central density, and investigate how this ratio depends on
different quantities. One natural choice of independent variable is
$C_{*}=\frac{M_{*}}{R_{*}}\frac{km}{M_{\odot}}$, the non-dimensional
compactness of the non-rotating star in a sequence. (Note that many authors
define compactness as the ratio $GM/Rc^{2}$ which differs from this definition
by a constant factor.) Similarly we could introduce the star’s angular
velocity normalized by $\Omega_{*}=\left(GM_{*}/R_{*}^{3})\right)^{1/2}$ or
alternatively, the dimensionless spin parameter $\chi=cJ/GM^{2}$, where $J$ is
the star’s angular momentum. (We use the following notation: the star’s spin
frequency, as measured at infinity is $\nu$, and the angular velocity,
$\Omega$, is $\Omega=2\pi\nu$.) However, we have found that the behaviour of
quantities such as $M/M_{*}$ for multiple EOS depends more closely on how
close the angular velocity is to the Kepler limit than either the quantity
${\Omega}/\Omega_{*}$or $\chi$. For this reason, we begin our exploration of
universality with the Kepler limit.

###  3.1 The Kepler limit

The Kepler limit occurs when a star spins at the same rate that a test
particle orbits at the star’s equator. We use the symbol $\Omega_{K}$ to
denote a star spinning at the Kepler limit. To investigate the Kepler limit,
we first choose an EOS, and find all of the stellar models that are spinning
at the Kepler limit. For instance, in the case of EOS PP0, this would
correspond to the set of 31 points that terminate the constant-density spin-up
sequences on the heavy dashed curve shown in Figure 1. Each of these spin-up
sequences has different non-rotating reference values for $M_{*}$ and $R_{*}$.
The angular velocity for each of these stars is then normalized by
$\Omega_{*}$ and plotted versus the compactness of the zero spin member of its
sequence. This procedure is then repeated for each of the EOS in our library.
The result is shown in Figure 2.

Figure 2: The Kepler angular velocity, $\Omega_{K}$ of neutron stars as a
function of the compactness of the zero spin star with the same central
density as the rotating star. Stars computed with a PP or CS EOS are shown
with blue or gold dots. The best-fit curve, Equation (1) is shown as a solid
curve.

Each point shown in Figure 2 corresponds to a star, constructed with any of
the EOS, that is spinning at the Kepler limit. The angular velocity of the
star, normalized by $\Omega_{*}$ follows a monotonic trend when plotted versus
the compactness $C_{*}$ of the non-rotating reference star. The best-fit curve
for this approximate relation for the Kepler angular velocity is given by a
4th order polynomial,

$\displaystyle\Omega_{K}=\Omega_{*}\times\sum_{i=0}^{4}\;a_{i}C_{*}^{i},$ (1)

where the coefficients $a_{i}$ and the $R^{2}$ statistic are shown in Table 1.
This best-fit curve and all other fits in this paper were made using the
Python packages numpy.polyfit (Harris et al., 2020) for polynomial fits and
scipy.optimize.curve_fit (Virtanen et al., 2020) for more complicated
functions.

Equation | Symbol | $a_{0}$ | $a_{1}$ | $a_{2}$ | $a_{3}$ | $a_{4}$ | $a_{5}$ | A | B | $R^{2}$
---|---|---|---|---|---|---|---|---|---|---
(1) | $\Omega_{K}$ | 0.552 | 3.304 | -35.211 | 180.61 | -326.48 | - | - | - | 0.9895
(3) | $R_{e}/R_{*}$ | - | -15.496 | 442.60 | -4945.62 | 23458.06 | -40544.25 | 0.203 | 0.1611 | 0.9939
(4) | $M/M_{*}$ | -0.0160 | 3.123 | -20.721 | 41.202 | -6.464 | - | 1.127 | - | 0.9808
|  | $b_{0}$ | $b_{1}$ | $b_{2}$ | $b_{3}$ | $b_{4}$ | $b_{5}$ | $A_{2}$ | $B_{2}$ | $R^{2}$
(9) | $\Omega_{K2}$ | 0.9933 | 0.4183 | -2.328 | -17.907 | 79.377 | - | - | - | 0.9121
(11) | $R_{e}/R_{*}$ | - | -8.8981 | 216.845 | -2298.43 | 10211.98 | -16587.60 | 0.39998 | -0.0058314 | 0.9980
(12) | $M/M_{*}$ | - | 1.765 | -10.985 | 11.069 | 34.996 | - | - | - | 0.9722
|  |  | $d_{1}$ | $d_{2}$ | $d_{3}$ | $d_{4}$ |  |  |  |
(12) | $M/M_{*}$ |  | 1.0 | -2.924 | 15.305 | -9.908 |  |  |  |
Table 1: Coefficients and the $R^{2}$ statistic for the best-fit equations.

Given a solution of the TOV equations, $M_{*}$ and $R_{*}$, the empirical
relation (1) predicts the largest possible angular velocity of a star that has
the same central density. We will examine the mass and radius of this star in
a later section. However, this formula does not predict whether or not this
star would instead belong to a sequence that terminates at the quasi-radial
instability, so our Kepler frequency formula has some limitations. Typical
values for the Kepler frequency are shown in Table 2. For realistic ranges of
mass and radius, the maximum spin frequency for a neutron star ranges over a
factor of 2, from about 900 to 1900 Hz.

$M_{*}$ | $R_{*}$ | $C_{*}$ | $\Omega_{K}/2\pi$ | $\nu$ | $\Omega_{n}$ | $\frac{M-M_{*}}{M_{*}}$ | $\frac{R_{e}-R_{*}}{R_{*}}$ | $\nu$ | $\Omega_{n}$ | $\frac{M-M_{*}}{M_{*}}$ | $\frac{R_{e}-R_{*}}{R_{*}}$ | $\nu$ | $\Omega_{n}$ | $\frac{M-M_{*}}{M_{*}}$ | $\frac{R_{e}-R_{*}}{R_{*}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$M_{\odot}$ | km |  $M_{\odot}/$km | Hz | Hz |  |  |  | Hz |  |  |  | Hz |  |  |
1.3 | 13 | 0.100 | 957 | 200 | 0.209 | 0.007 | 0.007 | 400 | 0.418 | 0.028 | 0.025 | 600 | 0.627 | 0.072 | 0.041
1.4 | 11 | 0.127 | 1295 | 200 | 0.154 | 0.003 | 0.004 | 400 | 0.309 | 0.015 | 0.014 | 600 | 0.463 | 0.035 | 0.027
1.6 | 13 | 0.123 | 1075 | 200 | 0.186 | 0.005 | 0.005 | 400 | 0.372 | 0.022 | 0.020 | 600 | 0.558 | 0.055 | 0.035
1.8 | 12 | 0.150 | 1309 | 200 | 0.153 | 0.003 | 0.003 | 400 | 0.305 | 0.013 | 0.013 | 600 | 0.458 | 0.033 | 0.025
2.0 | 10 | 0.200 | 1885 | 200 | 0.106 | 0.001 | 0.001 | 400 | 0.212 | 0.005 | 0.005 | 600 | 0.318 | 0.012 | 0.011
Table 2: Spin corrections for typical neutron star parameters.

The empirical relation for the Kepler frequency can be compared with earlier
calculations. Haensel & Zdunik (1989) and Friedman et al. (1989) calculated
the maximum angular velocity, $\Omega_{max}$ for a particular EOS, compared
with the mass, $M_{TOV}$ and radius, $R_{TOV}$ of the maximum mass non-
rotating star with the same EOS, and found an approximately universal scaling
of $\Omega_{max}=K(GM_{TOV}/R_{TOV}^{3})^{1/2}$. Lasota et al. (1996)
generalized this relationship so that $K$ is a function of the maximum
compactness, $C_{TOV}$ of the non-rotating star. A somewhat different approach
is to compare the maximum rotation rate of normal neutron stars with the non-
rotating neutron star with the same baryonic mass, as was done by Cook et al.
(1994) (see their Figure 27). This mapping using the star’s baryon mass does
not yield a relation that appears as universal as the mapping using the
central density shown in Figure 2. Similar universal relations for the maximum
rotation rates of neutron stars have been explored by a number of authors
(e.g. Haensel et al. (2009), Haskell et al. (2018), and Koliogiannis &
Moustakidis (2020)). We have verified that our results for the Kepler
frequency agree (within reasonable error) with earlier work.

###  3.2 Equatorial compactness ratio

The shape of a rotating neutron star is oblate, with an approximately
universal shape (Morsink et al., 2007) that depends on the equatorial
compactness $C_{e}=\frac{M}{R_{e}}\frac{km}{M_{\odot}}$ and the angular
velocity of the rotating star. The hot spots on a rotation-powered X-ray
pulsar can be at any latitude of the neutron star, but the universal relations
for the shape allow the determination of the equatorial values of radius and
compactness, through pulse profile modelling (Bogdanov et al., 2019) using
X-ray timing telescopes such as NICER.

It has been noted elsewhere (Miller et al., 2019; Raaijmakers et al., 2019;
Annala et al., 2022) that the equatorial compactness appears to be almost
constant if the star’s central density is kept constant as the star is spun
up. This property can be seen, for example, for the case of slow rotation in
Figure 5 of Raaijmakers et al. (2019), for stars spinning at the Kepler limit
in Figure 4 of Annala et al. (2022) , and for a wide range of stellar rotation
in Figure 3 of Mallick et al. (2022).

In the example EOS PP0 shown in the right-hand panel of Figure 1, the
compactness is the slope of the curve formed by connecting the sequences of
constant central density. The slopes of these sequences are approximately
constant between the TOV limit (bold black curve) and the gold curve. The gold
curve connects neutron stars spinning at an angular velocity of 0.95
$\Omega_{K}$, as predicted by Equation (1). To the right of the gold curve the
slope decreases somewhat as the sequences approach the Kepler limit, with
departures from a constant slope that are more pronounced for less compact
stars. This suggests that the equatorial compactness depends on both $C_{*}$
and the “closeness” of the star’s rotation rate to the Kepler limit. We define
a simple normalized angular velocity $\Omega_{n}$ of a particular rotating
star by

$\Omega_{n}(\Omega;C_{*},\Omega_{*})=\frac{\Omega}{\Omega_{K}(C_{*},\Omega_{*})},$
(2)

where $\Omega_{K}$ is given by the empirical fit shown in Equation (1). Each
member of a constant-central density spin-up sequence is described by the same
non-rotating reference values of $C_{*}$ and $\Omega_{*}$, so the normalized
angular velocity $\Omega_{n}$ ranges from 0 to a value close to 1 for
sequences that terminate at the Kepler limit. Note that due to the nature of
the empirical fit, the highest spin star in any sequence could be slightly
larger or smaller than 1. We use this same normalization for the sequences
that terminate in the quasi-radial instability, which occurs at values of the
normalized angular velocity that can be much smaller than 1. Typical values of
$\Omega_{n}$ for realistic masses, radii, and spin frequencies are shown in
Table 2.

Figure 3: Average normalized equatorial compactness $C_{e}/C_{*}$ for rotating
stars as a function of the compactness of the non-rotating star with the same
central density and the normalized angular velocity. The colour in each bin is
the average (over all 32 EOS) value of equatorial compactness of stars with
the same value of $C_{*}$ and $\Omega_{n}$.

The dimensionless ratio of the equatorial compactness of a rotating neutron
star to the compactness of a non-rotating star with same central density,
$C_{e}/C_{*}$, is shown in Figure 3. This plot shows the dependence of
$C_{e}/C_{*}$ on the normalized angular velocity $\Omega_{n}$ (horizontal
axis) and the non-rotating compactness (vertical axis). The colour in each bin
represents the average over all stars computed from the 32 EOS with the given
values of $\Omega_{n}$ and $C_{*}$. Although the ratio $C_{e}/C_{*}$ is not
exactly 1.0, the contours show that for stars rotating at less than about 90%
of the Kepler frequency, the decrease in compactness (ie $1-C_{e}/C_{*}$) is
less than 1% for neutron stars with astrophysical masses and radii ($C_{*}\geq
0.1$). Although the stars with $C_{*}<0.1$ have unrealistically low masses,
these stars are included on the figure in order to illustrate the dependence
of this feature on the strength of the gravitational field.

###  3.3 Equatorial radius and mass

Although $C_{e}$ is roughly constant along constant-central-density spin-up
sequences, neither $M$ nor $R_{e}$ is constant along these sequences. Typical
changes in the mass or radius are at least an order of magnitude larger than
the change in their ratio.

Other authors have examined the universality in the increase in mass for some
special cases. Breu & Rezzolla (2016) calculated the values of mass, called
the critical mass, along the quasi-radial instability line (the thin dashed
line shown in Figure 1). They found a universal dependence of the critical
masses normalized by the TOV mass on the dimensionless spin parameter $\chi$.
Shao et al. (2020) added to this relation the dependence of the critical mass
on the compactness of the TOV model, as well as a universal relation for the
radii of the critical models. These relations for the critical masses were
shown to have some EOS dependence when phase transitions and differential
rotation are included in the computations (Bozzola et al., 2019). In this
paper, we examine a more general universality in the increase in mass (and
radius) for all uniformly rotating neutron stars, not just the critical stars.

The rotationally-induced increase in the neutron star’s equatorial radius and
mass is shown in Figures 4 and 5 respectively. The average values and standard
deviations of $R_{e}/R_{*}$ and $M/M_{*}$ (over all 32 EOS) are shown in each
bin of constant $C_{*}$ and $\Omega_{n}$. Figure 4 shows that on average (left
panel), the radius increases by about 15% when the neutron star is rotating at
about 80% of the Kepler limit. The standard deviation in the radius increase
is less than 1% of $R_{*}$ for the same range. Since the increase in radius is
almost independent of the EOS, and only depends on $C_{*}$ and $\Omega_{n}$,
the increase appears to be “universal”.

Figure 4: Average (left) and standard deviation (right) of the increase in
equatorial radius due to rotation as a function of $\Omega_{n}$ and $C_{*}$.
The colour bar represents the average (or standard deviation) of the value of
$R_{e}/R_{*}$ for all 32 EOS in each bin.

Similarly, the average and standard deviation in the increase in the mass due
to rotation is shown in Figure 5. The mass increases by 15% at a somewhat
larger normalized angular velocity, and there is more dependence on
compactness than is seen in Figure 4 for the increase in radius. The standard
deviation in $M/M_{*}$ grows as large as 1% at spins of 80% of the Kepler
limit. While this is a larger variation than seen in the increase in radius,
this is still almost independent of the EOS.

Figure 5: Average (left) and standard deviation (right) in the increase in
mass due to rotation as a function of $\Omega_{n}$ and $C_{*}$. The colour bar
represents the average (or standard deviation) of the value of $M/M_{*}$ for
all 32 EOS.

Based on these results, we made the hypothesis that the changes in mass and
equatorial radius along a constant-density spin-up sequence can be
approximated by a surface in a space spanned by $C_{*}$ and $\Omega_{n}$. The
best-fit equations for the surfaces representing these quantities are

$\displaystyle\frac{R_{e}}{R_{*}}$ $\displaystyle=$ $\displaystyle
1+(e^{A_{r}\Omega_{n}^{2}}-1+B_{r}\left[\ln(1-(\frac{\Omega_{n}}{1.1})^{4})\right]^{2})\times\left(1+\sum_{i=1}^{5}a_{r,i}C_{*}^{i}\right)$
(3) $\displaystyle\frac{M}{M_{*}}$ $\displaystyle=$ $\displaystyle
1+(e^{A_{m}\Omega_{n}^{2}}-1)\times\left(\sum_{i=0}^{4}a_{m,i}C_{*}^{i}\right),$
(4)

where the best-fit coefficients and the $R^{2}$ statistics for the fits are
shown in Table 1. The functional forms of these equations were found by first
plotting $R_{e}$ and $M$ as surfaces depending on $C_{*}$ and $\Omega_{n}^{2}$
and noting the approximately exponential dependence on the square of the
angular velocity by eye. The logarithmic term was added to model the rapid
change in radius as the Kepler limit is approached. The factor of 1.1 was
required since the empirical nature of the fits leads to some stars having
values of $\Omega_{n}>1.$

In order to quantify the quality of these best-fit equations, we introduce the
absolute deviation of some quantity $Q$. The absolute deviation of $Q$,
$Dev(Q)$ for a particular rotating neutron star model $Q_{i}$ from the best-
fit surface $Q_{fit}$ is defined by

$Dev(Q)=\frac{|Q_{i}-Q_{fit}|}{Q_{i}},$ (5)

which is the fractional error incurred by using the best-fit function.
Examination of the absolute deviations is useful for flagging possible outlier
EOS that are not described well by the best-fit functions.

The largest absolute deviations introduced by the best-fit value of
$R_{e}/R_{*}$, Equation (3) for any of the EOS are shown in Figure 6 (upper
panels). The colour in each pixel represents the largest absolute deviation of
the fitting function for the scaled equatorial radius $R/R_{*}$, from the
value computed for any star with the given value of $C_{*}$ and $\Omega_{n}$.
The left-hand panel shows the largest absolute deviation within the 32
randomly generated EOS while the right-hand panel shows the largest deviation
for the smaller set of tabulated hadronic and hybrid EOS described in the last
paragraph of Section 2.1. Similarly, the largest absolute deviation of
$M/M_{*}$ for the same groups of EOS is shown in the lower panels of 6.

Figure 6 shows that for values of $\Omega_{n}\leq 0.8$ the error introduced by
the formula for the rotating star’s radius given by Equation (3) is typically
less than 2% hadronic and hybrid EOS. The deviations for $M/M_{*}$ are larger
than the radius deviations, and more equation of state dependence can be seen
in this Figure. Stars with $\Omega_{n}=0.8$ can have deviations in the
rotating mass that are as large as 4%, depending on the value of $C_{*}$ and
the EOS. The largest deviations come about mainly from the hyperonic EOS,
which has a very low maximum mass of 1.5 $M_{\odot}$. Bare quark stars show a
larger deviation, as large as 5% in mass at 60% of Kepler. As has been seen
for other universal relations, the bare quark EOS probably should be treated
separately from hadronic EOS. We speculate that the larger EOS dependence seen
in the $M/M_{*}$ plots (for all EOS) is due to some dependence on how close
the star’s mass is to the maximum mass.

Figure 6: Largest absolute deviation in the radius (upper) and mass (lower).
The colour in each $C_{*}$-$\Omega_{n}$ bin represents the largest absolute
deviation in the ratio $R_{e}/R_{*}$ or $M/M_{*}$ for any of the 32 randomly
generated EOS (left) or tabulated hadronic or hybrid EOS (right).

The absolute deviations in the best-fit formulae can rise to about 8% error
for stars rotating close to the Kepler limit. This is not surprising, since it
has been seen that other universal relations, for instance the I-Love-Q
relations break down for very rapid rotation (Doneva et al., 2014).

###  3.4 Approximate universality

The approximately constant value for the equatorial compactness on a constant-
density spin-up sequence can be understood intuitively through a simple
Newtonian model. In Euclidean space, the volume of an oblate spheroid with
equatorial radius $R_{e}$ and polar radius $R_{p}$ is $V=4\pi
R_{e}^{2}R_{p}/3$. Suppose that this ellipsoid is a member of the same
sequence of models as a non-rotating sphere with radius $R_{*}$ and volume
$V_{*}$, with $R_{e}=R_{*}+\Delta R_{e}$ and $R_{p}=R_{*}-\Delta R_{p}$. For
small changes in radius, the change in volume is

$\Delta V=V_{*}\frac{R_{e}}{R_{*}}\left(2\frac{\Delta R_{e}}{R*}-\frac{\Delta
R_{p}}{R*}\right)\simeq V_{*}\frac{\Delta R_{e}}{R_{*}},$ (6)

if $\Delta R_{p}\simeq\Delta R_{e}$.

In the case of a uniform density star, if the central density is kept
constant, then the logarithmic change in the total mass of the star will
increase at the same rate as the logarithmic change in the volume, so

$\frac{\Delta M}{M_{*}}=\frac{\Delta V}{V_{*}}.$ (7)

In the situations where $\Delta R_{p}\simeq\Delta R_{e}$, this leads to the
logarithmic change in the mass being roughly equal to the logarithmic change
in the equatorial radius, leading to constant compactness along constant
central-density sequences.

In reality, the density is not constant inside a neutron star. However, in the
case of Newtonian polytropes, the average density of a star is equal to a
constant (for a particular EOS) times the central density. If this property is
preserved as a star spins faster, then a sequence with constant central-
density will also have constant average density, leading to equation (7).

In order to quantitatively show how the polar radius of a star changes
relative to the equatorial radius, we introduce a quantity, $\alpha$, that
measures the asymmetry in the rotational deformation, defined by

$\alpha=\frac{2R_{*}-R_{e}-R_{p}}{2R_{*}}.$ (8)

The rotational deformation asymmetry, $\alpha$, vanishes for zero rotation by
definition. Rotating stars with $\Delta R_{e}=\Delta R_{p}$ will also have
vanishing $\alpha$. Note that a star can be highly elliptical even if $\alpha$
vanishes. The definition of $\alpha$ is not a coordinate-independent
definition. In our application, we compute $\alpha$ using the Schwarzschild-
like radial coordinate, since it is most closely connected to physically
meaningful quantities like the circumference. Figure 7 shows the largest
absolute value of $\alpha$ for each value of $C_{*}$ and $\Omega_{n}$. The
values of $\alpha$ are small compared to the values of $\Delta R_{e}/R_{*}$
shown in Figure 4, and trends to smaller values with increasing compactness.
The vanishing of the rotational deformation asymmetry in the limit of very
compact stars is most likely related to the iso-density self-similarity
discussed by Yagi & Yunes (2017).

Figure 7: Neutron star deformation asymmetry. The colour scale shows the
largest absolute value (over the 32 EOS) of the dimensionless rotational
deformation asymmetry, $(2R_{*}-R_{e}-R_{p})/(2R_{*})$, for each bin of
$C_{*}$ and $\Omega_{n}$.

When we examine the constant-density contours inside a rotating neutron star,
we see a general qualitative trend that the deformation of the star is
dominated by the deformation of the crust. This suggests a simple toy model of
a dense, incompressible core that stays roughly spherical and a low density
crust that is easily deformed by the centrifugal force. If we think of the two
regions as two springs in series, with a stiff spring attached to a wall
representing the core, and a soft spring representing the crust, the overall
effective spring constant will have a value similar to the soft spring’s
constant (in the limit that the crust’s spring constant is much smaller than
the core’s spring constant). The applied centrifugal force, proportional to
$\Omega_{n}$, then induces a strain in the crust that depends on $\Omega_{n}$,
but is almost independent of the length of the spring (or thickness of the
crust). Since most of the EOS in our libraries (including the tabulated
hadronic and hybrid EOS) have similar crust EOS, it is natural that the change
in equatorial radius with spin should show the almost universal behaviour, as
seen in Figure 4. The bare quark star EOS does not include a crust, and with
this toy model it is not surprising that bare quarks stars are not as well
described by universal relations.

###  3.5 Application to EOS inference

The empirical relations for the mass and radius of a rotating neutron star
given the values found from the spherically symmetric TOV equations with the
same central density lead to a simple set of corrections that can be used in
EOS inference codes.

We imagine that a number of neutron stars’ masses, equatorial radii, and spins
have been measured, for instance using pulse-profile modelling as is done by
NICER. If we have $N$ measurements, the $i$-th neutron star has spin frequency
$\nu_{i}$, and the mass and equatorial radius will be described by a
probability distribution. The spin-corrected EOS inference procedure is to
first choose an EOS, and then solve the TOV equations for a range of central
densities to construct a zero-spin mass-radius curve. For each value of
central density use the zero-spin values $M_{*}$, $R_{*}$, and the spin
frequency of the $i$th star to compute the predicted mass and radius using
equations (4) and (3). The likelihood of this empirically predicted mass and
radius can then be computed using the $i$th mass-radius distribution, and then
the procedure repeated for other central densities until the central density
that maximizes the $i$th likelihood function has been found (for this specific
choice of EOS). Since the empirical relations are simple algebraic
expressions, they are computationally inexpensive to implement. This makes it
simple to add rotation to an EOS inference code that simply solves the TOV
equations.

Whether or not these empirical corrections to the rotating star’s mass and
radius should be added to an EOS inference code depends on the precision of
the measurements, and the relative increase in the star’s mass and radius due
to rotation. For instance, NICER measured the mass and radius of the pulsar
J0030 with a precision of about 10% in both quantities (Miller et al., 2019;
Riley et al., 2019). This pulsar’s spin frequency is 205 Hz so comparing with
values in Table 2, the increases in mass and radius are less than 1% for
typical mass and radius values. The small size of these increases means that
the simple use of the TOV equations in the EOS inference results of Miller et
al. (2019) and Raaijmakers et al. (2019) is valid. It is only when the
precision of mass and radius parameter estimation for the NICER pulsars is
similar to the typical magnitudes of the changes given in Table 2 will the
rotational corrections presented in this paper be significant.

##  4 Inverse mapping

Given that we have found a mapping between the mass and radius of non-rotating
neutron stars to rotating neutron stars, it is natural to look at inverse
mappings. The inverse mapping is potentially useful if many neutron stars,
with different spins, have precise mass and equatorial radius measurements. An
inverse mapping allows the construction of an effective zero-spin mass-radius
relation using multiple observations.

The inverse map from the measurement of a rotating star’s mass and radius to a
non-rotating star with the same central energy density can be constructed
using a method similar to Section 3. The first step is to construct a new
normalized angular velocity, since the previous definition relies on knowledge
of $M_{*}$ and $R_{*}$. A natural choice is to compute the Kepler angular
velocity normalized by $\sqrt{GM/R_{e}^{3}}$, where $M$ and $R_{e}$ are the
mass and equatorial radius of a neutron star rotating at the Kepler limit.
Figure 8 shows this normalized Kepler frequency plotted versus the equatorial
compactness for all 32 EOS. Each blue or orange point represents a maximally
rotating star from one of the constant density sequences computed using either
a PP or CS EOS. The solid curve shows a best-fit curve for these points, which
remarkably, is within a couple percent of 1. The equation of this curve is

$\Omega_{K2}=\sqrt{\frac{GM}{R_{e}^{3}}}\times\sum_{i=0}^{4}b_{i}C_{e}^{i}.$
(9)

The best-fit coefficients for this equation can be found in Table 1.

Figure 8: The normalized Kepler limit angular velocity versus the compactness
ratio of the rotating neutron star. Blue and orange points represent neutron
stars spinning at the Kepler limit computed with either a PP or CS EOS
respectively. The best-fit equation (9) to these points is shown as a solid
curve.

The Kepler angular velocity is, by definition, both the angular velocity of
the star spinning at the mass-shed limit and the angular velocity of a
particle orbiting the star at the equator (in the same direction as the star’s
rotation). The orbital angular velocity of a particle orbiting a non-rotating
star is $\Omega_{s,orb}=\sqrt{GM/R^{3}}$. Since the formula for the orbital
frequency at the mass-shed limit is almost the same, this suggests that the
orbital frequency at intermediate spins be given by a similar formula
$\Omega_{orb}(\Omega)=F\sqrt{GM/R_{e}^{3}}$, where $F\sim 1$. When we plot
(not shown) the orbital frequencies for all models, we find that while there
is more scatter, the constant is always within 3% of 1. This does not depend
on the stability of the orbit. We have not attempted to fit this relationship
for the orbital angular velocity for intermediate spin rates, but note that
approximating $F\sim 1$ may be good enough for other applications.

The spin parameter $\chi$ is approximately proportional to the star’s angular
velocity. In the case of stars rotating at the Kepler frequency, it has been
argued (Breu & Rezzolla, 2016; Shao et al., 2020) that since
$\Omega_{K}\sim\sqrt{GM/R^{3}}$, the value of the spin parameter at the Kepler
limit should be $\chi_{K}\sim I/(MR^{2})\times\sqrt{Rc^{2}/GM}$, where $I$ is
the moment of inertia, and $I,M,$ and $R$ are all evaluated at the Kepler
limit. When combined with the universal relations for $I$ (Koliogiannis &
Moustakidis, 2020) at the Kepler limit, Shao et al. (2020) showed that this
leads to a very simple formula for the largest value of $\chi$ allowed by an
equation of state, $\chi_{max}=1.24(\sqrt{GM/Rc^{2}})_{TOV}$. This leads to a
prediction of about a 10% variation in the maximum possible value of $\chi$,
depending on the EOS. It would be interesting to investigate the dependence of
mass and radius on dimensionless parameters such as $\chi/\chi_{K}$. However,
our preliminary investigations found that the dependence on $\chi/\chi_{K}$
introduced more scatter than occurs with a quantity such as
$\Omega/\Omega_{K2}$ that traces the relative magnitudes of the centrifugal
and gravitational forces.

Given an observation of a rotating neutron star’s mass, equatorial radius, and
angular velocity, we can compute the quantities $C_{e}$ and a normalized
angular velocity

$\Omega_{n2}=\frac{\Omega}{\Omega_{K2}}.$ (10)

Based on our library of neutron star models, we find the following best-fit
equations for the mass and radius of a non-rotating neutron star with the same
central density as the rotating star,

$\displaystyle\frac{R_{e}}{R_{*}}$ $\displaystyle=$ $\displaystyle
1+(e^{A_{2}\Omega_{n}^{2}}-1+B_{2}\left[\ln(1-(\frac{\Omega_{n2}}{1.1})^{4})\right]^{2})\times\left(1+\sum_{i=1}^{5}b_{r,i}C_{e}^{i}\right)$
(11) $\displaystyle\frac{M}{M_{*}}$ $\displaystyle=$ $\displaystyle
1+\left(\sum_{i=1}^{4}d_{i}\Omega_{n2}^{i}\right)\times\left(\sum_{i=1}^{4}b_{m,i}C_{e}^{i}\right),$
(12)

where the coefficients appearing in these equations can be found in Table 1.
If we have a few neutron star mass-radius measurements, these equations can be
used to produce an effective zero-spin mass-radius curve for this group of
neutron stars.

Figure 1 demonstrates the use of the inverse map on the example EOS PP0. Each
blue dot representing a rotating neutron star has been mapped to a non-
rotating neutron star with the same central density, represented with an
orange dot, using Equations (9) - (12). If the inverse map were perfect, all
of the orange dots for a particular central density would overlap and lie
exactly on the solid black TOV curve. It can be seen in Figure 1 that there is
a small spread in the orange dots, however, the spread is much smaller than
the typical uncertainty in the measurements.

NICER observations of PSR J0030+0451, which spins with a frequency of
$\nu=205$ Hz, were analyzed by Miller et al. (2019) and Riley et al. (2019)
leading to estimates for this pulsar’s mass and equatorial radius with
precisions (at one sigma) of about 10%. For example, Miller et al. (2019)
found values M = $1.44^{+0.15}_{-0.14}$ $M_{\odot}$ and $R_{e}$ =
$13.02^{+1.24}_{-1.06}$ km for the rotating star. Making use of Equations (9)
- (12) with $M=1.44M_{\odot}$, $R_{e}=13.02$ km, and $\nu=200$ Hz, we find
that the equivalent non-rotating mass and radius are $M_{*}=1.43M_{\odot}$ and
$R_{*}=12.94$ km. These changes in mass and radius due to rotation are about
an order of magnitude smaller than the measurement precision, so there is no
loss of accuracy if the rotational correction terms are ignored.

The precision of the mass and radius measurements made by NICER will improve
as more data is accumulated with longer observations. However, NICER’s target
precision of 5% for the radius measurement is still larger than the $~{}1\%$
typical change in radius predicted for neutron stars spinning at the
relatively slow rate of 200 Hz. The other pulsars observed by NICER spin at
similar rates, so the order of magnitudes for their spin corrections are
similar. However, future missions, such as Strobe-X (Ray et al., 2019) and
eXTP (Watts et al., 2019) will have a 10x larger collecting area, which will
allow for more precise radius measurements as well as observations of fainter
pulsars. In addition, these missions are likely to target accreting ms-period
X-ray pulsars that spin more rapidly. The universal formulae for the changes
in mass and radius will be more important for these future X-ray missions.

##  5 Discussion

Neutron star properties such as their masses and radii are decidedly non-
universal, since they depend strongly on the EOS of supranuclear density
matter. However, given a value of a neutron star’s mass and radius, many
secondary properties only appear to depend on the given mass and radius, and
are almost independent of the choice of EOS. In this paper we investigate the
universality of the change in mass and radius caused by rapid rotation.

We computed equilibrium sequences of rapidly rotating neutron stars with
constant central density to explore the universality of the increase in mass
and equatorial radius. We found that the fractional increase of mass and
radius can be fit with an empirical formula that depends only on the mass and
radius of the non-rotating star in the sequence and the rotating star’s spin
frequency. As part of this method, we have also derived a formula for the
fastest spin frequency of a neutron star given the mass and radius of a non-
rotating neutron star with the same central density. These empirical formulae
can then be applied as corrections to EOS inference codes, since EOS inference
codes compute the mass and radius of a non-rotating star with a given central
density.

We also explored the surprising property that the equatorial compactness of a
neutron star is very close to constant along a constant central-density
sequence. This property appears to be related to an internal symmetry that
keeps the relative decrease in the polar radius close to equal in magnitude to
the increase in the equatorial radius.

It is also interesting to consider whether an effective zero-spin mass-radius
relation can be constructed from observations of the masses and radii of many
neutron stars with a wide variety of spin frequencies. We constructed a set of
equations that allow the mapping of a rotating star’s mass and radius to a
corresponding non-rotating neutron star with the same central density.

The rotational corrections derived in this paper are generally smaller than
the typical uncertainties in determining the masses and radii of rotating
neutron stars. So at present, these correction factors are not required in EOS
inference codes used with results obtained from pulse-profile modelling done
by NICER. However, we expect these formulae to be useful in the future when
more precise measurements are possible, and when radius measurements of more
rapidly rotating neutron stars are possible.

##  6 Acknowledgements

We thank Jorge Calderon Noguez, Rodrigo Fernandez, and Craig Heinke for
helpful discussions, and the anonymous referee for useful comments. This
research was supported by a grant from NSERC to SMM.

## References

  * Akmal et al. (1998) Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Physical Review C, 58, 1804, doi: 10.1103/PhysRevC.58.1804
  * Alcock et al. (1986) Alcock, C., Farhi, E., & Olinto, A. 1986, ApJ, 310, 261, doi: 10.1086/164679
  * Alford et al. (2005) Alford, M., Braby, M., Paris, M., & Reddy, S. 2005, ApJ, 629, 969, doi: 10.1086/430902
  * AlGendy & Morsink (2014) AlGendy, M., & Morsink, S. M. 2014, ApJ, 791, 78, doi: 10.1088/0004-637X/791/2/78
  * Annala et al. (2022) Annala, E., Gorda, T., Katerini, E., et al. 2022, Physical Review X, 12, 011058, doi: 10.1103/PhysRevX.12.011058
  * Baldo et al. (1997) Baldo, M., Bombaci, I., & Burgio, G. F. 1997, A&A, 328, 274.  https://arxiv.org/abs/astro-ph/9707277
  * Bogdanov et al. (2019) Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJ, 887, L26, doi: 10.3847/2041-8213/ab5968
  * Bozzola et al. (2019) Bozzola, G., Espino, P. L., Lewin, C. D., & Paschalidis, V. 2019, European Physical Journal A, 55, 149, doi: 10.1140/epja/i2019-12831-2
  * Breu & Rezzolla (2016) Breu, C., & Rezzolla, L. 2016, MNRAS, 459, 646, doi: 10.1093/mnras/stw575
  * Cook et al. (1992) Cook, G. B., Shapiro, S. L., & Teukolsky, S. A. 1992, ApJ, 398, 203, doi: 10.1086/171849
  * Cook et al. (1994) —. 1994, ApJ, 424, 823, doi: 10.1086/173934
  * Doneva et al. (2014) Doneva, D. D., Yazadjiev, S. S., Stergioulas, N., & Kokkotas, K. D. 2014, ApJ, 781, L6, doi: 10.1088/2041-8205/781/1/L6
  * Friedman et al. (1989) Friedman, J. L., Ipser, J. R., & Parker, L. 1989, Phys. Rev. Lett., 62, 3015, doi: 10.1103/PhysRevLett.62.3015
  * Gendreau et al. (2012) Gendreau, K. C., Arzoumanian, Z., & Okajima, T. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8443, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, ed. T. Takahashi, S. S. Murray, & J.-W. A. den Herder, 844313, doi: 10.1117/12.926396
  * Glendenning (1997) Glendenning, N. K. 1997, Compact stars. Nuclear physics, particle physics, and general relativity
  * Greif et al. (2019) Greif, S. K., Raaijmakers, G., Hebeler, K., Schwenk, A., & Watts, A. L. 2019, MNRAS, 485, 5363, doi: 10.1093/mnras/stz654
  * Hachisu (1986) Hachisu, I. 1986, ApJS, 61, 479, doi: 10.1086/191121
  * Haensel & Zdunik (1989) Haensel, P., & Zdunik, J. L. 1989, Nature, 340, 617, doi: 10.1038/340617a0
  * Haensel et al. (2009) Haensel, P., Zdunik, J. L., Bejger, M., & Lattimer, J. M. 2009, A&A, 502, 605, doi: 10.1051/0004-6361/200811605
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Haskell et al. (2018) Haskell, B., Zdunik, J. L., Fortin, M., et al. 2018, A&A, 620, A69, doi: 10.1051/0004-6361/201833521
  * Hebeler et al. (2013) Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11, doi: 10.1088/0004-637X/773/1/11
  * Kojo et al. (2021) Kojo, T., Baym, G., & Hatsuda, T. 2021, arXiv e-prints, arXiv:2111.11919.  https://arxiv.org/abs/2111.11919
  * Koliogiannis & Moustakidis (2020) Koliogiannis, P. S., & Moustakidis, C. C. 2020, Phys. Rev. C, 101, 015805, doi: 10.1103/PhysRevC.101.015805
  * Komatsu et al. (1989) Komatsu, H., Eriguchi, Y., & Hachisu, I. 1989, MNRAS, 237, 355, doi: 10.1093/mnras/237.2.355
  * Krüger & Kokkotas (2020) Krüger, C. J., & Kokkotas, K. D. 2020, Phys. Rev. Lett., 125, 111106, doi: 10.1103/PhysRevLett.125.111106
  * Landry & Essick (2019) Landry, P., & Essick, R. 2019, Phys. Rev. D, 99, 084049, doi: 10.1103/PhysRevD.99.084049
  * Lasota et al. (1996) Lasota, J.-P., Haensel, P., & Abramowicz, M. A. 1996, ApJ, 456, 300, doi: 10.1086/176650
  * Lattimer & Prakash (2001) Lattimer, J. M., & Prakash, M. 2001, ApJ, 550, 426, doi: 10.1086/319702
  * Legred et al. (2022) Legred, I., Chatziioannou, K., Essick, R., & Landry, P. 2022, Phys. Rev. D, 105, 043016, doi: 10.1103/PhysRevD.105.043016
  * Lindblom (1992) Lindblom, L. 1992, ApJ, 398, 569, doi: 10.1086/171882
  * Lindblom & Indik (2012) Lindblom, L., & Indik, N. M. 2012, Phys. Rev. D, 86, 084003, doi: 10.1103/PhysRevD.86.084003
  * Lindblom & Indik (2014) —. 2014, Phys. Rev. D, 89, 064003, doi: 10.1103/PhysRevD.89.064003
  * Mallick et al. (2022) Mallick, R., Kuzur, D., & Nandi, R. 2022, European Physical Journal C, 82, 512, doi: 10.1140/epjc/s10052-022-10468-w
  * Miller et al. (2020) Miller, M. C., Chirenti, C., & Lamb, F. K. 2020, ApJ, 888, 12, doi: 10.3847/1538-4357/ab4ef9
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24, doi: 10.3847/2041-8213/ab50c5
  * Miller et al. (2021) —. 2021, ApJ, 918, L28, doi: 10.3847/2041-8213/ac089b
  * Morsink et al. (2007) Morsink, S. M., Leahy, D. A., Cadeau, C., & Braga, J. 2007, ApJ, 663, 1244, doi: 10.1086/518648
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, Annual Review of Astronomy and Astrophysics, 54, 401, doi: 10.1146/annurev-astro-081915-023322
  * Pandharipande & Smith (1975) Pandharipande, V. R., & Smith, R. A. 1975, Nuclear Physics, Section A, 237, 507, doi: 10.1016/0375-9474(75)90415-7
  * Raaijmakers et al. (2019) Raaijmakers, G., Riley, T. E., Watts, A. L., et al. 2019, ApJ, 887, L22, doi: 10.3847/2041-8213/ab451a
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, ApJ, 918, L29, doi: 10.3847/2041-8213/ac089a
  * Ravenhall & Pethick (1994) Ravenhall, D. G., & Pethick, C. J. 1994, ApJ, 424, 846, doi: 10.1086/173935
  * Ray et al. (2019) Ray, P. S., Arzoumanian, Z., Ballantyne, D., et al. 2019, arXiv e-prints, arXiv:1903.03035.  https://arxiv.org/abs/1903.03035
  * Read et al. (2009) Read, J. S., Lackey, B. D., Owen, B. J., & Friedman, J. L. 2009, Physical Review D, 79, doi: 10.1103/physrevd.79.124032
  * Rezzolla et al. (2018) Rezzolla, L., Most, E. R., & Weih, L. R. 2018, ApJ, 852, L25, doi: 10.3847/2041-8213/aaa401
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21, doi: 10.3847/2041-8213/ab481c
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27, doi: 10.3847/2041-8213/ac0a81
  * Shao et al. (2020) Shao, D.-S., Tang, S.-P., Sheng, X., et al. 2020, Phys. Rev. D, 101, 063029, doi: 10.1103/PhysRevD.101.063029
  * Silva et al. (2021) Silva, H. O., Pappas, G., Yunes, N., & Yagi, K. 2021, Phys. Rev. D, 103, 063038, doi: 10.1103/PhysRevD.103.063038
  * Stergioulas & Friedman (1995) Stergioulas, N., & Friedman, J. L. 1995, ApJ, 444, 306, doi: 10.1086/175605
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Völkel & Krüger (2022) Völkel, S. H., & Krüger, C. J. 2022, Phys. Rev. D, 105, 124071, doi: 10.1103/PhysRevD.105.124071
  * Watts et al. (2016) Watts, A. L., Andersson, N., Chakrabarty, D., et al. 2016, Reviews of Modern Physics, 88, 021001, doi: 10.1103/RevModPhys.88.021001
  * Watts et al. (2019) Watts, A. L., Yu, W., Poutanen, J., et al. 2019, Science China Physics, Mechanics, and Astronomy, 62, 29503, doi: 10.1007/s11433-017-9188-4
  * Yagi & Yunes (2013) Yagi, K., & Yunes, N. 2013, Science, 341, 365, doi: 10.1126/science.1236462
  * Yagi & Yunes (2017) —. 2017, Phys. Rep., 681, 1, doi: 10.1016/j.physrep.2017.03.002