# Mass spectra of hidden heavy-flavor tetraquarks with two and four heavy
quarks

Ting-Qi Yan1<EMAIL_ADDRESS>Wen-Xuan Zhang1<EMAIL_ADDRESS>Duojie Jia1,2 111Corresponding author<EMAIL_ADDRESS>1Institute of
Theoretical Physics, College of Physics and Electronic Engineering, Northwest
Normal University, Lanzhou 730070, China
2Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou, 730000,
China

###### Abstract

Inspired by the observation of the $X(6900)$ by LHCb and the $X(6600)$ (with
mass $6552\pm 10$ $\pm 12$ MeV) recently by CMS and ATLAS experiments of the
LHC in the di-$J/\Psi$ invariant mass spectrum, we systemically study masses
of all ground-state configurations of the hidden heavy-flavor tetraquarks
$q_{1}Q_{2}\bar{q}_{3}\bar{Q}_{4}$ and $Q_{1}Q_{2}\bar{Q}_{3}\bar{Q}_{4}$
($Q=c,b$;$q=u,d,s$) contaning two and four heavy quarks in the MIT bag model
with chromomagnetic interaction and enhanced binding energy. Considering
color-spin mixing due to chromomagnetic interaction, our mass computation
indicates that the observed $X(6600)$ is likely to be the $0^{++}$ ground
states of hidden-charm tetraquark $cc\bar{c}\bar{c}$ with computed masses
$6572$ MeV, which has a $0^{++}$ color partner around $6469$ MeV. The fully
bottom system of tetraquark $bb\bar{b}\bar{b}$ has masses of 19685 MeV and
19717 MeV for the the $0^{++}$ ground states. Further computation is given to
the tetraquark systems $sc\bar{s}\bar{c}$, $sb\bar{s}\bar{b}$,
$cb\bar{c}\bar{b}$, $nc\bar{n}\bar{c}$ and $nb\bar{n}\bar{b}$, suggesting that
the $Z_{c}(4200)$ is the tetraquark $nc\bar{n}\bar{c}$ with $J^{PC}=1^{+-}$.
All of these tetraquarks are above their lowest thresholds of two mesons and
unstable against the strong decays.

PACS number(s):12.39Jh, 12.40.Yx, 12.40.Nn

Key Words: Heavy pentaquark, Spectroscopy, Quantum number

##  I Introduction

All known strongly interacting particles (mesons and baryons) could be
classified as bound states made of a quark-antiquark pair or three quarks for
a long time based on the conventional scheme of the quark model by Gell-
MannGell-Mann:1964ewy  and ZweigZweig:1964ruk . Meanwhile, they also suggested
possible existence of the hadron states of multiquarks like tetraquarks (with
quark configuration $q^{2}\bar{q}^{2}$) and pentaquarks ($q^{4}\bar{q}$). In
the 1970s, multiquark states (the exotic light mesons like the $a_{0}$ and
$f_{0}$) are calculated by Jaffe based on the dynamical framework of the MIT
bag model Jaffe:1976ig ; Jaffe:1976ih . Despite that multiquarks are
considered to be exotic in the sense that they go beyond the conventional
scheme of quark model, they are, in principle, allowed by the quantum
chromodynamics (QCD), the theory of the strong force that binds quarks into
hadrons.

Since observation of the first exotic hadron $X(3872)$ Belle:2003nnu  in 2003
by the Belle, many (more than 20) tetraquark candidates have been observed
among charmonium-like or bottomonium-like $XYZ$ states, which include the
charmonium-like states the $Z_{c}(3900)$ BESIII:2013ris , the $Z_{c}(4200)$
Belle:2014nuw , the $Z_{c}(4430)$ Belle:2007hrb ; Belle:2013shl ; LHCb:2014zfx
; LHCb:2015sqg . Some of the observed $XYZ$ states, like the charged state
$Z_{c}(3900)$ BESIII:2013ris , are undoubtedly exotic. In 2020, a candidate of
fully charm tetraquark, the $X(6900)$, has been observed by LHCb in the
di-$J/\Psi$ invariant mass spectrum around the mass of $6905$ MeV, which is
later confirmed by CMS and ATLAS of the LHC at CERN Bouhova-Thacker:2022vnt ;
Zhang:2022toq ; LHCb:2020bwg . Meanwhile in the same di-$J/\Psi$ invariant
mass spectrum, a new structure, the $X(6600)$, are also found by CMS with mass
of $6552\pm 10$ $\pm 12$ MeV, which is very likely to be the fully charm
tetraquark.

The purpose of this work is to use the MIT bag model with enhanced binding
energy to systemically study the ground-state masses of the hidden heavy-
flavor tetraquarks containing two or four heavy quarks. Based on color-spin
wavefunctions constructed for the hidden heavy-flavor tetraquarks, we solve
the bag model and diagonalize the chromomagnetic interaction (CMI) to take
into account the possible color-spin mixing of the states with same quantum
numbers. We find that the computed masses of the fully charmed tetraquark
$cc\bar{c}\bar{c}$ is in a good agreement with the mass measurement by the CMS
experimentZhang:2022toq . Further mass computation is performed for hidden
heavy-flavor systems of the tetraquarks $bb\bar{b}\bar{b}$,
$cb\bar{c}\bar{b}$, $sc\bar{s}\bar{c}$, $sb\bar{s}\bar{b}$,
$nc\bar{n}\bar{c}$, $nb\bar{n}\bar{b}$, with a suggestion that the particle
$Z_{c}(4200)$ reported by Belle:2014nuw  is likely to be the hidden-charm
tetraquark $nc\bar{n}\bar{c}$ with $J^{PC}=1^{+-}$.

In Section 2, we present the allowed wavefunctions of hidden heavy-flavor
tetraquarks with two or four heavy quarks. In Section 3, We describe the
framework of MIT bag model to be used in this work. The mass matrix evaluation
for the CMI and its diagonalization are detailed in section 4. The masses of
the hidden heavy-flavor tetraquarks are computed numerically for the systems
($cc\bar{c}\bar{c}$, $bb\bar{b}\bar{b}$), $cb\bar{c}\bar{b}$,
($sc\bar{s}\bar{c}$, $sb\bar{s}\bar{b}$) and ($nc\bar{n}\bar{c}$,
$nb\bar{n}\bar{b}$) in Section 5. We end with conclusions and remarks in
Section 6.

##  II Wavefunctions of hidden-flavor tetraquarks

We consider hidden heavy-flavor tetraquarks containing two or four heavy
quarks($q_{1}Q_{2}\bar{q}_{3}\bar{Q}_{4}$ and
$Q_{1}Q_{2}\bar{Q}_{3}\bar{Q}_{4}$, $Q=c,b$, $q=u,d,s$), which include seven
flavor combinations of four quark systems: $cc\bar{c}\bar{c}$,
$bb\bar{b}\bar{b}$, $sc\bar{s}\bar{c}$, $sb\bar{s}\bar{b}$,
$cb\bar{c}\bar{b}$, $nc\bar{n}\bar{c}$, $nb\bar{n}\bar{b}$, with $n=u$, $d$.
In this section, we describe the wavefunctions of the hidden-flavor
tetraquarks in the flavor and the color-spin space.

In the flavor space, we utilized $\delta_{12}^{S}$ if $q_{1}q_{2}$ is
symmetric and $\delta_{12}^{A}\equiv 1-\delta_{12}^{A}$ if $q_{1}q_{2}$ is
antisymmetric to restrict the wavefunction. If the wavefunction has no flavor
symmetry (beyond the isospin symmetry $SU(2)_{I}$) under the exchange of
$q_{1}$ and $q_{2}$, then $\delta_{12}^{S}=\delta_{12}^{A}=1$.

In color space, the hidden heavy-flavor tetraquark
$q_{1}q_{2}\bar{q}_{3}\bar{q}_{4}$ can be in two color states:
$6_{c}\otimes\bar{6}_{c}$ and $\bar{3}_{c}\otimes 3_{c}$, with the respective
wave functions (superscript stands for color representation),

$\phi_{1}^{T}=\left|{\left(q_{1}q_{2}\right)}^{6}{\left(\bar{q}_{3}\bar{q}_{4}\right)}^{\bar{6}}\right\rangle,\quad\phi_{2}^{T}=\left|{\left(q_{1}q_{2}\right)}^{\bar{3}}{\left(\bar{q}_{3}\bar{q}_{4}\right)}^{3}\right\rangle.$
(1)

With the help of the color $SU(3)_{c}$ symmetry, one can write the two
configurations $\phi_{1,2}^{T}$ here in terms of the fundamental
representations, i.e., of the color bases $c_{n}=|r\rangle$, $|b\rangle$ and
$|g\rangle$ of the $SU(3)_{c}$ group (see Appendix A).

In the spin space, there are six states of a tetraquark state allowed
(Appendix A), with the wavefunctions (subscript stands for spin),

$\displaystyle\chi_{1}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}\right\rangle}_{2},\quad\chi_{2}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}\right\rangle}_{1},$
(2)
$\displaystyle\chi_{3}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}\right\rangle}_{0},\quad\chi_{4}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{0}\right\rangle}_{1},$
$\displaystyle\chi_{5}^{T}={\left|{\left(q_{1}q_{2}\right)}_{0}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}\right\rangle}_{1},\quad\chi_{6}^{T}={\left|{\left(q_{1}q_{2}\right)}_{0}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{0}\right\rangle}_{0}.$

Based on the Pauli’s principle, one can construct twelve color-spin
wavefunctions for the lowest S-wave (in coordinate space) tetraquarks:

$\displaystyle\phi_{1}^{T}\chi_{1}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{6}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{\bar{6}}\right\rangle}_{2}\delta_{12}^{A}\delta_{34}^{A},$
(3)
$\displaystyle\phi_{2}^{T}\chi_{1}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{\bar{3}}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{3}\right\rangle}_{2}\delta_{12}^{S}\delta_{34}^{S},$
$\displaystyle\phi_{1}^{T}\chi_{2}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{6}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{\bar{6}}\right\rangle}_{1}\delta_{12}^{A}\delta_{34}^{A},$
$\displaystyle\phi_{2}^{T}\chi_{2}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{\bar{3}}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{3}\right\rangle}_{1}\delta_{12}^{S}\delta_{34}^{S},$
$\displaystyle\phi_{1}^{T}\chi_{3}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{6}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{\bar{6}}\right\rangle}_{0}\delta_{12}^{A}\delta_{34}^{A},$
$\displaystyle\phi_{2}^{T}\chi_{3}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{\bar{3}}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{3}\right\rangle}_{0}\delta_{12}^{S}\delta_{34}^{S},$
$\displaystyle\phi_{1}^{T}\chi_{4}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{6}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{0}^{\bar{6}}\right\rangle}_{1}\delta_{12}^{A}\delta_{34}^{S},$
$\displaystyle\phi_{2}^{T}\chi_{4}^{T}={\left|{\left(q_{1}q_{2}\right)}_{1}^{\bar{3}}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{0}^{3}\right\rangle}_{1}\delta_{12}^{S}\delta_{34}^{A},$
$\displaystyle\phi_{1}^{T}\chi_{5}^{T}={\left|{\left(q_{1}q_{2}\right)}_{0}^{6}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{\bar{6}}\right\rangle}_{1}\delta_{12}^{S}\delta_{34}^{A},$
$\displaystyle\phi_{2}^{T}\chi_{5}^{T}={\left|{\left(q_{1}q_{2}\right)}_{0}^{\bar{3}}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{1}^{3}\right\rangle}_{1}\delta_{12}^{A}\delta_{34}^{S},$
$\displaystyle\phi_{1}^{T}\chi_{6}^{T}={\left|{\left(q_{1}q_{2}\right)}_{0}^{6}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{0}^{\bar{6}}\right\rangle}_{0}\delta_{12}^{S}\delta_{34}^{S},$
$\displaystyle\phi_{2}^{T}\chi_{6}^{T}={\left|{\left(q_{1}q_{2}\right)}_{0}^{\bar{3}}{\left(\bar{q}_{3}\bar{q}_{4}\right)}_{0}^{3}\right\rangle}_{0}\delta_{12}^{A}\delta_{34}^{A}.$

We choose these wavefunctions to be the bases (the first approximation) of the
tetraquark eigenstates for which the chromomagnetic interaction (CMI) are
ignored. We are going to employ these bases to take into account the
chromomagnetic mixing due to the CMI. For example, for the $J^{PC}=0^{++}$
state of the $cc\bar{c}\bar{c}$ tetraquark, one can write two bases of the
wavefunctions $\phi_{2}^{T}\chi_{3}^{T}$, $\phi_{1}^{T}\chi_{6}^{T}$ in Eq.
(3) as a zero-order approximation, which satisfy the required symmetry in the
color-spin space and can lead to mixing of the color-spin states when the CMI
added.

Table 1: Allowed state mixing of the hidden heavy-flavor tetraquarks due to chromomagnetic interaction. State | $J^{PC}$ | Allowed states for mixing
---|---|---
$QQ\bar{Q}\bar{Q}$ | $0^{++}$ | $\left(\phi_{2}^{T}\chi_{3}^{T},\phi_{1}^{T}\chi_{6}^{T}\right)$
| $1^{+-}$ | $\left(\phi_{2}^{T}\chi_{2}^{T}\right)$
| $2^{++}$ | $\left(\phi_{2}^{T}\chi_{1}^{T}\right)$
$Q{Q}^{\prime}\bar{Q}\bar{Q}^{\prime}$,$qQ\bar{q}\bar{Q}$ | $0^{++}$ | $\left(\phi_{2}^{T}\chi_{3}^{T},\phi_{2}^{T}\chi_{6}^{T},\phi_{1}^{T}\chi_{3}^{T},\phi_{1}^{T}\chi_{6}^{T}\right)$
| $1^{++}$ | ($\frac{1}{\sqrt{2}}\left(\phi_{2}^{T}\chi_{4}^{T}+\phi_{2}^{T}\chi_{5}^{T}\right),\frac{1}{\sqrt{2}}\left(\phi_{1}^{T}\chi_{4}^{T}+\phi_{1}^{T}\chi_{5}^{T}\right)$)
| $1^{+-}$ | ($\phi_{2}^{T}\chi_{2}^{T},\frac{1}{\sqrt{2}}\left(\phi_{2}^{T}\chi_{4}^{T}-\phi_{2}^{T}\chi_{5}^{T}\right),\phi_{1}^{T}\chi_{2}^{T},\frac{1}{\sqrt{2}}\left(\phi_{1}^{T}\chi_{4}^{T}-\phi_{1}^{T}\chi_{5}^{T}\right)$)
| $2^{++}$ | $\left(\phi_{2}^{T}\chi_{1}^{T},\phi_{1}^{T}\chi_{1}^{T}\right)$

With respect to given flavor compositions of the tetraquarks, one can write
the allowed color-spin states that may mix due to the CMI for each choice of
the quantum number $J^{PC}$ in Table 1, where $Q^{\prime}$ denoting heavy
quark differing with $Q$. Note that for the flavor composition
$QQ\bar{Q}\bar{Q}$ with quantum numbers $J^{PC}=1^{+-}$ and $2^{++}$, there is
only one color-spin state for each of them, that is, the
$\phi_{2}^{T}\chi_{2}^{T}$ associated with $1^{+-}$ and
$\phi_{2}^{T}\chi_{1}^{T}$ associated with $2^{++}$, for which not mixing
occurs in reality.

