Commun. Theor. Phys.

Angular momentum and parity projected multidimensionally constrained
relativistic Hartree-Bogoliubov model

††footnotetext<EMAIL_ADDRESS>
Kun Wang a,b), Bing-Nan Lu c)†

a) CAS Key Laboratory of Theoretical Physics, Institute of Theoretical
Physics, Chinese Academy of Sciences, Beijing 100190, China
b) School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China
c) Graduate School of China Academy of Engineering Physics, Beijing 100193,
China

###### Abstract

The nuclear deformations are of fundamental importance in nuclear physics.
Recently we developed a multi-dimensionally constrained relativistic Hartree-
Bogoliubov (MDCRHB) model, in which all multipole deformations respecting the
$V_{4}$ symmetry can be considered self-consistently. In this work we extend
this model by incorporating the angular momentum projection (AMP) and parity
projection (PP) to restore the rotational and parity symmetries broken in the
mean-field level. This projected-MDCRHB (p-MDCRHB) model enables us to connect
certain nuclear spectra to exotic intrinsic shapes such as triangle or
tetrahedron. We present the details of the method and an exemplary calculation
for 12C. We develop a triangular moment constraint to generate the triangular
configurations consisting of three $\alpha$ clusters arranged as an
equilateral triangle. The resulting 12C spectra are consistent with that from
a triangular rigid rotor for large separations between the $\alpha$ clusters.
We also calculate the $B(E2)$ and $B(E3)$ values for low-lying states and find
good agreement with the experiments.

##  1\. Introduction

Most atomic nuclei are deformed in their ground states [1, 2]. The nuclear
shapes can be described by the surface multipole expansion

$R(\theta,\varphi)=R_{0}\left[1+\beta_{00}+\sum_{\lambda=2}^{\infty}\sum_{\mu=-\lambda}^{\lambda}\beta_{\lambda\mu}^{*}Y_{\lambda\mu}(\theta,\varphi)\right],$
(1)

where $R$ is the distance from the center of mass to the nuclear surface,
$\lambda$ and $\mu$ are integers characterizing different kinds of multipole
deformations. For example, $\beta_{\lambda 0}$’s describe axially symmetric
deformations, the triaxial deformations correspond to $\beta_{\lambda\mu}$
with even $\lambda$’s, and the tetrahedral shape means $\beta_{32}\neq 0$
while all other $\beta_{\lambda\mu}$’s vanish. These shape parameters can be
viewed as generalized coordinates, from which we can build dynamical theories
describing the rotations and vibrations of the nuclei. These collective
motions are foundations of the nuclear spectroscopy.

The nuclear spectra are the fingerprints of the underlying nuclear shapes. The
mostly studied deformations are the axially symmetric quadrupole deformation
with only $\beta_{20}\neq 0$, which manifest itself in the typical spectra of
a symmetric rigid rotor, $E\propto I(I+1)$. A more interesting case is the
triaxial deformation with nonzero $\beta_{20}$ and $\beta_{22}$, which results
in the wobbling motion [2] and the chiral doublet bands [3, 4, 5]. The
multiple chiral doublet (M$\chi$D) in nuclei with different triaxial
configurations were predicted by Meng et al. [6] and recently confirmed
experimentally in several nuclei [7, 8, 9, 10, 11]. A number of studies have
shown that the low-spin signature inversion indicates the existence of
triaxiality in certain regions of the nuclide chart [12, 13, 14, 15, 16, 17].
On the other hand, the parity-breaking octupole deformations $\beta_{3\mu}\neq
0$ are also found in many experiments [18, 19, 20]. For example, the low-lying
alternative parity bands found in actinides and some rare-earth nuclei hint
the possible existence of reflection asymmetric shapes [21, 22, 23, 24, 25,
26, 27, 28]. Direct evidences of a static reflection asymmetric shape were
found in 224Ra [29], 144Ba [30] and 228Th [31].

In recent years some exotic nuclear deformations have attracted a lot of
theoretical interests. For example, the reflection asymmetric deformations
combined with the triaxial deformations have been proposed to influence the
nuclear fission barrier [32, 33, 34, 35, 36]. The nuclear clustering has been
proved to be important in light nuclei, especially in the states near
$n\alpha$ binding threshold [37, 38, 39, 40]. The nuclear clustering
configurations with various arrangement of $\alpha$ clusters are naturally
accompanied by exotic shape deformations. It implies that a proper theoretical
description requires simultaneous inclusion of various shape degrees of
freedom. It has been demonstrated that the density functional theories are
able to describe the clustering structure in light nuclei [41, 42, 43, 44, 45,
46, 47, 48, 49]. The study in Ref. [41] discussed the clustering structure of
light nuclei with $N=Z$ by using the axially deformed relativistic mean field
model. By using the cranking covariant density functional theory (CDFT), the
rod-shaped carbon isotopes at high spins were studied [48, 50, 51]. The
studies based on the algebraic model [52] and the antisymmetrized molecular
dynamics (AMD) model [53] suggested that the $\alpha$ clustering can appear in
the low-lying states of 12C and other light even-even nuclei [54, 55, 56]. In
order to study how the nuclei cluster and how the clustering influences the
nuclear properties, we need theories that can describe both the excited states
and the ground state, and include as many shape degrees of freedom as
possible.

The CDFT gives a universal description of the nuclei all over the nuclide
chart [57, 58, 59, 60, 61, 62, 63, 64, 65]. In CDFT, the Lorentz invariance
nature of the Lagrangian automatically determines the spin-orbit interaction,
and reveals the origin of the pseudospin symmetry [66, 67, 68, 69, 70, 71, 64]
in the Fermi sea and the spin symmetry in the Dirac sea [72, 73]. An advantage
of studying nuclear clustering with CDFT is that it does not presume the
existence of $\alpha$ clusters a priori but generates such configurations from
microscopic interactions. However, the mean field approximation used in CDFT
brings the problem that the intrinsic symmetries and conservation laws such as
the translational invariance, rotational invariance, parity and particle
number conservations, are usually broken [74]. One solution is to restore
these symmetries with the projection method, examples are the angular momentum
projection (AMP) and the parity projection (PP) methods [62]. Combining AMP
and PP with the generator coordinate method, we can calculate the rotational-
vibrational spectra, which can be directly compared with the experiments [62,
75, 76, 77, 78].

In recent years, the multidimensionally constrained covariant density
functional (MDCCDFT) theories [79, 33, 80, 81] have been established and
applied successfully in several topics. In these models, the reflection
asymmetric and axially asymmetric deformations can be considered
simulteneously, with only a spatial $V_{4}$ symmetry remains unbroken for the
densities and potentials. Depending on how to deal with the pairing
correlation, two variants was developed: multidimensionally constrained
relativistic mean field model (MDCRMF) [79, 33] and multidimensionally
constrained relativistic Hartree-Bogoliubov model (MDCRHB) [80]. The MDCCDFTs
have been applied to study the potential energy surfaces and fission barriers
in actinides [79, 32, 33, 34] and superheavy nuclei [82], the third minima on
PESs of light actinides [35], the higher order deformations in superheavy
nuclei [83], the non-axial octupole $Y_{32}$ correlations in $N=150$ isotones
[84] and Zr isotopes [80], and the axial octupole $Y_{30}$ correlations in
M$\chi$D [11]. The fission dynamics in actinides was studied based on the PESs
calculated with the MDCCDFTs [85, 86, 87, 88, 89]. Furthermore, by including
the strangeness degrees of freedom, the MDCCDFTs was employed to study the
shape evolution [90, 91] and hyperon pairing correlation [92] in hypernuclei,
and build new nucleon-hyperon effective interactions [93].

In this work we focus on the MDCRHB model and its extensions. Mean field state
$\ket{\phi(q)}$ can be obtained by solving the relativistic Hartree Bogoliubov
equation in the intrinsic frame. Here $q$ is a collection of quantum numbers
specifying the deformation of the state or other intrinsic properties. Several
intrinsic symmetries may be broken due to the mean field approximation and
need to be restored using the symmetry projection methods. For example, the
angular momentum projection for a relativistic mean field model has been
implemented upon a triaxially deformed ground state [94] and an axially
deformed reflection-asymmetric ground state [95, 96]. The generating
coordinate method has also been included to treat soft potential energy
surfaces [75, 76]. In this work, we incorporate the AMP and PP methods into
the MDCRHB model. Our implementation differs from the previous works in two
aspects. First, we apply a unified framework to treat both the triaxial and
octupole deformations. While in Ref. [94] the triaxial deformations are
considered with an anisotropic 3-dimensional harmonic oscillator (3DHO) basis,
in this work we employ an axially symmetric harmonic oscillator basis in all
cases. This setting allows us to simultaneously break as many symmetries as
possible and consider exotic shapes like the triangle or tetrahedron. Second,
we use the Bogoliubov transformation together with a separable pairing force
of finite range to treat the pairing effect. This makes our results more
reliable for extreme deformations for which the pairing force can play key
roles in cluster formations.

With the projected MDCRHB model (p-MDCRHB), we are able to explore various
exotic nuclear shapes. One direct application is to search for possible
$\alpha$-clustering structure in light nuclei. As the first step toward this
direction, in this work we only consider the AMP and PP methods. The particle
number projection method and generating coordinate method also play essential
roles in some circumstances and will be considered in future works. In Sec.
2., we present the theoretical framework of the p-MDCRHB model. In Sec. 3., we
present the numerical check and the results for the $\alpha$-conjugate nucleus
12C. We summarize this work in Sec. 4..

##  2\. Theoretical Framework

The CDFTs solve the nuclear many-body problems in the framework of the density
functional theories. We design a general energy density functional
$E[\rho(\bm{r})]$ with $\rho(\bm{r})$ the intrinsic nucleon densities, then
the ground state energy and densities are given by the variational principle
$\delta E/\delta\rho=0$. In most CDFTs this framework is formally equivalent
to a mean field approximation for an effective Lagrangian, in which the
nucleons are Dirac fermions interacting indirectly via meson exchange or
directly through zero-range point-coupling interactions [97, 98]. To give
correct saturation properties of the symmetric nuclear matter, the nonlinear
self-coupling terms [99, 100, 101] or the density dependence of the coupling
constants [102, 103] were introduced and applied extensively. Accordingly, the
framework of CDFT can be formulated in four representations: effective
interactions with meson exchange (ME) or point-coupling (PC) combined with
nonlinear self-coupling (NL) or density dependence of coupling constants (DD).
In MDCCDFT we implemented all four possibilities. In this work, we illustrate
the formalism of the projected MDCRHB model using the density-dependent point-
coupling interaction DD-PC1 [104]. The formalisms for other parameters are
similar and omitted for brevity.

###  2.1. Relativistic Hartree-Bogoliubov model

In this section we briefly introduce the theoretical framework of the MDCRHB
model with density-dependent point-coupling interactions. We start with the
effective Lagrangian

$\begin{split}\mathcal{L}=&~{}\bar{\psi}({\rm
i}\gamma_{\mu}\partial^{\mu}-M)\psi\\\
&-\frac{1}{2}\alpha_{S}(\rho)\rho_{S}^{2}-\frac{1}{2}\alpha_{V}(\rho)j_{V}^{2}-\frac{1}{2}\alpha_{TV}(\rho)\vec{j}_{TV}^{2}\\\
&-\frac{1}{2}\delta_{S}(\partial_{\nu}\rho_{S})(\partial^{\nu}\rho_{S})\\\
&-eA_{\mu}j^{\mu}_{C}-\frac{1}{4}F^{\mu\nu}F_{\mu\nu},\end{split}$ (2)

where $M$ is the nucleon mass, $\alpha_{S}(\rho)$, $\alpha_{V}(\rho)$ and
$\alpha_{TV}(\rho)$ are density dependent coupling constants in scalar, vector
and isospin-vector channels, respectively. $\rho_{S},~{}j_{V}$ and
$\vec{j}_{TV}$ are the corresponding densities and four-currents. $\delta_{S}$
is the coupling constant for the derivative term. $j^{\mu}_{C}$ is the proton
current, $A_{\mu}$ is the electromagnetic field and
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the field strength
tensor. For DD-PC1 functional [104] the functions $\alpha_{S}(\rho)$,
$\alpha_{V}(\rho)$ and $\alpha_{TV}(\rho)$ share the same ansatz

$\alpha_{i}(\rho)=a_{i}+(b_{i}+c_{i}x)e^{-d_{i}x}\quad(i=S,V,TV),$ (3)

