#  Proximity effects of vortices in neutron ${}^{3}P_{2}$ superfluids in
neutron stars:
Vortex core transitions and covalent bonding of vortex molecules

Michikazu Kobayashi  School of Environmental Science and Engineering, Kochi
University of Technology, Miyanoguchi 185, Tosayamada, Kami, Kochi 782-8502,
Japan  Muneto Nitta  Department of Physics, and Research and Education Center
for Natural Sciences, Keio University, Hiyoshi 4-1-1, Yokohama, Kanagawa
223-8521, Japan

###### Abstract

Neutron ${}^{3}P_{2}$ superfluids consisting of neutron pairs with the total
angular momentum $J=2$ with spin-triplet and $P$-wave are believed to be
realized in neutron star cores. Within the Ginzburg-Landau theory it was
previously found that a singly quantized vortex is split into two half-
quantized non-Abelian vortices connected by one (or three) soliton(s) forming
a vortex molecule with the soliton bond(s), in the absence (presence) of
magnetic field parallel to them. In this paper, we investigate proximity
effects of two vortex molecules by exhausting all possible two vortex molecule
states consisting of four half-quantized vortices and determining the phase
diagram spanned by the magnetic field and rotation speed. As the rotation
speed is increased, the distance between the two vortex molecules becomes
shorter. In the magnetic field below the critical value, we find that as the
rotation speed is increased, the two separated vortex molecules transit to a
dimerized vortex molecule, where the two vortex molecules are bridged by two
solitons that we call “covalent bonds” in analogy with chemical molecules. We
also find that the orders of the constituent half-quantized vortex cores
transit from a ferromagnetic order to a cyclic order as the vortex molecules
come closer. On the other hand, no dimerization occurs in the magnetic field
above the critical value. Instead, we find a transition for the polarization
direction of the vortex molecules from a configuration parallel to the
separation to one perpendicular to the separation as they come closer. We also
show some examples of three and four vortex molecule states.

##  I Introduction

Recently, there have been great progresses in astrophysical observations of
neutron stars (pulsars), such as massive neutron stars with masses about twice
as large as the solar mass [1, 2], detection of gravitational waves from a
binary neutron star merger [3, 4], and the Neutron star Interior Composition
Explorer (NICER) mission [5, 6]. These are providing us astrophysical
laboratories for exploring nuclear and QCD matter under extreme condition:
extremely high density, with strong magnetic fields and under rapid rotations
[7, 8].

The interior of neutron stars exhibit neutron superfluidity and proton
superconductivity [9], see Refs. [10, 11, 12, 13, 7, 14] as a review. Such
supercomponents provide low-energy excitations affecting several processes and
properties of neutron stars: neutrino emissivities and specific heats relevant
to the long relaxation time after pulsar glitches (sudden speed-up events of
neutron stars) [15, 16, 17], and the enhancement of neutrino emission around
the critical point of the superfluid transition [18, 19, 20, 21, 22, 23]. When
superfluids are rotating, there appear a large number of quantum vortices
along the rotation axis, and they form a vortex lattice. The circulation
around each quantum vortex is quantized due to the Feynman-Onsager’s
quantization in order for the wave function to be single-valued. In neutron
stars, there should exist typically $10^{17}$ quantum vortices, considered to
play significant roles of neutron star dynamics. For instance, the origin of
pulsar glitches was suggested to be explained by sudden unpinning of a large
number of quantum vortices [24, 25].

At lower density corresponding to outer cores of neutron stars, neutron
superfluids are realized by the attraction by the ${}^{1}S_{0}$ channel
between two neutrons [26]. On the other hand, at higher density corresponding
to the inner cores of neutron stars, the ${}^{1}S_{0}$ channel is repulsive
due to the strong short-range repulsion. Instead, the ${}^{3}P_{2}$ channel
originating from a strong spin-orbit force at large scattering energy becomes
more dominant, where neutron Cooper pairs possess spin-triplet and $P$-wave
with the total angular momentum $J=2$ [27, 28, 29, 28, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44]. Furthermore, the ${}^{3}P_{2}$ channel is
tolerant against the strong magnetic field such as $10^{15}-10^{18}$ G for
magnetars, because aligned Cooper pairs with the spin-triplet pairing are not
broken by the Zeeman effect, in contrast to the $S$-wave Cooper pairs which
can survive at most around the magnetic field $10^{17}$ G [45]. In
astrophysical observations, the rapid cooling of the neutron star in
Cassiopeia A was proposed to be explained by the enhancement of neutrino
emissivities due to the formation and dissociation of neutron ${}^{3}P_{2}$
Cooper pairs [46, 22, 23].

Theoretically, there are two frameworks to deal with the ${}^{3}P_{2}$
superfluids: a microscopic theory known as the Bogoliubov-de Gennes (BdG)
equation describing fermion degrees of freedom, and the Ginzburg-Landau (GL)
theory for the order parameters conveniently describing bosonic excitations.
The latter can be obtained from the former by integrating out fermion degrees
of freedom as an expansion of the order parameters and spatial derivatives.
Thus, the GL theory is the low-energy effective theory describing large
distance behaviors, which is valid only in the region close to the critical
temperature. Among superfluid states with $J=2$ classified into nematic,
cyclic, and ferromagnetic phases etc [47], the GL theory for ${}^{3}P_{2}$
superfluids [32, 33, 48, 49, 50, 51, 52, 53, 54, 55, 56] predicts that the
ground state is in the nematic phase at least in the weak coupling limit [48,
49, 50]. The nematic phase consists of almost degenerate three different
states with different unbroken symmetries: uniaxial nematic (UN), $D_{2}$
biaxial nematic ($D_{2}$BN), and $D_{4}$ biaxial nematic ($D_{4}$BN) phases
with unbroken groups $O(2)$, $D_{2}$ and $D_{4}$, respectively. Here, $D_{n}$
is a dihedral group of order $n$ (see Table 1(a) and (b)). Nematic phases were
also known in spin-2 Bose-Einstein condensates (BECs) [57, 58, 59, 60, 61,
62]. Depending on the magnetic field and temperature, UN, $D_{2}$BN or
$D_{4}$BN state is realized as the ground state for zero magnetic field,
nonzero one below the critical value $B_{c}$, or nonzero one above $B_{c}$,
respectively [52, 54, 55, 56] [see Table 1(a)]. Apart from nematic phases,
more general uniform states (which do not have to be realized as the ground
states) of ${}^{3}P_{2}$ superfluids were classified according to symmetries
[62]. Beyond the quasi-classical approximation, the ferromagnetic phase was
found in the region close to the critical temperature [63]. The GL approach is
useful not only to determine the ground states but also to deal with bosonic
collective excitations [64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76]
relevant for cooling process of neutron stars, and various topological
excitations and defects, such as vortices (as explained below in more detail),
domain walls [77] and the boundary defect (boojums) [78].

On the other hand, the BdG approach offers a microscopic description with
fermion degrees of freedom valid at short distances and all ranges of
temperatures including zero temperature [79, 80, 63, 81, 82]. It was applied
to the phase diagram of ${}^{3}P_{2}$ superfluids in the plane of the
temperature and magnetic field, which is valid even at zero temperature [79],
including a tri-critical point connecting first and second order phase
transition lines between $D_{4}$ and $D_{2}$BN phases [79, 80]. Furthermore,
${}^{3}P_{2}$ superfluids were shown to be topological superfluids of a class
DIII in the classification of topological insulators and superconductors [83,
84], ensuring a topologically protected gapless Majorana fermion on its
boundary [79] and inside vortex cores [81, 82] as explained below.

Since it is promising that quantum vortices in $S$-wave superfluids play
significant roles in neutron star dynamics, the same should be expected for
the ${}^{3}P_{2}$ superfluids as well. In fact, quantum vortices were
investigated in the case of ${}^{3}P_{2}$ superfluids both in the GL theory
[49, 50, 32, 52, 85, 53, 86] (coreless vortices [87]) and in the BdG theory
[81, 82]. The first homotopy group classifies type of vortices in each phase
[52] as in Table 1(d). Singly quantized vortices in ${}^{3}P_{2}$ superfluids
were studied in the GL theory [49, 50, 32, 52, 85] and in the BdG theory [81]
with topologically protected Majorana fermions in the vortex cores. Vortices
more peculiar to the ${}^{3}P_{2}$ superfluids are half-quantized non-Abelian
vortices [53, 82, 86], which have a half of the Feynman-Onsager’s quantized
circulations and are characterized by a non-Abelian first homotopy group
$D_{4}^{*}$, thus giving non-commutativity when exchanging two vortices.
Isolated half-quantized non-Abelian vortices are topologically allowed only in
the $D_{4}$BN phase. The existence of half-quantized vortices was proposed to
explain a scaling law of pulsar glitches [88].

| $|\boldsymbol{B}|=0$ | $0<|\boldsymbol{B}|<B_{c}$ | $B_{c}<|\boldsymbol{B}|$
---|---|---|---
(a) Phase | UN |  $D_{2}$BN |  $D_{4}$BN
(b) Symmetry | $O(2)$ | $D_{2}$ | $D_{4}$
(c) OPM | $S^{1}\times{\mathbb{R}}P^{2}$ | $U(1)\times\frac{SO(3)}{D_{2}}$ | $\frac{U(1)\times SO(3)}{D_{4}}$
(d) $\pi_{1}$(OPM) | ${\mathbb{Z}}\oplus{\mathbb{Z}}_{2}$ | ${\mathbb{Z}}\oplus{\mathbb{Q}}$ | ${\mathbb{Z}}\rtimes_{h}D_{4}^{*}$
(e) Vortex core order | Ferro | Cyclic | Cyclic
(Cyclic)
(f) $\\#$ of solitons | 1 (3) | 3 | 3
(g) Soliton core order |  $D_{4}$BN |  $D_{4}$BN |  $D_{2}$BN
Table 1:  In the absence and presence of the magnetic field below and above
the critical magnetic field $B_{c}$, (a) bulk phases, (b) unbroken symmetries,
(c) OPMs, (d) the first homotopy groups of the OPMs, (e) the order of half-
quantized vortex cores, (f) the number of solitons connecting two half-
quantized vortices, and (g) the order of soliton cores are summarized. Here,
${\mathbb{Q}}=D_{2}^{*}$ is a quaternion group with $*$ implying the universal
covering group. $\rtimes_{h}$ is a product defined in Ref. [89], ensuring the
existence of isolated half-quantized non-Abelian vortices. The blackets in the
rows (e) and (f) and the column $|\boldsymbol{B}|=0$ represent for a
metastable solution. The most parts of this table are taken from our previous
paper [86].

