# A low-crosstalk double-side addressing system using acousto-optic deflectors
for atomic ion qubits

Rui-Rui Li  Yi-Long Chen  Ran He  Shu-Qian Chen  Wen-Hao Qi  CAS Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei 230026, China  CAS Center For Excellence in Quantum Information
and Quantum Physics, University of Science and Technology of China, Hefei
230026, China  Jin-Ming Cui  Electronic mail<EMAIL_ADDRESS>Yun-Feng
Huang  Electronic mail<EMAIL_ADDRESS>Chuan-Feng Li  Guang-Can Guo  CAS Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei 230026, China  CAS Center For Excellence in Quantum Information
and Quantum Physics, University of Science and Technology of China, Hefei
230026, China

###### Abstract

The ability to individually and agilely manipulate qubits is crucial for the
scalable trapped-ion quantum information processing. A plethora of challenging
proposals have been demonstrated with the utilization of optical addressing
systems, in which single ions is addressed exclusively by individual laser
beam. However, crosstalk error in optical addressing systems limits the gate
fidelity, becoming an obstacle to quantum computing, especially quantum error
correction. In this work, we demonstrate a low-crosstalk double-side
addressing system based on a pair of acousto-optic deflectors (AODs). The AODs
addressing method can flexibly and parallelly address arbitrary ions between
which the distance is variable in a chain. We employ two 0.4 NA objective
lenses in both arms of the Raman laser and obtain a beam waist of 0.95
$\upmu\mathrm{m}$, resulting in a Rabi rate crosstalk as low as $6.32\times
10^{-4}$ when the neighboring ion separation is about 5.5 $\upmu\mathrm{m}$.
This agile and low-crosstalk double-side addressing system is promising for
higher-fidelity gates and the practical application of the quantum error
correction.

††preprint: AIP/123-QED

Trapped ion is one of the most promising platforms of the quantum information
process (QIP) due to its longest coherent time Wang  _et al._ (2017, 2021),
high-fidelity state initialization, detection An  _et al._ (2022) and logic
manipulations Harty  _et al._ (2014); Gaebler  _et al._ (2016); Ballance  _et
al._ (2016); Clark  _et al._ (2021) among all candidates. The ability to
individually and agilely manipulate an arbitrary qubit in trapped-ion systems
brings a plethora of challenging proposals to reality, e.g., the parallel
entangling operations Figgatt  _et al._ (2019), the quantum Grover search
Figgatt  _et al._ (2017) algorithms, the quantum simulation in many-body
systems Smith  _et al._ (2016); Martinez  _et al._ (2016), and the quantum
error correction (QEC) Nguyen  _et al._ (2021); Schindler  _et al._ (2011);
Erhard  _et al._ (2021); Egan  _et al._ (2021). However, the crosstalk error
in individually addressing (IA) systems limits the fidelity of gate operation
Schmidt-Kaler  _et al._ (2003), resulting in an obstacle for the QIP,
especially the QEC, which has more stringent requirements Parrado-Rodríguez
_et al._ (2021). Therefore, a low-crosstalk IA system is the key to realizing
a universal error-corrected QIP platform.

The intensity crosstalk, resulting from the IA beam’s intensity distribution
at its neighboring sites, is strongly related to the spacing of ions
(typically 3 to 6 microns in the practical experiment) and spot sizes of IA
beams (much smaller than the ion spacing). We list the intensity crosstalk in
several previously reported IA systems, shown in Tab. 1. The intensity
crosstalks in those IA systems with 0.4 NA objective lens are typically on the
level at $10^{-4}$. To reduce the intrinsic intensity crosstalk of IA systems,
one can employ a higher NA objective lens to reduce the spot sizes of IA
beams. Some IA systems with 0.6 NA objective lens acquire an intensity
crosstalk on the level at $10^{-6}$. However, the NA of the objective lens
can’t be increased unboundedly and the employment of such a high NA objective
lens often causes extra aberrations resulting from misalignment. Hence the
intrinsic intensity crosstalk of IA systems is difficult to be further
improved and is limited by optical engineering.

Table 1: Parameters of different IA systems. The best SSA system that employs 0.6 NA objective lens acquire a Rabi rate crosstalk as low as 0.2%, yielding a low intensity crosstalk of $4\times 10^{-6}$. Our DAS system employs a lower NA objective lens to acquire the lowest Rabi rate crosstalk, relaxing the requirement of intensity crosstalk. S. Crain, et. al. have reported a MEMS-based IA system with a Rabi rate crosstalk of $1.3\times 10^{-4}$ but the ion spacing is 7.4 $\upmu\mathrm{m}$ Crain  _et al._ (2014) hence we don’t list it here.  IA system | NA | Intensity crosstalk | Rabi crosstalk | ion spacing
---|---|---|---|---
MCAOMs Egan (2021) | 0.37 | 1e-4 | 1%-2% | 5 $\upmu\mathrm{m}$
GLIAS Binai-Motlagh  _et al._ (2023) | 0.37 | 1e-4 | 1%(estimate) | 4 $\upmu\mathrm{m}$
AODs Pogorelov  _et al._ (2021) | 0.6 | 4e-6 | 0.2% | 3.5 $\upmu\mathrm{m}$
MEMS Wang  _et al._ (2020) | 0.6 | 4e-6 | 0.2%-0.6% | 5 $\upmu\mathrm{m}$
This work | 0.37$\times$2 | <1e-3 | 0.06%-0.1% | 5.5 $\upmu\mathrm{m}$

Intensity crosstalk eventually leads to the crosstalk error, which is the
unwanted Rabi drive of neighboring ions caused by the addressing laser. For
those single-side addressing (SSA) systems, i.e. one configures a global beam
on one side and tightly focusing IA beams on the other side, the Rabi rate
crosstalk is proportional to the square root of the intensity crosstalk of IA
beams, shown in Fig. 1(a). The SSA system with the best intensity crosstalk
performance (on the level at $10^{-6}$) has a Rabi rate crosstalk of 0.2%,
shown in Tab. 1. As a comparison, the Rabi rate crosstalk scales linearly with
the intensity crosstalk of IA beams in the case of double-side addressing
(DSA) methods, i.e. one configures tightly focusing IA beams on the two sides.
Consequently, the DSA system can reach much lower crosstalk, breaking the
limitation of optical engineering suffered by SSA systems. From Tab. 1, in our
DSA system, one can employ a lower NA objective lens and moderate optical
system to achieve lower Rabi rate crosstalk than the best SSA system.

Figure 1: Schematic diagrams of different IA systems. (a) The single-side
addressing (SSA) methods typically configure a global beam and IA beams to
drive the SRT. Rabi frequency is proportional to the square root of light
field intensity of two Raman lasers, $\Omega\propto\sqrt{I_{1}}\sqrt{I_{2}}$.
The spot size of the global beam is sufficiently large to cover all the ion
qubits hence its intensity distribution is uniform at different positions. In
this case, we have
$C_{\mathrm{SSA}}\equiv\sqrt{I_{\mathrm{IA}}\left(x_{1}\right)/I_{\mathrm{IA}}\left(x_{0}\right)}$.
In particular, for optical qubits, the above equation also holds when using
SSA methods. (b) The double-side addressing (DSA) method reported in this work
configure tightly focusing IA beams on both sides of the Raman laser. In this
case, we have $C_{\mathrm{DSA}}\equiv
I_{\mathrm{IA}}\left(x_{1}\right)/I_{\mathrm{IA}}\left(x_{0}\right)$, assuming
the identical intensity distribution of IA beams from two sides. In this
figure, the position of the addressed ion and neighboring ion is denoted as
$x_{0}$ and $x_{1}$, respectively.

The previous demonstratations of the AOD-based SSA system have been used to
manipulate ${{}^{40}}\mathrm{Ca}^{+}$ optical qubit Pogorelov  _et al._
(2021). In this work, we firstly demonstrate a low-crosstalk DSA system using
a pair of 0.4 NA objective lenses and AODs to manipulate
${{}^{171}}\mathrm{Yb}^{+}$ hyperfine qubit via stimulated Raman transition
(SRT) Bergmann, Theuer, and Shore (1998); Vitanov  _et al._ (2017); Sørensen
_et al._ (2006). This free space AOD-based DSA system enables a flexible,
reprogramming and low error-rate manipulation for atomic ion qubits. It
addresses ions that are not evenly distributed in a long ion chain by flexibly
tuning the driven frequencies of AODs. By moderate optical configuration and
simple alignment procedures, we obtain IA beams with beam waists of 0.95
$\upmu\mathrm{m}$. The intensity crosstalk measured by a high-dynamic range
camera is below the level of $10^{-3}$ on the neighboring ion sites, and the
Rabi-rate crosstalk of this DSA system is as low as $6\times 10^{-4}$. The
intensity crosstalk can be further suppressed by eliminating the aberrations
or employing a higher NA objective lens (up to 0.66 for our system), leading
to a much lower Rabi crosstalk.

Figure 2: (a) Schematic of optical setup of IA system. Only the trajectories
of multiple laser beams are shown, while the profiles of the laser beams are
not shown. (b) The relevant energy level of ${{}^{171}}\mathrm{Yb}^{+}$
qubits. The transition wavelength for Doppler cooling is 369.5 nm, which is
marked with purple lines. The energy levels associated with the SRT are marked
with green lines and the wavelength is 532 nm in the visible domain. (c)
Schematic of optical setup of the imaging system. The distances between each
optical element are not drawn as the true scale. The optic axis of the imaging
system is perpendicular to the plane of the Raman beam. HOAT: High optical
access trap, HWP: Half-wave plate, PBS: Polarization beam splitter, QWP:
Quarter-wave plate, AL: Achromatic lens, AOM: Acousto-optic modulator, AOD:
Acousto-optic deflector, RLB: Raman laser beam, OL: objective lens, BL:
Biconvex lens, MCPMT: Multi-channel photomultiplier.

The optical setup of our experimental apparatus is shown in Fig. 2(a). For our
experiment, the ${{}^{171}}\mathrm{Yb}^{+}$ qubits are confined in a glass
cell Paul trap with high optical access (HOAT) in multi-directions He  _et
al._ (2021). The Raman laser is a continuous wave (CW) laser modulated by a
fiber electro-optical modulator (FEOM) and its wavelength is 532 nm in the
visible domain Li  _et al._ (2022). This visible Raman laser is fiber (a
single mode and polarization-maintaining fiber) delivered to a beam splitter
comprised of polarizing elements and a mirror. Then the Raman laser is split
on the beam splitter into two beams and the power of each beam can be adjusted
optionally via a half-wave plate. The two portions of the Raman laser
propagate through two symmetrically configured optical paths and are then
recombined at the ion position, forming a Mach–Zehnder interferometer.
Symmetrical configuration can reduce undesirable crosstalk caused by
misalignment of the Raman beams from two sides.

To make the tightly focusing beams reach the diffraction limit, the beams are
magnified by lens groups, which consist of achromatic lenses to suppress the
aberration. Given that the optical paths on both sides are symmetrically
configured, we would like to describe the optical path on one side here. The
acousto-optic modulator (AOM) is placed between the first pair of lens groups.
By tuning the driving frequencies of AOMs, we can shift the relative frequency
of the Raman laser propagated on two sides to drive the carrier, blue-sideband
or red-sideband transition. Then the Raman laser passes through a total
reflection prism installed on a high-precision displacement stage (not shown
in the figure) to precisely tune the length of the optical path without
changing the orientation of the laser beam. The AOD is placed after the prism
and driven by an arbitrary waveform generator (AWG).

Figure 3: (a) Image of diffractive beam spots of AOD. The beams are focused on
the plane of ions. The arrows in the figure indicate the orientation of the
ion chain. The distance between beam spots is equal to the distance between
the ions. This picture is recorded by a high-dynamic range camera. (b) The
profiles of the corresponding beam spots. The result shows that the intensity
crosstalk is below $10^{-3}$ when the addressing beam is illuminating on the
nearest neighboring site.

By reprogramming the amplitude, frequency and phase of the driving signal
generated by AWG, the IA system allows for individual control of any ion qubit
in the chain.

Diffractive beams from AOD are reflected by a piezoelectric mirror and are
magnified by the second pair of lens groups. The piezoelectric mirror is
utilized to fine-tune the position of tightly focusing beams to align them
with the ions. The beams are finally demagnified by a 0.4 NA objective lens
and focused on the plane of ions. The half and quarter-wave plate is utilized
to rotate the polarization of the Raman laser.

The axis of the optical system for imaging is perpendicular to the plane of
the Raman beam, shown in Fig. 2(c). Another 0.4 NA objective lens is employed
to collect the fluorescence from the ion. The multi-channel photomultiplier
(MCPMT) has 32 segmented channels. The interval between each channel is 200
$\upmu\mathrm{m}$ while the width of each channel is 800 $\upmu\mathrm{m}$.
The optical system for imaging can magnify the image of ions by 200 to 1000
times, depending on the distance between the last achromatic lens and MCPMT
mounted on the moveable optical stage. For individual ion in the ion chain,
its fluorescence will be steering to the corresponding channel of MCPMT after
the magnification, enabling parallel state detection of multiple ions. The
function of apertures in the optical path is optical filtering, which improves
the imaging quality. Crosstalk from the imaging system is negligible in the
case of good imaging quality and high magnification.

