#  Core structures of vortices in Ginzburg-Landau theory for neutron
${}^{3}P_{2}$ superfluids

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

We investigate vortex solutions in the Ginzburg-Landau theory for neutron
${}^{3}P_{2}$ superfluids relevant for neutron star cores in which neutron
pairs possess the total angular momentum $J=2$ with spin-triplet and $P$ wave,
in the presence of the magnetic field parallel to the angular momentum of
vortices. The ground state is known to be in the uniaxial nematic (UN) phase
in the absence of magnetic field, while it is in the $D_{2}$ ($D_{4}$) biaxial
nematic (BN) phase in the presence of the magnetic field below (above) the
critical value. We find that a singly quantized vortex always splits into two
half-quantized non-Abelian vortices connected by soliton(s) as a vortex
molecule with any strength of the magnetic field. In the UN phase, two half-
quantized vortices with ferromagnetic cores are connected by a linear soliton
with the $D_{4}$ BN order. In the $D_{2}$ ($D_{4}$) BN phase, two half-
quantized vortices with cyclic cores are connected by three linear solitons
with the $D_{4}$ ($D_{2}$) BN order. The energy of the vortex molecule
monotonically increases and the distance between the two half-quantized
vortices decreases with the magnetic field increases, except for a
discontinuously increasing jump of the distance at the critical magnetic
field. We also construct an isolated half-quantized non-Abelian vortex in the
$D_{4}$ BN phase.

##  I Introduction

Pulsars or rapidly rotating neutron stars are dense and compact stars under
extreme conditions, thereby serving as astrophysical laboratories for studying
nuclear and QCD matter at high density, with a strong magnetic field and under
rapid rotation [1, 2]. The recent progresses in astrophysical observations
promote us to study the interiors of neutron stars more precisely: the
observation of massive neutron stars whose masses are about twice as large as
the solar mass [3, 4], the gravitational waves from a binary neutron star
merger [5, 6], and the Neutron Star Interior Composition Explorer (NICER)
mission [7, 8], expected to reveal interior states of neutron stars.

The interior of neutron stars is believed to exhibit neutron superfluidity and
proton superconductivity as first predicted by Migdal [9] (see Refs. [10, 11,
12, 13, 1, 14] for recent reviews). Such superfluid and superconducting
components can alter low-energy excitation modes compared with the normal
phase, and thus their existence can affect several processes and properties of
neutron stars, such as neutrino emissivities and specific heats relevant to
the long relaxation time after the sudden speed-up events, that is pulsar
glitches, 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]. The neutron superfluids are realized by the attraction between
two neutrons in the ${}^{1}S_{0}$ channel at the low density [9]. It was,
however, shown in Ref. [24] that this channel is repulsive at higher densities
as a consequence of the strong short-range repulsion. Thus, it was proposed
that neutron ${}^{3}P_{2}$ superfluids, in which neutron pairs possess the
total angular momentum $J=2$ with spin triplet and $P$ wave, are more relevant
at higher density [25, 26, 27, 26, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42]. The ${}^{3}P_{2}$ interaction originates from a strong spin-
orbit force at large scattering energy, and thus the neutron ${}^{3}P_{2}$
superfluids are expected to be realized in the high-density regions in the
inner cores of neutron stars. They also can survive in neutron stars with
strong magnetic fields, such as magnetars with the magnetic field
$10^{15}$–$10^{18}$ G, since they are tolerant against the strong magnetic
field due to the fact that the aligned pairs of Cooper pairs with the spin-
triplet pairing are not broken by the Zeeman effect. In the $S$-wave case, it
has also predicted that Cooper pairs can survive at around the magnetic field
$10^{17}$ G [43].

In astrophysical observations, the possibility of the existence of neutron
${}^{3}P_{2}$ superfluids inside neutron stars are investigated; the rapid
cooling of the neutron star in Cassiopeia A might be explained by the
enhancement of neutrino emissivities due to the formation and dissociation of
neutron ${}^{3}P_{2}$ Cooper pairs [21, 22, 23].

On the other hand, as the theoretical aspects are concerned, neutron
${}^{3}P_{2}$ superfluids have rich topological structures both in bosonic and
fermionic excitations. There are basically two related approaches for
theoretical study of superfluids. The most fundamental theory is the
Bogoliubov-de Gennes (BdG) equation offering a framework to deal with fermion
degrees of freedom. The other is the Ginzburg-Landau (GL) approach as the low-
energy effective theory obtained by integrating out fermion degrees of
freedom, which is an expansion of both the order parameters and spatial
derivatives, and thus is valid only in the region close to the critical
temperature. The BdG approach was applied to ${}^{3}P_{2}$ superfluids and the
phase diagram in the plane of the temperature and magnetic field was obtained
[44]. Furthermore, ${}^{3}P_{2}$ superfluids were shown to be topological
superfluids of a class DIII in the classification of topological insulators
and superconductors [45, 46], allowing a topologically protected gapless
Majorana fermion on its boundary [44]. On the other hand, within the GL
theory, superfluid states with $J=2$ are in general classified into nematic,
cyclic, and ferromagnetic phases etc [47]. The GL theory for ${}^{3}P_{2}$
superfluids was obtained [30, 31, 48, 49, 50, 51, 52, 53, 54, 55, 56], and in
the weak coupling limit, the nematic phase was found to be the ground state of
${}^{3}P_{2}$ superfluids [48, 49, 50]. The nematic phase consists of three
subphases 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, where $D_{n}$
is a dihedral group of order $n$ [see Table 1(a) and 1(b)]. Corresponding
order parameter manifolds (OPMs) are $U(1)\times SO(3)/O(2)\simeq
S^{1}\times{\mathbb{R}}P^{2}$, $U(1)\times SO(3)/D_{2}$, and $[U(1)\times
SO(3)]/D_{4}$, respectively [see Table 1(c)]. These are continuously
degenerated in the absence of magnetic field in the GL expansion up to the
fourth order 111 More precisely these are connected by a parameter of
continuous degeneracy called a quasi-Nambu-Goldstone mode [57], and these OPMs
are submanifolds of an extended OPM $(S^{1}\times S^{4})/{\mathbb{Z}}_{2}$. .
In the presence of the magnetic field and/or with the inclusion of the sixth-
order term into the GL theory, the continuous degeneracy is lifted to pick up
either UN, $D_{2}$ BN, or $D_{4}$ BN state as the ground state for zero
magnetic field, nonzero one below the critical value $B_{c}$, and nonzero one
above $B_{c}$, respectively [52, 54, 55] [see Table 1(a)]. There is a subtle
problem on the instability of the ground states for large value of the order
parameter, which is cured by the expansion up to the eighth order [56]. The
phase diagram up to the eighth order captures the essential features of that
determined in the BdG equation [44], including a tricritical point connecting
first- and second-order phase transition lines between $D_{4}$ and $D_{2}$ BN
phases [44, 58]. Apart from nematic phases, more general uniform states of
${}^{3}P_{2}$ superfluids were classified [59], which is also useful to
identify local states such as vortex cores. As a uniform ground state, the
ferromagnetic state is in fact found to appear, beyond the quasi-classical
approximation, in the region close to the critical temperature [60]. The GL
approach is useful to deal with bosonic collective excitations and topological
defects. Bosonic excitations in the ${}^{3}P_{2}$ superfluids yield collective
modes [61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73] relevant, for
instance, for cooling process of neutron stars. Topological defects such as
domain walls [74] and the boundary defect (boojums) of ${}^{3}P_{2}$
superfluids [75] were investigated.