While an axisymmetric ansatz was employed in the previous studies of vortex
solutions [49, 50, 32, 52, 85, 53], it was shown in the BdG equation [82] that
a singly quantized vortex always splits into two half-quantized non-Abelian
vortices with any strength of the magnetic field, forming a vortex
molecule.111 Similar molecules of half-quantized vortices connected by a
linear soliton are present in the other multicomponent condensed matter
physics and high energy physics: multicomponent or multigap superconductors
[112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125],
coherently-coupled multicomponent BECs [126, 93, 94, 127, 128, 129, 130, 95,
131, 132, 133, 134, 135, 99, 105, 136], high density QCD [137], and two-Higgs
doublet models [138]. However, the unique feature of ${}^{3}P_{2}$ superfluids
is that constituent half-quantized vortices are non-Abelian vortices
associated with a non-Abelian first homotopy group.  (It was also found in
Ref. [82] that a Majorna fermion zero mode is trapped in each half-quantized
vortex.) Such a splitting of a singly quantized vortex was also confirmed in
the GL theory without enforcing axisymmetry [86]. In the GL theory, cores of
two half-quantized vortices exhibit a ferromagnetic order in the UN phase with
the zero magnetic field, and a cyclic order in the $D_{2}$ and $D_{4}$BN
phases in the presence of the magnetic field, as summarized in Table 1(e). In
the UN phase in the absence of the magnetic field, the most stable singly-
quantized vortex configuration consists of two half-quantized vortices with
the ferromagnetic cores connected by a single soliton of the $D_{4}$BN order,
and in addition to this, there is also a metastable configuration consisting
of those of the cyclic cores connected by three solitons of the $D_{4}$BN
order. On the other hand, in the $D_{2}$ ($D_{4}$)BN phase in the presence of
small (large) magnetic field, two half-quantized vortices are connected by
three linear solitons of the $D_{4}$ ($D_{2}$)BN order. See Table 1(f) and
(g). Even in the bulk UN and $D_{2}$BN phases, the $D_{4}$BN order locally
appears around the vortex cores as solitons, because isolated half-quantized
vortices can topologically exist only in the $D_{4}$BN state and thus
splitting into two half-quantized vortices is possible only inside the
$D_{4}$BN order.

In this paper, we investigate proximity effects of two vortex molecules. As
explained above, a singly quantized vortex is of the form of a vortex molecule
of two half-quantized non-Abelian vortices connected by a single soliton
(three solitons) in the absence (presence) of magnetic field. In the absence
of magnetic field, we find that a transition of half-quantized non-Abelian
vortex cores occurs from a ferromagnetic order (connected by a single soliton)
to a cyclic order (connected by three solitons), when two vortex molecules
come close to each other with increasing the rotation speed. Furthermore, one
of three solitons reconnects to one of the other vortex molecule and bridges
the two vortex molecules, forming a “covalent bond.” Thus, the four
constituent half-quantized vortices are connected by one or two soliton(s)
alternately. In the presence of magnetic field below the critical magnetic
field $B<B_{c}$, we find the same transition from two isolated vortex
molecules to a “dimerized molecule” with a covalent bonding of solitons when
the two vortex molecules come close to each other with increasing the rotation
speed 222 A similar reconnection of solitons was observed in numerical
simulations of two-component BECs [95, 99], and similar chemical bonds of
vortex molecules were found in multi-component superconductors with arbitrary
charges and Josephson couplings [124]. . However, we do not find such
phenomena in the $D_{4}$BN phase above the critical magnetic field. Instead,
we find another type of a transition for the polarization direction of the
vortex molecules from one parallel to the separation of the two vortex
molecules to the other perpendicular to the separation of the two vortex
molecules with increasing the external rotation. Then, we further study three
and four vortex molecules. In the $D_{2}$BN phase below the critical magnetic
field, as in the case of two vortex molecules, we find a transition from
isolated vortex molecules to trimerized and tetramerized vortex molecules
consisting of six and eight half-quantized non-Abelian vortices.

This paper is organized as follows. In Sec. II, we begin with formulations of
${}^{3}P_{2}$ superfluids within the GL approach in our notation, and shortly
summarize our previous results for a single vortex molecule state. In Sec.
III, we show our main results for two, three and four vortex molecule states.
Sec. IV is devoted to a summary and discussion.

##  II Ginzburg-Landau free energy and single-vortex molecule state for
${}^{3}P_{2}$ neutron superfluids

We start from a brief review of the GL theory for ${}^{3}P_{2}$ superfluids
[56] reformulated in the notation of Ref. [62] and single-vortex molecule
state within the GL formalism. The details were discussed in Refs. [56, 62]
for GL theory and Ref. [86] for single-vortex state.

###  II.1 Ginzburg-Landau theory

The effective GL Lagrangian density $f$ is given by

$\displaystyle\begin{aligned}
f&=K_{0}\left(f_{202}^{(0)}+f_{202}^{(1)}\right)+\alpha
f_{002}+\beta_{0}f_{004}+\gamma_{0}f_{006}\\\
&+\delta_{0}f_{008}+\beta_{2}f_{022}+\gamma_{2}f_{024}\\\
&+\sum_{4l+2m+n=10}\mathcal{O}(\nabla^{l}|\boldsymbol{B}|^{m}A^{n}),\end{aligned}$
(1)

where $f_{lmn}$ is the free energy part including $l$ spatial derivatives
$\nabla$, $m$-th order of the magnetic field $\boldsymbol{B}$, and $n$-th
order of spin-2 spinor order parameter
$\psi=(\psi_{2},\psi_{1},\psi_{0},\psi_{-1},\psi_{-2})^{T}$. The spatial
derivative term $f_{202}$ is further separated into current-spin independent
and dependent parts $f_{202}^{(0)}$ and $f_{202}^{(1)}$, respectively. Their
specific forms can be written as

$\displaystyle\begin{aligned}
&f_{202}^{(0)}=3\boldsymbol{j}^{\dagger}\cdot\boldsymbol{j},\quad
f_{202}^{(1)}=4\boldsymbol{j}^{\dagger}\cdot\boldsymbol{j}-\frac{i}{2}\boldsymbol{j}^{\dagger}\cdot\hat{\boldsymbol{S}}\times\boldsymbol{j}-\left(\boldsymbol{j}^{\dagger}\cdot\hat{\boldsymbol{S}}\right)\left(\hat{\boldsymbol{S}}\cdot\boldsymbol{j}\right),\quad
f_{002}=3\rho,\quad
f_{004}=6\rho^{2}+\frac{3}{4}\boldsymbol{S}^{2}-\frac{3}{2}|\Psi_{20}|^{2},\\\
&f_{022}=2\rho\boldsymbol{B}^{2}-\frac{1}{2}\psi^{\dagger}\hat{S}_{\boldsymbol{B}}\hat{S}_{\boldsymbol{B}}\psi,\quad
f_{006}=-324\rho^{3}-81\rho\boldsymbol{S}^{2}+162\rho|\Psi_{20}|^{2}+15|\Psi_{30}|^{2}-27|\Phi_{30}|^{2},\\\
&f_{024}=\left(-106\rho^{2}+\frac{9}{2}\boldsymbol{S}^{2}+31|\Psi_{20}|^{2}\right)\boldsymbol{B}^{2}\\\
&\phantom{f_{024}}+\left(22\rho\psi^{\dagger}\hat{S}_{\boldsymbol{B}}\hat{S}_{\boldsymbol{B}}\psi+\mathrm{Re}\left[\Psi_{20}^{\ast}\psi^{T}\hat{S}_{\boldsymbol{B}}^{T}\hat{T}\hat{S}_{\boldsymbol{B}}\psi\right]\vphantom{\frac{5}{4}}+\frac{5}{4}\Psi_{22}^{\dagger}\hat{S}_{\boldsymbol{B}}\hat{S}_{\boldsymbol{B}}\Psi_{22}+\frac{1}{2}\Phi_{22}^{T}\hat{S}_{\boldsymbol{B}}^{T}\hat{T}\hat{S}_{\boldsymbol{B}}\Phi_{22}\right),\\\
&f_{008}=6480\rho^{4}+1944\rho^{2}\boldsymbol{S}^{2}-5184\rho^{2}|\Psi_{20}|^{2}-864\rho|\Psi_{30}|^{2}+2592\rho|\Phi_{30}|^{2}+81\boldsymbol{S}^{4}+648|\Psi_{20}|^{4}-1296\Gamma_{4}\\\
&\boldsymbol{j}\equiv-i\nabla\psi,\quad\rho\equiv\psi^{\dagger}\psi,\quad\boldsymbol{S}\equiv\psi^{\dagger}\hat{\boldsymbol{S}}\psi,\quad\Gamma_{4}\equiv\mathrm{Re}\left[\Psi_{20}\Phi_{30}^{\ast
2}\right],\quad\Psi_{20}\equiv\sqrt{5}C^{00}_{2m_{1},2m_{2}}\psi_{m_{1}}\psi_{m_{2}},\\\
&\Psi_{30}\equiv-\sqrt{\frac{35}{2}}C^{00}_{JM,2m_{3}}C^{JM}_{2m_{1},2m_{2}}\psi_{m_{1}}\psi_{m_{2}}\psi_{m_{3}},\quad\Phi_{30}\equiv-\sqrt{\frac{35}{2}}C^{00}_{JM,2m_{3}}C^{JM}_{2m_{1},2m_{2}}\psi_{m_{1}}\psi_{m_{2}}(-1)^{m_{3}}\psi^{\ast}_{-m_{3}},\end{aligned}$
(2)

where $\hat{S}_{i}$ ($i=x,y,z$) are $5\times 5$ spin-2 matrices. For
$\Psi_{20}$, $\Psi_{30}$, and $\Phi_{30}$, we take the Einstein summation
notation for $-2\leq m_{1,2,3}\leq 2$, $0\leq J\leq 4$, and $-J\leq M\leq J$
with the Clebsch-Gordan coefficients $C^{JM}_{s_{1}m_{1},s_{2}m_{2}}$. The GL
coefficients can be obtained in the weak coupling limit within the quasi-
classical approximation starting from the non-relativistic spin-1/2 fermion
field theory as [56]