##  III The MIT bag model

We use the MIT bag model which includes enhanced binding energy and the CMI in
the interaction correction $\Delta M$. The mass formula for the MIT bag model
isZhang:2021yul

$M\left(R\right)=\sum_{i}\omega_{i}+\frac{4}{3}\pi
R^{3}B-\frac{Z_{0}}{R}+\Delta M,$ (4)
$\omega_{i}=\left(m_{i}^{2}+\frac{x_{i}^{2}}{R^{2}}\right)^{1/2},$ (5)

with the first term describes (relativistic) kinetic motion of each quark $i$
in tetraquark, the second is the volume energy of bag with bag constant $B$,
the third is the zero-point-energy with coefficient $Z_{0}$ and $R$ the bag
radius to be determined variationally. In Eq. (5), the dimensionless
parameters $x_{i}=x_{i}(mR)$ are related to $R$ through an transcendental
equation

$\tan
x_{i}=\frac{x_{i}}{1-m_{i}R-\left(m_{i}^{2}R^{2}+x_{i}^{2}\right)^{1/2}}.$ (6)

In Eq. (4), we denote the sum of the first three terms to be $M(T)$. The
interaction correction $\Delta M$ includes the enhanced binding energy $M_{B}$
among the quarks in tetraquark and the mass splitting $M_{CMI}$ corresponding
to the CMI:

$\Delta M=M_{B}+M_{CMI}=\sum_{i<j}B_{ij}+\langle H_{CMI}\rangle,$ (7)

where $B_{ij}$ stands for the binding energyKarliner:2014gca ;
Karliner:2017elp  between quarks $i$ and $j$, described below at the end of
this section, and the chromomagnetic interaction $H_{CMI}$ is given by

$H_{CMI}=-\sum_{i<j}\left(\mathbf{\lambda}_{\mathbf{i}}\cdot\mathbf{\lambda}_{j}\right)\left(\mathbf{\sigma}_{\mathbf{i}}\cdot\mathbf{\sigma}_{j}\right)C_{ij}.$
(8)

where $\mathbf{\lambda}\bm{{}_{i}}$ and $\mathbf{\sigma}\bm{{}_{i}}$ are the
Gell-Mann and Pauli matrices of the quark $i$, respectively, and $C_{ij}$ the
CMI coupling parameters, given byDeGrand:1975cf

$C_{ij}=3\frac{\alpha_{s}\left(R\right)}{R^{3}}\bar{\mu}_{i}\bar{\mu}_{j}I_{ij},$
(9)

with $\alpha_{s}(R)$ is the running coupling given in Ref. Zhang:2021yul ,
$\bar{\mu}_{i}$ the reduced magnetic moment of quark $i$,

$\alpha_{s}(R)=\frac{0.296}{ln\left[1+{\left(0.281R\right)}^{-1}\right]}.$
(10)
$\bar{\mu}_{i}=\frac{R}{6}\frac{4{\omega_{i}R}+2\lambda_{i}-3}{2{\omega_{i}R}\left({\omega_{i}R}-1\right)+\lambda_{i}},$
(11)

and

$I_{ij}=1+2\int_{0}^{R}\frac{dr}{r^{4}}\bar{\mu}_{i}\bar{\mu}_{j}=1+F\left(x_{i},x_{j}\right).$
(12)

where $\lambda_{i}\equiv{m_{i}R}$. The function $F\left(x_{i},x_{j}\right)$ is
given by

$\displaystyle
F\left(x_{i},x_{j}\right)=\left(x_{i}\sin^{2}x_{i}-\frac{3}{2}y_{i}\right)^{-1}\left(x_{j}\sin^{2}x_{j}-\frac{3}{2}y_{j}\right)^{-1}$
(13)
$\displaystyle\left\\{-\frac{3}{2}y_{i}y_{j}-2x_{i}x_{j}sin^{2}x_{i}sin^{2}x_{j}+\frac{1}{2}x_{i}x_{j}\left[2x_{i}Si\left(2x_{i}\right)\right.\right.$
$\displaystyle\left.\left.+2x_{j}Si\left(2x_{j}\right)-\left(x_{i}+x_{j}\right)Si\left(2\left(x_{i}+x_{j}\right)\right)\right.\right.$
$\displaystyle\left.\left.-\left(x_{i}-x_{j}\right)Si\left(2\left(x_{i}-x_{j}\right)\right)\right]\right\\}$

where $y_{i}=x_{i}-\cos\left(x_{i}\right)\sin\left(x_{i}\right)$, $x_{i}$ is
the solution of Eq. (6), and

$Si(x)=\int\limits_{0}^{x}\frac{\sin(t)}{t}dt.$ (14)

Note that the functional of the running coupling $\alpha_{s}(R)$ in Eq. (9)
and other parameters (the quark mass $m_{i}$, zero-point energy coefficient
$Z_{0}$, bag constant $B$) are evaluated in Ref. Zhang:2021yul  via mapping
the model mass prediction to the ground-state masses of the observed mesons
and baryons. The obtained values for these model parameters areZhang:2021yul

$\begin{Bmatrix}m_{n}=0\,\text{GeV,}&m_{s}=0.279\,\text{GeV,}\\\
m_{c}=1.641\,\text{GeV,}&m_{b}=5.093\,\text{GeV,}\\\
Z_{0}=1.84,&B^{1/4}=0.145\,\text{GeV.}\end{Bmatrix}$ (15)

We will use these parameters to analyze the heavy tetra- quarks in this work,
with the bag radius $R$ determined variationally via the MIT bag model.

The binding energy $M_{B}$ in Eq. (7) measures the short-range chromoelectric
interaction between quarks and/or antiquarks. For the massive quarks of $i$
and $j$, this energy, which scales like $-\alpha_{s}(r_{ij})/r_{ij}$, becomes
sizable when both quarks($i$ and $j$) are massive, moving nonrelativistically.
We treat this energy as the sum of the pair binding energies,
$B_{QQ^{\prime}}(B_{Qs})$, between heavy quarks ($Q$ and $Q^{\prime}$) and
between heavy quarks $Q$ and the strange quarks $s$Karliner:2014gca ;
Karliner:2017elp . This leads to five binding energies $B_{cs}$, $B_{cc}$,
$B_{bs}$, $B_{bb}$, and $B_{bc}$ for any quark pair in the color configuration
$\bar{3}_{c}$, which are extractable from heavy mesons and can be scaled to
other color configurations.

Assuming two quarks $QQ^{\prime}$ to be in the color anti-triplet
$\bar{3}_{c}$ inside baryon, the binding energy $B_{QQ^{\prime}}\equiv
B_{QQ^{\prime}}[\bar{3}_{c}]$ are extracted in the MIT bag modelZhang:2021yul
(appendix A) for the combination of $QQ^{\prime}=cc$, $bb$, $bc$, $bs\ $and
$cs$, so that a unified parameter setup was established for the ground states
of meson, baryons and heavy hadrons (including doubly baryon and tetraquar-
ks). The results areZhang:2021yul

$\begin{Bmatrix}B_{cs}=-0.025\,\text{GeV,}&B_{cc}=-0.077\,\text{GeV,}\\\
B_{bs}=-0.032\,\text{GeV,}&B_{bb}=-0.128\,\text{GeV,}\\\
B_{bc}=-0.101\,\text{GeV.}&\end{Bmatrix}$ (16)

##  IV Color and spin factors for tetraquarks

To determine the mass splitting $M_{CMI}=\langle H_{CMI}\rangle$ via the CMI
Hamiltonian $H_{CMI}$ in Eq. (8), one has to evaluate the chromomagnetic
matrices $H_{CMI}$ of the tetraquarks $T$ for a given quantum number $J^{PC}$.
For this, one can firstly work out the color factors
${\left\langle\bm{\lambda}_{i}\cdot\bm{\lambda}_{j}\right\rangle}$ and spin
factors ${\left\langle\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}\right\rangle}$ as
matrices over the color and spin bases, respectively,the allowed states of
tetraquarks with given $J^{PC}$ in Table 1. In this section, we present the
color and spin factors as a matrix elements in the color and spin space, and
give an unified expressions for binding energy
$M_{B}=\sum_{i<j}F_{c}B_{ij}(\bar{3}_{c})$ for the both color :

Color factor in the color states $|n\rangle$ and $|m\rangle$:

${\left\langle\bm{\lambda}_{i}\cdot\bm{\lambda}_{j}\right\rangle}_{nm}=\sum_{\alpha=1}^{8}Tr\left(c_{in}^{\dagger}\lambda^{\alpha}c_{im}\right)Tr\left(c_{jn}^{\dagger}\lambda^{\alpha}c_{jm}\right),$
(17)

and spin factor in the spin states $|x\rangle$ and $|y\rangle$:

${\left\langle\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}\right\rangle}_{xy}=\sum_{\alpha=1}^{3}Tr\left(\chi_{ix}^{\dagger}\sigma^{\alpha}\chi_{iy}\right)Tr\left(\chi_{jx}^{\dagger}\sigma^{\alpha}\chi_{jy}\right),$
(18)

where $c_{in}$ stands for color basis (three colors $r$, $g$, and $b$) of a
given quark $i$, and $\chi_{ix}$ represents its spin basis (with two spin
components of $\uparrow$ and $\downarrow$).

In color-spin wavefunction of the tetraquark $T$, one can compute explicitly
the expectation values of $H_{CMI}$,

${\left\langle
T|H_{CMI}|T\right\rangle=}-\sum_{i<j}{\left\langle\bm{\lambda}_{i}\cdot\bm{\lambda}_{j}\right\rangle}_{TT}{\left\langle\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}\right\rangle}_{TT}C_{ij}\text{,}$
(19)

to obtain the color and spin factor, writing the mass formula for $M_{CMI}$ in
terms of the CMI couplings $C_{ij}$, which are given further by Eq. (9) in the
MIT bag model. Here the state of $T$ are the mixed states listed in Table 1,
with the mixed weight $w=(w_{1},w_{2},\cdots,w_{f})$ solved (as eigenvector
during the CMI diagonalization) numerically in Table 2,5-7 in the section 5.

Given the two formula (17) and (18), one can compute the color factors
${\left\langle\phi_{1}^{T},\phi_{2}^{T}|\bm{\lambda}_{i}\cdot\bm{\lambda}_{j}|\phi_{1}^{T},\phi_{2}^{T}\right\rangle}$
as $2$ by $2$ matrix in the color subspace of
$\left(\phi_{1}^{T},\phi_{2}^{T}\right)$, via applying Eqs. (30) and (31) in
appendix A. The result are obtained to be

$\displaystyle\left\langle\bm{\lambda_{1}}\cdot\bm{\lambda_{2}}\right\rangle=\left\langle\bm{\lambda_{3}}\cdot\bm{\lambda_{4}}\right\rangle=\begin{bmatrix}\frac{4}{3}&0\\\
0&-\frac{8}{3}\end{bmatrix},$ (20)
$\displaystyle\left\langle\bm{\lambda_{1}}\cdot\bm{\lambda_{3}}\right\rangle=\left\langle\bm{\lambda_{2}}\cdot\bm{\lambda_{4}}\right\rangle=\begin{bmatrix}-\frac{10}{3}&2\sqrt{2}\\\
2\sqrt{2}&-\frac{4}{3}\end{bmatrix},$
$\displaystyle\left\langle\bm{\lambda_{1}}\cdot\bm{\lambda_{4}}\right\rangle=\left\langle\bm{\lambda_{2}}\cdot\bm{\lambda_{3}}\right\rangle=\begin{bmatrix}-\frac{10}{3}&-2\sqrt{2}\\\
-2\sqrt{2}&-\frac{4}{3}\end{bmatrix}.$

From the above matrices, we see that the color conifgurations $\phi_{1}^{T}\
$and $\phi_{2}^{T}$ may mix for a tetraquark state $T$ due to the
chromomagnetic interaction.

We further consider the binding energy $M_{B}$ based on Eq. (16), which
corresponds to the binding energy $B_{ij}\equiv B_{ij}[\bar{3}_{c}]$ in
baryons with the quark pair ($i,j$) in $\bar{3}_{c}$. Let us then consider the
binding energy $M_{B}$ for a given color configurations of the tetraquark
$T=(q_{1}q_{2})^{R}(\bar{q}_{3}\bar{q}_{4})^{\bar{R}}$(with representation
$R=$ $6_{c}$ and $\bar{3}_{c}$). First of all, one can scale the pair binding
energy $B_{ij}\equiv B_{ij}[\bar{3}_{c}]$ of the pair in baryon to
$F_{c}[R]B_{ij}[\bar{3}_{c}]$ of the pair in tetraquark $T$, where $F_{c}[R]$
is the ratio of the color factor in Eq. (20) to the color factor
${\left\langle\bm{\lambda}_{i}\cdot\bm{\lambda}_{j}\right\rangle}_{B}=-8/3$
for baryon with each of quark pair ($i,j$) in $\bar{3}_{c}$. At last, applying
to all quark pair ($i,j$) of the tetraquark $T$ with configurations
$\phi_{1}^{T}\ $and $\phi_{2}^{T}$, one can obtain the pair binding energies
$F_{c}[R]B_{ij}$, whose sums are,

$M_{B}(\phi_{1}^{T})=-\frac{1}{2}B_{12}-\frac{1}{2}B_{34}+\frac{5}{4}B_{13}+\frac{5}{4}B_{14}+\frac{5}{4}B_{23}+\frac{5}{4}B_{24},$
(21)
$M_{B}(\phi_{2}^{T})=B_{12}+B_{34}+\frac{1}{2}B_{13}+\frac{1}{2}B_{14}+\frac{1}{2}B_{23}+\frac{1}{2}B_{24},$
(22)

for the tetraquark $T$, respectively, where $B_{ij}$ is the binding energy
with ($i,j$) in $\bar{3}_{c}$.

For color sextets of the pair ($1,2$) and ($3,4$), for instance, the binding
energy is $-B_{12}/2$ and $-B_{34}/2$, respectively, with
$F_{c}=(4/3)/(-8/3)=-1/2$. For any of representation of the quark $i$ and
antiquark $j$, the binding energies in $T$ are either $-5B_{ij}/4$ or
$B_{ij}/2$. We note that $B_{ij}$ vanishes if both of quark $i$ and $j$ are
light quarks or one of them is non-strange light quark ($B_{nQ}=0$,
$B_{nn}=0$, $B_{n\bar{n}}=0$, $B_{s\bar{s}}=0$) since the short range
interactions between ($i,j$) quarks are small and thereby ignorable averagely
for quark pair $(i,j)=(n,n)$ or $(i,j)=(s,s)$, due to their relativistic
motion.

We come to consider the spin factors, which is given by
${\left\langle\chi^{T}|\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}|\chi^{T}\right\rangle}$.
In the subspace spanned by $\\{\chi_{1\sim 6}^{T}\\}$ in Eq. (32), the direct
computation yields the following matrices,