One of the most salient features of superfluidity is the fact that
circulations of vortices are quantized so that the wave function is single
valued around the vortices (the Feynman-Onsager’s quantization), yielding the
existence of quantized vortices. When a superfluid is rotating, a vortex
lattice is formed as observed in helium superfluids and ultracold atomic
gasses. In the context of superfluids in neutron stars, the origin of pulsar
glitches was proposed to be explained by sudden releases of a large number of
quantized vortices [76, 77]. In the case of ${}^{3}P_{2}$ superfluids,
quantized vortices were investigated both in the GL theory [49, 50, 30, 52,
78, 53] (see also Ref. [79] for coreless vortices), and in the BdG theory [80,
81]. The first homotopy group classifying vortices is given in Ref. [52] [see
Table 1(d)]. Singly quantized vortices in ${}^{3}P_{2}$ superfluids were
studied in the GL theory [49, 50, 30, 52, 78], and topologically protected
Majorana fermions in the vortex core were found in the BdG theory [80].
Vortices more peculiar to the ${}^{3}P_{2}$ superfluids are half-quantized
non-Abelian vortices which are allowed only in the $D_{4}$ BN phase [53, 81].
Their circulations are a half of the Feynman-Onsager’s quantized circulations,
and the first homotopy group characterizing these vortices is non-Abelian,
thus giving non-commutativity of exchanging vortices. The existence of half-
quantized vortices was proposed to explain a scaling law of pulsar glitches
[82]. In Ref. [81], it was found in the BdG formalism that a singly quantized
vortex is split into two half-quantized non-Abelian vortices. It was also
found that a Majorna fermion zero mode is trapped in each half-quantized
vortex.

Apart from ${}^{3}P_{2}$ superfluids, spin-2 spinor ultracold atomic Bose-
Einstein condensates (BECs) are also $J=2$ condensates whose ground states are
possibly nematic phase [83, 84, 85, 57, 86, 59] sharing almost the same
bosonic properties with ${}^{3}P_{2}$ superfluids, thus admitting the same
order parameter manifold and non-Abelian half-quantized vortices [57, 86].
Therefore, studying bosonic properties of ${}^{3}P_{2}$ superfluids is also
applicable to spin-2 BECs, which can be experimentally testable in principle,
although the current experiments of 87Rb atoms imply their ground state to be
in the cyclic or nematic phase [87, 88, 89, 90, 91, 92].

$|\boldsymbol{B}|$ | $|\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
(f) $\\#$ of solitons | 1 | 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 value $B_{c}$, (a) phases, (b) unbroken symmetries, (c) OPMs, (d)
the first homotopy groups of the OPMs, (e) the order inside half-quantized
vortex cores, (f) the number of solitons connecting the two half-quantized
vortices, and (g) the order inside 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. [93], ensuring the
existence of isolated half-quantized non-Abelian vortices. (e), (f), and (g)
are new results presented in this paper.

In this paper, we investigate vortex solutions, namely singly quantized
vortices and half-quantized non-Abelian vortices, in neutron ${}^{3}P_{2}$
superfluids within the GL approach. The orientation of the magnetic field is
fixed to be parallel to the angular momentum of vortices. In the previous
studies of vortices in the GL theory, an axial symmetry around the vortex axis
was assumed [49, 50, 30, 52, 78, 53]. In contrast, imposing no axial symmetry,
we find that a singly quantized vortex always splits into two half-quantized
non-Abelian vortices with any strength of the magnetic field. An advantage to
use the GL theory compared with the BdG equation employed in Ref. [80] is that
we do not have to consider the direction of splitting a priori. In the UN
phase with the zero magnetic field, cores of two half-quantized vortices are
found to be filled with the ferromagnetic states, while they are filled with
the cyclic states 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, the two half-
quantized vortices are connected by a single soliton of the $D_{4}$ BN order
while they are connected by three linear solitons with the $D_{4}$ ($D_{2}$)
BN order in the $D_{2}$ ($D_{4}$) BN phase, as summarized in Table 1(f) and
1(g). The appearance of the $D_{4}$ BN order around the vortex core in the UN
and $D_{2}$ BN phases is a consequence of the fact that isolated half-
quantized vortices can topologically exist only in the $D_{4}$ BN state. We
also show that the energy of the vortex molecule monotonically increases as
the magnetic field increases, which is continuous everywhere including the
critical magnetic field separating $D_{2}$ and $D_{4}$ BN states. On the other
hand, the distance between the two half-quantized vortices decreases with the
magnetic field increases, except for a discontinuously increasing jump at the
critical magnetic field. We also construct an isolated half-quantized non-
Abelian vortex in the $D_{4}$ BN phase above the critical magnetic field.

A molecule of half-quantized vortices connected by a soliton or domain wall
can be found in various systems such as multicomponent or multigap
superconductors [94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106],
coherently coupled multicomponent BECs [107, 108, 109, 110, 111, 112, 113,
114, 115, 116, 117, 118, 119, 120, 121, 122], dense QCD of quark matter [123],
and the two-Higgs doublet model [124] as a model beyond the standard model of
elementary particles. Compared with these systems, the unique feature of
${}^{3}P_{2}$ superfluids is that constituent half-quantized vortices are non-
Abelian vortices, that is, characterized by a non-Abelian first homotopy
group.

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. Section III
shows our numerical results for vortex states in the ${}^{3}P_{2}$
superfluids. Section IV is devoted to a summary and discussion.

##  II Ginzburg-Landau free energy 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. [59]. 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}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},\\\
&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),\\\
&f_{002}=3\rho,\\\
&f_{004}=6\rho^{2}+\frac{3}{4}\boldsymbol{S}^{2}-\frac{3}{2}|\Psi_{20}|^{2},\\\
&f_{006}=-324\rho^{3}-81\rho\boldsymbol{S}^{2}+162\rho|\Psi_{20}|^{2}+15|\Psi_{30}|^{2}-27|\Phi_{30}|^{2},\\\
&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}\\\
&f_{022}=2\rho\boldsymbol{B}^{2}-\frac{1}{2}\psi^{\dagger}\hat{S}_{\boldsymbol{B}}\hat{S}_{\boldsymbol{B}}\psi,\\\
&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),\end{aligned}$
(2)

where $\hat{S}_{i}$ ($i=x,y,z$) are $5\times 5$ spin-2 matrices, $\hat{T}$ is
the time reversal operator defined by
$(\hat{T}\psi)_{m}\equiv(-1)^{m}\psi_{-m}$, and
$\hat{S}_{\boldsymbol{B}}\equiv\boldsymbol{\hat{S}}\cdot\boldsymbol{B}$. The
$\psi$-dependent terms $\boldsymbol{j}$, $\rho$, $\boldsymbol{S}$,
$\Psi_{20}$, $\Psi_{30}$, $\Phi_{30}$, $\Gamma_{4}$ are defined by

