# Möbius strip microlasers: a testbed for non-Euclidian photonics.

Yalei Song  Laboratoire Lumière, Matière et Interfaces (LuMIn) CNRS, ENS
Paris-Saclay, Universié Paris-Saclay, CentraleSupélec, 91190 Gif-sur-Yvette
School of Physical Science and Technology, and Key Laboratory for Magnetism
and Magnetic Materials of MOE, Lanzhou University, Lanzhou, Gansu 730000,
China  Yann Monceaux  Laboratoire Lumière, Matière et Interfaces (LuMIn) CNRS,
ENS Paris-Saclay, Universié Paris-Saclay, CentraleSupélec, 91190 Gif-sur-
Yvette  Stefan Bittner  Chair in Photonics, LMOPS EA-4423 Laboratory,
CentraleSupélec and Université de Lorraine, 2 rue Edouard Belin, Metz 57070,
France  Kimhong Chao  Laboratoire Lumière, Matière et Interfaces (LuMIn) CNRS,
ENS Paris-Saclay, Universié Paris-Saclay, CentraleSupélec, 91190 Gif-sur-
Yvette  Héctor M. Reynoso de la Cruz  Science and Engineering Division of the
University of Guanajuato, Department of Physical Engineering, Academic Body of
Statistical Mechanics, León, Gto., México  Laboratoire Lumière, Matière et
Interfaces (LuMIn) CNRS, ENS Paris-Saclay, Universié Paris-Saclay,
CentraleSupélec, 91190 Gif-sur-Yvette  Clément Lafargue  Laboratoire Lumière,
Matière et Interfaces (LuMIn) CNRS, ENS Paris-Saclay, Universié Paris-Saclay,
CentraleSupélec, 91190 Gif-sur-Yvette  Dominique Decanini  Centre de
Nanosciences et de Nanotechnologies, CNRS, Université Paris-Saclay, 10
Boulevard Thomas Gobert, Palaiseau, 91120, France.  Barbara Dietz  School of
Physical Science and Technology, and Key Laboratory for Magnetism and Magnetic
Materials of MOE, Lanzhou University, Lanzhou, Gansu 730000, China  Joseph
Zyss  Laboratoire Lumière, Matière et Interfaces (LuMIn) CNRS, ENS Paris-
Saclay, Universié Paris-Saclay, CentraleSupélec, 91190 Gif-sur-Yvette  Alain
Grigis  Laboratoire d’Analyse, Géométrie et Applications, CNRS UMR 7539,
Université Sorbonne Paris Cité, Université Paris 13, Institut Galilée, 99
avenue Jean-Baptiste Clément, 93430 Villetaneuse, France  Xavier Checoury
Centre de Nanosciences et de Nanotechnologies, CNRS, Université Paris-Saclay,
10 Boulevard Thomas Gobert, Palaiseau, 91120, France.  Melanie Lebental
Laboratoire Lumière, Matière et Interfaces (LuMIn) CNRS, ENS Paris-Saclay,
Universié Paris-Saclay, CentraleSupélec, 91190 Gif-sur-Yvette  lebental@ens-
paris-saclay.fr

###### Abstract

We report on experiments with Möbius strip microlasers which were fabricated
with high optical quality by direct laser writing. A Möbius strip looks like a
ring with a twist and exhibits the fascinating property that it has only one
boundary and a one-sided nonorientable surface. Hence, in contrast to
conventional ring or disk resonators, a Möbius strip cavity cannot sustain
Whispering Gallery Modes (WGMs). Comparison between experiments and FDTD
simulations evidenced that the resonances are indeed not of whispering gallery
type, but localized along periodic geodesics. This finding is supported, on
the one hand, by an extension of the effective index approximation to curved
layers, and on the other hand, by an algorithm to systematically identify the
periodic geodesics.

###### pacs:

42.55.Sa, 03.65.Sq, 05.45.Mt, 42.60.Da

The emergence of three-dimensional (3D) optics is now enabled by the
availability of advanced printing technologies. It is in particular possible
to fabricate 3D surfaces with thicknesses down to 1 $\mu$m, allowing for light
to be confined and to propagate within a curved layer. Many questions then
arise regarding the structure of the electromagnetic field in non-Euclidean
photonics. It remains an emergent domain with only a few pioneering works prb-
goldstone ; peschel-PRL ; acoustic-3D ; arie-polariton . The present Letter
considers the question, how the propagation of light in such a thin but curved
structure influences the properties of microcavities and -lasers. The most
striking differences with respect to flat cavities is that light rays no
longer propagate along straight lines, but instead follow _geodesics_.
Actually we can assume that the Maupertuis principle is still valid and that
light follows the shortest paths. A geodesic is indeed the shortest path
between two points on a curved surface.

A Möbius strip is formed by connecting two ends of a strip after twisting one
end by 180∘, whereas connecting them without twist results in a simple ring
resonator. In spite of its simplicity, it exhibits peculiar topological
properties: it has a single, non-orientable surface and a single boundary.
Hence Möbius strips have captured the attention of generations of scientists
to this day, from pure mathematics to physical, chemical and engineering
sciences nature-cristal ; isolant-topo ; prb-electron ; lpr-plasmon ; epl-
hentschel . In particular, a Möbius microcavity cannot display whispering
gallery modes (WGMs) because its boundary possesses both concave and convex
parts, in distinction to conventional ring resonators. The central question we
address is what kind of structure its lasing modes exhibit and how it is
related to its non-Euclidean geometry.

Microlasers provide an ideal testbed for such studies, since we can lean on a
deep experimental and theoretical background expertise which has been built
over more than two decades in the domain of planar microlasers. Their typical
size ranges from 50 to 200 $\mu$m for wavelengths around 1 $\mu$m to get
enough gain, while maintaining a high spectral resolution. The microlasers
then operate in the semiclassical regime brack ; stoeckmann-livre . In
particular, for flat layers, the resonance spectrum is sustained by periodic
orbits, which are classical periodic ray trajectories. In a ring cavity, the
modes with the highest quality factor are WGMs which propagate along the outer
circular boundary. Their corresponding ray periodic orbits are regular
polygons with reflections at the outer circle. Such classical trajectories do
not exist in a Möbius strip. In this letter, we show that the lasing modes of
a Möbius strip are located along periodic geodesics acoustic-3D ;
colindeverdiere ; sebbah .

Figure 1: SEM image of a Möbius strip microlaser with $R=50~{}\mu$m,
$W=15~{}\mu$m and $h=3~{}\mu$m.

##  I Möbius fabrication and parametrization

The 3D design of the Möbius strip was drawn with the OpenSCAD software by
rotating a rectangle by 180∘ while its center travels along a circle of radius
$R$. The sidelengths of the rectangle are $W$, which is fixed to 15 $\mu$m,
and $h$, which corresponds to the Möbius strip thickness, and varies from 1 to
5 $\mu$m. The radius $R$ varies from 40 to 60 $\mu$m. The notations are
reported in Fig. 1.
The microlasers were then fabricated by Direct Laser Writing (DLW) lithography
using a Photonic Professional GT system with negative resist IP-G780 both
developed and manufactured by the company Nanoscribe 2PP . The resist was
doped by $0.5$ wt% pyrromethene 597 laser dye (by the company Exciton) which
is homogeneously distributed in the bulk host resist pyramid . Each cavity is
supported by six circular pylons to avoid coupling with the glass substrate. A
Scanning Electronic Microscope (SEM) image of a Möbius strip microlaser is
shown in Fig. 1.
A standard parametrization of a Möbius strip surface
$\vec{r}=\vec{F}(\varphi,w)$ is

$\displaystyle x$ $\displaystyle=$ $\displaystyle
F_{x}(\varphi,w)=(R+w\cos\frac{\varphi}{2})\,\cos\varphi$ (1) $\displaystyle
y$ $\displaystyle=$ $\displaystyle
F_{y}(\varphi,w)=(R+w\cos\frac{\varphi}{2})\,\sin\varphi$ (2) $\displaystyle
z$ $\displaystyle=$ $\displaystyle F_{z}(\varphi,w)=w\sin\frac{\varphi}{2}$
(3)

with $w\in[-W/2;W/2]$ and $\varphi\in[0;2\pi]$, which looks similar to the SEM
image of Fig. 1 for vanishing thickness $h$. A Möbius strip can be twisted the
other way, which corresponds to a minus on line (2). We checked that the
experiments and their interpretation were consistent for both chiralities.
Note that this parametrization differs from a twisted paper sheet, which can
be unfolded on a flat plane, and is then called a _flat_ Möbius strip, see
Appendix A. This property leads to important consequences on geodesics and
mode localization, which will be discussed below.

The Möbius strip microlasers were pumped individually and with uniform
intensity by a beam perpendicular to the substrate from a frequency doubled
Nd:YAG laser (532 nm, 500 ps, 10 Hz). Their emission was analyzed by a
spectrometer coupled to a CCD camera with an overall resolution of 0.03 nm.
The experiments were carried out at room pressure and temperature.
The microlaser light was mostly emitted parallel to the substrate plane. The
inset of Fig. 2 shows a laser threshold at 5 MW/cm2. As expected, this value
decreases with increasing size of the Möbius strip and is similar to the laser
threshold of 3D cavities made of the same laser dye with an equivalent gain
volume pyramid .
The photograph of Fig. 3(b) shows a Möbius strip lasing, whose emission mainly
originates from the lateral boundaries of the strip at a near-grazing angle.
This is observed from all directions, which is consistent with modes confined
inside the Möbius strip and propagating within this curved and twisted guiding
layer.
A laser emission spectrum is shown in Fig. 2(a). It does not significantly
depend on the direction of observation, but depends on the thickness $h$, and
becomes comb-like for $h=1\mu$m. Actually, the Möbius strip behaves like a
guiding layer and an in depth analysis of the spectrum yields information on
the mode localization.

Figure 2: (a) Experimental laser spectrum of a Möbius strip with $R=50$ $\mu$m
and $h=1$ $\mu$m. (b) Normalized Fourier transform of the spectrum in (a)
considered as intensity versus wavenumber. Inset: Experimental laser threshold
curve for a Möbius microcavity with $R=50$ $\mu$m and $h=3$ $\mu m$. Figure 3:
Camera image of a Möbius strip illuminated with a low intensity white light
(a) without pumping, and (b) under pumping. We can see the cavity and its
reflection on the glass substrate. The resist IP-G780 is transparent. The
green pump light is removed by a notch filter, while the laser emission is
yellow. (c) Effective refractive index for a flat layer (Eq. (4) in Ref. pra-
simple ) and a curved layer for $n=1.515$, an outer radius $b=50~{}\mu$m, and
an inner radius $a=49~{}\mu$m in the experimental range $\lambda=2\pi/k=560$
nm to $\lambda=600$ nm.

##  II Helmholtz equation

In Cartesian coordinates, the electric field fulfills the Helmholtz equation

$(\Delta+n^{2}k^{2})\vec{E}=\vec{0}$ (4)

where $k$ is the wavenumber in vacuum and $n$ the refractive index, which is 1
outside of the strip and $n\simeq 1.5$ inside it. As the fabrication method
does not generate stress in the resist, $n$ is assumed to be homogeneous
within the strip. Although it is experimentally evidenced that the laser
emission of a Möbius microcavity is polarized, for the sake of simplicity, we
will only consider a single electric field component $\psi$.
Following Ref. dacosta , the Cartesian coordinate system is changed to a local
coordinate system ($q_{1}$, $q_{2}$, $q_{3}$) defined by the local frame
($\vec{u}_{1}$, $\vec{u}_{2}$, $\vec{u}_{3}$) with

$\displaystyle\vec{u}_{1}$ $\displaystyle=$
$\displaystyle\frac{\partial\vec{F}}{\partial\varphi}$ (5)
$\displaystyle\vec{u}_{2}$ $\displaystyle=$
$\displaystyle\frac{\partial\vec{F}}{\partial w}$ (6)
$\displaystyle\vec{u}_{3}$ $\displaystyle=$
$\displaystyle\frac{\vec{u}_{1}\times\vec{u}_{2}}{||\vec{u}_{1}\times\vec{u}_{2}||}$
(7)

By construction, $\vec{u}_{1}$ and $\vec{u}_{2}$ lie within the strip, whereas
$\vec{u}_{3}$ is orthogonal to it. For the Möbius strip considered here
$(\varphi,w)=(q_{1},q_{2})$. The Laplace operator $\Delta$ in a curvilinear
coordinate system is given by

$\Delta\psi=\frac{1}{\sqrt{g}}\sum_{i,j=1}^{3}\frac{\partial}{\partial
q_{i}}\left(\sqrt{g}\,g^{ij}\frac{\partial\psi}{\partial q_{j}}\right)$ (8)