with $x=\rho/\rho_{\rm sat}$, where $\rho_{\rm sat}$ denotes the saturation
density of the symmetric nuclear matter. The constants $a_{i}$, $b_{i}$ and
$c_{i}$ are parameters fitted to the experiments.

The pairing correlations are included with the Bogoliubov transformation.
According to the mean field approximation, we take the Bogoliubov vacuum as
the trial wave function and solve the variational problem to find the ground
state. In order to prevent the system from collapsing to the negative energy
Dirac sea, we apply the “no-sea” approximation, i.e., only the positive energy
states are kept in calculating the densities and currents. The resulting RHB
equation [58, 105] in coordinate space is

$\int{\rm d}^{3}\bm{r}^{\prime}\left(\begin{array}[]{cc}h-\lambda&\Delta\\\
-\Delta^{*}&-h+\lambda\end{array}\right)\left(\begin{array}[]{c}U_{k}\\\
V_{k}\end{array}\right)=E_{k}\left(\begin{array}[]{c}U_{k}\\\
V_{k}\end{array}\right),$ (4)

where $E_{k}$ is the quasi-particle energy and $\lambda$ is the chemical
potential. $U_{k}=U_{k}(\bm{r}^{\prime}\sigma^{\prime})$ and
$V_{k}=V_{k}(\bm{r}^{\prime}\sigma^{\prime})$ are quasi-particle wave
functions. $\Delta=\Delta(\bm{r}\sigma,\bm{r}^{\prime}\sigma^{\prime})$ is the
pairing potential and $h$ is the single particle Hamiltonian

$h=\bm{\alpha}\cdot[\bm{p}-\bm{V}(\bm{r})]+\beta[M+S(\bm{r})]+V_{0}(\bm{r})+\Sigma_{R}(\bm{r}),$
(5)

where the scalar potential $S(\bm{r})$, vector potential
$V^{\mu}(\bm{r})=(V_{0}(\bm{r}),\bm{V}(\bm{r}))$ and the rearrangement term
$\Sigma_{R}(\bm{r})$ are expressed with various nuclear densities and
currents,

$\displaystyle S$
$\displaystyle=\alpha_{S}(\rho)\rho_{S}+\delta_{S}\Delta\rho_{S},$ (6a)
$\displaystyle V^{\mu}$
$\displaystyle=\alpha_{V}(\rho)j_{V}^{\mu}+\alpha_{TV}(\rho)j_{TV}^{\mu}\cdot\vec{\tau}+e\frac{1-\tau_{3}}{2}A^{\mu},$
(6b) $\displaystyle\Sigma_{R}$
$\displaystyle=\frac{1}{2}\frac{\partial\alpha_{S}}{\partial\rho}\rho_{S}^{2}+\frac{1}{2}\frac{\partial\alpha_{V}}{\partial\rho}j_{V}^{2}+\frac{1}{2}\frac{\partial\alpha_{TV}}{\partial\rho}\vec{j}_{TV}^{2}.$
(6c)

The pairing potential $\Delta$ is written as

$\begin{split}&\Delta_{p_{1}p_{2}}(\bm{r}_{1}\sigma_{1},\bm{r}_{2}\sigma_{2})\\\
&=\int{\rm d}^{3}\bm{r}^{\prime}_{1}{\rm
d}^{3}\bm{r}^{\prime}_{2}\sum_{\sigma^{\prime}_{1}\sigma^{\prime}_{2}}^{p^{\prime}_{1}p^{\prime}_{2}}V_{p_{1}p_{2},p^{\prime}_{1}p^{\prime}_{2}}^{\mathrm{pp}}(\bm{r}_{1}\sigma_{1},\bm{r}_{2}\sigma_{2},\bm{r}^{\prime}_{1}\sigma^{\prime}_{1},\bm{r}^{\prime}_{2}\sigma^{\prime}_{2})\times\kappa_{p^{\prime}_{1}p^{\prime}_{2}}(\bm{r}^{\prime}_{1}\sigma^{\prime}_{1},\bm{r}^{\prime}_{2}\sigma^{\prime}_{2})\end{split}$
(7)

with the pairing tensor
$\kappa=\langle\hat{a}(\bm{r}_{1}\sigma_{1})\hat{a}(\bm{r}_{2}\sigma_{2})\rangle$.
In Eq. (7) and what follows we use the symbols with and without prime to
denote the quantities for initial and final states, respectively. As a
standard procedure in RHB models [106], we only include the large component of
the Dirac spinors in calculating the pairing tensor.

In this work we use a separable pairing force of finite range [107, 108]

$V(\bm{r}_{1}\sigma_{1},\bm{r}_{2}\sigma_{2},\bm{r}^{\prime}_{1}\sigma^{\prime}_{1},\bm{r}^{\prime}_{2}\sigma^{\prime}_{2})=-G\delta(\bm{R}-\bm{R}^{\prime})P(r)P(r^{\prime})\frac{1}{2}(1-P_{\sigma}),$
(8)

where $G$ is the pairing strength, $(1-P_{\sigma})/2$ is the projector onto
the total spin $S=0$, $\bm{R}=(\bm{r}_{1}+\bm{r}_{2})/2$ and
$\bm{r}=\bm{r}_{1}-\bm{r}_{2}$ are center of mass coordinate and relative
coordinate, respectively. $P(r)$ is a normalized Gaussian function $P(r)=(4\pi
a^{2})^{-3/2}\exp\left(-r^{2}/(4a^{2})\right)$ with $a$ an adjustable
parameter specifying the range of the pairing force. The details of
calculating the pairing tensor and pairing potential in the ADHO basis can be
found in Ref. [80]. In this illustrative work we fix the pairing strength to
$G=1.1G_{0}$ with $G_{0}=728.0$ MeV$\cdot$fm3 and the range to $a=0.644$ fm,
in which $G_{0}$ and $a$ are determined by fitting to the pairing gap in the
nuclear matter [107, 109].

The RHB equation (4) can be solved self-consistently starting from a properly
chosen initial state. The resulting minima represent the ground state or the
shape isomers. To obtain the potential energy surfaces $E(\beta)$ we use a
modified linear constraint method [33]. The Routhian with deformation
constraints reads

$E^{\prime}=E_{\mathrm{MF}}+\sum_{\lambda\mu}\frac{1}{2}C_{\lambda\mu}\beta_{\lambda\mu},$
(9)

where $E_{\rm{MF}}$ is the energy functional obtained with the Lagrangian (2)
and the pairing force (8). The variable $C_{\lambda\mu}$ varies during the
iteration. $\beta_{\lambda\mu}$ is the nuclear deformation

$\beta_{\lambda\mu}=\frac{4\pi}{3AR^{\lambda}}\langle Q_{\lambda\mu}\rangle,$
(10)

where $R=1.2A^{1/3}$ fm is the approximate nuclear radius, $A$ is the number
of nucleons and $Q_{\lambda\mu}=r^{\lambda}Y_{\lambda\mu}(\Omega)$, where
$Y_{\lambda\mu}$ is the spherical harmonics, is the multipole moment operator.

###  2.2. Axially symmetric harmonic oscillator basis

In the MDCRHB model, we solve the RHB equation (4) by expanding the spinor
wave functions $U_{k}(\bm{r}\sigma)$ and $V_{k}(\bm{r}\sigma)$ in an axially
deformed harmonic oscillator (ADHO) basis [110, 33, 80],

$U_{k}(\bm{r}\sigma)=\left(\begin{array}[]{c}\sum_{\alpha}f_{U}^{k\alpha}\Phi_{\alpha}(\bm{r}\sigma)\\\
\sum_{\alpha}g_{U}^{k\alpha}\Phi_{\alpha}(\bm{r}\sigma),\end{array}\right)$
(11a)
$V_{k}(\bm{r}\sigma)=\left(\begin{array}[]{c}\sum_{\alpha}f_{V}^{k\alpha}\Phi_{\alpha}(\bm{r}\sigma)\\\
\sum_{\alpha}g_{V}^{k\alpha}\Phi_{\alpha}(\bm{r}\sigma),\end{array}\right)$
(11b)

where $\Phi_{\alpha}(\bm{r}\sigma)$ is the basis wave function generated by
solving the Schrodinger equation with an ADHO potential. Here we impose the
$V_{4}$ symmetry, which means that the densities and potentials are invariant
under the mirror reflections

$\begin{split}\hat{S}_{x}\phi(x,y,z)&=\phi(-x,y,z),\\\
\hat{S}_{y}\phi(x,y,z)&=\phi(x,-y,z),\end{split}$ (12)

where $\phi$ represents any scalar field or spatial distribution in the
theory. It immediately follows that the rotation
$\hat{S}_{z}(\pi)=\hat{S}_{x}\hat{S}_{y}$ around the $z$-axis by $\pi$ is also
a symmetry operation.

In the ADHO basis, the reflection asymmetric shapes can be included by mixing
the basis states with different parities. While the axially deformed shapes
preserve the $z$-axis angular momentum projection $K$, the triaxial shapes mix
different $K$ blocks and we need to calculate the corresponding mixing matrix
elements. In this case, we expand the densities and potentials in terms of the
Fourier series of the azimuthal coordinate $\varphi$ in cylindrical
coordinates. Under the restriction of the $V_{4}$ symmetry, the expansion for
any field $f$ writes

$f(\rho,\varphi,z)=f_{0}(\rho,z)\frac{1}{2\sqrt{\pi}}+\sum_{n=1}^{\infty}f_{n}(\rho,z)\frac{1}{\sqrt{\pi}}\cos(2n\varphi).$
(13)

The ADHO basis are also eigenfunctions of the simplex operator
$\hat{S}_{z}(\pi)$ with eigenvalues $S=\pm 1$. The basis with $S=-1$ can be
calculated with a time-reversal operation on those with $S=1$. For a system
with time-reversal symmetry, the single particle energy levels have the
Krammer degeneracy and we only need to diagonalize the Hamiltonian in the
$S=1$ subspace.

In mean field calculations we need to calculate the potentials on a
cylindrical lattice. The number of mesh points in $z$\- and $\rho$\-
directions are $n_{z}=24$ and $n_{\rho}=12$, respectively. In MDCCDFT, the
triaxial deformation is included by the Fourier transformation given in Eq.
(13). Choosing $z$ as the principal axis, only the $n=0$ Fourier component
survives in the axially symmetric case, which simplifies the calculation
regardless of $\varphi$. A triangular deformed system lying on the $x$-$z$
plane is not axially symmetric about the intrinsic $z$-axis. In such a system,
$n>1$ components in the Fourier transformation needs to be taken into account.
Due to the $V_{4}$ symmetry restriction, we only need to calculate the
components with even $n$ in the Fourier transformation. The $V_{4}$ symmetry
can simplify the calculation into $1/4$ of the full space with respect to
$\varphi$ and we choose the number of points for $\varphi$ integral as
$n_{\varphi}=12$ for $0<\varphi<\pi/2$.

###  2.3. Angular momentum and parity projections

The angular momentum projection operator for a general wave function breaking
both the axial symmetry and reflection symmetry can be written as

$\hat{P}_{MK}^{I}=\frac{2I+1}{8\pi^{2}}\int{\rm d}\Omega
D_{MK}^{I*}(\Omega)\hat{R}(\Omega),$ (14)

where $D_{MK}^{I}(\Omega)={\rm e}^{-iM\alpha}d_{MK}^{I}(\beta){\rm
e}^{-iK\gamma}$ is the Wigner function with Euler angle
$\Omega\equiv\\{\alpha,\beta,\gamma\\}$ and $\hat{R}(\Omega)={\rm
e}^{-i\alpha\hat{I}_{z}}{\rm e}^{-i\beta\hat{I}_{y}}{\rm
e}^{-i\gamma\hat{I}_{z}}$ is the rotational operator applied to the many-body
wave function. $I$ is the total angular momentum, $M$ and $K$ are angular
momentum projections along the $z$-axis in the laboratory frame and intrinsic
frame, respectively. For nuclei with reflection asymmetric shapes, the parity
projection should also be included to build states with definite parities
[74]. The parity projection operator writes

$\hat{P}^{\pi}=\frac{1}{2}(1+\pi\hat{P}),$ (15)

where $\hat{P}$ is the spatial reflection operation and $\pi=\pm 1$ is the
parity.

With $\hat{P}^{I}_{MK}$ and $\hat{P}^{\pi}$, we can project out states with
good angular momentum $I$, $M$ and parity $\pi$ from a RHB mean field wave
function $|\Phi\rangle$,

$\ket{\Psi_{M}^{I\pi}}=\sum_{K}f_{K}^{I\pi}\hat{P}^{\pi}\hat{P}_{MK}^{I}\ket{\Phi}=\sum_{K}f_{K}^{I\pi}\ket{IMK\pi},$
(16)