$\displaystyle\begin{aligned}
&\boldsymbol{j}=-i\nabla\psi,\quad\rho=\psi^{\dagger}\psi,\quad\boldsymbol{S}=\psi^{\dagger}\hat{\boldsymbol{S}}\psi,\\\
&\Psi_{20}=\psi^{T}\hat{T}\psi=\sqrt{5}\sum_{m_{1},m_{2}=-2}^{2}C^{00}_{2m_{1},2m_{2}}\psi_{m_{1}}\psi_{m_{2}},\\\
&\Psi_{30}=-\sqrt{\frac{35}{2}}\sum_{J=0}^{4}\sum_{M=-J}^{J}\sum_{m_{1},m_{2},m_{3}=-2}^{2}\\\
&\phantom{\Psi_{30}}\times
C^{00}_{JM,2m_{3}}C^{JM}_{2m_{1},2m_{2}}\psi_{m_{1}}\psi_{m_{2}}\psi_{m_{3}},\\\
&\Phi_{30}=-\sqrt{\frac{35}{2}}\sum_{J=0}^{4}\sum_{M=-J}^{J}\sum_{m_{1},m_{2},m_{3}=-2}^{2}\\\
&\phantom{\Phi_{30}}\times
C^{00}_{JM,2m_{3}}C^{JM}_{2m_{1},2m_{2}}\psi_{m_{1}}\psi_{m_{2}}\varphi^{\ast}_{m_{3}},\\\
&\Gamma_{4}=\mathrm{Re}\left[\Psi_{20}\Phi_{30}^{\ast 2}\right].\end{aligned}$
(3)

The GL coefficients can be obtained in the weak coupling limit within the
quasiclassical approximation starting from the nonrelativistic 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}},$ (4)
$\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}$.

The spin-2 spinor order parameter is often written by the $3\times 3$
traceless symmetric matrix $A$ given by

$\displaystyle\begin{aligned}
&[A]_{11}=\frac{\sqrt{3}}{2}(\psi_{2}+\psi_{-2})-\frac{1}{\sqrt{2}}\psi_{0},\\\
&[A]_{12}=[A]_{21}=\frac{\sqrt{3}i}{2}(\psi_{2}-\psi_{-2}),\\\
&[A]_{13}=[A]_{31}=-\frac{\sqrt{3}}{2}(\psi_{1}-\psi_{-1}),\\\
&[A]_{22}=-\frac{\sqrt{3}}{2}(\psi_{2}+\psi_{-2})-\frac{1}{\sqrt{2}}\psi_{0},\\\
&[A]_{23}=[A]_{32}=-\frac{\sqrt{3}i}{2}(\psi_{1}+\psi_{-1}),\\\
&[A]_{33}=\sqrt{2}\psi_{0}.\end{aligned}$ (5)

All candidates for uniform ground states were classified in Ref. [59], and
characterized by $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$,
$|\Psi_{20}|^{2}$, and $|\Psi_{30}|^{2}$. The five characteristic symmetric
states are ferromagnetic ($\boldsymbol{S}^{2}/\rho^{2}=4\rho^{2}$,
$|\Psi_{20}|^{2}=|\Psi_{30}|^{2}=0$), uniaxial nematic (UN)
($\boldsymbol{S}^{2}=0$, $|\Psi_{20}|^{2}=\rho^{2}$,
$|\Psi_{30}|^{2}=\rho^{3}$), $D_{4}$ biaxial nematic (BN)
($\boldsymbol{S}^{2}=0$, $|\Psi_{20}|^{2}=\rho^{2}$, $|\Psi_{30}|^{2}=0$),
$D_{2}$ BN ($\boldsymbol{S}^{2}=0$, $|\Psi_{20}|^{2}=\rho^{2}$,
$0<|\Psi_{30}|^{2}<\rho^{3}$), and cyclic ($\boldsymbol{S}^{2}=0$,
$|\Psi_{20}|^{2}=0$, $|\Psi_{30}|^{2}=2\rho^{3}$) states.

For the effective Lagrangian density $f$ in Eq. (1), the UN, $D_{2}$ BN, and
$D_{4}$ BN states are predicted to be realized at $|\boldsymbol{B}|=0$,
$0<|\boldsymbol{B}|<B_{\rm c}$, and $|\boldsymbol{B}>B_{c}$, in ${}^{3}P_{2}$
superfluids [44] 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$.

##  III Vortex solutions

###  III.1 Ansatz

Next, we consider vortex solutions with vortex cores placed at $r=0$ in the
cylindrical coordinate $(r,\theta,z)$ and the boundary $\psi|_{r\to\infty}$
far from vortex cores. The vortices and the angular momentum are parallel to
the $z$ axis. For the so-called integer vortices, the order parameters behave
as $\psi_{-2\leq m\leq 2}|_{r\to\infty}\propto e^{i\theta}$ around which the
overall phase of $\psi$ winds by $2\pi$. In this section, we determine the
boundary conditions around the vortices for the cases of $\boldsymbol{B}=0$
and $\boldsymbol{B}\neq 0$,

For $\boldsymbol{B}=0$, the uniform ground state is degenerate within the
possible UN state

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

or

$\displaystyle\begin{aligned}
&\tilde{A}=\frac{1}{\sqrt{2}}R_{z}^{a}R_{y}^{b}\begin{pmatrix}-1&0&0\\\
0&-1&0\\\ 0&0&2\end{pmatrix}R_{y}^{-b}R_{z}^{-a},\\\
&R_{y}^{b}\equiv\begin{pmatrix}\cos b&0&\sin b\\\ 0&1&0\\\ -\sin b&0&\cos
b\end{pmatrix},\\\ &R_{z}^{a}\equiv\begin{pmatrix}\cos a&-\sin a&0\\\ \sin
a&\cos a&0\\\ 0&0&1\end{pmatrix},\end{aligned}$ (7)

for $0\leq a,2b,\phi\leq 2\pi$. Here $\tilde{\psi}$ is defined as
$\tilde{\psi}\equiv\psi/\sqrt{\rho}$, and $\tilde{A}$ is defined by replacing
$A$ and $\psi$ with $\tilde{A}$ and $\tilde{\psi}$, respectively in Eq. (5).
$a$, $b$, and $\phi$ represent overall spin rotations along $z$ axis and $y$
axis, and overall phase shift, respectively. Under the spatial phase gradient
$e^{i\theta}$ for the vortex solution, however, the current-spin-dependent
free energy density $f_{202}^{(2)}$ in Eq. (1) favors $b=0$, giving

$\displaystyle\tilde{\psi}\xrightarrow{r\to\infty}\left(0,0,e^{i\theta},0,0\right)^{T},$
(8)

or

$\displaystyle\tilde{A}\xrightarrow{r\to\infty}\frac{e^{i\theta}}{\sqrt{2}}\begin{pmatrix}-1&0&0\\\
0&-1&0\\\ 0&0&2\end{pmatrix},$ (9)

We next consider the case of $|\boldsymbol{B}|>0$. In this paper, the simplest
situation $\boldsymbol{B}\parallel\boldsymbol{z}$ is considered, i.e., the
magnetic field is parallel to the angular momentum for the vortex. For
$0<|\boldsymbol{B}|<B_{\rm c}$, field-dependent free energy density
$\beta_{2}f_{022}+\gamma_{2}f_{024}$ favors $\psi_{\pm 2}$, leading

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

or

$\displaystyle\begin{aligned}
&\tilde{A}\xrightarrow{r\to\infty}\sqrt{2}e^{i\theta}R_{z}^{a}\begin{pmatrix}\sin
g_{-}&0&0\\\ 0&-\sin g_{+}&0\\\ 0&0&\cos g\end{pmatrix}R_{z}^{-a},\\\
&g_{\pm}\equiv g\pm\frac{\pi}{6},\end{aligned}$ (11)

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 $z$ axis as well as that in Eq. (6). Without
current-spin dependent free energy $f_{202}^{(1)}$, $a$ takes the arbitrary
value, but is fixed with the finite $f_{202}^{(1)}$ to minimize this. In the
limit of $|\boldsymbol{B}|\searrow 0$, $g$ becomes $g\to\pi/3$ giving

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