with $g_{ij}=(\vec{u}_{i}\,\cdot\,\vec{u}_{j})$ the metric tensor, $(g^{ij})$
its inverse, and $g=\det(g_{ij})$. From Eqs. (5-7) follow $g_{33}=1$ and
$g_{13}=g_{31}=g_{23}=g_{32}=0$. Consequently the Laplace operator within the
strip $\Delta_{s}$ can be seperated from the Laplacian perpendicular to it,

$\Delta_{s}\psi+\frac{\partial^{2}\psi}{\partial
q_{3}^{2}}+n^{2}k^{2}\psi=0\,.$ (9)

Using the separation Ansatz $\psi=\psi_{s}(q_{1},q_{2})\,\psi_{3}(q_{3})$ and
assuming that $n$ depends only on $q_{3}$, the Helmholtz equation can be split
into the vertical direction

$\frac{\partial^{2}\psi_{3}}{\partial
q_{3}^{2}}+(n^{2}-n^{2}_{\mathrm{eff}})\,k^{2}\psi_{3}=0$ (10)

and the propagation within the surface of the strip:

$\Delta_{s}\psi_{s}+n_{\mathrm{eff}}^{2}k^{2}\psi_{s}=0$ (11)

The effective index $n_{\mathrm{eff}}$ was originally introduced for flat
layers. In the ($q_{1}$, $q_{3}$) plane, a section of the Möbius strip locally
looks like a ring. Following Ref. jap-indice-courbe , we extend the effective
index approach to a curved layer, which is assumed to be a segment of a ring.
The full derivation is detailed in Ref. anneau . As shown in Fig. 3(c), the
agreement of the effective index for a flat and a curved layer is excellent in
the range of $k$ corresponding to the experiments. This approximation allows
us to investigate the propagation of light in the surface defined by
$(\vec{u}_{1},\vec{u}_{2})$, without considering the vertical $\vec{u}_{3}$
direction. In a geometric optics picture, this corresponds to rays traveling
along the curved surface of the Möbius strip. Since the radius $R$ and the
width $W$ of the strip are large enough compared to the wavelength, a
semiclassical approach is valid. Hence we surmise that the resonant modes are
related to periodic geodesics.

Figure 4: Examples of periodic geodesics in a Möbius strip with $R=50~{}\mu$m.
(a, b) This periodic geoedesic has a length $L=7.03R$ and each consecutive
reflection is on the opposite side. (c, d) This periodic geoedesic has a
length $L=6.73R$. The rotation symmetry axis $y=z=0$ is indicated as green
dashed lines in panels (a, c) and as green dots in panels (b, d). The red
numbers are the indices of the vertices. In (d), the blue line is the mean
position calculated from Eq. (15) for the wave function in Fig. 6.

##  III Geodesics

A classical trajectory $\gamma$ on a surface can be represented by a path
$q(s)$ with the arclength $s\in[0,1]$; $q(0)$ is the initial point and $q(1)$
the final point. This trajectory $\gamma$ is a geodesic if its length is
stationary with respect to variations of the path $q(s)$ while maintaining
$q(0)$ and $q(1)$. The length is given by Kuehnel

$L_{\gamma}=\int_{0}^{1}ds\left[\sum\limits_{i,j}g_{ij}[q(s)]\,\frac{dq_{i}}{ds}(s)\,\frac{dq_{j}}{ds}(s)\right]^{1/2}$
(12)

with $i,j=1,2$ and $q=(q_{1},q_{2})$. A periodic geodesic is such that
$q(0)=q(1)$. Since the geodesics on the Möbius strip cannot be derived
analytically, we calculated them numerically using this variational principle
(see Appendix B.2).
In contrast to the flat Möbius billiard, which is described in Appendix A,
each periodic geodesic has at least one reflection on the boundary, and is
isolated111A periodic geodesic is called _isolated_ if there is no periodic
geodesic in its vicinity.. Another striking difference to the flat Möbius
billiard is the sequence of reflections on the boundary. In the flat Möbius
billiard, the geodesic bounces alternatively opposite sides222The notion of
same or opposite side makes sense from a local perspective only, since a
Möbius strip actually has only one side., namely $w=+W/2$, then $w=-W/2$, then
$w=+W/2$, etc.. Such periodic geodesics exist in the Möbius strip as well, but
only if they have at least 9 reflections, see Fig. 4(a) and Appendix C. Figure
4(c) shows a different kind of periodic geodesics which exhibits also
consecutive reflections on the same side, a feature specific to the 3D Möbius
strip. It experiences 1 revolution and 5 bounces in one period. We call it
_geodesic 5a_ to distinguish it from other periodic geodesics with 5 bounces
and similar length. Hereafter we will compare it to experiments and numerical
simulations.

##  IV Spectrum analysis

The comb-like spectrum of Fig. 2(a) indicates that the lasing modes are
longitudinal modes located along a single classical trajectory (or a few
trajectories with similar lengths). To verify this conjecture, we analyze the
spectra.
The free spectral range or distance between two adjacent peaks $\Delta k$ is
inversely proportional to the geometric length $L$ of this classical
trajectory APL-FP3D

$\Delta k=\frac{2\pi}{n_{g}L}$ (13)

where $n_{g}$ is the group refractive index. $\Delta k$ can be inferred
systematically from the Fourier transform, as illustrated in Fig. 2(b), which
is peaked at the position $n_{g}L$ and its harmonics. Hence we can determine
the geometric length $L$ of the trajectory, if the group refractive index
$n_{g}$ is known. For 3D surfaces like the Möbius strip, the group refractive
index is deduced from the effective index and contributions from the material
dispersion of $n$ and modal dispersions due to $n_{\mathrm{eff}}$
approximation:

$n_{g}=n_{\mathrm{eff}}+k\frac{\partial n}{\partial
k}|_{\mathrm{mat}}+k\frac{\partial n}{\partial k}|_{\mathrm{eff}}$ (14)

The bulk refractive index is n=1.515 $\pm$ 0.005 indice . Its dispersion is
low, thus we assume that $n$ is constant. But $k$ is large, hence
$k\frac{\partial n}{\partial k}|_{\mathrm{mat}}$ cannot be neglected.
Following Ref. APL-FP3D , it is inferred from the laser spectra of cuboids,
since their lasing modes are located on a well-known periodic orbit in the
bulk of the cavity. From these data, the term due to the material dispersion
is estimated as $k\frac{\partial n}{\partial k}|_{\mathrm{mat}}=$0.05 $\pm$
0.04.
The effective index and its dispersion are calculated, as discussed above. In
the frequency range of interest and for a 1 $\mu$m thick strip, there are four
branches (vertical excitations) for each polarization with significantly
different index and dispersion values. Since the spectrum in Fig. 2(a)
exhibits a single comb, we consider only the branch with the highest effective
index, which actually has the best confinement factor, and thus the highest
laser gain. The corresponding effective index is $n_{\mathrm{eff}}=1.49$ and
its dispersion equals $k\frac{\partial n}{\partial
k}|_{\mathrm{eff}}=0.037\,\pm\,0.002$. Finally, taking into account the three
terms of Eq. (14) yields a group index of $n_{g}=1.58\,\pm\,0.05$.

Figure 5: Comparison of experimental optical lengths for different Möbius
microlasers ($R=50\,\mu$m, $h=1\,\mu$m) (red circles) with the expected
optical length of geodesic 5a (continuous blue lines) and the half perimeter
(dotted black lines). The lower (resp. upper) line corresponds to $n_{g}=1.53$
(resp. $n_{g}=1.63$), which is the lower (resp. upper) limit of the group
index $n_{g}$. Inset: The dashed line indicates the perimeter.

Dozens of Möbius microlasers with various radii and thicknesses were
investigated. For illustration, Fig. 5 shows the measured optical lengths for
Möbius microlasers with radius $R=50~{}\mu$m and thickness $h=1~{}\mu$m
compared with half the perimeter. The obvious disagreement evidences that the
lasing modes do not propagate close to the cavity boundary, as WGMs would do.
On the contrary, the agreement is very good for the geodesic 5a. A few
geodesics are also consistent with experiments. All the corresponding rays are
confined in the Möbius strip by total internal reflection. Some geodesics can
be excluded because their length is too short or too long. This includes for
instance all the geodesics with consecutive reflections on opposite sides,
like the one in Fig. 4(a).

##  V Wave functions

Numerical simulations were performed with a home-made 3D Finite Difference
Time Domain (FDTD) code for $R=10\,\mu$m, $W=3\,\mu$m, and $h=150$ nm. The
aspect ratio $W/R$ is the same as in experiments with $R=50\,\mu$m and
$W=15\,\mu$m. For such parameters, there is a single vertical excitation in
the thickness direction ($\vec{u}_{3}$) but there exist several modes in the
transverse $w$ direction. Some wave functions are presented in Appendix D. A
typical wave function is plotted in Fig. 6. It exhibits maximal intensity
close to the boundary at positions coinciding with the vertices of the
geodesic 5a which is superimposed. To quantify the agreement, the mean
position of the wave function along the strip width is calculated for each
$\varphi$ value:

$\langle w\rangle(\varphi)=\frac{\iint\limits_{\textrm{Half-
plane}}w\,\rho(w,\varphi,q_{3})\,dw\,dq_{3}}{\iint\limits_{\textrm{Half-
plane}}\,\rho(w,\varphi,q_{3})\,dw\,dq_{3}}$ (15)

where
$\rho=\frac{1}{2}\epsilon_{0}\epsilon_{r}|\vec{E}|^{2}+\frac{1}{2}\mu_{0}|\vec{H}|^{2}$
is the energy density of the electromagnetic field for a given wave function.
This mean position is plotted in Fig. 4(d). The agreement is very good, in
particular for the reflections on the boundary, further supporting the
hypothesis that the modes are localized on periodic geodesics.

Figure 6: Different views of the same wave function with
$\mathrm{Re}(kR)=88.0$ and $\mathrm{Im}(kR)=-2.18\cdot 10^{-4}$. Plots of
$|\vec{E}|^{2}$. The geodesic 5a is superimposed as red lines.

##  VI Conclusion

Organic Möbius strip microlasers are investigated as non-trivial examples of
non-Euclidean photonic structures. Actually direct laser writing (DLW)
techniques enable the fabrication of almost arbitrary 3D photonic polymer
structures with sub-micron feature sizes. We evidence that Möbius strip
microlasers do not exhibit WGMs, contrary to conventional ring cavities.
Instead their lasing modes are located on periodic geodesics. These findings
are based on 3D FDTD numerical simulations, and on a dedicated algorithm to
systematically identify the periodic geodesics. They are substantiated by a
step-by-step derivation of the wave propagation in a curved layer using the
effective index approximation, which implies reduction of the vectorial
electric field to a scalar wave function. An objective for the future will be
the derivation of a vectorial equation for the electromagnetic wave
propagation in curved surfaces and an in-depth investigation of the non-
trivial polarization features of various surface microlasers. This work opens
the way to the exploration of non-Euclidian photonic devices.

###### Acknowledgements.

This work was supported by the French RENATECH network, the CNano IdF DIM
Nano-K, the Labex NanoSaclay (ANR-10-LABX-0035), and the Laboratoire
International Associé ImagiNano. Y. S. and B. D. are supported by the National
Natural Science Foundation of China under Grants No. 11775100 and No.
11775101. S. B. acknowledges support for the Chair in Photonics from Ministère
de l’Enseignement Supérieur, de la Recherche et de l’Innovation, Région Grand-
Est, Département Moselle, European Regional Development Fund (ERDF), Metz
Métropole, GDI Simulation, CentraleSupélec, and Fondation CentraleSupélec. H.
M. R is supported by Consejo Nacional de Ciencia y Tecnología - México under
Grant No. 2019-000016-01NACF-03530 of Movilidad Extranjera 2019, and the
Dirección de Apoyo a la Investigación y el Posgrado (DAIP) Guanajuato under
the Convocatoria de apoyo a posgrados-2019 Formato-B.
J.-F. Audibert and H. K. Warner are acknowledged for the fabrication of
centimeter-scale Möbius strips with a 3D printer, N. Pavloff for pointing Ref.
prb-goldstone , P. Pansu, F. Jean, and G. Bossard for advices on geodesics,
and E. Bogomolny for his derivation of the effective index approximation.

##  Appendix A Periodic orbits of the flat Möbius billiard

The conventional Moebius strip parametrized by the set of Eqs. (1-3) is
hereafter referred to as _3D Möbius strip_ to distinguish it from the _flat
Möbius billiard_. Please note that all calculations presented in the
appendices are made for a Möbius strip with chirality corresponding to an
additional minus sign in Eq. (2) for the $y$ coordinate.