$\left\langle\bm{\sigma_{1}}\cdot\bm{\sigma_{2}}\right\rangle_{\chi_{1}^{T}}=\begin{bmatrix}1&0&0&0&0&0\\\
0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&1&0&0\\\ 0&0&0&0&-3&0\\\
0&0&0&0&0&-3\end{bmatrix},$ (23)
$\left\langle\bm{\sigma_{1}}\cdot\bm{\sigma_{3}}\right\rangle_{\chi_{2}^{T}}=\begin{bmatrix}1&0&0&0&0&0\\\
0&-1&0&\sqrt{2}&-\sqrt{2}&0\\\ 0&0&-2&0&0&-\sqrt{3}\\\ 0&\sqrt{2}&0&0&1&0\\\
0&-\sqrt{2}&0&1&0&0\\\ 0&0&-\sqrt{3}&0&0&0\end{bmatrix},$ (24)
$\left\langle\bm{\sigma_{1}}\cdot\bm{\sigma_{4}}\right\rangle_{\chi_{3}^{T}}=\begin{bmatrix}1&0&0&0&0&0\\\
0&-1&0&-\sqrt{2}&-\sqrt{2}&0\\\ 0&0&-2&0&0&\sqrt{3}\\\ 0&-\sqrt{2}&0&0&-1&0\\\
0&-\sqrt{2}&0&-1&0&0\\\ 0&0&\sqrt{3}&0&0&0\end{bmatrix},$ (25)
$\left\langle\bm{\sigma_{2}}\cdot\bm{\sigma_{3}}\right\rangle_{\chi_{4}^{T}}=\begin{bmatrix}1&0&0&0&0&0\\\
0&-1&0&\sqrt{2}&\sqrt{2}&0\\\ 0&0&-2&0&0&\sqrt{3}\\\ 0&\sqrt{2}&0&0&-1&0\\\
0&\sqrt{2}&0&-1&0&0\\\ 0&0&\sqrt{3}&0&0&0\end{bmatrix},$ (26)
$\left\langle\bm{\sigma_{2}}\cdot\bm{\sigma_{4}}\right\rangle_{\chi_{5}^{T}}=\begin{bmatrix}1&0&0&0&0&0\\\
0&-1&0&-\sqrt{2}&\sqrt{2}&0\\\ 0&0&-2&0&0&-\sqrt{3}\\\ 0&-\sqrt{2}&0&0&1&0\\\
0&\sqrt{2}&0&1&0&0\\\ 0&0&-\sqrt{3}&0&0&0\end{bmatrix},$ (27)
$\left\langle\bm{\sigma_{3}}\cdot\bm{\sigma_{4}}\right\rangle_{\chi_{6}^{T}}=\begin{bmatrix}1&0&0&0&0&0\\\
0&1&0&0&0&0\\\ 0&0&1&0&0&0\\\ 0&0&0&-3&0&0\\\ 0&0&0&0&1&0\\\
0&0&0&0&0&-3\end{bmatrix}.$ (28)

Combing the spin factors in Eqs. (23)-(28) with Eqs. (20), we are the position
to use Eqs. (19), Eq. (17) and Eq. (18) to compute the mass splitting
$M_{CMI}$ duo to chromomagnetic interaction. Using Eqs. (21) and (22), one can
compute the mass sum $\Delta M=M_{B}+M_{CMI\text{ }}$in Eq. 7 and further
obtain, via adding mass of the bag $M_{bag}=\sum_{i}\omega_{i}+(4/3)\pi
R^{3}B-Z_{0}/R$, a complete mass formula for the hidden heavy-flavor
tetraquark systems $T$ addressed in this work,

$M(T)=M_{bag}+M_{B}+M_{CMI\text{ }}(C_{ij}),$ (29)

in which $M_{CMI\text{ }}(C_{ij})$ are linear functions of the CMI couplings
$C_{ij}$, with the linear coefficients given by the color and spin factors
shown in this section.

##  V Masses of hidden heavy-flavor tetraquarks

Given the input parameters in Eqs. (15), one can numerically solve Eq. (4)
variationally, with the mass splitting $M_{CMI}$ and the CMI couplings
$C_{ij}$ given by Eqs. (9), (10), (11) and (12), to obtain bag radius $R$ and
numerically give the masses $M(T)$ of the hidden heavy-flavor tetraquarks $T$.
Meanwhile, we show the numerical corresponding results for the bag radius
$R_{0}$, the mixing weights (eigenvectors of the CMI matrix $H_{CMI\text{
}}$), the tetraquark masses $M(T)$ and thresholds of two mesons as a final
states in the Tables 2, 5-7. In the following, we present the results and
discussions with respect to the tetraquark systems addressed below in order.

###  V.1 Fully heavy tetraquark systems

In the case of fully charmed systems of the tetraquarks $cc\bar{c}\bar{c}$, we
show the numerical results for $R_{0}$, the state-mixing weights (eigenvectors
of $H_{CMI\text{ }}$), the tetraquark masses $M(T)$ and thresholds of two
mesons final states in the Table 2, with the later two plotted in Fig. 1. We
see that for $J^{PC}=0^{++}$ there are two states of the tetraquarks
$cc\bar{c}\bar{c}$ with the masses of $6572\,$MeV and $6469\,$MeV, splitted by
$103\,$MeV. The tetraquark ($cc\bar{c}\bar{c}$) states with $J^{PC}=1^{+-}$
and $J^{PC}=2^{++}$ have the masses within a similar mass region, as shown in
Fig 1. We find that all these $cc\bar{c}\bar{c}$ states relatively far above
their two mesons thresholds shown explicitly. For instance, the $0^{++}$ state
are all above the thresholds of the $J/\psi J/\psi$ and $\eta_{c}\eta_{c}$,
about $275-605\,$MeV, indicating that they are not stable against strong
decays through quark rearrangement to the final state of $J/\psi J/\psi$ as
well as $\eta_{c}\eta_{c}$. For the $1^{+-}$ state, there is one state, and
its mass is above the thresholds of the two mesons $\eta_{c}J/\psi$ and
$J/\psi J/\psi$ about $325-440\,$MeV, unstable against the strong decay to the
later. In the case of the $2^{++}$ state, there is one state with the mass
above the threshold($J/\psi$$J/\psi$) about $350\,$MeV, also strongly
unstable. We also compare our calculations with other works cited and list the
results in Table 3.

Figure 1: The computed masses (MeV the solid lines) of the $cc\bar{c}\bar{c}$ tetraquark system in their ground-states, as well as two meson thresholds (MeV the dotted lines).  Figure 2: The computed masses (MeV the solid lines) of the $bb\bar{b}\bar{b}$ tetraquark system in their ground-states, as well as two meson thresholds (MeV the dotted lines).  Table 2: The numerical results for the bag radius $R_{0}$, the state-mixing weights (eigenvectors of $H_{CMI}$), the tetraquark masses $M(T)$ and thresholds of two mesons final states for the hidden heavy-flavor tetraquarks($cc\bar{c}\bar{c}$, $bb\bar{b}\bar{b}$) State | $J^{PC}$ | Eigenvector |  $R_{0}$(GeV-1) |  $M(T)$(MeV) | Threshold (MeV)
---|---|---|---|---|---
$cc\bar{c}\bar{c}$ | $0^{++}$ | $(0.58,0.81)$ | $4.74$ | $6572$ |  $J/\psi J/\psi=6194$;$\eta_{c}\eta_{c}=5967$
|  | $(-0.81,0.58)$ | $4.44$ | $6469$ |
| $1^{+-}$ | $1.00$ | $4.59$ | $6519$ |  $\eta_{c}J/\psi=6080$;$J/\psi J/\psi=6194$
| $2^{++}$ | $1.00$ | $4.66$ | $6545$ | $J/\psi J/\psi=6194$
$bb\bar{b}\bar{b}$ | $0^{++}$ | $(0.58,0.81)$ | $3.15$ | $19717$ |  $\Upsilon\Upsilon=18921$;$\eta_{b}\eta_{b}=18798$
|  | $(-0.81,0.58)$ | $2.99$ | $19685$ |
| $1^{+-}$ | $1.00$ | $3.07$ | $19700$ |  $\eta_{b}\Upsilon=18859$;$\Upsilon\Upsilon=18921$
| $2^{++}$ | $1.00$ | $3.11$ | $19708$ | $\Upsilon\Upsilon=18921$
Table 3: Comparision of our results for the $cc\bar{c}\bar{c}$ systems with other calculations cited.All masses are in unit of MeV. State | $J^{PC}$ | This work | Liu:2019zuc  | Lloyd:2003yc  | Wu:2016vtq  | Chen:2016jxd  | Ader:1981db  | Iwasaki:1975pv  | Karliner:2016zzc  | Barnea:2006sd
---|---|---|---|---|---|---|---|---|---|---
${(cc\bar{c}\bar{c})}$ | $0^{++}$ | $6469$ | $6487$ | $6477$ | $6797$ | $6440-6820$ | $6437$ | $6200$ | $6192$ | $6038-6115$
| $0^{++}$ | $6572$ | $6518$ | $6695$ | $7016$ | $6440-6820$ | $6383$ | $...$ | $...$ | $...$
| $1^{+-}$ | $6519$ | $6500$ | $6528$ | $6899$ | $6370-6510$ | $6437$ | $...$ | $...$ | $6101-6176$
| $2^{++}$ | $6545$ | $6524$ | $6573$ | $6956$ | $6370-6510$ | $6437$ | $...$ | $...$ | $6172-6216$

For fully bottom systems of the tetraquarks $bb\bar{b}\bar{b}$, the solved
results of the model are shown in Table 2. We find that all these
$bb\bar{b}\bar{b}$ states (with $J^{PC}=0^{++}$, $1^{+-}$ and $2^{++}$ ) are
close to each other and strongly unstable as they are far above their two
mesons final states shown. For instance, two of the $0^{++}$ states have the
masses of $19717\,$MeV and $19685\,$MeV (with mass splitting $32\,$MeV). As
seen in Fig 2, the two of the $0^{++}$ states are above thresholds
($\Upsilon\Upsilon$, $\eta_{b}\eta_{b}$) about $764-919\,$MeV. For the
$1^{+-}$ state of $bb\bar{b}\bar{b}$, its mass is higher than the
threshold($\eta_{b}\Upsilon$ and $\Upsilon\Upsilon$) about $780-840\,$MeV. For
the $2^{++}$ state, it is above the threshold ($\Upsilon\Upsilon$), about
$787\,$MeV. By the way, our results for the $bb\bar{b}\bar{b}$ systems are
also compared to other works cited, as shown in Table 4.

Table 4: Comparision of our results for the $bb\bar{b}\bar{b}$ systems with other calculations cited.All masses are in unit of MeV.  State | $J^{PC}$ | This work | Liu:2019zuc  | Wu:2016vtq  | Wang:2017jtz ; Wang:2018poa  | Karliner:2016zzc  | Berezhnoy:2011xn
---|---|---|---|---|---|---|---
${(bb\bar{b}\bar{b})}$ | $0^{++}$ | 19685 | 19322 | 20155 | 18840 | 18826 | 18754
| $0^{++}$ | 19717 | 19338 | 20275 | … | … | …
| $1^{+-}$ | 19700 | 19329 | 20212 | 18840 | … | 18808
| $2^{++}$ | 19708 | 19341 | 20243 | 18850 | … | 18916

###  V.2 The bottom-charmed system($cb\bar{c}\bar{b}$)

Figure 3: Computed masses (MeV the solid lines) of the $cb\bar{c}\bar{b}$
systems of tetraquarks in their ground states, and the thresholds (MeV the
dotted lines) of the two meson final states.

For bottom-charmed systems of the tetraquarks $cb\bar{c}\bar{b}$, we show in
Table 5 the computed results for $R_{0}$, the mixing weights (the CMI
eigenvectors), the tetraquark masses $M(T)$ and thresholds (two mesons), with
the later two plotted in Fig 3. We find that there are four $J^{PC}=0^{++}$
states for the $cb\bar{c}\bar{b}$ systems, all above the thresholds
($B_{c}^{\ast}B_{c}^{\ast}$, $\Upsilon J/\psi$, $B_{c}B_{c}$ and
$\eta_{b}\eta_{c}$) about $424-794\,$MeV. There are four states of the
$cb\bar{c}\bar{b}$ systems with $J^{PC}=1^{+-}$, all highly above the
thresholds ($B_{c}^{\ast}B_{c}^{\ast}$, $B_{c}^{\ast}B_{c}$, $\eta_{b}J/\psi$
and $\Upsilon\eta_{c}$) about $1289-1567\,$MeV, and two states of the
$cb\bar{c}\bar{b}$ systems with $J^{PC}=1^{++}$, all highly above the
thresholds ($B_{c}^{\ast}B_{c}$ and $\Upsilon J/\psi$) about $1347-1417\,$MeV.
There are also two states with $J^{PC}=2^{++}$, both above the thresholds
($B_{c}^{\ast}B_{c}^{\ast}$ and $\Upsilon J/\psi$) about $506-608\,$MeV. This
indicates that the $cb\bar{c}\bar{b}$ systems are unstable against strong
decay to the final states of the mesons.

Table 5: Computed results for the bottom-charmed tetraquark states $cb\bar{c}\bar{b}$. The thresholds of two mesons are also listed. State | $J^{PC}$ | Eigenvector |  $R_{0}$(GeV-1) |  $M(T)$(MeV) | Threshold (MeV)
---|---|---|---|---|---
$cb\bar{c}\bar{b}$ | $0^{++}$ | (-0.21,-0.52,0.82,0.16) | 3.76 | 13076 |  $B_{c}^{*}B_{c}^{*}=12652$;$\Upsilon J/\psi=12557$
|  | (-0.77, 0.24,-0.16,0.58) | 3.90 | 13117 |  $B_{c}B_{c}=12550$;$\eta_{b}\eta_{c}=12382$
|  | (-0.14,-0.82,-0.55,0.01) | 4.02 | 13147 |
|  | (0.59,-0.06,-0.05,0.80) | 4.12 | 13176 |
| $1^{+-}$ | (-0.43,0.39,0.14,0.80) | 3.95 | 13941 |  $B_{c}^{*}B_{c}^{*}=12652$;$B_{c}^{*}B_{c}=12597$
|  | (0.63,0.74,-0.25,0.02) | 4.0 | 13959 |  $\eta_{b}J/\psi=12496$;$\Upsilon\eta_{c}=12444$
|  | (0.57,-0.26,0.71,0.32) | 4.04 | 13966 |
|  | (0.30,-0.48,-0.65,0.51) | 4.18 | 14011 |
| $1^{++}$ | (-0.58,0.82) | 3.94 | 13944 |  $B_{c}^{*}B_{c}=12597$;$\Upsilon J/\psi=12557$
|  | (0.82,0.58) | 4.07 | 13974 |
| $2^{++}$ | (0.74,0.67) | 4.06 | 13158 |  $B_{c}^{*}B_{c}^{*}=12652$;$\Upsilon J/\psi=12557$
|  | (-0.68,0.74) | 4.08 | 13165 |

###  V.3 The strange-heavy systems ($sc\bar{s}\bar{c}$ and $sb\bar{s}\bar{b}$)

Figure 4: Computed masses (Mev the solid lines) of the $sc\bar{s}\bar{c}$
tetraquarks and corresponding two meson thresholds (MeV the dotted lines)
Figure 5: Computed masses (MeV the solid lines) of the $sb\bar{s}\bar{b}$
tetraquarks and corresponding two meson thresholds (MeV the dotted lines)