or

$\displaystyle\tilde{A}\xrightarrow{r\to\infty}\frac{e^{i\theta}}{\sqrt{2}}R_{z}^{a}\begin{pmatrix}1&0&0\\\
0&-2&0\\\ 0&0&1\end{pmatrix}R_{z}^{-a}.$ (13)

This solution belongs to the uniaxial nematic state in Eq. (6) with $b=\pi/2$,
but is different from that for $\boldsymbol{B}=0$ shown in Eq. (8), which
leads the discontinuity between $\boldsymbol{B}=0$ and
$|\boldsymbol{B}|\searrow 0$.

In the limit of $|\boldsymbol{B}|\nearrow B_{\rm c}$, $g$ becomes $g\to\pi/2$,
giving

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

or

$\displaystyle\tilde{A}\xrightarrow{r\to\infty}\sqrt{\frac{3}{2}}e^{i\theta}R_{z}^{a}\begin{pmatrix}1&0&0\\\
0&-1&0\\\ 0&0&0\end{pmatrix}R_{z}^{-a}.$ (15)

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

For $|\boldsymbol{B}|\geq B_{\rm c}$, the vortex state is the same as that in
Eq. (14) that belongs to $D_{4}$ BN state.

For the half-quantized vortex, the overall phase of $\psi$ winds by $\pi$. To
keep the single-valued property of the order parameter, the spin also rotates.
Only the $D_{4}$ BN state under $|\boldsymbol{B}|>B_{\rm c}$ enables the half-
quantized vortex with the order parameter giving

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

or

$\displaystyle\begin{aligned}
&\tilde{A}\xrightarrow{r\to\infty}\frac{\sqrt{3}}{2\sqrt{2}}R_{z}^{a}\begin{pmatrix}\theta_{+}&i\theta_{-}&0\\\
i\theta_{-}&-\theta_{+}&0\\\ 0&0&0\end{pmatrix}R_{z}^{-a},\\\
&\theta_{\pm}\equiv e^{i\theta}\pm 1\end{aligned}$ (17)

in the case of $\boldsymbol{B}\parallel\boldsymbol{z}$. For the vortex
solution (16), the $z$ component of the spin rotates by $\pi/2$ around the
vortex. The other solution with the spin rotation by $-\pi/2$ is

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

or

$\displaystyle\tilde{A}\xrightarrow{r\to\infty}\frac{\sqrt{3}}{2\sqrt{2}}R_{z}^{a}\begin{pmatrix}\theta_{+}&-i\theta_{-}&0\\\
-i\theta_{-}&-\theta_{+}&0\\\ 0&0&0\end{pmatrix}R_{z}^{-a}.$ (19)

###  III.2 Numerical results

Figure 1:  The vortex state in the UN phase at the zero magnetic field
$\boldsymbol{B}=0$. 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 1.33\times 10^{4}$ fm. The two half-quantized vortices are
connected by a single line soliton of the $D_{4}$ BN order.

Figure 2:  The vortex state in the $D_{2}$ BN phase at the magnetic field
$\boldsymbol{B}=0.5B_{\rm c}\hat{\boldsymbol{z}}$. 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 1.33\times 10^{4}$ fm. The two half-
quantized vortices are connected by three line solitons with the $D_{4}$ BN
order in their cores.

In this subsection, we show the numerical results for the overall vortex state
by minimizing the free-energy $f$ under the boundary condition

$\displaystyle\psi_{m}(\theta+\pi/2)=i\psi_{m}(\theta),$ (20)

or

$\displaystyle[A]_{ij}(\theta+\pi/2)=i[A]_{ij}(\theta),$ (21)

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

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

The solution of Eq. (22) can be obtained by introducing the relaxation time
$t$ and the dependence of the order parameter $\psi_{m}$ on $t$, and solving

$\displaystyle\dot{\psi}_{m}=-\frac{\delta f}{\delta\psi_{m}^{\ast}}.$ (23)

After the long time evolution of Eq. (23), we attain the solution of Eq. (22).
With an appropriate scaling of the time $t$, Eq. (23) is nothing but the time-
dependent GL equation that is often used in the research field of
superconductivity. However, it has not yet derived from the microscopic theory
for the ${}^{3}P_{2}$ superfluids.

Here, $L$ is the system size fixed to be $L=128p_{\rm F}/[\pi^{2}N(0)T_{\rm
c}^{2}]$. The temperature is fixed to be $T=0.854T_{\rm c}\approx 0.171$ MeV
for which $B_{\rm c}$ takes the maximal value $B_{\rm c}=7.06\times
10^{-2}\pi(1+F_{0}^{a})T_{\rm c}/\gamma_{\rm n}\approx 7.36\times 10^{15}$ G
with the critical temperature $T_{\rm c}\approx 0.2$ MeV and Landau parameter
$F_{0}^{a}\approx 1$. The magnetic field
$\boldsymbol{B}(\parallel\boldsymbol{z})$ changes from 0 to $1.5B_{\rm c}$.

####  III.2.1 The vortex state with ferromagnetic core in the UN phase

We start from the case of the UN phase at zero magnetic field. Figure 1 shows
the squared modulus $|\psi_{m}|^{2}$, argument $\mathrm{Arg}[\psi_{m}]$ of the
order parameter, $U(1)\times SO(3)$ invariants $\boldsymbol{S}^{2}$,
$|\Psi_{20}|^{2}$, and $|\Psi_{30}|^{2}$, and the free energy density $f$ at
the zero magnetic field $\boldsymbol{B}=0$. The radius of circles in figures
is fixed to be $64\hbar^{2}k_{\rm F}/(\pi M_{\rm n}^{\ast}T_{\rm c})\approx
1.33\times 10^{4}$ fm. Where $k_{\rm F}=(3\pi^{2}n)^{1/3}\approx 2.20$ fm-1 is
the neutron Fermi wave number with the neutron number density $n\approx
2.25n_{0}\approx 0.36$ fm-3 for the saturation density $n_{0}\approx 0.16$
fm-3 of nuclear matter. $M_{\rm n}^{\ast}$ is the effective neutron mass
$M_{\rm n}^{\ast}c^{2}\approx 0.7M_{\rm n}c^{2}\approx 658$ MeV for the
neutron vacuum mass $M_{\rm n}c^{2}\approx 940$ MeV. The critical temperature
is set to be $T_{\rm c}\approx 0.2$ MeV. The order parameter $\psi$ satisfies
Eq. (8) near the edge of the system, where the state belongs to the uniaxial
nematic state $\boldsymbol{S}^{2}=0$,
$|\Psi_{20}|^{2}/\rho^{2}=|\psi_{30}|^{2}/\rho^{3}=1$. At the center of the
system, there are two holes of $|\Psi_{20}|^{2}$ implying the breakdown of the
UN order. These two holes correspond to the vortex cores and each of them
carries half circulations. This result clearly shows that a singly quantized
vortex splits into two half-quantized non-Abelian vortices around the vortex
core. 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. 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 along a line structure bridging two
vortex cores, as can be seen in the plot of $|\Psi_{30}|^{2}$ locally inducing
the $D_{4}$ BN order and half-quantized vortices. On the other hand, at the
vortex cores, the $U(1)\times SO(3)$ invariants are
$\boldsymbol{S}^{2}/\rho^{2}\simeq 4$ and
$|\Psi_{20}|^{2}/\rho^{2}=|\psi_{30}|^{2}/\rho^{3}=0$, implying the appearance
of the ferromagnetic order. Therefore, we characterize this vortex by the
ferromagnetic core.