To characterize the intensity crosstalk between multiple beams diffracted from
AOD, the laser beams focused on the ion plane are re-imaged onto a high-
dynamic range camera. The AOD generates a single diffractive beam spot by
applying a monochromatic drive signal. Then we sequentially increase the
driving signal frequency to obtain ten diffractive spots. The driving
frequencies array here equals those applied to the AOD when addressing
multiple ions. We then eliminate the background noise and stitch the images of
each diffractive beam spot together, shown in Fig. 3(a). The arrows in Fig.
3(a) indicate the orientation of the ion chain and define a cross-section
across the center of beam spots. Fig. 3(b) shows the profiles of IA beams in
this cross-section. From Fig. 3(b), the intensity crosstalk of the neighboring
site of ions is below the level of $10^{-3}$. We also notice the stray spots
distributed along the ion chain direction in the optical system, which is
highly likely to cause an increase in intensity crosstalk.

Figure 4: Scan of the addressing beam over a ion chain containing two ions.
The fluorescence scattered by ion A and B is collected by two individually
channels of MCPMT, respectively. Dashed lines are the Gaussian fits to data
from corresponding channels. The Rabi frequencies of ion A and B are slightly
different but they can be fine-tuned through changing the amplitude of the
driving signal of AODs.

The schematic diagram of energy level of ${{}^{171}}\mathrm{Yb}^{+}$ ion is
shown in Fig. 2(b). The qubits are defined by the two hyperfine ground states
of ${{}^{171}}\mathrm{Yb}^{+}$, namely,
$|0\rangle\equiv{{}^{2}}S_{1/2}|F=0,m_{F}=0\rangle$ and
$|1\rangle\equiv{{}^{2}}S_{1/2}|F=1,m_{F}=0\rangle$. Doppler cooling,
electromagnetic induction transparency (EIT) cooling Feng  _et al._ (2020);
Qiao  _et al._ (2021), optical pumping and state-dependent fluorescence
detection are implemented via 369.5 nm CW laser while the repump from
${{}^{2}}D_{3/2}$ state is implemented via 935 nm CW laser Olmschenk  _et al._
(2007).

Figure 5: The Rabi rate crosstalk between ion A and B. (a)$\&$(b) The ion A
(B) is illuminated by the Raman beam within short period to measure the
accurate Rabi frequency (or the period of Rabi cycle). When ion A (B) is
addressed, the period of Rabi cycle is 22.84 us (30.85 us). (c)$\&$(d) The ion
A (B) is illuminated by the Raman beam within long period to measure the
crosstalk caused by addressing beam. When ion A (B) is addressed, ion B (A)
will oscillate between $|1\rangle$ and $|0\rangle$ state with a very low Rabi
frequency due to the crosstalk.

The qubits are driven by the SRT using two coherent laser fields of which the
frequency gap between them equals the frequency splitting of the hyperfine
ground states, namely, $\omega_{\mathrm{HF}}=$12.64 GHz. When an ion is
illuminated by the Raman laser beams, its quantum state will flip between
$|0\rangle$ and $|1\rangle$. The Raman laser system used in this experiment
has been detailedly described in our previous work Li  _et al._ (2022).

In the AOD-based IA system, each specific diffractive beam of AODs matches its
corresponding addressed ion. The center frequency of AOD (DTSX-400, AA Opto
Electronic) is 77.5 MHz and the -1st order diffractive beam’s angular
direction is colinear with the incident laser beam. In preparation for the
following experiments, we confine a single ion in the trap and set the drive
frequency of AODs equal to the AODs’ center frequency. Taking advantage of
single Raman beam flipping, we align the diffractive beam of both sides to the
single ion by adjusting the voltage applied to the piezoelectric mirrors Li
_et al._ (2022).

For the first experiment, we confine two ions in the trap and sweep the drive
frequency of AOD to investigate the corresponding frequency of each ion. For
convenience of description, we denote the two ions in the experiment as ion A
and B, respectively. With the sweep of the drive frequency, the orientation of
movement of the diffractive beam is parallel with the orientation of the ion
chain. When an ion is overlapped with the scanning beam, its quantum state
will oscillate between $|1\rangle$ and $|0\rangle$ state due to the SRT
process. Depending on the overlap between the ion and the diffractive beam,
the ion experiences different laser intensity, resulting in a different Rabi
frequency. In the case of that the interaction time is fixed, the excitation
probability of $|1\rangle$ state is positively proportional to the intensity
experienced by the ion. Hence, the relationship between the drive frequency of
AODs and the excitation probability is expected to be the intensity
distribution curve of the diffractive. Fig. 4 shows a Gaussian distribution
curve.

We utilize the function

$A\left(x\right)=\sum_{i}{A_{i}e^{\left[-2\left(\left(x-x_{i}\right)/w_{i}\right)^{2}\right]}}$
(1)

to fit the data, where $A\left(i\right)$ is the peak intensity of each beam
spot, $x_{i}$ is the center position of each beam spot, $w_{i}$ is the waist
of each beam spot. The fit curve indicates that the profile of the addressing
beam is Gaussian distribution but with stray spots, which is consistent with
the result shown in Fig. 3. The ion A and B can be addressed by setting the
drive frequency of AODs as 77.07 MHz and 78.03 MHz respectively, according to
the fitting result in Fig. 4. Given that the distance between ion A and B is
approximately 5.5 $\upmu\mathrm{m}$, the beam waist is approximately 0.95
$\upmu\mathrm{m}$, which doesn’t reach the diffraction limit (0.8113
$\upmu\mathrm{m}$) of the optical system. For every date point we obtain in
Fig. 4, the complete experimental sequence consists of Doppler cooling laser
with 1 ms, EIT cooling lasers with 1 ms, pumping laser with 20 us, Raman
addressing laser with a fixed interaction time and detecting laser with 1 ms.
The sequence is repeated 100 times to obtain the probability.

For the further experiment, we investigate the Rabi rate crosstalk, which is
defined by the ratio of Rabi frequency of adjacent ion and addressed ion,
namely, $\varOmega_{\mathrm{adjacent}}/\varOmega_{\mathrm{addressed}}$. To
characterize the accurate Rabi frequency in the case of that the ion A or B is
individually addressed, we set the drive frequency of AODs as the
corresponding frequency of ion A or B and keep the time interval of
interaction time short. For instance, by setting the drive frequency of AODs
as 77.07 MHz, the Raman beam will only focus on the position of ion A, and
vice versa. Both ion A and B are cooled to the ground motional state and
initialized to the $|0\rangle$ state by optical pumping. Then ion A is
illuminated by the Raman beam during a varying interaction time, leaving ion B
being affected by the addressing beam of ion A. The final states of ion A and
B versus interaction time are recorded as Fig. 5(a) and (b). The Rabi
frequency of ion A or B when they are individually addressed is $2\pi\times
43.78$ kHz or $2\pi\times 32.42$ kHz respectively.

To characterize the crosstalk of the Rabi rate, we implement a similar
experiment but we extend the maximum interaction time to 5 ms. The results are
shown in Fig. 5(c) and (b). It is worth noting that the intensity of the Raman
laser hasn’t been stabilized during our experiment, resulting in the decay of
Rabi flopping. When ion A is illuminated by its addressing beam, ion B
transitions from $|0\rangle$ state to $|1\rangle$ state with a very low Rabi
frequency, and vice versa. The final excitation probability of $|1\rangle$
state of unaddressed ion A and B after 5 ms is 0.38 and 0.11, respectively.
The expression of excitation probability of $|1\rangle$ state is

$P_{i,|1\rangle}\left(t\right)=\frac{1}{2}\left[1-\cos\left(\varOmega_{i}t\right)\right],$
(2)

where $i$ denotes ion A or B. With Eq. 2, the crosstalk of Rabi rate for ion A
and B is calculated as $6.32\times 10^{-4}$ and $9.66\times 10^{-4}$,
respectively.

In this work, we demonstrate a low-crosstalk individual optical addressing
system using a pair of AODs and investigate that the crosstalk of Rabi rate on
the site of neighboring ion is on the level at $6\times 10^{-4}$ to $9\times
10^{-4}$. Taking advantage of the high optical access of our glass cell trap,
we employ two 0.4 NA objective lenses on the double side of Raman laser and
obtain a beam waist of 0.95 $\upmu\mathrm{m}$ at the ion location. The
crosstalk of our double-side addressing system is substantially lower than
those of single-side addressing systems. The intensity crosstalk of our DSA
system is below $10^{-3}$ and it is mainly caused by undesirable optical
aberrations. Higher NA objective lenses (up to 0.66 NA for our system) and
aberration elimination techniques can further reduce the crosstalk. By
optionally reprogramming the frequency, amplitude and phase of the driving
signal of AODs, the DSA system can offer an agile, flexible and low error-rate
manipulation of atomic ion qubits. This low error-rate DSA system enables the
implementation of high-fidelity two-qubit gates in a scalable trapped ion
system and supports the realization of the error-corrected quantum computer.

###### Acknowledgements.