For strange-charmed systems of the tetraquarks $sc\bar{s}\bar{c}$, we show in
Table 6 the computed results for $R_{0}$, the mixing weights, the masses
$M(T)$ and thresholds (two mesons), with the later two plotted in Fig 4. We
find that there are four $J^{PC}=0^{++}$ states of the $sc\bar{s}\bar{c}$
systems, all below the threshold of $D_{s1}^{\ast}D_{s1}^{\ast}$, in which
three states with masses ($4492,4378,4254$) MeV are above the thresholds of
$D_{s}D_{s}$ and $\phi\left(1020\right)J/\psi$ about $137-556\,$ MeV and
unstable against strong decay to them. The lowest state with mass of
$4091\,$MeV is above the threshold of $D_{s}D_{s}$ about $155\,$MeV while it
is near to the threshold of $\phi(1020)J/\psi$, far below the threshold of
$D_{s1}^{\ast}D_{s1}^{\ast}$. It is uncertain whether the lowest state is
above or below the threshold of $\phi\left(1020\right)J/\psi$ as the model
uncertainty is as large as $\pm 40\,$MeVZhang:2021yul . In the case of the
$J^{PC}=1^{+-}$ states, there are four states, with three of them having the
mass of ($4529,4596,4638$) MeV and all all below the thresholds of
$D_{s}D_{s1}^{\ast}$ and $D_{s1}^{\ast}D_{s1}^{\ast}$ about 110-1031 MeV and
one state, with mass of $4843\,$MeV, above the threshold of
$D_{s}D_{s1}^{\ast}$ about $95\,$MeV but below the threshold of
$D_{s1}^{\ast}D_{s1}^{\ast}$ about $717\,$MeV. There are two
$sc\bar{s}\bar{c}$ systems with $J^{PC}=1^{++}$, both of which are above the
threshold of $\phi\left(1020\right)J/\psi$ about $455-538\,$MeV and below the
threshold of $D_{s}D_{s1}^{\ast}$ about $93-176\,$MeV. There are two
$sc\bar{s}\bar{c}$ systems with $J^{PC}=2^{++}$, both above the threshold of
$\phi\left(1020\right)J/\psi$ about $304-333\,$MeV and below the threshold of
$D_{s1}^{\ast}D_{s1}^{\ast}$ about $1110-1139\,$MeV, unstable to strong decay
to $\phi\left(1020\right)J/\psi$.

For strange-bottom systems $sb\bar{s}\bar{b}$, we show in Table 6 the computed
results for $R_{0}$, the mixing weights, the masses $M(T)$ and thresholds,
with the later two plotted in Fig 5. Similarly, there are four states for each
of $J^{PC}=0^{++}$ and $J^{PC}=1^{+-}$, and two states for each of
$J^{PC}=1^{++}$ and $J^{PC}=2^{++}$. All of the $sb\bar{s}\bar{b}$ systems are
above the thresholds except for the lowest one with mass of $10843$ MeV which
is near to thresholds ($13\,$MeV) of the $B_{s}^{\ast}B_{s}^{\ast}$. The
$0^{++}$ states of the $sb\bar{s}\bar{b}$ systems are above the thresholds of
$B_{s}^{0}B_{s}^{0}$, $B_{s}^{\ast}B_{s}^{\ast}$ and
$\phi\left(1020\right)\Upsilon$. Among them, the minimum mass of $10843\,$MeV
can be strongly decayed into $B_{s}^{0}B_{s}^{0}$ and
$\phi\left(1020\right)\Upsilon$. Because of the error in the model, it is
uncertain whether it is above or below the threshold of
$B_{s}^{\ast}B_{s}^{\ast}$. The four $J^{PC}=1^{+-}$ states are all highly
above the thresholds of $B_{s}^{0}B_{s}^{\ast}$, $B_{s}^{\ast}B_{s}^{\ast}$
(about $562-945\,$MeV and $514-897\,$MeV, respectively).There are two states
$J^{PC}=1^{++}$, which are higher than $B_{s}^{0}B_{s}^{\ast}$ and
$\phi\left(1020\right)\Upsilon$(about $626-666\,$MeV and $928-968\,$MeV,
respectively). $J^{PC}=2^{++}$ has two states, which are higher than
$B_{s}^{\ast}B_{s}^{\ast}$ and $\phi\left(1020\right)\Upsilon$ thresholds
(about $263-324\,$MeV and $613-674\,$MeV), indicating they are unstable.

Table 6: Computed results for the strange-heavy tetraquark states $sc\bar{s}\bar{c}$ and $sb\bar{s}\bar{b}$. The thresholds of two mesons are also listed.  State | $J^{PC}$ | Eigenvector |  $R_{0}$(GeV-1) |  $M(T)$(MeV) | Threshold (MeV)
---|---|---|---|---|---
$sc\bar{s}\bar{c}$ | $0^{++}$ | $(-0.18,-0.51,0.83,0.14)$ | $4.70$ | $4091$ |  $D_{s}D_{s}=3936$;$D_{s1}^{\ast}D_{s1}^{\ast}=5560$
|  | $(-0.76,0.21,-0.15,0.60)$ | $4.93$ | $4254$ | $\phi\left(1020\right)J/\psi=4117$
|  | $(-0.12,-0.84,-0.54,0.01)$ | $5.11$ | $4378$ |
|  | $(0.62,-0.06,-0.04,0.78)$ | $5.33$ | $4492$ |
| $1^{+-}$ | $(-0.33,0.57,0.03,0.75)$ | $5.21$ | $4529$ |  $D_{s}D_{s1}^{\ast}=4748$;$D_{s1}^{\ast}D_{s1}^{\ast}=5560$
|  | $(0.87,0.49,0.08,0.02)$ | $5.30$ | $4596$ |
|  | $(0.18,-0.46,0.77,0.41)$ | $5.38$ | $4638$ |
|  | $(0.32,-0.47,-0.64,0.52)$ | $5.46$ | $4843$ |
| $1^{++}$ | $(-0.58,0.82)$ | $5.22$ | $4572$ |  $\phi\left(1020\right)J/\psi=4117$;$D_{s}D_{s1}^{\ast}=4748$
|  | $(0.82,0.58)$ | $5.41$ | $4655$ |
| $2^{++}$ | $(0.55,0.83)$ | $5.39$ | $4421$ |  $D_{s1}^{\ast}D_{s1}^{\ast}=5560$;$\phi\left(1020\right)J/\psi=4117$
|  | $(-0.83,0.55)$ | $5.39$ | $4450$ |
$sb\bar{s}\bar{b}$ | $0^{++}$ | $(-0.36,-0.35,0.80,0.33)$ | $4.43$ | $10843$ |  $B_{s}^{0}B_{s}^{0}=10734$;$B_{s}^{\ast}B_{s}^{\ast}=10830$
|  | $(-0.58,0.39,-0.36,0.62)$ | $4.60$ | $11023$ | $\phi\left(1020\right)\Upsilon=10480$
|  | $(0.35,0.81,0.47,0.10)$ | $4.74$ | $11111$ |
|  | $(0.64,-0.30,-0.13,0.70)$ | $4.86$ | $11158$ |
| $1^{+-}$ | $(-0.41,0.67,-0.20,0.58)$ | $4.88$ | $11344$ |  $B_{s}^{0}B_{s}^{\ast}=10782$;$B_{s}^{\ast}B_{s}^{\ast}=10830$
|  | $(0.76,0.39,0.45,0.26)$ | $4.94$ | $11388$ |
|  | $(-0.31,-0.47,0.63,0.53)$ | $4.88$ | $11457$ |
|  | $(0.38,-0.43,-0.59,0.56)$ | $5.11$ | $11727$ |
| $1^{++}$ | $(0.82,0.58)$ | $4.98$ | $11408$ |  $B_{s}^{0}B_{s}^{\ast}=10782$;$\phi\left(1020\right)\Upsilon=10480$
|  | $(-0.58,0.82)$ | $4.85$ | $11448$ |
| $2^{++}$ | $(0.51,0.86)$ | $4.97$ | $11093$ |  $B_{s}^{\ast}B_{s}^{\ast}=10830$;$\phi\left(1020\right)\Upsilon=10480$
|  | $(-0.86,0.51)$ | $4.99$ | $11154$ |

###  V.4 The heavy-light(non-strange) systems ($nc\bar{n}\bar{c}$ and
$nb\bar{n}\bar{b}$)

Figure 6: Computed masses(MeV the solid lines) of the hidden-bottom
tetraquarks $nc\bar{n}\bar{c}$ and the two meson thresholds (MeV the dotted
lines).  Figure 7: Computed masses(MeV the solid lines) of the hidden-bottom
tetraquark $nb\bar{n}\bar{b}$ and the two meson thresholds (MeV the dotted
lines).

For hidden charmed systems of the tetraquarks $nc\bar{n}\bar{c}$, we show the
computed results for $R_{0}$, the mixing weights, the masses $M(T)$ and
thresholds in Table 7, with the later two plotted in Fig 6. There are four
states for each of $J^{PC}=0^{++}$ and $J^{PC}=1^{+-}$, and two states for
each of $J^{PC}=1^{++}$ and $J^{PC}=2^{++}$. For the $0^{++}$ states, two
higher states ($4259$ MeV, $4127\,$MeV) are all above the thresholds of
$D^{0}D^{0}$, $D^{\ast}D^{\ast}$, $\omega\left(782\right)J/\psi$ and
$\pi^{0}\eta_{c}$ (about $397-529\,$MeV, $110-242\,$MeV, $248-380\,$MeV and
$1008-1140\,$MeV). The lower state with mass $3954\,$MeV, which is above the
thresholds of $D^{0}D^{0}$, $\omega\left(782\right)J/\psi$, $\pi^{0}\eta_{c}$
and below the threshold of $D^{\ast}D^{\ast}$, can strongly decay to the three
former final states. The lowest state, which is below the thresholds of
$D^{\ast}D^{\ast}$, $\omega\left(782\right)J/\psi$,$D^{0}D^{0}$ and above the
thresholds of $\pi^{0}\eta_{c}$, can decay to two final states of
$\pi^{0}\eta_{c}$ . Further, all states with $J^{PC}=1^{+-}$, $1^{++}$ and
$2^{++}$ are above the thresholds of $D^{0}D^{\ast}$, $D^{\ast}D^{\ast}$,
$\pi^{0}J/\psi$, $\omega\left(782\right)J/\psi$,
$\omega\left(782\right)\eta_{c}$, can decay to the laters with same quantum
numbers. For instance, the $1^{++}$ states can decay to
$\omega\left(782\right)J/\psi$ and $D^{0}D^{\ast}$, the $2^{++}$ states can
decay to $D^{\ast}D^{\ast}$, $\omega\left(782\right)$$J/\psi$.

For hidden-bottom systems of tetraquarks $nb\bar{n}\bar{b}$, we show the
computed results for $R_{0}$, the mixing weights, the masses $M(T)$ and
thresholds in Table 7, with the later two plotted in Fig 7. We find from Fig 7
that all states of $nb\bar{n}\bar{b}$ systems are above the thresholds of
their final states of two mesons, except for the lowest state ($10484\,$MeV),
which is below the thresholds of $B^{0}B^{0}$ and $B^{\ast}B^{\ast}$ only and
it can decay into $\omega\left(782\right)\Upsilon$, $\pi^{0}\eta_{b}$. The
possible decays are, for instance, the $nb\bar{n}\bar{b}(0^{++})$ to
$B^{0}B^{0}$, $B^{\ast}B^{\ast}$, $\omega\left(782\right)\Upsilon$ and
$\pi^{0}\eta_{b}$, the $nb\bar{n}\bar{b}(1^{+-})$ to $B^{0}B^{\ast}$,
$B^{\ast}B^{\ast}$, $\pi^{0}\Upsilon$ and $\omega\left(782\right)\eta_{b}$,
the $nb\bar{n}\bar{b}(1^{++})$ to $B^{0}B^{\ast}$,
$\omega\left(782\right)\Upsilon$, the $nb\bar{n}\bar{b}(2^{++})$ to
$B^{\ast}B^{\ast}$, $\omega\left(782\right)\Upsilon$.

Table 7: Computed results for the hidden heavy-flavor tetraquark $nc\bar{n}\bar{c}$ and $nb\bar{n}\bar{b}$, with respective thresholds of two mesons shown also.  State | $J^{PC}$ | Eigenvector |  $R_{0}$(GeV-1) |  $M(T)$(MeV) | Threshold (MeV)
---|---|---|---|---|---
$nc\bar{n}\bar{c}$ | $0^{++}$ | $(-0.24,-0.45,0.83,0.20)$ | $4.90$ | $3715$ |  $D^{0}D^{0}=3730$;$D^{\ast}D^{\ast}=4017$
|  | $(-0.71,0.26,-0.22,0.62)$ | $5.01$ | $3954$ | $\omega\left(782\right)J/\psi=3879$
|  | $(0.17,0.85,0.50,0.00)$ | $5.17$ | $4127$ | $\pi^{0}\eta_{c}=3119$
|  | $(0.65,-0.11,-0.06,0.75)$ | $5.38$ | $4259$ |
| $1^{+-}$ | $(-0.31,0.65,-0.06,0.69)$ | $5.15$ | $4079$ |  $D^{0}D^{\ast}=3874$;$D^{\ast}D^{\ast}=4017$
|  | $(0.87,0.33,0.34,0.11)$ | $5.26$ | $4172$ |  $\pi^{0}J/\psi=3264$;$\omega\left(782\right)\eta_{c}=3767$
|  | $(-0.15,-0.50,0.70,0.48)$ | $5.33$ | $4263$ |
|  | $(0.35,-0.46,-0.61,0.54)$ | $5.32$ | $4626$ |
| $1^{++}$ | $(-0.58,0.82)$ | $5.18$ | $4210$ |  $D^{0}D^{\ast}=3874$;$\omega\left(782\right)J/\psi=3879$
|  | $(0.82,0.58)$ | $5.36$ | $4233$ |
| $2^{++}$ | $(0.46,0.89)$ | $5.34$ | $4152$ |  $D^{\ast}D^{\ast}=4017$;$\omega\left(782\right)J/\psi=3879$
|  | $(-0.89,0.46)$ | $5.31$ | $4219$ |
$nb\bar{n}\bar{b}$ | $0^{++}$ | $(-0.37,-0.31,0.80,0.36)$ | $4.65$ | $10484$ |  $B^{0}B^{0}=10560$;$B^{\ast}B^{\ast}=10650$
|  | $(-0.51,0.37,-0.40,0.67)$ | $4.70$ | $10760$ | $\omega\left(782\right)\Upsilon=10242$
|  | $(0.42,0.78,0.43,0.14)$ | $4.82$ | $10887$ | $\pi^{0}\eta_{b}=9533$
|  | $(0.65,-0.40,-0.14,0.63)$ | $4.90$ | $10943$ |
| $1^{+-}$ | $(-0.43,0.67,-0.25,0.54)$ | $4.79$ | $10865$ |  $B^{0}B^{\ast}=10605$;$B^{\ast}B^{\ast}=10650$
|  | $(0.73,0.41,0.47,0.29)$ | $4.84$ | $10925$ |  $\pi^{0}\Upsilon=9594$;$\omega\left(782\right)\eta_{b}=10182$
|  | $(-0.36,-0.45,0.61,0.55)$ | $4.81$ | $11097$ |
|  | $(0.39,-0.42,-0.59,0.57)$ | $4.95$ | $11509$ |
| $1^{++}$ | $(0.82,0.58)$ | $4.89$ | $10949$ |  $B^{0}B^{\ast}=10605$;$\omega\left(782\right)\Upsilon=10242$
|  | $(-0.58,0.82)$ | $4.78$ | $11090$ |
| $2^{++}$ | $(0.45,0.89)$ | $4.88$ | $10835$ |  $B^{\ast}B^{\ast}=10650$;$\omega\left(782\right)\Upsilon=10242$
|  | $(-0.89,0.45)$ | $4.88$ | $10944$ |

