# Phase diagram of the Hubbard model on a honeycomb lattice:
A cluster slave-spin study

Ming-Huan Zeng  Department of Physics, Beijing Normal University, Beijing
100875, China  Tianxing Ma<EMAIL_ADDRESS>Department of Physics, Beijing
Normal University, Beijing 100875, China  Y.-J. Wang<EMAIL_ADDRESS>Department of Physics, Beijing Normal University, Beijing 100875, China

###### Abstract

The cluster slave-spin method is implemented to research the ground state
properties of the honeycomb lattice Hubbard model with doping $\delta$ and
coupling $U$ being its parameters. At half-filling, a single direct and
continuous phase transition between the semi-metal and antiferromagnetic (AFM)
insulator is found at $U_{\text{AFM}}=2.43t$ that is in the Gross-Neveu-Yukawa
universality class, precluding the existence of the intermediate spin liquid
state. In the vicinity of this transition, a relation between the staggered
magnetization $M$ and the AFM energy gap $\Delta_{\text{AFM}}$ is established
as $M\propto\Delta_{\text{AFM}}$, compared to
$M\propto\Delta_{\text{AFM}}(\ln{\Delta_{\text{AFM}}})^{2}$ in the square
lattice case. For the underlying paramagnetic (PM) state, a first-order semi-
metal to insulator Mott transition is corroborated at $U_{\text{Mott}}=8.36t$,
characterized by discontinuities and hystereses, which is responsible for a
broad crossover around $U_{c}=5.4t$ between the weak- and strong-coupling
regimes in the AFM state. In the doped system, unlike the square lattice
Hubbard model, this crossover coupling strength $U_{c}$ increases with doping
$\delta$. The compressibility $\kappa$ near the van Hove singularity at
$\delta=1/4$ is suppressed substantially by the interaction before the semi-
metal to AFM transition occurs, whereas $\kappa$ near the Dirac points is very
close to the noninteracting one, indicating that the Dirac cone structure of
the energy dispersion is rather robust. Finally, an overall phase diagram in
the $U$-$\delta$ plane is presented, which consists of four regimes: the AFM
insulator at $\delta=0$ for $U>U_{\text{AFM}}$, the AFM metal with
compressibility $\kappa>0$ or $\kappa<0$, and the PM semi-metal, and the AFM
metal with $\kappa<0$ only exists in an extremely small area near the phase
boundary between the AFM and PM state and at intermediate coupling strengths.

###### pacs:

71.27.+a, 71.10.-w, 75.10.Jm

##  I INTRODUCTION

The Hubbard Hamitonian Hubbard (1963); *Hubbard1964 has been acting over
decades as a prototypical model for the description of interacting electrons.
In spite of its seeming simplicity, this model captures a rich phenomenology
of strongly-correlated electrons such as metallic-insulating, nonmagnetic-
antiferromagnetic and normal-superconducting phase transitions and can not be
solved exactly in more than one dimensions, which necessitates some
nonperturbative approaches to deal with the strong-coupling aspect of the
model White  _et al._ (1989); Qin  _et al._ (2021).

In this paper, we will focus on the one-band Hubbard model defined on a
honeycomb network, which is bipartite and admits the antiferromagnetism (AFM)
in the strong coupling limit. This model has a linear free electron energy
dispersion with nodal gapless points at the corners of the Brillouin zone,
leading to the so-called Dirac semi-metal. Due to the gapless Dirac points,
there is a nontrivial semi-metal to antiferromagnetic insulator (AFMI)
transition at a finite coupling strength at half-filling, which makes this
model an ideal playground to research the interaction-driven semi-metal to
AFMI transition. Up till now, many numerical and analytical methods have been
applied to the half-filled system to study this transition and its critical
behaviors. Large scale quantum Monte Carlo (QMC) simulations of $648$ sites
predict a spin liquid state in a range of interaction $3.5t<U<4.3t$, beyond
which the AFM sets in Meng  _et al._ (2010), and this argument was supported
by some numerical works Hohenadler  _et al._ (2011, 2012a, 2012b); Zheng  _et
al._ (2011). Nevertheless, this picture was disputed by many other numerical
studies Otsuka  _et al._ (2016); Raczkowski  _et al._ (2020); Assaad and
Herbut (2013); Paiva  _et al._ (2005); Yamada (2016); Sorella and Tosatti
(1992); Sorella  _et al._ (2012); Ostmeyer  _et al._ (2020, 2021), especially
those using the same method containing up to $2592$ sites Sorella  _et al._
(2012); Otsuka  _et al._ (2016) and $20808$ sites Ostmeyer  _et al._ (2020,
2021). By means of cluster dynamical mean-field theory, variational cluster
approximation, and cluster dynamical impurity approximation, Hassan _et al._
Hassan and Sénéchal (2013) showed that the results are dependent on the shape
and size of the clusters, and they claimed that only the system with two bath
orbits per cluster boundary site is able to describe the correct behavior and
found that the Mott transition for the spin liquid state is actually preempted
by the AFM long-range order. Though the early variational cluster calculations
Seki and Ohta (2012) argued that the single-particle gap opens at an
infinitesimal value of $U$, recent dynamical cluster approximation study found
that this spurious excitation gap is due to the violation of the translation
symmetry of the system and the cluster with one bath orbital per cluster site
is sufficient for the description of the short-range correlations within the
honeycomb unit cell Liebsch and Wu (2013a). Recent study by density matrix
embedding theory revealed a paramagnetic insulating state with possible
honeycomb cluster state at intermediate coupling strength whose stability is
highly cluster and lattice size dependent, and this state is nonexistent in
the thermodynamic limit, signaling no intermediate state in the half-filled
Hubbard model on a honeycomb lattice Chen  _et al._ (2014). In addition, a
two-particle self-consistent study presented a semi-metal to AFMI transition
and proved that the transition from a semi-metal to spin liquid phase is
forestalled by this transition Arya  _et al._ (2015). The functional
renormalization group theory predicts a critical interaction strength $U=3.8t$
that is consistent with the results from the methods mentioned above,
supporting that there is no spin liquid state at intermediate coupling
strengths Honerkamp (2008); Raghu  _et al._ (2008).

Based on the charge-spin separation theory Wang  _et al._ (1993); Feng  _et
al._ (1993); *Feng_1994; *Feng_2003; *Feng_2004; *Feng_2015; Hassan and de’
Medici (2010), the $U(1)$ slave-spin method has been proposed to cope with the
Mott transition in multi-orbital systems Yu and Si (2012). Then, a cluster
slave-spin approach was developed to address strongly correlated systems Lee
and Lee (2017) and has been employed to solve the square lattice Hubbard model
to obtain an overall ground state phase diagram in the parameter space of
doping $\delta$ and interaction $U$ Zeng  _et al._ (2021). In the present
work, we apply this method with the Lanczos exact-diagonalization as the
cluster solver to the honeycomb lattice Hubbard model to study its ground
state properties, including the quantum critical behavior in the vicinity of
the interaction-driven semi-metal to AFMI transition at half-filling and an
overall phase diagram in the whole $U$-$\delta$ plane. Our motivation is two-
fold: (i) Because the results of this model is shown to be highly dependent on
the size of the lattice adopted for QMC simulations, as well as the size and
shape of the clusters used within various cluster approximations, more results
from different approaches ought to be included and compared with each other.
(ii) Away from half-filling, much attention was paid to the 1/4-doping, where
the free density of states shows a van Hove singularity of logarithmic type,
favoring an instability towards superconductivity in the weak interaction
regime Gu  _et al._ (2013); Wang  _et al._ (2012); Nandkishore  _et al._
(2012), whereas an overall $U$-$\delta$ phase diagram pertaining to the
magnetism is still absent.

In the honeycomb lattice Hubbard model, we find that the first-order Mott
transition occurs at $U_{\text{Mott}}=8.36t$ in the half-filled PM state,
characterized by discontinuities and hystereses in all quantities, and
transforms into a broad crossover in the AFM state because of long range AFM
correlations. Besides, the phase separation, manifested by a negative
compressibility, has been observed in a region near the phase boundary
$\delta_{M}(U)$ between the AFM and PM state and at intermediate couplings,
whose area is much smaller compared to the square lattice Hubbard model Zeng
_et al._ (2021). Finally, a phase diagram in the $U$-$\delta$ plane is
presented, consisting of four regimes: AFMI, AFM metal with positive and
negative compressibility, and the PM semi-metal.

The rest of this paper is organized as follows. In Sec. II, we reintroduce the
cluster slave-spin mean-field theory Lee and Lee (2017); Zeng  _et al._ (2021)
and implement it in the honeycomb lattice Hubbard model by making use of two-
and six-site cluster approximations. In Sec. III.1, for the half-filled
system, an analytical relation between the staggered magnetization $M$ and the
AFM energy gap $\Delta_{\text{AFM}}$ in the vicinity of the semi-metal to AFMI
transition is established, and the first-order Mott transition at
$U_{\text{Mott}}=8.36t$ is observed in the PM state. In Sec. III.2, the
results of finite doping cases obtained by two- and six-site clusters are
discussed thoroughly, and we find that the two-site cluster is inadequate to
capture the AFM transition appropriately because it violates the symmetry of
the honeycomb lattice. In Sec. IV, the properties of $M$,
$\Delta_{\text{AFM}}$, and the compressibility $\kappa$ are combined to show a
phase diagram of the model in the $U$-$\delta$ plane.

##  II Formalism

The standard one-band fermionic Hubbard model Hubbard (1963, 1964) reads

$H=-t\sum_{\langle
i,j\rangle\sigma}(c_{i\sigma}^{\dagger}c_{j\sigma}+\mbox{\sc
h.c.})+U\sum_{i}n_{i\uparrow}n_{i\downarrow}-\mu\sum_{i\sigma}n_{i\sigma}\;,$
(1)

where $t$, $U$, $\mu$ are the nearest hopping constant, the on-site Coulomb
repulsion energy and the chemical potential, respectively. The sum $\langle
i,j\rangle$ runs over all pairs of nearest-neighbor sites on a honeycomb
lattice, and $c_{i\sigma}^{\dagger}$ is the creation operator of the electron
at site $i$ with spin $\sigma=\uparrow,\,\downarrow$, and the number operator
$n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$. Hereafter, we use $t=1$ as the
unit of energy.

In the U(1) slave-spin method Yu and Si (2012), an electron operator is
factorized into a slave-spin operator and a fermionic spinon operater

$c^{\dagger}_{\alpha}\equiv S^{\dagger}_{\alpha}f^{\dagger}_{\alpha}\;.$ (2)

With the constraint
$a^{\dagger}_{\alpha}a_{\alpha}+b^{\dagger}_{\alpha}b_{\alpha}=1$, the slave-
spin operator is rewritten in the Schwinger boson representation

$S^{\dagger}_{\alpha}=a^{\dagger}_{\alpha}b_{\alpha}\;,\;\;\;\;S^{z}_{\alpha}=\frac{1}{2}\left(a^{\dagger}_{\alpha}a_{\alpha}-b^{\dagger}_{\alpha}b_{\alpha}\right)\;.$
(3)

On account of the charge-spin separation (2), the original Hillbert space with
basis $\\{|0\rangle,|1\rangle\\}$ is enlarged to
$\\{|n^{f}_{\alpha},S^{z}_{\alpha}\rangle\\}=\\{|0,-\tfrac{1}{2}\rangle,|1,\tfrac{1}{2}\rangle,|0,\tfrac{1}{2}\rangle,|1,-\tfrac{1}{2}\rangle\\}$.
Thus, an extra constraint needs to be imposed to restrict the Hillbert space
to the physical one:
$\\{|n^{f}_{\alpha},S^{z}_{\alpha}\rangle\\}=\\{|0,-\tfrac{1}{2}\rangle,|1,\tfrac{1}{2}\rangle\\}$,

$S^{z}_{\alpha}=f^{\dagger}_{\alpha}f_{\alpha}-\tfrac{1}{2}\;.$ (4)

To ensure the correct non-interacting behaviors, the slave-boson operators
need to be dressed as follows Kotliar and Ruckenstein (1986)

$\displaystyle\tilde{S}^{\dagger}_{\alpha}$ $\displaystyle=$ $\displaystyle
P^{\dagger}_{\alpha}a^{\dagger}_{\alpha}b_{\alpha}P_{\alpha}\;,$ (5a)
$\displaystyle P^{\pm}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{1/2\pm S^{z}_{\alpha}}}\;.$ (5b)

Following the recipe of Lee and Lee Lee and Lee (2017), with the local
constraints (4) being ensured roughly by two global Lagrange multipliers
$\lambda_{I\sigma}$ on sublattices $I=A$ and $B$, Hamiltonian (1) can be cast
into the form

$\displaystyle H^{f}_{\text{MF}}$ $\displaystyle=$ $\displaystyle-
tZ\sum_{\langle
i,j\rangle\sigma}(a_{i\sigma}^{\dagger}b_{j\sigma}+\textsc{h.c.})$ (6a)
$\displaystyle-\sum_{i\sigma}\big{[}(\mu+\lambda_{A\sigma}-\tilde{\mu}_{A\sigma})a_{i\sigma}^{\dagger}a_{i\sigma}$
$\displaystyle\hskip
23.00006pt+(\mu+\lambda_{B\sigma}-\tilde{\mu}_{B\sigma})b_{i\sigma}^{\dagger}b_{i\sigma}\big{]}\;,$
$\displaystyle H_{n_{c}\text{-site}}^{S}$ $\displaystyle=$ $\displaystyle
H^{\lambda}_{n_{c}\text{-site}}+H^{U}_{n_{c}\text{-site}}+H^{K}_{n_{c}\text{-site}}\;,$
(6b)