This work was supported by the National Key Research and Development Program
of China (No. 2017YFA0304100), National Natural Science Foundation of China
(Nos. 11774335, 11734015), the Key Research Program of Frontier Sciences, CAS
(No. QYZDY-SSWSLH003), Innovation Program for Quantum Science and Technology
(Nos. 2021ZD0301604, 2021ZD0301200). This work was partially carried out at
the USTC Center for Micro and Nanoscale Research and Fabrication.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Wang  _et al._ (2017) Y. Wang, M. Um, J. Zhang, S. An, M. Lyu, J.-N. Zhang, L.-M. Duan, D. Yum, and K. Kim, “Single-qubit quantum memory exceeding ten-minute coherence time,” Nat. Photonics 11, 646–650 (2017).
  * Wang  _et al._ (2021) P. Wang, C.-Y. Luan, M. Qiao, M. Um, J. Zhang, Y. Wang, X. Yuan, M. Gu, J. Zhang, and K. JKim, “Single ion qubit with estimated coherence time exceeding one hour,” Nat. Commun. 12, 233 (2021).
  * An  _et al._ (2022) F. A. An, A. Ransford, A. Schaffer, L. R. Sletten, J. Gaebler, J. Hostetter, and G. Vittorini, “High fidelity state preparation and measurement of ion hyperfine qubits with $i>\frac{1}{2}$,” Phys. Rev. Lett. 129, 130501 (2022).
  * Harty  _et al._ (2014) T. P. Harty, D. T. C. Allcock, C. J. Ballance, L. Guidoni, H. A. Janacek, N. M. Linke, D. N. Stacey, and D. M. Lucas, “High-fidelity preparation, gates, memory, and readout of a trapped-ion quantum bit,” Phys. Rev. Lett. 113, 220501 (2014).
  * Gaebler  _et al._ (2016) J. P. Gaebler, T. R. Tan, Y. Lin, Y. Wan, R. Bowler, A. C. Keith, S. Glancy, K. Coakley, E. Knill, D. Leibfried, and D. J. Wineland, “High-fidelity universal gate set for ${{}^{9}\mathrm{Be}}^{+}$ ion qubits,” Phys. Rev. Lett. 117, 060505 (2016).
  * Ballance  _et al._ (2016) C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, and D. M. Lucas, “High-fidelity quantum logic gates using trapped-ion hyperfine qubits,” Phys. Rev. Lett. 117, 060504 (2016).
  * Clark  _et al._ (2021) C. R. Clark, H. N. Tinkey, B. C. Sawyer, A. M. Meier, K. A. Burkhardt, C. M. Seck, C. M. Shappert, N. D. Guise, C. E. Volin, S. D. Fallek, H. T. Hayden, W. G. Rellergert, and K. R. Brown, “High-fidelity bell-state preparation with ${}^{40}{\mathrm{ca}}^{+}$ optical qubits,” Phys. Rev. Lett. 127, 130505 (2021).
  * Figgatt  _et al._ (2019) C. Figgatt, A. Ostrander, N. M. Linke, K. A. Landsman, D. Zhu, D. Maslov, and C. Monroe, “Parallel entangling operations on a universal ion-trap quantum computer,” Nature 572, 368–372 (2019).
  * Figgatt  _et al._ (2017) C. Figgatt, D. Maslov, K. A. Landsman, N. M. Linke, S. Debnath, and C. Monroe, “Complete 3-qubit grover search on a programmable quantum computer,” Nat. Commun. 8, 1918 (2017).
  * Smith  _et al._ (2016) J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, “Many-body localization in a quantum simulator with programmable random disorder,” Nat. Phys. 12, 907–911 (2016).
  * Martinez  _et al._ (2016) E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt, “Real-time dynamics of lattice gauge theories with a few-qubit quantum computer,” Nature 534, 516–519 (2016).
  * Nguyen  _et al._ (2021) N. H. Nguyen, M. Li, A. M. Green, C. Huerta Alderete, Y. Zhu, D. Zhu, K. R. Brown, and N. M. Linke, “Demonstration of shor encoding on a trapped-ion quantum computer,” Phys. Rev. Appl. 16, 024057 (2021).
  * Schindler  _et al._ (2011) P. Schindler, J. T. Barreiro, T. Monz, V. Nebendahl, D. Nigg, M. Chwalla, M. Hennrich, and R. Blatt, “Experimental repetitive quantum error correction,” Science 332, 1059–1061 (2011), https://www.science.org/doi/pdf/10.1126/science.1203329 .
  * Erhard  _et al._ (2021) A. Erhard, H. Poulsen Nautrup, M. Meth, L. Postler, R. Stricker, M. Stadler, V. Negnevitsky, M. Ringbauer, P. Schindler, H. J. Briegel, R. Blatt, N. Friis, and T. Monz, “Entangling logical qubits with lattice surgery,” Nature 589, 220–224 (2021).
  * Egan  _et al._ (2021) L. Egan, D. M. Debroy, C. Noel, A. Risinger, D. Zhu, D. Biswas, M. Newman, M. Li, K. R. Brown, M. Cetina, and C. Monroe, “Fault-tolerant control of an error-corrected qubit,” Nature 598, 281–286 (2021).
  * Schmidt-Kaler  _et al._ (2003) F. Schmidt-Kaler, H. Häffner, M. Riebe, S. Gulde, G. P. T. Lancaster, T. Deuschle, C. Becher, C. F. Roos, J. Eschner, and R. Blatt, “Realization of the cirac–zoller controlled-not quantum gate,” Nature 422, 408–411 (2003).
  * Parrado-Rodríguez  _et al._ (2021) P. Parrado-Rodríguez, C. Ryan-Anderson, A. Bermudez, and M. Müller, “Crosstalk Suppression for Fault-tolerant Quantum Error Correction with Trapped Ions,” Quantum 5, 487 (2021).
  * Crain  _et al._ (2014) S. Crain, E. Mount, S. Baek, and J. Kim, “Individual addressing of trapped 171Yb+ ion qubits using a microelectromechanical systems-based beam steering system,” Appl. Phys. Lett. 105 (2014).
  * Egan (2021) L. N. Egan,  _Scaling Quantum Computers With Long Chains Of Trapped Ions_ , Ph.D. thesis, University of Maryland (2021).
  * Binai-Motlagh  _et al._ (2023) A. Binai-Motlagh, M. Day, N. Videnov, N. Greenberg, C. Senko, and R. Islam, “A guided light system for agile individual addressing of ba+ qubits with $10^{-4}$ level intensity crosstalk,” (2023), arXiv:2302.14711 [quant-ph] .
  * Pogorelov  _et al._ (2021) I. Pogorelov, T. Feldker, C. D. Marciniak, L. Postler, G. Jacob, O. Krieglsteiner, V. Podlesnic, M. Meth, V. Negnevitsky, M. Stadler, B. Höfer, C. Wächter, K. Lakhmanskiy, R. Blatt, P. Schindler, and T. Monz, “Compact ion-trap quantum computing demonstrator,” PRX Quantum 2, 020343 (2021).
  * Wang  _et al._ (2020) Y. Wang, S. Crain, C. Fang, B. Zhang, S. Huang, Q. Liang, P. H. Leung, K. R. Brown, and J. Kim, “High-fidelity two-qubit gates using a microelectromechanical-system-based beam steering system for individual qubit addressing,” Phys. Rev. Lett. 125, 150505 (2020).
  * Bergmann, Theuer, and Shore (1998) K. Bergmann, H. Theuer, and B. W. Shore, “Coherent population transfer among quantum states of atoms and molecules,” Rev. Mod. Phys. 70, 1003–1025 (1998).
  * Vitanov  _et al._ (2017) N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, “Stimulated raman adiabatic passage in physics, chemistry, and beyond,” Rev. Mod. Phys. 89, 015006 (2017).
  * Sørensen  _et al._ (2006) J. L. Sørensen, D. Møller, T. Iversen, J. B. Thomsen, F. Jensen, P. Staanum, D. Voigt, and M. Drewsen, “Efficient coherent internal state transfer in trapped ions using stimulated raman adiabatic passage,” New J. Phys. 8, 261 (2006).
  * He  _et al._ (2021) R. He, J.-M. Cui, R.-R. Li, Z.-H. Qian, Y. Chen, M.-Z. Ai, Y.-F. Huang, C.-F. Li, and G.-C. Guo, “An ion trap apparatus with high optical access in multiple directions,” Rev. Sci. Instrum. 92 (2021), 10.1063/5.0043985, 073201, https://pubs.aip.org/aip/rsi/article-pdf/doi/10.1063/5.0043985/16016833/073201_1_online.pdf .
  * Li  _et al._ (2022) R.-R. Li, R. He, J.-M. Cui, Y. Chen, W.-R. Ye, Y.-L. Chen, Y.-F. Huang, C.-F. Li, and G.-C. Guo, “A versatile multi-tone laser system for manipulating atomic qubits based on a fiber mach–zehnder modulator and second harmonic generation,” Opt. Express 30, 30098–30107 (2022).
  * Feng  _et al._ (2020) L. Feng, W. L. Tan, A. De, A. Menon, A. Chu, G. Pagano, and C. Monroe, “Efficient ground-state cooling of large trapped-ion chains with an electromagnetically-induced-transparency tripod scheme,” Phys. Rev. Lett. 125, 053001 (2020).
  * Qiao  _et al._ (2021) M. Qiao, Y. Wang, Z. Cai, B. Du, P. Wang, C. Luan, W. Chen, H.-R. Noh, and K. Kim, “Double-electromagnetically-induced-transparency ground-state cooling of stationary two-dimensional ion crystals,” Phys. Rev. Lett. 126, 023604 (2021).
  * Olmschenk  _et al._ (2007) S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz, and C. Monroe, “Manipulation and detection of a trapped ${\mathrm{yb}}^{+}$ hyperfine qubit,” Phys. Rev. A 76, 052314 (2007).


# Implications of comprehensive nuclear and astrophysics data
on the equations of state of neutron star matter

Sk Md Adil Imam1,2<EMAIL_ADDRESS>Tuhin Malik3<EMAIL_ADDRESS>Constança
Providência3<EMAIL_ADDRESS>B. K. Agrawal1,2<EMAIL_ADDRESS>1Saha
Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata 700064, India.  2Homi
Bhabha National Institute, Anushakti Nagar, Mumbai 400094, India.  3CFisUC,
Department of Physics, University of Coimbra, 3004-516 Coimbra, Portugal

###### Abstract

The equations of state (EoSs) governing neutron star (NS) matter obtained for
both non-relativistic and relativistic mean-field models are systematically
confronted with a diverse set of terrestrial data and astrophysical
observations within the Bayesian framework. The terrestrial data, spans from
bulk properties of finite nuclei to the heavy-ion collisions, constrain the
symmetric nuclear matter EoS and the symmetry energy up to twice the
saturation density ($\rho_{0}$= 0.16 fm-3). The astrophysical observations
encompass the NS radius, the tidal deformability, and the lower bound on
maximum mass. Three distinct posterior distributions of EoSs are generated by
gradually updating the priors with different constraints: (i) only the maximum
NS mass, (ii) incorporating additional terrestrial data, (iii) combining both
the terrestrial data and astrophysical observations. These distributions are
then compared using the Kullback-Liebler divergence which highlights the
significant constraints imposed on the EoSs by the currently available lower
bound of NS maximum mass and terrestrial data. The remaining astrophysical
observations marginally refine the EoS within the density range $\sim$
2-3$\rho_{0}$. It is observed that the non-relativistic mean field model
yields stiffer EoS for the symmetric nuclear matter while softer density
dependent symmetry energy compared to the relativistic one when all the
data/observations are considered.

##  I Introduction

Determining the unified equation of state (EoS) accurately requires inputs
from both nuclear physics experiments and astrophysical observations. The
experimental data on heavy-ion collisions (HIC) and bulk properties of finite
nuclei confine the behavior of the EoS at low densities ($\rho\leq
2\rho_{0}$). These data specifically constrain the EoS for symmetric nuclear
matter (SNM) and the density-dependent symmetry energy, crucial elements for
determining the EoS of neutron star matter Lattimer and Prakash (2001);
Lattimer (2021); Fattoyev  _et al._ (2010); Tsang  _et al._ (2012); Fortin
_et al._ (2016); Margueron  _et al._ (2018); Fattoyev  _et al._ (2018); Greif
_et al._ (2020); Reed  _et al._ (2023).The astrophysical observations such as
neutron star radius, tidal deformability, and maximum mass further constrain
the EoS at high densities Raaijmakers  _et al._ (2021); Annala  _et al._
(2022); Ghosh  _et al._ (2022); Huang  _et al._ (2023); Char  _et al._ (2020);
Most  _et al._ (2018).The radius and tidal deformability parameter of neutron
stars are sensitive to the EoS. They have been deduced from the gravitational
wave events GW170817 and GW190425 observed using the Advanced-LIGO Aasi  _et
al._ (2015) and Advanced-Virgo detectors Acernese  _et al._ (2015). In
particular, the gravitational waves associated to the event GW170817 stemmed
from a binary neutron star (BNS) merger with component masses ranging from
1.17 to 1.6 M⊙ Abbott et al. (2019a, b) has sparked numerous theoretical
investigations into the neutron star (NS) properties Abbott et al. (2017,
2018); De  _et al._ (2018); Fattoyev  _et al._ (2018); Landry and Essick
(2019); Piekarewicz and Fattoyev (2019); Abbott et al. (2020); Dietrich  _et
al._ (2020). Future observations of coalescing BNS events by detectors like
LIGO-Virgo-KAGRA, Einstein Telescope Punturo  _et al._ (2010), and Cosmic
Explorer Reitze  _et al._ (2019) are likely to occur more frequently, enabling
a more precise determination of the EoS Iacovelli  _et al._ (2023); Huxford
_et al._ (2023); Rose  _et al._ (2023); Walker  _et al._ (2024). The Neutron
star Interior Composition Explorer (NICER) also offers complementary
information on the NS properties Watts et al. (2016); Psaltis  _et al._
(2014).

The nuclear matter parameters are the key quantities that determine the EoS.
The combined constraints from the bulk properties of finite nuclei and NS
observable has been employed to constrain the nuclear matter parameters which
are the expansion coefficients in the Taylor model for the EoS Zhang  _et al._
(2018). The finite nuclei constraints are encoded through a few low order
nuclear matter parameters treating them independent of each other. The NS
observable considered are the radius, 10.62 km $<R_{1.4}<$ 12.83 km Lattimer
and Steiner (2014a), the dimensionless tidal deformability
$\Lambda_{1.4}\leq$800 Abbott et al. (2017); Abbott  _et al._ (2017a) and
lower limit of the maximum mass imposed by the NS PSR J0348 \- 0432 with a
mass 2.01 $\pm$ 0.04 M⊙ Antoniadis  _et al._ (2013); Lattimer and Steiner
(2014b). It is demonstrated that these constraints limit the multidimensional
space of few higher order nuclear matter parameters. It may be emphasized that
the astrophysical inputs employed in Ref.Zhang  _et al._ (2018) have undergone
significant revision. Two different groups of NICER have reported neutron
star’s mass and radius simultaneously for PSR J0030+0451 with radius R
$=13.02^{+1.24}_{-1.06}$km for mass M $=1.44^{+0.15}_{-0.14}$ M⊙ Miller et al.
(2019) and R $=12.71^{+1.14}_{-1.19}$ km for M $=1.34^{+0.15}_{-0.16}$ M⊙
Riley et al. (2019). Recently, another (heavier) pulsar PSR J0740+6620, R
$=13.7^{+2.6}_{-1.5}$ km with M $=2.08\pm 0.07$ M⊙ Miller  _et al._ (2021) and
R $=12.39^{+1.30}_{-0.98}$ km with M $=2.072^{+0.067}_{-0.066}$ M⊙ Riley et
al. (2021) have been reported. The improvement in the precision of
astrophysical observation and the terrestrial data has triggered the search
for the universal EoS .

A Bayesian inference technique has been employed Huth  _et al._ (2022) to
combine data from astrophysical multi-messenger observations of neutron stars
Abbott et al. (2017); Abbott  _et al._ (2020, 2017b); Coughlin  _et al._
(2017); Miller et al. (2019); Riley et al. (2019); Miller  _et al._ (2021);
Riley et al. (2021); Dietrich  _et al._ (2020), high-energy heavy-ion
collisions involving gold nuclei Le Fèvre  _et al._ (2016); Russotto and et.
al. (2016) and microscopic nuclear theory calculations Hebeler  _et al._
(2013); Tews  _et al._ (2013); Lynn  _et al._ (2016); Drischler  _et al._
(2019, 2020); Huth  _et al._ (2021) to constrain the EoS of dense matter. The
microscopic input was extended up to 1.5 times the nuclear saturation density
($\rho_{0}$) for the $\beta$-EoS derived from chiral effective field theory
($\chi$EFT). The nuclear experimental data includes the symmetry energy and
symmetric nuclear matter pressure extracted from the FOPI Le Fèvre  _et al._
(2016) and the ASY-EoS Russotto and et. al. (2016) experiments, respectively.
The astrophysical data included mass measurements of massive neutron stars PSR
J0348+0432 Antoniadis  _et al._ (2013) and PSR J1614-2230 Arzoumanian  _et
al._ (2018). Constraints on the maximum mass of neutron stars were also
derived from the binary neutron-star collision GW170817 Margalit and Metzger
(2017); Rezzolla  _et al._ (2018). The NS radii of PSR J0030+0451 and PSR
J0740+6620 utilizing NICER and XMM-Newton data Miller et al. (2019); Miller
_et al._ (2021); Riley et al. (2021) was also included. A Bayesian inference
technique was employed to analyze the GW information from the GW170817 Abbott
et al. (2017) and GW190425 Abbott  _et al._ (2020) events, aligning observed
GW data with theoretical models reliant on neutron-star properties. The
analysis relied on a GW model Dietrich  _et al._ (2019), an enhanced version
of the primary waveform model utilized by the LIGO/Virgo Collaboration in the
exploration of GW170817 Abbott et al. (2019a) and GW190425 Abbott  _et al._
(2020). Information pertaining to the kilonova AT2017gfo Abbott  _et al._
(2017b), linked with the GW signal, was also incorporated. Electromagnetic
observations underwent comprehensive radiative transfer simulations Bulla
(2019) to extract insights from observed light curves and spectra Coughlin
_et al._ (2017). However, the EoS employed in this analysis extended the
$\chi$EFT calculation to $\beta$-equilibrium with a crust Tews (2017), while
beyond this density range, a collection of six randomly distributed points
within the speed of sound (cs) plane at baryon densities spanning 1.5 to
12$\rho_{0}$ were sampled, enforcing 0 $\leq$c${}_{s}\leq$ 1 at each point.