##  VI Summary

Stimulated by observations of the $X(6900)$ by LHCb and the recent
observations of the $X(6600)$ by CMS and ATLAS experiments of the LHC, we have
systematically investigated the ground-state masses of hidden heavy-flavor
tetraquarks with two and four hidden heavy-flavor within a unified framework
of MIT bag model which incorporates chromomagnetic interactions and enhanced
binding energy. Based on color-spin wavefunctions constructed for the hidden
heavy-flavor tetraquarks, we solve the MIT bag model and diagonalize the
chromomagnetic interaction (CMI) to predict masses of the color-spin
multiplets of hidden heavy-flavor tetraquarks in their ground states with
spin-parity quantum numbers $J^{PC}=0^{++}$, $1^{++}$, $2^{++}$, and $1^{+-}$.
We find that the fully charmed tetraquark $cc\bar{c}\bar{c}$ with
$J^{PC}=0^{++}$ has mass about $6572$ MeV and is very likely to be the
$X(6600)$ reported by CMS and ATLAS experiments of the LHC, with the measured
mass $6552\pm 10\pm 12$ MeV. We further computed masses of the tetraquark
systems $bb\bar{b}\bar{b}$, $cb\bar{c}\bar{b}$, $sc\bar{s}\bar{c}$,
$sb\bar{s}\bar{b}$, $nc\bar{n}\bar{c}$ and $nb\bar{n}\bar{b}$ in their color-
spin multiplets and suggested that the particle $Z_{c}(4200)$ reported by
Belle:2014nuw  is likely to be the hidden-charm tetraquark made of
$nc\bar{n}\bar{c}$ with $J^{PC}=1^{+-}$.

Compared to two-meson thresholds determined via the final states in details,
the most-likely strong decay channels are noted. Our mass computation shows
that all of these hidden heavy-flavor tetraquarks are above the thresholds of
the lowest two-mesons final states and unstable against strong decay to these
final states. For the doubly heavy systems of the tetraquarks
$sb\bar{s}\bar{b}$, $sc\bar{s}\bar{c}$, $nb\bar{n}\bar{b}$ and
$nc\bar{n}\bar{c}$, there are a few states below thresholds except for their
lowest final states, indicating that they may have longer lifetime compared to
the fully heavy tetraquarks. We also find some near-threshold states for which
coupled channel effects are possible. We hope that upcoming LHCb experiments
with increased data can test the prediction in this work.

Acknowledgments

D. J. is supported by the National Natural Science Foundation of China under
the no. 12165017.

## Appendix A

Based on the color $SU(3)_{c}$ symmetry, one can obtain two components of
color singlets $6_{c}\otimes\bar{6}_{c}$ and $\bar{3}_{c}\otimes 3_{c}$ for
the hidden-flavor tetraquarks,

$\displaystyle\phi_{1}^{T}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{6}}\left(rr\bar{r}\bar{r}+gg\bar{g}\bar{g}+bb\bar{b}\bar{b}\right)$
(30)
$\displaystyle+\frac{1}{2\sqrt{6}}\left(rb\bar{b}\bar{r}+br\bar{b}\bar{r}+gr\bar{g}\bar{r}+rg\bar{g}\bar{r}+gb\bar{b}\bar{g}+bg\bar{b}\bar{g}\right.$
$\displaystyle+gr\bar{r}\bar{g}+rg\bar{r}\bar{g}+gb\bar{g}\bar{b}\left.+bg\bar{g}\bar{b}+rb\bar{r}\bar{b}+br\bar{r}\bar{b}\right),$
$\displaystyle\phi_{2}^{T}$ $\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{3}}\left(rb\bar{b}\bar{r}-br\bar{b}\bar{r}-gr\bar{g}\bar{r}+rg\bar{g}\bar{r}+gb\bar{b}\bar{g}-bg\bar{b}\bar{g}\right.$
(31)
$\displaystyle\left.+gr\bar{r}\bar{g}-rg\bar{r}\bar{g}-gb\bar{g}\bar{b}+bg\bar{g}\bar{b}-rb\bar{r}\bar{b}+br\bar{r}\bar{b}\right),$

which corresponds to two color configurations in Eq. (1).

For six states $\chi_{1\sim 6}^{T}$ (2) of heavy tetraquarks, one can
construct their spin wave functions via writing the $CG$ coefficients
explicitly:

$\displaystyle\chi_{1}^{T}$ $\displaystyle=\uparrow\uparrow\uparrow\uparrow,$
$\displaystyle\chi_{2}^{T}$
$\displaystyle=\frac{1}{2}\left(\uparrow\uparrow\uparrow\downarrow+\uparrow\uparrow\downarrow\uparrow-\uparrow\downarrow\uparrow\uparrow-\downarrow\uparrow\uparrow\uparrow\right),$
$\displaystyle\chi_{3}^{T}$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\uparrow\uparrow\downarrow\downarrow+\downarrow\downarrow\uparrow\uparrow\right),$
$\displaystyle-\frac{1}{2\sqrt{3}}\left(\uparrow\downarrow\uparrow\downarrow+\uparrow\downarrow\downarrow\uparrow+\downarrow\uparrow\uparrow\downarrow+\downarrow\uparrow\downarrow\uparrow\right),$
$\displaystyle\chi_{4}^{T}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\uparrow\uparrow\uparrow\downarrow-\uparrow\uparrow\downarrow\uparrow\right),$
$\displaystyle\chi_{5}^{T}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\uparrow\downarrow\uparrow\uparrow-\downarrow\uparrow\uparrow\uparrow\right),$
$\displaystyle\chi_{6}^{T}$
$\displaystyle=\frac{1}{2}\left(\uparrow\downarrow\uparrow\downarrow-\uparrow\downarrow\downarrow\uparrow-\downarrow\uparrow\uparrow\downarrow+\downarrow\uparrow\downarrow\uparrow\right),$
(32)

in which notations $\uparrow$ and $\downarrow$ represent the third component
of the quark’s spin. Alternatively, one can also use the $CG$ coefficients
given in Ref. Zhang:2021yul  to examine Eq. (32) for the different spin
states. Note that the results of spin factors in Ref. Zhang:2021yul  are shown
as matrix form in the spanned space of the states $\chi_{1\sim 6}^{T}$ as the
spin multiplets (2) indicated. Combining with two color configurations
$\phi_{1\sim 2}^{T}$ in Eq. (1) and six spin configurations $\chi_{1\sim
6}^{T}$ in Eq. (2), one can then construct their color-spin wavefunctions (3).
The allowed states of the hidden-flavor tetraquarks to be mix due to
chromomagnetic interaction are listed in Table 1.

## References

  * (1) Murray Gell-Mann.  A Schematic Model of Baryons and Mesons.  Phys. Lett., 8:214–215, 1964.
  * (2) G. Zweig.  An SU(3) model for strong interaction symmetry and its breaking. Version 1.  1 1964.
  * (3) Robert L. Jaffe.  Multi-Quark Hadrons. 1. The Phenomenology of (2 Quark 2 anti-Quark) Mesons.  Phys. Rev. D, 15:267, 1977.
  * (4) Robert L. Jaffe.  Multi-Quark Hadrons. 2. Methods.  Phys. Rev. D, 15:281, 1977.
  * (5) S. K. Choi et al.  Observation of a narrow charmonium-like state in exclusive $B^{\pm}\to K^{\pm}\pi^{+}\pi^{-}J/\psi$ decays.  Phys. Rev. Lett., 91:262001, 2003.
  * (6) M. Ablikim et al.  Observation of a Charged Charmoniumlike Structure in $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ at $\sqrt{s}$ =4.26 GeV.  Phys. Rev. Lett., 110:252001, 2013.
  * (7) K. Chilikin et al.  Observation of a new charged charmoniumlike state in $\bar{B}^{0}\to J/\psi K^{-}\pi^{+}$ decays.  Phys. Rev. D, 90(11):112009, 2014.
  * (8) S. K. Choi et al.  Observation of a resonance-like structure in the $pi^{\pm}\psi^{\prime}$ mass distribution in exclusive $B\to K\pi^{\pm}\psi^{\prime}$ decays.  Phys. Rev. Lett., 100:142001, 2008.
  * (9) K. Chilikin et al.  Experimental constraints on the spin and parity of the $Z$(4430)+.  Phys. Rev. D, 88(7):074026, 2013.
  * (10) Roel Aaij et al.  Observation of the resonant character of the $Z(4430)^{-}$ state.  Phys. Rev. Lett., 112(22):222002, 2014.
  * (11) Roel Aaij et al.  Model-independent confirmation of the $Z(4430)^{-}$ state.  Phys. Rev. D, 92(11):112009, 2015.
  * (12) Evelina Bouhova-Thacker.  ATLAS results on exotic hadronic resonances.  PoS, ICHEP2022:806, 2022.
  * (13) Jingqing Zhang and Kai Yi.  Recent CMS results on exotic resonances.  PoS, ICHEP2022:775, 11 2022.
  * (14) Roel Aaij et al.  Observation of structure in the $J/\psi$ -pair mass spectrum.  Sci. Bull., 65(23):1983–1993, 2020.
  * (15) Wen-Xuan Zhang, Hao Xu, and Duojie Jia.  Masses and magnetic moments of hadrons with one and two open heavy quarks: Heavy baryons and tetraquarks.  Phys. Rev. D, 104(11):114011, 2021.
  * (16) Marek Karliner and Jonathan L. Rosner.  Baryons with two heavy quarks: Masses, production, decays, and detection.  Phys. Rev. D, 90(9):094007, 2014.
  * (17) Marek Karliner and Jonathan L. Rosner.  Quark-level analogue of nuclear fusion with doubly-heavy baryons.  Nature, 551:89, 2017.
  * (18) Thomas A. DeGrand, R. L. Jaffe, K. Johnson, and J. E. Kiskis.  Masses and Other Parameters of the Light Hadrons.  Phys. Rev. D, 12:2060, 1975.
  * (19) Ming-Sheng Liu, Qi-Fang Lü, Xian-Hui Zhong, and Qiang Zhao.  All-heavy tetraquarks.  Phys. Rev. D, 100(1):016006, 2019.
  * (20) Richard J. Lloyd and James P. Vary.  All charm tetraquarks.  Phys. Rev. D, 70:014009, 2004.
  * (21) Jing Wu, Yan-Rui Liu, Kan Chen, Xiang Liu, and Shi-Lin Zhu.  Heavy-flavored tetraquark states with the $QQ\bar{Q}\bar{Q}$ configuration.  Phys. Rev. D, 97(9):094015, 2018.
  * (22) Wei Chen, Hua-Xing Chen, Xiang Liu, T. G. Steele, and Shi-Lin Zhu.  Hunting for exotic doubly hidden-charm/bottom tetraquark states.  Phys. Lett. B, 773:247–251, 2017.
  * (23) J. P. Ader, J. M. Richard, and P. Taxil.  DO NARROW HEAVY MULTI - QUARK STATES EXIST?  Phys. Rev. D, 25:2370, 1982.
  * (24) Yoichi Iwasaki.  A Possible Model for New Resonances-Exotics and Hidden Charm.  Prog. Theor. Phys., 54:492, 1975.
  * (25) Marek Karliner, Shmuel Nussinov, and Jonathan L. Rosner.  $QQ\bar{Q}\bar{Q}$ states: masses, production, and decays.  Phys. Rev. D, 95(3):034011, 2017.
  * (26) N. Barnea, J. Vijande, and A. Valcarce.  Four-quark spectroscopy within the hyperspherical formalism.  Phys. Rev. D, 73:054004, 2006.
  * (27) Zhi-Gang Wang.  Analysis of the $QQ\bar{Q}\bar{Q}$ tetraquark states with QCD sum rules.  Eur. Phys. J. C, 77(7):432, 2017.
  * (28) Zhi-Gang Wang and Zun-Yan Di.  Analysis of the vector and axialvector $QQ\bar{Q}\bar{Q}$ tetraquark states with QCD sum rules.  Acta Phys. Polon. B, 50:1335, 2019.
  * (29) A. V. Berezhnoy, A. V. Luchinsky, and A. A. Novoselov.  Tetraquarks Composed of 4 Heavy Quarks.  Phys. Rev. D, 86:034004, 2012.


# Theory of Josephson current on a lattice model of grain boundary in $d$-wave
superconductors

Takashi Sakamori1, Satoshi Kashiwaya1, Rikizo Yano1,Yukio Tanaka1,Takafumi
Hatano2 and Keiji Yada1 1 Department of Applied Physics, Nagoya University,
Nagoya 464-8603, Japan  2 Department of Materials Physics, Nagoya University,
Nagoya 464-8603, Japan<EMAIL_ADDRESS>
###### Abstract

Identifying the origins of suppression of the critical current at grain
boundaries of high-critical-temperature superconductors, such as cuprates and
iron-based superconductors, is a crucial issue to be solved for future
applications with polycrystalline materials. Although the dominant factor of
current suppression might arise during material fabrication and/or processing,
investigating it due to an internal phase change of the pair potential is an
important issue in understanding the threshold of the critical current. In
this paper, we study the Josephson current on a symmetric [001]-tilt grain
boundary (GB) of a $d$-wave superconductor on a lattice model. In addition to
the suppression of the maximum Josephson current associated with the internal
phase change of the $d$-wave pair potential which has been predicted in
continuum models, we find a unique phase interference effect due to folding of
the Fermi surface in the lattice model. In particular, the resultant maximum
Josephson current at low-tilting-angle regions tends to be suppressed more
than that in preexisting theories. Because similar suppressions of the
critical current at GBs have been reported in several experimental works, the
present model can serve as a guide to clarify the complicated transport
mechanism in GBs.

††: Supercond. Sci. Technol.

  * March 2023

Keywords: Josephson junction, HTSC, grain boundary, $d$-wave supercondunctor

##  1 Introduction

The application of high-critical-temperature superconductors (HTSCs), such as
cuprates and iron-based superconductors (IBSs), to superconducting cables and
joints has been the subject of research. It requires realizing a stable
persistent current even if the circuit includes the grain boundaries (GBs) of
the crystals because these boundaries inevitably exist inside cables and at
their joints. Despite the successful development of superconducting cables
based on low-critical-temperature (low-$T_{c}$) superconductors, the
realization of high-current-density cables based on HTSCs requires additional
techniques for crystalline orientation control, such as ion-beam-assisted
deposition. The critical current flow through GBs has been experimentally
confirmed to be seriously suppressed when the tilting angle of the GB exceeds
approximately 5∘ [1, 2, 3, 4] for cuprates and 10∘ for IBSs [5, 6, 7, 8]. Two
origins of this severe suppression of the critical current have been proposed.
One proposed origin is related to material issues, such as crystal growth,
doping inhomogeneities, and impurity scattering around GBs [9, 10, 11]. The
other proposed origin is related to the anisotropic pairing states of
superconductors. Unlike conventional low-$T_{c}$ superconductors, the pair
potentials of cuprates and IBSs exhibit phase changes in the momentum space.
According to the current–phase relation (CPR) of the Josephson current, the
internal phase change contributes to the reversal of the current direction.
Although the dominant origin of the suppression might be the former in actual
cables and joints, the influences of the latter should also be clarified in
detail.[12]