In this phase, there is also a metastable vortex molecule state with the
cyclic cores similar to the $D_{2}$ phase discussed below.

####  III.2.2 The vortex state with the cyclic core in the $D_{2}$ BN phase

At $|\boldsymbol{B}|=0.5B_{\rm c}$ as shown in Fig. 2, the order parameter
drastically changes from that in Eq. (8) to that in Eq. (10) where $\psi_{\pm
2}$ become finite at $r\to\infty$. As well as the case for $\boldsymbol{B}=0$,
a singly quantized vortex splits into two half-quantized vortices with two
holes of $|\Psi_{20}|^{2}$. In contrast to the vortex molecule in the UN
phase, there are three soliton lines of $D_{4}$ BN order with
$|\Psi_{30}|^{2}/\rho^{3}=0$ bridging two half-quantized vortices. At the
vortex cores, we have $\boldsymbol{S}^{2}/\rho^{2}\simeq 0$ and
$|\psi_{30}|^{2}/\rho^{3}\simeq 2$ supporting the cyclic order. We
characterize this vortex by the cyclic core. With turning off the magnetic
field, it becomes a metastable state with higher energy than the lowest energy
state of the vortex molecule discussed in Sec III III.2 III.2.1.

Figure 3:  The vortex state in the $D_{4}$ BN phase at the magnetic field
$\boldsymbol{B}=1.3B_{\rm c}\hat{\boldsymbol{z}}$. 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$ are
shown. The radius of the figure shown here is $64p_{\rm F}/(\pi m_{\rm
n}T_{\rm c})\approx 1.33\times 10^{4}$ fm. The two half-quantized vortices are
connected by three line solitons of the $D_{2}$ BN order in their cores.

Figure 4:  The vortex state at the boundary between the $D_{2}$ BN and $D_{4}$
BN phases at the magnetic field $\boldsymbol{B}=B_{\rm
c}\hat{\boldsymbol{z}}$. 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$ are shown. The radius of
the figure shown here is $64p_{\rm F}/(\pi m_{\rm n}T_{\rm c})\approx
1.33\times 10^{4}$ fm. The two half-quantized non-Abelian vortices are
connected by three line solitons asymmetrically.

####  III.2.3 The vortex state with the cyclic core in the $D_{4}$ BN phase

At $|\boldsymbol{B}|>B_{\rm c}$, the order parameter in the bulk becomes the
$D_{4}$ BN state shown in Eq. (14) at $r\to\infty$. Figure 4 shows the $SO(3)$
invariants and the free-energy density at the magnetic field
$|\boldsymbol{B}|=1.3B_{\rm c}$. Although one half-quantized vortex is
topologically stable in this state, two half-quantized vortices form a bound
state to be an integer vortex. The cores of half-quantized vortices are filled
with the cyclic order having $\boldsymbol{S}^{2}/\rho^{2}\simeq 0$ and
$|\psi_{30}|^{2}/\rho^{3}\simeq 2$ as well as the case of the $D_{2}$ BN phase
($0<|\boldsymbol{B}|<B_{\rm c}$). In contrast to the case of the $D_{2}$ BN
phase, the three line solitons bridging two half-quantized vortices are
characterized by the $D_{2}$ BN order with $0<|\Psi_{30}|^{2}/\rho^{3}<1$.

Our results suggest that the order of three line solitons connecting the two
half-quantized vortices are exchanged between $D_{2}$ and $D_{4}$ BN orders
for $|\boldsymbol{B}|<B_{\rm c}$ and $|\boldsymbol{B}|>B_{\rm c}$,
respectively. Another characteristic feature of the case of
$|\boldsymbol{B}|>B_{\rm c}$ is a fact that $\psi_{\pm 1}$ components
completely vanish.

####  III.2.4 The vortex state at the boundary between the $D_{2}$ and $D_{4}$
BN phases

Figure 4 shows the vortex state at the critical magnetic field
$|\boldsymbol{B}|=B_{\rm c}$. Especially, we can see the asymmetric shape of
$|\Psi_{30}|^{2}/\rho^{3}$ as an intermediate state between those below and
above the critical magnetic field $B_{\rm c}$, in which two $D_{4}$ and
$D_{2}$ BN orders compete as candidates of the line solitons connecting the
two half-quantized vortices. This asymmetric structure of the vortex core soon
vanishes as the magnetic field $\boldsymbol{B}$ becomes away from the critical
magnetic field $B_{\rm c}$. In our case, the vortex core becomes symmetric at
$|\boldsymbol{B}|=0.9B_{\rm c}$ and $|\boldsymbol{B}|=1.1B_{\rm c}$.

####  III.2.5 Energetics and distance between half-quantized vortices

Figure 5:  (a): Free energy $F$ and (b): distance $\Delta$ between two half-
quantized non-Abelian vortices as a function of the magnetic field
$|\boldsymbol{B}|$.

Figure 5 (a) shows the free energy $F=\int drd\theta\>rf$ as a function of the
magnetic field $|\boldsymbol{B}|$. The free energy $F$ monotonically increases
with magnetic field, and it is continuous when the magnetic field across the
critical one $B_{\rm c}$. At the zero magnetic field $\boldsymbol{B}=0$, the
vortex solution with ferromagnetic cores with the boundary condition in Eq.
(8) has the lower free energy $F$ than that with the cyclic cores and the
boundary shown in Eq. (10), and discontinuously connects at
$\boldsymbol{B}=0$. However, a vortex solution with cyclic cores with the
boundary in Eq. (10) can also exist as the metastable solution, continuously
connecting to the solutions in the $D_{2}$ BN phase $|\boldsymbol{B}|>0$.

Figure 5 (b) shows the distance $\Delta$ between two half-quantized vortices
as a function of the magnetic field $|\boldsymbol{B}|$. $\Delta$ monotonically
decreases with the magnetic field $|\boldsymbol{B}|$, discontinuously
increases at $|\boldsymbol{B}|=B_{\rm c}$, and again monotonically decreases
for $|\boldsymbol{B}|>B_{\rm c}$. At the zero magnetic field
$\boldsymbol{B}=0$, the half-quantized vortices having the ferromagnetic cores
have smaller $\Delta$ than that for those having the cyclic cores.

Figure 6:  A single isolated half-quantized non-Abelian vortex in the $D_{4}$
BN phase at the magnetic field $\boldsymbol{B}=1.3B_{\rm
c}\hat{\boldsymbol{z}}$. 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$ are shown. The radius of
the figure shown here is $64p_{\rm F}/(\pi m_{\rm n}T_{\rm c})\approx
1.33\times 10^{4}$ fm.

####  III.2.6 Isolated half-quantized vortices in the $D_{4}$ BN phase

Finally, we also obtain the solution for a single half-quantized vortex at
$|\boldsymbol{B}|>B_{\rm c}$. We put the boundary condition as

$\displaystyle\begin{aligned} &\psi_{2}(\theta+\pi/2)=i\psi_{2}(\theta),\\\
&\psi_{1,0,-1}=0,\\\ &\psi_{-2}(\theta+\pi)=\psi_{-2}(\theta),\end{aligned}$
(24)