where $f^{I\pi}_{K}$ are some undetermined constants. Taking Eq. (16) as the
trial wave function and solving the variational problem against the variables
$f^{I\pi}_{K}$, we end up with the Hill-Wheeler-Griffin equation [74]

$\sum_{K^{\prime}}[\mathcal{H}_{KK^{\prime}}^{I\pi}-E_{\alpha}^{I\pi}\mathcal{N}_{KK^{\prime}}^{I\pi}]f_{K^{\prime}\alpha}^{I\pi}=0,$
(17)

where the Hamiltonian kernel $\mathcal{H}$ and the norm kernel $\mathcal{N}$
are defined as

$\begin{split}\mathcal{H}_{KK^{\prime}}^{I\pi}&=\bra{\Phi}HP^{\pi}P_{KK^{\prime}}^{I}\ket{\Phi}\\\
&=\frac{2I+1}{16\pi^{2}}\int\left[D_{KK^{\prime}}^{I^{*}}(\Omega)+\pi(-1)^{I}D_{K-K^{\prime}}^{I}(\Omega)\right]\bra{\Phi}\hat{H}\hat{R}(\Omega)\ket{\Phi}{\rm
d}\Omega,\end{split}$ (18a)
$\begin{split}\mathcal{N}_{KK^{\prime}}^{I\pi}&=\bra{\Phi}P^{\pi}P_{KK^{\prime}}^{I}\ket{\Phi}\\\
&=\frac{2I+1}{16\pi^{2}}\int\left[D_{KK^{\prime}}^{I^{*}}(\Omega)+\pi(-1)^{I}D_{K-K^{\prime}}^{I}(\Omega)\right]\bra{\Phi}\hat{R}(\Omega)\ket{\Phi}{\rm
d}\Omega.\end{split}$ (18b)

Note that the expressions of the kernels are derived by considering the
$V_{4}$ symmetry. By discretizing the group space, the problem turns to the
calculation of the norm overlap $n(\Omega)$ and the Hamiltonian overlap
$h(\Omega)$ as functions of the Euler angle $\Omega$,

$\displaystyle n(\Omega)$ $\displaystyle=\bra{\Phi}\hat{R}(\Omega)\ket{\Phi},$
(19a) $\displaystyle h(\Omega)$
$\displaystyle=\bra{\Phi}\hat{H}\hat{R}(\Omega)\ket{\Phi}/n(\Omega).$ (19b)

We adopt the method introduced in Ref. [94] to calculate the overlaps in the
ADHO basis using the generalized Wick’s theorem.

We follow the standard way of solving the Hill-Wheeler-Griffin equation (17)
as discussed in Ref. [74]. The first step is to diagonalize the norm kernel

$\sum_{K^{\prime}}\mathcal{N}_{KK^{\prime}}^{I\pi}u_{m}^{IK^{\prime}\pi}=n_{m}^{I\pi}u_{m}^{IK\pi}.$
(20)

The normalized vectors can be used to build orthonormalized basis

$\ket{m}=\frac{1}{\sqrt{n_{m}^{I\pi}}}\sum_{K}u_{m}^{IK\pi}\hat{P}^{\pi}\hat{P}_{MK}^{I}\ket{\Phi}.$
(21)

Eigen vectors with small eigen values $n^{I\pi}_{m}$ may induce numerical
instabilities. To solve this problem, we introduce a cutoff $\chi=10^{-3}$ and
only consider eigenstates with eigenvalues $n_{m}^{I\pi}>\chi$. In the
othornormalized basis $|m\rangle$ the Hill-Wheeler-Griffin equation becomes a
standard eigenvalue problem,

$\sum_{m^{\prime}}\braket{m}{\hat{H}}{m^{\prime}}g_{m^{\prime}}^{I\pi}=E^{I\pi}g_{m}^{I\pi}.$
(22)

The eigen energies $E_{\alpha}^{I\pi}$ can be solved by matrix diagonalization
and the wave function components $f^{I\pi}_{K}$ can be obtained from the
corresponding eigenvectors $g^{I\pi}_{m}$,

$f_{K}^{I\pi}=\sum_{m,n_{m}^{I\pi}\neq
0}\frac{g_{m}^{I\pi}}{\sqrt{n_{m}^{I\pi}}}u_{m}^{IK\pi}.$ (23)

With the projected wave functions, the reduced transition probabilities
$B(E\lambda)$ from $J_{i}$ state to $J_{f}$ state writes

$\begin{split}B(E\lambda;J_{i}\rightarrow
J_{f})&=\frac{1}{2J_{i}+1}\sum_{M_{i}M_{f}\mu}\big{|}\bra{J_{f},M_{f}}\hat{Q}_{\lambda\mu}\ket{J_{i},M_{i}}\big{|}^{2}\\\
&=\frac{1}{2J_{i}+1}\big{|}\bra{J_{f}|}\hat{Q}_{\lambda}\ket{|J_{i}}\big{|}^{2}.\end{split}$
(24)

The reduced matrix element has the following expression

$\bra{J_{f}|}\hat{Q}_{\lambda}\ket{|J_{i}}=\sum_{K_{i}K_{f}}f_{K_{f}}^{I_{f}\pi_{f}}f_{K_{i}}^{I_{i}\pi_{i}}\left[Q_{\lambda}\right]_{K_{f}K_{i}}^{J_{f}J_{i}},$
(25)

where

$\left[Q_{\lambda}\right]_{K_{f}K_{i}}^{J_{f}J_{i}}=\frac{(2J_{i}+1)(2J_{f}+1)}{8\pi^{2}}(-1)^{J_{f}-K_{f}}\sum_{\mu
M}\left(\begin{array}[]{ccc}J_{f}&\lambda&J_{i}\\\
-K_{f}&\mu&M\end{array}\right)Q_{\lambda\mu},$ (26)

in which

$Q_{\lambda\mu}=\int{\rm d}\Omega
D_{MK_{i}}^{J_{i}*}(\Omega)\bra{\Phi}\hat{Q}_{\lambda\mu}\hat{R}(\Omega)\ket{\Phi}.$
(27)

##  3\. Results and Discussions

In this section, we apply the p-MDCRHB model to study the low-lying
excitations of the nucleus 12C, which is a typical nucleus with possible
$\alpha$ clustering structures. Among the proposed clustering states, the
second $0^{+}$ state of 12C, namely the Hoyle state, is the most interesting
because it plays a crucial role in the stellar nuclear synthesis [111]. The
Hoyle state is a resonance consisting of three interacting $\alpha$ clusters
lying above the three-$\alpha$ breakup threshold [112]. The energy and width
of this resonance play an essential role in synthesizing three $\alpha$’s into
12C in the main sequence stars. The abundance of the carbon elements would be
insufficient for supporting life in our universe if the Hoyle state has a
different energy or does not exist at all. Furthermore, several recent studies
of the ground state band indicate that the 3-$\alpha$ structure may already
appear in the ground state of 12C [113, 53]. A $5^{-}$ state was measured in
the ground state band [54], providing an evidence of a triangular $D_{3h}$
symmetry.

In MDCRHB model and its extensions, the system respects the mirror symmetries
about the $x$-$z$ and $y$-$z$ planes. One interesting case which these models
specialize is the clustering structure in the light $\alpha$-conjugate nuclei.
Here the $\alpha$ clusters can form various molecular-like structures and the
possible shapes include rod, triangle, tetrahedron, octahedron, etc. Most of
these shapes can not be described by conventional multipole expansion method
that only includes either triaxial or octupole deformations. Further, the
structure might be different for ground state and excited state, a mean field
description is too simplified and we must take into account the symmetry
restorations with projection techniques.

In this section, we apply the p-MDCRHB model to study the ground state and
low-lying excitations of 12C nucleus with an intrinsic triangular shape. To
generate such a configuration, we define a triangular deformation constraint
in the mean field calculations. In analogy with the usual multipole
deformation constraint Eq. (10), we use the ansatz

$\lambda_{3}=\frac{4\pi}{3NR^{3}}\langle S_{3}\rangle,$ (28)

where the “triangular moment operator” is defined as

$S_{3}=(z^{2}+\rho^{2}\cos^{2}\varphi)^{3/2}\cos(3\theta),$ (29)

with $\rho$, $\varphi$ and $z$ the cylindrical coordinates and $\theta$ is the
polar angle on the $x$-$z$ plane. We can verify that Eq. (29) has the desired
symmetry. If rotated around $y$-axis, the operator $S_{3}$ has a period
$2\pi/3$ and the maxima point to the directions along which the nucleus is
stretched. The parameter $\lambda_{3}$ is the counterpart of the multipole
deformation parameters $\beta_{\lambda\mu}$. When constrained to large
$\lambda_{3}$, the operator $S_{3}$ elongates the nucleus in directions
$\theta=0,2\pi/3,4\pi/3$ and compress it in other directions. The resulting
shape would be a regular triangle with three 4He at the vertices. In MDCRHB
model, the triangular constraint can be implemented in exactly the same way
with the multipole moment constraint. In Fig. 1 we show the density profiles
of 12C constrained to $\lambda_{3}=0.8$. We present the contours in all three
coordinate planes $x$=0, $y$=0, $z$=0 and a 3-D isosurface. The densities are
consistent with what we expect from the structure of the triangular moment
constraint. The energetically preferred directions are distributed evenly in
the $x$-$z$ plane. With much stronger constraining potential individual
$\alpha$ clusters can form and be squeezed out of the nucleus.

Figure 1: Density profiles of 12C constrained for $\lambda_{3}=0.8$. The
contour plots display the total density $\rho_{\mathrm{total}}$ in $x$-$z$,
$y$-$z$ and $x$-$y$ planes. The 3-D plot shows the isosurface of
$\rho_{\mathrm{total}}$=0.08 fm-3.

In mean field calculations, the Dirac equation is solved by expanding the wave
functions in an ADHO basis with a major shell truncation $N_{f}$. For Dirac
spinors we use $N_{f}$ shells to expand the large component and $N_{f}+1$
shells to expand the small component in order to remove the spurious states.
Firstly, we examine the convergence of the mean field results against the
basis truncation $N_{f}$. In Fig. 2 we present the potential energy surfaces
(PESs) of 12C calculated with the triangular constraint using different values
of $N_{f}$. We found that all energy minima locate at $\lambda_{3}=0$
corresponding to a spherical shape. There is no explicit triangular
correlation in the mean field level. As $N_{f}$ increases, differences between
adjacent curves become smaller, signifying a clear convergence pattern. At the
ground state $\lambda_{3}=0$, the calculated binding energy already converges
at $N_{f}=8$, with the truncation error smaller than 0.1 MeV estimated by
comparing to the result at $N_{f}=12$. For larger $\lambda_{3}$, the mean
field potential assumes a triangular shape, differing from that of the
harmonic oscillator basis. In this case we observe a slower convergence. For
example, the mean field energies for $N_{f}=10$ and $N_{f}=12$ differ by about
1 MeV for $\lambda_{3}\approx 3.0$. Nevertheless, the qualitative behaviours
of the PES’s are not changed when $N_{f}$ is large enough. Concerning the big
computational cost of the 3-D angular momentum projection, in this
illustrative work we use $N_{f}=8$ in all following calculations. For more
realistic calculations we need to either employ a triangular harmonic
oscillator basis or increase $N_{f}$ until all truncation errors are smaller
than the desired precision.

Figure 2: The energy of 12C as functions of the triangular shape parameter
$\lambda_{3}$ calculated with different basis truncations in MDCRHB model. The
squares, circles, up triangles and down triangles denote the results obtained
with major shell truncation $N_{f}$=6, 8, 10 and 12, respectively. The center
of mass correction to the binding energy is not included.