Figure 7:  Periodic orbits in the flat Möbius billiard. (a) Geometry of the
flat Möbius billiard of length $D$ and width $W$. The boundaries at $x=0$ and
$x=D$ are identified with each other with an inversion of the $y$ direction,
hence the points marked by o and x are identified with each other,
respectively. (b) The isolated periodic orbit $(n_{t},n_{l})=(0,1)$ (solid
line) and one PO of the $(0,2)$ family (dashed line). (c) Two examples of the
$(1,1)$ PO family. (d) Two examples of the $(3,1)$ family. (e) Two examples of
the $(4,2)$ family. (f) Two examples of the $(1,3)$ PO family. The numbers
$1,2,3$ in panels (b), (e) and (f) indicate the connections of the different
PO segments.

In this section we discuss the flat Möbius billiard333Note that the flat
Möbius billiard is not equivalent to a developable Möbius strip Gravesen .
which serves as a two-dimensional analog of the actual 3D Möbius strip. The
comparison of their periodic orbits / geodesics demonstrates how the curvature
and torsion of the Möbius strip change the possible structures of periodic
trajectories. The flat Möbius billiard is defined as a rectangle with length
$D=2\pi R$ and width $W$ where the transverse boundaries, defined by
$x\in[0,D]$ at $y=\pm W/2$, are hard walls with specular reflections. The
diagonally opposite corners are identified with each other, that is, periodic
boundary conditions in $x$ direction with an inversion of the $y$ coordinate
are imposed to obtain the topology of a Möbius strip, that is

$\begin{array}[]{lcl}x=D&\rightarrow&x=0\\\ y|_{x=0}&=&-y|_{x=D}\\\
\dot{y}|_{x=0}&=&-\dot{y}|_{x=D}\end{array}$ (16)

as illustrated in Fig. 7(a).

The periodic orbits (POs) of the flat Möbius billiard are classified by the
index pair $(n_{t},n_{l})$ where $n_{t}$ is the number of reflections at the
transverse boundaries and $n_{l}$ is the number of roundtrips in longitudinal
($x$) direction. Several examples are shown in Fig. 7. Since $\dot{y}$, the
transverse component of the velocity vector, must return to its initial value
after one roundtrip along a PO, and each reflection as well as each roundtrip
in longitudinal direction change the sign of $\dot{y}$, the sum of $n_{t}$ and
$n_{l}$ must be even. The only exception is the $(0,1)$ periodic orbit (and
its odd repetitions) because $\dot{y}=0$ in this case [see Fig. 7(b)].

All POs are part of a family that covers the whole billiard, with the
exception of the $(0,1)$ orbit which is isolated. The lengths of the POs are
given by

$L(n_{t},n_{l})=[(n_{t}W)^{2}+(n_{l}D)^{2}]^{1/2}$ (17)

and their momentum vectors $\vec{k}=(k_{x},k_{y})$ fulfill
$|k_{x}|=k\cos(\alpha)$ and $|k_{y}|=k\sin(\alpha)$ with

$\alpha(n_{t},n_{l})=\arctan[n_{t}W/(n_{l}D)]\,.$ (18)

Any trajectory with a momentum vector fulfilling Eq. (18) is a PO, independent
of its initial position.

##  Appendix B Calculation of geodesics

###  B.1 Local approach

Figure 8:  Metric tensor element $g_{11}$ for the 3D Möbius strip as function
of the generalized coordinates $q=(q_{1},q_{2})$.

A curve $\gamma(s)$ along $\vec{r}\,[q_{1}(s),q_{2}(s)]$ on a curved surface
parametrized by, e.g., the arclength $s$, is a geodesic if it fulfills

$\frac{d^{2}q^{k}}{ds^{2}}+\sum\limits_{i,j}\Gamma_{ij}^{k}\frac{dq^{i}}{ds}\frac{dq^{j}}{ds}=0$
(19)

with $i,j,k=1,2$ (see Eq. 7.19b in Ref. Nakahara ). The Christoffel symbols
are given by

$\Gamma_{ij}^{k}=\frac{1}{2}\sum\limits_{m}g^{km}\left(\frac{\partial
g_{mi}}{\partial q^{j}}+\frac{\partial g_{mj}}{\partial q^{i}}-\frac{\partial
g_{ij}}{\partial q^{m}}\right)\,$ (20)

where $(g_{ij})$ is the metric tensor and $(g^{ij})$ its inverse. For the 3D
Möbius strip

$g_{11}=R^{2}+\frac{q_{2}^{2}}{4}\,[3+2\cos(q_{1})]+2R\,q_{2}\cos\left(\frac{q_{1}}{2}\right)\,$
(21)

with the coordinates $q=(q_{1},q_{2})$ parametrizing the surface of a 3D
Möbius strip as defined in Eqs. (1-3). The other tensor elements are
$g_{22}=1$ and $g_{12}=g_{21}=0$. The element $g_{11}$ is shown as function of
$(q_{1},q_{2})$ in Fig. 8.

Solving this initial conditions problem evidently requires knowledge of the
initial values of $q(s)$ and of $\dot{q}(s)=dq(s)/ds$. For the Möbius strip
the Eqs. (19) must be solved numerically with an algorithm such as the fourth-
order Runge-Kutta method.

Finding periodic geodesics by solving an initial conditions problem is
difficult since it is a priori unknown which initial position and velocity
lead to a periodic trajectory. Hence, a large parameter space needs to be
explored with high resolution. Moreover, it is difficult to judge whether a
numerical solution of these equations is indeed a periodic geodesic, in view
of the unavoidable numerical inaccuracies. A variational approach is simpler
and more efficient.

###  B.2 Variational approach

Since a geodesic is the locally shortest path between two points — which is
essentially Fermat’s principle — we can find a geodesic segment connecting two
points with coordinates $Q^{(0)}=q(0)$ and $Q^{(1)}=q(1)$ by searching for a
trajectory $\gamma(s)$ for which the length $L_{\gamma}$ given by

$L_{\gamma}=\int_{0}^{1}ds\left[\sum\limits_{i,j}g_{ij}[q(s)]\,\frac{dq_{i}}{ds}(s)\,\frac{dq_{j}}{ds}(s)\right]^{1/2}$
(22)

is stationary. This formulation as a boundary condition problem no longer
requires a priori knowledge of the initial velocity $\dot{q}$. However, the
two vertices $Q^{(0,1)}$ must be given beforehand and are fixed. Constructing
a periodic geodesic requires calculating several such segments that are
connected in a closed loop with vertices on the boundary of the Möbius strip.
Then all the vertices, including the initially fixed vertices $Q^{(0,1)}$, are
varied. This will be described in more detail further below.

We start with the algorithm to calculate a single geodesic segment. Since an
analytic solution of this infinitely-dimensional variational problem is not
possible for the 3D Möbius strip, we solve it numerically using a discretized
trajectory $\gamma=\\{q(n)\\}$ with $n=1\dots N$ where $q(1)$ corresponds to
$Q^{(0)}$ and $q(N)$ to $Q^{(1)}$. The integral for $L_{\gamma}$ is then
approximated by a sum over the $N$ discretization points using the trapezoidal
rule where the derivatives of $q(s)$ in Eq. (22) are replaced by central
finite differences. The trajectory $\gamma$ is stationary when all the
derivatives $\partial L_{\gamma}/\partial q_{i}(n)$ vanish. Here, variation is
performed perpendicular to the curve $\gamma(s)$, implying that the two
coordinates $[q_{1}(n),q_{2}(n)]$ at discretization point $n$ are not varied
independently. For this, the derivative $v_{n}=\partial L_{\gamma}/\partial
q_{\perp}(n)$ is computed in the direction $q_{\perp}(n)$ which is
perpendicular to the trajectory in $q$ space at discretization point $n$. Note
that a variation tangential to the trajectory would correspond to a
perturbation along the trajectory which implies a change in $s$. Therefore, as
is general practice, we do not consider this variation.

The trajectory segment is initialized with a straight line in $q$ space. The
correct solution for which $\vec{v}$ vanishes is found using a Newton
iteration, that is, $q(n)\rightarrow q(n)-d_{n}\,q_{\perp}(n)$ where $\vec{d}$
– a $N$ component vector – is the solution of $H\vec{d}=\vec{v}$. The matrix
elements of $H$ are

$H_{mn}=\frac{\partial^{2}L_{\gamma}}{\partial q_{\perp}(m)\,\partial
q_{\perp}(n)}\,,$ (23)

so $\vec{v}$ and $H$ are the first and second derivatives (gradient and
Hessian), respectively, of $L_{\gamma}$ with respect to the $q_{\perp}(n)$.
Since the first and last point of the segment are kept fixed, only the
elements with indices in the range $2\dots N-1$ of vectors $\vec{v}$ and
$\vec{d}$ and matrix $H$ are considered in the above equations. The Newton
iteration is continued till $\vec{v}$ converges to $\vec{0}$ within machine
precision.

In order to find a periodic geodesic consisting of $M$ segments, $M+1$
vertices $Q^{(m)}$ on the boundary of the Möbius strip are chosen with values
$Q_{2}^{(m)}=\pm W/2$ selected to correspond to the symbolic code of interest
and the $Q_{1}^{(m)}$ chosen randomly in ascending order. In a first step, the
$M$ individual segments are iterated as described above to find the geodesics
connecting the different vertices which are kept fixed. In a second step, one
step of the Newton iteration for all discretization points including all of
the vertices $Q^{(m)}$ is calculated, where the vertices are only allowed to
vary in the $q_{1}$ direction to keep them on the boundary of the Möbius
strip. Steps one and two are repeated till the trajectory becomes stationary
within machine precision. To ensure periodicity, $Q^{(0)}$ and $Q^{(M)}$ are
initialized such that they represent the same point in configuration space,
that is, $Q_{1}^{(M)}=Q_{1}^{(0)}+2\pi n_{l}$ and
$Q_{2}^{(M)}=(-1)^{n_{l}}Q_{2}^{(0)}$ where $n_{l}$ is the number of
revolutions around the Möbius strip, and the same variation is applied to both
of them.

A large number of random initial conditions for the vertices $Q^{(m)}$ was
tested for each considered symbolic code since the algorithm only converges to
a physical trajectory if the initial conditions are sufficiently close to a
correct solution. A vanishing gradient $\vec{v}$ proves that the calculated
trajectories are numerically exact geodesics. Typically, $10,000$
discretization points per revolution were used in the calculations, and it was
checked that increasing the number of discretization points did not
significantly change the length $L_{\gamma}$ or the shape of the trajectory.

Furthermore, it was verified that the law of elastic reflection is respected
at all vertices, i.e., that the ingoing angle is equal to the outgoing angle.
It should be noted that the law of elastic reflection is not a priori assumed
in any way during the numerical calculations of the periodic geodesics but
results when a correct solution is found. In addition, several periodic
geodesics were verified independently by solving the equations (19) with the
correct initial conditions.

##  Appendix C Periodic geodesics of the 3D Möbius strip

Figure 9:  Further examples of periodic geodesics (red solid lines) in
configuration (a, c, e) and $q$ space (b, d, f). The dashed green line
indicates the symmetry axis. (a, b) Periodic geodesic $(3,1,sos)$ with length
$L=6.53R$ on a Möbius strip with aspect ratio $W/R=0.3$. (c, d) Periodic
geodesic $(4,1,sos^{2})$, called 4a in the following, with length $L=6.66R$ on
a Möbius strip with $W/R=0.3$. (e, f) Periodic geodesic $(5,1,so^{3}s)$ with
length $L=6.67R$ on a Möbius strip with $W/R=0.375$. Its symmetry partner
obtained by rotation around the axis $y=z=0$ (green dashed line) is indicated
as cyan dash-dotted line.

In this Appendix, further examples of periodic geodesics of the 3D Möbius
strip are shown and their structure is compared to the periodic orbits in the
flat Möbius billiard and conventional ring cavities.

Figure 9 shows three typical examples. Some periodic geodesics are symmetric
with respect to a rotation of $180^{\circ}$ around the symmetry axis $x=z=0$
of the Möbius strip, indicated as green dashed line in the configuration space
plots, such as the ones in Figs. 9(a-d). These trajectories all pass through
the point $Q=(\pi,0)$ indicated as green dot in the $q$ space plots. Other
orbits such as the one shown in Figs. 9(e, f) are not, and they become a
distinct periodic geodesic under the symmetry rotation.