or

$\displaystyle\begin{aligned}
&{[A]}_{11}(\theta+\pi/2)=-[A]_{22}(\theta+\pi/2)\\\
&\phantom{{[A]}_{11}(\theta+\pi/2)}=\frac{e^{i\pi/4}}{\sqrt{2}}\left\\{[A]_{11}(\theta)-[A]_{12}(\theta)\right\\},\\\
&[A]_{12}(\theta+\pi/2)=[A]_{21}(\theta+\pi/2)\\\
&\phantom{{[A]}_{12}(\theta+\pi/2)}=\frac{e^{i\pi/4}}{\sqrt{2}}\left\\{[A]_{11}(\theta)+[A]_{12}(\theta)\right\\},\\\
&[A]_{13}=[A]_{23}=[A]_{31}=[A]_{32}=[A]_{33}=0\\\ \end{aligned}$ (25)

at the boundary which induces half-quantized vortex solutions. Figure 6 shows
a single half-quantized vortex state at $|\boldsymbol{B}|=1.3B_{\rm c}$. The
system has a three-fold rotational symmetry around the vortex core which can
be seen in $|\psi_{-1}|^{2}$, $|\psi_{0}|^{2}$, and $|\Psi_{30}|^{2}$. The
existence of the three-fold symmetry in $|\Psi_{30}|^{2}$ is closely related
to the existence of the three soliton lines between two half-quantized
vortices for the vortex shown in Fig. 4. Such a three-fold symmetry at the
vortex core has been also observed in ultracold spin-2 atomic BECs [125].

One of the main differences of the isolated half-quantized vortex in Fig. 6
from the constituent one in the molecule in Fig. 4 is that the vortex core in
Fig. 6 is filled with the state having $\boldsymbol{S}^{2}/\rho^{2}>0$ and
$|\Psi_{30}|^{2}/\rho^{3}>0$, which corresponds to the mixed state [59] being
intermediate state between the ferromagnetic and cyclic states. Another
difference is that $\psi_{\pm 1}$ takes finite values near the vortex core.

##  IV Summary and Discussion

In this paper, we have presented vortex solutions, that is, singly quantized
vortices and half-quantized non-Abelian vortices, in the neutron ${}^{3}P_{2}$
superfluids in the case that the external magnetic field is parallel to the
angular momentum of the vortices. We have found that a singly quantized vortex
splits into two half-quantized non-Abelian vortices connected by soliton(s)
forming a vortex molecule, at any strength of the magnetic field. The main
results are summarized in Table 1. In the absence of the magnetic field in
which the UN state is the bulk ground state, the cores of the half-quantized
vortices are filled with the ferromagnetic states, and a single linear soliton
with the $D_{4}$ BN state connects the two half-quantized vortices, as shown
in Fig. 1. At the finite magnetic field below the critical one separating
$D_{2}$ and $D_{4}$ BN states, the bulk ground state is the $D_{2}$ BN state.
In this case, the cores of the half-quantized vortices are filled with the
cyclic state and are connected by three line solitons composed of the $D_{4}$
BN order, as shown in Fig. 2. Above the critical magnetic field for which the
bulk ground state is the $D_{4}$ BN state, a single half-quantized vortex is
topologically allowed stably, as constructed in Fig. 6. Nevertheless, two
half-quantized vortices are confined with three line solitons composed of
$D_{2}$ BN order, still forming a vortex molecule as shown in Fig. 4. At the
critical magnetic field, the vortex core becomes asymmetric as in Fig. 4, as
an intermediate state between the two kinds of molecules of two half-quantized
vortices connected by the three $D_{4}$ BN solitons and those connected by the
three $D_{2}$ BN solitons. We have also found that the energy of the vortex
molecule monotonically increases as the magnetic field increases, which is
continuous including the critical magnetic field as in Fig. 5(a). The distance
between the two half-quantized vortices decreases as the magnetic field
increases, except for a discontinuous jump with an increase at the critical
magnetic field in Fig. 5(b).

Except for the case of the zero magnetic field, vortex cores are always filled
with the cyclic state, which also appears even at the zero magnetic field as a
metastable state as in Fig. 5(a). Our results contradict with the BdG approach
[80, 81, 126] in which vortex cores are filled with the mixed state for
$|\boldsymbol{B}|<B_{\rm c}$ [126] without $\psi_{\pm 1}$ components and
completely separated as two isolated half-quantized vortices without linear
soliton [81] for $|\boldsymbol{B}|>B_{\rm c}$. A main possible reason for this
contradiction comes from the difference of the GL and BdG approaches. The
latter approach has an advantage in describing the microscopic structure such
as vortex cores and fermion degrees of freedom, and the cyclic state inside
cores that we have obtained could be an artifact of the low energy theory. We
should study this point in more detail with, for example, a GL expansion to
higher order. Another minor possible reason is the difference of treatments of
the boundary condition and positions of half-quantized vortices. In the
previous study, the boundary is fixed with the bulk integer vortex state [Eqs.
(10) and (14)]. The positions of the half-quantized vortices are also fixed
and treated as a differential parameter. On the other hand, neither boundary
state determined in the boundary condition (20) nor the positions of half-
quantized vortices are fixed and automatically determined to minimize the
whole free-energy density in our study. This subtle difference may affect the
difference of the vortex-core states.

Here we address further discussions for future studies. We here have studied
only the case for the magnetic field parallel to the angular momentum (the
direction of vortices). We will report the case for an arbitrary angle between
them elsewhere, which should be important for study of neutron stars in more
general situations.

Further studies should be done for multiple vortex states such as a vortex
lattice under rapid rotation relevant for neutron star interiors. In
particular, it is important to study whether, in a vortex lattice, constituent
half-quantized vortices are still tightly bound as a singly quantized vortex
as found in this paper or they are separated by distances of the same order as
two-component BECs [108, 109, 114] and whether a lattice is triangular or
square.

Another important subject is a collision dynamics of vortices in three spatial
dimensions. It is important whether two vortices reconnect in collision or a
formation of a rung between them occurs as the case of non-Abelian vortices in
the cyclic phase of spin-2 spinor BECs [127]. The presence or absence of such
a vortex reconnection is crucial for states of the quantum turbulence. With
this regards, vortex reconnection was reported in the nematic phase of a
spin-2 BEC [86], a superfluid similar to ${}^{3}P_{2}$ superfluids. Collision
of two vortex molecules may be accompanied by swapping partners as found in
vortex molecules in two-component BECs [120].

The coexistence of ${}^{1}S_{0}$ and ${}^{3}P_{2}$ superfluids drastically
changes the phase diagram [128], and vortex states in this case will be also
one direction to be explored. In this case, vortices having winding only in
either of ${}^{1}S_{0}$ and ${}^{3}P_{2}$ condensates would be further
fractionalized. In addition, vortex states in the ferromagnetic phase
appearing without quasi-classical approximation in the region close to the
critical temperature [60] are also interesting to be investigated.

Recently, it has been proposed that in the deep inside of neutron star cores,
the quark-hadron continuity for two-flavor quarks holds; the ${}^{3}P_{2}$
superfluid (nuclear matter) is continuously connected through crossover to a
two-flavor quark matter called the 2SC+$dd$ phase [129]; in addition to the
conventional 2SC phase, a $P$-wave condensation $\left<d\nabla d\right>$ is
suggested. Vortex structures in the 2SC+$dd$ phase were studied in Refs. [130,
131, 132] in which an exotic vortex called a non-Abelian Alice string was
found. In particular, in Ref. [131], how vortices in nuclear and quark matter
should be connected along the crossover. It will be interesting whether core
structures found in this paper are preserved or deformed along this
connection.