For angular momentum and parity projections, the Hamiltonian kernel Eq. (18a)
and the norm kernel Eq. (18b) are calculated by integrating over the SO(3)
group space parameterized by the Euler angles $(\alpha,\beta,\gamma)$. These
integrals are performed numerically by the Gauss-Legendre integral with
numbers of mesh points $(n_{\alpha},n_{\beta},n_{\gamma})$ in $\alpha$,
$\beta$ and $\gamma$ directions, respectively. As most of the computational
cost comes from the calculation of the Hamiltonian overlap and norm overlap
for each Euler angle, it is desirable to reduce the total number of integral
mesh points for a given precision. Here we make a convergence check against
the number of mesh points in the group space. In Fig. 3, we show the projected
energies as functions of $n_{\alpha}$, $n_{\beta}$ and $n_{\gamma}$. On the
left panel, we set $n_{\beta}=8$ and study the projected energies for $6\leq
n_{\alpha}=n_{\gamma}\leq 14$, where $E^{J\pi}$ denotes the energy of the
lowest $\ket{J^{\pi}}$ state with corresponding Euler angle mesh size, and
$E^{J\pi}_{0}$ represents the energy with the largest mesh size. we see that
the influence of $n_{\alpha}$ and $n_{\gamma}$ on the projected energies for
$J<6$ is less than 10 keV, which is negligible compared with the excitation
energies of the order of MeV. In contrast, the right panel shows that the
number $n_{\beta}$ has a rather large influence on the projected energies. All
projected energies up to $6^{+}$ achieve relative precisions less than
$10^{-6}$ when $n_{\beta}\geq 10$. In this work, we choose the mesh size
$n_{\alpha}=n_{\gamma}=12$ and $n_{\beta}=14$, with which the numerical
integral in the group space can be evaluated rather precisely for $J\leq 6$.

Figure 3: The projected energies of 12C calculated with different mesh size
for SO(3) group integral. For each curve, the energies are relative to the
corresponding converged values $E_{0}^{J}$ calculated at large mesh sizes. All
calculations are performed with a mean field wave function constrained to
$\lambda_{3}=1.2$. $n_{\alpha}$, $n_{\beta}$ and $n_{\gamma}$ are number of
mesh points for integrating Euler angles $\alpha$, $\beta$ and $\gamma$,
respectively. Squares, circles, up triangles, down triangles, diamonds and
left triangles denote results projected to $J^{\pi}$=$2^{+}$, $4^{+}$,
$6^{+}$, $3^{-}$, $4^{-}$, $5^{-}$, respectively.

The constraint Eq. (29) produces a triangular shape with the 3-fold axis
pointing along the $y$-axis and one of the 2-fold axes pointing along the
$z$-axis. For angular momentum projection we rotate the mean field wave
function to Euler angles $(\alpha,\beta,\gamma)$ and calculate its overlap
with the original wave function. These norm overlaps should exhibit
periodicity if the rotations are about the symmetry axes of the intrinsic
shape, which can be examined numerically. In Fig. 4 we show the norm overlaps
for different Euler angles. All results are calculated with a mean field wave
function constrained to $\lambda_{3}=1.2$. In left panels we set
$\beta=\gamma=0$ and rotate the nucleus around the $z$-axis by increasing
$\alpha$ from $0^{\circ}$ to $360^{\circ}$. As the $z$-axis is a 2-fold axis
of the triangle, we found a period of $180^{\circ}$ for both the norm overlap
and Hamiltonian overlap. In right panels we set $\alpha=\gamma=0$, then
variation of $\beta$ corresponds to a rotation around the $y$-axis. In this
case we observe a period of $120^{\circ}$ for both overlaps, which is
consistent with the fact that $y$-axis is a 3-fold axis of the triangle. In
all subplots, the overlaps return to its original values at
$\alpha,\beta,\gamma=0$ after rotating by the corresponding period. This is
the direct consequence of the invariance of the nuclear shape under these
symmetry operations.

Figure 4: The norm and Hamiltonian overlaps of 12C as functions of the Euler
angles $\alpha$ and $\beta$, calculated with a mean field wave function
constrained to $\lambda_{3}=1.2$. $\beta=\gamma=0$ on the left panels and
$\alpha=\gamma=0$ on the right panels.

The spectrum of a triangular rotor can be mostly determined by the group
theory analysis. The triangle shape corresponds to a point group $D_{3h}$,
which has 12 elements and 6 irreducible representations (irreps)
$A_{1}^{\prime}$, $A_{2}^{\prime}$, $E^{\prime}$, $A_{1}^{\prime\prime}$,
$A_{2}^{\prime\prime}$ and $E^{\prime\prime}$. As $D_{3h}$ is a subgroup of
the $O(3)$ group, we can decompose the irreps of $O(3)$ group into irreps of
the $D_{3h}$ group [114, 115, 116]. The $O(3)$ group has infinitely many
irreps corresponding to different angular momentum $I$ and parity $\pi$. In
Table 1 we present the decomposition of the $O(3)$ irreps into the $D_{3h}$
irreps. In each column we show the number of each $D_{3h}$ irrep appearing in
the decomposition. For example, the identity representation corresponds to
$0^{+}$ for $O(3)$ group and is marked by $A_{1}^{\prime}$ for $D_{3h}$ group,
while the irrep of $O(3)$ with $I^{\pi}=1^{+}$ can be written as a direct sum
of $A_{2}^{\prime}$ and $E^{\prime\prime}$ representations.

For any state of the triangular rotor with quantum numbers $I^{\pi}$, the wave
function transforms under the corresponding irrep of the $O(3)$ group. If we
rotate the rotor with an element of the $D_{3h}$ group in the intrinsic frame,
nothing is actually changed, thus the wave functions should be scalars for
these rotations. In the language of the group theory, the allowed states all
belong to the $A_{1}^{\prime}$ irrep of the intrinsic $D_{3h}$ group. If an
angular momentum and parity combination $I^{\pi}$ does not contain any
$A_{1}^{\prime}$ irrep, it should not appear in the spectrum. In other words,
the $I^{\pi}$ states with non-zero entries in the first row of Table 1 are all
states allowed by the symmetry analysis. In summary, the lowest energy levels
of a triangular rigid rotor should be
$I^{\pi}=0^{+},2^{+},3^{-},4^{+},4^{-},5^{-},6^{+},6^{-},...$. Other angular
momentum and parity combinations should be consequences of the intrinsic
excitations.

Table 1: The decomposition of the irreps of the $O(3)$ group with quantum numbers $I^{\pi}$ into irreps of the $D_{3h}$ group. In each column the numbers are the degeneracies.  $I^{\pi}$ | $0^{+}$ | $1^{+}$ | $2^{+}$ | $3^{+}$ | $4^{+}$ | $5^{+}$ | $6^{+}$ | $7^{+}$ | $0^{-}$ | $1^{-}$ | $2^{-}$ | $3^{-}$ | $4^{-}$ | $5^{-}$ | $6^{-}$ | $7^{-}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$A^{\prime}_{1}$ | 1 | 0 | 1 | 0 | 1 | 0 | 2 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1
$A^{\prime}_{2}$ | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 2 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1
$E^{\prime}$ | 0 | 0 | 1 | 1 | 2 | 2 | 2 | 2 | 0 | 1 | 1 | 1 | 1 | 2 | 2 | 3
$A^{\prime\prime}_{1}$ | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 2 | 1
$A^{\prime\prime}_{2}$ | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 2
$E^{\prime\prime}$ | 0 | 1 | 1 | 1 | 1 | 2 | 2 | 3 | 0 | 0 | 1 | 1 | 2 | 2 | 2 | 2

On the one hand, the group theory analysis is strict and can be generalized to
any other intrinsic shapes. For example, we can predict the low-lying states
of a tetrahedral rotor and compare the results with that of the 16O nucleus,
which was predicted to exhibit a four-$\alpha$ cluster structure. On the other
hand, there are also less abstract ways of understanding the essentiality of
the nuclear rotational spectrum. For a triangular rotor, a stable rotation can
only occur around its symmetry axes. The wave function for angular momentum
$I$ is $\psi(\varphi)=\exp(iI\varphi)$ with $\varphi$ the rotational angle.
For a $n$-fold axis we must have $\psi(\varphi+2\pi/n)=\psi(\varphi)$, thus
among all integer values of $I$, only integer multiples of $n$ are allowed.
For example, the rotations about a 2-fold axis generate angular momenta
$I=0,2,4,6,...$, while the rotations about a 3-fold axis give the angular
momenta $I=0,3,6,...$. Further, the parity of the rotational wave function can
be determined by $\pi=(-1)^{I}$. The total rotational energy of a triangular
rotor is the sum of that around the 2-fold and 3-fold axes. If there is no
rotation around the 3-fold axis, we have the spectrum
$I^{\pi}=0^{+},2^{+},4^{+},6^{+},...$ If there is one unit of quantized
angular momentum around the third axis, we have states
$I^{\pi}=3^{-},4^{-},5^{-},6^{-},...$. Here the states $4^{-}$ and $5^{-}$ can
be understood as a rotation around the axis in between the 2-fold axis and the
3-fold axis, which is already a complex motion in the classical level.

Next we turn to the projected potential energy surfaces (PES’s). We first
generate mean field configurations with $\lambda_{3}$ varying from 0.0 to 5.0.
The results describe the disintegration of 12C into three $\alpha$ clusters.
In Fig. 5 we show the projected energies as well as mean field energies for
all $\lambda_{3}$ in this interval with a step size $\Delta\lambda_{3}=0.2$.
Here we only show the lowest energy levels with $I^{\pi}=0^{+}$, $2^{+}$,
$3^{-}$, $4^{\pm}$ and $5^{-}$. The mean field PES gives a ground state at the
spherical shape $\lambda_{3}=0$. The energy increases quickly as three
$\alpha$’s are separated. For $\lambda_{3}=1.0$ the deformation energy is more
than 5 MeV, preventing the dissociation of the nucleus. By projecting out the
spurious rotational energies, we see that the PES for $0^{+}$ state is quite
different from the mean field PES. The $0^{+}$ PES is lower than the mean
field PES and becomes very flat for $\lambda_{3}<2.0$. We find a minimum at
$\lambda_{3}=1.2$, which is very shallow that the energy at $\lambda_{3}=0$ is
only about 1 MeV higher. In the insets of Fig. 5 we show the density profiles
from the mean field calculations at $\lambda_{3}=1.2$ and $\lambda_{3}=3.8$.
The first inset corresponds to the ground state identified with the $0^{+}$
PES. Although at this point there is no emerging clustering structure, we can
see a clear triangular shape from the total densities. At $\lambda_{3}=3.8$
the deformation is so strong that the clustering structure is very pronounced.
Theoretically for extremely large $\lambda_{3}$ the 12C nucleus would break up
into three $\alpha$ particles and the energy should converge to three times
the 4He binding energy. However, for such light nuclei we might need a very
large center of mass correction and the mean field model might not be a good
approximation. Even though there are still a lot of difficulties, the
extremely large $\lambda_{3}$ deformation still deserves some further
investigations. This three-$\alpha$ disintegration can be viewed as an inverse
reaction of the famous $3\alpha\rightarrow^{12}$C reaction and may be related
to the Hoyle state. Whether the mean field model is adequate for describing
such a process and what is lack is thus an interesting topic. Although in
recent years the ab initio calculations are more popular in studying the
clustering structures, mean field models have a clear physical pictures and
have proven successful on this topic. The p-MDCRHB model we presented here can
describe various $\alpha$ cluster configurations and can serve as a starting
point for beyond-mean-field calculations in this direction.

The PES’s for higher angular momenta tell us how the $\alpha$ clusters move
about the rotational axes. At $\lambda_{3}=0$ the system is spherical and
there is no rotation. For larger $\lambda_{3}$, the $2^{+}$ and $4^{+}$ states
describe rotations around the 2-fold axes, the corresponding PES’s have mimina
at relatively small $\lambda_{3}$. The $3^{-}$ and $4^{-}$ states contain
rotations around the 3-fold axis, thus the centrifugal force prefers a larger
$\lambda_{3}$ which also means a larger moment of inertia.