The periodic geodesics are labeled by indices $(n_{t},n_{l})$ and a symbolic
code, where $n_{t}$ is the number of vertices at the boundary and $n_{l}$ is
the number of revolutions around the Möbius strip. The symbolic code consists
of $n_{t}$ letters o (s) which designate trajectory segments that connect
vertices with a different (the same) sign of $q_{2}$. This means that the
segment appears to connect opposite boundaries (the same boundary), though the
notion of ”same” and ”opposite” side only makes sense from a local perspective
given that the Möbius strip actually has only one boundary444The labels are
not unambiguous since several distinct periodic geodesics with the same
indices and code can exist.. The $(11,1,o^{11})$ orbit in Figs. 4(a, b), for
example, has a zig-zag structure where all segments connect opposite sides,
whereas the $(5,1,s^{2}os^{2})$ orbit in Figs. 4(c, d) starts with two
segments connecting the same side, followed by a segment connecting to the
opposite side after half a roundtrip followed by two segments connecting the
same side again.

The existence of specific periodic geodesics depends on the aspect ratio
$W/R$. For example, the $(5,1,so^{3}s)$ geodesic in Figs. 9(e, f) does not
exist for smaller aspect ratios such as $W/R=0.3$. As the aspect ratio is
decreased and the strip becomes narrower, the segment connecting vertices $1$
and $2$ eventually touches the lower boundary $q_{2}=-W/2$ and the geodesic
hence ceases to exist. In contrast, the periodic orbits of the flat Möbius
billiard exist independently of its aspect ratio.

The three periodic geodesics in Fig. 9 all have structures that cannot be
found for periodic orbits of the flat Möbius billiard. First, all of them
feature segments that connect to the same boundary (code $s$), however, such
segments are impossible in the flat Möbius billiard. Furthermore, some
periodic geodesics such as the $(4,1,sos^{2})$ in Figs. 9(c, d) have
$n_{t}+n_{l}$ odd, whereas it must always be even for the flat Möbius. These
differences can be directly attributed to the curvature and torsion of the 3D
Möbius strip. First, due to the curvature, some parts of the boundary are
concave and thus permit $s$ trajectory segments, whereas the straight boundary
of the flat Möbius billiard does not allow these. Second, $n_{t}+n_{l}$ must
be even in the flat Möbius since $\dot{y}$ can only change its sign at a
reflection but stays constant in between, whereas $\dot{q}_{2}$ changes
continuously during the propagation on the curved surface of the 3D Möbius
strip in addition to the sign changes at the vertices.

Another difference is that the POs of the flat Möbius billiard are all part of
a continuous family whereas the periodic geodesics of the 3D Möbius strip are
all isolated. The existence of families of POs is due to the continuous
translational symmetry in $x$ direction of the flat Möbius billiard. While the
3D Möbius strip has a similar periodic boundary condition as the flat Möbius
billiard, the metric tensor $g_{ij}$ does not exhibit a continuous
translational symmetry, see Fig. 8. Hence the periodic geodesics are isolated.
This reduction of symmetry is a consequence of the non-uniform curvature of
the surface.

Figure 10:  Periodic geodesics of $(n_{t},1,o^{n_{t}})$ type (red solid lines)
on a 3D Möbius strip with apsect ratio $W/R=0.375$. They resemble the
$(n_{t},1)$ periodic orbits of the flat Möbius billiard which are indicated as
blue dashed lines. (a) Geodesic $(9,1,o^{9})$ with length $L=7.03R$. The
length of the $(9,1)$ PO is $7.13R$. (b) Geodesic $(11,1,o^{11})$ with length
$L=7.44R$. The length of the $(11,1)$ PO is $7.52R$. (c) Geodesic
$(13,1,o^{13})$ with length $L=7.89R$. The length of the $(13,1)$ PO is
$7.95R$.

There is one series of periodic geodesics, the $(n_{t},1,o^{n_{t}})$ orbits
with $n_{t}$ odd, that resembles the $(n_{t},1)$ periodic orbits of the flat
Möbius billiard. Three examples of such periodic geodesics are shown as red
lines in Fig. 10. The analog POs of the flat Möbius billiard are indicated as
blue dashed lines, where the coordinate $x$ ($y$) of the flat Möbius billiard
is identified with $q_{1}$ ($q_{2}$). As the number of reflections $n_{t}$
increases, the segments of the periodic geodesics become less and less curved
in $q$ space and increasingly resemble the straight segments of the POs of the
flat Möbius billiard, where the segments in the region around $q_{1}=0$,
$2\pi$ exhibit generally more curvature than those in the region around
$q_{1}=\pi$. Also the lengths of the periodic geodesics becomes closer to the
lengths of the corresponding POs in the flat Möbius.

These two effects can be understood by considering the metric tensor. For the
flat Möbius billiard, $g_{ij}=\delta_{ij}$, and hence all geodesics are
straight lines. For the 3D Möbius strip, $g_{22}=1$ as well, but $g_{11}$
generally deviates from $1$ (see Fig. 8) and causes the geodesics to be not
straight lines in $q$ space. However, the more the velocity vector $\dot{q}$
of a trajectory becomes parallel to the $q_{2}$ direction, which is
necessarily the case when $n_{t}$ increases, the less important the
contribution of $g_{11}$ is, which results in less curvature. In fact, any
trajectory with $q_{1}=\mathrm{const.}$ is a geodesic of the 3D Möbius strip.
Therefore, the segments of the periodic geodesics resemble more and more
straight lines as $n_{t}$ increases. Moreover, the deviation of $g_{11}$ from
$1$ is smallest in the region around $q_{1}=\pi$, and hence the trajectory
segments are less curved there.

So while the $(n_{t},1,o^{n_{t}})$ periodic geodesics slowly converge towards
the corresponding POs of the flat Möbius billiard, some differences remain.
First, the periodic geodesics remain isolated since the three-dimensional
Möbius lacks the translational symmetry of its flat counterpart. Second, these
periodic geodesics only exist for $n_{t}\geq 9$ for the aspect ratio
$W/R=0.375$ considered here ($n_{t}\geq 11$ for $W/R=0.3$) because the convex
parts of the boundary make some of the longer segments needed for smaller
$n_{t}$ impossible.

While the flat Möbius billiard features POs without reflections at the
boundaries, i.e., the $(0,n_{l})$ orbits in Fig. 7(b), no periodic geodesics
without reflections at the boundary were found for the 3D Möbius strip. In
particular, trajectories with $q_{2}=\mathrm{const.}$ are not geodesics since
$g_{11}\neq 1$, and hence $\dot{q}_{2}\neq 0$ for almost all points of a
geodesic. Therefore, at least one reflection at the boundary is needed to
compensate the change of sign of $\dot{q}_{2}$ after one roundtrip due to the
periodic boundary conditions with inversion of $q_{2}$.

It is also interesting to compare the Möbius strip with conventional
whispering gallery resonators exhibiting azimuthal rotational symmetry such as
rings, spheres or toroids. Such resonators are of interest since they feature
modes with ultra-high quality ($Q$) factors due to the existence of whispering
gallery trajectories that propagate along the outer resonator wall with a
constant angle of incidence above the critical angle for total internal
reflection. High-$Q$ modes exist even for small deformations of the rotational
symmetry that break angular momentum conservation.

Evidently, the 3D Möbius strip does not have an azimuthal symmetry that would
preserve the angular momentum. More importantly, whispering gallery type
trajectories, which would have a symbolic code $s^{n_{t}}$, cannot exist since
parts of the boundary are convex and thus locally prohibit $s$-type segments.
Moreover, the boundary of the Möbius strip is not a geodesic itself. Hence,
all periodic geodesics must feature at least one $o$-segment that crosses the
strip, and in many cases an $o$-segment passing through or near the symmetry
point $Q=(\pi,0)$ is found. As a consequence, the 3D Möbius strip cannot
support whispering gallery modes. Nonetheless, resonances with high quality
factors exist since many periodic geodesics such as the ones shown above are
confined by total internal reflection at all reflections.

##  Appendix D FDTD simulations of the modes of a Möbius strip

In order to better understand the relation between the geodesics and the
resonant modes of the Möbius strip, several mode patterns have been calculated
in addition to the one shown in Fig. 6. The spectrum of a Möbius strip and two
further examples of mode patterns are presented here.

###  D.1 Spectrum of the Möbius strip

To obtain the spectrum of a Möbius strip, a 3D Finite Difference Time Domain
(FDTD) simulation is performed with a temporally short (spectrally broad)
excitation. The time evolution of the field at one point in the Möbius strip
is recorded. To get an accurate estimate of the frequencies and lifetimes from
the recorded time data, a spectral estimation technique known as Prony’s
method is used Marple . The calculated spectrum in Fig. 11 features several
series of resonant modes.

Figure 11:  Spectrum calculated by 3D FDTD simulations for a Möbius strip with
$R=10~{}\mu$m, $W=3~{}\mu$m, and $h=150$ nm. The arrow indicates the wave
function plotted in Fig. 6.

###  D.2 Fundamental transverse mode

To calculate the mode pattern of a particular resonance, a narrow band
excitation at the targeted frequency is used so that only the targeted mode is
excited in the FDTD simulation.

In addition to the pattern depicted in Fig. 6, a low-loss mode at the lower
frequency of $\mathrm{Re}(kR)=66.54$ [$\mathrm{Im}(kR)=-2.46\cdot 10^{-2}$],
but belonging to the same series of resonances, is shown in Fig. 12. This
series consists of the fundamental transverse modes. The corresponding
periodic geodesic (called 4a) is superimposed. This geodesic has only four
vertices and is hence different from the geodesic 5a shown in the Fig. 6 as
can be seen in Fig. 13. The mean position calculated from Eq. (16) for the
wave function in Fig. 12 is also shown for comparison.

Figure 12:  Different views of the same wave function (plot of
$|\vec{E}|^{2}$) with $\mathrm{Re}(kR)=66.54$ and $\mathrm{Im}(kR)=-2.46\cdot
10^{-2}$. The geodesic 4a is superimposed as red line.

###  D.3 Higher order transverse mode

As in other waveguides, there also exist higher order modes with several
antinodes in the transverse ($q_{2}$) direction in the Möbius strip. One such
mode is plotted in Fig. 14. The real and imaginary part of its frequency are
$\mathrm{Re}(kR)=87.73$ and $\mathrm{Im}(kR)=-1.15\cdot 10^{-3}$,
respectively. This mode is almost at the same frequency as the mode presented
in Fig. 6 but has higher losses. In contrast to previous modes which belong to
the resonance series with the lowest losses, this mode belongs to the second
series with larger losses consisting of higher order transverse modes. These
modes cannot be easily related to a known geodesic.

Figure 13:  Calculated periodic geodesic 4a with four vertices (red) and mean
position calculated from the wave function of Fig. 12 (blue). Figure 14:
Different views of the same wave function (plot of $|\vec{E}|^{2}$) with
$\mathrm{Re}(kR)=87.73$ and $\mathrm{Im}(kR)=-1.15\cdot 10^{-3}$.