$\displaystyle K_{0}=\frac{7\zeta(3)N(0)p_{\rm F}^{4}}{240\pi^{2}m_{\rm
n}^{2}T^{2}},\quad$ $\displaystyle\alpha=\frac{N(0)p_{\rm
F}^{2}}{3}\log\frac{T}{T_{\rm c}},$ (3)
$\displaystyle\beta_{0}=\frac{7\zeta(3)N(0)p_{\rm
F}^{4}}{60\pi^{2}T^{2}},\quad$
$\displaystyle\gamma_{0}=\frac{31\zeta(5)N(0)p_{\rm
F}^{6}}{13440\pi^{4}T^{4}},$
$\displaystyle\delta_{0}=\frac{127\zeta(7)N(0)p_{\rm
F}^{8}}{387072\pi^{6}T^{6}},\quad$
$\displaystyle\beta_{2}=\frac{7\zeta(3)N(0)p_{\rm F}^{2}\gamma_{\rm
n}^{2}}{48(1+F_{0}^{a})^{2}\pi^{2}T^{2}},$
$\displaystyle\gamma_{2}=\frac{31\zeta(5)N(0)p_{\rm F}^{4}\gamma_{\rm
n}^{2}}{3840(1+F_{0}^{a})^{2}\pi^{4}T^{4}},$

with the temperature $T$, the critical temperature $T_{\rm c}$, the neutron
mass $m_{\rm n}$, the neutron gyromagnetic ratio $\gamma_{\rm n}$, the Fermi
momentum $p_{\rm F}$, the state-number density $N(0)=m_{\rm n}p_{\rm
F}/(2\pi)^{2}$ at the Fermi surface, and the Landau parameter $F_{0}^{a}$.

Phase | $\boldsymbol{S}^{2}/\rho^{2}$ | $|\Psi_{20}|^{2}/\rho^{2}$ | $|\Psi_{30}|^{2}/\rho^{3}$
---|---|---|---
F | 4 | 0 | 0
UN | 0 | 1 | 1
$D_{4}$BN | 0 | 1 | 0
$D_{2}$BN | 0 | 1 | $(0,1)$
C | 0 | 0 | 2
Table 2:  $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$,
$|\Psi_{20}|^{2}$, and $|\Psi_{30}|^{2}$ for ferromagnetic (F), uniaxial
nematic (UN), $D_{4}$ biaxial nematic (BN), $D_{2}$BN, and cyclic (C) states.

All candidates for uniform ground states were classified in Ref. [62]. The
five characteristic symmetric states are ferromagnetic (F), uniaxial nematic
(UN), $D_{4}$ biaxial nematic (BN), $D_{2}$BN, and cyclic (C) states. Each
ground state is characterized by $U(1)\times SO(3)$ invariants
$\boldsymbol{S}^{2}$, $|\Psi_{20}|^{2}$, and $|\Psi_{30}|^{2}$ as summarized
in Table 2.

For the effective Lagrangian density $f$ in Eq. (1), the UN, $D_{2}$BN, and
$D_{4}$BN states are predicted to be realized as the ground states of
${}^{3}P_{2}$ superfluids at $|\boldsymbol{B}|=0$, $0<|\boldsymbol{B}|<B_{\rm
c}$, and $|\boldsymbol{B}|>B_{c}$, [79] as shown in Table 1(a). The critical
magnetic field $B_{\rm c}$ separating the $D_{2}$BN and $D_{4}$BN states
depends on the temperature and takes the maximum value $B_{\rm c}=7.06\times
10^{-2}\pi(1+F_{0}^{a})T_{\rm c}/\gamma_{\rm n}$ at $T\simeq 0.854T_{\rm c}$.
With an estimation for the critical temperature $T_{\rm c}\approx 0.2$ MeV and
the Landau parameter $F_{0}^{a}\approx 1$, this critical magnetic field can be
estimated as $B_{\rm c}\approx 7.36\times 10^{15}$ G. At $T\lesssim
0.796T_{\rm c}$, we obtain $B_{\rm c}=0$.

###  II.2 Single vortex molecule solutions

Figure 1:  Singly quantized vortices as molecules of two half-quantized
vortices. $\boldsymbol{S}^{2}$ (left panels), $|\Psi_{20}|^{2}$ (middle
panels), and $|\Psi_{30}|^{2}$ (right panels) for the single-vortex solution
at (a) $T=0.854T_{\rm c}$ and $\boldsymbol{B}=0$, (b)
$\boldsymbol{B}=0.5B_{\rm c}\hat{\boldsymbol{z}}$ and (c)
$\boldsymbol{B}=1.3B_{\rm c}\hat{\boldsymbol{z}}$. The radius of the figure
shown here is $64p_{\rm F}/(\pi m_{\rm n}T_{\rm c})\approx 9.29\times 10^{3}$
fm.

Here, we briefly summarize our previous results [86] for singly quantized
vortex states as molecules of two half-quantized vortices.

First, we consider the ansatz for vortex solutions with vortex cores placed at
$r=0$ in the cylindrical coordinates $(r,\theta,z)$ and the boundary
$\psi_{r\to\infty}$ far from vortex cores. For singly quantized vortices, the
order parameters behave as $\psi_{-2\leq m\leq 2}|_{r\to\infty}\propto
e^{i\theta}$. For $\boldsymbol{B}=0$, the uniform ground state is degenerate
within the possible UN state

$\displaystyle\begin{aligned} &\psi_{\pm 2}=\frac{e^{i\phi\mp
2a}\sqrt{3}\sin^{2}b}{2\sqrt{2\rho}},\quad\psi_{\pm 1}=\frac{e^{i\phi\mp
a}\sin(2b)}{2\sqrt{2\rho}},\\\
&\psi_{0}=\frac{e^{i\phi}\\{1+3\cos(2b)\\}}{4\sqrt{\rho}},\end{aligned}$ (4)

for $0\leq a,2b,\phi<2\pi$. Here $a$, $b$, and $\phi$ represent overall spin
rotations along the $z$ axis and $y$ axis, and overall phase shift,
respectively. Under the spatial phase gradient $e^{i\theta}$ for the vortex
solution, the current-spin dependent free energy density $f_{202}^{(2)}$ in
Eq. (1) favors $b=0$ giving

$\displaystyle\psi\xrightarrow{r\to\infty}\frac{e^{i\theta}}{\sqrt{\rho}}(0,0,1,0,0)^{T}.$
(5)

In the case of $0<|\boldsymbol{B}|<B_{\rm c}$, we obtain

$\displaystyle\psi\xrightarrow{r\to\infty}\frac{e^{i\theta}}{\sqrt{\rho}}\left(\frac{e^{-2ia}\sin
g}{\sqrt{2}},0,\cos g,0,\frac{e^{2ia}\sin g}{\sqrt{2}}\right)^{T},$ (6)

where $g$ depends on $|\boldsymbol{B}|$ and satisfies $\pi/3<g<\pi/2$, making
$\psi|_{r\to\infty}$ to be the $D_{2}$BN state. $a$ also represents the
overall spin rotation along the $z$ axis and takes arbitrary (fixed) value
without (with) the current-spin dependent free energy $f_{202}^{(1)}$. In the
limit of $|\boldsymbol{B}|\searrow 0$, $g$ becomes $g\to\pi/3$ giving

$\displaystyle\psi\xrightarrow{r\to\infty}\frac{e^{i\theta}}{\sqrt{\rho}}\left(\frac{e^{-2ia}\sqrt{3}}{2\sqrt{2}},0,\frac{1}{2},0,\frac{e^{2ia}\sqrt{3}}{2\sqrt{2}}\right)^{T}.$
(7)

This solution belongs to the UN state in Eq. (4) with $b=\pi/2$, but is
different from that for $\boldsymbol{B}=0$ shown in Eq. (5), which leads the
discontinuity between $\boldsymbol{B}=0$ and $|\boldsymbol{B}|\searrow 0$ and
the metastable solution at $\boldsymbol{B}=0$. For the case of
$|\boldsymbol{B}|\geq B_{\rm c}$, $g$ becomes $g=\pi/2$, giving

$\displaystyle\psi\xrightarrow{r\to\infty}\frac{e^{i\theta}}{\sqrt{\rho}}\left(\frac{e^{-2ia}}{\sqrt{2}},0,0,0,\frac{e^{2ia}}{\sqrt{2}}\right)^{T},$
(8)

which belongs to the $D_{4}$BN state.

Next, we show our numerical results for singly quantized vortex solutions.
They are obtained by minimizing the free-energy density $f$ under the boundary
conditions with the cylindrical coordinate $(r,\theta,z)$:

$\displaystyle\left.\psi_{m}(\theta+\pi)\right|_{r=L}=-\left.\psi_{m}(\theta)\right|_{r=L},$
(9)

at the boundary $r=L$, which induces the integer vortex solution. The
minimization of the free energy density $f$ can be done by finding the
stationary solution of the GL equation

$\displaystyle\frac{\delta f}{\delta\psi_{m}^{\ast}}=0.$ (10)

The solution of Eq. (10) can be obtained by the Nesterov’s method with
introducing the relaxation time $t$ and the dependence of the order parameter
$\psi_{m}$. The time dependences of $\psi_{m}$ is given by

$\displaystyle\begin{aligned} \ddot{\psi}_{m}=-\frac{\delta
f}{\delta\psi_{m}^{\ast}}-\frac{3}{t}\dot{\psi}_{m},\quad\dot{\psi}_{m}(t=0)=0.\end{aligned}$
(11)

After the long time evolution of Eq. (11), we obtain the solution of Eq. (10).
We note that Eq. (11) is just one of methods to effectively obtain the
solution of Eq. (10).