Figure 5:  The mean field and the projected potential energies of 12C versus the triangular deformation $\lambda_{3}$. The squares, circles, down triangles, up triangles, diamonds and left triangles represent $I^{\pi}=0^{+}$, $2^{+}$, $3^{-}$, $4^{+}$, $4^{-}$ and $5^{-}$, respectively. The mean field energies are displayed as black solid line. The insets show the mean field density profiles at $\lambda_{3}=1.2$ and $3.8$. Table 2: Reduced $E2$ transition probabilities from states $I^{\pi}_{i}$ to states $I^{\pi}_{f}$ and the corresponding excitation energies of 12C. The calculation is made for $\lambda_{3}$=1.2 corresponding to the energy minimum of the $0^{+}$ PES shown in Fig. 5. The experimental $B(E2)$ value is taken from Ref. [117]. $I^{\pi}_{i}$ |  $E_{xi}^{\mathrm{th}}$ (MeV) |  $E_{xi}^{\mathrm{expt}}$ (MeV) | $I^{\pi}_{f}$ |  $E_{xf}^{\mathrm{th}}$ (MeV) |  $B(E2)^{\mathrm{th}}$ (e2 fm4) |  $B(E2)^{\mathrm{expt}}$ (e2 fm4)
---|---|---|---|---|---|---
$2^{+}$ | 1.951 | 4.440 | $0^{+}$ | 0.000 | 9.72 | 7.6 $\pm$ 0.4 [117]
$4^{+}$ | 9.015 | 14.079 | $2^{+}$ | 1.951 | 16.28 |
$4^{-}$ | 15.253 |  | $3^{-}$ | 12.864 | 20.95 |
$5^{-}$ | 19.915 |  | $3^{-}$ | 12.864 | 6.78 |
|  |  | $4^{-}$ | 15.253 | 18.90 |
$6^{+}_{1}$ | 26.734 |  | $4^{+}$ | 9.015 | 26.04 |
$6^{-}$ | 25.131 |  | $4^{-}$ | 15.253 | 13.38 |
|  |  | $5^{-}$ | 19.915 | 15.01 |
Table 3: Reduced $E3$ transition probabilities from states $I^{\pi}_{i}$ to states $I^{\pi}_{f}$ and the corresponding excitation energies of 12C. The calculation is made for $\lambda_{3}$=1.2 corresponding to the energy minimum of the $0^{+}$ PES shown in Fig. 5. The experimental $B(E3)$ value is taken from Ref. [117]. $I^{\pi}_{i}$ |  $E_{xi}^{\mathrm{Th}}$ (MeV) |  $E_{xi}^{\mathrm{Expt}}$ (MeV) | $I^{\pi}_{f}$ |  $E_{xf}^{\mathrm{Th}}$ (MeV) |  $B(E3)^{\mathrm{Th}}$ (e3 fm6) |  $B(E3)^{\mathrm{Expt}}$ (e3 fm6)
---|---|---|---|---|---|---
$3^{-}$ | 12.864 | 9.641 | $0^{+}$ | 0.000 | 145.81 | 103 $\pm$ 17 [117]
|  |  | $2^{+}$ | 1.951 | 260.25 |
|  |  | $4^{+}$ | 9.015 | 53.53 |
$4^{-}$ | 15.253 | 13.35 | $2^{+}$ | 1.951 | 268.69 |
|  |  | $4^{+}$ | 9.015 | 198.55 |
$5^{-}$ | 19.915 | 22.4(2) | $2^{+}$ | 1.951 | 96.28 |
|  |  | $4^{+}$ | 9.015 | 318.35 |
$6^{+}_{1}$ | 26.734 |  | $4^{-}$ | 15.253 | 33.95 |
|  |  | $5^{-}$ | 19.915 | 129.86 |
|  |  | $6^{-}$ | 25.131 | 250.19 |
$6^{+}_{2}$ | 27.234 |  | $3^{-}$ | 12.864 | 480.28 |
|  |  | $4^{-}$ | 15.253 | 219.26 |
|  |  | $5^{-}$ | 19.915 | 90.68 |
|  |  | $6^{-}$ | 25.131 | 30.99 |
$6^{-}$ | 25.131 |  | $4^{+}$ | 9.015 | 274.90 |

The energy spectra of 12C at large $\lambda_{3}$ can be understood by
comparing them to the spectra of a triangular rigid rotor discussed above.
Taking such a rotor with equal masses on the three vertices, we can compute
the moment of inertia $\mathcal{M}$ along every principle axis. For the
orientation shown in Fig. 1, the principle axes are simply the three
coordinate axes $x$, $y$ and $z$. It is straightforward to prove that
2$\mathcal{M}_{x}$ = 2$\mathcal{M}_{z}$ = $\mathcal{M}_{y}$ = $\mathcal{M}$.
Only considering rotations about the center of mass, we can write the
Hamiltonian as

$\hat{H}=\frac{\hat{I}_{x}^{2}}{\mathcal{M}}+\frac{\hat{I}_{y}^{2}}{2\mathcal{M}}+\frac{\hat{I}_{z}^{2}}{\mathcal{M}}=\frac{\hat{I}^{2}}{\mathcal{M}}-\frac{\hat{I}_{y}^{2}}{2\mathcal{M}},$
(30)

with $\hat{I}_{x}$, $\hat{I}_{y}$ and $\hat{I}_{z}$ the angular momentum
components in the intrinsic frame and $\hat{I}^{2}$ the total angular
momentum. As $\hat{H}$, $\hat{I}^{2}$ and $\hat{I}_{y}^{2}$ commute with each
other, we can take the total energy $E$, total angular momentum $I$ and the
angular momentum component $K$ as a group of good quantum numbers specifying
the eigenstates. Here $I(I+1)$ is the eigenvalue of $\hat{I}^{2}$ and $K$ is
the eigenvalue of $\hat{I}_{y}$. Rotation around the intrinsic $y$-axis by
$2\pi/3$ takes the nucleus back to its original position while produces a
phase $\exp(i2\pi K/3)$. This phase should be unit one for any physically
allowed state, thus we have $K=3n$ with $n$ an integer. For a state
$|IK\alpha\rangle$ with $\alpha$ quantum numbers other than $I$ and $K$, we
can calculate the energy according to Eq. (30),

$E=\langle
IK\alpha|\hat{H}|IK\alpha\rangle=\frac{1}{\mathcal{M}}[I(I+1)-K^{2}/2].$ (31)

The spatial reflection $P$ combined with a rotation around the 3-fold axis by
$180^{\circ}$ reset the system. The operator $P$ produce a factor of parity
$\pi$ and the rotation gives an additional phase of $(-1)^{K}$, thus we have
$\pi=(-1)^{K}$ for all physically allowed states.

The energy spectrum of Eq. (31) can be classified according to the quantum
number $K$. For $K=0$ we have the ground state band $0^{+}$, $2^{+}$, $4^{+}$,
…, with positive parity. For $K=3$ we have a negative parity band $3^{-}$,
$4^{-}$, $5^{-}$, … The excitation energies relative to the $0^{+}$ state are

$E_{x}(4^{+}):E_{x}(4^{-}):E_{x}(3^{-}):E_{x}(2^{+})=3.333:2.583:1.25:1$ (32)

up to an overall multiplicative factor. One important feature of this spectrum
is that the $4^{+}$ energy is higher than the $4^{-}$ energy. The reason is
that the $4^{+}$ state is a pure rotation around the 2-fold axes, while
$4^{-}$ state contains rotations around the 3-fold axis. The moment of inertia
for the latter rotation is relatively larger. However, in Fig. 5 we see that
the $4^{+}$ state is much lower than $3^{-}$ and $4^{-}$ for small
$\lambda_{3}$, only for extremely large $\lambda_{3}$ the projected energies
can be explained by the rigid rotor model. The reason is twofold. First, the
relation $\mathcal{M}_{y}$ = 2$\mathcal{M}_{z}$ used in deriving the rigid
rotor spectrum is only approximately satisfied for large deformations where
the masses of $\alpha$ clusters are concentrated at the vertices. For small
$\lambda_{3}$ we have a triangular density distribution without clustering and
the moments of inertia along different axis are not so different. Second, the
angular momentum and parity projections capture the full quantum mechanical
elements. The impact of the intrinsic states are reflected by the calculated
norm kernel and Hamiltonian kernel. Only for large $\lambda_{3}$ we can omit
the intrinsic structures and use the Gaussian overlap approximation (GOA) to
treat the $\alpha$ clusters as point particles.

In Table 2 and Table 3 we compare the low-lying spectrum and reduced
transition probabilities $B(E2)$ and $B(E3)$ of 12C with the experiments. The
results are calculated by projecting a mean field state constrained to
$\lambda_{3}=1.2$, which corresponds to the minimum of the $0^{+}$ PES in Fig
5. The experimental values are well reproduced at this level of
approximations. These results can be improved from various aspects. First, as
the PES’s are soft at small $\lambda_{3}$, in Fig. 5 we see that the
$\lambda_{3}$ value minimizing the projected energy has a strong dependence on
$I^{\pi}$, which should be considered in calculating the energy spectrum.
Second, the soft PES’s suggest we introduce the generating coordinate method
(GCM) to mix different deformations. Lastly, other configurations like the rod
shape might also play an important role in the low-lying excitations and
should be considered in building the low-lying states with good quantum
numbers.

In Ref. [52], the $\alpha$ arrangement in 12C is assumed to be an equilateral
triangle. In Ref. [118], the author performed calculations with the
orthogonality condition model (OCM) and the shallow potential model, and
obtained different configurations for each model. Due to insufficient
experimental data, whether the intrinsic configuration of the ground state of
12C is a regular triangle remains an open questions [118]. It is necessary to
explore this point by performing multi-dimensional constrained calculations
that include all important shape configurations. In Fig. 6 we show the
projected $0^{+}$ PES of 12C with respect to the deformation parameters
$\beta_{20}$ and $\lambda_{3}$. Here we constrain $\lambda_{3}$ from 0.0 to
3.0 with a step size $\Delta\lambda_{3}=0.2$, constrain $\beta_{20}$ from
$-0.2$ to $1.0$ with a step size $\Delta\beta_{20}=0.1$. The projected
energies are much lower than the mean field energies for non-zero
deformations. The rotational energy correction defined as the difference
between the $0^{+}$ energy and the mean field energy plays a key role in
determining the ground state configurations. On the $0^{+}$ PES, we found two
energy minima. The lowest one locates at $\lambda_{3}=0.8$ and
$\beta_{20}=0.6$, corresponding to a isosceles acute triangular shape, while
the shape isomer at $\lambda_{3}=1.2$ and $\beta_{20}=-0.1$ corresponds to an
isosceles obtuse triangular shape.

Figure 6: The 2-D potential energy surface of 12C projected to $0^{+}$ state.
The $\lambda_{3}$ and $\beta_{20}$ constraints are imposed with
$0<\lambda_{3}<3.0,~{}\Delta\lambda_{3}=0.2$ and
$-0.2<\beta_{20}<1.0,~{}\Delta\beta_{20}=0.1$.

##  4\. Summary

We have developed a projected multi-dimensionally constrained relativistic
Hartree-Bogoliubov (p-MDCRHB) model. In this model, the shape of the system is
only restricted as mirror reflection symmetric about $x$-$z$ and $y$-$z$
planes. In the mean field calculation, the RHB equation is solved in the
axially deformed harmonic oscillator (ADHO) basis. The parity and angular
momentum projection is performed to restore the intrinsic symmetries.

This model is well suited for exploring exotic shapes like triangle and
tetrahedron found in both the nuclear ground state and low-lying excited
states. One important scenario where these shapes play essential roles is the
nuclear clustering. The most important stellar nuclear reactions such as the
production of carbon and oxygen depend heavily on the clustering structures of
the reactant. In this work we benchmark the pMDCRHB model by studying the
clustering states in 12C. We define a triangular moment constraint and study
the corresponding triangular deformation $\lambda_{3}$. This constraint
generates a clustering state with three $\alpha$ clusters arranged on the
vertices of a regular triangle. We found that with the growth of
$\lambda_{3}$, the spectra approach that of a regular triangular rigid rotor.
This was explained by the symmetry analysis and the rigid rotor Hamiltonian.
We also performed the 2D constrained calculations of 12C. Two minima are found
on the projected $0^{+}$ PES. The softness of the PES’s indicates that the
configuration mixing GCM calculation is needed for a deeper understanding of
the structure of 12C and the related nuclear reactions occurring in the
stellar environment.

##  5\. Acknowledgments

Helpful discussions with Xiang-Xiang Sun and Yu-Ting Rong are gratefully
acknowledged. We thank Xiang-Xiang Sun and Shan-Gui Zhou for reading the
manuscript and valuable suggestions. This work has been supported by NSAF
(Grant No. U1930403), the National Key R&D Program of China (Grant No.
2018YFA0404402), the National Natural Science Foundation of China (Grants No.
11525524, No. 12070131001, No. 12047503, and No. 11961141004), the Key
Research Program of Frontier Sciences of Chinese Academy of Sciences (Grant
No. QYZDB-SSWSYS013) and the Strategic Priority Research Program of Chinese
Academy of Sciences (Grants No. XDB34010000 and No. XDPB15). The results
described in this paper are obtained on the High-performance Computing Cluster
of ITP-CAS and the ScGrid of the Supercomputing Center, Computer Network
Information Center of Chinese Academy of Sciences.

## Appendix: Calculation of the Rotation Matrix

Since the rotation matrix is naturally written as Wigner $D$-function with
spherical harmonic oscillators, the calculation of the rotation matrix can be
simplified by expanding the ADHO basis with the spherical harmonic oscillator
(SHO) basis.