Focusing on cuprates, we note that the Josephson current formula for $d$-wave
superconductor junctions has been developed mainly on the basis of the
continuum theory. Sigrist and Rice [13] investigated the influence of the
anisotropic pair potential on the Josephson current and proposed a theoretical
reversal of the Josephson current direction due to the internal phase change.
A reversal of the current’s direction accompanies the transition of a
0-junction to a $\pi$-junction. This formula successfully explains the
anomalous magnetic field dependence of superconducting quantum interference
devices and the spontaneous magnetic flux in tri-crystals [14, 15, 16, 17].
These experimental results have been accepted as strong evidence for $d$-wave
pairing symmetry of cuprates.

Another critical effect expected for $d$-wave Josephson junctions is the
formation of zero-energy Andreev bound states (ZEABSs) at the interface [18,
19, 20, 21]. The presence of ZEABSs is predicted to enhance the higher
harmonics of the CPR of the Josephson current [22, 23, 24]. The nonmonotonic
temperature dependence of the critical current [25, 26], the nonsinusoidal CPR
[27], and the higher harmonics of the CPR [28] due to the ZEABS have been
detected experimentally in GB junctions of cuprates, supporting the validity
of the theoretical models. The character of the CPR determines the free-energy
minimum because the CPR is given by the derivative of the free energy with
respect to the phase difference $\phi$. When the free-energy minimum locates
at $\phi=0$ ($\pi$), it is referred to as a 0 ($\pi$)-junction. In the
presence of higher harmonics, the free-energy minimum might appear at neither
$\phi=0$ nor $\pi$ but at $\phi=\phi_{m}$; this junction is classified as a
$\phi$-junction[22, 23]

The dependence of the maximum Josephson current on the tilting angle in HTSC
GBs is an unresolved issue. The experimentally detected suppression of the
critical current at GBs [2] is more drastic than the suppression predicted by
the continuum models. The influences of the short coherence length have been
introduced in the calculation using lattice models [29, 30]. Although the
nonmonotonic temperature dependence of the critical current and the existence
of the $\varphi$-junction was reproduced in a lattice model [31], the tilting-
angle dependence was inconsistent with the experimental observations. These
facts indicate that the suppression mechanism needs to be clarified beyond the
mechanisms treated in the preexisting models to explain the Josephson current
through the GBs.

In the present work, we study symmetric [001]-tilt GB Josephson junctions
using a lattice model in which the sites at the interface are shared by the
lattices on the left and right sides (see figure 1). In this model, the
periodicity of the unit cell $m$ in the direction parallel to the interface
changes in accordance with the tilting angle and the corresponding Brillouin
zone (BZ) parallel to the interface is folded. When the tilting angle is
reduced, folding of the BZ becomes significant with increasing $m$. This
effect strongly enhances the interference of the internal phase change of the
anisotropic pair potential. We found that, as a result of the enhanced
interference, the Josephson current is substantially suppressed with
increasing $m$.

The remainder of this paper is organized as follows. In section 2, we explain
the model and formulation. In section 3, the results of numerical calculations
are presented. In section 4, we summarize the obtained results and clarify the
physics in the proposed model.

##  2 Model and formulation

We consider symmetric [001]-tilt GBs as shown in figure 1, where ($m10$) and
($\bar{m}10$) surfaces are connected by sharing the most end sites. Hereafter,
this structure is referred to as the ($m10$) GB. The ($m10$) and ($\bar{m}10$)
structures are tilted by $\pm\theta/2$ from the untilted square lattice with
coordinates $(ia,ja)$ at the lattice points, where $i$ and $j$ are integers
and $a$ is the lattice constant. The coordinate of the lattice point $x<0$ is
$ia\left(\cos\frac{\theta}{2},\sin\frac{\theta}{2}\right)+ja\left(-\sin\frac{\theta}{2},\cos\frac{\theta}{2}\right)$;
for $x>0$, it is
$ia\left(\cos\frac{\theta}{2},-\sin\frac{\theta}{2}\right)+ja\left(\sin\frac{\theta}{2},\cos\frac{\theta}{2}\right)$.
We consider the ($m10$) structure for $x>0$. As shown in figure 1, the ($m10$)
structure has discrete translational symmetry along the $y$ direction by
$\sqrt{m^{2}+1}a$. That is, the $x$ coordinate of the lattice point at
$i=i_{0}$ and $j=j_{0}$ is equivalent that at $i=i_{0}-p$ and $j=j_{0}-mq$
with integers $p$ and $q$ that index the unit cell. Thus, we consider the unit
cell composed of $m$ lined lattice points along the $y_{R(L)}$ direction. Two
translation vectors $\hat{u}_{R(L)}$ and $\hat{v}$ are then given by

$\displaystyle\hat{u}_{R(L)}$
$\displaystyle=\left(a\left(\sin\frac{\theta}{2}+\cos\frac{\theta}{2}\right),(-)a\left(\sin\frac{\theta}{2}-\cos\frac{\theta}{2}\right)\right)$
$\displaystyle\hat{v}$ $\displaystyle=\left(0,a\sqrt{m^{2}+1}\right).$ (1)

We label the lattice point by three integers $o$, $p$, and $q$; its coordinate
$r(o,p,q)$ is

$\displaystyle\hat{r}(o,p,q)=o\left(\sin\frac{\theta}{2},-\mathrm{sgn}(p)\cos\frac{\theta}{2}\right)+p\hat{u}+q\hat{v},$
(2)

where $o$ is the sublattice index, which takes values from 0 to $m$ for $p\geq
1$, and from 0 to $-m$ for $p\leq-1$. For $p=0$, $o$ takes values from $-m$ to
$m$. The corresponding lattice sites for the (210) GB are shown in figure 1.
To consider a superconductor/normal metal/superconductor Josephson junction,
we assume that the left (blue) and right (red) sites are in the
superconducting states with the phases of $\phi_{L}$ and $\phi_{R}$,
respectively, and that the center (green) sites with $p=0$ are in the normal
state. The phase difference between the left and right $\phi$ is given by
$\phi_{L}-\phi_{R}$. For $p=0$, the ($m10$) structure at $x<0$ and the
($\bar{m}10$) structure at $x>0$ share the same lattice point with $o=0$. The
resultant unit cell of the normal state has a dog-legged shape with $2m+1$
lattice points.

Figure 1: Lattice structure of the (210) GB. Blue, red, and green solid
circles denote the lattice points of the left superconductor, right
superconductor, and the normal metal site, respectively. The
$x_{L,(R)},y_{L(R)}$ axes are in a tilted geometry.

For this lattice structure, we consider the tight-binding model. The original
bulk Hamiltonian without tilting is a single-band model with a $d$-wave pair
potential given by

$\displaystyle\mathcal{H}_{R(L)}=$
$\displaystyle\sum_{k}(\varepsilon_{k}-\mu)c_{k}^{\dagger}c_{k}+\sum_{k}\Delta_{k}c_{k}^{\dagger}c_{-k}^{\dagger}+h.c.,$
(3) $\displaystyle\varepsilon_{k}=$ $\displaystyle-2t(\cos k_{x}+\cos k_{y}),$
(4) $\displaystyle\Delta_{k}^{R(L)}=$ $\displaystyle\Delta_{0}(\cos k_{x}-\cos
k_{y})e^{i\phi_{R(L)}},$ (5)

where $t$ and $\Delta_{0}$ are the hopping integral and the pair potential,
respectively. The matrix components of the Hamiltonian in Bogoliubov–de Gennes
form in real space for $x\geq 0$ and $1\leq o\leq m$ are then

$\displaystyle\langle o,p,q|\mathcal{H}|o-1,p+1,q\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-t&\Delta(p)\\\ \Delta^{*}(p)&t\end{pmatrix},$
(6) $\displaystyle\langle o-1,p,q|\mathcal{H}|o,p,q\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-t&-\Delta(p)\\\ -\Delta^{*}(p)&t\end{pmatrix},$
(7)

and

$\displaystyle\langle o=0,p,q|\mathcal{H}|o^{\prime}=m,p,q+1\rangle$
$\displaystyle=$ $\displaystyle\begin{pmatrix}-t&\Delta(p)\\\
\Delta^{*}(p)&t\end{pmatrix},$ (8) $\displaystyle\langle
o=m,p,q|\mathcal{H}|o^{\prime}=0,p+1,q-1\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-t&-\Delta(p)\\\ -\Delta^{*}(p)&t\end{pmatrix}.$
(9)

Similarly, we can consider corresponding matrix components for the
($\bar{m}10$) structure in $x<0$. Here, because the layer with $p=0$ is in the
normal state,