Of late Tsang  _et al._ (2023), the EoS has been constrained using a diverse
set of data from nuclear experiments and astrophysical observations. The
nuclear inputs consist of experimental data on bulk properties of finite
nuclei and HIC and astrophysical observations on NS properties. The properties
of finite nuclei and HIC are encoded through the symmetry energy, pressure for
symmetry energy and SNM spanning the density up to 2$\rho_{0}$. The
experimental data on nuclear masses, iso-vector giant dipole resonance,
neutron skin thickness and isobaric analog states constrain the symmetry
energy and the symmetry energy pressure at sub-saturation densities Brown and
Schwenk (2014); Lynch and Tsang (2022); Kortelainen  _et al._ (2012);
Danielewicz  _et al._ (2017); Zhang and Chen (2015); Adhikari  _et al._
(2021); Reed  _et al._ (2021). The HIC data from isospin diffusion, neutron-
proton ratios, charged pion spectra are incorporated to constrain the
behaviour of SNM and symmetry energy from sub-saturation to supra-saturation
densities Tsang  _et al._ (2009); Lynch and Tsang (2022); Morfouace  _et al._
(2019); Estee  _et al._ (2021); Cozma (2018); Russotto  _et al._ (2011);
Russotto and et. al. (2016); Danielewicz  _et al._ (2002); Le Fèvre  _et al._
(2016). The recent astrophysical observations for the radius of a NS with mass
$\sim$ 1.4 and 2.1 M⊙ Riley et al. (2019); Miller et al. (2019); Antoniadis
_et al._ (2013) and the dimensionless tidal deformability for 1.4 M⊙ NSAbbott
et al. (2018) have also been considered. These data are used in a Bayesian
framework to constrain the EoS obtained using Taylor expansion.

We use a Bayesian framework for the systematic analysis of the implications of
the diverse set of data Tsang  _et al._ (2023) in constraining the EoS of NS
matter obtained within a Skyrme based non-relativistic mean field and non-
linear variant of relativistic mean field (RMF) models. Our analysis involves
three distinct sets of calculations
(i) Priors updated with the NS maximum mass, $M_{\rm max}\geq$ 2.1 M⊙ .
(ii) Further inclusion of 12 experimental data.
(iii) Inclusion of the experimental data as well as astrophysical observations
corresponding to the NS properties.

The lower bound on maximum mass constraint considered was chosen taking into
account: the mass of the pulsar PSR J0740+6620 with M=2.08$\pm
0.07\,M_{\odot}$ Fonseca  _et al._ (2021); the mass of the pulsar PSR
J0952-0607, M = 2.35 $\pm 0.14\,M_{\odot}$ Romani  _et al._ (2022), which is
an analysis together with information from other black widow and red back
pulsars gave a minimum value for the maximum NS mass of 2.19$M_{\odot}$ (2.09
$M_{\odot}$) at 1$\sigma$ (3$\sigma$) confidence Romani  _et al._ (2022). In
the present analysis we will take as the lower bound on maximum NS mass
2.1$M_{\odot}$, compatible with these two results, but will also discuss the
effect of taking as the lower bound on maximum mass 2.17$M_{\odot}$, i.e. the
mass of pulsar PSR J0952-0607 at 1$\sigma$.

This paper is organised as follows. In Sec.II a brief description of mean
field models and Bayesian framework are presented. Sec.III summarizes the
experimental data and astrophysical observations considered. The results are
discussed in Sec.IV. Finally, the conclusions are drawn in Sec.V.

##  II Methodology

We use non-relativistic and the relativistic mean field models for the EoSs.
The non-relativistic mean field models are based on Skyrme type effective
interaction. The relativistic mean field models are derived from the effective
Lagrangian which describes the interactions of nucleons through scalar-
isoscalar meson $\sigma$, vector-isoscalar meson $\omega$ and vector-isovector
meson $\varrho$ which also includes non-linear self interaction terms for
$\sigma$ mesons and mixed interaction term involving $\omega$ and $\rho$
mesons. Once we have these EoSs we can calculate the energy of the SNM,
e${}_{SNM}(\rho)$ and symmetry energy, e${}_{\rm sym}(\rho)$ and use them to
obtain the pressure of the symmetric nuclear matter (PSNM) and the symmetry
energy (Psym). Empirical values of these quantities are deduced from
experimental data on finite nuclei and HIC which are employed to constrain the
required EoS.

###  II.1 SKYRME EoS

The Hamiltonian density for the Skyrme type effective interaction contributing
to uniform nuclear matter can be written as

${\cal H}={\cal K}+{\cal H}_{0}+{\cal H}_{3}+{\cal H}_{\rm eff}$ (1)

where, ${\cal K}=\frac{\hbar^{2}}{2m}\tau$ is the kinetic energy term, ${\cal
H}_{0}$ is the zero-range term, ${\cal H}_{3}$ the density dependent term,
${\cal H}_{\rm eff}$ an effective-mass term. For the Skyrme interaction we
have,

${\cal
H}_{0}=\frac{1}{4}t_{0}\left[(2+x_{0})\rho^{2}-(2x_{0}+1)(\rho_{p}^{2}+\rho_{n}^{2})\right],$
(2) ${\cal
H}_{3}=\frac{1}{24}t_{3}\rho^{\sigma}\left[(2+x_{3})\rho^{2}-(2x_{3}+1)(\rho_{p}^{2}+\rho_{n}^{2})\right],$
(3) $\displaystyle{\cal H}_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{1}{8}\left[t_{1}(2+x_{1})+t_{2}(2+x_{2})\right]\tau\rho$
(4) $\displaystyle+$
$\displaystyle\frac{1}{8}\left[t_{2}(2x_{2}+1)-t_{1}(2x_{1}+1)\right](\tau_{p}\rho_{p}+\tau_{n}\rho_{n}),$

Here, $\rho=\rho_{p}+\rho_{n},$ and $\tau=\tau_{p}+\tau_{n},$ are the particle
number density and kinetic energy density with $p$ and $n$ denoting the
protons and neutrons, respectively.

We can express the Skyrme parameters in terms of the nuclear matter parameters
as Chabanat  _et al._ (1997); Margueron  _et al._ (2002); Gómez  _et al._
(1992); Agrawal  _et al._ (2005)

$\displaystyle t_{0}$ $\displaystyle=$
$\displaystyle(\frac{8}{\rho_{0}})*(\frac{f_{1}}{f_{2}}+f_{3})$ (5)
$\displaystyle\sigma$ $\displaystyle=$ $\displaystyle\frac{f_{5}}{f_{6}}$ (6)
$\displaystyle t_{3}$ $\displaystyle=$
$\displaystyle(\frac{16}{\rho_{0}^{\sigma+1}})(\frac{f_{4}}{f_{2}})$ (7)
$\displaystyle\theta_{SNM}$ $\displaystyle=$
$\displaystyle\frac{1}{8}\left(3t_{1}+(5+4x_{2})t_{2}\right)$ (8)
$\displaystyle=$ $\displaystyle\frac{\hbar^{2}}{m\rho_{0}}(\frac{m}{m^{*}}-1)$

Where,

$\displaystyle f_{1}$ $\displaystyle=$
$\displaystyle(e_{0}+(\frac{2m}{m^{*}}-3)(\frac{\hbar^{2}}{10m})k_{f}^{2})$
$\displaystyle(\frac{K_{0}}{27}-(1-\frac{6m^{*}}{5m})(\frac{\hbar^{2}}{9m^{*}})k_{f}^{2})$
$\displaystyle f_{2}$ $\displaystyle=$ $\displaystyle
e_{0}+\frac{K_{0}}{9}-(\frac{4m^{*}}{3m}-1)(\frac{\hbar^{2}}{10m})k_{f}^{2}$
$\displaystyle f_{3}$ $\displaystyle=$
$\displaystyle(1-\frac{5m}{3m^{*}})(\frac{\hbar^{2}}{10m})k_{f}^{2}$
$\displaystyle f_{4}$ $\displaystyle=$
$\displaystyle(e_{0}+(\frac{2m^{*}}{m}-3)(\frac{\hbar^{2}}{10m})k_{f}^{2})^{2}$
$\displaystyle f_{5}$ $\displaystyle=$ $\displaystyle-
e_{0}-\frac{K_{0}}{9}+(\frac{4m^{*}}{3m}-1)(\frac{\hbar^{2}}{10m})k_{f}^{2}$
$\displaystyle f_{6}$ $\displaystyle=$ $\displaystyle
e_{0}+(\frac{2m^{*}}{m}-3)(\frac{\hbar^{2}}{10m})k_{f}^{2}$

Similarly, we can also express a few of the Skyrme parameters that determine
the symmetry energy in a linear combination of the nuclear matter parameters
J0, L0 and Ksym0 of the symmetry energy with help of the Skyrme parameters of
SNM as

$\displaystyle\left(\begin{matrix}x_{0}\\\ \theta_{\rm sym}\\\ x_{3}\\\
\end{matrix}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{matrix}a_{01}&a_{11}&a_{21}\\\
a_{02}&a_{12}&a_{22}\\\ a_{03}&a_{13}&a_{23}\\\
\end{matrix}\right)^{-1}\left(\begin{matrix}J_{\rm 01}\\\ L_{\rm 01}\\\ K_{\rm
sym01}\\\ \end{matrix}\right),$ (9)

Where,

$\displaystyle\theta_{sym}$ $\displaystyle=$ $\displaystyle
3t_{1}x_{1}-t_{2}(4+5x_{2})$ $\displaystyle J_{\rm 01}$ $\displaystyle=$
$\displaystyle
J_{0}-a_{s}\rho_{0}^{\frac{2}{3}}+\frac{t_{0}\rho_{0}}{8}+\frac{t_{3}}{48}\rho_{0}^{\sigma+1}$
$\displaystyle L_{\rm 01}$ $\displaystyle=$
$\displaystyle\frac{L_{0}}{3\rho_{0}}-\frac{2}{3}a_{s}\rho_{0}^{-\frac{1}{3}}+\frac{t_{0}}{8}+\frac{t_{3}}{48}(\sigma+1)\rho_{0}^{\sigma}$
$\displaystyle K_{\rm sym01}$ $\displaystyle=$ $\displaystyle\frac{K_{\rm
sym0}}{3\rho_{0}^{2}}+\frac{2}{9}a_{s}\rho_{0}^{-\frac{4}{3}}+\frac{t_{3}}{48}\sigma(\sigma+1)\rho_{0}^{\sigma-1}$
$\displaystyle a_{01}$ $\displaystyle=$
$\displaystyle-\frac{t_{0}\rho_{0}}{4}$ $\displaystyle a_{11}$
$\displaystyle=$
$\displaystyle-\frac{1}{24}(\frac{3\pi^{2}}{2})^{2/3}{\rho_{0}}^{5/3}$
$\displaystyle a_{21}$ $\displaystyle=$
$\displaystyle-\frac{t_{3}}{24}\rho_{0}^{\sigma+1}$ $\displaystyle a_{02}$
$\displaystyle=$ $\displaystyle-\frac{t0}{4}$ $\displaystyle a_{12}$
$\displaystyle=$
$\displaystyle-\frac{5}{72}(\frac{3\pi^{2}}{2})^{2/3}\rho_{0}^{2/3}$
$\displaystyle a_{22}$ $\displaystyle=$
$\displaystyle-(\sigma+1)\frac{t_{3}}{24}\rho_{0}^{\sigma}$ $\displaystyle
a_{03}$ $\displaystyle=$ $\displaystyle 0$ $\displaystyle a_{13}$
$\displaystyle=$
$\displaystyle-\frac{5}{108}(\frac{3\pi^{2}}{2})^{2/3}\rho_{0}^{-1/3}$
$\displaystyle a_{23}$ $\displaystyle=$
$\displaystyle-\sigma(\sigma+1)\frac{t_{3}}{24}\rho_{0}^{\sigma-1},$

with as = $\frac{\hbar^{2}}{6m}(\frac{3\pi^{2}}{2})^{2/3}$.

Once these Skyrme parameters are known the energy density for a given
$\rho_{n}$ and $\rho_{p}$ can be computed Chabanat  _et al._ (1997).

Table 1: The empirical values of symmetry energy (esym), symmetry energy
pressure (Psym) and symmetric nuclear matter pressure (PSNM) from experimental
data on the bulk properties of finite nuclei and HIC. The astrophysical
observational constraints of radii and tidal deformability of neutron star.
See RefTsang  _et al._ (2023) for details.