An ADHO basis $\ket{\alpha}$ with quantum numbers
$\\{n_{z},n_{\rho},m_{l},m_{s}\\}$ reads

$\ket{\alpha}=\ket{n_{z}n_{\rho}m_{l}m_{s}}=C_{L}\phi_{n_{z}}(z)R_{m_{l}}^{n_{\rho}}(\rho,\varphi)\chi_{m_{s}},$
(33)

with

$\displaystyle\phi_{n_{z}}(z)$
$\displaystyle=\frac{1}{\sqrt{b_{z}}}\frac{1}{\pi^{1/4}\sqrt{2^{n_{z}}n_{z}!}}H_{n_{z}}\left(\frac{z}{b_{z}}\right){\rm
e}^{-\frac{z^{2}}{2b_{z}}},$ (34a) $\displaystyle
R_{m_{l}}^{n_{\rho}}(\rho,\varphi)$
$\displaystyle=\frac{1}{b_{\rho}}\sqrt{\frac{2n_{\rho}!}{(n_{\rho}+|m|)!}}\frac{1}{\sqrt{2\pi}}{\rm
e}^{-\rho^{2}/(2b_{\rho}^{2})}\left(\frac{\rho}{b_{\rho}}\right)^{|m|}L_{n_{\rho}}^{|m|}\left(\frac{\rho^{2}}{b_{\rho}^{2}}\right){\rm
e}^{im\varphi}.$ (34b)

$b_{z}$ and $b_{\rho}$ are the harmonic oscillator lengths along $z$ and
$\rho$ directions which are defined as $b_{z}=1/\sqrt{M\omega_{z}}$ and
$b_{\rho}=1/\sqrt{M\omega_{\rho}}$.

The SHO basis function reads

$\ket{n_{r}ljm}=\sum_{m_{s}}C_{m_{l}m_{s}jm}^{l\frac{1}{2}}\phi_{lm_{l}}^{n_{r}}(r,\Omega)\chi_{m_{s}},$
(35)

with $C_{m_{l}m_{s}jm}^{l\frac{1}{2}}$ the Clebsch-Gordon coefficient and

$\phi_{lm_{l}}^{n_{r}}(r,\Omega)=\frac{1}{b_{r}^{3/2}}\sqrt{\frac{2n_{r}!}{\Gamma(n_{r}+l+3/2)}}\left(\frac{r}{b_{r}}\right)^{l}\exp\left(-\frac{r^{2}}{2b_{r}^{2}}\right)L_{n_{r}}^{l+1/2}\left(\frac{r^{2}}{b_{r}^{2}}\right)Y_{lm}(\Omega).$
(36)

The transformation matrix can be written as

$\begin{split}\braket{n_{z}n_{\rho}m_{l}m_{s}}{n_{r}ljm_{j}}&=C_{L}^{*}C_{m_{l}m_{s}jm_{j}}^{l\frac{1}{2}}\delta_{m_{l}+m_{s},m_{j}}\braket{\phi_{n_{z}}R_{m_{l}}^{n_{\rho}}}{\phi_{lm_{l}}^{n_{r}}}\\\
&=C_{L}^{*}C_{m_{l}m_{s}jm_{j}}^{l\frac{1}{2}}\delta_{m_{l}+m_{s},m_{j}}\braket{n_{z}n_{\rho}m_{l}}{n_{r}lm_{l}}.\end{split}$
(37)

From Eq. (34) and Eq. (36), one can calculate the transformation matrix
element, which has also been given in Ref. [119]

$\braket{n_{r}lm}{n_{z}n_{\rho}m}=\sum_{\lambda}\frac{n_{\rho}!(2l+1)2^{n_{\rho}-n_{r}-\lambda}(-1)^{\lambda}(n_{r}+2\lambda)!}{\lambda!(n_{\rho}-\lambda)!(\lambda-
n_{\rho}+n_{r})!(2n_{r}+2\lambda-2n_{\rho}+2l+1)!!}.$ (38)

Thus, one obtains the relation between ADHO and SHO bases,

$\ket{\alpha}=\sum_{n_{r}ljm}C_{n_{r}ljm}^{\alpha}\ket{n_{r}ljm}.$ (39)

The time reversal states of the spherical harmonic oscillator basis are
obtained by applying the time reversal operator $\hat{T}$

$\hat{T}\ket{n_{r}ljm_{j}}=(-1)^{l+m_{j}+j}\ket{n_{r}lj\bar{m}_{j}}.$ (40)

The rotation matrix in the ADHO basis is calculated as

$\displaystyle\bra{\beta}\hat{R}(\Omega)\ket{\alpha}$
$\displaystyle=\sum_{n_{r}lj}\sum_{mm^{\prime}}C_{n_{r}ljm}^{\beta*}D_{mm^{\prime}}^{j}(\Omega)C_{n_{r}ljm^{\prime}}^{\alpha},$
(41a) $\displaystyle\bra{\bar{\beta}}\hat{R}(\Omega)\ket{\alpha}$
$\displaystyle=\sum_{n_{r}lj}\sum_{mm^{\prime}}(-1)^{l+m-j}C_{n_{r}ljm}^{\beta}D_{mm^{\prime}}^{j}(\Omega)C_{n_{r}ljm^{\prime}}^{\alpha}$
(41b)

By applying the properties of Wigner $D$-function, one finds the relation

$R_{\bar{\beta}\bar{\alpha}}=R_{\beta\alpha}^{*},\quad
R_{\beta\bar{\alpha}}=-R_{\bar{\beta}\alpha}^{*},\quad\mathrm{for}~{}\alpha,\beta>0,$
(42)

that means only half of $R$ matrix elements are needed in Eq. (41), while the
other half can be calculated directly from the relation (42).