Figure 1 shows numerical solutions for single vortex molecule states [86] at
the temperature $T=0.854T_{\rm c}$ and the magnetic field $\boldsymbol{B}=0$,
$0.5B_{\rm c}\hat{\boldsymbol{z}}$, and $1.3B_{\rm c}\hat{\boldsymbol{z}}$. In
any cases, a singly quantized vortex splits into two half-quantized vortices
with holes of $|\Psi_{20}|^{2}$ forming a vortex molecule. The fact that
isolated half-quantized vortices can topologically exist only in the $D_{4}$BN
state implies that the $D_{4}$BN order should appear around the vortex core
for even for $|\boldsymbol{B}|<B_{\rm c}$ in which the ground states are
either UN or $D_{2}$BN states (Figs. 1 (a), (b), and (c)). In fact, we can
confirm that the $D_{4}$BN order characterized by $\boldsymbol{S}^{2}=0$,
$|\Psi_{20}|^{2}/\rho^{2}=1$, and $|\Psi_{30}|^{2}/\rho^{3}=0$ appears as one
(three) soliton(s) bridging two vortex cores for Fig. 1 (a) at
$\boldsymbol{B}=0$ (Fig. 1 (b) at $\boldsymbol{B}=0$ and Fig. 1 (c) at
$\boldsymbol{B}=0.5B_{\rm c}\hat{\boldsymbol{z}}$), as can be seen in the plot
of $|\Psi_{30}|^{2}$ locally inducing the $D_{4}$BN order. At
$\boldsymbol{B}=1.3B_{\rm c}\hat{\boldsymbol{z}}$ shown in Fig. 1 (d), where
$D_{4}$BN state becomes the ground state, two half-quantized vortices also
form a vortex molecule bridged by three $D_{2}$BN solitons characterized by
$0<|\Psi_{30}|^{2}/\rho^{3}<1$. The cores of the half-quantized vortices are
filled by F, C, C, and C orders for Figs. 1 (a), (b), (c), and (d),
respectively. At the numerical boundary $r=L$, the order parameters for Figs.
1 (a), (b), (c), and (d) satisfy Eqs. (5), (7), (6), and (8), respectively. In
particular, the two solutions at $\boldsymbol{B}=0$ shown in Figs. 1 (a) and
(b) correspond to the most stable and metastable ones, respectively.

##  III Multi-vortex state for ${}^{3}P_{2}$ Neutron superfluids

In this section, we show our main results for multi-vortex state for
${}^{3}P_{2}$ superfluids. Because two vortices with the same circulation have
long-range repulsion, they cannot be stabilized by the boundary condition such
as Eq. (9). We thus introduce an external rotation term by rewriting the free-
energy density $f$ as $f-\mathrm{Re}[\boldsymbol{\Omega}\cdot\boldsymbol{L}]$
[92], where $\boldsymbol{\Omega}=(0,0,\Omega)$ is the rotation vector parallel
to the $z$-axis, and $\boldsymbol{L}=m_{\rm n}^{2}p_{\rm
F}\psi^{\dagger}(-i\boldsymbol{r}\times\nabla)\psi$ is the angular momentum.
This is nothing but a physical situation realized in rotating neutron stars.
As well as the case of the single-vortex solutions, we obtain the stationary
solution of the GL equation by the Nesterov’s method in Eq. (11) under the
boundary condition

$\displaystyle\left.\left(\nabla\psi\right)_{\parallel}\right|_{r=L}=0,$ (12)

where $(\ )_{\parallel}$ means the component parallel to
$\hat{\boldsymbol{r}}$ in the cylindrical coordinate, and $L$ is the numerical
boundary set to be $L=128p_{\rm F}/(\pi m_{\rm n}T_{\rm c})\approx 1.86\times
10^{4}$ fm.

###  III.1 Two-vortex molecule states

####  III.1.1 The stable solution in the case of zero magnetic field
$\boldsymbol{B}=0$

Figure 2:  The most stable two-vortex molecule state in the UN phase at
$\boldsymbol{B}=0$ and $\Omega^{2}=0.8\Omega_{0}^{2}$ with
$\Omega_{0}=160\pi^{2}m_{n}T_{\rm c}^{2}/(256^{2}p_{\rm F}^{2})\approx
7.15\times 10^{15}$ rad$\cdot$fs-1. The squared modulus $|\psi_{m}|^{2}$ (top
row), argument $\mathrm{Arg}[\psi_{m}]$ (middle row) of the order parameter,
the $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$, $|\Psi_{20}|^{2}$, and
$|\Psi_{30}|^{2}$, and the free-energy density $f$ (bottom row) are shown. The
radius of the figure shown here is $64p_{\rm F}/(\pi m_{\rm n}T_{\rm
c})\approx 9.29\times 10^{3}$ fm.

For the $\boldsymbol{B}=0$ case in which the UN phase is the ground state, we
expect the F-core vortex molecules connected by a single $D_{4}$BN soliton as
shown in Fig. 1 (a). Figure 2 shows the two-vortex molecule state at
$\boldsymbol{B}=0$ and $\Omega^{2}=0.8\Omega_{0}^{2}$, where $\Omega_{0}$ is
defined as $\Omega_{0}=160\pi^{2}m_{n}T_{\rm c}^{2}/(256^{2}p_{\rm
F}^{2})\approx 7.30$ rad$\cdot$fs-1. Compared to the singly quantized vortex
state, the two $D_{4}$BN solitons between half-quantized vortex cores repel
each other to be bend, as can be clearly seen in $|\Psi_{30}|^{2}$. At the
boundary $r=L$, the order parameter approximately satisfies

$\displaystyle\psi|_{r=L}=\frac{e^{2i\theta}}{\sqrt{\rho}}(0,0,1,0,0)^{T},$
(13)

which is just the double winding of the singly quantized vortex solution Eq.
(5) at $\boldsymbol{B}=0$.

Figure 3:  The most stable two-vortex molecule state in the UN phase at
$\boldsymbol{B}=0$ and $\Omega^{2}=1.0\Omega_{0}^{2}$. The squared modulus
$|\psi_{m}|^{2}$ (top row), argument $\mathrm{Arg}[\psi_{m}]$ (middle row) of
the order parameter, the $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$,
$|\Psi_{20}|^{2}$, and $|\Psi_{30}|^{2}$, and the free-energy density $f$
(bottom row) are shown. The radius of the figure shown here is $64p_{\rm
F}/(\pi m_{\rm n}T_{\rm c})\approx 9.29\times 10^{3}$ fm. The two vortex
molecules share two solitons.

With increasing the rotation, the “dimerization” of two vortex molecules
occurs while keeping the boundary state unchanged as shown in Eq. (13). Fig. 3
shows the two-vortex molecule state consisting of four half-quantized vortices
at $\Omega^{2}=1.0\Omega_{0}^{2}$. The details of the dimerization is as
follows: First, the F-core vortex molecules having one $D_{4}$BN soliton
changes to C-core vortex molecules having three $D_{4}$BN solitons as in the
case of the vortex molecule at $0<|\boldsymbol{B}|<B_{\rm c}$. Then, one of
three $D_{4}$BN solitons reconnect to one of the other molecule to be shared
by the two vortex molecules to form a “covalent bond”. In Fig. 3, the upper
and lower $D_{4}$BN solitons in $|\Psi_{30}|^{2}$ are covalent bonds between
left and right vortex molecules. The transition from the two isolated F-core
vortex molecules to a dimerized C-core vortex molecule occurs at
$\Omega^{2}\simeq\Omega_{0}^{2}$.

####  III.1.2 The metastable solution in the case of zero magnetic field
$\boldsymbol{B}=0$

Figure 4:  The metastable two-vortex molecule state in the UN phase at
$\boldsymbol{B}=0$ and $\Omega^{2}=0.8\Omega_{0}^{2}$. The squared modulus
$|\psi_{m}|^{2}$ (top row), argument $\mathrm{Arg}[\psi_{m}]$ (middle row) of
the order parameter, the $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$,
$|\Psi_{20}|^{2}$, and $|\Psi_{30}|^{2}$, and the free-energy density $f$
(bottom row) are shown. The radius of the figure shown here is $64p_{\rm
F}/(\pi m_{\rm n}T_{\rm c})\approx 9.29\times 10^{3}$ fm.

As well as singly quantized vortex states, there are metastable solutions
having the other boundary

$\displaystyle\psi|_{r=L}=\frac{e^{2i\theta}}{\sqrt{\rho}}\left(\frac{e^{-2ia}\sqrt{3}}{2\sqrt{2}},0,\frac{1}{2},0,\frac{e^{2ia}\sqrt{3}}{2\sqrt{2}}\right)^{T},$
(14)

which is double winding of Eq. (7). In this case, we expect the C-core vortex
molecule connected by three $D_{4}$BN solitons as shown in Fig. 1 (b). Figure
4 shows the metastable two-vortex molecule state at $\boldsymbol{B}=0$ and
$\Omega^{2}=0.8\Omega_{0}^{2}$. There are two isolated vortex molecules in
which the three $D_{4}$BN solitons connect half-quantized vortices. Compared
to the single vortex state shown in Fig. 1 (b), the $D_{4}$BN soliton in the
side of the other vortex molecule is shorten. Oh the other hand, the $D_{4}$BN
soliton in the opposite side is enlarged and bent. As a result, a symmetric
structure of one vortex molecule as shown in Fig. 1 (b) is strongly distorted.

Figure 5:  The metastable two-vortex molecule state in the UN phase at
$\boldsymbol{B}=0$ and $\Omega^{2}=1.0\Omega_{0}^{2}$. The squared modulus
$|\psi_{m}|^{2}$ (top row), argument $\mathrm{Arg}[\psi_{m}]$ (middle row) of
the order parameter, the $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$,
$|\Psi_{20}|^{2}$, and $|\Psi_{30}|^{2}$, and the free-energy density $f$
(bottom row) are shown. The radius of the figure shown here is $64p_{\rm
F}/(\pi m_{\rm n}T_{\rm c})\approx 9.29\times 10^{3}$ fm.

As well as the stable case, increasing $\Omega$ causes the dimerization of
vortex molecules under the same boundary condition Eq. (14). The shape of the
dimerized vortex molecule is similar to that shown in Fig. 3, i.e., there are
two $D_{4}$BN covelent bond shared by the two vortex molecules.

####  III.1.3 The case of a small magnetic field $0<|\boldsymbol{B}|<B_{\rm
c}$

Figure 6:  The two-vortex molecule state in the $D_{2}$BN phase at
$\boldsymbol{B}=0.5B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=0.8\Omega_{0}^{2}$. The squared modulus $|\psi_{m}|^{2}$ (top
row), argument $\mathrm{Arg}[\psi_{m}]$ (middle row) of the order parameter,
the $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$, $|\Psi_{20}|^{2}$, and
$|\Psi_{30}|^{2}$, and the free-energy density $f$ (bottom row) are shown. The
radius of the figure shown here is $64p_{\rm F}/(\pi m_{\rm n}T_{\rm
c})\approx 9.29\times 10^{3}$ fm.