## References

  * (1) J. Goldstone and R. L. Jaffe, _Bound states in twisting tubes_ , Phys. Rev. B, 45, 14100 (1992).
  * (2) V. H. Schultheiss, S. Batz, A. Szameit, F. Dreisow, S. Nolte, A. Tünnermann, S. Longhi, and U. Peschel, _Optics in curved space_ , Phys. Rev. Lett. , 105, 143901 (2010).
  * (3) M. Avlund, C. Ellegaard, M. Oxborrow, T. Guhr, and N. Sondergaard, _Observation of Periodic Orbits on Curved Two-Dimensional Geometries_ , Phys. Rev. Lett. , 104, 164101 (2010).
  * (4) A. Libster-Hershko, R. Shiloh, and A. Arie, _Surface plasmon polaritons on curved surfaces_ , Optica, 6, 115 (2019).
  * (5) S. Tanda, T. Tsuneta, Y. Okajima, K. Inagaki, K. Yamaya, and N. Hatakenaka, _A Möbius strip of single crystals_ , Nature, 417, 397 (2002).
  * (6) L.-T. Huang and D.-H. Lee, _Topological insulators on a Möbius strip_ , Phys. Rev. B, 84, 193106 (2011).
  * (7) Zehao Li and L. R. Ram-Mohan, _Quantum mechanics on a Mobius ring: Energy levels, symmetry, optical transitions, and level splitting in a magnetic field_ , Phys. Rev. B, 85, 195438 (2012).
  * (8) Y. Yin, S. Li, V. Engemaier, E. S. Ghareh Naz, S. Giudicatti, L. Ma, and O. G. Schmidt, _Topology induced anomalous plasmon modes in metallic Mobius nanorings_ , Laser and Photonics Review, 11, 1600219 (2017).
  * (9) J. Kreismann and M. Hentschel, _The optical Möbius strip cavity: Tailoring geometric phases and far fields_ , Europhysics Lett., 121, 24001 (2018).
  * (10) M. Brack and R. K. Bhaduri, _Semiclassical physics_ , Addison-Wesley publishing company (1997).
  * (11) H.-J. Stöckmann, _Quantum chaos, an introduction_ , Cambridge University Press (1999).
  * (12) Y. Colin de Verdière, _Spectrum of the Laplace operator and periodic geodesics: thirty years after_ , Annales de l’Institut Fourier, Tome 57, 2429 (2007).
  * (13) C. Xu, I. Dana, L.-G. Wang, and P. Sebbah, _Classical and wave chaos in the transformation from curved to flat surfaces_ , Arxiv:2010.12220v1 (2020).
  * (14) Joachim Fischer and Martin Wegener, _Three-dimensional optical laser lithography beyong the diffraction limit_ , Laser and Photonics Review, 7, 22 (2013).
  * (15) M. A. Guidry, Y. Song, C. Lafargue, R. Sobczyk, D. Decanini, S. Bittner, B. Dietz, L. Huang, J. Zyss, A. Grigis, and M. Lebental, _Three-dimensional micro-billiard lasers: The square pyramid_ , Europhysics Lett., 126, 64004 (2019).
  * (16) M. Lebental, N. Djellali, C. Arnaud, J.-S. Lauret, J. Zyss, R. Dubertrand, C. Schmit, E. Bogomolny, Phys. Rev. A, 76, 023830 (2007).
  * (17) R. C. T. da Costa, _Quantum mechanics of a constrained particle_ Phys. Rev. A, 23, 1982 (1981).
  * (18) R. S. Elliott, _Azimuthal surface waves on circular cylinders_ , J. App. Phys., 26, 368 (1955).
  * (19) T. Gissibl, S. Wagner, J. Sykora, M. Schmid, and H. Giessen, _Refractive index measurements of photoresists for three-dimensional direct laser writing_ Optical Materials Express, 7, 2293 (2017).
  * (20) Kimhong Chao, X. Checoury, and M. Lebental, to be submitted.
  * (21) W. Kühnel, _Differential geometry. Curves - surfaces - manifolds_ , Third edition, American Mathematical Society (2015).
  * (22) N. Sobeshchuk, M. A. Guidry, C. Lafargue, R. Gashemi, D. Decanini, J. Zyss, and M. Lebental, _Out-of-plane modes in three-dimensional Fabry-Perot microlasers_ , App. Phys. Lett., 112, 261102 (2018).
  * (23) J. Gravesen and M. Willatzen, _Eigenstates of Möbius nanostructures including curvature effects,_ Phys. Rev. A, 72, 032108 (2005).
  * (24) M. Nakahara, _Geometry, topology and Physics_ , second edition, Taylor and Francis (2003).
  * (25) S. Lawrence Marple, Jr., _Digital Spectral Analysis_ , second edition, Courier Dover Publications (2019).


††thanks: Corresponding author††thanks: Corresponding author

# Bayesian Inference of the Symmetry Energy and the Neutron Skin in 48Ca and
208Pb from CREX and PREX-2

Zhen Zhang<EMAIL_ADDRESS>Sino-French Institute of Nuclear
Engineering and Technology, Sun Yat-sen University, Zhuhai 519082, China  Lie-
Wen Chen<EMAIL_ADDRESS>School of Physics and Astronomy, Shanghai Key
Laboratory for Particle Physics and Cosmology, and Key Laboratory for Particle
Astrophysics and Cosmology (MOE), Shanghai Jiao Tong University, Shanghai
200240, China

###### Abstract

Using the recent model-independent determination of the charge-weak form
factor difference $\Delta F_{\rm CW}$ in [48]Ca and [208]Pb by the CREX and
PREX-2 collaborations, we perform Bayesian inference of the symmetry energy
$E_{\rm sym}(\rho)$ and the neutron skin thickness $\Delta r_{\rm np}$ of 48Ca
and 208Pb within the Skyrme energy density functional (EDF). We find the
inferred $E_{\rm sym}(\rho)$ and $\Delta r_{\rm np}$ separately from CREX and
PREX-2 are compatible with each other at $90\%$ C.L., although they exhibit
strong tension at $68.3\%$ C.L. with CREX (PREX-2) favoring a very soft
(stiff) $E_{\rm sym}(\rho)$ and rather small (large) $\Delta r_{\rm np}$. By
combining the CREX and PREX-2 data, we obtain $E_{\rm
sym}(\rho_{0})=30.2_{-3.0}^{+4.1}$ MeV, the symmetry energy slope parameter
$L=15.3^{+46.8}_{-41.5}$ MeV at saturation density $\rho_{0}$, $\Delta r_{\rm
np}(^{48}{\rm Ca})=0.149_{-0.032}^{+0.035}$ fm and $\Delta r_{\rm
np}(^{208}{\rm Pb})=0.139_{-0.063}^{+0.070}$ fm at $90\%$ C.L.. The
$E_{\rm{sym}}(\rho)$ and $\Delta r_{\rm np}$ from combining CREX and PREX-2
are closer to the corresponding results from CREX alone, implying the PREX-2
is less effective to constrain the $E_{\rm sym}(\rho)$ and $\Delta r_{\rm np}$
due to its lower precision of $\Delta F_{\rm CW}$. Furthermore, we find the
Skyrme EDF results inferred by combining the CREX and PREX-2 data at $90\%$
C.L. nicely agree with the measured dipole polarizabilities $\alpha_{D}$ in
[48]Ca and [208]Pb as well as the neutron matter equation of state from
microscopic calculations.

_Introduction._ — The CREX Adhikari  _et al._ (2022) and PREX-2 Adhikari  _et
al._ (2021) collaborations recently reported precise measurements of the
parity-violating asymmetry $A_{\rm PV}$ in the elastic scattering of
longitudinally polarized electrons from [48]Ca and [208]Pb, respectively.
These measurements lead to model-independent extractions of the difference
between the charge form factor $F_{C}$ and the weak form factor $F_{W}$, i.e.,
$\Delta F_{{\rm{CW}}}(q)\equiv F_{C}(q)-F_{W}(q)=0.0277\pm 0.0055$ at
$q=0.8733~{}\rm{fm}^{-1}$ for [48]Ca and $\Delta F_{\rm{CW}}(q)=0.041\pm
0.013$ at a smaller four-momentum transfer $q=0.3977~{}\rm{fm}^{-1}$ for
[208]Pb Adhikari  _et al._ (2022). Since these extractions are free from the
strong interaction uncertainties, they allow to determine with minimal model-
dependence the neutron skin thickness $\Delta r_{\rm np}\equiv r_{n}-r_{p}$
[$r_{n(p)}$ is the neutron(proton) rms radius of the nucleus] and further to
constrain the density dependence of the symmetry energy $E_{\rm sym}(\rho)$
Brown (2000); Typel and Brown (2001); Horowitz and Piekarewicz (2001);
Furnstahl (2002); Chen  _et al._ (2005); Centelles  _et al._ (2009); Warda
_et al._ (2009); Roca-Maza  _et al._ (2011); Agrawal  _et al._ (2012); Zhang
and Chen (2013). The $E_{\rm sym}(\rho)$ encodes the isospin dependence of
nuclear matter equation of state (EOS) and plays an important role in both
nuclear physics and astrophysics Baran  _et al._ (2005); Steiner  _et al._
(2005); Lattimer and Prakash (2007); Li  _et al._ (2008).

From the PREX-2 data, the $\Delta r_{\rm np}$ of [208]Pb is extracted to be
$0.283\pm 0.071~{}\rm{fm}$ Adhikari  _et al._ (2021). An analysis based on a
relativistic energy density functional (EDF) indicates that such a thick
neutron skin in [208]Pb leads to a very stiff $E_{\rm sym}(\rho)$ with a
rather large symmetry energy slope parameter of $L=106\pm 37$ MeV at
saturation density $\rho_{0}$ Reed  _et al._ (2021), which seriously
challenges our current knowledge on the $E_{\rm sym}(\rho)$ Lattimer and
Steiner (2014); Li and Han (2013); Oertel  _et al._ (2017). Many studies have
been devoted to understanding the PREX-2 result and exploring its implication
in nuclear physics and astrophysics Yue  _et al._ (2022); Reed  _et al._
(2021); Piekarewicz (2021); Reinhard  _et al._ (2021a). In particular, the
tension between the PREX-2 data and the measured electric dipole
polarizabilities $\alpha_{D}$ in [48]Ca and [208]Pb at RCNP in Osaka Tamii
_et al._ (2011); Roca-Maza  _et al._ (2015); Birkhan  _et al._ (2017) is
observed with the latter favoring a much softer $E_{\rm sym}(\rho)$
Piekarewicz (2021); Reinhard  _et al._ (2021a).

Very remarkably, the CREX collaboration adopt the same experimental approach
as PREX-2 and recently report a thin neutron skin of $\Delta r_{\rm
np}=0.121\pm 0.026(\rm exp)\pm 0.024(\rm model)$ fm in [48]Ca Adhikari  _et
al._ (2022). Analyses within modern nonrelativistic and relativistic EDFs
Reinhard  _et al._ (2022); Yüksel and Paar (2022) conclude a significant
tension between the CREX and PREX-2 results, calling for further critical
theoretical and experimental investigations.

In this Letter, we employ the Bayesian inference method, which provides a
consistent probabilistic approach to extract quantitative information from
experimental measurements von Toussaint (2011), to analyze the CREX and PREX
results on the $\Delta F_{{\rm{CW}}}(q)$ together with other well-known
ground-sate and collective-excitation data of [48]Ca and [208]Pb based on the
Skyrme EDF. Our Bayesian inference indicates the CREX and PREX-2 results are
compatible at $90\%$ confidence level (C.L.), although they exhibit tension at
$68.3\%$ C.L., and furthermore the PREX-2 is less effective to constrain the
$E_{\rm sym}(\rho)$ and $\Delta r_{\rm np}$ due to its lower precision of
$\Delta F_{\rm CW}$ compared to the CREX. By combining the CREX and PREX-2
results at $90\%$C.L., we find the Skyrme EDF can nicely describe the measured
$\alpha_{D}$ in [48]Ca and [208]Pb as well as the neutron matter EOS from
microscopic many-body calculations.

_Model and method._ — The nuclear properties are calculated within the widely
used standard Skyrme EDF. Since we focus on the doubly magic [48]Ca and 208Pb,
pairing interaction is not taken into account. The Skyrme EDF can then be
characterized by ten parameters: the nuclear saturation density $\rho_{0}$,
the binding energy per nucleon of symmetric nuclear matter $E_{0}(\rho_{0})$,
the incompressibility $K_{0}$, $E_{\mathrm{sym}}(\rho_{0})$, $L$, the
isoscalar effective mass $m_{s,0}^{*}$ and the isovector effective mass
$m_{v,0}^{*}$ at $\rho_{0}$, the gradient coefficient $G_{S}$, the symmetry-
gradient coefficient $G_{V}$, and the spin-orbit coupling constant $W_{0}$
Chen  _et al._ (2010); Chen and Gu (2012); Kortelainen  _et al._ (2010). Based
on the Skyrme EDF, once given a parameter set

$\displaystyle\bm{p}$ $\displaystyle=$
$\displaystyle\\{\rho_{0},~{}E_{0}(\rho_{0}),~{}K_{0},~{}E_{\mathrm{sym}}(\rho_{0}),~{}L,$
(1)
$\displaystyle~{}G_{S},~{}G_{V},~{}W_{0},~{}m_{s,0}^{\ast},~{}m_{v,0}^{\ast}\\},$

the ground-state properties of [208]Pb and [48]Ca are calculated with the
Hartree-Fock (HF) method, and the breading mode energy is obtained from the
constrained HF (CHF) calculation.

Table 1: Prior ranges of the ten parameters used, together with the posterior median values and $68.3\%$($90\%$) credible intervals from B-All and the parameter values of the Skyrme interaction SkREx. Quantity | prior | posterior (B-All) | SkREx
---|---|---|---
$\rho_{0}~{}(\mathrm{fm}^{-3})$ | $[0.155,0.165]$ | $0.1611_{-0.0017(0.0028)}^{+0.0017(0.0027)}$ | 0.1590
$E_{0}~{}(\mathrm{MeV})$ | $[-16.5,-15.5]$ | $-16.05_{-0.06(0.10)}^{+0.06(0.10)}$ | -16.03
$K_{0}~{}(\mathrm{MeV})$ | $[210,250]$ | $225.4_{-3.0(4.9)}^{+3.1(5.2)}$ | 233.1
$E_{\mathrm{sym}}(\rho_{0})~{}(\mathrm{MeV})$ | $[22,55]$ | $30.2_{-1.9(3.0)}^{+2.3(4.1)}$ | 29.0
$L~{}(\mathrm{MeV})$ | $[-90,240]$ | $15.3_{-25.8(41.5)}^{+27.5(46.8)}$ | 12.8
$G_{S}~{}(\mathrm{MeV}\cdot\mathrm{fm}^{5})$ | $[110,170]$ | $126.6_{-7.0(10.8)}^{+9.4(17.4)}$ | 115.2
$G_{V}~{}(\mathrm{MeV}\cdot\mathrm{fm}^{5})$ | $[-70,70]$ | $-20.1_{-35.3(45.3)}^{+50.5(74.5)}$ | -46.6
$W_{0}~{}(\mathrm{MeV}\cdot\mathrm{fm}^{5})$ | $[90,140]$ | $118.6_{-15.9(23.6)}^{+13.9(18.9)}$ | 95.3
$m_{s,0}^{\ast}/m$ | $[0.7,1.0]$ | $0.928_{-0.085(0.150)}^{+0.052(0.066)}$ | 0.979
$m_{v,0}^{\ast}/m$ | $[0.6,0.9]$ | $0.726_{-0.091(0.115)}^{+0.113(0.155)}$ | 0.601