Finally, a novel type of the Berezinskii-Kosterlitz-Thouless (BKT) transition
of vortex molecules in two-component systems was reported in Ref. [121], and
thus the BKT transition should be investigated in neutron ${}^{3}P_{2}$
superfluids.

## Acknowledgment

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

## References

  * 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] .
  * 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] .
  * Migdal [1959] A. B. Migdal, Zh. Eksp. Teor. Fiz. 37, 249 (1959), [Sov. Phys. JETP 10, 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, Nat. 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 [2010] C. O. Heinke and W. C. G. Ho, Astrophys. J. 719, L167 (2010), 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, Mon. Not. Roy. Astron. Soc. Lett. 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] .
  * 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, Prog. Theor. Phys. 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).
  * 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] .
  * 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).
  * 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] .
  * Uchino  _et al._ [2010a] S. Uchino, M. Kobayashi, M. Nitta, and M. Ueda, Phys. Rev. Lett. 105, 230406 (2010a), arXiv:1010.2864 [cond-mat.quant-gas] .
  * Mizushima  _et al._ [2020] T. Mizushima, S. Yasui, and M. Nitta, Phys. Rev. Research 2, 013194 (2020).
  * 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).
  * Yasui and Nitta [2020] S. Yasui and M. Nitta, Phys. Rev. C 101, 015207 (2020).
  * Yasui  _et al._ [2020a] S. Yasui, C. Chatterjee, and M. Nitta, Phys. Rev. C 101, 025204 (2020a), arXiv:1905.13666 [nucl-th] .
  * 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).
  * Chatterjee  _et al._ [2017] C. Chatterjee, M. Haberichter, and M. Nitta, Phys. Rev. C96, 055807 (2017), arXiv:1612.05588 [nucl-th] .
  * Leinson [2020] L. B. Leinson, Mon. Not. Roy. Astron. Soc. 498, 304 (2020).
  * Masaki  _et al._ [2020] Y. Masaki, T. Mizushima, and M. Nitta, Phys. Rev. Research 2, 013193 (2020).
  * Masaki  _et al._ [2021] Y. Masaki, T. Mizushima, and M. Nitta, (2021), arXiv:2107.02448 [cond-mat.supr-con] .
  * Marmorini  _et al._ [2020] G. Marmorini, S. Yasui, and M. Nitta, (2020), arXiv:2010.09032 [astro-ph.HE] .
  * 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._ [2010b] S. Uchino, M. Kobayashi, and M. Ueda, Phys. Rev. A 81, 063632 (2010b).
  * 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] .
  * 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).
  * 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] .
  * 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. Carlström, 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, J. High Energy Phys. 04, 109, arXiv:1912.02685 [hep-th] .
  * Son and Stephanov [2002] D. T. Son and M. A. Stephanov, Phys. Rev. A 65, 063621 (2002), arXiv:cond-mat/0103451 .
  * 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).
  * 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] .
  * Cipriani and Nitta [2013] M. Cipriani and M. Nitta, Phys. Rev. Lett. 111, 170401 (2013), arXiv:1303.2592 [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._ [2014] M. Nitta, M. Eto, and M. Cipriani, J. Low Temp. Phys. 175, 177 (2014), 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] .
  * Eto  _et al._ [2020] M. Eto, K. Ikeno, and M. Nitta, Phys. Rev. Research 2, 033373 (2020), arXiv:1912.09014 [cond-mat.quant-gas] .
  * Kobayashi  _et al._ [2019] M. Kobayashi, M. Eto, and M. Nitta, Phys. Rev. Lett. 123, 075303 (2019), arXiv:1802.08763 [cond-mat.stat-mech] .
  * 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, J. High Energy Phys. 02, 099, arXiv:2111.13345 [hep-ph] .
  * Kobayashi  _et al._ [2011] M. Kobayashi, Y. Kawaguchi, and M. Ueda, (2011), arXiv:0907.3716 [cond-mat] .
  * [126] Y. Masaki, T. Mizushima, and M. Nitta, In preparation.
  * 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] .
  * 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, J. High Energy Phys. 09, 192, arXiv:2103.15185 [hep-ph] .


#  Angular upsampling in diffusion MRI using contextual _HemiHex_ sub-sampling
in _q_ -space

Abrar Faiyaz  Department of ECE
University of Rochester
Rochester, US
<EMAIL_ADDRESS>Md Nasir Uddin  Department of Neurology
University of Rochester
Rochester, US
<EMAIL_ADDRESS>Giovanni Schifitto  Department of Neurology
University of Rochester
Rochester, US
<EMAIL_ADDRESS>
###### Abstract

Artificial Intelligence (Deep Learning(DL)/ Machine Learning(ML)) techniques
are widely being used to address and overcome all kinds of ill-posed problems
in medical imaging which was or in fact is seemingly impossible. Reducing
gradient directions but harnessing high angular resolution(HAR) diffusion data
in MR that retains clinical features is an important and challenging problem
in the field. While the DL/ML approaches are promising, it is important to
incorporate relevant context for the data to ensure that maximum prior
information is provided for the AI model to infer the posterior. In this
paper, we introduce _HemiHex_(HH) subsampling to suggestively address training
data sampling on q-space geometry, followed by a nearest neighbor regression
training on the HH-samples to finally upsample the dMRI data. Earlier studies
has tried to use regression for up-sampling dMRI data but yields performance
issues as it fails to provide structured geometrical measures for inference.

Our proposed approach is a geometrically optimized regression technique which
infers the unknown q-space thus addressing the limitations in the earlier
studies.

###### Index Terms:

HemiHex, angular, upsampling, diffusion, MRI, deep learning, machine learning,
Subsampling, q-space, delaunay, triangulation, low angular resolution, lar,
har

##  I Introduction

Harnessing high angular resolution(HAR) with reduced number of gradients in
diffusion data in MR that retains clinical features is an important and
challenging problem in the field[1]. Interpolation techniques to estimate DWI
signals for Fiber oreientation density(fod) estimation and tractography with
dMRI has shown promising result in simulation and _in-vivo_[2, 3, 4]. We
hypothesized that if the upsampling procedure can be guided through a Delanuy
Triangulation sampling scheme, the nearest neighbor geometry could potentially
be generalized in terms of diffusion signals.

The dataset provided in the MICCAI Challenge held in the CDMRI workshop 2022
(QuaD22) has geometrically optimized sampling scheme for both upsampled and
low angular resolution q-space data and this enables us to generalize
subsamples of the q-space and finally estimate unknowns in the geometry.

##  II Method

###  II-A Optimized sampling in q-space

We have utilized Iterative Maximum Overlap Construction (IMOC) and 1 Opt
greedy method for optimizing gradient directions.[5, 6] The gradient nodes in
LAR were considered in such a way that they have the best angular coverage.
Then the HAR gradients can be optimized considering LAR gradients to be a
subset of the HAR protocol. The protocol optimization was done using dMRItool
in MATLAB.

###  II-B HemiHex Subsampling in q-space

HemiHex (HH) Subsampling can be regarded as subsampling the q-space centering
on an unknown node in such a way that the known and unknown q-space nodes fall
on approximated hexagon nodes, alternating knowns and unknowns. Hemi-Hex
interpolation refers to the regression of unknown center node on a hemi-hex
sample. The subsampling scheme results in more training data points per
subject and thus requires less subjects to be trained contextually. (Figure-1)