Symmetric matter |  |  |  |
---|---|---|---|---
Constraints |  $n$ (fm-3) | P${}_{\text{SNM}}$ (MeV/fm3) |  | Ref.
HIC(DLL) | 0.32 | $10.1\pm 3.0$ |  | Danielewicz  _et al._ (2002)
HIC(FOPI) | 0.32 | $10.3\pm 2.8$ |  | Le Fèvre  _et al._ (2016)
Asymmetric matter |  |  |  |
Constraints |  $n$ (fm-3) | S($n$) (MeV) | P${}_{\text{\rm sym}}$ (MeV/fm3) | Ref.
Nuclear structure |  |  |  |
$\alpha_{D}$ | 0.05 | $15.9\pm 1.0$ |  |  Zhang and Chen (2015)
PREX-II | 0.11 |  | $2.38\pm 0.75$ |  Adhikari  _et al._ (2021); Reed  _et al._ (2021); Lynch and Tsang (2022)
Nuclear masses |  |  |  |
Mass(Skyrme) | 0.101 | $24.7\pm 0.8$ |  |  Brown and Schwenk (2014); Lynch and Tsang (2022)
Mass(DFT) | 0.115 | $25.4\pm 1.1$ |  |  Kortelainen  _et al._ (2012); Lynch and Tsang (2022)
IAS | 0.106 | $25.5\pm 1.1$ |  |  Danielewicz  _et al._ (2017); Lynch and Tsang (2022)
Heavy-ion collisions |  |  |  |
HIC(Isodiff) | 0.035 | $10.3\pm 1.0$ |  |  Tsang  _et al._ (2009); Lynch and Tsang (2022)
HIC(n/p ratio) | 0.069 | $16.8\pm 1.2$ |  |  Morfouace  _et al._ (2019); Lynch and Tsang (2022)
HIC($\pi$) | 0.232 | $52\pm 13$ | $10.9\pm 8.7$ |  Estee  _et al._ (2021); Lynch and Tsang (2022)
HIC(n/p flow) | 0.240 |  | $12.1\pm 8.4$ |  Cozma (2018); Russotto  _et al._ (2011); Russotto and et. al. (2016); Lynch and Tsang (2022)

Astronomical |  |  |  |
---|---|---|---|---
Constraints | $M_{\odot}$ | R (km) | $\Lambda_{1.36}$ | Ref.
LIGO | 1.36 |  | $300_{-230}^{+420}$ |  Abbott et al. (2019a)
*Riley PSR J0030+0451 | 1.34 | $12.71^{+1.14}_{-1.19}$ |  |  Riley et al. (2019)
*Miller PSR J0030+0451 | 1.44 | $13.02^{+1.24}_{-1.06}$ |  |  Miller et al. (2019)
*Riley PSR J0740+6620 | 2.07 | 12.39 |  |  Riley et al. (2021)
*Miller PSR J0740+6620 | 2.08 | $13.7~{}^{+2.6}_{-1.5}$ |  |  Miller  _et al._ (2021)

###  II.2 RMF EoS

We also create a collection of EoSs within a model of nuclear matter that is
based on a relativistic field theory approach. The nuclear interaction between
nucleons is represented by the exchange of the scalar-isoscalar meson
$\sigma$, the vector-isoscalar meson $\omega$, and the vector-isovector meson
$\varrho$.

The Lagrangian density is given by (Fattoyev  _et al._ , 2010; Dutra  _et al._
, 2014; Malik  _et al._ , 2023)

$\mathcal{L}=\mathcal{L}_{N}+\mathcal{L}_{M}+\mathcal{L}_{NL}+\mathcal{L}_{leptons}$
(10)

where

$\mathcal{L}_{N}=\bar{\Psi}\Big{[}\gamma^{\mu}\left(i\partial_{\mu}-g_{\omega}\omega_{\mu}-\frac{1}{2}g_{\varrho}{\bm{t}}\cdot\bm{\varrho}_{\mu}\right)-\left(m_{N}-g_{\sigma}\sigma\right)\Big{]}\Psi$

denotes the Dirac equation for the nucleon doublet (neutron and proton) with
bare mass $m_{N}$, $\Psi$ is a Dirac spinor, $\gamma^{\mu}$ are the Dirac
matrices, and $\bm{t}$ is the isospin operator. $\mathcal{L}_{M}$ is the
Lagrangian density for the mesons, given by

$\displaystyle\mathcal{L}_{M}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\partial_{\mu}\sigma\partial^{\mu}\sigma-
m_{\sigma}^{2}\sigma^{2}\right]-\frac{1}{4}F_{\mu\nu}^{(\omega)}F^{(\omega)\mu\nu}+\frac{1}{2}m_{\omega}^{2}\omega_{\mu}\omega^{\mu}$
$\displaystyle-$
$\displaystyle\frac{1}{4}\bm{F}_{\mu\nu}^{(\varrho)}\cdot\bm{F}^{(\varrho)\mu\nu}+\frac{1}{2}m_{\varrho}^{2}\bm{\varrho}_{\mu}\cdot\bm{\varrho}^{\mu}$

where
$F^{(\omega,\varrho)\mu\nu}=\partial^{\mu}A^{(\omega,\varrho)\nu}-\partial^{\nu}A^{(\omega,\varrho)\mu}$
are the vector meson tensors, and

$\displaystyle\mathcal{L}_{NL}$ $\displaystyle=$
$\displaystyle-\frac{1}{3}b~{}m_{N}~{}g_{\sigma}^{3}(\sigma)^{3}-\frac{1}{4}c(g_{\sigma}\sigma)^{4}+\frac{\xi}{4!}g_{\omega}^{4}(\omega_{\mu}\omega^{\mu})^{2}$
$\displaystyle+$
$\displaystyle\Lambda_{\omega}g_{\varrho}^{2}\bm{\varrho}_{\mu}\cdot\bm{\varrho}^{\mu}g_{\omega}^{2}\omega_{\mu}\omega^{\mu}$

contains the non-linear mesonic terms with parameters $b$, $c$, $\xi$,
$\Lambda_{\omega}$ to take care of the high-density behavior of nuclear
matter. The parameters $g_{i}$’s are the couplings of the nucleons to the
meson fields $i=\sigma,\omega,\varrho$, with masses $m_{i}$. Finally, the
Lagrangian density for the leptons is given as
$\mathcal{L}_{leptons}=\bar{\Psi_{l}}\Big{[}\gamma^{\mu}\left(i\partial_{\mu}-m_{l}\right)\Psi_{l}\Big{]}$,
where $\Psi_{l}~{}(l=e^{-},\mu^{-})$ denotes the lepton spinor for electrons
and muons; leptons are considered non-interacting. The equations of motion for
the meson fields are obtained from the Euler-Lagrange equations:

$\displaystyle{\sigma}$ $\displaystyle=$
$\displaystyle\frac{g_{\sigma}}{m_{\sigma,{\rm eff}}^{2}}\sum_{i}\rho^{s}_{i}$
(11) $\displaystyle{\omega}$ $\displaystyle=$
$\displaystyle\frac{g_{\omega}}{m_{\omega,{\rm eff}}^{2}}\sum_{i}\rho_{i}$
(12) $\displaystyle{\varrho}$ $\displaystyle=$
$\displaystyle\frac{g_{\varrho}}{m_{\varrho,{\rm
eff}}^{2}}\sum_{i}t_{3}\rho_{i},$ (13)

where $\rho^{s}_{i}$ and $\rho_{i}$ are, respectively, the scalar density and
the number density of nucleon $i$, and the effective meson masses $m_{i,{\rm
eff}}$ are defined as

$\displaystyle m_{\sigma,{\rm eff}}^{2}$ $\displaystyle=$ $\displaystyle
m_{\sigma}^{2}+{b~{}m_{N}~{}g_{\sigma}^{3}}{\sigma}+{cg_{\sigma}^{4}}{\sigma}^{2}$
(14) $\displaystyle m_{\omega,{\rm eff}}^{2}$ $\displaystyle=$ $\displaystyle
m_{\omega}^{2}+\frac{\xi}{3!}g_{\omega}^{4}{\omega}^{2}+2\Lambda_{\omega}g_{\varrho}^{2}g_{\omega}^{2}{\varrho}^{2}$
(15) $\displaystyle m_{\varrho,{\rm eff}}^{2}$ $\displaystyle=$ $\displaystyle
m_{\varrho}^{2}+2\Lambda_{\omega}g_{\omega}^{2}g_{\varrho}^{2}{\omega}^{2},.$
(16)

The meson equations are solved using the relativistic mean-field approximation
for a given $\rho_{n}$ and $\rho_{p}$. Note, that the effect of the nonlinear
terms on the magnitude of the meson fields enters through the effective meson
masses, $m_{i,{\rm eff}}$. The energy density of the baryons and leptons are
given by the following expressions:

$\displaystyle\epsilon$
$\displaystyle=\sum_{i=n,p,e,\mu}\frac{1}{\pi^{2}}\int_{0}^{k_{Fi}}\sqrt{k^{2}+{m_{i}^{*}}^{2}}\,k^{2}\,dk$
(17)
$\displaystyle+\frac{1}{2}m_{\sigma}^{2}{\sigma}^{2}+\frac{1}{2}m_{\omega}^{2}{\omega}^{2}+\frac{1}{2}m_{\varrho}^{2}{\varrho}^{2}$
$\displaystyle+\frac{b}{3}(g_{\sigma}{\sigma})^{3}+\frac{c}{4}(g_{\sigma}{\sigma})^{4}+\frac{\xi}{8}(g_{\omega}{\omega})^{4}+\Lambda_{\omega}(g_{\varrho}g_{\omega}{\varrho}{\omega})^{2},$

where $m_{i}^{*}=m_{i}-g_{\sigma}\sigma$ for protons and neutrons and
$m_{i}^{*}=m_{i}$ for electrons and muons, and $k_{Fi}$ is the Fermi moment of
particle $i$.

Once we have the energy density for a given EoS model, we can compute the
chemical potential of neutron ($\mu_{n}$) and proton ($\mu_{p}$). The chemical
potential of electron ($\mu_{e}$) and muon ($\mu_{\mu}$) can be computed using
the condition of $\beta$-equilibrium : $\mu_{n}-\mu_{p}=\mu_{e}$ and $\mu_{e}$
= $\mu_{\mu}$ and the charge neutrality: $\rho_{p}=\rho_{e}+\rho_{\mu}$. Where
$\rho_{e}$ and $\rho_{\mu}$ are the electron and muon number density. The
pressure is then determined from the thermodynamic relation:

$p=\sum_{i}\mu_{i}\rho_{i}-\epsilon.$ (18)

###  II.3 Bayesian Likelihood

The Bayesian likelihood is a fundamental concept in Bayesian statistics,
providing a way to update the probability estimate for a hypothesis as more
evidence or information becomes available. In the context of Bayesian
inference, the posterior distribution measures the plausibility of a set of
parameters values given the observed data.

Bayes’ Theorem:- Bayes’ theorem is used to update the probability estimate of
a hypothesis as more evidence is available. It is formulated as:

$P(\theta|{\mathcal{D}})=\frac{{\mathcal{L}}(\mathcal{D}|{\theta})P(\theta)}{P({\mathcal{D}})}$
(19)

In this equation:

  * •

$P(\theta|{\mathcal{D}})$ is the posterior probability of the parameters
$\theta$ given the data, $\mathcal{D}$.

  * •

${\mathcal{L}}(\mathcal{D}|{\theta})$ is the likelihood of the data under the
parameters (likelihood function).

  * •

$P(\theta)$ is the prior probability of the parameters $\theta$.

  * •

$P({\mathcal{D}})$ is the marginal likelihood or evidence, often acting as a
normalizing constant.

Likelihood:- The likelihood function in Bayesian analysis is defined as the
probability of the observed data under a specific statistical model,
parameterized by a set of parameters $\theta$. To include the asymmetric
uncertainties (${\mathcal{D}}_{\sigma_{-}}^{\sigma^{+}}$) in the NS properties
as listed in Tab1 the likelihood function considered is

$\displaystyle ln({\mathcal{L}}({\mathcal{D}|}\theta))$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{i}\left(w_{i}(\frac{R_{i}}{\sigma_{i}})^{2}+ln(2\pi\sigma_{i}^{2})\right)$
(20)

Here $R_{i}$ is the residual between the the data, $\mathcal{D}_{i}$ and model
value $M_{i}(\theta)$ which is defined as,

$\displaystyle R_{i}$ $\displaystyle=$
$\displaystyle{\mathcal{D}_{i}}-M_{i}(\theta)$ (21)

$\sigma_{i}$ is the standard deviation either $\sigma_{-}$ or $\sigma_{+}$
depending on whether $\frac{R_{i}}{\mid R_{i}\mid}$ is +1 or -1. $w_{i}$ is
the weighting factor which is always 1 except $\frac{1}{2}$ in the likelihood
evaluation of NS radii (Since there are two sets of measurements for the same
objects).

##  III Terrestrial and astrophysics constraints on EoSs