The calibration and uncertainty quantification of the ten parameters is
carried out using a Bayesian approach. According to Bayes’ theorem, the
posterior distribution of model parameters $\bm{p}$, given experimental data
$\mathcal{O}^{\rm{exp}}$ for a set of observables $\mathcal{O}$, can be
evaluated as

$\displaystyle
P\left(\bm{p}\mid\mathcal{M},\mathcal{O}^{\exp}\right)=\frac{P\left(\mathcal{O}^{\exp}\mid\mathcal{M},\bm{p}\right)P(\bm{p})}{\int
P\left(\mathcal{O}^{\exp}\mid\mathcal{M},\bm{p}\right)P(\bm{p})\mathrm{d}\bm{p}},$
(2)

where $\mathcal{M}$ is the given model, $P(\bm{p})$ is the prior probability
density of model parameters $\bm{p}$ before being confronted with the data
$\mathcal{O}^{\exp}$, and $P\left(O^{\exp}\mid\mathcal{M},\bm{p}\right)$
denotes the likelihood or the conditional probability distribution of
observing $\mathcal{O}^{\exp}$ with given model $\mathcal{M}$ predictions at
$\bm{p}$. The prior distribution of model parameters are normally chosen to be
uniform in their empirical ranges listed in Tab. 1. In particular, given the
rather thick (thin) neutron skin in [208]Pb ([48]Ca) from PREX-2 (CREX), the
prior ranges of $E_{\mathrm{sym}}(\rho_{0})$ and $L$ are taken to be as large
as $22\sim 55$ MeV and $-90\sim 250$ MeV, respectively.

The likelihood function is taken to be the commonly used Gaussian form

$\displaystyle
P\left(\mathcal{M},O_{i}^{\exp}\mid\bm{p}\right)\propto\exp\left\\{-\sum_{i}\frac{\left[\mathcal{O}_{i}(\bm{p})-\mathcal{O}_{i}^{\exp}\right]^{2}}{2\sigma_{i}^{2}}\right\\},$
(3)

where $\mathcal{O}_{i}(\bm{p})$ is the model prediction on the $i$-th
observable for a given parameter set $\bm{p},\mathcal{O}_{i}^{\exp}$ is the
corresponding data, and $\sigma_{i}$ is the adopted error.

To estimate the posterior distribution given by Eq. (2), the Markov chain
Monte Carlo (MCMC) process is carried out using the Methropolis-Hasting
algorithm. We first run $5\times 10^{5}$ burn-in MCMC steps to allow the chain
to reach equilibrium, and then generate $10^{6}$ MCMC steps in parameter
space. The posterior distributions of model parameters and observables are
estimated from 15 parallel MCMC process, i.e., $1.5\times 10^{7}$ MCMC
samples.

The normalized nuclear form factors $F_{C}(q)$ and $F_{W}(q)$ are calculated
by folding the nucleon form factor $F_{t}(q)~{}(t=n,p)$ and the spin-orbit
current form factor ($F_{t}^{ls}$) with the intrinsic nucleon electromagnetic
form factor $G_{E/M,t}$ and weak form factor $G^{(W)}_{E/M,t}$ by Reinhard
_et al._ (2021a)

$\displaystyle F_{C}(q)$ $\displaystyle=$
$\displaystyle\frac{1}{Z}\sum_{t=p,n}\left[G_{E,t}(q)F_{t}(q)+G_{M,t}(q)F_{t}^{(ls)}(q)\right],$
(4) $\displaystyle F_{W}(q)$ $\displaystyle=$
$\displaystyle\sum_{t=p,n}\frac{\left[G_{E,t}^{(W)}(q)F_{t}(q)+G_{M,t}^{(W)}(q)F_{t}^{(ls)}(q)\right]}{ZQ_{p}^{(W)}+NQ_{n}^{(W)}},$
(5)

where $N(Z)$ is the neutron(proton) number, and $Q_{p}^{(W)}=0.0713$ and
$Q_{n}^{(W)}=-0.9888$ are proton and neutron weak charges, respectively. The
$G_{E/M,t}$ are derived from the isospin-coupled Sachs form factors, the
relativistic Darwin correction has been included and the center-of-mass
corrections are taken into account by simply renormalizing the nucleon mass
$m_{N}$ to $(1-1/A)m_{N}$ in the HF calculation (see, e.g., Ref. Bender  _et
al._ (2003) for details). From $G_{E/M,t}$, the $G^{(W)}_{E/M,t}$ can then be
determined by further considering the contribution of the strange-quark
electromagnetic form factors $G_{E/M,s}$ (see, Ref. Reinhard  _et al._ (2021a)
for details).

_Selected observables._ — The key observables in this work are the model-
independent $\Delta F_{\rm{CW}}$ in [48]Ca and [208]Pb, i.e., $\Delta
F_{\rm{CW}}^{48}\equiv\Delta F_{\rm{CW}}(q=0.8733~{}\rm{fm}^{-1})$ from CREX
Adhikari  _et al._ (2022) and $\Delta F_{\rm{CW}}^{208}\equiv\Delta
F_{\rm{CW}}(q=0.3977~{}\rm{fm}^{-1})$ from PREX-2 Adhikari  _et al._ (2021).
Besides, we also include three quantities defined according to the well-
determined $F_{C}$ Friedrich and Voegler (1982); Klupfel  _et al._ (2009),
i.e., the charge rms radius

$r_{C}=\sqrt{-\left.\frac{3}{F_{C}(0)}\frac{d^{2}}{dq^{2}}F_{C}(q)\right|_{q=0}},$
(6)

the diffraction radius $R_{d}=\frac{4.493}{q_{0}}$ determined from the first
zero of $F_{C}[q_{0}]=0$, and the surface thickness

$\sigma=\sqrt{\frac{2}{q_{m}}\log\left[\frac{3j_{1}(q_{m}R_{d})}{q_{m}R_{d}F_{C}\left(q_{m}\right)}\right]},\quad
q_{m}=5.6/R,$ (7)

where $j_{1}(x)=\frac{\sin x}{x^{2}}-\frac{\cos x}{x}$ is the spherical Bessel
function of the first kind.

To guarantee the MCMC process is carried out in reasonable parameter space, we
further include in our analysis the total binding energies $E_{B}$ and the
neutron-proton Fermi energy difference $\Delta\epsilon_{F}$ in both [48]Ca and
[208]Pb, along with the neutron $3p_{1/2}-3p_{1/3}$ splitting
$\Delta\epsilon_{ls}$ and the breathing mode energy $E_{\rm{GMR}}$ in [208]Pb.
We note that the $\Delta\epsilon_{F}$ provides a sensitive probe of the
$E_{\rm sym}(\rho)$ at $\rho\approx 2/3\rho_{0}$ Wang  _et al._ (2013).

Table 2: The data and adopted errors $\sigma_{i}$ of $E_{B}$ Wang  _et al._ (2017), $r_{C}$ Angeli and Marinova (2013), $R_{d}$ Klupfel  _et al._ (2009), $\sigma$ Klupfel  _et al._ (2009), $\Delta F_{\rm{CW}}$ Adhikari  _et al._ (2022, 2021) and $\Delta\epsilon_{F}$ Wang  _et al._ (2013); Schwierz  _et al._ (2007) for [208]Pb and [48]Ca. The superscript “48”(“208”) represents [48]Ca([208]Pb). Also listed are the $E_{\mathrm{GMR}}$ measured by RCNP Patel  _et al._ (2013); Howard  _et al._ (2020) and TAMU Youngblood  _et al._ (1999) groups, and the neutron $3p_{1/2}$-$3p_{2/3}$ splitting $\Delta\epsilon_{ls}$ Klupfel  _et al._ (2009) in [208]Pb.  | value | $\sigma_{i}$ |  | value | $\sigma_{i}$
---|---|---|---|---|---
$E_{B}^{208}$ (MeV) | -1636.43 | 1 |  $E_{B}^{48}$ (MeV) | -416.001 | 1
$r_{C}^{208}$ (fm) | 5.5012 | 0.02 |  $r_{C}^{48}$ (fm) | 3.479 | 0.04
$R_{d}^{208}$ (fm) | 6.776 | 0.04 |  $R_{d}^{48}$ (fm) | 3.964 | 0.04
$\sigma^{208}$ (fm) | 0.913 | 0.04 |  $\sigma^{48}$ (fm) | 0.881 | 0.04
$\Delta\epsilon_{F}^{208}$ (MeV) | 0.64 | 0.13 |  $\Delta\epsilon_{F}^{48}$ (MeV) | 6.10 | 1.22
$\Delta F_{\rm{CW}}^{208}$ | 0.041 | 0.013 | $\Delta F_{\rm{CW}}^{48}$ | 0.0277 | 0.0055
$E_{\mathrm{GMR}}^{\rm{(RCNP)}}$ (MeV) | 13.1 | 0.1 | $E_{\mathrm{GMR}}^{\rm{(TAMU)}}$ | 14.18 | 0.11
$\Delta\epsilon_{ls}$ (MeV) | 0.90 | 0.18 |  |  |

Table 2 summarizes the values and the adopted errors of the selected
observables. Since the $E_{B}$, $r_{C}$ and $\Delta\epsilon_{F}$ in [208]Pb
and [48]Ca as well as the neutron $3p_{1/2}$-$3p_{3/2}$ splitting in [208]Pb
are very accurately measured, the adopted errors are taken to be the usual
values in the fitting of Skyrme EDF Klupfel  _et al._ (2009). As for the
$E_{\rm{GMR}}$ of [208]Pb, the two independently measured values by RCNP Patel
_et al._ (2013); Howard  _et al._ (2020) and TAMU Youngblood  _et al._ (1999)
are both used.

_Results and discussions._ — Using the data listed in Tab. 2, we conduct
Bayesian analyses in four different cases. The Bayesian analysis without
including the $\Delta F_{\rm{CW}}^{208}$ and $\Delta F_{\rm{CW}}^{48}$ is
considered as the base point labeled with “B-Bas”. The $\Delta
F_{\rm{CW}}^{48}$ and $\Delta F_{\rm{CW}}^{208}$ data are then separately
added into the analysis to quantify the tension between the CREX and PREX
data, and the results are accordingly labeled by “B-$\Delta F_{\rm{CW}}^{48}$”
and “B-$\Delta F_{\rm{CW}}^{208}$”. A Bayesian analysis including all data
listed in Tab. 2 labeled by “B-All” is further carried out to constrain the
$E_{\rm sym}(\rho)$ by combining the CREX and PREX results. The posterior
median values and $68.3\%(90\%)$ credible intervals of parameters $\bm{p}$
obtained with B-All are listed in Tab. 1.

Figure 1 shows the obtained posterior joint distributions between $\Delta
F_{\rm{CW}}^{48}$ and $\Delta F_{\rm{CW}}^{208}$ with $68.3\%$ and $90\%$
credible regions from B-$\Delta F_{\rm{CW}}^{48}$, B-$\Delta
F_{\rm{CW}}^{208}$ and B-All, together with the corresponding experimental
joint distribution of $90\%$ credible region as well as the CREX and PREX
individual measurements with $90\%$ uncertainties. One clearly sees the
tension between the $\Delta F_{\rm{CW}}^{48}$ and $\Delta F_{\rm{CW}}^{208}$
data at $1\sigma$ ($68.3\%$) C.L., with the corresponding credible regions
having no overlap, confirming the conclusion in Refs. Reinhard  _et al._
(2022); Yüksel and Paar (2022). However, at $90\%$ C.L., the obtained credible
regions are compatible with each other as well as the experimental one. By
combining the $\Delta F_{\rm{CW}}^{48}$ and $\Delta F_{\rm{CW}}^{208}$ data,
the B-All prediction is between the results with B-$\Delta F_{\rm{CW}}^{48}$
and B-$\Delta F_{\rm{CW}}^{208}$, and has larger overlap with the experimental
region. From the MCMC samples, we find a Skyrme interaction (named “SkREx”)
that is consistent with both the CREX and PREX data at $90\%$ C.L. (as
indicated by star in Fig. 1), the measured $\alpha_{D}$ in [48]Ca and [208]Pb,
and the neutron matter EOS from microscopic calculations as shown later. The
parameter values of SkREx are listed in the last column of Tab. 1.