Figure 7:  The two-vortex molecule state in the $D_{2}$BN phase at
$\boldsymbol{B}=0.5B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=1.0\Omega_{0}^{2}$. The squared modulus $|\psi_{m}|^{2}$ (top
row), argument $\mathrm{Arg}[\psi_{m}]$ (middle row) of the order parameter,
the $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$, $|\Psi_{20}|^{2}$, and
$|\Psi_{30}|^{2}$, and the free-energy density $f$ (bottom row) are shown. The
radius of the figure shown here is $64p_{\rm F}/(\pi m_{\rm n}T_{\rm
c})\approx 9.29\times 10^{3}$ fm.

For the $0<|\boldsymbol{B}|<B_{\rm c}$ case in which the $D_{2}$BN phase is
the ground state, we expect the C-core vortex molecule connected by three
$D_{4}$BN solitons as shown in Fig. 1 (c). Figure 6 (Figure 7) shows two-
vortex molecule state at $\boldsymbol{B}=0.5B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=0.8\Omega_{0}^{2}$ ($\Omega^{2}=1.0\Omega_{0}^{2}$). Two
characteristic structures for two isolated vortex molecules shown in Fig. 6
and the dimerized vortex molecule shown in Fig. 7 are almost the same as those
for the metastable solutions at $\boldsymbol{B}=0$ for Figs. 4 and 5, while
the order parameter at the boundary satisfy

$\displaystyle\psi|_{r=L}=\frac{e^{2i\theta}}{\sqrt{\rho}}\left(\frac{e^{-2ia}\sin
g}{\sqrt{2}},0,\cos g,0,\frac{e^{2ia}\sin g}{\sqrt{2}}\right)^{T},$ (15)

which is double winding of Eq. (6). As well as the singly quantized vortex
solution, the two-vortex molecule state in the $D_{2}$BN phase at
$0<|\boldsymbol{B}|<B_{\rm c}$ continuously changes to the metastable state in
the UN phase at $\boldsymbol{B}=0$ ($g\to\pi/3$ limit in Eq. (15) at the
boundary).

####  III.1.4 The case of large magnetic field $|\boldsymbol{B}|\geq B_{\rm
c}$

Figure 8:  The two-vortex molecule state in the $D_{4}$BN phase at
$\boldsymbol{B}=1.3B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=0.4\Omega_{0}^{2}$. The squared modulus $|\psi_{m}|^{2}$ (top
row), argument $\mathrm{Arg}[\psi_{m}]$ (middle row) of the order parameter,
the $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$, $|\Psi_{20}|^{2}$, and
$|\Psi_{30}|^{2}$, and the free-energy density $f$ (bottom row) are shown. The
radius of the figure shown here is $64p_{\rm F}/(\pi m_{\rm n}T_{\rm
c})\approx 9.29\times 10^{3}$ fm.

Finally, for the $|\boldsymbol{B}|\geq B_{\rm c}$ case in which the $D_{4}$BN
phase is the ground state, we expect the C-core vortex molecule connected by
three $D_{2}$BN solitons as shown in Fig. 1 (c). Figure 8 shows the two-vortex
molecule state at $\boldsymbol{B}=1.3B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=0.4\Omega_{0}^{2}$. There are two isolated vortex molecules in
each of which the three $D_{2}$BN solitons connect two half-quantized
vortices. At the boundary, the order parameter satisfies

$\displaystyle\psi|_{r=L}=\frac{e^{2i\theta}}{\sqrt{\rho}}\left(\frac{e^{-2ia}}{\sqrt{2}},0,0,0,\frac{e^{2ia}}{\sqrt{2}}\right)^{T},$
(16)

which is double winding of Eq. (8). We define the molecule angle $\theta_{\rm
molecule}$ as the angle between the polarization direction of each vortex
molecule and the direction of sepration of the centers of the two vortex
molecules. In Fig. 8, the molecule angle takes $\theta_{\rm rel}=0$. On the
other hand, for the $|\boldsymbol{B}|<B_{\rm c}$ cases shown in Figs. 2, 3,
and 6, the molecule angle takes $\theta_{\rm molecule}=\pi/2$.

Figure 9:  A transition of the polarization direction of the two-vortex
molecules in the $D_{4}$BN phase. $|\Psi_{30}|^{2}$ for two-vortex state in
the $D_{4}$BN phase at $\boldsymbol{B}=1.3B_{\rm c}\hat{\boldsymbol{z}}$ and
(a) $\Omega^{2}=0.4\Omega_{0}^{2}$, (b) $\Omega^{2}=0.6\Omega_{0}^{2}$, and
(c) $\Omega^{2}=0.8\Omega_{0}^{2}$. The radius of the figure shown here is
$64p_{\rm F}/(\pi m_{\rm n}T_{\rm c})\approx 9.29\times 10^{3}$ fm.

Being different from the $|\boldsymbol{B}|<B_{\rm c}$ cases, the dimerization
of vortex molecules never occurs even with increasing the rotation frequency
$\Omega$. On the other hand, a rapid change of $\theta_{\rm molecule}$ from
$0$ to $\pi/2$ occurs. Figure 9 shows $|\Psi_{30}|^{2}$ profiles for the two-
vortex molecule states at the rotation frequencies
$\Omega^{2}/\Omega_{0}^{2}=0.4$ (panel (a)), $0.6$ (panel (b)), and $0.8$
(panel (c)), and the molecule angle $\theta_{\rm molecule}$ takes $\theta_{\rm
molecule}\approx 0$, $0.22\pi$, and $\pi/2$, respectively.

####  III.1.5 Phase diagram

Figure 10:  The $\Omega$–$|\boldsymbol{B}|$ phase diagram for two-vortex
molecule states. ⚫: Coexistence of isolated stable F-core molecules and
isolated metastable C-core molecules. $\times$: Coexistence of isolated stable
F-core molecules, isolated metastable C-core molecules, dimerized stable
molecule, and dimerized metastable molecule. $-$: Coexistence of isolated
metastable C-core molecules, dimerized stable molecule, and dimerized
metastable molecule. $\blacksquare$: isolated C-core molecules. $+$:
Coexistence of isolated C-core molecules and dimerized molecule.
$\blacktriangle$: Dimerized molecule. $\blacklozenge$: Molecules with
$\theta_{\rm molecule}=0$. $\blacktriangledown$: Molecules with $\theta_{\rm
molecule}=\pi/2$. $\bigstar$: Molecules with $0<\theta_{\rm molecule}<\pi/2$.

Figure 10 shows the $\Omega$-$|\boldsymbol{B}|$ phase diagram for two-vortex
molecule states. The dimerized molecules for $|\boldsymbol{B}|<B_{\rm c}$ and
molecules with $\theta_{\rm molecule}=\pi/2$ can exist at large $\Omega$ and
small $|\boldsymbol{B}|$. This result supports that these states originate
from a proximity effect between two vortex molecules, because the distance
between molecules decreases with increasing $\Omega$ and the size of a
molecule decreases with increasing $|\boldsymbol{B}$ (see Fig. 5 (b) in Ref.
[86]). At $|\boldsymbol{B}|<B_{\rm c}$, the dimerized molecule state coexists
with the isolated molecule state in a wide region of the phase diagram
($\times$, $-$, and $+$ symbols). Here, the coexistence of several states
means that one state is stable and others are metastable, or that some states
are degenerate and equally stable.

Figure 11:  Two-dimensional free energy $F=\int d^{2}x\>f$ at
$|\boldsymbol{B}|=0$ as a function of $\Omega^{2}$ for the isolated stable
F-core molecule state, the isolated metastable C-core molecule state, the
dimerized stable molecule state, and the dimerized metastable molecule state.

To check this, we calculate the two-dimensional free energy $F=\int
d^{2}x\>f$. Figure 11 shows the two-dimensional free energy $F$ as a function
of $\Omega^{2}$ at $|\boldsymbol{B}|=0$ where the isolated stable F-core
molecule state (see Fig. 2), the isolated metastable C-core molecule state
(see Fig. 4), the dimerized stable molecule state (see Fig. 3), and the
dimerized metastable molecule state (see Fig. 5) are stabilized. The isolated
stable F-core molecule state and dimerized stable molecule state are stable
states at small and large rotation frequencies $\Omega$, respectively, having
lower free-energy than those of the isolated metastable C-core molecule state
and the dimerized metastable molecule state. On the other hand, the isolated
stable F-core (metastable C-core) molecule state and the dimerized stable
(metastable) molecule state have almost degenerate energies in a wide range of
$\Omega$, which does not change in the case of finite values of
$|\boldsymbol{B}|$ ($+$ symbols in Fig. 10).

###  III.2 Three and four-vortex molecule states

For more than two vortex molecule states, we expect various stable
configurations such as completely isolated molecules, partially polymerized
molecules, fully polymerized molecule at $|\boldsymbol{B}|<B_{\rm c}$, and
clusters of molecule pairs with various $\theta_{\rm molecule}$ at
$|\boldsymbol{B}|>B_{\rm c}$. Actually, it is quite difficult to exhaust all
(meta-)stable states.

Figure 12:  $|\Psi_{30}|^{2}$ for three-vortex states (upper panels) and four-
vortex states (lower panels) at (a) $\boldsymbol{B}=0.2B_{\rm
c}\hat{\boldsymbol{z}}$ and $\Omega^{2}=1.9\Omega_{0}^{2}$, (b)
$\boldsymbol{B}=0.2B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=5.0\Omega_{0}^{2}$, (c) $\boldsymbol{B}=1.3B_{\rm
c}\hat{\boldsymbol{z}}$ and $\Omega^{2}=5.0\Omega_{0}^{2}$, (d)
$\boldsymbol{B}=0.5B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=2.0\Omega_{0}^{2}$, and (e) $\boldsymbol{B}=0.4B_{\rm
c}\hat{\boldsymbol{z}}$ and $\Omega^{2}=8.0\Omega_{0}^{2}$, and (f)
$\boldsymbol{B}=1.5B_{\rm c}\hat{\boldsymbol{z}}$ and
$\Omega^{2}=8.0\Omega_{0}^{2}$. The radius of the figure shown here is
$64p_{\rm F}/(\pi m_{\rm n}T_{\rm c})\approx 9.29\times 10^{3}$ fm.