We use relativistic and non-relativistic mean-field models, as briefly
outlined in the previous section, to calculate the properties of nuclear
matter at several densities and the EoSs for the neutron star matter. The
properties of nuclear matter considered are the pressure for the symmetric
nuclear matter ($P_{\rm SNM}$), symmetry energy pressure ($P_{\rm sym}$) and
the symmetry energy ($e_{\rm sym}$) which are constrained empirically over a
range of densities by the experimental data on bulk properties of finite
nuclei such as the nuclear masses (2), neutron skin thickness in 208Pb (1) and
the dipole polarizability (1), isobaric analog states (1) as well as the HIC
data (7) spanning the density range $0.03$ \- $0.32$ fm-3, where the number of
data points for a given quantity has been indicated in the parentheses. The
astrophysics data considered are the radii of neutron stars with masses around
the canonical one and $\sim$ 2.1 M⊙, the dimensionless tidal deformability for
the neutron star around the canonical mass and the current bound on the
maximum mass of the neutron star. The EoSs are subjected to these constraints
within the Bayesian framework. The terrestrial and astrophysics data Tsang
_et al._ (2023) considered are listed in Tab 1. It should noted that the tidal
deformability constraints considered in this work is specifically for NS with
a mass of 1.36 M⊙. This represents the limitations on the combined tidal
deformability, denoted as $\tilde{\Lambda}$, for the GW170817 binary
companion. When the mass ratio, q, is equal to 1, the combined tidal
deformability is equal to $\Lambda_{1.36}$. It is worth noting that the
$\tilde{\Lambda}$ value is almost constant, almost insensitive to the NS mass
ratio q Raithel  _et al._ (2018); Malik and Pais (2022).

Table 2:  The median and 68$\%$ confidence interval for the nuclear matter parameters obtained for three different cases as labelled by $M_{\rm max}$, $M_{\rm max}$ \+ EXPT and ALL described in the text. Each case is studied with 2.10 and 2.17 M⊙ as the lower limit on the neutron star maximum mass. NMPs | Models | Mmax |  $M_{\rm max}$ \+ EXPT | ALL
---|---|---|---|---
(in MeV) |  $\geq$2.10 |  $\geq$2.17 |  $\geq$2.10 |  $\geq$2.17 |  $\geq$2.10 |  $\geq$2.17
K0 | Skyrme | 246.73${}_{-20.89}^{+22.03}$ | 257.42${}_{-14.89}^{+18.13}$ | 242.46${}_{-10.61}^{+9.88}$ | 257.69${}_{-10.79}^{+12.91}$ | 244.10${}_{-12.38}^{+12.47}$ | 264.47${}_{-10.18}^{+5.51}$
RMF | $231.77_{-25.47}^{+26.76}$ | 231.70${}_{-27.76}^{+25.88}$ | 202.94${}_{-16.16}^{+18.26}$ | 199.71${}_{-16.82}^{+17.71}$ | 195.28${}_{-16.50}^{+18.38}$ | 206.92${}_{-13.72}^{+16.38}$
J0 | Skyrme | 33.24${}_{-5.29}^{+4.92}$ | 33.96${}_{-5.32}^{+4.16}$ | 35.52${}_{-1.75}^{+1.43}$ | 35.32${}_{-1.63}^{+1.70}$ | 34.42${}_{-1.45}^{+1.45}$ | 35.22${}_{-1.67}^{+1.07}$
RMF | $29.80_{-3.84}^{+11.06}$ | 29.47${}_{-3.8}^{+10.97}$ | 38.06${}_{-0.98}^{+1.08}$ | 38.25${}_{-1.06}^{+1.05}$ | 38.16${}_{-1.06}^{+1.07}$ | 38.39${}_{-1.02}^{+0.99}$
L0 | Skyrme | 77.18${}_{-35.29}^{+42.93}$ | 78.09${}_{-32.09}^{+38.64}$ | 94.07${}_{-13.91}^{+13.45}$ | 93.04${}_{-14.09}^{+13.98}$ | 81.14${}_{-11.10}^{+10.87}$ | 82.50${}_{-10.38}^{+8.62}$
RMF | $46.80_{-12.27}^{+23.63}$ | 47.57${}_{-11.37}^{+22.08}$ | 107.50${}_{-5.59}^{+4.74}$ | 108.92${}_{-5.72}^{+4.53}$ | 107.28${}_{-6.19}^{+5.17}$ | 107.92${}_{-5.82}^{+4.54}$
Ksym0 | Skyrme | 72.78${}_{-112.04}^{+107.35}$ | 92.94${}_{-120.33}^{+93.84}$ | 76.08${}_{-73.94}^{+82.60}$ | 83.00${}_{-78.04}^{+78.10}$ | 6.50${}_{-55.93}^{+63.73}$ | 10.74${}_{-43.81}^{+42.05}$
RMF | $-71.00_{-92.54}^{+55.46}$ | -57.33${}_{-90.68}^{+58.98}$ | 2.75${}_{-32.12}^{+20.36}$ | 10.30${}_{-33.48}^{+20.78}$ | 0${}_{-37.61}^{+25.45}$ | -0.40${}_{-38.69}^{+22.35}$
Table 3:  The median and 68$\%$ confidence interval for the radius of neutron star of mass 1.4, 2.08 M⊙ and tidal deformability for 1.36 M⊙ NS obtained for three different cases as labelled by $M_{\rm max}$, $M_{\rm max}$ \+ EXPT and ALL described in the text. NS properties | Models | Mmax |  $M_{\rm max}$ \+ EXPT | ALL
---|---|---|---|---
$\geq$2.10 |  $\geq$2.17 |  $\geq$2.10 |  $\geq$2.17 |  $\geq$2.10 |  $\geq$2.17
R1.4 (km) | Skyrme | 13.37${}_{-1.69}^{+1.10}$ | 13.42${}_{-0.92}^{+1.51}$ | 13.79${}_{-0.43}^{+0.47}$ | 13.77${}_{-0.46}^{+0.43}$ | 13.28${}_{-0.40}^{+0.39}$ | 13.18${}_{-0.26}^{+0.32}$
RMF | $12.35_{-0.34}^{+0.54}$ | 12.44${}_{-0.23}^{+0.54}$ | 13.10${}_{-0.23}^{+0.19}$ | 13.17${}_{-0.21}^{+0.18}$ | 13.05${}_{-0.19}^{+0.22}$ | 13.05${}_{-0.20}^{+0.21}$
R2.08(km) | Skyrme | 12.30${}_{-0.68}^{+0.70}$ | 12.65${}_{-0.59}^{+0.63}$ | 12.06${}_{-0.36}^{+0.34}$ | 12.47${}_{-0.30}^{+0.28}$ | 11.75${}_{-0.33}^{+0.36}$ | 11.99${}_{-0.18}^{+0.29}$
RMF | $11.69_{-0.40}^{+0.47}$ | 12.02${}_{-0.23}^{+0.54}$ | 11.86${}_{-0.31}^{+0.36}$ | 12.20${}_{-0.25}^{+0.27}$ | 12.10${}_{-0.26}^{+0.23}$ | 11.97${}_{-0.19}^{+0.18}$
$\Lambda_{1.36}$ (…) | Skyrme | 824${}_{-255}^{+317}$ | 885${}_{-246}^{+286}$ | 805${}_{-129}^{+133}$ | 838${}_{-133}^{+130}$ | 660${}_{-97}^{+106}$ | 664${}_{-64}^{+85}$
RMF | 508${}^{+104}_{-81}$ | 548${}_{-87}^{+102}$ | $622_{-74}^{+71}$ | 645${}_{-72}^{+67}$ | 607${}_{-66}^{+70}$ | 610${}_{-69}^{+70}$

##  IV Results and Discussion

We have applied two distinct models outlined in Section II to construct the
EoSs. These EoS models used in this analysis are (i) Skyrme based non-
relativistic mean field and (ii) RMF models. These EoSs are employed to
compute esym, Psym, and PSNM corresponding to the experimental data derived
from finite nuclei and HIC experiments, as detailed in the preceding section.

The radius (R), the dimensionless tidal deformability ($\Lambda$) and the
maximum mass of NS ($M_{\rm max}$) for a given EoS are obtained from the
solutions of the Tolman–Oppenheimer–Volkoff (TOV) equations Oppenheimer and
Volkoff (1939); Tolman (1939); Hinderer (2008). The EoS for the outer crust up
to 0.0016$\rho_{0}$ is taken to be the one given by Baym-Pethick-Sutherland
Baym  _et al._ (1971). The EoS for the inner crust corresponding to the
density range $0.0016\rho_{0}<\rho<0.5\rho_{0}$ is obtained by using a
polytropic form Carriere  _et al._ (2003); Patra  _et al._ (2022). The core
EoSs are constructed using non-relativistic and relativistic mean field models
as described in SecII. These EoSs are subjected to (i) thermodynamic
stability, (ii) causality (speed of sound below 1), (iii) positive semi-
definiteness of the symmetry energy, and (iv) a maximum mass of the stable
non-rotating neutron stars exceeding 2.1 M⊙Romani  _et al._ (2022) as
discussed in the Introduction. These EoSs are further constrained within the
Bayesian framework by experimental data and astrophysical observations. The
prior distributions of the parameters of these models are provided in Tabs.A1
and A2 of the appendix.

Figure 1:  Distributions of the value of the log likelihood for the case of
all the data considered within the Skyrme (top) and RMF (bottom) models. The
contributions to the logarithmic likelihood of each distinct group of data are
also displayed for comparison (see text for details). Figure 2:  The 68%
confidence interval for $e_{\rm sym}$, $P_{\rm sym}$, $P_{\rm SNM}$ at fixed
densities (left) and radius and tidal deformability of NS (right) for the case
“ALL”. The fixed values of densities and NS masses are specified along the
abscissa as listed in Tab1. The results are obtained using the posterior
distributions of nuclear matter parameters as presented in Fig. (3). The
dashed vertical lines with hollow circles representing the
experimental/observational data and the corresponding results from Skyrme(top)
and RMF(bottom) models are depicted by the solid vertical lines with solid
circles. The scale along the ordinate only represents the numerical values,
while the associated units for the presented quantities are provided in the
top panel. Figure 3:  Corner plots obtained for the Skyrme and RMF models with
the 12 experimental data and 3 astrophysical observations within the Bayesian
framework. The marginalized posterior distributions of the nuclear matter
parameters and selected neutron star properties are plotted along the
diagonal. The confidence ellipses for two-dimensional posterior distributions
are plotted with $1\sigma$, $2\sigma$ and $3\sigma$ intervals along the off-
diagonal.

In Fig.(1), we display the distributions of the logarithm of the likelihood
obtained for both the EoS models considered. The curves labeled as “ALL”
represent the results derived from employing all the data considered.
Additionally, to allow a comparison, we depict contributions from various data
groups. Specifically, the data groups denoted as PREX, HIC, and Astro
correspond to single data from PREX-II, seven data from HIC, and five data
from astrophysical observations, respectively. The group labeled as “Rest”
encompasses the remaining four data points from nuclear masses, dipole
polarizability, and IAS. Looking at the figure, it is evident that the peak
positions of the dashed curves for both the models are very close to each
other. The major contribution to the likelihood stems from the HIC data (in
red) since it is closest to the total likelihood (in dashed green) for both
the models. The constraints imposed on the symmetry energy from finite nuclei
data (in blue) appear to align somewhat better with the Skyrme model.
Meanwhile, the astrophysical observations related to neutron stars seem to be
fitted better marginally with the RMF model (in black). The constraint on Psym
from PREX-II seems to be fitted equally well (in brown) for both the models.
In Fig. (2), we display the values of $e_{\rm sym}$, $P_{\rm sym}$, and
$P_{\rm SNM}$ at specific fixed densities, along with R1.4, R2.08, and
$\Lambda_{1.36}$ within a 68$\%$ confidence interval. These values are derived
from the Skyrme (top) and RMF (bottom) models within the Bayesian framework.
Additionally, for comparison, we include corresponding values obtained from
experimental data and astrophysical observations as used in the present work
to constrain the EoS (see Tab1). The value of PSNM (in magenta) aligns well
with the RMF model, whereas, esym (in green) from the nuclear masses are
fitted well with the Skyrme model. The astrophysical observations related to
NS properties show similar agreement for both the models.

To analyse systematically the implications of various data in constraining the
EoS, we perform our calculations in the following steps, (i) Priors are
updated with the current lower bound for the maximum mass of neutron stars
($M_{\rm max}\geq$2.1 M⊙). (ii) Next, we further include 12 experimental data
points. (iii) All data considered in Tab1 are incorporated. The outcomes for
the distributions of nuclear matter parameters and neutron star properties
under these three scenarios, identified as “$M_{\rm max}$”, “$M_{\rm max}$ \+
EXPT”, and “ALL”, are detailed in Tabs.2 and 3, respectively. From Tab2 it is
clear that for the Skyrme model the K0 is essentially constrained by the lower
limit of NS maximum mass, while the experimental data and the other
astrophysical observations have a smaller effect. In the case of the RMF model
the maximum mass as well as experimental data have a great impact on K0. The
other parameters are also significantly constrained by the experimental data
in the RMF case. However, for the Skyrme model the Ksym0 is mainly constrained
by the astrophysical data. Finally from the distributions of the nuclear
matter parameters in the case of ALL we can conclude that the Skyrme presents
a stiffer SNM EoS but a softer symmetry energy compared to RMF in the vicinity
of $\rho_{0}$.