where

$\displaystyle H^{\lambda}_{n_{c}\text{-site}}$ $\displaystyle=$
$\displaystyle\sum_{i_{c}=1,\sigma}^{n_{c}}\lambda_{I\sigma}S_{i_{c}\sigma}^{z}\;,$
(7a) $\displaystyle H^{U}_{n_{c}\text{-site}}$ $\displaystyle=$
$\displaystyle\sum_{i_{c}=1}^{n_{c}}U(S_{i_{c}\sigma}^{z}+\tfrac{1}{2})(S_{i_{c}\bar{\sigma}}^{z}+\tfrac{1}{2})\;,$
(7b) $\displaystyle H^{K}_{\text{2-site}}$ $\displaystyle=$
$\displaystyle\sum_{\sigma}\Big{\\{}\epsilon_{\sigma}^{\delta_{1}}(\tilde{z}_{A\sigma}^{\dagger}\tilde{z}_{B\sigma}+\tilde{z}_{B\sigma}^{\dagger}\tilde{z}_{A\sigma})+(\epsilon_{\sigma}^{\delta_{2}}+\epsilon_{\sigma}^{\delta_{3}})$
(7c) $\displaystyle\hskip
20.00003pt\times\big{[}\tilde{z}_{A\sigma}^{\dagger}\langle\tilde{z}_{B\sigma}\rangle+\tilde{z}_{B\sigma}^{\dagger}\langle\tilde{z}_{A\sigma}\rangle+\textsc{h.c.}\big{]}\Big{\\}}\;,$
$\displaystyle H^{K}_{\text{6-site}}$ $\displaystyle=$
$\displaystyle\sum_{\sigma}\Big{\\{}\epsilon_{\sigma}^{\delta_{1}}(\tilde{z}_{1\sigma}^{\dagger}\tilde{z}_{2\sigma}+\tilde{z}_{4\sigma}^{\dagger}\tilde{z}_{5\sigma})$
(7d)
$\displaystyle+\epsilon_{\sigma}^{\delta_{2}}(\tilde{z}_{1\sigma}^{\dagger}\tilde{z}_{6\sigma}+\tilde{z}_{3\sigma}^{\dagger}\tilde{z}_{4\sigma})+\epsilon_{\sigma}^{\delta_{3}}(\tilde{z}_{2\sigma}^{\dagger}\tilde{z}_{3\sigma}+\tilde{z}_{5\sigma}^{\dagger}\tilde{z}_{6\sigma})$
$\displaystyle+\epsilon_{\sigma}^{\delta_{1}}(\tilde{z}_{3\sigma}^{\dagger}\langle\tilde{z}_{6\sigma}\rangle+\langle\tilde{z}_{3\sigma}^{\dagger}\rangle\tilde{z}_{6\sigma})$
$\displaystyle+\epsilon_{\sigma}^{\delta_{2}}(\tilde{z}_{2\sigma}^{\dagger}\langle\tilde{z}_{5\sigma}\rangle+\langle\tilde{z}_{2\sigma}^{\dagger}\rangle\tilde{z}_{5\sigma})$
$\displaystyle+\epsilon_{\sigma}^{\delta_{3}}(\tilde{z}_{1\sigma}^{\dagger}\langle\tilde{z}_{4\sigma}\rangle+\langle\tilde{z}_{1\sigma}^{\dagger}\rangle\tilde{z}_{4\sigma})+\textsc{h.c.}\Big{\\}}\;.$

The slave-spin Hamiltonian, describing the electron’s charge degrees of
freedom, will be treated with $n_{c}$-site clusters ($n_{c}=$ 2, 6) marked by
the red color geometry in Fig. 1. The cluster will be solved by using the
Lanczos exact-diagonalization method. The parameters $Z$,
$\tilde{\mu}_{I\sigma}$, and $\epsilon_{\sigma}^{\delta}$ in Eqs. (6) and (7)
are calculated as follows:

$\displaystyle Z$ $\displaystyle=$
$\displaystyle\langle\tilde{z}_{A\sigma}^{\dagger}\rangle\langle\tilde{z}_{B\sigma}\rangle,\;\hskip
20.00003pt\epsilon_{\sigma}^{\delta_{1/2/3}}=-t\langle
a_{i\sigma}^{\dagger}b_{i+\hat{\delta}_{1/2/3}\sigma}\rangle,$ (8)
$\displaystyle\tilde{\mu}_{I\sigma}=\frac{2Z\langle
S_{I\sigma}^{z}\rangle(\epsilon_{\sigma}^{\delta_{1}}+\epsilon_{\sigma}^{\delta_{2}}+\epsilon_{\sigma}^{\delta_{3}})}{(\tfrac{1}{2})^{2}-\langle
S_{I\sigma}^{z}\rangle^{2}}\;.$

Figure 1: Schematic illustration of the (a) two- and (b) six-site cluster
configurations. The lattice constant $a$ is set to unity, and the positions of
three nearest neighbors of site A are $\delta_{1}=(1,0)$,
$\delta_{2}=(-1/2,\sqrt{3}/2)$, $\delta_{3}=-(1/2,\sqrt{3}/2)$.

Moreover, the fermionic spinon Hamiltonian can be Fourier transformed into
momentum space:

$\displaystyle H^{f}_{\text{MF}}$ $\displaystyle=$
$\displaystyle\sum_{\bm{k},\sigma}\big{(}\varepsilon_{A\sigma}a_{\bm{k}\sigma}^{\dagger}a_{\bm{k}\sigma}+\varepsilon_{B\sigma}b_{\bm{k}\sigma}^{\dagger}b_{\bm{k}\sigma}$
(9) $\displaystyle\hskip
20.00003pt+\Gamma_{\bm{k}}a_{\bm{k}\sigma}^{\dagger}b_{\bm{k}\sigma}+\Gamma_{\bm{k}}^{*}b_{\bm{k}\sigma}^{\dagger}a_{\bm{k}\sigma}\big{)}$

with

$\displaystyle\varepsilon_{A/B\sigma}=\tilde{\mu}_{A/B\sigma}-\mu-\lambda_{A/B\sigma}\;,$
$\displaystyle\Gamma_{\bm{k}}=-tZ\gamma_{\bm{k}}\;,\;\;\;\;\gamma_{\bm{k}}=\sum_{\bm{\delta}}e^{i\bm{k}\cdot\bm{\delta}}\;.$
(10)

Diagnalization of the spinon Hamiltonian (9) gives rise to the eigenenergies
as

$\displaystyle E_{\bm{k}}^{\pm}$ $\displaystyle=$
$\displaystyle\pm\sqrt{|\Gamma_{\bm{k}}|^{2}+\Delta_{\sigma}^{2}}-\mu_{\text{eff}}\;,$
(11a) $\displaystyle\mu_{\text{eff}}$ $\displaystyle=$
$\displaystyle\mu-\frac{1}{2}(\tilde{\mu}_{A\sigma}-\lambda_{A\sigma}+\tilde{\mu}_{B\sigma}-\lambda_{B\sigma})\;,$
(11b) $\displaystyle\Delta_{\sigma}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\tilde{\mu}_{A\sigma}-\lambda_{A\sigma}-\tilde{\mu}_{B\sigma}+\lambda_{B\sigma})\;.$
(11c)

Here, AFM energy gap $\Delta_{\text{AFM}}=|\Delta_{\sigma}|$ is identical in
form to that in the square lattice caseLee and Lee (2017); Zeng  _et al._
(2021).

In most occasions, it proves effective to adopt the density of states (DOS) of
the non-interacting electrons to calculate the physical quantities in the
thermodynamic limit. On the hexagon lattice, it is defined as