## References

  * [1] Bohr A and Mottelson B 1969 _Nuclear Structure_ volume 1 Benjamin Inc., New York 1 edition.
  * [2] Bohr A and Mottelson B 1975 _Nuclear Structure_ volume 2 Benjamin Inc., New York 1 edition.
  * [3] Frauendorf S and Meng J 1997 _Nucl. Phys. A_ 617 131–147 http://dx.doi.org/10.1016/S0375-9474(97)00004-3.
  * [4] Ødegård S W _et al_ 2001 _Phys. Rev. Lett._ 86 5866–5869 http://dx.doi.org/10.1103/PhysRevLett.86.5866.
  * [5] Starosta K _et al_ 2001 _Phys. Rev. Lett._ 86 971–974 http://dx.doi.org/10.1103/PhysRevLett.86.971.
  * [6] Meng J _et al_ 2006 _Phys. Rev. C_ 73 037303 http://dx.doi.org/10.1103/PhysRevC.73.037303.
  * [7] Ayangeakaa A D _et al_ 2013 _Phys. Rev. Lett._ 110 172504 http://dx.doi.org/10.1103/PhysRevLett.110.172504.
  * [8] Lieder E O _et al_ 2014 _Phys. Rev. Lett._ 112 202502 http://dx.doi.org/10.1103/PhysRevLett.112.202502.
  * [9] Kuti I _et al_ 2014 _Phys. Rev. Lett._ 113 032501 http://dx.doi.org/10.1103/PhysRevLett.113.032501.
  * [10] Tonev D _et al_ 2014 _Phys. Rev. Lett._ 112 052501 http://dx.doi.org/10.1103/PhysRevLett.112.052501.
  * [11] Liu C _et al_ 2016 _Phys. Rev. Lett._ 116 112501 http://dx.doi.org/10.1103/PhysRevLett.116.112501.
  * [12] Bengtsson R _et al_ 1984 _Nucl. Phys. A_ 415 189–214 http://dx.doi.org/10.1016/0375-9474(84)90620-1.
  * [13] Liu Y Z _et al_ 1995 _Phys. Rev. C_ 52 2514–2519 http://dx.doi.org/10.1103/PhysRevC.52.2514.
  * [14] Liu Y Z _et al_ 1996 _Phys. Rev. C_ 54 719–730 http://dx.doi.org/10.1103/PhysRevC.54.719.
  * [15] Zhou S G _et al_ 1996 _J. Phys. G: Nucl. Part. Phys._ 22 415–420 http://dx.doi.org/10.1088/0954-3899/22/3/014.
  * [16] Riedinger L _et al_ 1997 _Prog Part Nucl Phys_ 38 251–261 http://dx.doi.org/10.1016/S0146-6410(97)00035-5 4 pi High Resolution Gamma Ray Spectroscopy and Nuclear Structure.
  * [17] Liu Y Z _et al_ 1998 _Phys. Rev. C_ 58 1849–1852 http://dx.doi.org/10.1103/PhysRevC.58.1849.
  * [18] Butler P A and Nazarewicz W 1996 _Rev. Mod. Phys._ 68 349–421 http://dx.doi.org/10.1103/RevModPhys.68.349.
  * [19] Möller P _et al_ 2008 _At. Data. Nucl. Data Tables_ 94 758–780 http://dx.doi.org/https://doi.org/10.1016/j.adt.2008.05.002.
  * [20] Agbemava S E, Afanasjev A V and Ring P 2016 _Phys. Rev. C_ 93 044304 http://dx.doi.org/10.1103/PhysRevC.93.044304.
  * [21] Shneidman T M _et al_ 2003 _Phys. Rev. C_ 67 014313 http://dx.doi.org/10.1103/PhysRevC.67.014313.
  * [22] Shneidman T M _et al_ 2006 _Phys. Rev. C_ 74 034316 http://dx.doi.org/10.1103/PhysRevC.74.034316.
  * [23] Wang S _et al_ 2005 _Phys. Rev. C_ 72 024317 http://dx.doi.org/10.1103/PhysRevC.72.024317.
  * [24] Yang D _et al_ 2009 _Chinese Phys. Lett._ 26 082101 http://dx.doi.org/10.1088/0256-307x/26/8/082101.
  * [25] Robledo L M and Bertsch G F 2011 _Phys. Rev. C_ 84 054302 http://dx.doi.org/10.1103/PhysRevC.84.054302.
  * [26] Zhu S J _et al_ 2012 _Phys. Rev. C_ 85 014330 http://dx.doi.org/10.1103/PhysRevC.85.014330.
  * [27] Nomura K, Vretenar D and Lu B N 2013 _Phys. Rev. C_ 88 021303 http://dx.doi.org/10.1103/PhysRevC.88.021303.
  * [28] Nomura K, Rodríguez-Guzmán R and Robledo L M 2015 _Phys. Rev. C_ 92 014312 http://dx.doi.org/10.1103/PhysRevC.92.014312.
  * [29] Gaffney L P _et al_ 2013 _Nature_ 497 199–204 http://dx.doi.org/10.1038/nature12073.
  * [30] Bucher B _et al_ 2016 _Phys. Rev. Lett._ 116 112503 http://dx.doi.org/10.1103/PhysRevLett.116.112503.
  * [31] Chishti M M R _et al_ 2020 _Nat. Phys._ 16 853–856 http://dx.doi.org/10.1038/s41567-020-0899-4.
  * [32] Lu B N _et al_ 2012 _EPJ Web Conf._ 38 05003 http://dx.doi.org/10.1051/epjconf/20123805003.
  * [33] Lu B N _et al_ 2014 _Phys. Rev. C_ 89 014323 http://dx.doi.org/10.1103/PhysRevC.89.014323.
  * [34] Lu B N _et al_ 2014 _Phys. Scr._ 89 054028 http://stacks.iop.org/1402-4896/89/i=5/a=054028.
  * [35] Zhao J _et al_ 2015 _Phys. Rev. C_ 91 014321 http://dx.doi.org/10.1103/PhysRevC.91.014321.
  * [36] Chai Q Z, Zhao W J and Wang H L 2019 _Comm. Theo. Phys._ 71 067 http://dx.doi.org/10.1088/0253-6102/71/1/67.
  * [37] Freer M 2007 _Rep. Progr. Phys._ 70 2149–2210 http://dx.doi.org/10.1088/0034-4885/70/12/r03.
  * [38] von Oertzen W, Freer M and Kanada-En’yo Y 2006 _Phys. Rep._ 432 43–113 http://dx.doi.org/10.1016/j.physrep.2006.07.001.
  * [39] Freer M _et al_ 2018 _Rev. Mod. Phys._ 90 035004 http://dx.doi.org/10.1103/RevModPhys.90.035004.
  * [40] Feng J _et al_ 2018 _Science China Physics, Mechanics & Astronomy_ 62 12011 https://doi.org/10.1007/s11433-018-9258-7.
  * [41] Arumugam P _et al_ 2005 _Phys. Rev. C_ 71 064308 http://dx.doi.org/10.1103/PhysRevC.71.064308.
  * [42] Maruhn J A _et al_ 2006 _Phys. Rev. C_ 74 044311 http://dx.doi.org/10.1103/PhysRevC.74.044311.
  * [43] Maruhn J _et al_ 2010 _Nucl. Phys. A_ 833 1–17 http://dx.doi.org/10.1016/j.nuclphysa.2009.12.005.
  * [44] Reinhard P G _et al_ 2011 _Phys. Rev. C_ 83 034312 http://dx.doi.org/10.1103/PhysRevC.83.034312.
  * [45] Ichikawa T _et al_ 2011 _Phys. Rev. Lett._ 107 112501 http://dx.doi.org/10.1103/PhysRevLett.107.112501.
  * [46] Ebran J P _et al_ 2012 _Nature_ 487 341–344 http://dx.doi.org/10.1038/nature11246.
  * [47] Ebran J P _et al_ 2014 _Phys. Rev. C_ 90 054329 http://dx.doi.org/10.1103/PhysRevC.90.054329.
  * [48] Zhao P W, Itagaki N and Meng J 2015 _Phys. Rev. Lett._ 115 022501 http://dx.doi.org/10.1103/PhysRevLett.115.022501.
  * [49] Zhou E _et al_ 2016 _Phys. Lett. B_ 753 227–231 http://dx.doi.org/https://doi.org/10.1016/j.physletb.2015.12.028.
  * [50] Ren Z _et al_ 2019 _Sci. China Phys. Mech._ 62 112062 https://doi.org/10.1007/s11433-019-9412-3.
  * [51] Ring P 2019 _Sci. China Phys. Mech._ 62 112063 https://doi.org/10.1007/s11433-019-9440-9.
  * [52] Bijker R and Iachello F 2002 _Ann. Phys._ 298 334–360 http://dx.doi.org/10.1006/aphy.2002.6255.
  * [53] Kanada-En’yo Y 2007 _Prog. Theo. Phys._ 117 655–680 http://dx.doi.org/10.1143/PTP.117.655.
  * [54] Marín-Lámbarri D J _et al_ 2014 _Phys. Rev. Lett._ 113 012502 http://dx.doi.org/10.1103/PhysRevLett.113.012502.
  * [55] Bijker R and Iachello F 2014 _Phys. Rev. Lett._ 112 152501 http://dx.doi.org/10.1103/PhysRevLett.112.152501.
  * [56] Bijker R and Iachello F 2021 _Nucl. Phys. A_ 1006 122077 http://dx.doi.org/10.1016/j.nuclphysa.2020.122077.
  * [57] Reinhard P G 1989 _Rep. Prog. Phys._ 52 439 http://stacks.iop.org/0034-4885/52/i=4/a=002.
  * [58] Ring P 1996 _Prog. Part. Nucl. Phys._ 37 193–263 http://dx.doi.org/10.1016/0146-6410(96)00054-3.
  * [59] Vretenar D _et al_ 2005 _Phys. Rep._ 409 101–259 http://dx.doi.org/10.1016/j.physrep.2004.10.001.
  * [60] Meng J _et al_ 2006 _Prog. Part. Nucl. Phys._ 57 470–563 http://dx.doi.org/10.1016/j.ppnp.2005.06.001.
  * [61] Meng J _et al_ 2011 _Prog. Phys._ 31 199–336 http://pip.nju.edu.cn/Home/ShowArticle/743.
  * [62] Nikšić T, Vretenar D and Ring P 2011 _Prog. Part. Nucl. Phys._ 66 519–548 http://dx.doi.org/10.1016/j.ppnp.2011.01.055.
  * [63] Meng J and Zhou S G 2015 _J. Phys. G: Nucl. Phys._ 42 093101 http://stacks.iop.org/0954-3899/42/i=9/a=093101.
  * [64] Liang H Z, Meng J and Zhou S G 2015 _Phys. Rep._ 570 1–84 http://dx.doi.org/10.1016/j.physrep.2014.12.005.
  * [65] Meng J, ed. 2016 _Relativistic Density Functional for Nuclear Structure_ World Scientific, Singapore http://dx.doi.org/10.1142/9872.
  * [66] Ginocchio J N 1997 _Phys. Rev. Lett._ 78 436–439 http://dx.doi.org/10.1103/PhysRevLett.78.436.
  * [67] Sugawara-Tanabe K and Arima A 1998 _Phys. Rev. C_ 58 R3065–R3068 http://dx.doi.org/10.1103/PhysRevC.58.R3065.
  * [68] Meng J _et al_ 1998 _Phys. Rev. C_ 58 R628–R631 http://dx.doi.org/10.1103/PhysRevC.58.R628.
  * [69] Meng J _et al_ 1999 _Phys. Rev. C_ 59 154–163 http://dx.doi.org/10.1103/PhysRevC.59.154.
  * [70] Lu B N, Zhao E G and Zhou S G 2012 _Phys. Rev. Lett._ 109 072501 http://dx.doi.org/10.1103/PhysRevLett.109.072501.
  * [71] Lu B N, Zhao E G and Zhou S G 2013 _Phys. Rev. C_ 88 024323 http://dx.doi.org/10.1103/PhysRevC.88.024323.
  * [72] Zhou S G, Meng J and Ring P 2003 _Phys. Rev. Lett._ 91 262501 http://dx.doi.org/10.1103/PhysRevLett.91.262501.
  * [73] He X T _et al_ 2006 _Eur. Phys. J. A_ 28 265–269 http://dx.doi.org/10.1140/epja/i2006-10066-0.
  * [74] Ring P and Schuck P 1980 _The Nuclear Many-Body Problem_ Springer-Verlag New York Inc., New York. http://dx.doi.org/10.1063/1.2915762.
  * [75] Yao J M _et al_ 2010 _Phys. Rev. C_ 81 044311 http://dx.doi.org/10.1103/PhysRevC.81.044311.
  * [76] Yao J M _et al_ 2011 _Phys. Rev. C_ 83 014308 http://dx.doi.org/10.1103/PhysRevC.83.014308.
  * [77] Sun X X and Zhou S G 2021 _Phys. Rev. C_ 104 064319 https://doi.org/10.1103/physrevc.104.064319.
  * [78] Sun X X and Zhou S G 2021 _Sci. Bull._ 66 2072–2078 http://dx.doi.org/https://doi.org/10.1016/j.scib.2021.07.005.
  * [79] Lu B N, Zhao E G and Zhou S G 2012 _Phys. Rev. C_ 85 011301 http://dx.doi.org/10.1103/PhysRevC.85.011301.
  * [80] Zhao J _et al_ 2017 _Phys. Rev. C_ 95 014320 http://dx.doi.org/10.1103/PhysRevC.95.014320.
  * [81] Zhou S G 2016 _Phys. Scr._ 91 063008 http://stacks.iop.org/1402-4896/91/i=6/a=063008.
  * [82] Meng X, Lu B N and Zhou S G 2019 _Sci. China Phys. Mech._ 63 212011 http://dx.doi.org/10.1007/s11433-019-9422-1.
  * [83] Wang X Q, Sun X X and Zhou S G 2021 _Chinese Phys. C_ in press http://dx.doi.org/10.1088/1674-1137/ac3904.
  * [84] Zhao J _et al_ 2012 _Phys. Rev. C_ 86 057304 http://dx.doi.org/10.1103/PhysRevC.86.057304.
  * [85] Zhao J _et al_ 2015 _Phys. Rev. C_ 92 064315 http://dx.doi.org/10.1103/PhysRevC.92.064315.
  * [86] Zhao J _et al_ 2016 _Phys. Rev. C_ 93 044315 http://dx.doi.org/10.1103/PhysRevC.93.044315.
  * [87] Zhao J _et al_ 2019 _Phys. Rev. C_ 99 014618 http://dx.doi.org/10.1103/PhysRevC.99.014618.
  * [88] Zhao J _et al_ 2019 _Phys. Rev. C_ 99 054613 http://dx.doi.org/10.1103/PhysRevC.99.054613.
  * [89] Zhao J _et al_ 2020 _Phys. Rev. C_ 101 064605 http://dx.doi.org/10.1103/PhysRevC.101.064605.
  * [90] Lu B N, Zhao E G and Zhou S G 2011 _Phys. Rev. C_ 84 014328 http://dx.doi.org/10.1103/PhysRevC.84.014328.
  * [91] Lu B N _et al_ 2014 _Phys. Rev. C_ 89 044307 http://dx.doi.org/10.1103/PhysRevC.89.044307.
  * [92] Rong Y T, Zhao P W and Zhou S G 2020 _Phys. Lett. B_ 807 135533 http://dx.doi.org/10.1016/j.physletb.2020.135533.
  * [93] Rong Y T, Tu Z H and Zhou S G 2021 _Phys. Rev. C_ 104 054321 http://dx.doi.org/10.1103/PhysRevC.104.054321.
  * [94] Yao J M _et al_ 2009 _Phys. Rev. C_ 79 044312 http://dx.doi.org/10.1103/PhysRevC.79.044312.
  * [95] Yao J M, Zhou E F and Li Z P 2015 _Phys. Rev. C_ 92 041304 http://dx.doi.org/10.1103/PhysRevC.92.041304.
  * [96] Xia H _et al_ 2018 _Sci. China Phys. Mech._ 62 42011 https://doi.org/10.1007/s11433-018-9308-0.
  * [97] Nikolaus B A, Hoch T and Madland D G 1992 _Phys. Rev. C_ 46 1757–1781 http://dx.doi.org/10.1103/PhysRevC.46.1757.
  * [98] Bürvenich T _et al_ 2002 _Phys. Rev. C_ 65 044308 http://dx.doi.org/10.1103/PhysRevC.65.044308.
  * [99] Boguta J and Bodmer A 1977 _Nucl. Phys. A_ 292 413–428 http://dx.doi.org/10.1016/0375-9474(77)90626-1.
  * [100] Brockmann R and Toki H 1992 _Phys. Rev. Lett._ 68 3408–3411 http://dx.doi.org/10.1103/PhysRevLett.68.3408.
  * [101] Sugahara Y and Toki H 1994 _Nucl. Phys. A_ 579 557–572 http://dx.doi.org/10.1016/0375-9474(94)90923-7.
  * [102] Fuchs C, Lenske H and Wolter H H 1995 _Phys. Rev. C_ 52 3043–3060 http://dx.doi.org/10.1103/PhysRevC.52.3043.
  * [103] Nikšić T _et al_ 2002 _Phys. Rev. C_ 66 024306 http://dx.doi.org/10.1103/PhysRevC.66.024306.
  * [104] Nikšić T, Vretenar D and Ring P 2008 _Phys. Rev. C_ 78 034318 http://dx.doi.org/10.1103/PhysRevC.78.034318.
  * [105] Kucharek H and Ring P 1991 _Z. Phys. A_ 339 23–35 http://dx.doi.org/10.1007/BF01282930.
  * [106] Serra M and Ring P 2002 _Phys. Rev. C_ 65 064324 http://dx.doi.org/10.1103/PhysRevC.65.064324.
  * [107] Tian Y, Ma Z and Ring P 2009 _Phys. Lett. B_ 676 44–50 http://dx.doi.org/https://doi.org/10.1016/j.physletb.2009.04.067.
  * [108] Tian Y, Ma Z y and Ring P 2009 _Phys. Rev. C_ 79 064301 http://dx.doi.org/10.1103/PhysRevC.79.064301.
  * [109] Tian Y, Ma Z y and Ring P 2009 _Phys. Rev. C_ 80 024313 http://dx.doi.org/10.1103/PhysRevC.80.024313.
  * [110] Ring P, Gambhir Y and Lalazissis G 1997 _Comput. Phys. Commun._ 105 77–97 http://dx.doi.org/10.1016/S0010-4655(97)00022-2.
  * [111] Hoyle F 1954 _Astrophys. J. Suppl. Ser._ 1 121 http://dx.doi.org/10.1086/190005.
  * [112] Freer M and Fynbo H 2014 _Prog. Part. Nucl. Phys._ 78 1–23 http://dx.doi.org/10.1016/j.ppnp.2014.06.001.
  * [113] Bijker R 2016 _Phys. Scr._ 91 073005 http://stacks.iop.org/1402-4896/91/i=7/a=073005.
  * [114] Hamermesh M 1962 _Group Theory and Its Application to Physical Problems_ Dover Publications, Inc., New York.
  * [115] Góźdź A and Dudek J 1996 _Acta Phys. Pol. B_ 27 469–474 https://www.actaphys.uj.edu.pl/R/27/1/469.
  * [116] Tagami S, Shimizu Y R and Dudek J 2013 _Phys. Rev. C_ 87 054306 http://dx.doi.org/10.1103/PhysRevC.87.054306.
  * [117] Ajzenberg-Selove F 1990 _Nucl. Phys. A_ 506 1–158 http://dx.doi.org/10.1016/0375-9474(90)90271-M.
  * [118] Moriya H _et al_ 2021 _Few-Body Systems_ 62 46 http://dx.doi.org/10.1007/s00601-021-01631-2.
  * [119] Talman J D 1970 _Nucl. Phys. A_ 141 273–288 http://dx.doi.org/10.1016/0375-9474(70)90847-X.