Figure 1: Hemi-hex subsampling on q-space, the center node is an unknown
surrounded with 3 known and unknowns.

###  II-C HemiHex Subsampling Prerequisites

  * •

Both the LAR and HAR images must have optimized sampled schemes in q space for
diffusion protocol. Optimization makes sure that the gradients are equidistant
and the hexagonal formation is possible between the known and unknown samples.

  * •

The LAR gradients used must be common to HAR gradients and should be a subset
of HAR grads. Number of gradients in LAR must be 3 times the number of
gradients in the HAR.

Figure 2: Schematics of angular upsampling in diffusion MR using contextual
HemiHex sub-sampling in q-space

###  II-D Mahcine Learning Model

A fully connected network (FCN) is used with 27x3 input layer and one node for
diffusion weighted signal output for each gradient. The schematic diagram of
the approach is shown in Figure-2. HemiHex subsampling generates training and
validation data required for training the FCN. Spatial neighbors are
incorporated in the input domain of the learner.

The network is trained in three different phases with stochastic gradient
descent with momentum, adam and rmsprop algorithms respectively. Successive
iteration and details are reported on Figure-4,5,6.

ROIs were selected from GM, WM and merging areas from three subjects to obtain
the training data. Validation data was similarly extracted from five other
subjects. In total, eight subjects are used for training and validating the
network for DWI signal prediction. Since the training takes place in
regressing the fourth dimension (the dwi signals), sufficient training data is
sampled from eight healthy subjects.

The input samples are 27 dimensional as each centering voxel posses 26
neighbours. Based on Delaneuy triangulation geometry, for each unknown signal
we always have 3 knowns on the sphere. The architecture is expected to be
rotation invariant because the training being generalized on the sphere due to
the inherent sampling geometry of the provided data. For each voxel, diffusion
gradient signals for 40 unknown gradients are estimated.

###  II-E Estimation Philosophy

The network estimates DW signals for the unknown directions trained on the
variable yet close geometric patterns. The final metric of evaluation are the
DTI metrics. Once the data is upsampled, DTI is applied on the data and the
final result is generated. The results are not shown in this article.

###  II-F Dataset

Two datasets in the study were obtained through MICCAI Challenge 2022 - QuaD22
Website

  * •

Training dataset: Diffusion-weighted MRI dataset acquired with 61 gradient
directions at b=1000 s/mm2 coming from 60 healthy controls. The sampling
scheme allows the 61 gradient directions to be easily subsampled to 21
gradients. (Resolution- 1.875x1.875x2 mm3)

  * •

Migraine dataset: A set of 50 Chronic migraine (CM) and 50 episodic migraine
(EM) patients, all acquired in a subsampled scenario with 21 gradient
directions and b=1000 s/mm2.(Resolution- 1.875x1.875x2 mm3)

_Note: Although 60 healthy controls were provided from the challenge website,
only 3 was used for training purposes. Others were used for validation and
testing._

##  III Optimization

For training the nearest neighbor regression network, we have used three
different optimization algorithms. Stochastic gradient descent with momentum
(sgdm), adam and rmsprop algorithms successively were applied on independent
training data points to minimize the loss function. Three phase training had
the same objective function for minimization,

$MSE=\sum_{i=1}^{D}(x_{i}-y_{i})^{2}$

where, $D$ is the number of total gradient directions for a voxel in x-space.
$x_{i}$ denotes the original gradient signal and $y_{i}$ is the predicted
gradient signal inferred from the subsampled input.

Figure 3: Three step optimization reduces the error in logscale.

##  IV Discussion

Angular Upsampling is an ill-posed inverse problem which can be made stable
with the help of tessellated geometries designed in this work through
sampling. Tessellated geometries ensure rotation invariant learning on the
sphere (q-space). And the neighboring voxels provide local spatial prior to
learn the tissue property for a specific gradient. This makes sure high
resolution counterpart can retain clinically relevant information through
upsampling.

The implemented approach requires significantly less computational resources
than general CNN based architectures. (We have used core i7 CPU with 16GB
memory used for training). Which took less than 10 minutes for training.

The approach is highly efficient with training resources which yields
important advantage for clinical applications and research. As less as 2/3
subjects is enough for training and acquire the contextual information for the
gradients.

##  V Conclusion

We proposed a fully connected regression network integrating the philosophy of
geometric distribution and interpolation of DWI data. The low resource need
and ease of computation makes the approach more suitable to be used in dMRI
clinical studies.

Figure 4: Training Phase-1 with stochastic gradient descent with momentum
details highlights initial loss and RMSE for training and validation
iterations. Figure 5: Training Phase-2 with ADAM highlights secondary loss and
RMSE for training and validation iterations. Figure 6: Training Phase-3 with
rmsprop highlights final loss and RMSE for training and validation iterations.

## Acknowledgements & Note

The work is a participant of MICCAI Challenge 2022, organized by CDMRI
workshop, QuaD22. The provided dataset result or the competition result is not
released or shown in this work in agreement with the organizers. Only the
developed method is presented in the article.

## References

  * [1] G. Chen, B. Dong, Y. Zhang, W. Lin, D. Shen, and P.-T. Yap, “Angular Upsampling in Infant Diffusion MRI Using Neighborhood Matching in x-q Space,” _Frontiers in Neuroinformatics_ , vol. 12, 2018, publisher: Frontiers Media SA. [Online]. Available: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6137306/
  * [2] I. B. Alaya, M. Jribi, F. Ghorbel, M. Mars, and T. Kraiem, “Quantitative evaluation of fiber tractography with a Delaunay triangulation-based interpolation approach,” _Medical & Biological Engineering & Computing_, vol. 57, no. 4, pp. 925–938, Apr. 2019.
  * [3] I. Ben Alaya, M. Jribi, F. Ghorbel, D. Sappey-Marinier, and T. Kraiem, “Fast and Accurate Estimation of the HARDI Signal in Diffusion MRI Using a Nearest-Neighbor Interpolation Approach,” _IRBM_ , vol. 38, no. 3, pp. 156–166, Jun. 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1959031817300404
  * [4] I. Ben Alaya, M. Jribi, F. Ghorbel, and T. Kraiem, “Comparison analysis of local angular interpolation methods in diffusion MRI,” _Computer Methods in Biomechanics and Biomedical Engineering: Imaging & Visualization_, vol. 0, no. 0, pp. 1–10, Jan. 2022, publisher: Taylor & Francis _eprint: https://doi.org/10.1080/21681163.2021.2024089. [Online]. Available: https://doi.org/10.1080/21681163.2021.2024089
  * [5] J. Cheng, D. Shen, and P.-T. Yap, “Designing Single- and Multiple-Shell Sampling Schemes for Diffusion MRI Using Spherical Code,” _Medical image computing and computer-assisted intervention : MICCAI … International Conference on Medical Image Computing and Computer-Assisted Intervention_ , vol. 17, no. Pt 3, pp. 281–288, 2014. [Online]. Available: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8167438/
  * [6] J. Cheng, D. Shen, P.-T. Yap, and P. J. Basser, “Single- and Multiple-Shell Uniform Sampling Schemes for Diffusion MRI Using Spherical Codes,” _IEEE transactions on medical imaging_ , vol. 37, no. 1, pp. 185–199, Jan. 2018.