Instead, here we just show several characteristic examples for three and four-
vortex states. Figure 12 shows six examples of three and four vortex molecule
states. At $|\boldsymbol{B}|<B_{\rm c}$, we have various kinds of
“polymerization” of vortex molecules. In Fig. 12 (a) and (d), we have one
dimerized molecule and one isolated molecule, and two dimerized molecules, as
three and four vortex molecule states, respectively. We further obtain
“trimerized” and “tetramerized” molecules in Figs. 12 (b) and (d),
respectively. At $|\boldsymbol{B}|\geq B_{\rm c}$, there are also various
kinds of clusters of vortex molecules. We show symmetric examples of them in
Figs. 12 (c) and (f).

##  IV Summary and Discussion

In this paper, we have worked out solutions for two vortex molecules
consisting of four half-quantized non-Abelian vortices in the neutron
${}^{3}P_{2}$ superfluids in the presence of the external magnetic field
parallel to the angular momentum of the vortices. We have found two
characteristic transitions at $|\boldsymbol{B}|<B_{\rm c}$ and
$|\boldsymbol{B}|\geq B_{\rm c}$ as a result of the proximity effect between
the two vortex molecules. At $0<|\boldsymbol{B}|<B_{\rm c}$ where the C-core
vortex molecule with three $D_{4}$BN solitons is stabilized in the $D_{2}$BN
ground state, two isolated vortex molecules dimerize through covalent bonds of
$D_{4}$BN solitons namely the solitons shared by the two molecules. At
$\boldsymbol{B}=0$ where the F-core vortex molecule with a $D_{4}$BN soliton
is the most stable state in the UN ground state, not only dimerization but
also the transition of the vortex core state from F state to C state
simultaneously occurs. There are also metastable isolated C-core vortex
molecules and the dimerized vortex molecule which are continuously changed
from states at $0<|\boldsymbol{B}|<B_{\rm c}$. Our results suggest that a
dimerized F-core vortex molecule never exist. This situation is similar to the
chemical dimerization in which two molecules having only a $\sigma$ bond
cannot dimerize and those having more than two bonds consisting of a $\sigma$
bond and one or two $\pi$ bonds can dimerize. Because the F-core vortex
molecule has only one $D_{4}$BN soliton, the vortex core needs to change from
F state to C state with three $D_{4}$BN solitons for isolated F-core vortex
molecules to dimerize. At $|\boldsymbol{B}|\geq B_{\rm c}$ where the C-core
vortex molecule with three $D_{2}$BN soliton is stabilized in $D_{4}$BN ground
state, no dimerization occurs. Instead, a transition of the molecule angle
from $\theta_{\rm molecule}=0$ to $\theta_{\rm molecule}=\pi/2$ occurs. The
main results are summarized as the $|\boldsymbol{B}|$-$\Omega$ phase diagram
shown in Fig. 10. We have also obtained some trimerized and tetramerized
vortex molecule states as three and four vortex molecule solutions,
respectively.

Here, let us give discussions for future studies. We have exhausted possible
two-vortex molecule states, and have briefly given some examples of three and
four vortex molecule states since these states are quite nontrivial. Further
studies should be done for a vortex-molecule lattice under a rapid rotation
relevant for neutron star interiors. In particular, it is important to
investigate which is realized among a symmetric vortex-molecule lattice as
found in two-component BECs [93, 94, 95] or a non-symmetric glassy state
realized under the complicated polimerization of vortex molecules. We here
have studied only the case for the magnetic field parallel to the the
direction of vortices. It will be one of important directions for neutron star
interiors to investigate the case of an arbitrary angle between them.

Apart from steady states under rotation studied thus far, a dynamics of
vortices in three spatial dimensions is also crucial for neutron star dynamics
such as the pulsar glitch phenomenon. The most important feature of vortex
dynamics in ${}^{3}P_{2}$ superfluids would be the non-Abelian property of
half-quantized vortices implying that they are noncommutative under exchange.
It is important whether two vortices reconnect in collision or a formation of
a rung between them occurs as the cases of non-Abelian vortices in the cyclic
phase of a spin-2 BEC [96] and the nematic phase of a spin-2 BEC [97, 61]. For
instance, a vortex reconnection is crucial for states of the quantum
turbulence [98]. Collision of two vortex molecules may be accompanied by
swapping partners as in vortex molecules in two-component BECs [99]. In
addition, the dimerization of vortex molecules found in this paper may occur
in a certain part of two vortex lines, thus forming for instance a X junction.
Thus, multiple vortex molecules are expected to be entangled in general.

In this paper, we have not considered a ${}^{1}S_{0}$ paring, but the phase
diagram under the coexistence of ${}^{1}S_{0}$ and ${}^{3}P_{2}$ superfluids
is quite different from the one of ${}^{3}P_{2}$ superfluids [100]. Vortex
states in this case will be also one direction to be explored. In addition,
vortex states in the ferromagnetic phase, that was found in the region close
to the critical temperature without quasi-classical approximation [63], are
also worth to be studied.

Quark matter consisting of diquark condensations exhibiting color
superconductivity may exist in the region deeper than nuclear matter in
neutron star cores. A quark-hadron continuity for two-flavor quarks was
suggested to continulusly connect the ${}^{3}P_{2}$ superfluid (nuclear
matter) to a two-flavor quark matter called the 2SC+dd phase through crossover
(rather than a phase transition) [101], and vortex structures in the 2SC+dd
phase were studied in Refs. [102, 103, 104]. It will be interesting to study
whether vortex molecule structures found in this paper is preserved or
deformed through the quark-hadron continuity.

A phase transition due to unbinding of a vortex and an anti-vortex is known as
the Berezinskii-Kosterlitz-Thouless (BKT) transition in two spatial dimension.
Recently, a novel type of BKT transition due to vortex molecules was found in
two-component systems [105]. Therefore, the BKT transition of neutron
${}^{3}P_{2}$ superfluids confined in a quasi-two dimensional plane should be
an interesting subject.

Finally, apart from ${}^{3}P_{2}$ superfluids, spin-2 spinor ultracold atomic
BECs are also $J=2$ condensates whose ground states are possibly nematic phase
[57, 58, 59, 60, 61, 62] although the current experiments of 87Rb atoms imply
their ground state to be in the cyclic or nematic phase [106, 107, 108, 109,
110, 111]. Nematic spin-2 BECs share almost the same bosonic properties with
${}^{3}P_{2}$ superfluids, and thus admit the same order parameter manifold
and non-Abelian half-quantized vortices [60, 61]. Therefore, our present
results for the dimerization of vortex molecules are also applicable to spin-2
nematic BECs which can be experimentally testable in principle.

## Acknowledgment

We would like to thank Yusuke Masaki for helpful discussions and comments.
This work is supported in part by JSPS KAKENHI (Grant Numbers JP22H01221 (M.K.
and M.N.), 20K03765(M.K.), 19KK0066(M.K.)), and by Osaka City University
Advanced Mathematical Institute (MEXT Joint Usage/Research Center on
Mathematics and Theoretical Physics JPMXP0619217849).