$\displaystyle\Delta(p)=\left\\{\begin{matrix}e^{i\phi_{R}}\Delta_{0}/2&(p\geq
1)\\\ 0&(p=0)\\\ e^{i\phi_{L}}\Delta_{0}/2&(p\leq 1)\\\ \end{matrix},\right.$
(10)

Using these matrices, we calculate the Green’s function at the center sites
and obtain the Josephson current through the normal-state sites.

To calculate the Josephson current, we define the current operator $J$ from
the equation of continuity in the lattice model:

$\displaystyle
e\frac{d}{dt}\\{c^{\dagger}_{ky,0,0}(t)c_{ky,0,0}(t)\\}=J_{LC}(k_{y})-J_{CR}(k_{y}).$
(11)

where $c_{ky,p,o}$ ($c_{ky,p,o}^{\dagger}$) is the annihilation (creation)
operator at the site of sublattice $o$ in the $p$th layer with momentum
$k_{y}$. The left-hand side of Eq. (11) is obtained by the Heisenberg
equation, and we obtain the current operator $J_{CR}(k_{y})$:

$\displaystyle
J_{CR}(k_{y})=\frac{iet}{\hbar}(c^{\dagger}_{ky,0,1}c_{ky,0,0}-c^{\dagger}_{ky,0,0}c_{ky,0,1}+e^{ik_{y}}c^{\dagger}_{ky,0,m}c_{ky,0,0}-e^{-ik_{y}}c^{\dagger}_{ky,0,0}c_{ky,0,m}),$
(12)

The expectation value of $J$ is calculated by the Green’s function at
$\tau=0$, where $\tau$ is the imaginary time. Because we consider the
stationary state, $J_{ky}(\phi)=\langle J_{LC}(k_{y})\rangle=\langle
J_{CR}(k_{y})\rangle$ and $\langle J_{CR}(k_{y})\rangle$ is given by

$\displaystyle\langle
J_{CR}(k_{y})\rangle=\frac{iet}{\hbar}(G_{C,0,1}(\tau=0,k_{y})-G_{C,1,0}(\tau=0,k_{y})$
$\displaystyle+e^{ik_{y}}G_{C,0,m}(\tau=0,k_{y})-e^{-ik_{y}}G_{C,m,0}(\tau=0,k_{y})).$
(13)

The Josephson current through the normal sites $J(\phi)$ is given by

$\displaystyle J(\phi)=\frac{1}{\sqrt{m^{2}+1}a}\int
J_{ky}(k_{y},\phi)dk_{y},$ (14)

because the sites at $x=0$ are aligned in the $y$ direction with a period of
$\sqrt{m^{2}+1}a$.

The Green’s function at $\tau=0$ is obtained by the Matsubara Green’s
function:

$\displaystyle
G_{C}(\tau=0,k_{y})=\frac{1}{k_{B}T}\sum_{n}G_{C}(i\epsilon_{n},k_{y}),$ (15)

where $\epsilon_{n}$ is the Matsubara frequency. To calculate the Matsubara
Green’s function in this junction system, we first calculate the functions in
the semi-infinite system: $G_{L}(i\epsilon_{n})$ for $p\leq-1$ and
$G_{R}(i\epsilon_{n})$ for $p\geq 1$. The details of the calculation of
$G_{L}(i\epsilon_{n})$ and $G_{R}(i\epsilon_{n})$ are given in the appendix.
The Green’s function at $p=0$ in the ($m$10) GB is then given by

$\displaystyle
G_{C}(i\epsilon_{n})=(i\epsilon_{n}I-H_{0,0}-H_{0,1}G_{R}(i\epsilon_{n})H_{1,0}-H_{0,-1}G_{L}(i\epsilon_{n})H_{-1,0})^{-1},$
(16)

where $H_{i,j}$ is the matrix $H_{i,j}=\langle p=i|\mathcal{H}|p=j\rangle$.

We also calculate the Josephson current for the continuum model using the
following equations for $d$-wave superconductor junctions [22, 23]:

$\displaystyle
R_{N}J=\frac{\pi\overline{R_{N}}k_{B}T}{e}\sum_{\epsilon_{n}}\int_{-\pi/2}^{\pi/2}F(\varphi,i\epsilon_{n},\phi)\sin(\phi)\sigma_{N}\cos(\varphi)d\varphi,$
(17) $\displaystyle\epsilon_{n}=2\pi k_{B}T\left(n+1/2\right),$ (18)
$\displaystyle
F(\varphi,i\epsilon_{n},\phi)=\frac{2\Delta(\varphi_{+})\Delta(\varphi_{-})}{\Omega_{n,+}\Omega_{n,-}+\epsilon^{2}_{n}+\\{1-2\sigma_{N}\sin^{2}(\phi/2)\\}\Delta(\varphi_{+})\Delta(\varphi_{-})},$
(19)
$\displaystyle\Omega_{n,\pm}=\rm{sgn}(\epsilon_{n})\sqrt{\Delta^{2}(\varphi_{\pm})+\epsilon^{2}_{n}},$
(20)
$\displaystyle\Delta(\varphi_{\pm})=\Delta\cos(2(\varphi\mp\theta/2))e^{i\phi},$
(21)
$\displaystyle\sigma_{N}=\frac{\cos^{2}(\varphi)}{\cos^{2}(\varphi)+Z^{2}},$
(22)
$\displaystyle\overline{R_{N}}^{-1}=\int_{-\pi/2}^{\pi/2}\sigma_{N}\cos(\varphi)d\varphi,$
(23)

and $Z$ is the barrier parameter.

The corresponding formula for $s$-wave superconductor junctions is [32]

$\displaystyle
R_{N}J=\frac{\pi\Delta_{0}}{2e}\frac{1}{\sqrt{1-\sigma_{N}\sin^{2}(\phi/2)}}\tanh\left(\frac{\Delta_{0}}{2k_{B}T}\sqrt{1-\sigma_{N}\sin^{2}(\phi/2)}\right)\sin\phi.$
(24)

##  3 Results

We used $\mu=-t,k_{B}T_{c}=10^{-3}t,T=0.05T_{c}$; $\Delta_{0}$ was taken to be
Bardeen–Cooper–Schrieffer (BCS)-like in the calculations given by

$\displaystyle\Delta_{0}(T)=\Delta_{0}(T=0)\tanh\left(1.74\sqrt{\frac{T_{c}-T}{T}}\right).$
(25) $\displaystyle\Delta_{0}(T=0)=1.76k_{B}T_{c}.$ (26)

The maximum value of $\Delta_{k}^{L,R}$ is 1.5 $\Delta_{0}$.

Figure 2: Current-phase relation (CPR) for a (100) grain boundary (GB). The
horizontal axis shows the phase difference and the vertical axis shows the
Josephson current. Figure 3: CPR for the (110) GB. The horizontal axis shows
the phase difference, and the vertical axis shows the Josephson current.
Figure 4: CPR for the (210) GB. The horizontal axis shows the phase
difference, and the vertical axis shows the Josephson current.

By changing $\phi$, we obtained the CPR (figures 2 to 4). Because the present
superconducting state does not break the time-reversal symmetry, the CPR is
expanded as

$J(\phi)=\sum J_{n}\sin(n\phi).$ (27)

The CPR for a (100) GB (figure 2) exhibits standard sinusoidal 0-junction
behavior with $J_{1}>0$, whereas it exhibits $\pi$-junction behavior with
$J_{1}<0$ for the (110) GB (figure 3). Because $J_{1}$ changes sign, a region
exists where $|J_{1}|$ becomes smaller than the other higher harmonics. In
this case, the CPR shows so-called $\phi$-junction behavior [22, 23], where
the higher harmonics of the CPR are enhanced (figure 4). The maximum Josephson
current $J_{max}$ is obtained from the maximum value of $J(\phi)$. We then
calculate $J_{max}$ for ($m10$) GBs for $m=1$ to 17. We also calculate
$J_{max}$ for the case of $m=\infty$ (i.e., the (100) GB). Among these
structures, the (110) and (100) GBs have a perfect crystalline square lattice
structure without defects.

Figure 5: Maximum Josephson current $J_{max}$ for $s$-wave (cross) and
$d$-wave (circle) junctions with tilting angle $\theta$.

Figure 5 shows the $J_{max}$ of ($m10$) GBs for the $s$-wave and $d$-wave
cases. First, $J_{max}$ for the $s$-wave pair potential decreases with
decreasing tilting angle $\theta$, except in the case of $\theta=0^{\circ}$,
even though the pair potential does not exhibit the phase change. This
suppression of $J_{max}$ depends on the character of the present lattice
model. In the present model, ($m10$) GBs are connected by sharing the sites at
$x=0$ with periodicity $\sqrt{m^{2}+1}a$ in the $y$ direction. The number of
conducting channels per length along the $y$ direction is then proportional to
$1/\sqrt{m^{2}+1}=1/\sqrt{\tan^{-2}(\theta/2)+1}$. The decrease in the number
of conducting channels leads to suppression of the magnitude of $J_{max}$.

Figure 6: The ratio of $J_{max}$ for a $d$-wave junction to that of an
$s$-wave junction. The dotted line shows the result predicted by Sigrist and
Rice’s theory [13]. We also plot the corresponding results obtained by the
continuum model for $Z=0$, $Z=1$, and $Z=4$. Figure 7: Momentum-resolved
Josephson current $J_{ky}$ for the (210) GB. We chose $\phi$ where $J$ gives
the maximum Josephson current; $k_{1,2,3}$ are defined in figure 8.

Next, we compare the $J_{max}$ of $s$-wave and $d$-wave pair potentials. The
$J_{max}$ for the $d$-wave cases are found to be much smaller than those for
the $s$-wave cases. To focus on the effect of the internal phase change, we
plot $J_{max}$ of $d$-wave superconductor junctions normalized by that of
$s$-wave junctions in figure 6. We also show the $\theta$ dependence of the
ratio of $J_{max}$ obtained by the phenomenological continuum model for
comparison [13]. In the case of the phenomenological model, the suppression of
$J_{max}$ originates from the cancellation of the 0-phase and $\pi$-phase
current components where the sign of $J_{1}$ is positive and negative,
respectively. Similar cancellation occurs in the present lattice model.

To confirm this calcellation, we show the $k_{y}$ dependence of the momentum-
resolved Josephson current $J_{ky}$ in figure 7, where $k_{y}$ is the parallel
momentum. We find that the sign of $J_{ky}$ changes at $k_{y}=k_{1}$ and
$k_{3}$. These momenta correspond to the positions of the nodes. That is, the
sign change of the pair potential causes the sign change of the first-order
Josephson coupling of $J_{ky}$. The cancellation of the first-order Josephson
coupling then occurs in the integration by $k_{y}$, and the resulting
Josephson current is suppressed. Thus, the magnitude of $J_{max}$ for $d$-wave
junctions is smaller than that for $s$-wave junctions.

Figure 8: Fermi surfaces and Brillouin zone for the (210) GB. The blue (red)
lines correspond to the Fermi surfaces with positive (negative) pair
potential; $k_{1},k_{2}$ are the wavelengths at nodal points of the $d$-wave
pair potential, and $k_{2}$ is the wavelength where the number of the Fermi
surface changes.

Although the reason for the sign change of the momentum-resolved current in
the lattice model is the same as that in the continuum model, the $\theta$
dependence of the Josephson current shown in figure 6 exhibits different
behaviors. For example, at lower $\theta$ values, the continuum model shows a
single minimum value, whereas the lattice model shows oscillating behavior. We
find two reasons for this difference. One reason is the folding of the Fermi
surface (FS), as shown in figure 8. This folding causes a change of the
position of the nodes and leads to overlap of the FS. The other reason is the
change of the effective barrier potential. As $m$ increases, the number of
conducting channels per unit length decreases. The corresponding effective
barrier potential then becomes pronounced. The contribution to the 0-phase is
then suppressed relative to the contribution of the $\phi$-phase, and the
current is suppressed at lower $\theta$. To observe these effects, we analyze
the physical properties of the Josephson current in the lattice model in
greater detail.

Figure 9: Position of the free-energy minimum $\phi_{m}$ plotted as a function
of tilting angle $\theta$.

First, we show the phase difference with the minimum free energy, which is
obtained by integrating the CPR with respect to $\phi$, in figure 9. A
0-junction appears in the lower-$\theta$ region, and a $\pi$-junction appears
at higher-$\theta$ values. We also observe $\phi$-junctions in certain places.
In the case of the continuum model, a $\phi$-junction only appears between the
0-junction at low $\theta$ and the $\pi$-junction at high $\theta$. In the
lattice model, a $\phi$-junction appears for the (410) GB at $\theta\simeq
28^{\circ}$ between the 0-junction and $\pi$-junction regions. However, a
$\phi$-junction also appears in the cases of the (210) and (610) GBs. This
result is distinct from that of the continuum model.

To clarify this point, we examine figure 7 in greater detail. In figure 7, a
sign change of the momentum-resolved Josephson current occurs at the position
of the nodes, as previously noted. However, the position of nodes in the
lattice model differs from that of the continuum model. In the continuum
model, a node appears at $k_{y}=k_{F}\sin(\pi/4\pm\theta/2)$ because the Fermi
surface for a symmetric-tilt GB is just given by the rotation of the $k_{x}$\-
and $k_{y}$-axes. However, in the lattice model, when
$k_{y}=k_{F}\sin(\pi/4\pm\theta/2)$ is greater than the momentum of the BZ
boundary, the FS is folded (figure 8). In the case of the (210) GB, the nodal
position is closer than in the case of the continuum model because the nodes
at $k_{y}=k_{1}$ are on the folded FSs. The area of the $\pi$-phase
($k_{3}<|k_{y}|<k_{1}$) then becomes small, and the negative contribution to
the current decreases. The $\phi$-phase then appears at $T=0.05T_{c}$, and the
resultant Josephson current is smaller than that of the continuum model.
Similarly, a $\phi$-junction also appears in the case of the (610) GB between
the 0-junction of the (510) and (710) GBs. This effect is highly sensitive to
the position of the nodal position in the BZ of the tilted structure; thus, it
strongly depends on the chemical potential $\mu$.

A further difference between the lattice model and continuum model is the
interference by the folded FSs. Considering the (210) GB as an example, two
FSs exist at $|k_{y}|>k_{2}$, as shown in figure 8, and the signs of the pair
potentials on these two FSs are opposite. Thus, when we consider the incident
and transmitted quasiparticles in the (210) GB with $|k_{y}|>k_{2}$, the first
order of the Josephson coupling is canceled because the quasiparticles’
contributions to the current have opposite signs. Actually, as shown in figure
7, the magnitude of $J_{ky}$ at $|k_{y}|>k_{2}$ is much smaller than that of
$|k_{y}|<k_{3}$, where there is no folded FS.

Figure 10: Temperature dependence of $J_{max}$ for $d$-wave junctions
normalized by $J_{max}(T=0.05T_{c})$ for an $s$-wave. Tilting angle $\theta$
for the (210), (310), and (610) GBs is $\theta\simeq$ 53∘, 37∘, and 19 ∘,
respectively. Figure 11: Spectral function
$A(k_{y},\omega+i\delta)=-\frac{1}{\pi}\sum_{i=0,1,2}\mathrm{Im}G_{C,i,i}(k_{y},\omega+i\delta)$
of the (210) GB. We chose $\phi=0,\delta=10^{-3}\Delta_{0}$.

Another feature of the proposed lattice model is that the normal conductance
depends on its tilting angle. In the present model, we do not introduce the
barrier potential or the change of the hopping integral $t$ at the interface.
In the case of the continuum model without a barrier potential ($Z=0$), a
0–$\pi$ transition occurs at $\theta=45^{\circ}$; however, the $\pi$-phase
appears at $\theta\simeq 37^{\circ}$ for the (310) GB. Thus, an effective $Z$
exists in the present model because of the decreases of the conducting channel
density. We can also confirm this result from the temperature dependence of
$J_{max}$ in figure 10. The temperature dependence of $J_{max}$ for the (210)
GB exhibits saturation behavior at low temperatures and is similar to the
previous result reported by Kulik-Omel’yanchuk [33]. However, the temperature
dependence of $J_{max}$ for the (310) and (610) GBs shows a rapid increase at
low temperatures and nonmonotonic behaviors. The rapid increase of $J_{max}$
is known to originate from the existence of the ZEABS at the interface. The
results also suggest that this ZEABS induces the $\pi$-phase. As shown in
figures 7 and 11, the 0-phase and $\pi$-phase coexist in the momentum-resolved
Josephson current and the area of the $\pi$-phase corresponds to that of the
ZEABS [22, 23]. When the effective $Z$ increases with decreasing $\theta$, the
contribution of the $\pi$-phase increases at lower temperatures. The
nonmonotonic temperature dependence due to the 0–$\pi$ transition and the
rapid increase of the Josephson current at lower temperatures then appear.

Finally, we note that a clear increase as decreasing $\theta$ in the $J_{max}$
at lower $\theta$ values is not observed in the lattice model. As shown in
figure 6, the maximum Josephson current is not suppressed for the (100) GB
($\theta=0^{\circ}$) because it is a perfect crystal structure. However, once
the tilting angle becomes nonzero, the large suppression of the $J_{max}$ is
observed. In the present calculation, the Josephson current always passes
through the sharing site at $x=0$ and its density decreases with decreasing
$\theta$. However, when $m$ becomes large, the edges of ($m10$) and
($\bar{m}10$) become nearly parallel and their interval is almost the same as
the lattice constant $a$. We then must consider the effect of direct hopping
along the $x$ direction at the interface. This hopping will make the (100) GB
effective, likely resulting in an increase of the $J_{max}$.

We here comment on the comparison between the present calculation results and
the results of actual experiments. In many experimental studies of the GBs of
HTSCs, the critical current has been reported to be severely suppressed with
increasing tilting angle [1, 2, 3, 4, 5, 6, 7, 8]. The origin of this
suppression cannot be explained using the previously reported continuum
models, as shown in figure 6. By contrast, the suppression in the small
$\theta$ region obtained in the present model suggests that the folding of the
FS and the defect formation successfully reproduced the experimental trends.
However, the recovery of $J_{max}$ with further increasing $\theta$ does not
coincide with the experimental trends. This observation indicates that
additional effects, such as diffusive scattering near the interface, occur in
actual experiments corresponding to large-$\theta$ regions. We believe the
present model is important because it provides a guide for experiments showing
that a mismatch at GBs strongly influences the Josephson current properties,
even if the mismatch angle is slight.

##  4 Summary

In this paper, we calculated the Josephson current on the symmetric [001]-tilt
GB of $d$-wave superconductor Josephson junctions with ($m10$)- and
($\bar{m}10$)-oriented surfaces on the lattice model. By changing the tilting
angle, we found a wide variety of CPRs, including $0$-, $\pi$-, and
$\phi$-junctions. In addition to the suppression of the maximum Josephson
current associated with the internal phase change of the pair potential, as
observed in the continuum models, we found that further phase interference
occurs because of the folding of the FS. The obtained maximum Josephson
current in the present lattice model is smaller than that in preexisting
theories based on the continuum model in low-angle regions. Because similar
suppression of the critical current corresponding to maximum Josephson current
at GBs has been reported in several experimental studies, the obtained results
can serve as a guide to clarify the complicated transport mechanism in GBs.

The roughness of the surface/interface is known to strongly influence the
charge transport behavior in $d$-wave superconductor junctions [34, 35, 36].
The diffusive scattering near the interface destructively influences the
contribution from ZEABSs because of the hidden odd-frequency odd-parity spin-
singlet pairing [37, 38, 39, 40, 41]. Clarifying how the grain angle
dependence of the critical current due to the interference of the phase from
the folding of the BZ predicted in the present paper is influenced by the
diffusive scattering would be an interesting topic for a future study.

In the present work, we have focused on the GB effect in $d$-wave
superconductors. Studies of iron-based superconductors in which $s_{\pm}$
pairing is a promising symmetry have also shown interesting results [42, 43].
Because there are several theoretical works related to surface Andreev bound
states, quasiparticle tunneling, and Josephson effects in $s_{\pm}$-wave
superconductors [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55], it is
interesting to calculate the Josephson current in GB $s_{\pm}$ superconductor
junctions.

This work was supported by Scientific Research (A) (KAKENHIGrant No.
JP20H00131), Scientific Research (B) (KAKENHIGrant No. JP20H01857), and
Scientific Research (Early-Career Scientists) (KAKENHIGrant No. 21K13854).

##  Appendix A Recursive Green’s function method for ($m10$) surface

In this appendix, we derive the Green’s function in the junction of ($m10$)
GB. To calculate this, we first calculate the surface Green’s function at
($m10$) and ($\bar{m}10$) surfaces. For that purpose, we use the recursive
Green’s function method. As given in the main text, the matrix elements for
($m10$) surface is given by

$\displaystyle\langle p,q|\mathcal{H}|p+1,q\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&0&0&\cdots&0\\\ t_{x}&0&0&\cdots&0\\\
0&t_{x}&0&\ddots&0\\\ \vdots&\ddots&\ddots&\ddots&0\\\
0&0&0&t_{x}&0\end{pmatrix},$ (28) $\displaystyle\langle
p,q|\mathcal{H}|p,q\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}h_{0}&t_{y}&0&\cdots&0\\\
t_{y}^{\dagger}&h_{0}&t_{y}&\cdots&0\\\ 0&t_{y}^{\dagger}&h_{0}&\ddots&0\\\
\vdots&\ddots&\ddots&\ddots&t_{y}\\\
0&0&0&t_{y}^{\dagger}&h_{0}\end{pmatrix},$ (29) $\displaystyle\langle
p,q|\mathcal{H}|p,q+1\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&\cdots&0&t_{x}\\\ 0&\ddots&0&0\\\
\vdots&\ddots&\ddots&\vdots\\\ 0&\cdots&0&0\end{pmatrix},$ (30)
$\displaystyle\langle p,q|\mathcal{H}|p+1,q-1\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&\cdots&\cdots&0\\\
\vdots&\ddots&\ddots&\vdots\\\ 0&0&\ddots&\vdots\\\
t_{y}&0&\cdots&0\end{pmatrix},$ (31)