Figure 1: (Color online). Joint distributions of $\Delta
F_{\rm{CW}}^{48}$-$\Delta F_{\rm{CW}}^{208}$ from B-$\Delta F_{\rm{CW}}^{48}$
(blue, dash-dotted line), B-$\Delta F_{\rm{CW}}^{208}$ (green, dashed line)
and B-All (red, solid line). The shaded regions and the lines correspond to
$68.3\%$ and $90\%$ credible regions, respectively. The $90\%$ credible
regions of the experimental joint and marginal distributions by CREX and
PREX-2 are indicated as the dotted ellipse and rectangle, respectively. Figure
2: (Color online). Posterior distributions of $E_{\mathrm{sym}}(2\rho_{0}/3)$
(a), $L(2\rho_{0}/3)$ (b), $E_{\mathrm{sym}}(\rho_{0})$ (c) and $L$ (d) in the
four cases (see text for details).

Shown in Fig. 2 are the posterior distributions of
$E_{\mathrm{sym}}(2\rho_{0}/3)$, $L(2\rho_{0}/3)$,
$E_{\mathrm{sym}}(\rho_{0})$ and $L$ from the Bayesian analysis in the four
cases. It is seen from Fig. 2 (a) that for B-Bas, which includes the
$\Delta\epsilon_{F}$ of [208]Pb and [48]Ca in the analysis, the
$E_{\mathrm{sym}}(2\rho_{0}/3)$ has already been well constrained, and further
including the data of $\Delta F_{\rm{CW}}$ has minor effects on the posterior
distribution of $E_{\mathrm{sym}}(2\rho_{0}/3)$. This observation confirms the
sensitivity of $\Delta\epsilon_{F}$ to the $E_{\rm{sym}}(2\rho_{0}/3)$ Wang
_et al._ (2013). With the $E_{\rm{sym}}(2\rho_{0}/3)$ tightly constrained by
the $\Delta\epsilon_{F}^{208}$ and $\Delta\epsilon_{F}^{48}$, the
$L(2\rho_{0}/3)$, $E_{\rm{sym}}(\rho_{0})$ and $L$ become highly correlated,
leading to very similar shapes for their posterior distributions as shown in
Figs. 2 (b), (c) and (d). Unlike the $E_{\rm{sym}}(2\rho_{0}/3)$, the
$L(2\rho_{0}/3)$, $E_{\mathrm{sym}}(\rho_{0})$ and $L$ are all weakly
constrained with B-Bas, with their posterior distributions exhibiting large
distribution widths. Comparing the results from B-$\Delta F_{\rm{CW}}^{48}$
and B-$\Delta F_{\rm{CW}}^{208}$, one sees the tension between CREX and
PREX-2. For example, the measured $\Delta F_{\rm{CW}}^{208}$ by PREX-2
collaboration favors a larger $L$ [i.e. $111.7_{-46.2(75.6)}^{+48.2(76.6)}$
MeV at $68.3\%$($90\%$) C.L.], while the $\Delta F_{\rm{CW}}^{48}$ by CREX
prefers a much smaller $L$ [i.e., $-10.0_{-26.7(42.9)}^{+29.2(49.1)}$ at
$68.3\%$($90\%$) C.L.]. The $68.3\%$ credible intervals obtained from the CREX
and PREX data are incompatible, whereas the $90\%$ credible intervals overlap
within the region of $L=36.1\sim 59.1$ MeV. The overlap region of the $L$
distributions extracted from CREX and PREX-2 amounts for only about $11\%$.

From the posterior distributions of $E_{\mathrm{sym}}(2\rho_{0}/3)$,
$L(2\rho_{0}/3)$, $E_{\mathrm{sym}}(\rho_{0})$ and $L$ from B-All in Fig. 2,
we obtain $E_{\mathrm{sym}}(2\rho_{0}/3)=26.7_{-0.5(0.8)}^{+0.6(1.1)}$ MeV,
$L(2\rho_{0}/3)=34.1_{-10.6(17.0)}^{+11.7(20.0)}$ MeV,
$E_{\mathrm{sym}}(\rho_{0})=30.2_{-1.9(3.0)}^{+2.3(4.1)}$ MeV and
$L=15.3_{-25.8(41.5)}^{+27.6(46.8)}$ MeV at $68.3\%$($90\%$) C.L.. The
obtained $E_{\rm{sym}}(2\rho_{0}/3)=26.7_{-0.8}^{+1.1}$ MeV at $90\%$ C.L. is
consistent with $E_{\rm{sym}}(2\rho_{0}/3)\approx 26~{}\mathrm{MeV}$ Horowitz
and Piekarewicz (2001) and $E_{\mathrm{sym}}(0.1~{}\mathrm{fm}^{-3})=25.4\pm
0.8~{}\mathrm{MeV}$ Brown (2013) obtained respectively from relativistic and
nonrelativistic EDFs constrained by nuclear masses, as well as
$E_{\mathrm{sym}}(0.11~{}\mathrm{fm}^{-3})=26.2\pm 1.0$ MeV extracted from
$\Delta\epsilon_{F}$ in doubly magic nuclei Wang  _et al._ (2013). In
particular, the present $E_{\mathrm{sym}}(2\rho_{0}/3)$ is in very good
agreement with the value extracted from the binding energy difference of heavy
isotope pairs Zhang and Chen (2013), i.e.,
$E_{\mathrm{sym}}(\rho_{c})=26.65\pm 0.20$ MeV at
$\rho_{c}=0.11~{}\mathrm{fm}^{-3}\approx 2\rho_{0}/3$ at $95\%$ C.L.. The
obtained $E_{\mathrm{sym}}(\rho_{0})$ and $L$ from B-All are also consistent
with many previous works within the $68.3\%$ uncertainties. For example,
although the obtained median value of $L=15.3$ MeV indicates a very soft
$E_{\rm{sym}}(\rho)$, the upper limit of $L=43$ MeV at $68.3\%$ C.L. still
lies in the constraint of $L=53_{-15}^{+14}$ MeV extracted by combining
astrophysical data, PREX-2 result and chiral effective theory constraints
Essick  _et al._ (2021), and the constraint of $L=54\pm 8$ MeV from EDF
analysis of the parity-violating asymmetry and the $\alpha_{D}$ in [208]Pb
Reinhard  _et al._ (2021a).

Figure 3: (Color online). Posterior distributions of $\Delta r_{\rm{np}}^{48}$
(a) and $\Delta r_{\rm{np}}^{208}$ (b) for B-$\Delta F_{\rm{CW}}^{48}$,
B-$\Delta F_{\rm{CW}}^{208}$ and B-All. Dots and bars indicate the median
values, along with the $68.3\%$ and $90\%$. uncertainties

Shown in Fig. 3 are the posterior distributions of $\Delta r_{\rm{np}}^{48}$
and $\Delta r_{\rm{np}}^{208}$ for B-$\Delta F_{\rm{CW}}^{48}$, B-$\Delta
F_{\rm{CW}}^{208}$ and B-All. As expected, the CREX(PREX-2) data result in
thinner(thicker) $\Delta r_{\rm{np}}$. Again, their $90\%$ credible intervals
overlap with each other. The extracted $\Delta
r_{\rm{np}}^{208}=0.281_{-0.069}^{+0.069}~{}\mathrm{fm}$ agrees well with
$0.283\pm 0.071~{}\rm{fm}$ reported by PREX-2 Adhikari  _et al._ (2021), and
the extracted $\Delta r_{\rm{np}}^{48}=0.130_{-0.021}^{+0.023}~{}\mathrm{fm}$
is also consistent with $0.121\pm 0.026(\rm{exp})\pm
0.024(\rm{model})~{}\rm{fm}$ by CREX Adhikari  _et al._ (2022). Combining the
CREX and PREX data results in $\Delta
r_{\rm{np}}^{208}=0.139_{-0.039(0.063)}^{+0.041(0.070)}$ fm and $\Delta
r_{\rm{np}}^{48}=0.149_{-0.020(0.032)}^{+0.021(0.035)}$ fm at $68.3\%$
($90\%$) C.L.. The predicted $\Delta r_{\rm{np}}$ is relatively thin, but is
still consistent with many previous experimental and theoretical studies Zhang
and Chen (2013); Roca-Maza  _et al._ (2013, 2015), e.g., the very recent ab
initio predictions of $\Delta r_{\rm{np}}^{208}=0.14-0.20~{}\mathrm{fm}$ and
$\Delta r_{\rm{np}}^{48}=0.14-0.19~{}\mathrm{fm}$ at $68.3\%$ C.L. Hu  _et
al._ (2021), and the $\Delta r_{\rm{np}}^{48}=0.15-0.21$ fm from [54]Ni-[54]Fe
charge radius difference Pineda  _et al._ (2021). Overall, our Bayesian
analyses indicate that the CREX and PREX data are compatible with each other
at $90\%$ C.L., although they exhibit tension at $68.3\%$ C.L.. Furthermore,
the inferred results of $E_{\mathrm{sym}}(\rho)$ and $\Delta r_{\rm{np}}$ by
combining the CREX and PREX data much favor the results from CREX alone,
implying the PREX-2 is less effective to constrain the $E_{\rm sym}(\rho)$ and
$\Delta r_{\rm np}$ due to its lower precision of $\Delta F_{\rm CW}$ compared
to the CREX, consistent with the conclusion in Refs. Reinhard  _et al._
(2021a, 2022).

It is instructive to see the Bayesain inference on the neutron matter EOS
$E_{\mathrm{PNM}}(\rho)$ which has been well constrained by microscopic
calculations. Fig. 4 shows the inferred $E_{\mathrm{PNM}}(\rho)$ by combining
the CREX and PREX data at $68.3\%$ and $90\%$ C.L., together with the
predictions from many-body perturbation theory using N3LO chiral interactions
by Tews et al. Tews  _et al._ (2013), Wellenhofer et al. Wellenhofer  _et al._
(2015) and Drishcler et al. Drischler  _et al._ (2020), the quantum Monte
Carlo methods by Gandolfi et al. Gandolfi  _et al._ (2012), Wlazłowski et al.
Wlazłowski  _et al._ (2014) and Roggero et al. Roggero  _et al._ (2014), the
variational calculations by Akmal-Pandharipande-Ravenhall (APR) Akmal  _et
al._ (1998), the Bethe-Bruckner-Goldstone calculations (BBG-QM 3h-gap and BBG-
QM 3h-con) Baldo and Fukukawa (2014), and the self-consistent Green’s function
approach (SCGF-N3LO+N2LOdd) Carbone  _et al._ (2014). One sees that the
inferred $E_{\mathrm{PNM}}(\rho)$ is in good agreement with the microscopic
calculations.

Figure 4: (Color online). The Bayesian inferred $E_{\mathrm{PNM}}(\rho)$ by
combining CREX and PREX-2 data. The results from microscopic calculations and
the Skyrme interaction SkREx are also included for comparison (see text for
detals).

Also included in Fig. 4 is the prediction from SkREx, which well agrees with
the microscopic calculations. We note the SkREx can reproduce the data of the
binding energy better than $0.05\%$, the charge rms radius by $0.25\%$ for
[208]Pb and $0.91\%$ for [48]Ca. In addition, the $\alpha_{{D}}$ is known as
an important isovector indicator Reinhard and Nazarewicz (2010); Piekarewicz
_et al._ (2012); Roca-Maza  _et al._ (2013, 2015); Zhang and Chen (2015); Xu
_et al._ (2020). Analyses based on modern nuclear EDFs suggest a strong
correlation between $\alpha_{D}^{208}E_{\mathrm{sym}}(\rho_{0})$ and $L$ Roca-
Maza  _et al._ (2013, 2015), and further analysis indicates the
$\alpha_{D}^{208}$ is strongly correlated with $E_{\rm{sym}}(\rho_{0}/3)$
Zhang and Chen (2015); Xu  _et al._ (2020). Experimentally, the $\alpha_{D}$
in [208]Pb and [48]Ca have been determined to be $19.6\pm 0.6~{}\rm{fm^{3}}$
Tamii  _et al._ (2011); Roca-Maza  _et al._ (2015) and $2.07\pm
0.22~{}\rm{fm^{3}}$ Birkhan  _et al._ (2017) via forward-angle proton elastic
scattering experiments. From CHF calculation Reinhard  _et al._ (2021b) with
SkREx, we obtain $\alpha_{D}^{48}=2.24~{}\rm{fm^{3}}$ and
$\alpha_{D}^{208}=19.5~{}\rm{fm^{3}}$, agreeing well with the data.