It can be seen from the Tab3 that the “$M_{\rm max}$” and “$M_{\rm max}$ \+
EXPT” scenarios, for the two types of models used — Skyrme and RMF — show
noticeable differences in predicting the radius for neutron stars with 1.4 M⊙
and the tidal deformability ($\Lambda$) with 1.36 M⊙ measured in GW170817.
However, for higher NS masses, like 2.08 M⊙, the predicted radii $R_{2.08}$
are very similar. When astrophysical observations on NS radii and tidal
deformability are included, both the models predict similar radii for neutron
stars of both 1.4 and 2.08 M⊙, as well as a similar tidal deformability. It is
worth noting , though, agreement of Skyrme and RMF models with the
experimental data on PSNM and esym are different, the EoS for the NS matter
are such that the differences in the NS properties level off. It may be
pointed out that when all the data/observations are considered the shift on
the lower bound on Mmax from 2.1 $M_{\odot}$ to 2.17 $M_{\odot}$ has only a
small effect on the NMP giving rise a slightly larger $K_{0}$ (by 20 MeV for
the Skyrme models and 10 MeV for the RMF models).All the other properties do
not differ much more than 1 MeV, the $K_{\rm sym0}$ with the Skyrme model
being the only exception with a difference of $\sim 3$ MeV. So, the increase
in the lower limit of the maximum mass M${}_{\rm max}\geq$ to 2.17 M⊙ has a
negligible impact on the outcomes related to the NS properties. We have also
repeated our all the calculations with variable transition density between
0.04 -0.08 fm-3 depending on the values of L0 Malik  _et al._ (2024) to asses
the robustness of our results. The values of radius and tidal deformability
change minimally only by less than 1$\%$ and 5$\%$, respectively.

The corner plots in Fig.(3) display nuclear matter parameters and selected NS
properties under the ‘ALL’ scenario for both Skyrme and RMF models. Along the
diagonal, marginalised distributions are presented, while off-diagonal
sections depict confidence ellipses. The distributions of the nuclear matter
parameters for the Skyrme model are relatively wider compared to those for the
RMF models. It may be emphasized that the values of nuclear matter parameters
obtained from the Skyrme and RMF models differ significantly, despite being
inferred from the same data set. For instance, the median value for K0 is
higher in the Skyrme model by about 20$\%$ compared to the RMF model, while
the median value of L0 in the Skyrme model is lower by the same amount. This
suggests that the Skyrme model yields a stiffer EoS for symmetric nuclear
matter around $\sim$ $\rho_{0}$, but with a softer density-dependent symmetry
energy. Further L0 exhibits correlations with J0 and Ksym0 for both the
models. However, the strength of correlations between the remaining pairs of
nuclear matter parameters are relatively weaker. In the case of the Skyrme
model, the correlations between L0 and Ksym0 with R1.4 and $\Lambda_{1.36}$
appear relatively stronger compared to those in the RMF modelFortin  _et al._
(2016); Tsang  _et al._ (2020); Thi  _et al._ (2021); Patra  _et al._ (2022);
Patra et al. (2023); Imam  _et al._ (2023); Patra  _et al._ (2023). It may be
pointed out that R1.4 is better constrained for the RMF model. This may be
associated to the narrow distribution of L0.

Figure 4:  The $95\%$ confidence intervals for the radius and tidal
deformability as a function of neutron star mass evaluated using the posterior
distributions. The astrophysical observations used in the Bayesian framework
are also shown. The error bars in black and red colors indicate the
constraints from two different analysis for NS radii for each of the pulsar
PSR J0030+0451 and PSR J0740+6620 with masses $\sim$ 1.4 and 2.0 M⊙,
respectively.

We examine the impact of three different scenarios on the distributions of
radius, tidal deformability of NS and the underlying EoSs. In Fig.(4), the
$95\%$ confidence intervals display the NS mass-radius and mass-tidal
deformability relationships. Results from both the Skyrme and the RMF model
exhibit similar trends. The incorporation of experimental data (magenta band)
significantly reduces the spread in NS properties at a given mass compared to
the scenario involving only the maximum mass constraint (empty blue band). The
values of R1.4, R2.08, and $\Lambda_{1.36}$ are already reasonably constrained
by $M_{\rm max}$ and experimental data, aligning well with current
astrophysical observations. Interestingly, the inclusion of constraints from
astrophysical observations (green band) does not significantly alter the NS
properties. This suggests a need for more precise observations to further
refine and constrain the underlying EoS.

In Fig.(5), the $95\%$ confidence intervals represent the updated EoSs
achieved by sequentially incorporating the constraints. The original priors
(grey band) are also depicted for comparison purposes. The maximum NS mass
significantly constrains the EoSs (empty blue band) which further narrow down
with the inclusion of experimental data (magenta band). The introduction of
constraints from astrophysical observations does not notably impact the EoSs.

The Kullback-Liebler divergence (DKL) provides a quantitative comparison
between two different distributions. The DKL comparing distributions P and Q
is defined as :
DKL(P$\mid\mid$Q) = $\int_{-\infty}^{+\infty}$
P(x)log2($\frac{P(x)}{Q(x)}$)dx. For a more detailed comparison, we evaluate
the DKL as a function of baryon density for two distributions of the EoSs
corresponding to different scenarios. In Fig.(6), the solid lines represent
the value of DKL obtained by comparing the distributions of EoSs between
“$M_{\rm max}$” and “$M_{\rm max}$ \+ EXPT” scenarios. Similarly, the dashed
lines represent the value of DKL obtained by comparing the distributions of
EoSs between “$M_{\rm max}$” and “ALL” scenarios. The values of DKL for the
solid lines indicate that the experimental data provides tighter constraints
on the EoS in the vicinity of $\rho_{0}$. The proximity between solid and
dashed lines indicates that including astrophysical observations does not
impose stricter constraints beyond those set by the experimental data.
However, the presence of a minor peak around 2-3$\rho_{0}$ in both models
suggests that astrophysical observations on NS radii and tidal deformability
introduce additional constraints within this density range.

Figure 5:  The $95\%$ confidence intervals for the posterior distributions of
the pressure for the neutron star matter as a function of baryon density for
the RMF (top) and Skyrme (bottom) mean field models corresponding to different
scenarios as indicated (see text).

Finally, we would like to compare our results on the radius of the
1.4$M_{\odot}$ star with other results existing in the literature. In Ref.Pang
_et al._ (2023), the authors made an estimation of the radius of a
1.4$M_{\odot}$ star, $R_{1.4}=11.98^{+0.35}_{-0.40}$ km, using an extension of
the Nuclear-physics and Multi-Messenger Astrophysics framework NMMA, which
includes data from gravitational-wave signal from the GW170817 event, the
kilonova AT2017gfo and the gamma-ray burst GRB170817A afterglow, together with
nuclear-physics constraints at low densities and X-ray and radio observations
of isolated neutron stars. This predicted value is in agreement with several
other predictions listed in Tab1 of Ref.Pang  _et al._ (2023) which have been
obtained essentially from observational data and nuclear-physics constraints
at low densities for $\chi$EFT calculations for neutron matter, except for the
study Huth  _et al._ (2022) where data from HIC Russotto and et. al. (2016);
Le Fèvre  _et al._ (2016); Danielewicz  _et al._ (2002) has also been
considered. In Ref.Huth  _et al._ (2022) the predicted radius of a
1.4$M_{\odot}$ NS is 12.01${}^{+0.78}_{-0.77}$ km considering Astro and HIC
constraints, however taking only the HIC constraints the radius is larger,
12.06${}^{+1.13}_{-1.18}$ km. The quite large radius predicted in our study
13.28${}_{-0.40}^{+0.39}$ km (Skyrme) and 13.05${}_{-0.19}^{+0.22}$ km (RMF)
is compatible with the prediction from Huth  _et al._ (2022) from the HIC
data. Note, however, that in our study we have considered a larger set of HIC
and experimental nuclear data. The HIC and nuclear data seem to give rise to
larger radii than the astrophysics observations. These are clearly seen from
Table 3: nuclear data increase the NS radius with respect just to the maximum
mass constraint, and the rest of the astrophysical data pushes down the radius
coming from the nuclear data plus maximum mass constraints. In Ref.Tsang  _et
al._ (2023) the R1.4 = 12.9${}_{-0.5}^{+0.4}$ km has been obtained using a
diverse set of data from finite nuclei, HIC, GW170817 and NICER which compare
well with R1.4 $\sim$ 13 km as obtained in the present work (see Tab3). These
results suggest that lower values of R1.4 obtained in Pang  _et al._ (2023)
and references therein may be due to the constraints from $\chi$EFT for pure
neutron matter at low density. In a recent study Breschi  _et al._ (2024) a
combined analysis has been performed using GW170817, kilonova and NICER data
which predict radius of 1.4 M⊙ NS 12.30${}_{-0.56}^{+0.81}$ km and
13.20${}_{-0.91}^{+0.90}$ km for two different hot spot models used in the
NICER analysis.

Figure 6:  The variations of Kullback–Leibler divergence (DKL) as a function
of baryon density. The values of DKL shown by solid (dashed) line is obtained
for the scenario Mmax \+ EXPT (ALL) with respect to the prior updated with
maximum neutron star mass.

##  V Conclusions

We have performed a systematic investigation to assess the implications of a
comprehensive set of data, encompassing nuclear physics experiments and
astrophysical observations concerning neutron star properties, on the EoSs
obtained within the non-relativistic and relativistic mean field models. The
non-relativistic model relies on a Skyrme-type effective interaction, while
the relativistic mean field model is derived from an effective Lagrangian
which incorporates nonlinear terms for $\sigma$, $\omega$, and $\rho$ mesons.
The experimental data set includes bulk properties from finite nuclei and
heavy-ion collisions, contributing with information on the symmetry energy,
the symmetry energy pressure, and the pressure for symmetric nuclear matter up
to twice the saturation density. Additionally, astrophysical observations
including radii for approximately 1.4 M⊙ and 2.1 M⊙, tidal deformability for
1.36 M⊙, and the maximum neutron star mass are considered. We have employed a
Bayesian framework, to investigate the constraints imposed by these data on
nuclear matter properties. Three scenarios were explored: (i) priors updated
with the maximum neutron star mass constraint only, (ii) maximum neutron star
mass plus experimental data, and (iii) all the data considered.

Our determined values largely align within 68$\%$ confidence intervals with
most of the fitted data/observations (see Fig.2). The non-relativistic model
underestimates the pressure for symmetric nuclear matter, while the RMF model
tends to overestimate the symmetry energy at densities of 0.05 and 0.115 fm-3
derived from iso-vector dipole resonance and nuclear masses, respectively. The
dimensionless tidal deformability from the RMF model aligns slightly better
with the corresponding astrophysical observations for a neutron star of 1.36
M⊙. The underlying EoSs initially constrained by the maximum neutron star
mass, narrow further with the inclusion of experimental data. Interestingly,
there is a minimal impact on the EoSs by further inclusion of astrophysical
observations , especially evident in the distributions of neutron star mass-
radius and mass-tidal deformability. These trends are as-lo quantitatively
reflected from the values of KL divergence obtained by comparing the
distributions of the EoS from different scenarios. It implies more precise
observations require to confine the EoSs in narrower bounds. Finally, we
remark that there is a significant discrepancy in the values of nuclear matter
parameters obtained from the Skyrme and RMF models, despite both originating
from the same data set. In particular, as alluded above, Skyrme forces predict
a harder symmetric nuclear matter EoS while the RMF predicts a harder symmetry
energy. It is also observed that the value of R${}_{1.4}\sim$ 12 km in the
previous analysis may be due to the low density constraints from $\chi$EFT for
pure neutron matter.

##  VI Acknowledgements

TM and CP would like to appreciate the support of national funds from FCT
(Fundação para a Ciência e a Tecnologia, I.P, Portugal) under projects
UIDP/04564/2020 and No. UIDB/04564/2020, respectively with DOI identifiers
10.54499/UIDB/04564/2020 and 10.54499/UIDP/04564/2020, and Project
2022.06460.PTDC with the DOI identifier 10.54499/2022.06460.PTDC. They also
like to acknowledge the Laboratory for Advanced Computing at the University of
Coimbra for providing HPC resources that have contributed to the research
results reported within this paper, URL: https://www.uc.pt/lca. BKA
acknowledges partial support from the SERB, Department of Science and
Technology, Government of India with CRG/2021/000101. The authors sincerely
acknowledge the usage of the analysis software BILBY Ashton et al. (2019);
Romero-Shaw et al. (2020).