$\displaystyle D(\gamma)$ $\displaystyle=$ $\displaystyle\frac{1}{{\cal
N}_{\text{triangle}}}\sum_{\bm{k}}\delta(\gamma-|\gamma_{\bm{k}}|)$ (14)
$\displaystyle=$
$\displaystyle\Big{\\{}\begin{array}[]{ll}N(\gamma)\;,&\;\;\;(0\leq\gamma<1)\\\
\tilde{N}(\gamma)\;,&\;\;\;(1<\gamma\leq 3)\end{array}$

where

$\Bigg{\\{}\begin{array}[]{l}N(\gamma)=\frac{4}{\pi^{2}}\frac{\gamma}{\sqrt{(3-\gamma)(1+\gamma)^{3}}}\,K\Big{(}\sqrt{\frac{16\gamma}{(3-\gamma)(1+\gamma)^{3}}}\,\Big{)}\;,\\\
\tilde{N}(\gamma)=\frac{1}{\pi^{2}}\sqrt{\gamma}\,K\Big{(}\sqrt{\frac{(3-\gamma)(1+\gamma)^{3}}{16\gamma}}\,\Big{)}\;,\end{array}$
(15)

with $K(x)$ being the first kind complete elliptical integral. In comparison
to the self-dual situation on a square lattice, we now have a duality
transformation

$\tilde{\gamma}=\frac{3-\gamma}{1+\gamma}$ (16)

to connect these two parts, under which

$\displaystyle\frac{(3-\tilde{\gamma})(1+\tilde{\gamma})^{3}}{16\tilde{\gamma}}$
$\displaystyle=$ $\displaystyle\frac{16\gamma}{(3-\gamma)(1+\gamma)^{3}}\;,$
(17a)
$\displaystyle\frac{(3-\tilde{\gamma})(1+\tilde{\gamma})}{4\tilde{\gamma}}$
$\displaystyle=$ $\displaystyle\frac{4\gamma}{(3-\gamma)(1+\gamma)}\;,$ (17b)

and

$\displaystyle\tilde{N}(\gamma)$ $\displaystyle=$
$\displaystyle\frac{(3-\tilde{\gamma})(1+\tilde{\gamma})}{4\tilde{\gamma}}N(\tilde{\gamma})$
(18) $\displaystyle=$
$\displaystyle\frac{4\gamma}{(3-\gamma)(1+\gamma)}N\Big{(}\frac{3-\gamma}{1+\gamma}\Big{)}\;.$

Then, the self-consistent quantities $\epsilon^{\delta}_{\sigma}=\epsilon$ and
$n_{(A/B)\sigma}\equiv\langle a^{\dagger}_{i\sigma}a_{i\sigma}\rangle/\langle
b^{\dagger}_{i\sigma}b_{i\sigma}\rangle$ can be calculated through

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\int_{0}^{3}d\gamma\,D(\gamma)\frac{(tZ\gamma)^{2}}{6Z\sqrt{(tZ\gamma)^{2}+\Delta^{2}}}$
(19)
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\times\Big{\\{}\theta[-E^{+}(\gamma)]-\theta[-E^{-}(\gamma)]\Big{\\}}\;,$
$\displaystyle n_{(A/B)\sigma}$ $\displaystyle=$
$\displaystyle\int_{0}^{3}d\gamma\,D(\gamma)\bigg{[}\frac{1}{2}\Big{(}1\pm\frac{\Delta_{\sigma}}{\sqrt{(tZ\gamma)^{2}+\Delta^{2}}}\Big{)}\theta(-E_{\gamma}^{+})$
(20) $\displaystyle\hskip
10.00002pt+\frac{1}{2}\Big{(}1\mp\frac{\Delta_{\sigma}}{\sqrt{(tZ\gamma)^{2}+\Delta^{2}}}\Big{)}\theta(-E_{\gamma}^{-})\bigg{]}\;,$

where $E^{\pm}(\gamma)=\pm\sqrt{(tZ\gamma)^{2}+\Delta^{2}}-\mu_{\text{eff}}$
and $\Delta_{\sigma}^{2}=\Delta^{2}$.

##  III RESULTS AND DISCUSSIONS

###  III.1 HALF-FILLED SYSTEM

In this case, the particle-hole symmetry implies $\mu^{\text{eff}}=0$ and
$E^{+}(\gamma)>0$, and Eqs. (19) and (20) are simply

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle-\frac{t\lambda}{6}\int_{0}^{3}d\gamma\,D(\gamma)\frac{\gamma^{2}}{\sqrt{(\gamma\lambda)^{2}+1}}$
(21) $\displaystyle=$
$\displaystyle-\frac{t\lambda}{6}I_{\epsilon}(\lambda)\;,$
$\displaystyle\langle S_{(A/B)\sigma}^{z}\rangle$ $\displaystyle=$
$\displaystyle(-/+)\frac{\,{\rm
sgn}({\Delta}_{\sigma})\,}{2}\int_{0}^{3}d\gamma\,D(\gamma)\frac{1}{\sqrt{(\gamma\lambda)^{2}+1}}$
(22) $\displaystyle=$ $\displaystyle(-/+)\frac{\,{\rm
sgn}({\Delta}_{\sigma})\,}{2}I_{S}(\lambda)\;,$

where $\lambda=tZ/\Delta_{\text{AFM}}$ and, by relation (4), $\langle
S_{(A/B)\sigma}^{z}\rangle=n_{(A/B)\sigma}-\frac{1}{2}$.

Figure 2: (a) The AFM energy gap $\Delta_{\text{AFM}}$ and (b) the staggered
magnetization $M$ as function of $U$, where the insets show the same data in
the vicinity of the critical coupling together with the fitting data (blue
lines). The results from two-site cluster approximation are presented as well
(black).

For the half-filled square lattice Hubbard model at $T=0$, the AFM order
emerges for any nonzero $U$ because of the perfect nesting of the free Fermi
surface, whereas the honeycomb lattice is known to have a semi-metal phase at
small $U$ due to the low coordination number which allows more fluctuations
and an antiferromagnetic phase at large $U$. It is well established that there
is a single direct and continuous phase transition from semi-metal to AFMI at
a finite critical interaction strength $U_{\text{AFM}}$ for the half-filled
honeycomb lattice Hubbard model. This $U_{\text{AFM}}$ from large-scale QMC
simulations mainly locates around $U\approx 3.8$ Sorella  _et al._ (2012);
Otsuka  _et al._ (2016); Ostmeyer  _et al._ (2021); Raczkowski  _et al._
(2020); Assaad and Herbut (2013); Ma  _et al._ (2018), whereas the results
from various cluster scenarios, such as cluster dynamical impurity
approximation, variational cluster approximation, dynamical cluster
approximation and density matrix embedding theory, are strongly cluster
dependent and the $U_{\text{AFM}}$’s are in a wide range of $1.5\lesssim
U\lesssim 4.0$ Seki and Ohta (2012); Chen  _et al._ (2014); Liebsch and Wu
(2013b); Yamada (2016); Hassan and Sénéchal (2013). As shown in Fig. 2, we
find that $U_{\text{AFM}}=2.8$ or 2.43 in our two- or six-site cluster
approximation, which is larger than that from the Hartree-Fock approximation
of 2.235 Sorella and Tosatti (1992); Raczkowski  _et al._ (2020), but smaller
than those from QMC simulations Otsuka  _et al._ (2016); Assaad and Herbut
(2013); Sorella  _et al._ (2012); Ostmeyer  _et al._ (2021); Ma  _et al._
(2018), reflecting the fact that fluctuations have been incorporated in the
six-site cluster, but not enough to give the accurate value. This shortcoming
may be remedied by enlarging the cluster size and strictly dealing with the
constraint $S^{z}_{\alpha}=f^{\dagger}_{\alpha}f_{\alpha}-\tfrac{1}{2}$
locally. However, the two-site cluster value of $U_{\text{AFM}}$ is larger
than that from the six-site, necessitating more investigations on the
dependence of $U_{\text{AFM}}$ upon the cluster size. To extract the critical
information around $U_{\text{AFM}}$, we fit our self-consistent data from the
six-site cluster using the $U$ dependent form of $M$ and $\Delta_{\text{AFM}}$
that have been verified by QMC simulations Otsuka  _et al._ (2016); Assaad and
Herbut (2013); Sorella  _et al._ (2012); Ostmeyer  _et al._ (2021) and density
matrix embedding theory Chen  _et al._ (2014)

$\Delta_{\text{AFM}}/M=\alpha_{1/2}|U-U_{\text{AFM}}|^{\beta_{1/2}}\;,$ (23)

and we obtain

$\Big{\\{}\begin{array}[]{ll}\alpha_{1}=0.74021\pm
0.00836\;,&\beta_{1}=0.77086\pm 0.01597\;;\\\ \alpha_{2}=0.19739\pm
0.000753\;,&\beta_{2}=0.755\pm 0.0049\;.\end{array}$ (24)

The critical exponent $\beta_{2}=0.755\pm 0.0049$ for $M$ is very close to
those from the large-scale QMC simulations $\beta_{2}=0.75\pm 0.06$ Otsuka
_et al._ (2016) and 0.79 Assaad and Herbut (2013), and the critical exponent
for single particle gap $\Delta_{\text{AFM}}$ is slightly larger than that of
$M$ Assaad and Herbut (2013); Ostmeyer  _et al._ (2021).

We now expand asymptotically the integrals $I_{\epsilon}(\lambda)$ and
$I_{S}(\lambda)$ defined in Eqs. (21) and (22) as $\lambda\to\infty$
Bleinstein and Handelsman (1986):

$\epsilon=-\frac{t}{6}\Big{(}1.574597-0.448221\tilde{\Delta}^{2}+\frac{4}{3\sqrt{3}\pi}\tilde{\Delta}^{3}+\cdots\Big{)}\;,$
(25) $\displaystyle\langle S_{(A/B)\sigma}^{z}\rangle$ $\displaystyle=$
$\displaystyle(-/+)\frac{\,{\rm
sgn}({\Delta}_{\sigma})\,}{2}\Big{(}0.896441\tilde{\Delta}-\frac{2}{\sqrt{3}\pi}\tilde{\Delta}^{2}$
(26)
$\displaystyle-0.014\tilde{\Delta}^{3}+\frac{4}{9\sqrt{3}\pi}\tilde{\Delta}^{4}+\cdots\Big{)}\;,$

where $\tilde{\Delta}=\lambda^{-1}=\Delta_{\text{AFM}}/tZ$. The relation
between $M$ and $\Delta_{\text{AFM}}$ in the honeycomb lattice around
$U=U_{\text{AFM}}$ is established as

$\displaystyle M$ $\displaystyle=$ $\displaystyle\tfrac{1}{2}|\langle
S_{A\sigma}^{z}\rangle-\langle S_{B\sigma}^{z}\rangle|$ $\displaystyle\approx$
$\displaystyle\frac{1}{2}\Big{(}0.896441\tilde{\Delta}-\frac{2}{\sqrt{3}\pi}\tilde{\Delta}^{2}-0.014\tilde{\Delta}^{3}+\frac{4}{9\sqrt{3}\pi}\tilde{\Delta}^{4}\Big{)}\;.$

To the leading order, $M\propto\Delta_{\text{AFM}}$, compared to
$M\propto\Delta_{\text{AFM}}(\ln{\Delta_{\text{AFM}}})^{2}$ in the square
lattice, supporting the AFM at small $U$ in the latter case is driven by the
perfect nesting of its free Fermi surface.

On the other hand, $\Delta_{\text{AFM}}$ reaches its maximum around the
crossover coupling strength $U_{c}=5.4$ that separates the weak- and strong-
coupling regimes, which is consistent with the traditional mean-field behavior
$\Delta_{\text{AFM}}\sim U$ at small $U$, and $\Delta_{\text{AFM}}\sim
4t^{2}/U$ in the large $U$ limit supported by the super-exchange mechanism. It
ought to be mentioned that $M$ drops abruptly when $U$ is larger than a
certain value where the quasi-particle weight happens to drop to zero as shown
in Fig. 3(a), which may imply that the cluster slave-spin method breaks down
for the half-filled system in the large $U$ limit.

The quasi-particle residue $Z$, the generalized Gutzwiller factor $g_{t}$
Ogawa  _et al._ (1975); Abram  _et al._ (2013); Lee and Lee (2017), the holon-
doublon correlators $C_{ij}$ Lee and Lee (2017) between the nearest neighbors
$C_{12}$, the next-nearest neighbors $C_{13}$ and the next-next-nearest
neighbors $C_{14}$, the ground state energy of the slave-spin Hamiltonian per
site $\langle H_{n_{c}\text{-site}}^{S}\rangle/n_{c}$ with $n_{c}$ being the
cluster size, and the double occupancy $\langle D\rangle$ as function of $U$
at half-filling obtained from two- and six-site cluster approximations are
presented in Fig. 3, where $C_{ij}$ is defined as

$C_{ij}=\frac{\langle N_{i}D_{j}\rangle-\langle N_{i}\rangle\langle
D_{j}\rangle}{\langle N_{i}\rangle\langle D_{j}\rangle}$ (28)

with the holon operator $N_{i}=(1-n_{i\sigma})(1-n_{i\bar{\sigma}})$ and
doublon operator $D_{j}=n_{j\sigma}n_{j\bar{\sigma}}$.

Figure 3: (a) The quasiparticle weight $Z$ and the generalized Gutzwiller
factor $g_{t}$. (b)–(d) The holon-doublon correlators between the nearest
neighbors $C_{12}$, the next-nearest neighbors $C_{13}$ and the next-next-
nearest neighbors $C_{14}$. (e) The expectation value of the cluster slave-
spin Hamiltonian $\langle H_{n_{c}\text{-site}}^{S}\rangle/n_{c}$. (f) The
double occupancy $\langle D\rangle$ in the AFM state vs. $U$ obtained by the
2/6-site clusters (blue/red). All quantities in the PM state obtained by
6-site cluster are black lines.

The results are as follows: (i) All quantities in the PM state show
discontinuities and hystereses at the critical coupling strength
$U_{\text{Mott}}=8.36$ for the semi-metal to paramagnetic insulator transition
as the characteristics of the first-order Mott transition in the PM state
Hassan and Sénéchal (2013); Tran and Kuroki (2009). (ii) In Fig. 3(a),
compared to that in the PM state, $Z$ is largely suppressed as entering the
AFM phase. (iii) In Fig. 3(b)–(d), $C_{12}$ and $C_{13}$ are positive and
increase monotonically with $U$, signaling that the holon and doublon between
the nearest and next-nearest neighbors tend to attract each other, and this
tendency is enhanced by the coupling strength. However, $C_{14}$ presents a
minimum beyond the AFM transition or as $U$ approaches $U_{\text{Mott}}$ in
the PM state. On the other hand, $C_{14}$ is positive at small and large $U$
while negative for intermediate interaction strengths, suggesting that at
half-filling, the holon and doublon between the next-next-nearest neighbors
attract each other when $U$ is small or large, while behave repulsively at
intermediate $U$. (iv) In Fig. 3(e), $\langle
H_{n_{c}\text{-site}}^{S}\rangle/n_{c}$ in the AFM state is smaller than that
from the PM state, favoring an AFM ground state. (v) In Fig. 3(f), $\langle
D\rangle$ in the PM state decreases linearly with increasing $U$ when $U\ll
U_{\text{Mott}}$ Vollhardt (1984), whereas in the AFM state, its slope changes
abruptly as AFM sets in denoting a second-order transition from a semi-metal
to an AFMI. The double occupancy in the AFM state is smaller than that in the
PM state when $U_{\text{AFM}}<U\lesssim U_{\text{Mott}}$ or the opposite is
true when $U\gtrsim U_{\text{Mott}}$.

###  III.2 SYSTEMS WITH FINITE DOPING

Figure 4: (a) The quasi-particle weight $Z$ and the generalized Gutzwiller
factor $g_{t}$, (b) the AFM energy gap $\Delta_{\text{AFM}}/t$, (c) the
staggered magnetization $M$, (d) the holon-doublon correlators between the
nearest neighbors $C_{12}$, the next-nearest neighbors $C_{13}$ and the next-
next-nearest neighbors $C_{14}$, (e) the expectation value of the cluster
slave-spin Hamiltonian $\langle H_{n_{c}\text{-site}}^{S}\rangle/n_{c}$, and
(f) the double occupancy $\langle D\rangle$ as function of $U$ at
$\delta=0.02$ in the AFM state obtained by the two-site (blue) and six-site
(red) cluster, as well as the PM state (black) by six-site cluster.

In Fig. 4, we plot $Z$, $g_{t}$, $\Delta_{\text{AFM}}$, $M$, $C_{12/3/4}$,
$\langle H_{n_{c}\text{-site}}^{S}\rangle/n_{c}$, $\langle D\rangle$ as
function of $U$ at $\delta=0.02$ obtained from two- and six-site cluster
approximations. In Fig. 4(a), the quasi-particle residue from six-site cluster
is much smaller than that from two-site because of more quantum fluctuations,
and becomes flattened when $U>U_{\text{Mott}}$. In Fig. 4(b), the critical
coupling strength for AFM transition from six-site cluster is
$U_{\text{AFM}}=3.0$, while that from two-site locates at
$U_{\text{AFM}}=25.0$, which indicates that two-site cluster is inadequate to
describe the AFM transition in the honeycomb lattice because it keeps no track
of the lattice symmetry. Moreover, in six-site cluster, the coupling strength
for the crossover from weak to strong shifts from $U_{c}=5.4$ at half-filling
to $U_{c}=5.8$ at $\delta=0.02$. In Fig. 4(d), $C_{12}$ increases slowly when
$U<U_{\text{Mott}}$, then rises dramatically as $U$ approaches
$U_{\text{Mott}}$, and finally grows progressively as $U$ goes to infinity;
while $C_{13}$ shows a maximum near $U_{\text{Mott}}$, and the reason is that
the hopping probability between the next-nearest neighbors falls faster than
the one between the nearest neighbors when $U$ is increased as demonstrated in
our previous work on a square lattice Zeng  _et al._ (2021). Unlike $C_{12}$
and $C_{13}$, both are positive for all $U$’s, $C_{14}$ is negative at
$\delta=0.02$ and its magnitude grows monotonically with the coupling
strength, indicating that the holon and doublon between the next-next-nearest
neighbors repulse each other, and this tendency is strengthened as $U$
increases. The negativity of $C_{14}$ can be understood as follows: according
to our previous work Zeng  _et al._ (2021), $\langle N_{1}D_{4}\rangle$ is
connected to the hopping probability between site $1$ and $4$ that decreases
dramatically with the growth of $U$ and $\delta$, i.e., in comparison to
$\langle N_{1}D_{4}\rangle>\langle N_{1}\rangle\langle D_{4}\rangle$ at half-
filling, $\langle N_{1}D_{4}\rangle<\langle N_{1}\rangle\langle D_{4}\rangle$
at $\delta=0.02$ and their difference will be widened by the increase of the
coupling strength. In Fig. 4(e), as shown by the blue line with $0<U<25$
(where the system within two-site approximation is in the PM state.) and the
black line, the difference of $\langle H_{n_{c}\text{-site}}^{S}\rangle/n_{c}$
in PM state between two- and six-site cluster approximations is much smaller
when $U>U_{\text{Mott}}$, denoting that the cluster size’s effect on the
properties of the system in the PM state is less important at large $U$ as the
system becomes more localized and the inter-site correlations are much weaker
in contrast to the weak-coupling limit. In Fig. 4(f), there exists an
inflection in $\langle D\rangle$ in the PM state around $U\approx 10$, meaning
the first-order Mott transition in the half-filled system turns into a
continuous crossover at finite dopings. For $U<8.0$, the double occupancy in
the AFM state is smaller than that in the PM state while the opposite is true
for $U>8.0$, bespeaking that the AFM at small $U$ is triggered by the
interaction potential gain while that in the large $U$ limit is not driven by
this mechanism. This picture can also be seen in Fig. 5, where exists a region
($5.4<U<8.0$) with $\Delta E_{U}^{(S)}<0$ and $\Delta E_{K}^{(S)}<0$,
signaling that the AFM in this region is supported by both the kinetic energy
and interaction potential gain.

Figure 5: The difference of the kinetic energy $\Delta E_{K}^{(S)}$ (red),
interaction potential $\Delta E_{U}^{(S)}$ (black), and their summation
$\Delta E_{\text{Total}}^{(S)}$ (blue) of the six-site cluster slave-spin
Hamiltonian Eq. (7d) between the AFM and PM states as function of $U$ at
$\delta=0.02$. Figure 6: (a) The quasiparticle weight $Z$, (b) the AFM energy
gap $\Delta_{\text{AFM}}/t$, (c) the staggered magnetization $M$, (d)–(f) the
holon-doublon correlator between the nearest neighbors $C_{12}$, the next-
nearest neighbors $C_{13}$ and the next-next-nearest neighbors $C_{14}$, (g)
the expectation value of the cluster slave-spin Hamiltonian, and (h) the
double occupancy $\langle D\rangle$ as function of $U$ at a series of doping
concentrations $\delta=$ 0.02 (red), 0.0375 (blue), 0.1 (green), 0.1531 (dark
green), 0.1725 (violet) in the AFM state obtained by the six-site cluster.

The dependence of the quantities discussed above upon the interaction strength
at various dopings of $\delta=0.02,\,0.0375,\,0.1,\,0.1531,\,0.1725$ obtained
through six-site cluster approximation are presented in Fig. 6, where
$\delta=0.1531$ is the critical doping for the AFM to PM phase transition at
$U=40$, and $\delta=0.1725$ is the maximum of $\delta_{M}$, i.e., the boundary
between the AFM and PM state [see Figs. 8 and 11]. The following results are
concluded: (i) In Fig. 6(a), apart from its monotonous diminishing with the
growing $U$ at all dopings, $Z$ increases with $\delta$ for all coupling
strengths, suggesting that the system with the increasing doping tends to be
metallic. When $U>U_{\text{Mott}}$, the quasi-particle residues tend
approximately to constants, manifesting that the properties in the AFM state
are controlled by the underlying Mott transition of the half-filled PM state.
(ii) In Fig. 6(b), $\Delta_{\text{AFM}}$’s at all dopings exhibit a maximum
around the crossover coupling strength $U_{c}(\delta)$ that grows with
$\delta$. Moreover, $\Delta_{\text{AFM}}$ decreases with increasing $\delta$
when $U<U_{c}(\delta)$, while displays a nonmonotonic behavior for
$U>U_{c}(\delta)$. Interestingly, $\Delta_{\text{AFM}}$ diminishes to zero as
$U$ increases at large dopings, e.g., $\delta=0.1531,\,0.1725$, and this
property has also been found in the square lattice Hubbard model Zheng and
Chan (2016); Zeng  _et al._ (2021). The reason is that the AFM is mostly
itinerant when $\delta$ approaches the critical dopings $\delta_{M}$ at which
the AFM state transits to PM state, and the itinerant AFM is not favored by
the overly enhanced localization for large $U$. (iii) In Fig. 6(c), the
staggered magnetization $M$ also rises rapidly, but as $\delta$ approaches
$\delta_{M}$, it declines gradually to zero with the increasing $U$, same as
the situation in $\Delta_{\text{AFM}}$. (iv) In Fig. 6(d)–(f), $C_{12}$
increases monotonically with $U$ at all dopings and diminishes as $\delta$
goes up, which makes it eligible to be an indicator of the magnitude of
correlations. There exists a maximum for $C_{13}$ around $U_{\text{Mott}}$,
which has been identified as one of features characterizing the strong
correlations Zeng  _et al._ (2021). However, $C_{14}$’s are negative for
nonzero dopings, suggesting that the holon and doublon repulse each other when
$\delta>0$, and its magnitude increases with either $U$ or $\delta$. (v) In
Fig. 6(g), compared to the PM state, $\langle H_{\text{6-site}}^{S}\rangle/6$
is suppressed dramatically as soon as AFM sets in and this effect is weakened
by the increasing $\delta$. Moreover, this quantity grows monotonically with
$\delta$ when $U<U_{\text{Mott}}$ while displays a nonmonotonic behavior when
$U>U_{\text{Mott}}$. (vi) In Fig. 6(h), $\langle D\rangle$ decreases with $U$
and tends to zero as $U$ goes to infinity. Interestingly, this quantity at
small $U$’s diminishes with the increase of $\delta$, while behaves
nonmonotonically at large $U$, reflected by the intersection of the red and
blue lines.

The compressibility of the system is defined as $\kappa=n^{-2}\partial
n/\partial\mu$. At $U=0$, it is calculated by using the non-interacting DOS
$D(\gamma)$, Eq. (14), via

$\displaystyle\kappa(\mu)=\Bigg{\\{}\begin{array}[]{ll}\frac{N(|\mu|)}{[\frac{3}{4}+\int_{-1}^{\mu}N(|\gamma|)d\gamma]^{2}}\;,&\;\;(0<\delta<\frac{1}{4})\\\
\frac{\tilde{N}(|\mu|)}{[\int_{-3}^{\mu}\tilde{N}(|\gamma|)d\gamma]^{2}}\;,&\;\;(\frac{1}{4}<\delta<1)\end{array}$
(31)

which is proportional to the free DOS. For $U>0$, $\kappa$ should evolve
simultaneously with the quasi-particle DOS which makes it adequate to indicate
the dependence of this quantity upon interaction. To illustrate the
$\delta$-dependence of the compressibility at various couplings, we plot
$\kappa$ as function of $\delta$ at $U=0,\,2,\,4,\,6$, and 8 in Fig. 7. For
$U<U_{\text{AFM}}$, the compressibility near van Hove singularity is
suppressed most drastically by interaction, while that at low energy (near the
Dirac points) remains very close to the non-interacting one, reflecting that
the Dirac cone structure is very robust, and the DOS of the quasi-particles is
transferred away from the van Hove singularity as $U$ increases. Furthermore,
at $U=4,\,6$, there exhibit a discontinuity at $\delta=\delta_{M}$ where the
AFM-to-PM phase transition occurs, and the one-sided peak of $\kappa$ as
$\delta$ approaches $\delta_{M}$ manifests that the system now is an itinerant
AFM metal. However, at $U=8$, there exist two consecutive discontinuities: (i)
between positively and negatively divergent $\kappa$; (ii) between negatively
divergent and positive small $\kappa$.

Figure 7: The compressibility $\kappa$ obtained from the six-site cluster
approximation as function of $\delta$ at $U=0,\,2,\,4,\,6,\,8$ (red, orange,
magenta, olive, brown, respectively), and the exact one at $U=0$ calculated
through Eq. (31) (cyan).

##  IV PHASE DIAGRAM

The staggered magnetization $M$ with $U$ and $\delta$ being its parameters is
plotted in Fig. 8, where the phase boundary between the AFM and PM states is
delineated by $\delta_{M}(U)$. Obviously, $M$ maximizes at small dopings and
large couplings. The phase boundary $\delta_{M}(U)$ shows a nonmonotonic
behavior upon $U$ that may be connected to the crossover of
$\Delta_{\text{AFM}}$ as $U$ increases. We also notice that $M$ saturates when
$U>U_{\text{Mott}}$ at small dopings, reflecting that the physical properties
in the AFM state are dominated by the underlying Mott transition in the half-
filled PM state.

Figure 8: The staggered magnetization $M$ as function of $U$ and $\delta$
obtained by six-site cluster, where $\delta_{M}(U)$ is the phase boundary
between the AFM and PM states.

The AFM energy gap $\Delta_{\text{AFM}}$ in the same parameter space is
plotted in Fig. 9, where the phase boundary $\delta_{M}(U)$ still holds,
denoting that there is a single direct and continuous phase transition from a
PM semi-metal to an antiferromagnet. The phase transition is in the Gross-
Neveu-Yukawa universality class, Eq. (23). An overall crossover between the
weak- and strong-coupling regimes can be observed in $\Delta_{\text{AFM}}$
when $U$ grows, at which $\Delta_{\text{AFM}}$ reaches its maximum, and the
coupling $U_{c}$ for this crossover is highly $\delta$-dependent, in contrast
to that in the square lattice case Zeng  _et al._ (2021). For $U>U_{c}$, the
maximum of $\Delta_{\text{AFM}}$ occurs at $\delta\sim 0.075$, leading to an
interesting vertical re-entrance behavior as $\delta$ increases, same as the
square lattice case Zeng  _et al._ (2021).

Figure 9: The AFM gap $\Delta_{\text{AFM}}$ as function of $U$ and $\delta$
obtained by six-site cluster with $\delta_{M}(U)$ separating the AFM and PM
states. Figure 10: The compressibility $\kappa$ as function of $U$ and
$\delta$ obtained by six-site cluster. The staircase is an artifact because of
discrete $U$’s adopted to calculate $\kappa$, i.e., $\Delta U=t$ when $U\leq
10t$ and $\Delta U=2t$ when $10t<U\leq 20t$, which can only be eliminated in
the $\Delta U\to 0$ limit. The blue region between $\delta_{\kappa}^{1}(U)$
and $\delta_{\kappa}^{2}(U)$ is characterized by $\kappa<0$, where
$\delta_{\kappa}^{1}(U)$ and $\delta_{\kappa}^{2}(U)$ are delineated by the
midpoints of the blue and red steps, respectively. The discontinuities in
$\kappa$ at these two phase boundaries are reflected in the color jumpings.

Combining Figs. 8, 9, and 10, a phase diagram in the $U$-$\delta$ plane
emerges in Fig. 11. In contrast to the square lattice case Zeng  _et al._
(2021), the crossover $U_{c}$ in the AFM state at which $\Delta_{\text{AFM}}$
is maximized, symbolled by the red squares, is shown to be highly
$\delta$-dependent. On the other hand, $U_{c}$ at $\delta=0$ is smaller than
$U_{\text{Mott}}$ (red triangle), implying that at half-filling, the coupling
strength separating the weak- and strong-coupling regimes is suppressed by
long range AFM correlations Zeng  _et al._ (2021). The blue region in this
figure enclosed by $\delta^{1}_{\kappa}(U)$ and $\delta^{2}_{\kappa}(U)$ is
characterized by $M\neq 0$ and $\kappa<0$ with $\delta^{1}_{\kappa}(U)$ being
the phase boundary between $M\neq 0$, $\kappa>0$ and $M\neq 0$, $\kappa<0$,
and $\delta^{2}_{\kappa}(U)$ between $M\neq 0$, $\kappa<0$ and $M=0$,
$\kappa>0$. The region with $\kappa<0$ is extremely small compared to the
square lattice Hubbard model Zeng  _et al._ (2021), and exists only in the
vicinity of the phase boundary between the AFM and PM state and at
intermediate $U$.

Figure 11: The $U$-$\delta$ phase diagram of the hexagon lattice Hubbard model
within the six-site cluster scheme. The critical coupling
$U_{\text{Mott}}=8.36$ for the Mott transition in the half-filled PM state is
marked by the red triangle. The crossover coupling $U_{c}(\delta)$ in the AFM
state at $\delta=$ 0, 0.02, 0.0375, 0.1, 0.1531 and 0.1725 are symbolled by
red squares, at which $\Delta_{\text{AFM}}$ reaches its maximum. The half-
filling case with $U>2.43$ is highlighted by the heavy blue line, in which the
system is an AFM insulator with $\kappa=0$.

##  V CONCLUSION

We have exploited the cluster slave-spin method to explore extensively the
ground state properties of the one-band honeycomb lattice Hubbard model with
$U$ and $\delta$ as its parameters. At half-filling, the first-order semi-
metal to insulator Mott transition in the PM state is revealed, characterized
by discontinuities and hystereses in all quantities considered at
$U_{\text{Mott}}=8.36t$ Hassan and Sénéchal (2013); Chen  _et al._ (2014);
Tran and Kuroki (2009). In the AFM state, a single direct and continuous phase
transition between PM semi-metal and AFMI at $U_{\text{AFM}}=2.43t$ is
substantiated, which belongs to the Gross-Neveu-Yukawa universality class
Otsuka  _et al._ (2016); Raczkowski  _et al._ (2020); Assaad and Herbut
(2013); Paiva  _et al._ (2005); Yamada (2016); Sorella and Tosatti (1992);
Sorella  _et al._ (2012); Ostmeyer  _et al._ (2020); Hassan and Sénéchal
(2013), precluding the existence of intermediate phases such as a spin liquid
state. At finite dopings, there exists an extended crossover between the weak-
and strong-coupling regimes in the AFM state at which the AFM energy gap
$\Delta_{\text{AFM}}$ reaches its maximum, and the AFM within this crossover
is driven by both the kinetic energy and interaction potential gain. The
interaction $U_{c}$ for this crossover is smaller than $U_{\text{Mott}}$ at
small $\delta$ or otherwise when $\delta$ is large, which is distinct from the
square lattice system where $U_{c}$ remains almost unchanged with large
dopings Zeng  _et al._ (2021).

Moreover, for the half-filled system, by analytically calculating the relation
between $M$ and $\Delta_{\text{AFM}}$ in the vicinity of the PM semi-metal to
AFM insulator transition, Eq. (III.1), we found that to the leading order, $M$
is linearly dependent on $\Delta_{\text{AFM}}$, in contrast to the square
lattice result that is proportional to
$\Delta_{\text{AFM}}(\ln{\Delta_{\text{AFM}}})^{2}$ (Ref. Zeng  _et al._ ,
2021). This difference is consistent with the vanishing non-interacting DOS at
Dirac points in the honeycomb lattice, compared to the van Hove singularity of
the free electron DOS at the Fermi surface in a half-filled square lattice.

Finally, an overall phase diagram in the $U$-$\delta$ plane is presented in
Fig. 11, the phase boundary $\delta_{M}(U)$ separating the AFM and PM phases
shows a nonmonotonic behavior with the increasing $U$, which is consistent
with the crossover behavior of $\Delta_{\text{AFM}}$. The phase boundary
between the AFM metal with $\kappa>0$ and the AFM insulator with $\kappa=0$
locates exactly at $\delta=0$. The region with $\kappa<0$ only exists in the
vicinity of the phase boundary between the AFM and PM state and at
intermediate coupling strengths, whose area is extremely small compared to the
counterpart in the square lattice Hubbard model Zeng  _et al._ (2021).

###### Acknowledgements.

We thank Shiping Feng, Xiong Fan and Yu Ni for many helpful discussions. This
work was supported by NSFC (Nos. 11774033 and 11974049), Beijing Natural
Science Foundation (No. 1192011), and the HSCC program of Beijing Normal
University.

## References

  * Hubbard (1963) J. Hubbard, Proc. R. Soc. Lond. A 276, 238 (1963).
  * Hubbard (1964) J. Hubbard, Proc. R. Soc. Lond. A 277, 237 (1964).
  * White  _et al._ (1989) S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, Phys. Rev. B 40, 506 (1989).
  * Qin  _et al._ (2021) M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, “The hubbard model: A computational perspective,” (2021), arXiv:2104.00064 [cond-mat.str-el] .
  * Meng  _et al._ (2010) Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Muramatsu, Nature 464, 847–851 (2010).
  * Hohenadler  _et al._ (2011) M. Hohenadler, T. C. Lang, and F. F. Assaad, Phys. Rev. Lett. 106, 100403 (2011).
  * Hohenadler  _et al._ (2012a) M. Hohenadler, T. C. Lang, and F. F. Assaad, Phys. Rev. Lett. 109, 229902 (2012a).
  * Hohenadler  _et al._ (2012b) M. Hohenadler, Z. Y. Meng, T. C. Lang, S. Wessel, A. Muramatsu, and F. F. Assaad, Phys. Rev. B 85, 115132 (2012b).
  * Zheng  _et al._ (2011) D. Zheng, G.-M. Zhang, and C. Wu, Phys. Rev. B 84, 205121 (2011).
  * Otsuka  _et al._ (2016) Y. Otsuka, S. Yunoki, and S. Sorella, Phys. Rev. X 6, 011029 (2016).
  * Raczkowski  _et al._ (2020) M. Raczkowski, R. Peters, T. T. Phùng, N. Takemori, F. F. Assaad, A. Honecker, and J. Vahedi, Phys. Rev. B 101, 125103 (2020).
  * Assaad and Herbut (2013) F. F. Assaad and I. F. Herbut, Phys. Rev. X 3, 031010 (2013).
  * Paiva  _et al._ (2005) T. Paiva, R. T. Scalettar, W. Zheng, R. R. P. Singh, and J. Oitmaa, Phys. Rev. B 72, 085123 (2005).
  * Yamada (2016) A. Yamada, International Journal of Modern Physics B 30, 1650158 (2016).
  * Sorella and Tosatti (1992) S. Sorella and E. Tosatti, Europhysics Letters (EPL) 19, 699 (1992).
  * Sorella  _et al._ (2012) S. Sorella, Y. Otsuka, and S. Yunoki, Scientific Reports 2 (2012), 10.1038/srep00992.
  * Ostmeyer  _et al._ (2020) J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Lähde, T. Luu, and C. Urbach, Phys. Rev. B 102, 245105 (2020).
  * Ostmeyer  _et al._ (2021) J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Lähde, T. Luu, and C. Urbach, “The antiferromagnetic character of the quantum phase transition in the hubbard model on the honeycomb lattice,” (2021), arXiv:2105.06936 [cond-mat.str-el] .
  * Hassan and Sénéchal (2013) S. R. Hassan and D. Sénéchal, Phys. Rev. Lett. 110, 096402 (2013).
  * Seki and Ohta (2012) K. Seki and Y. Ohta, “Quantum phase transitions in the honeycomb-lattice hubbard model,” (2012), arXiv:1209.2101 [cond-mat.str-el] .
  * Liebsch and Wu (2013a) A. Liebsch and W. Wu, Phys. Rev. B 87, 205127 (2013a).
  * Chen  _et al._ (2014) Q. Chen, G. H. Booth, S. Sharma, G. Knizia, and G. K.-L. Chan, Phys. Rev. B 89, 165134 (2014).
  * Arya  _et al._ (2015) S. Arya, P. V. Sriluckshmy, S. R. Hassan, and A.-M. S. Tremblay, Phys. Rev. B 92, 045111 (2015).
  * Honerkamp (2008) C. Honerkamp, Phys. Rev. Lett. 100, 146404 (2008).
  * Raghu  _et al._ (2008) S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang, Phys. Rev. Lett. 100, 156401 (2008).
  * Wang  _et al._ (1993) Y. R. Wang, J. Wu, and M. Franz, Phys. Rev. B 47, 12140 (1993).
  * Feng  _et al._ (1993) S. Feng, J. B. Wu, Z. B. Su, and L. Yu, Phys. Rev. B 47, 15192 (1993).
  * Feng  _et al._ (1994) S. Feng, Z. B. Su, and L. Yu, Phys. Rev. B 49, 2368 (1994).
  * Feng (2003) S. Feng, Phys. Rev. B 68, 184501 (2003).
  * Feng  _et al._ (2004) S. Feng, J. Qin, and T. Ma, Journal of Physics: Condensed Matter 16, 343 (2004).
  * Feng  _et al._ (2015) S. Feng, Y. Lan, H. Zhao, L. Kuang, L. Qin, and X. Ma, International Journal of Modern Physics B 29, 1530009 (2015).
  * Hassan and de’ Medici (2010) S. R. Hassan and L. de’ Medici, Phys. Rev. B 81, 035106 (2010).
  * Yu and Si (2012) R. Yu and Q. Si, Phys. Rev. B 86, 085104 (2012).
  * Lee and Lee (2017) W.-C. Lee and T.-K. Lee, Phys. Rev. B 96, 115114 (2017).
  * Zeng  _et al._ (2021) M.-H. Zeng, T. Ma, and Y.-J. Wang, “Phase diagram of the hubbard model on a square lattice: A cluster slave-spin study,” (2021), arXiv:2106.05457 [cond-mat.str-el] .
  * Gu  _et al._ (2013) Z.-C. Gu, H.-C. Jiang, D. N. Sheng, H. Yao, L. Balents, and X.-G. Wen, Phys. Rev. B 88, 155112 (2013).
  * Wang  _et al._ (2012) W.-S. Wang, Y.-Y. Xiang, Q.-H. Wang, F. Wang, F. Yang, and D.-H. Lee, Phys. Rev. B 85, 035414 (2012).
  * Nandkishore  _et al._ (2012) R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Nature Physics 8, 158 (2012).
  * Kotliar and Ruckenstein (1986) G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57, 1362 (1986).
  * Ma  _et al._ (2018) T. Ma, L. Zhang, C.-C. Chang, H.-H. Hung, and R. T. Scalettar, Phys. Rev. Lett. 120, 116601 (2018).
  * Liebsch and Wu (2013b) A. Liebsch and W. Wu, Phys. Rev. B 87, 205127 (2013b).
  * Bleinstein and Handelsman (1986) N. Bleinstein and R. A. Handelsman,  _Asymptotic Expansions of Integrals_ (Dover Publications, New York, 1986).
  * Ogawa  _et al._ (1975) T. Ogawa, K. Kanda, and T. Matsubara, Progress of Theoretical Physics 53, 614 (1975).
  * Abram  _et al._ (2013) M. Abram, J. Kaczmarczyk, J. Jedrak, and J. Spalek, Phys. Rev. B 88, 094502 (2013).
  * Tran and Kuroki (2009) M.-T. Tran and K. Kuroki, Phys. Rev. B 79, 125125 (2009).
  * Vollhardt (1984) D. Vollhardt, Rev. Mod. Phys. 56, 99 (1984).
  * Zheng and Chan (2016) B.-X. Zheng and G. K.-L. Chan, Phys. Rev. B 93, 035126 (2016).


# The lowest excited states of 14C and 14O nuclei within a five-cluster model.

B. E. Grinyuk<EMAIL_ADDRESS>Bogolyubov Institute for Theoretical
Physics,
Kiev 03143, Ukraine  D. V. Piatnytskyi<EMAIL_ADDRESS>V. S. Vasilevsky
<EMAIL_ADDRESS>
###### Abstract

We study the ground and the first excited $0^{+}$ states of two mirror nuclei
14C and 14O within a five-cluster model (three alpha-particles and two extra
nucleons) with the use of high accuracy variational approach with Gaussian
bases. The first excited $0^{+}$ state of these nuclei is shown to be
connected with a change of the structure of the two-nucleon subsystem moving
in the field of the 12C cluster. The density distributions, electric form
factors (both the elastic and transition ones), pair correlation functions,
and the momentum distributions of the constituent particles are found. The
most probable shapes of five-cluster systems are revealed.

###### keywords:

Cluster model , cluster-cluster interaction , Coulomb interaction , five-
cluster model , mirror nuclei , resonance states

###### PACS:

21.60.Gx , 21.60.-n , 24.10.-i

††journal: Nuclear Physics A

##  1 Introduction

Cluster models are powerful tools to study a large verity of atomic nuclei and
nuclear reactions. They constantly demonstrate their ability to describe and
to interpret numerous phenomena in atomic nuclei observed both in scientific
laboratories and the Universe. They also allowed one to predict interesting
properties of nuclei or nuclear reactions. New developments of cluster models,
their effectiveness and new results obtained are thoroughly discussed in the
following reviews [1, 2, 3, 4, 5, 6, 7, 8].

To provide more adequate description of nuclear structure and nuclear
processes, cluster models are subject for a permanent improvement and
advancing, which, for instance, include more and more clusters and application
to study nuclei with a large excess of protons or neutrons. By including more
clusters and consequently more channels into consideration one faces with
bulky calculations, which are also complicated by the Pauli principle. There
were formulated some methods to reduce calculations associated with the Pauli
principle. One of the most famous methods is the orthogonality condition
method suggested by S. Saito [9, 10]. This method significantly reduces
computational efforts and provides a reasonable description of nuclear
structure and reactions. There are also further simplifications when clusters
are considered to be structureless particles and cluster-cluster interactions
are modelled by some effective interactions. Such a type of cluster models
does not require tremendous calculations and allows one to reveal some
interesting phenomena which can be confirmed in more sophisticated
calculations.

We consider the structure of 14C and 14O within a five-body model. These
nuclei are presented as a configuration $\alpha+\alpha+\alpha+N+N$. This five-
cluster configuration allows one, in principle, to analyze the structure of
14C and 14O nuclei as consisting from different subsystems. First, such a
cluster configuration allows one to consider these nuclei as a core 12C
comprised of three alpha-particles and two valence nucleons (neutrons or
protons). Thus our model accounts for a more flexible description of the 12C
inside the nuclei 14C and 14O, since the shape of triangle connecting the
centers of mass of three alpha-particles can be varied by valence neutrons
(protons). Second, one can consider the nucleus 14C (14O) as a “core” 6He(6Be)
surrounded by two alpha-particles. In this case, triangle connecting center of
mass of an alpha-particle and two valence neutrons can be essentially deformed
by two valence alpha-particles.

In the present paper we are going to study the ground and the first excited
$0^{+}$ states of 14C and 14O nuclei. Different types of the cluster-cluster
correlations will be calculated and analyzed in order to reveal the nature of
states under investigation.

The nucleus 14C attracts large attention from both theoretical (see Refs. [11,
12, 13, 14, 15, 16, 17, 18]) and experimental (see Refs. [19, 20, 21, 22, 23,
24, 25, 26]) points of view. Some of these investigations are devoted to
searching of linear configurations of clusters in these nuclei. The
theoretical investigations are mainly performed within the framework of the
antisymmetrized molecular dynamics (AMD).

Preliminary results of the present calculations carried out for the ground
state only, are reported in Refs. [15] and [18].

Experimental spectrum of the $0^{+}$ states of 14C is presented in Table 1 and
is taken from Ref. [19].

Table 1: Experimental spectrum of the lowest $0^{+}$ states and energies of two- and three-cluster thresholds in 14C and 14O. 14C |  $E^{*}\\!\left(0^{+}\\!\right)\\!-\\!E_{0}\\!\left(0^{+}\\!\right)$, MeV | 6.589 | 9.746 |  |
---|---|---|---|---|---
|  $E_{threshold}$, MeV | 8.177 | 12.012 | 13.123 | 19.424
| Decay channel |  13C$+n$ |  10Be$+\alpha$ |  12C$+n\\!+\\!n$ |  6He$+\alpha\\!+\\!\alpha$
14O |  $E^{*}\\!\left(0^{+}\\!\right)\\!-\\!E_{0}\\!\left(0^{+}\\!\right)$, MeV | 5.920 |  |  |
|  $E_{threshold}$, MeV | 4.628 | 10.118 | 6.572 |
| Decay channel |  13N$+p$ |  10C$+\alpha$ |  12C$+p\\!+\\!p$ |

In Figure 1, we show the spectrum of the $0^{+}$ states in 14C and the
dominant two- and three-cluster decay thresholds. This figure suggests certain
hierarchy of two- and three-cluster channels in 14C. The most important two-
body channel is the 13C$+n$ channel. One may conjecture that this channel is
responsible or, in other words, dominant for formation of the ground and first
excited $0^{+}$ states, and the $0^{+}$ resonance states which lies just above
the threshold of this channel. Second of importance is the 10Be$+\alpha$
channel and then the three-body 12C$\,+\,n+n$ channel. In the 14O nucleus,
there is only one excited $0^{+}$ state which lies between the two-body
13N$+p$ and three-body 12C$+p+p$ thresholds.

Figure 1: Spectrum of the lowest $0^{+}$ states of the nuclei 14C (on the
left) and 14O (on the right), and two- and three-cluster dominant decay
channels. Energies in MeV are accounted from the ground states of the nuclei.

The layout of the present paper is the following. In Section 2, we shortly
explain the main ideas and key elements of the five-cluster model used. The
key elements include an alpha-alpha potential, the nucleon-nucleon one in the
singlet state, and alpha-nucleon potential as well. In Section 3, we formulate
the method of finding the wave functions for the ground and excited states of
14C and 14O nuclei within the five-cluster approach. A detailed analysis of
the results obtained is also presented in this section.

##  2 Five-cluster model for 14O and 14C nuclei

In this section, we formulate our five-cluster model and describe the method
of calculations.

For 14O nucleus considered as three $\alpha$-particles and two extra protons,
we start from a model Hamiltonian [15]:

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum\limits_{i=1}^{2}\frac{\mathbf{p}_{i}^{2}}{2m_{p}}+\sum\limits_{i=3}^{5}\frac{\mathbf{p}_{i}^{2}}{2m_{\alpha}}+U_{pp}\left(r_{12}\right)+\sum\limits_{j>i=3}^{5}\hat{U}_{\alpha\alpha}\left(r_{ij}\right)$
(1) $\displaystyle+$
$\displaystyle\sum\limits_{i=1}^{2}\sum\limits_{j=3}^{5}\hat{U}_{p\alpha}\left(r_{ij}\right)+\sum\limits_{j>i=1}^{5}\frac{Z_{i}Z_{j}e^{2}}{r_{ij}}.$

Within our approach, we make some approximations to simplify the model. In
particular, considering the $0^{+}$ states, we take into account only central
part of interactions between particles, and this means that the total spin of
the nucleus is a good quantum number. We restrict ourselves with the zero
value of the total spin, which in our model coincides with the spin of two
valence nucleons. It is then naturally to involve nucleon-nucleon interaction
in the singlet state only. An account of an interaction of two nucleons in the
triplet state may be important for highly excited states of the nuclei under
consideration. We also neglect the spin-orbit part of the $\alpha–N$
interaction assuming that it plays minor role in formation of the $0^{+}$
states in 14C and 14O nuclei. Starting with the Hamiltonian (1) for “point-
like” particles, we nevertheless assume that the constituent particles have a
finite size, and that they have definite charge and density distributions
known from the experiment, and we take them into account within the Helm
approximation [37] after the solution of the five-body problem (1) to study
the charge and mass distributions, as well as the charge form factors of the
nuclei. Note also that the Coulomb repulsion in (1) could be chosen in a
modified form at short distances with an account of finite size of particles.
But, as it follows from comparison of mirror nuclei 10Be and 10C [28, 29], as
well as 14C and 14O [15], the Coulomb interaction between particles does not
play the decisive role in the structure of the nuclei as compared to the
nuclear forces. Moreover, due to the phase volume $\sim r^{2}dr$ in the
integrals, a contribution of the Coulomb repulsion at short distances into the
matrix elements is small. Thus, our choice of simple form of Coulomb repulsion
at short distances does not play essential role. More important is the fact
that, in our simple model, we ignore the relativistic effects in the kinetic
energy operator and tensor nuclear forces. We also do not take into account
other possible cluster configurations in the nuclei essentially important for
highly excited states. The role of all these effects needs a separate
consideration and analysis.

One can explicitly extract the kinetic energy of the center of mass of the
nucleus from expression (1), but there is no sense to carry out this trivial
procedure since we use the translation invariant wave functions (see below)
giving automatically the zero kinetic energy of the center of mass of the
system. The Hamiltonian for 14C nucleus within the five-cluster model (three
$\alpha$-particles plus two extra neutrons) is very similar to (1), but with
neutron mass $m_{n}$ instead of the proton one in the kinetic energy, and less
strong Coulomb repulsion (due to the fact that $Z_{1}=Z_{2}=0$ for extra
neutrons in 14C nucleus, while $Z_{1}=Z_{2}=1$ for extra protons in 14O, and
in both cases $Z_{3}=Z_{4}=Z_{5}=2$ for $\alpha$-particles). We assume the
nuclear potential between the extra protons $U_{pp}$ (in 14O) to be the same
as the one between extra neutrons $U_{nn}$ (in 14C). This assumption is in
concordance with the charge independence of nuclear forces. This interaction
potential in the singlet state proposed in [27] is chosen in the form

$U_{nn}\left(r\right)=U_{pp}\left(r\right)=\sum\limits_{k=1}^{3}V_{k}\cdot\exp\left(-\left(\frac{r}{R_{k}}\right)^{2}\right),$
(2)

where $V_{1}=952.15$ MeV, $R_{1}=0.44$ fm; $V_{2}=-79.39$ MeV, $R_{2}=0.959$
fm; and $V_{3}=-37.89$ MeV, $R_{3}=1.657$ fm. The potential (2) leads to a
description of the low-energy neutron-neutron phase shift in the singlet state
with reasonable accuracy, including the low-energy scattering parameters. It
was successfully used for extra nucleon interaction in studying 6He [27], 10Be
and 10C [28, 29] nuclei. The same potential was used for a study of the ground
state structure of 14C and 14O nuclei [15].

As to the interaction potential between an extra neutron and an
$\alpha$-particle, we use the generalized type potential with local
(attraction) and non-local (separable repulsion) terms. This type of
potentials [30, 31] enables one to concord the experimental phase shifts
(having difference in $\pi$ between their values at zero and infinite
energies) and the fact that neutron and $\alpha$-particle have no bound
states. This is possible due to a generalized Levenson’s theorem for such type
of potentials [31]. Parameters of this potential are given in [15]. The
$p\alpha$-interaction potential is almost the same as the $n\alpha$ one, but
with a little bit modified parameters [15]. The small difference in parameters
can be justified by the fact that, within an $\alpha$-particle, there is a
little bit different (of greater radius) distribution of protons as compared
to that of neutrons, mainly due to the Coulomb interaction [32].

For $\alpha\alpha$-interaction, we also use the generalized potentials with a
local part containing repulsion and attraction, and with additional non-local
(separable) repulsion. Parameters of the potential are given in [15]. The
total set of potentials used in Hamiltonian (1) is selected in such a way that
to reproduce experimental values (binding energy and charge radius) of the
ground state of 14C nucleus simultaneously with the corresponding phase shifts
(at least with qualitative accuracy) for corresponding pair collisions of
particles constituting the nucleus within the five-particle model. As to 14O
nucleus, its experimental charge radius is unknown (see our theoretical
estimation [15]), thus we adjusted the potentials to reproduce only the
binding energy of 14O.

The structure functions of the ground state of the both 14C and 14O nuclei
within the five-particle model were studied in detail in [15]. Now we are
going to study the first excited 0+ state of these nuclei within the same
variational method with Gaussian basis [33, 34, 35]. A study of higher excited
states needs some other approaches since they are higher than the first
13C$+n$ threshold. Moreover, treating the nuclei as consisting of
$\alpha$-particles ignores a contribution of their excitation or breakup
processes, and thus it has sense only for the lowest states of 14C and 14O
nuclei to be considered within our model.

The spectrum and wave functions of the system are found with high accuracy
within the variational method with the use of Gaussian basis [33, 34, 35, 36].
For $J^{\pi}=0^{+}$ states, the coordinate part of the wave functions of the
five-particle system can be found in the form

$\Phi=\hat{S}\sum_{k=1}^{K}C_{k}\varphi_{k}\equiv\hat{S}\sum_{k=1}^{K}C_{k}\,\exp\left(-\sum_{j>i=1}^{5}a_{k,ij}\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)^{2}\right),$
(3)

where $\hat{S}$ is the symmetrization operator with respect to identical
$\alpha$-particles, as well as to identical extra nucleons (since, according
to our assumption, they are in the singlet state, and their spin function and
thus the total function are antisymmetric with respect to the extra nucleons
permutation). The spin part of the total wave function is not indicated
explicitly for brevity. In all the calculations of matrix elements it gives
simply a unitary multiplier. In (3), $K$ is a dimension of the basis. The
coefficients $C_{k}$ can be found within the Galerkin method from the
algebraic system of equations

$\sum_{m=1}^{K}C_{m}\left\langle\hat{S}\varphi_{k}\left|\hat{H}-E\right|\hat{S}\varphi_{m}\right\rangle=0,\quad
k\,=\,0,1,...,K.$ (4)

Matrix elements for the Hamiltonian (1) are known to have an explicit form due
to the Gaussian functions of the basis. Within the complete basic set, at
$K\rightarrow\infty$, the system (4) is equivalent to the Schrödinger
equation. An advantage of the Gaussian bases lies in the fact that, at a
definite $K$, one can achieve better accuracy by varying the non-linear
parameters $a_{k,ij}$. This needs more time for calculations, but results in
the best accuracy for the energies and wave functions at a given $K$. The most
efficient method to save the time of calculations is the random choice of the
parameters $a_{k,ij}$ and the increasing of $K$ (“variational method with
stochastic Gaussian bases” [33, 34]). We used the both approaches. As a
result, we have energies and the wave functions of the ground and first
excited $0^{+}$ states of 14C and 14O nuclei within the five-particle model
with sufficiently high accuracy. After calculation of r.m.s. radii and density
distributions for “point-like” particles, we calculate the charge r.m.s. radii
and charge density distributions taking into account [15, 18] the size of
$\alpha$-particles and that of extra nucleons within the Helm approximation
[37]. Only after this correction, we compare the calculated charge radii with
the experimental data. The same Helm correction is used for calculation of
both elastic and transition electric form factors (see below).

##  3 Results and Discussion

In this section we present and discuss main results of our investigations.

In Table 2, results for energies and charge r.m.s. radii are given for both
14C and 14O nuclei calculated for the ground and first excited $0^{+}$ states.
The charge radii are found within the Helm approximation. Experimental values
of the ground state energies of the both nuclei, and charge r.m.s. radius of
14C nucleus, are reproduced with high accuracy [15] by adjusting the
parameters of potentials. The result $2.415$ fm for the charge radius of 14O
nucleus is our prediction [15].

The first excited $0^{+}$ state energy of 14O nucleus is calculated with
$K=500$ functions of the Gaussian basis, and it is seen to be in a good
coincidence with the experiment. The similar result for 14C nucleus is
obtained with $K=600$ Gaussian functions. It is in a qualitative agreement
with the experimental energy of the first $0^{+}$ excited state. We see that
the set of potentials proposed in [15] used for 14O nucleus appeared to be
more happy choice to treat the excited $0^{+}$ state than that used for 14C
one. The results for charge r.m.s. radii of the excited state are seen to be a
little bit greater than the ones for the ground state. Note also, that for the
ground state we have $R_{ch,0}(^{14}$O$)<R_{ch,0}(^{14}$C$)$, since the extra
protons of 14O nucleus move mainly inside the 12C cluster [15]. For the
considered excited states, one has more natural situation
$R_{ch,1}(^{14}$O$)>R_{ch,1}(^{14}$C$)$.

Table 2: Experimental (E) and calculated (T) energies of the ground and the first excited $0^{+}$ states and corresponding charge radii for 14C and 14O nuclei. |  14C (T) |  14C (E) |  14O (T) |  14O (E)
---|---|---|---|---
$E_{0}$, MeV | -20.398 | -20.3977 | -13.845 | -13.8445
$R_{ch}$, fm | 2.500 | 2.496 [38], 2.503 [39] | 2.415 | –
$E_{1}$, MeV | -13.223 | -13.808 | -7.932 | -7.925
$R_{ch}$, fm | 2.582 | – | 2.589 | –

More information about the spatial structure of the excited states under
consideration can be obtained from the analysis of r.m.s. relative distances
between constituting particles of the nuclei. These distances
($r_{\alpha\alpha}$ between $\alpha$-particles, $r_{NN}$ – between extra
nucleons, and $r_{N\alpha}$ – between an extra nucleon and an
$\alpha$-particle) are presented in Table 3. For comparison, the results [15]
for the ground state are given too. The general increase of the relative
distances between particles for the excited state is observed as compared to
the ones for the ground state. This fact is in concordance with the fact that,
generally speaking, the less is the binding energy, the greater are the
distances between particles. The greater relative distances mean better
substantiation of the five-particle model for the excited states considered,
since particles of finite size are less overlapping in these states then in
the ground state.

In the same Table 3, we present also the r.m.s. radii (the r.m.s. distances
from the center of mass of a nucleus) calculated for “point-like”
$\alpha$-particles ($R_{\alpha}$) and extra nucleons ($R_{N}$). These values
calculated directly are known to be connected with r.m.s. relative distances
by the identity [18]:

$R_{i}^{2}=\frac{1}{M^{2}}\left(\left(M-m_{i}\right)\sum_{\footnotesize{\begin{array}[]{c}{j\neq
i}\end{array}}}m_{j}r_{ij}^{2}-\sum_{\footnotesize{\begin{array}[]{c}{j<k}\\\
{\left(j\neq i,k\neq i\right)}\end{array}}}m_{j}m_{k}r_{jk}^{2}\right),$ (5)

where $M=\sum\limits_{j}m_{j}$ is the total mass of the system.

Table 3: Relative r.m.s. distances between particles, and r.m.s. radii in 14C and 14O nuclei (in fm), calculated within the five-particle model for “point-like” particles. The last column presents the mass radii $R_{m}$ with an account of the size of particles within the Helm approximation. | $r_{\alpha\alpha}$ | $r_{NN}$ | $r_{N\alpha}$ | $R_{\alpha}$ | $R_{N}$ | $R_{m}$
---|---|---|---|---|---|---
14C (ground $0^{+}$) | 3.189 | 2.621 | 2.667 | 1.852 | 1.786 | 2.433
14C (1st excited $0^{+}$) | 3.367 | 3.211 | 3.197 | 1.961 | 2.328 | 2.568
14O (ground $0^{+}$) | 3.239 | 2.732 | 2.750 | 1.882 | 1.865 | 2.461
14O (1st excited $0^{+}$) | 3.412 | 3.305 | 3.265 | 1.991 | 2.387 | 2.596

Now consider the density distributions of particles in the excited state in
comparison with those for the ground state of the nuclei 14C and 14O. In Fig.
2, density distributions multiplied by $r^{2}$ are shown. It is clearly seen
that, for the ground state of the both nuclei, the extra nucleons move mainly
inside the 12C cluster, while it is not so for the excited state. That is why
the charge density distribution calculated for the ground state of 14O nucleus
has maximum at shorter distances than the similar distribution for 14C does
(since there are charged extra protons in 14O nucleus moving more close to the
center of the nucleus than the $\alpha$-particles do, while in 14C nucleus,
only $\alpha$-particles make a contribution to the charge distribution). This
results in mentioned above inequality $R_{ch,0}(^{14}$O$)<R_{ch,0}(^{14}$C$)$.
And this effect is not revealed for the excited states.

Figure 2: Density distributions (multiplied by $r^{2}$) of $\alpha$-particles
and extra nucleons in the ground ($a$, $c$) and first excited ($b$, $d$)
$0^{+}$ states of 14C and 14O nuclei calculated for “point-like” particles.
The calculated charge density distribution (with an account of finite size of
constituent particles in Helm approximation) is depicted by thick line.

In Fig.3, we depict the pair correlation functions for the same nuclei in the
ground and the first excited states. The obtained distributions show that, in
the excited state, the average distances of extra nucleons from
$\alpha$-particles are enlarged as compared to those in the ground state. This
fact is in complete concordance with the results obtained for r.m.s. distances
(see Table 3).

Figure 3: Pair correlation functions for the ground ($a$, $c$) and first
excited ($b$, $d$) $0^{+}$ states of 14C and 14O nuclei.

Consider the momentum distributions of extra nucleons and $\alpha$-particles
in the ground and first excited $0^{+}$ states of the both nuclei. In Fig. 4,
we compare the momentum distributions of the 14C and 14O nuclei. The dashed
lines depict the momentum distributions for the ground state of these nuclei
[15]. It is clearly seen for the excited state momentum distributions to have
a somewhat greater low-energy (or low-momenta) values than those in the ground
state. This is natural for the momentum distributions of particles of less
bound systems as compared to more strongly bound ones. We also note that the
momentum distributions of extra nucleons in the excited state have dips
shifted towards the zero momenta as compared to the ground state. And the
change of regimes of decreasing of the momentum distribution profiles for
$\alpha$-particles is more clearly seen for the excited state than for the
ground state, and it is observed at lower momenta.

Figure 4: Momentum distributions in the ground (dashed lines) and first
excited (solid lines) $0^{+}$ states of 14C and 14O nuclei ($a$ and $b$,
respectively). Curves $1$ – for $\alpha$-particles, and curves $2$ – for extra
nucleons.

Generally speaking, the mirror nuclei 14C and 14O are seen to have very
similar spatial structure functions and the momentum distributions (at least
from qualitative point of view). That is why, for brevity, we consider below
some details of the calculated probability density distributions only for 14C
nucleus.

It is impossible to depict directly the squared wave function depending on too
large number of variables which is proper to the five-particle system. That is
why we introduce the values being the squared wave function partially
integrated over some of the variables and thus depending only on the rest few
of them. In particular, consider the following probability density
distribution

$P_{NN}\left(r,\rho,\theta\right)\equiv
r^{2}\rho^{2}\left\langle\Phi\left|\delta\left(\mathbf{r}-\mathbf{r}_{NN}\right)\delta\left(\mathbf{\rho}-\mathbf{\rho}_{\left(NN\right)\left(3\alpha\right)}\right)\right|\Phi\right\rangle,$
(6)

where $r$ is the distance between the extra nucleons, $\rho$ is the distance
between the centers of mass of the two-nucleon subsystem and the three-alpha
(12C) cluster, and $\theta$ is an angle between the vectors $\mathbf{r}$ and
$\mathbf{\rho}$. In Fig. 5, we show this value at fixed $\theta=0^{\circ}$ and
$\theta=90^{\circ}$ in case of 14C nucleus. Very similar pictures can be drawn
for 14O nucleus since an additional Coulomb repulsion between extra protons
themselves and that between the protons and $\alpha$-particles play minor role
as compared to the nuclear forces. It is obvious that two peaks [18] observed
for the probability density (6) at $\theta=0^{\circ}$ in the ground $0^{+}$
state change their positions in the first $0^{+}$ excited state. In
particular, the dineutron cluster starts to move around the 12C cluster at
greater distances of the order of $2$ fm, while a configuration corresponding
to two separate neutrons, vice versa, approaches to the center of mass of the
12C cluster. At $\theta=90^{\circ}$, the second peak does not reveal itself,
but it is seen once more that the two neutrons move further from the center of
mass of the 12C cluster, and it is seen that the size of the two extra
neutrons cluster in the excited state becomes greater.

We may consider another probability density,

$P_{\alpha\alpha}\left(r,\rho,\theta\right)\equiv
r^{2}\rho^{2}\left\langle\Phi\left|\delta\left(\mathbf{r}-\mathbf{r}_{\alpha\alpha}\right)\delta\left(\mathbf{\rho}-\mathbf{\rho}_{\left(\alpha\alpha\right)\left({}^{6}He\right)}\right)\right|\Phi\right\rangle,$
(7)

where $r$ is now a distance between two $\alpha$-particles, and $\rho$ is the
distance between the center of mass of the two $\alpha$-particle cluster and
the center of mass of the rest particles (of the 6He cluster). In Fig. 6, we
show this value for the same nucleus in the ground $0^{+}$ and the first
$0^{+}$ excited states. It is obvious from the pictures that the
configurations described by the value (7) almost do not change their positions
and form, on the contrary to the previous value (6) depicted in Fig.5. This
means that the 12C cluster is practically not excited in the $0^{+}$ first
excited state, while the excitation is connected with the two-nucleon
subsystem of 14C nucleus. We believe that first this effect was demonstrated
for 14C nucleus in [11] within an approach with the account of all nucleon
degrees of freedom.

Figure 5: Probability density distribution $P_{NN}\left(r,\rho,\theta\right)$
for the ground ($a$, $c$) and first excited ($b$, $d$) $0^{+}$ states of 14C
nucleus, at fixed $\theta=0^{\circ}$ ($a$,$b$) and $\theta=90^{\circ}$
($c$,$d$). Figure 6: Probability density distribution
$P_{\alpha\alpha}\left(r,\rho,\theta\right)$ for the ground ($a$, $c$) and
first excited ($b$, $d$) $0^{+}$ states of 14C nucleus.

We also present here the transition electric form factors of 14O and 14C
nuclei. The elastic form factors of these nuclei and their charge
distributions in the ground state are discussed in detail in [15]. Now we
calculate the electric transition form factors for both of these nuclei using
the wave functions of the ground $\Phi_{0}$ and the first excited $\Phi_{1}$
states found above.

We remind that, for a “point-like” charged particle $k$ of the system
described by the wave function $\Phi_{0}$, a contribution to the elastic form
factor $F_{k,00}\left(q\right)$ is known to be

$F_{k,00}\left(q\right)\equiv\int
e^{-i\left(\mathbf{q}\mathbf{r}\right)}n_{k}\left(r\right)d\mathbf{r}=\left\langle\Phi_{0}\left|e^{-i\left(\mathbf{q}\left(\mathbf{r}_{k}-\mathbf{R}_{c.m.}\right)\right)}\right|\Phi_{0}\right\rangle,$
(8)

where $n_{k}\left(r\right)$ is the probability density distribution of the
particle $k$. For the transition electric form factor between the states
$\Phi_{0}$ and $\Phi_{1}$, the same charged “point-like” particle makes a
contribution $F_{k,01}\left(q\right)$ calculated in a similar way:

$F_{k,01}\left(q\right)=\left\langle\Phi_{0}\left|e^{-i\left(\mathbf{q}\left(\mathbf{r}_{k}-\mathbf{R}_{c.m.}\right)\right)}\right|\Phi_{1}\right\rangle.$
(9)

An account of the contributions of all the charged particles results in a sum
over all such particles with weights proportional to their charges (the total
sum of weights being unity). In addition, one has to take into account that
the particles of the system are not “point-like”, and they have their own form
factors. Within the Helm approximation [37], a corresponding contribution of a
given particle into the total form factor should be multiplied by the own
charge form factor of these particle. In particular, for the transition
electric form factor of 14C nucleus we have (in the Helm approximation,
neglecting the contribution of extra neutrons)

$F_{{}^{14}\text{C},01}\left(q\right)=F_{\alpha,01}\left(q\right)\cdot
F_{\text{ch},^{4}\text{He}}\left(q\right),$ (10)

where $F_{\text{ch},^{4}\text{He}}\left(q\right)$ is the experimental form
factor of 4He nucleus [40]. For the transition form factor of the 14O nucleus,
one has two terms originating from $\alpha$-particles and extra protons
contributions:

$F_{{}^{14}\text{O},01}\left(q\right)=\frac{3}{4}F_{\alpha,01}\left(q\right)\cdot
F_{\text{ch},^{4}\text{He}}\left(q\right)+\frac{1}{4}F_{p,01}\left(q\right)\cdot
F_{\text{ch},p}\left(q\right),$ (11)

where $F_{\text{ch},p}\left(q\right)$ is the experimental charge form factor
of the proton [41].

An advantage of the Helm approximation lies in its simple and clear
formulation. But, unfortunately, one should take into account, that at large
momentum transfer corresponding to short distances in the coordinate space the
Helm approximation becomes not accurate enough, because it does not take into
account the Pauli principle and exchange effects between overlapping
particles. Thus our results for the form factors presented below (as well as
the results for elastic form factors obtained earlier [15]) are reliable at
not very high momentum transfer. In Fig. 7, the transition electric form
factor of 14C nucleus is presented calculated within approximation (10), while
in Fig. 8 we give the result for 14O obtained within (11). The dip observed in
Fig. 7 is situated not far from the dip [15] inherent to the elastic form
factor of 14C nucleus. The presence of such a dip is explained by the property
of the first multiplier in the formula (10) for the form factor of 14C
nucleus. On the other hand, the presence of the second term in the expression
(11) for the form factor of 14O nucleus “smooths out” the dip proper to the
first term. The same was observed for the elastic form factors [15].

Figure 7: Transition electric form factor for 14C nucleus. Figure 8:
Transition electric form factor for 14O nucleus.

##  4 Conclusions

To conclude, we note that the main structure functions (density distributions,
form factors, pair correlation functions, and momentum distributions) are
studied for the mirror nuclei 14C and 14O within a five-cluster model (three
alpha-particles and two extra nucleons). All the structure functions for the
ground and the excited states are compared between themselves, and a
conclusion is drawn that the first excited $0^{+}$ state of the both
considered nuclei looks like an excitation of the two-nucleon subsystem moving
in the field of the 12C cluster. The both nuclei are shown to have very
similar spatial structure, an additional Coulomb repulsion in 14O playing
secondary role. But the electric form factors of these nuclei, due to extra
protons present in 14O, are essentially different (this is true both for the
elastic and transition form factors).

## Acknowledgement

One of the authors (B.E.G.) would like to express appreciation to the European
Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT*) and
the Fondazione Bruno Kessler for financial support and hospitality within the
ECT* “Ukraine Fellowship Project”. This work was supported in part by the
Program of Fundamental Research of the Physics and Astronomy Department of the
National Academy of Sciences of Ukraine (Projects No. 0122U000886 and No.
0122U000889).

## References

  * [1] P. Descouvemont, D. Baye, The R-matrix theory, Rep. Prog. Phys. 73 (3) (2010) 036301.  arXiv:1001.0678, doi:10.1088/0034-4885/73/3/036301.
  * [2] H. Horiuchi, K. Ikeda, K. Katō, Recent Developments in Nuclear Cluster Physics, Prog. Theor. Phys. Suppl. 192 (2012) 1–238.  doi:10.1143/PTPS.192.1.
  * [3] Y. Funaki, H. Horiuchi, A. Tohsaki, Cluster models from RGM to alpha condensation and beyond, Prog. Part. Nucl. Phys. 82 (2015) 78–132.  doi:10.1016/j.ppnp.2015.01.001.
  * [4] C. Beck, From the stable to the exotic: clustering in light nuclei, arXiv e-prints (2016) arXiv:1603.02574arXiv:1603.02574.
  * [5] C. Beck, Clusters in light stable and exotic nuclei, arXiv e-prints (2018) arXiv:1812.08013arXiv:1812.08013.
  * [6] M. Freer, H. Horiuchi, Y. Kanada-En’yo, D. Lee, U.-G. Meißner, Microscopic clustering in light nuclei, Rev. Mod. Phys. 90 (3) (2018) 035004\.  arXiv:1705.06192, doi:10.1103/RevModPhys.90.035004.
  * [7] B. Zhou, Y. Funaki, H. Horiuchi, A. Tohsaki, Nonlocalized clustering and evolution of cluster structure in nuclei, Frontiers of Physics 15 (1) (2020) 14401.  arXiv:1905.00788, doi:10.1007/s11467-019-0917-0.
  * [8] T. Myo, K. Katō, Complex scaling: Physics of unbound light nuclei and perspective, Prog. Theor. Exp. Phys. 2020 (12) (2020) 12A101.  arXiv:2007.12172, doi:10.1093/ptep/ptaa101.
  * [9] S. Saito, Interaction between Clusters and Pauli Principle, Prog. Theor. Phys. 41 (3) (1969) 705–722.  doi:10.1143/PTP.41.705.
  * [10] S. Saito, Theory of Resonating Group Method and Generator Coordinate Method, and Orthogonality Condition Model, Prog. Theor. Phys. Suppl. 62 (1977) 11–89.  doi:10.1143/PTPS.62.11.
  * [11] T. Suhara, Y. Kanada-En’Yo, Cluster structures of excited states in 14C, Phys. Rev. C 82 (4) (2010) 044301.  arXiv:1004.4954, doi:10.1103/PhysRevC.82.044301.
  * [12] Y. Kanada-En’yo, Proton radii of Be, B, and C isotopes, Phys. Rev. C 91 (1) (2015) 014315.  arXiv:1411.0765, doi:10.1103/PhysRevC.91.014315.
  * [13] M. Kimura, T. Suhara, Y. Kanada-En’yo, Antisymmetrized molecular dynamics studies for exotic clustering phenomena in neutron-rich nuclei, Eur. Phys. J. A 52 (2016) 373.  arXiv:1612.09432, doi:10.1140/epja/i2016-16373-9.
  * [14] T. Baba, M. Kimura, Structure and decay pattern of the linear-chain state in 14C, Phys. Rev. C 94 (4) (2016) 044303.  arXiv:1605.05567, doi:10.1103/PhysRevC.94.044303.
  * [15] B. E. Grinyuk, D. V. Piatnytskyi, Structure of 14C and 14O nuclei calculated in the variational approach, Ukr. J. Phys. 61 (8) (2016) 674–680.
  * [16] T. Baba, M. Kimura, Structure and decay pattern of linear-chain states in neutron-rich Carbon isotopes, J. Phys. Conf. Ser. 863 (2017) 012033.  doi:10.1088/1742-6596/863/1/012033.
  * [17] T. Baba, M. Kimura, Three-body decay of linear-chain states in 14C, Phys. Rev. C 95 (6) (2017) 064318.  arXiv:1702.04874, doi:10.1103/PhysRevC.95.064318.
  * [18] B. E. Grinyuk, D. V. Piatnytskyi, Structure of 14N nucleus within a five-cluster model, Ukr. J. Phys. 62 (10) (2017) 835–844.
  * [19] D. R. Tilley, H. R. Weller, C. M. Cheves, Energy levels of light nuclei $A$ = 16-17, Nucl. Phys. A 564 (1993) 1–183.  doi:10.1016/0375-9474(93)90073-7.
  * [20] M. Freer, J. D. Malcolm, N. L. Achouri, N. I. Ashwood, D. W. Bardayan, S. M. Brown, W. N. Catford, K. A. Chipps, J. Cizewski, N. Curtis, K. L. Jones, T. Munoz-Britton, S. D. Pain, N. Soić, C. Wheldon, G. L. Wilson, V. A. Ziman, Resonances in 14C observed in the 4He(10Be ,$\alpha$ )10Be reaction, Phys. Rev. C 90 (5) (2014) 054324.  doi:10.1103/PhysRevC.90.054324.
  * [21] Z. Y. Tian, Y. L. Ye, Z. H. Li, C. J. Lin, Q. T. Li, Y. C. Ge, J. L. Lou, W. Jiang, J. Li, Z. H. Yang, J. Feng, P. J. Li, J. Chen, Q. Liu, H. L. Zang, B. Yang, Y. Zhang, Z. Q. Chen, Y. Liu, X. H. Sun, J. Ma, H. M. Jia, X. X. Xu, L. Yang, N. R. Ma, L. J. Sun, Cluster decay of the high-lying excited states in 14C, Chinese Physics C 40 (11) (2016) 111001.  arXiv:1607.00157, doi:10.1088/1674-1137/40/11/111001.
  * [22] S. Bedoor, A. H. Wuosmaa, M. Albers, M. Alcorta, S. Almaraz-Calderon, B. B. Back, P. F. Bertone, C. M. Deibel, C. R. Hoffman, J. C. Lighthall, S. T. Marley, D. G. Mcneel, R. C. Pardo, K. E. Rehm, J. P. Schiffer, D. V. Shetty, Structure of 14C and 14B from the 14,15C(d,3He)13,14B reactions, Phys. Rev. C 93 (4) (2016) 044323\.  doi:10.1103/PhysRevC.93.044323.
  * [23] A. Fritsch, S. Beceiro-Novo, D. Suzuki, W. Mittig, J. J. Kolata, T. Ahn, D. Bazin, F. D. Becchetti, B. Bucher, Z. Chajecki, X. Fang, M. Febbraro, A. M. Howard, Y. Kanada-En’yo, W. G. Lynch, A. J. Mitchell, M. Ojaruega, A. M. Rogers, A. Shore, T. Suhara, X. D. Tang, R. Torres-Isea, H. Wang, One-dimensionality in atomic nuclei: A candidate for linear-chain $\alpha$ clustering in 14C, Phys. Rev. C 93 (1) (2016) 014321.  doi:10.1103/PhysRevC.93.014321.
  * [24] D. Suzuki, T. Ahn, D. Bazin, F. D. Becchetti, S. Beceiro-Novo, A. Fritsch, J. J. Kolata, W. Mittig, AT-TPC Collaboration, Cluster structure of neutron-rich 10Be and 14C via resonant alpha scattering, Nuovo Cimento C Geophysics Space Physics C 39 (2017) 372.  doi:10.1393/ncc/i2016-16372-0.
  * [25] H. Yamaguchi, D. Kahl, S. Hayakawa, Y. Sakaguchi, K. Abe, T. Nakao, T. Suhara, N. Iwasa, A. Kim, D. H. Kim, S. M. Cha, M. S. Kwag, J. H. Lee, E. J. Lee, K. Y. Chae, Y. Wakabayashi, N. Imai, N. Kitamura, P. Lee, J. Y. Moon, K. B. Lee, C. Akers, H. S. Jung, N. N. Duy, L. H. Khiem, C. S. Lee, Experimental investigation of a linear-chain structure in the nucleus 14C, Phys. Lett. B 766 (2017) 11–16.  arXiv:1610.06296, doi:10.1016/j.physletb.2016.12.050.
  * [26] A. Di Pietro, J. P. Fernández-García, F. Ferrera, P. Figuera, M. Fisichella, M. Lattuada, S. Marletta, C. Marchetta, D. Torresi, M. Alcorta, M. J. G. Borge, T. Davinson, S. Heinitz, A. M. Laird, A. C. Shotter, D. Schumann, N. Soic, O. Tengblad, M. Zadro, Experimental investigation of exotic clustering in 13B and 14C using the resonance scattering method, J. Phys. Conf. Ser. 966 (2018) 012040.  doi:10.1088/1742-6596/966/1/012040.
  * [27] B. E. Grinyuk, I. V. Simenog, Structure of the 6He nucleus in the three-particle model, Phys. Atom. Nucl. 72 (2009) 6–19.  doi:10.1134/S1063778809010025.
  * [28] B. E. Grinyuk, I. V. Simenog, Structural properties of the 10Be and 10C four-cluster nuclei, Phys. Atom. Nucl. 77 (2014) 415–423.  doi:10.1134/S1063778814030090.
  * [29] B. E. Grinyuk, I. V. Simenog, Structure Peculiarities of Three- and Four-Cluster Nuclei 6He, 6Li, and 10Be, 10C, Nucl. Phys. Atom. En. 12 (1) (2011) 7–15.
  * [30] V. G. Neudatchin, V. I. Kukulin, V. L. Korotkikh, V. P. Korennoy, A microscopically substantiated local optical potential for $\alpha$\- $\alpha$\- scattering, Phys. Lett. B 34 (7) (1971) 581–583.  doi:10.1016/0370-2693(71)90142-0.
  * [31] V. I. Kukulin, V. G. Neudatchin, Y. F. Smirnov, An interaction of compound particles and the Pauli principle, Sov. J. Part. Nucl. 10 (1979) 1006\.
  * [32] B. E. Grinyuk, D. V. Piatnytskyi, I. V. Simenog, Structure Characteristics of a 4He Nucleus within the Microscopic Approach, Ukr. J. Phys. 52 (5) (2007) 424–435.
  * [33] V. I. Kukulin, V. M. Krasnopol’sky, A stochastic variational method for Few-Body Syst., J. Phys. G Nucl. Phys. 3 (1977) 795–811.  doi:10.1088/0305-4616/3/6/011.
  * [34] Y. Suzuki, K. Varga, Stochastic Variational Approach to Quantum Mechanical Few-Body Problems, Springer-Verlag, Berlin, Heidelberg, 1998.
  * [35] E. Hiyama, Y. Kino, M. Kamimura, Gaussian expansion method for few-body systems, Prog. Part. Nucl. Phys. 51 (1) (2003) 223–307.  doi:10.1016/S0146-6410(03)90015-9.
  * [36] E. Hiyama, M. Kamimura, Study of various few-body systems using Gaussian expansion method (GEM), Frontiers of Physics 13 (6) (2018) 132106.  arXiv:1809.02619, doi:10.1007/s11467-018-0828-5.
  * [37] R. H. Helm, Inelastic and Elastic Scattering of 187-Mev Electrons from Selected Even-Even Nuclei, Phys Rev. 104 (1956) 1466–1475.  doi:10.1103/PhysRev.104.1466.
  * [38] L. A. Schaller, L. Schellenberg, T. Q. Phan, G. Piller, A. Ruetschi, H. Schneuwly, Nuclear charge radii of the carbon isotopes 12C, 13C and 14C, Nucl. Phys. A 379 (3) (1982) 523–535.  doi:10.1016/0375-9474(82)90012-4.
  * [39] I. Angeli, K. P. Marinova, Table of experimental nuclear ground state charge radii: An update, Atomic Data and Nuclear Data Tables 99 (1) (2013) 69–95.  doi:10.1016/j.adt.2011.12.006.
  * [40] R. F. Frosch, J. S. McCarthy, R. E. Rand, M. R. Yearian, Structure of the He4 Nucleus from Elastic Electron Scattering, Phys. Rev. 160 (4) (1967) 874–879.  doi:10.1103/PhysRev.160.874.
  * [41] P. E. Bosted, L. Andivahis, A. Lung, L. M. Stuart, J. Alster, R. G. Arnold, C. C. Chang, F. S. Dietrich, W. Dodge, R. Gearhart, J. Gomez, K. A. Griffioen, R. S. Hicks, C. E. Hyde-Wright, C. Keppel, S. E. Kuhn, J. Lichtenstadt, R. A. Miskimen, G. A. Peterson, G. G. Petratos, S. E. Rock, S. Rokni, W. K. Sakumoto, M. Spengos, K. Swartz, Z. Szalata, L. H. Tao, Measurements of the electric and magnetic form factors of the proton from Q2=1.75 to 8.83 (GeV/c)2, Phys. Rev. Lett. 68 (26) (1992) 3841–3844.  doi:10.1103/PhysRevLett.68.3841.