_Conclusions._ — Based on the Bayesian inference method and the Skyrme EDF, we
have demonstrated that the CREX and PREX-2 data are compatible with each other
at $90\%$ C.L., although they exhibit strong tension at $68.3\%$ C.L.. We have
further obtained a new Skyrme interaction SkREx, which can describe the CREX
and PREX-2 data at $90\%$ C.L., the measured $\alpha_{D}$ in [48]Ca and
[208]Pb, and the microscopic neutron matter EOS. Our Bayesian analyses
indicate that the PREX-2 is less effective to constrain the $E_{\rm
sym}(\rho)$ and $\Delta r_{\rm np}$ due to its lower precision of $\Delta
F_{\rm CW}$ compared to the CREX, implying the future MESA experiment Becker
_et al._ (2018) is of particular importance to address the issue.

_Acknowledgements._ — This work was supported in part by the National Natural
Science Foundation of China under Grants No. 11905302 and No. 11625521, and
National SKA Program of China No. 2020SKA0120300.

## References

  * Adhikari  _et al._ (2022) D. Adhikari _et al._ (CREX), (2022), arXiv:2205.11593 [nucl-ex] .
  * Adhikari  _et al._ (2021) D. Adhikari _et al._ (PREX), Phys. Rev. Lett. 126, 172502 (2021), arXiv:2102.10767 [nucl-ex] .
  * Brown (2000) B. A. Brown, Phys. Rev. Lett. 85, 5296 (2000).
  * Typel and Brown (2001) S. Typel and B. A. Brown, Phys. Rev. C 64, 027302 (2001).
  * Horowitz and Piekarewicz (2001) C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001), arXiv:astro-ph/0010227 .
  * Furnstahl (2002) R. J. Furnstahl, Nucl. Phys. A 706, 85 (2002), arXiv:nucl-th/0112085 .
  * Chen  _et al._ (2005) L.-W. Chen, C. M. Ko, and B.-A. Li, Phys. Rev. C 72, 064309 (2005), arXiv:nucl-th/0509009 .
  * Centelles  _et al._ (2009) M. Centelles, X. Roca-Maza, X. Vinas, and M. Warda, Phys. Rev. Lett. 102, 122502 (2009), arXiv:0806.2886 [nucl-th] .
  * Warda  _et al._ (2009) M. Warda, X. Vinas, X. Roca-Maza, and M. Centelles, Phys. Rev. C 80, 024316 (2009), arXiv:0906.0932 [nucl-th] .
  * Roca-Maza  _et al._ (2011) X. Roca-Maza, M. Centelles, X. Vinas, and M. Warda, Phys. Rev. Lett. 106, 252501 (2011), arXiv:1103.1762 [nucl-th] .
  * Agrawal  _et al._ (2012) B. K. Agrawal, J. N. De, and S. K. Samaddar, Phys. Rev. Lett. 109, 262501 (2012), arXiv:1212.0292 [nucl-th] .
  * Zhang and Chen (2013) Z. Zhang and L.-W. Chen, Phys. Lett. B 726, 234 (2013), arXiv:1302.5327 [nucl-th] .
  * Baran  _et al._ (2005) V. Baran, M. Colonna, V. Greco, and M. Di Toro, Phys. Rep. 410, 335 (2005).
  * Steiner  _et al._ (2005) A. W. Steiner, M. Prakash, J. M. Lattimer, and P. J. Ellis, Phys. Rep. 411, 325 (2005).
  * Lattimer and Prakash (2007) J. M. Lattimer and M. Prakash, Phys. Rep. 442, 109 (2007).
  * Li  _et al._ (2008) B. A. Li, L. W. Chen, and C. M. Ko, Phys. Rep. 464, 113 (2008).
  * Reed  _et al._ (2021) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021), arXiv:2101.03193 [nucl-th] .
  * Lattimer and Steiner (2014) J. M. Lattimer and A. W. Steiner, Eur. Phys. J. A 50, 40 (2014).
  * Li and Han (2013) B. A. Li and X. Han, Phys. Lett. B 727, 276 (2013).
  * Oertel  _et al._ (2017) M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017), arXiv:1610.03361 [astro-ph.HE] .
  * Yue  _et al._ (2022) T.-G. Yue, L.-W. Chen, Z. Zhang, and Y. Zhou, Phys. Rev. Res. 4, L022054 (2022), arXiv:2102.05267 [nucl-th] .
  * Piekarewicz (2021) J. Piekarewicz, Phys. Rev. C 104, 024329 (2021), arXiv:2105.13452 [nucl-th] .
  * Reinhard  _et al._ (2021a) P.-G. Reinhard, X. Roca-Maza, and W. Nazarewicz, Phys. Rev. Lett. 127, 232501 (2021a), arXiv:2105.15050 [nucl-th] .
  * Tamii  _et al._ (2011) A. Tamii, I. Poltoratska, P. Von Neumann-Cosel, Y. Fujita, T. Adachi, C. A. Bertulani, J. Carter, M. Dozono, H. Fujita, K. Fujita, K. Hatanaka, D. Ishikawa, M. Itoh, T. Kawabata, Y. Kalmykov, A. M. Krumbholz, E. Litvinova, H. Matsubara, K. Nakanishi, R. Neveling, H. Okamura, H. J. Ong, B. Özel-Tashenov, V. Y. Ponomarev, A. Richter, B. Rubio, H. Sakaguchi, Y. Sakemi, Y. Sasamoto, Y. Shimbara, Y. Shimizu, F. D. Smit, T. Suzuki, Y. Tameshige, J. Wambach, R. Yamada, M. Yosoi, and J. Zenihiro, Phys. Rev. Lett. 107, 062502 (2011).
  * Roca-Maza  _et al._ (2015) X. Roca-Maza, X. Viñas, M. Centelles, B. K. Agrawal, G. Colo’, N. Paar, J. Piekarewicz, and D. Vretenar, Phys. Rev. C 92, 064304 (2015), arXiv:1510.01874 [nucl-th] .
  * Birkhan  _et al._ (2017) J. Birkhan _et al._ , Phys. Rev. Lett. 118, 252501 (2017), arXiv:1611.07072 [nucl-ex] .
  * Reinhard  _et al._ (2022) P.-G. Reinhard, X. Roca-Maza, and W. Nazarewicz, (2022), arXiv:2206.03134 [nucl-th] .
  * Yüksel and Paar (2022) E. Yüksel and N. Paar, (2022), arXiv:2206.06527 [nucl-th] .
  * von Toussaint (2011) U. von Toussaint, Rev. Mod. Phys. 83, 943 (2011).
  * Chen  _et al._ (2010) L. W. Chen, C. M. Ko, B. A. Li, and J. Xu, Phys. Rev. C 82, 024321 (2010).
  * Chen and Gu (2012) L.-W. Chen and J.-Z. Gu, J. Phys. G 39, 035104 (2012).
  * Kortelainen  _et al._ (2010) M. Kortelainen, T. Lesinski, J. Moré, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild, Phys. Rev. C 82, 024313 (2010).
  * Bender  _et al._ (2003) M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * Friedrich and Voegler (1982) J. Friedrich and N. Voegler, Nucl. Phys. A 373, 192 (1982).
  * Klupfel  _et al._ (2009) P. Klupfel, P. G. Reinhard, T. J. Burvenich, and J. A. Maruhn, Phys. Rev. C 79, 034310 (2009), arXiv:0804.3385 [nucl-th] .
  * Wang  _et al._ (2013) N. Wang, L. Ou, and M. Liu, Phys. Rev. C 87, 034327 (2013).
  * Wang  _et al._ (2017) M. Wang, G. Audi, F. G. Kondev, W. J. Huang, S. Naimi, and X. Xu, Chin. Phys. C 41, 030003 (2017).
  * Angeli and Marinova (2013) I. Angeli and K. Marinova, At. Data Nucl. Data Tables 99, 69 (2013).
  * Schwierz  _et al._ (2007) N. Schwierz, I. Wiedenhover, and A. Volya, (2007), arXiv:0709.3525 [nucl-th] .
  * Patel  _et al._ (2013) D. Patel, U. Garg, M. Fujiwara, T. Adachi, H. Akimune, G. Berg, M. Harakeh, M. Itoh, C. Iwamoto, A. Long, J. Matta, T. Murakami, A. Okamoto, K. Sault, R. Talwar, M. Uchida, and M. Yosoi, Phys. Lett. B 726, 178 (2013).
  * Howard  _et al._ (2020) K. B. Howard _et al._ , Phys. Lett. B 807, 135608 (2020), arXiv:2005.10172 [nucl-ex] .
  * Youngblood  _et al._ (1999) D. H. Youngblood, H. L. Clark, and Y. W. Lui, Phys. Rev. Lett. 82, 691 (1999).
  * Brown (2013) B. A. Brown, Phys. Rev. Lett. 111, 232502 (2013).
  * Essick  _et al._ (2021) R. Essick, I. Tews, P. Landry, and A. Schwenk, Phys. Rev. Lett. 127, 192701 (2021), arXiv:2102.10074 [nucl-th] .
  * Roca-Maza  _et al._ (2013) X. Roca-Maza, M. Centelles, X. Viñas, M. Brenna, G. Colò, B. K. Agrawal, N. Paar, J. Piekarewicz, and D. Vretenar, Phys. Rev. C 88, 024316 (2013), arXiv:1307.4806 [nucl-th] .
  * Hu  _et al._ (2021) B. Hu _et al._ , (2021), arXiv:2112.01125 [nucl-th] .
  * Pineda  _et al._ (2021) S. V. Pineda _et al._ , Phys. Rev. Lett. 127, 182503 (2021), arXiv:2106.10378 [nucl-ex] .
  * Tews  _et al._ (2013) I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 110, 032504 (2013), arXiv:1206.0025 [nucl-th] .
  * Wellenhofer  _et al._ (2015) C. Wellenhofer, J. W. Holt, and N. Kaiser, Phys. Rev. C 92, 015801 (2015), arXiv:1504.00177 [nucl-th] .
  * Drischler  _et al._ (2020) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125, 202702 (2020), arXiv:2004.07232 [nucl-th] .
  * Gandolfi  _et al._ (2012) S. Gandolfi, J. Carlson, and S. Reddy, Phys. Rev. C 85, 032801 (2012), arXiv:1101.1921 [nucl-th] .
  * Wlazłowski  _et al._ (2014) G. Wlazłowski, J. W. Holt, S. Moroz, A. Bulgac, and K. J. Roche, Phys. Rev. Lett. 113, 182503 (2014), arXiv:1403.3753 [nucl-th] .
  * Roggero  _et al._ (2014) A. Roggero, A. Mukherjee, and F. Pederiva, Phys. Rev. Lett. 112, 221103 (2014), arXiv:1402.1576 [nucl-th] .
  * Akmal  _et al._ (1998) A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998), arXiv:nucl-th/9804027 .
  * Baldo and Fukukawa (2014) M. Baldo and K. Fukukawa, Phys. Rev. Lett. 113, 242501 (2014), arXiv:1409.7206 [nucl-th] .
  * Carbone  _et al._ (2014) A. Carbone, A. Rios, and A. Polls, Phys. Rev. C 90, 054322 (2014), arXiv:1408.0717 [nucl-th] .
  * Reinhard and Nazarewicz (2010) P. G. Reinhard and W. Nazarewicz, Phys. Rev. C 81, 051303 (2010), arXiv:1002.4140 [nucl-th] .
  * Piekarewicz  _et al._ (2012) J. Piekarewicz, B. K. Agrawal, G. Colo, W. Nazarewicz, N. Paar, P. G. Reinhard, X. Roca-Maza, and D. Vretenar, Phys. Rev. C 85, 041302 (2012), arXiv:1201.3807 [nucl-th] .
  * Zhang and Chen (2015) Z. Zhang and L.-W. Chen, Phys. Rev. C 92, 031301 (2015), arXiv:1504.01077 [nucl-th] .
  * Xu  _et al._ (2020) J. Xu, J. Zhou, Z. Zhang, W.-J. Xie, and B.-A. Li, Phys. Lett. B 810, 135820 (2020), arXiv:2006.05217 [nucl-th] .
  * Reinhard  _et al._ (2021b) P. G. Reinhard, B. Schuetrumpf, and J. A. Maruhn, Comput. Phys. Commun. 258, 107603 (2021b).
  * Becker  _et al._ (2018) D. Becker _et al._ , Eur. Phys. J. A 54, 208 (2018), arXiv:1802.04759 [nucl-ex] .