where

$\displaystyle h_{0}$ $\displaystyle=$ $\displaystyle\begin{pmatrix}-\mu&0\\\
0&\mu\end{pmatrix},$ (32) $\displaystyle t_{x}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-t&\Delta(p)\\\ \Delta^{*}(p)&t\end{pmatrix},$
(33) $\displaystyle t_{y}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-t&-\Delta(p)\\\ -\Delta^{*}(p)&t\end{pmatrix}.$
(34)

By the Fourier transformation in $q$, we can diagonalize $\mathcal{H}$ in $q$,

$\displaystyle\langle p,k_{y}|\mathcal{H}|p+1,k_{y}\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&0&0&\cdots&0\\\ t_{x}&0&0&\cdots&0\\\
0&t_{x}&0&\ddots&0\\\ \vdots&\ddots&\ddots&\ddots&0\\\
t_{y}e^{iky}&0&0&t_{x}&0\end{pmatrix}$ (35) $\displaystyle\equiv$
$\displaystyle T_{x}(k_{y}),$ $\displaystyle\langle
p,k_{y}|\mathcal{H}|p,k_{y}\rangle$ $\displaystyle=$
$\displaystyle\begin{pmatrix}h_{0}&t_{y}&0&\cdots&t_{x}e^{-ik_{y}}\\\
t_{y}^{\dagger}&h_{0}&t_{y}&\cdots&0\\\ 0&t_{y}^{\dagger}&h_{0}&\ddots&0\\\
\vdots&\ddots&\ddots&\ddots&t_{y}\\\
t_{x}^{\dagger}e^{ik_{y}}&0&0&t_{y}^{\dagger}&h_{0}\end{pmatrix}$ (36)
$\displaystyle\equiv$ $\displaystyle H_{0}(k_{y}),$

where $T_{x}$ and $H_{0}$ are $2(m+1)\times 2(m+1)$ matrices which include
sublattice and particle-hole degree of freedom.

Here, we suppose that the surface Green’s function at $p=p_{0}+1$ in the
$n$-layer system from $p=p_{0}+1$ to $p=p_{0}+n$ is known. Then, the surface
Green’s function at $p=p_{0}$ in the $(n+1)$-layer system from $p=p_{0}$ to
$p=p_{0}+n$ is given by

$\displaystyle
G^{(n+1)}_{s}(i\epsilon_{n},k_{y})=\left(i\epsilon_{n}I-H_{0}-T_{x}G^{(n)}_{s}(i\epsilon_{n},k_{y})T_{x}^{\dagger}\right)^{-1},$
(37)

where $G^{(n+1)}_{s}(i\epsilon_{n},k_{y})$ stands for the surface Green’s
function with Matsubara Frequency $\epsilon_{n}$ at $p=p_{0}$ ($p=p_{0}+1$) in
the $(n+1)$-layer ($n$-layer) system. Since it is trivial to obtain the
Green’s function of the 1-layer system by
$G^{(1)}_{s}(i\epsilon_{n},k_{y})=(i\epsilon_{n}I-H_{0})^{-1}$, we can
calculate the surface Green’s function in the system with any number of layers
by the recursion relation in (37). By repeating this recursive equation,
surface Green’s function converges as
$G^{(n+1)}_{s}(i\epsilon_{n},k_{y})=G^{(n)}_{s}(i\epsilon_{n},k_{y})$.
However, to get the convergence, we have to consider the system much larger
than the coherence length. Thus, we can express the recursive equation in
terms of a Möbius transformation. As we will explain later, $T_{x}$ must have
the inverse matrix to do this. As seen from (35), $T_{x}$ is a lower
triangular matrix without any diagonal components. Then, $T_{x}$ does not have
the inverse matrix. Thus, we express (37) in the another form. The matrix form
of the total Hamiltonian is given by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\ddots&\ddots&\ddots&&\\\
\ddots&H_{0}&T_{x}&0&\\\ \ddots&T_{x}^{\dagger}&H_{0}&T_{x}&\ddots\\\
&0&T_{x}^{\dagger}&H_{0}&\ddots\\\ &&\ddots&\ddots&\ddots\\\ \end{pmatrix},$
(38)

where $H_{0}$ and $T_{x}$ are $2(m+1)\times 2(m+1)$ matrices. Consider the $m$
contiguous layers of this system. Then, we can express this Hamiltonian by the
$2m\times 2m$ matrices as shown in figure 12.

Figure 12: $2(m+1)\times 2(m+1)$ matrix $H_{0}$ and $T_{0}$ are represented by
$2m\times 2m$ matrix $\check{H}_{1,2,...,m+1},\check{T}_{12,23,...,(m+1,1)}$.

By using these matrices, we make the recursion relation between surface
Green’s functions for the $n$-layer system and $(n+m)$-layer system,

$\displaystyle\check{G}_{s}^{n+m}=X_{\bullet}\check{G}_{s}^{n},$ (39)

where $\check{G}_{s}$ is the submatrix of $G_{s}$ including the sublattice
$o=1$ to $m$ and therefore $\check{G}_{s}$ is $2m\times 2m$ matrix. Here, ∙
denotes the Möbius transformation defined by

$\displaystyle\begin{pmatrix}a&b\\\
c&d\end{pmatrix}_{\bullet}z\equiv(az+b)(cz+d)^{-1}.$ (40)

$X$ is given by

$\displaystyle X=$ $\displaystyle\prod_{i=1}^{m+1}X_{i},$ (41) $\displaystyle
X_{i}=$ $\displaystyle\begin{pmatrix}0&\check{T}_{m+i,i}^{-1}\\\
\check{T}_{m+i,i}^{\dagger}&(i\epsilon_{n}-\check{H}_{i})\check{T}_{m+i,i}^{-1}\end{pmatrix}$
(42)

where the subscript of $\check{H}$ and $\check{T}$ is integers in modulo
$m+1$. These equations rewrite the process of adding $m$ layers containing
$m+1$ sublattices into a process of adding $m$ sublattices $m+1$ times. Since
we can express the recursion relation in terms of a Möbius transformation in
(39), the problem to calculate the surface Green’s function for semiinfinite
layer is reduced into the eigenvalue problem of the matrix $X$.

$X$ is diagonalized as

$\displaystyle X=Q\begin{pmatrix}\lambda_{1}&&&{\huge O}\\\ &\lambda_{2}&&\\\
&&\ddots&\\\ {\huge O}&&&\lambda_{2m}\\\ \end{pmatrix}Q^{-1}$ (43)

$\lambda_{1,2,\ldots,2m}$ are the eigenvalue of the matrix $X$.

$\displaystyle|\lambda_{1}|<|\lambda_{2}|<\ldots<|\lambda_{2m}|$ (44)
$\displaystyle\check{G}_{s}^{\infty}=Q_{\bullet}I$ (45)

We can calculate $G_{s}^{\infty}$ by adding a single layer from
$\check{G}_{s}^{\infty}$.

## References

## References

  * [1] Held R, Schneider C W, Mannhart J, Allard L F, More K L and Goyal A 2009 Phys. Rev. B 79 014515
  * [2] Dimos D, Chaudhari P, Mannhart J and LeGoues F K 1988 Phys. Rev. Lett. 61(2) 219–222
  * [3] Hilgenkamp H and Manhart J 2002 Rev. Mov. Phys 74 485
  * [4] Ivanov Z G, Stepantsov E A, Claeson T, Wenger F, Lin S Y, Khare N and Chaudhari P 1998 Phys. Rev. B 57(1) 602–607
  * [5] Durrell J H, Eom C B, Gurevich A, Hellstrom E E, Tarantini C, Yamamoto A and Larbalestier D C 2011 Reports on Progress in Physics 74 124511
  * [6] Katase T, Ishimaru Y, Tsukamoto A, Hiramatsu H, Kamiya T, Tanabe K and Hosono H 2011 Nature Communications 2 1–6
  * [7] Si W, Zhang C, Shi X, Ozaki T, Jaroszynski J and Li Q 2015 Applied Physics Letters 106 032602
  * [8] Iida K, Omura T, Matsumoto T, Hatano T and Ikuta H 2019 Superconductor Science and Technology 32 074003
  * [9] Browning N D, Buban J P, Nellist P D, Norton D P, Chisholm M F and Pennycook S J 1998 Physica C 294 183–193
  * [10] JMannhart, Bielefeldt H, Goetz B, Hilgenkamp H, Schmehl A, Schneider C W and Schulz R R 2000 Philos. Mag. B 80 827–834
  * [11] Schmehl A, Goetz B, Schulz R R, Schneider C W, Bielefeldt H, Hilgenkamp H and Mannhart J 1999 Europhysics Letters 47 110
  * [12] Hammerl G, Schmehl A, Schulz R R, Götz B, Bielefeldt H, Schneider C W, Hilgenkamp H and Mannhart J 2000 Nature 407 162–164
  * [13] Sigrist M and M Rice T 1992 Journal of the Physical Society of Japan 61 4283–4286
  * [14] Wollman D A, Van Harlingen D J, Lee W C, Ginsberg D M and Leggett A J 1993 Phys. Rev. Lett. 71(13) 2134–2137
  * [15] Van Harlingen D J 1995 Rev. Mod. Phys. 67(2) 515–535
  * [16] Tsuei C C, Kirtley J R, Chi C C, Yu-Jahnes L S, Gupta A, Shaw T, Sun J Z and Ketchen M B 1994 Phys. Rev. Lett. 73(4) 593
  * [17] Tsuei C C and Kirtley J R 2000 Rev. Mod. Phys. 72(4) 969
  * [18] Hu C R 1994 Phys. Rev. Lett. 72(10) 1526
  * [19] Tanaka Y and Kashiwaya S 1995 Phys. Rev. Lett. 74(17) 3451
  * [20] Kashiwaya S and Tanaka Y 2000 Rep. Prog. Phys. 63 1641
  * [21] Löfwander T, Shumeiko V S and Wendin G 2001 Supercond. Sci. Technol. 14 R53
  * [22] Tanaka Y and Kashiwaya S 1996 Phys. Rev. B 53(18) R11957
  * [23] Tanaka Y and Kashiwaya S 1997 Phys. Rev. B 56(2) 892
  * [24] Barash Y S, Burkhardt H and Rainer D 1996 Phys. Rev. Lett. 77(19) 4070
  * [25] Testa G, Sarnelli E, Monaco A, Esposito E, Ejrnaes M, Kang D J, Mennema S H, Tarte E J and Blamire M G 2005 Phys. Rev. B 71(13) 134520
  * [26] Tafuri F and Kirtley J R 2005 Reports on Progress in Physics 68 2573
  * [27] Il’ichev E, Zakosarenko V, IJsselsteijn R P J, Schultze V, Meyer H G, Hoenig H E, Hilgenkamp H and Mannhart J 1998 Phys. Rev. Lett. 81(4) 894–897
  * [28] Il’ichev E, Zakosarenko V, IJsselsteijn R P J, Hoenig H E, Schultze V, Meyer H G, Grajcar M and Hlubina R 1999 Phys. Rev. B 60(5) 3096–3099
  * [29] Tanuma Y, Tanaka Y, Yamashiro M and Kashiwaya S 1998 Phys. Rev. B 57 7997
  * [30] Tanuma Y, Tanaka Y, Ogata M and Kashiwaya S 1999 Phys. Rev. B 60 9817
  * [31] Shirai S, Tsuchiura H, Asano Y, Tanaka Y, Inoue J i, Tanuma Y and Kashiwaya S 2003 Journal of the Physical Society of Japan 72 2299–2307
  * [32] Furusaki A and Tsukada M 1991 Solid State Commun. 78 299
  * [33] Kulik I O and Omelyanchuk A N 1978 Sov. J. Low Temp. Phys. 4 142
  * [34] Yamada K, Nagato Y, Higashitani S and Nagai K 1996 J. Phys. Soc. Jpn 65 1540
  * [35] Golubov A A and Kupriyanov M Y 1998 Journal of Experimental and Theoretical Physics Letters 67 501–507 ISSN 1090-6487
  * [36] Golubov A A and Kupriyanov M Y 1999 Journal of Experimental and Theoretical Physics Letters 69 262–267 ISSN 1090-6487
  * [37] Tanaka Y, Nazarov Y V and Kashiwaya S 2003 Phys. Rev. Lett. 90 167003
  * [38] Tanaka Y, Nazarov Y V, Golubov A A and Kashiwaya S 2004 Phys. Rev. B 69 144519
  * [39] Tanaka Y and Golubov A A 2007 Phys. Rev. Lett. 98 037003
  * [40] Tanaka Y, Golubov A A, Kashiwaya S and Ueda M 2007 Phys. Rev. Lett. 99 037005
  * [41] Tanaka Y, Tanuma Y and Golubov A A 2007 Phys. Rev. B 76 054522
  * [42] Mazin I I, Singh D J, Johannes M D and Du M H 2008 Phys. Rev. Lett. 101(5) 057003
  * [43] Kuroki K, Onari S, Arita R, Usui H, Tanaka Y, Kontani H and Aoki H 2008 Phys. Rev. Lett. 101(8) 087004
  * [44] Golubov A A, Brinkman A, Tanaka Y, Mazin I I and Dolgov O V 2009 Phys. Rev. Lett. 103(7) 077003
  * [45] Onari S and Tanaka Y 2009 Phys. Rev. B 79(17) 174526
  * [46] Ghaemi P, Wang F and Vishwanath A 2009 Phys. Rev. Lett. 102(15) 157002
  * [47] Nagai Y and Hayashi N 2009 Phys. Rev. B 79(22) 224508
  * [48] Linder J, Sperstad I B and Sudbø A 2009 Phys. Rev. B 80(2) 020503
  * [49] Ota Y, Nakai N, Nakamura H, Machida M, Inotani D, Ohashi Y, Koyama T and Matsumoto H 2010 Phys. Rev. B 81(21) 214511
  * [50] Huang W M and Lin H H 2010 Phys. Rev. B 81(5) 052504
  * [51] Bobkov A M and Bobkova I V 2011 Phys. Rev. B 84(13) 134527
  * [52] Seidel P 2011 Superconductor Science and Technology 24 043001
  * [53] V Burmistrova A, A Devyatov I, A Golubov A, Yada K and Tanaka Y 2013 Journal of the Physical Society of Japan 82 034716
  * [54] Burmistrova A V, Devyatov I A, Golubov A A, Yada K and Tanaka Y 2013 Superconductor Science and Technology 27 015010
  * [55] Burmistrova A V, Devyatov I A, Golubov A A, Yada K, Tanaka Y, Tortello M, Gonnelli R S, Stepanov V A, Ding X, Wen H H and Greene L H 2015 Phys. Rev. B 91(21) 214501