## References

  * Demorest  _et al._ [2010] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature 467, 1081 (2010), arXiv:1010.5788 [astro-ph.HE] .
  * Antoniadis  _et al._ [2013] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. W. Verbiest, and D. G. Whelan, Science 340, 1233232 (2013).
  * Abbott  _et al._ [2017] B. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2020] B. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020), arXiv:2001.01761 [astro-ph.HE] .
  * Riley  _et al._ [2019] T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ [2019] M. Miller _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Graber  _et al._ [2017] V. Graber, N. Andersson, and M. Hogg, Int. J. Mod. Phys. D 26, 1730015 (2017), arXiv:1610.06882 [astro-ph.HE] .
  * Baym  _et al._ [2018] G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Migdal [1960] A. B. Migdal, Zh. Eksp. Teor. Fiz. 37, 249 (1960), [Sov. Phys. JETP10, 176 (1960)].
  * Chamel and Haensel [2008] N. Chamel and P. Haensel, Living Rev. Rel. 11, 10 (2008), arXiv:0812.3955 [astro-ph] .
  * Chamel [2017] N. Chamel, J. Astrophys. Astron. 38, 43 (2017).
  * Haskell and Sedrakian [2018] B. Haskell and A. Sedrakian, Astrophys. Space Sci. Libr. 457, 401 (2018), arXiv:1709.10340 [astro-ph.HE] .
  * Sedrakian and Clark [2019] A. Sedrakian and J. W. Clark, Eur. Phys. J. A 55, 167 (2019), arXiv:1802.00017 [nucl-th] .
  * Andersson [2021] N. Andersson, Universe 7, 10.3390/universe7010017 (2021).
  * Baym  _et al._ [1969] G. Baym, C. Pethick, D. Pines, and M. Ruderman, Nature 224, 872 (1969).
  * Pines  _et al._ [1972] D. Pines, J. Shaham, and M. Ruderman, Nature Phys. Sci. 237, 83 (1972).
  * Takatsuka and Tamagaki [1988] T. Takatsuka and R. Tamagaki, Prog. Theor. Phys. 79, 274 (1988).
  * Yakovlev  _et al._ [2001] D. G. Yakovlev, A. D. Kaminker, O. Y. Gnedin, and P. Haensel, Phys. Rept. 354, 1 (2001), arXiv:astro-ph/0012122 [astro-ph] .
  * Potekhin  _et al._ [2015] A. Y. Potekhin, J. A. Pons, and D. Page, Space Sci. Rev. 191, 239 (2015), arXiv:1507.06186 [astro-ph.HE] .
  * Yakovlev  _et al._ [1999] D. G. Yakovlev, K. P. Levenfish, and Yu. A. Shibanov, Phys. Usp. 42, 737 (1999), arXiv:astro-ph/9906456 [astro-ph] .
  * Heinke and Ho [2010a] C. O. Heinke and W. C. G. Ho, Astrophys. J. 719, L167 (2010a), arXiv:1007.4719 [astro-ph.HE] .
  * Shternin  _et al._ [2011] P. S. Shternin, D. G. Yakovlev, C. O. Heinke, W. C. G. Ho, and D. J. Patnaude, Monthly Notices of the Royal Astronomical Society: Letters 412, L108 (2011).
  * Page  _et al._ [2011] D. Page, M. Prakash, J. M. Lattimer, and A. W. Steiner, Phys. Rev. Lett. 106, 081101 (2011), arXiv:1011.6142 [astro-ph.HE] .
  * Reichley and Downs [1971] P. E. Reichley and G. S. Downs, Nature Phys. Sci. 234, 48 (1971).
  * Anderson and Itoh [1975] P. W. Anderson and N. Itoh, Nature 256, 25 (1975).
  * Wolf [1966] R. A. Wolf, Astrophys. J. 145, 834 (1966).
  * Tabakin [1968] F. Tabakin, Phys. Rev. 174, 1208 (1968).
  * Hoffberg  _et al._ [1970] M. Hoffberg, A. E. Glassgold, R. W. Richardson, and M. Ruderman, Phys. Rev. Lett. 24, 775 (1970).
  * Tamagaki [1970] R. Tamagaki, Progress of Theoretical Physics 44, 905 (1970).
  * Takatsuka and Tamagaki [1971] T. Takatsuka and R. Tamagaki, Progress of Theoretical Physics 46, 114 (1971).
  * Takatsuka [1972] T. Takatsuka, Progress of Theoretical Physics 47, 1062 (1972).
  * Fujita and Tsuneto [1972] T. Fujita and T. Tsuneto, Progress of Theoretical Physics 48, 766 (1972).
  * Richardson [1972] R. W. Richardson, Phys. Rev. D 5, 1883 (1972).
  * Amundsen and Ostgaard [1985] L. Amundsen and E. Ostgaard, Nucl. Phys. A 442, 163 (1985).
  * Takatsuka and Tamagaki [1993] T. Takatsuka and R. Tamagaki, Prog. Theor. Phys. Suppl. 112, 27 (1993).
  * Baldo  _et al._ [1992] M. Baldo, J. Cugnon, A. Lejeune, and U. Lombardo, Nucl. Phys. A 536, 349 (1992).
  * Elgaroy  _et al._ [1996] O. Elgaroy, L. Engvik, M. Hjorth-Jensen, and E. Osnes, Nucl. Phys. A 607, 425 (1996), arXiv:nucl-th/9604032 [nucl-th] .
  * Khodel  _et al._ [1998] V. A. Khodel, V. V. Khodel, and J. W. Clark, Phys. Rev. Lett. 81, 3828 (1998), arXiv:nucl-th/9807034 [nucl-th] .
  * Baldo  _et al._ [1998] M. Baldo, O. Elgaroy, L. Engvik, M. Hjorth-Jensen, and H. J. Schulze, Phys. Rev. C 58, 1921 (1998), arXiv:nucl-th/9806097 [nucl-th] .
  * Khodel  _et al._ [2001] V. V. Khodel, V. A. Khodel, and J. W. Clark, Nucl. Phys. A 679, 827 (2001), arXiv:nucl-th/0001006 [nucl-th] .
  * Zverev  _et al._ [2003] M. V. Zverev, J. W. Clark, and V. A. Khodel, Nucl. Phys. A 720, 20 (2003), arXiv:nucl-th/0301028 [nucl-th] .
  * Maurizio  _et al._ [2014] S. Maurizio, J. W. Holt, and P. Finelli, Phys. Rev. C90, 044003 (2014), arXiv:1408.6281 [nucl-th] .
  * Bogner  _et al._ [2010] S. K. Bogner, R. J. Furnstahl, and A. Schwenk, Prog. Part. Nucl. Phys. 65, 94 (2010), arXiv:0912.3688 [nucl-th] .
  * Srinivas and Ramanan [2016] S. Srinivas and S. Ramanan, Phys. Rev. C 94, 064303 (2016), arXiv:1606.09053 [nucl-th] .
  * Stein  _et al._ [2016] M. Stein, A. Sedrakian, X.-G. Huang, and J. W. Clark, Phys. Rev. C 93, 015802 (2016).
  * Heinke and Ho [2010b] C. O. Heinke and W. C. G. Ho, The Astrophysical Journal Letters 719, L167 (2010b).
  * Mermin [1974] N. D. Mermin, Phys. Rev. A 9, 868 (1974).
  * Sauls and Serene [1978] J. A. Sauls and J. W. Serene, Phys. Rev. D 17, 1524 (1978).
  * Muzikar  _et al._ [1980] P. Muzikar, J. A. Sauls, and J. W. Serene, Phys. Rev. D 21, 1494 (1980).
  * Sauls  _et al._ [1982] J. A. Sauls, D. L. Stein, and J. W. Serene, Phys. Rev. D 25, 967 (1982).
  * Vulovic and Sauls [1984] V. Z. Vulovic and J. A. Sauls, Phys. Rev. D 29, 2705 (1984).
  * Masuda and Nitta [2016] K. Masuda and M. Nitta, Phys. Rev. C 93, 035804 (2016), arXiv:1512.01946 [nucl-th] .
  * Masuda and Nitta [2020] K. Masuda and M. Nitta, Prog. Theor. Exp. Phys. 2020, 013 (2020), arXiv:1602.07050 [nucl-th] .
  * Yasui  _et al._ [2019a] S. Yasui, C. Chatterjee, and M. Nitta, Phys. Rev. C 99, 035213 (2019a), arXiv:1810.04901 [nucl-th] .
  * Yasui  _et al._ [2019b] S. Yasui, C. Chatterjee, and M. Nitta,  _Proceedings, 8th International Conference on Quarks and Nuclear Physics (QNP2018): Tsukuba, Japan, November 13-17, 2018_ , JPS Conf. Proc. 26, 024022 (2019b), arXiv:1902.00674 [nucl-th] .
  * Yasui  _et al._ [2019c] S. Yasui, C. Chatterjee, M. Kobayashi, and M. Nitta, Phys. Rev. C 100, 025204 (2019c), arXiv:1904.11399 [nucl-th] .
  * Zhou and Semenoff [2006] F. Zhou and G. W. Semenoff, Phys. Rev. Lett. 97, 180411 (2006), arXiv:cond-mat/0607463 .
  * Semenoff and Zhou [2007] G. W. Semenoff and F. Zhou, Phys. Rev. Lett. 98, 100401 (2007), arXiv:cond-mat/0610162 [cond-mat] .
  * Uchino  _et al._ [2010a] S. Uchino, M. Kobayashi, and M. Ueda, Phys. Rev. A 81, 063632 (2010a).
  * Uchino  _et al._ [2010b] S. Uchino, M. Kobayashi, M. Nitta, and M. Ueda, Phys. Rev. Lett. 105, 230406 (2010b), arXiv:1010.2864 [cond-mat.quant-gas] .
  * Borgh and Ruostekoski [2016] M. O. Borgh and J. Ruostekoski, Phys. Rev. Lett. 117, 275302 (2016), [Erratum: Phys.Rev.Lett. 118, 129901(E) (2017)], arXiv:1611.09735 [cond-mat.quant-gas] .
  * Kobayashi and Nitta [2021] M. Kobayashi and M. Nitta, Phys. Rev. A 104, 053302 (2021), arXiv:2108.11676 [cond-mat.quant-gas] .
  * Mizushima  _et al._ [2021] T. Mizushima, S. Yasui, D. Inotani, and M. Nitta, Phys. Rev. C 104, 045803 (2021), arXiv:2108.01256 [nucl-th] .
  * Bedaque  _et al._ [2003] P. F. Bedaque, G. Rupak, and M. J. Savage, Phys. Rev. C 68, 065802 (2003), arXiv:nucl-th/0305032 [nucl-th] .
  * Leinson [2011a] L. B. Leinson, Phys. Lett. B 702, 422 (2011a), arXiv:1107.4025 [nucl-th] .
  * Leinson [2012] L. B. Leinson, Phys. Rev. C 85, 065502 (2012), arXiv:1206.3648 [nucl-th] .
  * Leinson [2013] L. B. Leinson, Phys. Rev. C 87, 025501 (2013), arXiv:1301.5439 [nucl-th] .
  * Bedaque and Nicholson [2013] P. F. Bedaque and A. N. Nicholson, Phys. Rev. C 87, 055807 (2013), [Erratum: Phys. Rev. C 89, 029902(E) (2014)], arXiv:1212.1122 [nucl-th] .
  * Bedaque and Sen [2014] P. Bedaque and S. Sen, Phys. Rev. C89, 035808 (2014), arXiv:1312.6632 [nucl-th] .
  * Bedaque and Reddy [2014] P. F. Bedaque and S. Reddy, Phys. Lett. B 735, 340 (2014), arXiv:1307.8183 [nucl-th] .
  * Bedaque  _et al._ [2015] P. F. Bedaque, A. N. Nicholson, and S. Sen, Phys. Rev. C 92, 035809 (2015), arXiv:1408.5145 [nucl-th] .
  * Leinson [2010a] L. B. Leinson, Phys. Rev. C 81, 025501 (2010a), arXiv:0912.2164 [astro-ph.SR] .
  * Leinson [2010b] L. B. Leinson, Phys. Lett. B 689, 60 (2010b), arXiv:1001.2617 [astro-ph.SR] .
  * Leinson [2010c] L. B. Leinson, Phys. Rev. C 82, 065503 (2010c), [Erratum: Phys. Rev.C84, 049901(E) (2011)], arXiv:1012.5387 [hep-ph] .
  * Leinson [2011b] L. B. Leinson, Phys. Rev. C 83, 055803 (2011b), arXiv:1007.2803 [hep-ph] .
  * Leinson [2011c] L. B. Leinson, Phys. Rev. C 84, 045501 (2011c), arXiv:1110.2145 [nucl-th] .
  * Yasui and Nitta [2020] S. Yasui and M. Nitta, Phys. Rev. C 101, 015207 (2020), arXiv:1907.12843 [nucl-th] .
  * Yasui  _et al._ [2020a] S. Yasui, C. Chatterjee, and M. Nitta, Phys. Rev. C 101, 025204 (2020a), arXiv:1905.13666 [nucl-th] .
  * Mizushima  _et al._ [2017] T. Mizushima, K. Masuda, and M. Nitta, Phys. Rev. B 95, 140503(R) (2017), arXiv:1607.07266 [cond-mat.supr-con] .
  * Mizushima  _et al._ [2020] T. Mizushima, S. Yasui, and M. Nitta, Phys. Rev. Research 2, 013194 (2020).
  * Masaki  _et al._ [2020] Y. Masaki, T. Mizushima, and M. Nitta, Phys. Rev. Research 2, 013193 (2020).
  * Masaki  _et al._ [2022] Y. Masaki, T. Mizushima, and M. Nitta, Phys. Rev. B 105, L220503 (2022), arXiv:2107.02448 [cond-mat.supr-con] .
  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
  * Ryu  _et al._ [2010] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New J. Phys. 12, 065010 (2010).
  * Chatterjee  _et al._ [2017] C. Chatterjee, M. Haberichter, and M. Nitta, Phys. Rev. C96, 055807 (2017), arXiv:1612.05588 [nucl-th] .
  * Kobayashi and Nitta [2022] M. Kobayashi and M. Nitta, Phys. Rev. C 105, 035807 (2022), arXiv:2203.09300 [nucl-th] .
  * Leinson [2020] L. B. Leinson, Mon. Not. Roy. Astron. Soc. 498, 304 (2020).
  * [88] G. Marmorini, S. Yasui, and M. Nitta, arXiv:2010.09032 [astro-ph.HE] .
  * Kobayashi  _et al._ [2012] S. Kobayashi, M. Kobayashi, Y. Kawaguchi, M. Nitta, and M. Ueda, Nucl. Phys. B 856, 577 (2012), arXiv:1110.1478 [math-ph] .
  * Note [1] Similar molecules of half-quantized vortices connected by a linear soliton are present in the other multicomponent condensed matter physics and high energy physics: multicomponent or multigap superconductors [112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125], coherently-coupled multicomponent BECs [126, 93, 94, 127, 128, 129, 130, 95, 131, 132, 133, 134, 135, 99, 105, 136], high density QCD [137], and two-Higgs doublet models [138]. However, the unique feature of ${}^{3}P_{2}$ superfluids is that constituent half-quantized vortices are non-Abelian vortices associated with a non-Abelian first homotopy group.
  * Note [2] A similar reconnection of solitons was observed in numerical simulations of two-component BECs [95, 99], and similar chemical bonds of vortex molecules were found in multi-component superconductors with arbitrary charges and Josephson couplings [124].
  * Salomaa and Volovik [1985] M. M. Salomaa and G. E. Volovik, Phys. Rev. B 31, 203 (1985).
  * Mueller and Ho [2002] E. J. Mueller and T.-L. Ho, Phys. Rev. Lett. 88, 180403 (2002).
  * Kasamatsu  _et al._ [2003] K. Kasamatsu, M. Tsubota, and M. Ueda, Phys. Rev. Lett. 91, 150406 (2003).
  * Cipriani and Nitta [2013] M. Cipriani and M. Nitta, Phys. Rev. Lett. 111, 170401 (2013), arXiv:1303.2592 [cond-mat.quant-gas] .
  * Kobayashi  _et al._ [2009] M. Kobayashi, Y. Kawaguchi, M. Nitta, and M. Ueda, Phys. Rev. Lett. 103, 115301 (2009), arXiv:0810.5441 [cond-mat.other] .
  * [97] M. Kobayashi, Y. Kawaguchi, and M. Ueda, arXiv:0907.3716 [cond-mat] .
  * [98] M. Kobayashi and M. Ueda, arXiv:1606.07190 [cond-mat] .
  * Eto  _et al._ [2020] M. Eto, K. Ikeno, and M. Nitta, Phys. Rev. Res. 2, 033373 (2020), arXiv:1912.09014 [cond-mat.quant-gas] .
  * Yasui  _et al._ [2020b] S. Yasui, D. Inotani, and M. Nitta, Phys. Rev. C 101, 055806 (2020b), arXiv:2002.05429 [nucl-th] .
  * Fujimoto  _et al._ [2020] Y. Fujimoto, K. Fukushima, and W. Weise, Phys. Rev. D 101, 094009 (2020), arXiv:1908.09360 [hep-ph] .
  * Fujimoto and Nitta [2021a] Y. Fujimoto and M. Nitta, Phys. Rev. D 103, 054002 (2021a), arXiv:2011.09947 [hep-ph] .
  * Fujimoto and Nitta [2021b] Y. Fujimoto and M. Nitta, Phys. Rev. D 103, 114003 (2021b), arXiv:2102.12928 [hep-ph] .
  * Fujimoto and Nitta [2021c] Y. Fujimoto and M. Nitta, Journal of High Energy Physics 2021, 192 (2021c), arXiv:2103.15185 [hep-ph] .
  * Kobayashi  _et al._ [2019] M. Kobayashi, M. Eto, and M. Nitta, Phys. Rev. Lett. 123, 075303 (2019), arXiv:1802.08763 [cond-mat.stat-mech] .
  * Schmaljohann  _et al._ [2004] H. Schmaljohann, M. Erhard, J. Kronjäger, M. Kottke, S. van Staa, L. Cacciapuoti, J. J. Arlt, K. Bongs, and K. Sengstock, Phys. Rev. Lett. 92, 040402 (2004).
  * Chang  _et al._ [2004] M.-S. Chang, C. D. Hamley, M. D. Barrett, J. A. Sauer, K. M. Fortier, W. Zhang, L. You, and M. S. Chapman, Phys. Rev. Lett. 92, 140403 (2004).
  * Kuwamoto  _et al._ [2004] T. Kuwamoto, K. Araki, T. Eno, and T. Hirano, Phys. Rev. A 69, 063604 (2004).
  * Widera  _et al._ [2006] A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and I. Bloch, New Journal of Physics 8, 152 (2006).
  * Tojo  _et al._ [2008] S. Tojo, A. Tomiyama, M. Iwata, T. Kuwamoto, and T. Hirano, Applied Physics B 93, 403 (2008).
  * Tojo  _et al._ [2009] S. Tojo, T. Hayashi, T. Tanabe, T. Hirano, Y. Kawaguchi, H. Saito, and M. Ueda, Phys. Rev. A 80, 042704 (2009).
  * Babaev [2002] E. Babaev, Phys. Rev. Lett. 89, 067001 (2002), arXiv:cond-mat/0111192 .
  * Babaev  _et al._ [2004] E. Babaev, A. Sudbo, and N. W. Ashcroft, Nature 431, 666 (2004), arXiv:cond-mat/0410408 .
  * Goryo  _et al._ [2007] J. Goryo, S. Soma, and H. Matsukawa, Europhys. Lett. 80, 17002 (2007).
  * Tanaka  _et al._ [2007] Y. Tanaka, A. Crisan, D. D. Shivagan, A. Iyo, K. Tokiwa, and T. Watanabe, Jpn. J. Appl. Phys. 46, 134 (2007).
  * Crisan  _et al._ [2007] A. Crisan, Y. Tanaka, D. D. Shivagan, A. Iyo, L. Cosereanu, K. Tokiwa, and T. Watanabe, Jpn. J. Appl. Phys. 46, L451 (2007).
  * Guikema  _et al._ [2008] J. W. Guikema, H. Bluhm, D. A. Bonn, R. Liang, W. N. Hardy, and K. A. Moler, Phys. Rev. B 77, 104515 (2008).
  * Nitta  _et al._ [2012] M. Nitta, M. Eto, T. Fujimori, and K. Ohashi, J. Phys. Soc. Jap. 81, 084711 (2012), arXiv:1011.2552 [cond-mat.supr-con] .
  * Garaud  _et al._ [2011] J. Garaud, J. Carlstrom, and E. Babaev, Phys. Rev. Lett. 107, 197001 (2011), arXiv:1107.0995 .
  * Garaud  _et al._ [2013] J. Garaud, J. Carlström, E. Babaev, and M. Speight, Phys. Rev. B 87, 014507 (2013), arXiv:1211.4342 .
  * Garaud and Babaev [2012] J. Garaud and E. Babaev, Phys. Rev. B 86, 060514(R) (2012).
  * Tanaka  _et al._ [2017] Y. Tanaka, H. Yamamori, T. Yanagisawa, T. Nishio, and S. Arisawa, Phys. C Supercond. its Appl. 538, 12 (2017).
  * Tanaka  _et al._ [2018] Y. Tanaka, H. Yamamori, T. Yanagisawa, T. Nishio, and S. Arisawa, Phys. C Supercond. its Appl. 548, 44 (2018).
  * Chatterjee  _et al._ [2020] C. Chatterjee, S. B. Gudnason, and M. Nitta, Journal of High Energy Physics 2020, 109 (2020), arXiv:1912.02685 [hep-th] .
  * [125] J. Garaud, arXiv:2203.12194[cond-mat.supr-con] .
  * Son and Stephanov [2002] D. T. Son and M. A. Stephanov, Phys. Rev. A 65, 063621 (2002), arXiv:cond-mat/0103451 .
  * Kasamatsu  _et al._ [2004] K. Kasamatsu, M. Tsubota, and M. Ueda, Phys. Rev. Lett. 93, 250406 (2004), arXiv:cond-mat/0406150 .
  * Kuopanportti  _et al._ [2012] P. Kuopanportti, J. A. M. Huhtamäki, and M. Möttönen, Phys. Rev. A 85, 043613 (2012).
  * Aftalion  _et al._ [2012] A. Aftalion, P. Mason, and J. Wei, Phys. Rev. A 85, 033614 (2012).
  * Eto and Nitta [2012] M. Eto and M. Nitta, Phys. Rev. A 85, 053645 (2012), arXiv:1201.0343 [cond-mat.quant-gas] .
  * Eto and Nitta [2013] M. Eto and M. Nitta, EPL 103, 60006 (2013), arXiv:1303.6048 [cond-mat.quant-gas] .
  * Nitta  _et al._ [2013] M. Nitta, M. Eto, and M. Cipriani, J. Low Temp. Phys. 175, 177 (2013), arXiv:1307.4312 .
  * Dantas  _et al._ [2015] D. S. Dantas, A. R. P. Lima, A. Chaves, C. A. S. Almeida, G. A. Farias, and M. V. Milošević, Phys. Rev. A 91, 023630 (2015).
  * Tylutki  _et al._ [2016] M. Tylutki, L. P. Pitaevskii, A. Recati, and S. Stringari, Phys. Rev. A 93, 043623 (2016), arXiv:1601.03695 [cond-mat.quant-gas] .
  * Eto and Nitta [2018] M. Eto and M. Nitta, Phys. Rev. A 97, 023613 (2018), arXiv:1702.04892 [cond-mat.quant-gas] .
  * Mencia Uranga and Lamacraft [2018] B. Mencia Uranga and A. Lamacraft, Phys. Rev. A 97, 043609 (2018).
  * Eto and Nitta [2021] M. Eto and M. Nitta, Phys. Rev. D 104, 094052 (2021), arXiv:2103.13011 [hep-ph] .
  * Eto  _et al._ [2022] M. Eto, Y. Hamada, and M. Nitta, JHEP 02, 099, arXiv:2111.13345 [hep-ph] .