## References

  * Lattimer and Prakash (2001) J. M. Lattimer and M. Prakash, Astrophys. J. 550, 426 (2001).
  * Lattimer (2021) J. Lattimer, Annual Review of Nuclear and Particle Science 71, 433 (2021).
  * Fattoyev  _et al._ (2010) F. J. Fattoyev, C. J. Horowitz, J. Piekarewicz, and G. Shen, Phys. Rev. C 82, 055803 (2010).
  * Tsang  _et al._ (2012) M. B. Tsang _et al._ , Phys. Rev. C 86, 015803 (2012).
  * Fortin  _et al._ (2016) M. Fortin, C. Providência, A. R. Raduta, F. Gulminelli, J. L. Zdunik, P. Haensel, and M. Bejger, Phys. Rev. C 94, 035804 (2016).
  * Margueron  _et al._ (2018) J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025806 (2018).
  * Fattoyev  _et al._ (2018) F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, Phys. Rev. Lett. 120, 172702 (2018).
  * Greif  _et al._ (2020) S. K. Greif, K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 901, 155 (2020).
  * Reed  _et al._ (2023) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, (2023), arXiv:2305.19376 [nucl-th] .
  * Raaijmakers  _et al._ (2021) G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, Astrophys. J. Lett. 918, L29 (2021).
  * Annala  _et al._ (2022) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Phys. Rev. X 12, 011058 (2022).
  * Ghosh  _et al._ (2022) S. Ghosh, B. K. Pradhan, D. Chatterjee, and J. Schaffner-Bielich, Frontiers in Astronomy and Space Sciences 9 (2022).
  * Huang  _et al._ (2023) C. Huang, G. Raaijmakers, A. L. Watts, L. Tolos, and C. Providência, (2023), arXiv:2303.17518 [astro-ph.HE] .
  * Char  _et al._ (2020) P. Char, S. Traversi, and G. Pagliara, Particles 3, 621 (2020).
  * Most  _et al._ (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018).
  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015).
  * Acernese  _et al._ (2015) F. Acernese _et al._ (VIRGO), Class. Quant. Grav. 32, 024001 (2015).
  * Abbott et al. (2019a) B. P. Abbott et al., Phys. Rev. X 9, 011001 (2019a).
  * Abbott et al. (2019b) B. P. Abbott et al., Astrophys. J. 882, L24 (2019b).
  * Abbott et al. (2017) B. P. Abbott et al., Phys. Rev. Lett. 119, 161101 (2017).
  * Abbott et al. (2018) B. P. Abbott et al., Phys. Rev. Lett. 121, 161101 (2018).
  * De  _et al._ (2018) S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Phys. Rev. Lett. 121, 091102 (2018).
  * Landry and Essick (2019) P. Landry and R. Essick, Phys. Rev. D 99, 084049 (2019).
  * Piekarewicz and Fattoyev (2019) J. Piekarewicz and F. J. Fattoyev, Phys. Rev. C 99, 045802 (2019).
  * Abbott et al. (2020) B. P. Abbott et al., Astrophys. J. 892, L3 (2020).
  * Dietrich  _et al._ (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020).
  * Punturo  _et al._ (2010) M. Punturo _et al._ , Class. Quant. Grav. 27, 194002 (2010).
  * Reitze  _et al._ (2019) D. Reitze _et al._ , Bull. Am. Astron. Soc. 51, 035 (2019).
  * Iacovelli  _et al._ (2023) F. Iacovelli, M. Mancarella, C. Mondal, A. Puecher, T. Dietrich, F. Gulminelli, M. Maggiore, and M. Oertel, Phys. Rev. D 108, 122006 (2023).
  * Huxford  _et al._ (2023) R. Huxford, R. Kashyap, S. Borhanian, A. Dhani, and B. S. Sathyaprakash, (2023), arXiv:2307.05376 [gr-qc] .
  * Rose  _et al._ (2023) H. Rose, N. Kunert, T. Dietrich, P. T. H. Pang, R. Smith, C. Van Den Broeck, S. Gandolfi, and I. Tews, Phys. Rev. C 108, 025811 (2023).
  * Walker  _et al._ (2024) K. Walker, R. Smith, E. Thrane, and D. J. Reardon, (2024), arXiv:2401.02604 [astro-ph.HE] .
  * Watts et al. (2016) A. L. Watts et al., Rev. Mod. Phys. 88, 021001 (2016).
  * Psaltis  _et al._ (2014) D. Psaltis, F. Özel, and D. Chakrabarty, Astrophys. J. 787, 136 (2014).
  * Zhang  _et al._ (2018) N.-B. Zhang, B.-A. Li, and J. Xu, Astrophys. J. 859, 90 (2018).
  * Lattimer and Steiner (2014a) J. M. Lattimer and A. W. Steiner, Eur. Phys. J. A 50, 40 (2014a).
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ , Astrophys. J. Lett. 848, L12 (2017a).
  * Antoniadis  _et al._ (2013) J. Antoniadis, P. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, _et al._ , Science 340, 6131 (2013).
  * Lattimer and Steiner (2014b) J. M. Lattimer and A. W. Steiner, Eur. Phys. J. A 50, 1 (2014b).
  * Miller et al. (2019) M. C. Miller et al., Astrophys. J. 887, L24 (2019).
  * Riley et al. (2019) T. Riley et al., Astrophys. J. Lett. 887, L21 (2019).
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021).
  * Riley et al. (2021) T. E. Riley et al., Astrophys. J. Lett. 918 (2021).
  * Huth  _et al._ (2022) S. Huth _et al._ , Nature 606, 276 (2022).
  * Abbott  _et al._ (2020) B. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020).
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (Virgo, Fermi-GBM, INTEGRAL, LIGO Scientific), Astrophys. J. 848, L13 (2017b).
  * Coughlin  _et al._ (2017) M. Coughlin, T. Dietrich, K. Kawaguchi, S. Smartt, C. Stubbs, and M. Ujevic, Astrophys. J. 849, 12 (2017).
  * Le Fèvre  _et al._ (2016) A. Le Fèvre, Y. Leifels, W. Reisdorf, J. Aichelin, and C. Hartnack, Nucl. Phys. A 945, 112 (2016).
  * Russotto and et. al. (2016) P. Russotto and et. al., Phys. Rev. C 94, 034608 (2016).
  * Hebeler  _et al._ (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013).
  * Tews  _et al._ (2013) I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 110, 032504 (2013).
  * Lynn  _et al._ (2016) J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis, K. E. Schmidt, and A. Schwenk, Phys. Rev. Lett. 116, 062501 (2016).
  * Drischler  _et al._ (2019) C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 122, 042501 (2019).
  * Drischler  _et al._ (2020) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125 (2020).
  * Huth  _et al._ (2021) S. Huth, C. Wellenhofer, and A. Schwenk, Phys. Rev. C 103 (2021).
  * Arzoumanian  _et al._ (2018) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Suppl. 235, 37 (2018).
  * Margalit and Metzger (2017) B. Margalit and B. D. Metzger, Astrophys. J. Lett. 850, L19 (2017).
  * Rezzolla  _et al._ (2018) L. Rezzolla, E. R. Most, and L. R. Weih, Astrophys. J. 852, L25 (2018).
  * Dietrich  _et al._ (2019) T. Dietrich, A. Samajdar, S. Khan, N. K. Johnson-McDaniel, R. Dudi, and W. Tichy, Phys. Rev. D 100, 044003 (2019).
  * Bulla (2019) M. Bulla, Mon. Not. Roy. Astron. Soc. 489, 5037 (2019).
  * Tews (2017) I. Tews, Phys. Rev. C 95, 015803 (2017).
  * Tsang  _et al._ (2023) C. Y. Tsang, M. B. Tsang, W. G. Lynch, R. Kumar, and C. J. Horowitz, (2023), arXiv:2310.11588 [nucl-th] .
  * Brown and Schwenk (2014) B. A. Brown and A. Schwenk, Phys. Rev. C 89, 011307 (2014).
  * Lynch and Tsang (2022) W. Lynch and M. Tsang, Physics Letters B 830, 137098 (2022).
  * Kortelainen  _et al._ (2012) M. Kortelainen, J. McDonnell, W. Nazarewicz, P.-G. Reinhard, J. Sarich, N. Schunck, M. V. Stoitsov, and S. M. Wild, Phys. Rev. C 85, 024304 (2012).
  * Danielewicz  _et al._ (2017) P. Danielewicz, P. Singh, and J. Lee, Nuclear Physics A 958, 147 (2017).
  * Zhang and Chen (2015) Z. Zhang and L.-W. Chen, Phys. Rev. C 92, 031301 (2015).
  * Adhikari  _et al._ (2021) D. Adhikari _et al._ (PREX), Phys. Rev. Lett. 126, 172502 (2021).
  * Reed  _et al._ (2021) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021).
  * Tsang  _et al._ (2009) M. B. Tsang, Y. Zhang, P. Danielewicz, M. Famiano, Z. Li, W. G. Lynch, and A. W. Steiner, Phys. Rev. Lett. 102, 122701 (2009).
  * Morfouace  _et al._ (2019) P. Morfouace _et al._ , Physics Letters B 799, 135045 (2019).
  * Estee  _et al._ (2021) J. Estee _et al._ ($\mathrm{S}\pi\mathrm{RIT}$ Collaboration), Phys. Rev. Lett. 126, 162701 (2021).
  * Cozma (2018) M. D. Cozma, Eur. Phys. J. A 54, 40 (2018).
  * Russotto  _et al._ (2011) P. Russotto, P. Wu, M. Zoric, M. Chartier, Y. Leifels, R. Lemmon, Q. Li, J. Łukasik, A. Pagano, P. Pawłowski, and W. Trautmann, Physics Letters B 697, 471 (2011).
  * Danielewicz  _et al._ (2002) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002).
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021).
  * Romani  _et al._ (2022) R. W. Romani, D. Kandel, A. V. Filippenko, T. G. Brink, and W. Zheng, The Astrophys. J. Lett. 934, L17 (2022).
  * Chabanat  _et al._ (1997) E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 627, 710 (1997).
  * Margueron  _et al._ (2002) J. Margueron, J. Navarro, and N. Van Giai, Phys. Rev. C 66, 014303 (2002).
  * Gómez  _et al._ (1992) J. Gómez, C. Prieto, and J. Navarro, Nuclear Physics A 549, 125 (1992).
  * Agrawal  _et al._ (2005) B. K. Agrawal, S. Shlomo, and V. K. Au, Phys. Rev. C 72, 014310 (2005).
  * Dutra  _et al._ (2014) M. Dutra, O. Lourenço, S. S. Avancini, B. V. Carlson, A. Delfino, D. P. Menezes, C. Providência, S. Typel, and J. R. Stone, Phys. Rev. C 90, 055203 (2014).
  * Malik  _et al._ (2023) T. Malik, M. Ferreira, M. B. Albino, and C. Providência, Phys. Rev. D 107, 103018 (2023).
  * Raithel  _et al._ (2018) C. Raithel, F. Özel, and D. Psaltis, Astrophys. J. Lett. 857, L23 (2018).
  * Malik and Pais (2022) T. Malik and H. Pais, Eur. Phys. J. A 58, 154 (2022).
  * Oppenheimer and Volkoff (1939) J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * Tolman (1939) R. C. Tolman, Phys. Rev. 55, 364 (1939).
  * Hinderer (2008) T. Hinderer, Astrophys. J. 677, 1216 (2008).
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Carriere  _et al._ (2003) J. Carriere, C. J. Horowitz, and J. Piekarewicz, Astrophys. J. 593, 463 (2003).
  * Patra  _et al._ (2022) N. K. Patra, S. M. A. Imam, B. K. Agrawal, A. Mukherjee, and T. Malik, Phys. Rev. D 106, 043024 (2022).
  * Malik  _et al._ (2024) T. Malik, H. Pais, and C. Providência, (2024), arXiv:2401.10842 [nucl-th] .
  * Tsang  _et al._ (2020) C. Y. Tsang, M. B. Tsang, P. Danielewicz, W. G. Lynch, and F. J. Fattoyev, Phys. Rev. C 102, 045808 (2020).
  * Thi  _et al._ (2021) H. D. Thi, C. Mondal, and F. Gulminelli, Universe 7, 373 (2021).
  * Patra et al. (2023) N. K. Patra et al., Phys. Rev. C 107, 055804 (2023).
  * Imam  _et al._ (2023) S. M. A. Imam, A. Mukherjee, B. K. Agrawal, and G. Banerjee, (2023), arXiv:2305.11007 [nucl-th] .
  * Patra  _et al._ (2023) N. K. Patra, P. Saxena, B. K. Agrawal, and T. K. Jha, Phys. Rev. D 108, 123015 (2023).
  * Pang  _et al._ (2023) P. T. H. Pang _et al._ , Nature Commun. 14, 8352 (2023).
  * Breschi  _et al._ (2024) M. Breschi, R. Gamba, G. Carullo, D. Godzieba, S. Bernuzzi, A. Perego, and D. Radice, (2024), arXiv:2401.03750 [astro-ph.HE] .
  * Ashton et al. (2019) G. Ashton et al., Astrophys. J. Suppl. Ser. 241, 27 (2019).
  * Romero-Shaw et al. (2020) I. M. Romero-Shaw et al., Mon. Not. Roy. Astron. Soc. 499, 3295 (2020).

## Appendix

Table A1:  Prior distributions of the NMPs for the Skyrme model. Here $\rho_{0}$, e0 and K0 are assumed as gaussian (G) distribution with median $\mu$ and standard deviation $\sigma$. Uniform (U) priors are used for the rest of the parameters with a minimum min and maximum max. The NMPs are listed in the units of MeV. Effective mass m∗ has the unit of nucleon mass, mN. |  $\rho_{0}$(G) | e0(G) | K0(G) | J0(U) | L0(U) | Ksym,0(U) | m∗(U)
---|---|---|---|---|---|---|---
min/$\mu$ | 0.16 | -16 | 230 | 24.7 | -11.4 | -328.5 | 0.5
max/$\sigma$ | 0.005 | 0.3 | 30 | 40.3 | 149.4 | 237.9 | 1
Table A2: The uniform prior is considered for the parameters of the RMF models (NL). Specifically, B and C are $b\times 10^{3}$ and $c\times 10^{3}$, respectively. It follows that “b” and “c” are going into Lagrangian density. The ’min’ and ’max’ entries denote the minimum and maximum values of the uniform distribution. All the parameters are dimensionless. |  $g_{\sigma}$ |  $g_{\omega}$ |  $g_{\varrho}$ |  B |  C |  $\xi$ |  $\Lambda_{\omega}$
---|---|---|---|---|---|---|---
min |  5.5 |  5.5 |  5.5 |  0.5 |  -10.0 |  0.0 |  0
max |  15.5 |  15.5 |  20.5 |  10.0 |  10.0 |  0.04 |  0.12