# Gaia Astrometry and MIKE+PFS Doppler Data Joint Analysis Reveals that HD
175167b is a Massive Cold Jupiter

Tianjun Gan Department of Astronomy, Tsinghua University, Beijing 100084,
People’s Republic of China

###### Abstract

HD 175167b is a cold ($P_{b}\sim 1200$ days) Jupiter with a minimum mass of
$M_{p}\sin i=7.8\pm 3.5\ M_{J}$ orbiting a Sun-like star, first discovered by
the Magellan Planet Search Program based on MIKE observations. Through a joint
analysis of the MIKE data and the Gaia two-body orbital solution, Winn (2022)
found a companion mass of $M_{p}=14.8\pm 1.8\ M_{J}$ and suggested that it
might be better designated as a brown dwarf. Additional publicly available
radial velocity data from Magellan/PFS better constrains the model, and
reveals that the companion is a massive cold Jupiter with a mass of
$M_{p}=10.2\pm 0.4\ M_{J}$ and a period of $P_{b}=1275.8\pm 0.4$ days. The
planet orbit is inclined by $i=38.6\pm 1.7^{\circ}$ with an eccentricity of
$0.529\pm 0.002$.

Giant planets; Radial Velocity; Astrometry

††facilities: Gaia, Magellan/MIKE, Magellan/PFS††software: radvel (Fulton et
al., 2018), emcee (Foreman-Mackey et al., 2013), Rebound (Rein & Liu, 2012)

##  1 Introduction

HD 175167b is a long-period massive gas giant around a G5V star reported by
Arriagada et al. (2010) using the MIKE radial velocity (RV) data taken under
the Magellan Planet Search Program (Bernstein et al., 2003). The host star has
a mass about $1\ M_{\odot}$ while the planet has a period of $1290\pm 22$ days
with a minimum mass of $M_{p}\sin i=7.8\pm 3.5\ M_{J}$ on an eccentric
($e=0.54\pm 0.09$) orbit. The minimum mass of the companion was later refined
to $9.0\pm 3.3\ M_{J}$ by Stassun et al. (2017) based on a better-constrained
empirical stellar mass. More recently, Feng et al. (2022) and Xiao et al.
(2023) both carried out full orbit fits with a combination of ground-based RV
measurements and astrometric data from Hipparcos (Perryman et al., 1997; van
Leeuwen, 2007) as well as Gaia (Gaia Collaboration et al., 2016, 2018, 2021),
enabling them to determine the true mass of the secondary. However, since the
time-series Gaia astrometric data have not been published and only two
absolute astrometry measurements are available now at the epoch of Year 1991
and 2016, the derived companion mass ($M_{p}=6.4\pm 1.4\ M_{J}$ and $9.8\pm
1.9\ M_{J}$) and other orbital parameters still have large uncertainties.

As part of Gaia Data Release 3 (Gaia DR3; Gaia Collaboration et al., 2023),
Halbwachs et al. (2023) presented a “two-body” catalog (nss_two_body_orbit)
based on Gaia observations spanning about 33 months. It includes orbital
solutions for about 165,500 systems, among which 1162 sources are potential
substellar objects (Holl et al., 2023). Each orbital solution contains the
best-fit orbital parameters like four Thiele-Innes coefficients (Binnendijk,
1960; Heintz, 1978), period, eccentricity, time of periastron and parallax
along with a correlation matrix between them. A recent work from Winn (2022)
proposed a method to conduct a joint analysis of these two-body solutions and
RV data, which has been used in many relevant studies (e.g., Unger et al.,
2023; Fitzmaurice et al., 2023). Basically, the Gaia two-body solution is
treated as a basis and the correlation matrix is transformed to a covariance
matrix using the nsstools111https://www.cosmos.esa.int/web/gaia/dr3-nss-tools
algorithm (Halbwachs et al., 2023) to construct a Gaia likelihood function.
MCMC sampling is then performed to maximize the combined Doppler$+$Gaia
likelihood function. Note that this approach is only valid when the system
contains a single companion. For multi-companion systems, the Gaia two-body
solution from a single Keplerian model fit would not be robust thus the joint
fit cannot be done.

With this method, Winn (2022) performed a joint analysis for the HD 175167
system using the published MIKE data and the Gaia two-body solution, leading
to a companion mass of $14.8\pm 1.8\ M_{J}$, a period of $1175\pm 25$ days and
an orbital inclination of $35.5\pm 2.3^{\circ}$. Therefore, Winn (2022)
suggested that the companion might be located in the brown dwarf mass range
(Burrows et al., 1997; Spiegel et al., 2011). However, two caveats were also
pointed out by Winn (2022): (1) the Gaia orbital solution has large
uncertainties, which may be due to the shorter time span of Gaia observations
compared with the companion’s period; (2) the periods and eccentricities are
discrepant between the Gaia and Doppler solutions, which is perhaps caused by
non-Gaussian uncertainties of Gaia orbital parameters.

##  2 Magellan/PFS data and Joint Analysis

In addition to the MIKE radial velocity data, HD 175167 was also monitored by
the Planet Finder Spectrograph (PFS; Crane et al., 2006, 2008, 2010) on the
6.5 m Magellan II (Clay) telescope between 2010 May 23 and 2019 August 14 (see
the data published in Feng et al. 2022). A total of 22 observations were
collected with a median RV uncertainty of 1.1 m/s.

Following the methodology proposed by Winn (2022), I perform a joint-fit
analysis including the PFS data. The joint model has 11 free parameters
including stellar mass $M$, companion mass $m$, $e\cos\omega$, $e\sin w$, the
cosine of the inclination $\cos i$, the longitude of the ascending node
$\Omega$, orbital period $P$, time of periastron $t_{p}$, parallax $\varpi$
and two RV offsets $\gamma_{\rm MIKE}$, $\gamma_{\rm PFS}$, where $e$ is the
orbital eccentricity and $\omega$ represents the argument of pericenter. The
Keplerian model is built by utilizing the code radvel (Fulton et al., 2018).
Since the companion is a dark body based on its mass measurement from previous
works, the flux ratio $\varepsilon$ between the companion and the primary star
is thus fixed at 0. Except for a Gaussian prior $\mathcal{N}$ ($1.00$ ,
$0.04^{2}$) adopted on the host star mass based on the value from literature
(Feng et al., 2022), all the other parameters are set with non-informative
uniform priors. The posterior sampling is carried out with emcee (Foreman-
Mackey et al., 2013). A total of 100 walkers are initialized and all of them
are run for 35,000 steps with the first 5,000 steps excluded as burnt-in
samples. Figure 1 presents the final joint-fit results, the orbital geometry
made by Rebound (Rein & Liu, 2012) as well as the posterior distributions.

##  3 Results and Discussions

After including additional RV data points, the uncertainty on the companion
mass is successfully reduced. The joint-fit of Doppler velocities from
Magellan/MIKE and Magellan/PFS along with the Gaia two-body solution reveals
that the companion has a mass of $10.2\pm 0.4\ M_{J}$, an eccentricity of
$e=0.529\pm 0.002$ and a period of $1275.8\pm 0.4$ days. The planet orbital
period constrained by the new joint-fit here is better in accord with the
values $1290.0\pm 22.0$, $1276.1\pm 0.7$ and $1289^{+8}_{-28}$ days from
Arriagada et al. (2010), Feng et al. (2022) and Xiao et al. (2023), compared
with $1175\pm 25$ days determined by Winn (2022). The best-fit planet orbit
inclination is $38.6\pm 1.7^{\circ}$, which is consistent with the result
$35.5\pm 2.3^{\circ}$ reported by Winn (2022) but different from $94.6\pm
23.7^{\circ}$ and $60\pm 17^{\circ}$ measured by Feng et al. (2022) and Xiao
et al. (2023). One possibility is that the small number (2 data points from
Hipparcos and Gaia) of absolute astrometry measurements limits the accuracy
and precision of the orbital inclination measurement in Feng et al. (2022) and
Xiao et al. (2023). Based on the updated mass measurement from the joint
analysis, HD 175167b turns out to be a massive cold gas giant planet rather
than a brown dwarf.

Figure 1: Top panels: Time-series (left) and phase-folded (right) RV data from
MIKE and PFS along with the best-fit model. Middle right panel: The orbit of
the HD 175167b constructed based on the best-fit physical parameters using
Rebound (Rein & Liu, 2012). The parent star is located at the origin. Bottom
panel: The posterior distributions of the joint-fit analysis. We follow the
Gaia convention where $t_{p}$ is the Julian date of pericenter minus
2,457,389.

##  4 Acknowledgments

The author thanks Drs. Shude Mao, Sharon X. Wang and Joshua Winn for the
useful discussions. This work is supported by the National Science Foundation
of China (Grant No. 12133005).

## References

  * Arriagada et al. (2010) Arriagada, P., Butler, R. P., Minniti, D., et al. 2010, ApJ, 711, 1229, doi: 10.1088/0004-637X/711/2/1229
  * Bernstein et al. (2003) Bernstein, R., Shectman, S. A., Gunnels, S. M., Mochnacki, S., & Athey, A. E. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1694–1704, doi: 10.1117/12.461502
  * Binnendijk (1960) Binnendijk, L. 1960, Properties of double stars; a survey of parallaxes and orbits.
  * Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856, doi: 10.1086/305002
  * Crane et al. (2006) Crane, J. D., Shectman, S. A., & Butler, R. P. 2006, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6269, The Carnegie Planet Finder Spectrograph, 626931, doi: 10.1117/12.672339
  * Crane et al. (2010) Crane, J. D., Shectman, S. A., Butler, R. P., et al. 2010, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, The Carnegie Planet Finder Spectrograph: integration and commissioning, 773553, doi: 10.1117/12.857792
  * Crane et al. (2008) Crane, J. D., Shectman, S. A., Butler, R. P., Thompson, I. B., & Burley, G. S. 2008, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, The Carnegie Planet Finder Spectrograph: a status report, 701479, doi: 10.1117/12.789637
  * Feng et al. (2022) Feng, F., Butler, R. P., Vogt, S. S., et al. 2022, ApJS, 262, 21, doi: 10.3847/1538-4365/ac7e57
  * Fitzmaurice et al. (2023) Fitzmaurice, E., Stefánsson, G., Kavanagh, R. D., et al. 2023, arXiv e-prints, arXiv:2310.07827, doi: 10.48550/arXiv.2310.07827
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504, doi: 10.1088/1538-3873/aaaaa8
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2, doi: 10.1051/0004-6361/201629512
  * Gaia Collaboration et al. (2018) —. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2021) —. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Halbwachs et al. (2023) Halbwachs, J.-L., Pourbaix, D., Arenou, F., et al. 2023, A&A, 674, A9, doi: 10.1051/0004-6361/202243969
  * Heintz (1978) Heintz, W. D. 1978, Double stars, Vol. 15
  * Holl et al. (2023) Holl, B., Sozzetti, A., Sahlmann, J., et al. 2023, A&A, 674, A10, doi: 10.1051/0004-6361/202244161
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Spiegel et al. (2011) Spiegel, D. S., Burrows, A., & Milsom, J. A. 2011, ApJ, 727, 57, doi: 10.1088/0004-637X/727/1/57
  * Stassun et al. (2017) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136, doi: 10.3847/1538-3881/aa5df3
  * Unger et al. (2023) Unger, N., Ségransan, D., Barbato, D., et al. 2023, arXiv e-prints, arXiv:2310.02758, doi: 10.48550/arXiv.2310.02758
  * van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653, doi: 10.1051/0004-6361:20078357
  * Winn (2022) Winn, J. N. 2022, AJ, 164, 196, doi: 10.3847/1538-3881/ac9126
  * Xiao et al. (2023) Xiao, G.-Y., Liu, Y.-J., Teng, H.-Y., et al. 2023, Research in Astronomy and Astrophysics, 23, 055022, doi: 10.1088/1674-4527/accb7e