††institutetext: Laboratoire de physique math matique et subatomique.
Fr res Mentouri university,
Costantine 1, Algeria

# A new anomaly free flipped 341 model

Meriem Djouala  Noureddine Mebarki and Habib Aissaoui
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A new flipped $SU(3)_{C}\otimes SU(4)_{L}\otimes U(1)_{X}$ model without
exotic electric charges is proposed. All the quarks families are arranged in
the same representation while leptons generations are in different
representations leading to a tree level FCNC. Moreover, It is shown that the
cancellation of the triangle chiral anomalies requires new additional leptons
10-plet and quadruplet. All fermions masses have been also discussed.
Furthermore, using the most recent experimental data of the branching ratios
of $\mu\longrightarrow ee\overline{e}$ and $\mu\longrightarrow e\gamma$ rare
decay modes, stringent bounds on the heavy neutral bosons masses and the muon-
electron mixing matrix element are obtained.

###### Keywords:

The flipped 341 model, fermions masses, FCNC

##  1 Introduction

Neutrinos oscillation phenomenon, dark matter, replication of quarks families,
charge quantization and many other fundamental questions reveal that the
Standard Model is an effective gauge field theory. Thus, going beyond the
Standard Model (BSM) becomes mandatory to explain these outstanding unsolved
problems. In this paper, we are interested in a specific BSM model based on
the Lie gauge group $SU(3)_{C}\otimes SU(4)_{L}\otimes U(1)_{X}$ (denoted by
341 for short) ref1 . Many interest has been devoted to this kind of models
especially those related to the LHC physics ref1 ; ref7 . The most attractive
feature of those models is the explanation of the family replication coming
from the triangle gauge anomaly cancellation which together with QCD
asymptotic freedom, requires that only three fermion families might be present
in the spectrum, as in the 331 models. Moreover, it reproduces all the SM
success including the global fit to electroweak precision data. In general,
requiring the gauge anomaly cancellation and avoiding the existence of
fractional leptons charges and since $SU(2)_{L}\otimes U(1)_{Y}$ of the
Standard Model (SM) can be embedded into $SU(4)_{L}\otimes U(1)_{X}$ in more
than one way, these models have to depend on two parameters $\beta$ and
$\gamma$ ref5 . In fact, in order to assign the fermionic content of the
model, it is useful to recall that the electric charge operator is written as
a linear combination of the diagonal generators $T_{3}$, $T_{8}$, $T_{15}$ and
X of the gauge group $SU(4)_{L}\otimes U(1)_{X}$. This allows us to write the
electric charge assignment for the fields in the fundamental representation of
$SU(4)_{L}$ as:

$\dfrac{Q}{e}=T_{3}+\beta T_{8}+\gamma T_{15}+X$ (1)

It is worth to mention that the 341 models were also classified according to
the scalars sectors ref5 ; ref7 . Moreover, the self consistency of any BSM
model is to be free from gauge anomalies. For the case of our interest and as
it was pointed out before, the cancellation of the $[SU(4)_{L}]^{3}$ gauge
anomaly in all versions of the 341 models is fulfilled if the three quarks
generations belong to different $SU(4)_{L}$ representations: two with a left
handed chirality $Q_{iL}$ (i=1,2) lie in the fundamental representation 4,
whereas , the third one $Q_{3L}$ together with the three leptons generations
$\psi_{lL}(l=e,\mu,\tau)$ have to transform under the conjugate fundamental
representation $\bar{4}$ (or vice versa). The right handed quarks transform as
singlets. For the scalar content with desirable properties concerning the
symmetry breakdown and fermion masses including neutrinos, one can take four
quartets $\phi_{k}$ (k=1,4) which are in the conjugate representation
$\bar{4}$. Table 1 shows particles content in all possible 341 models
parameterized by $\beta$ and $\gamma$, and the corresponding 331 sub-
representations. The 321 sub-representations embedded in the fundamental and
its conjugate representations of SU(3) are given by:

$\displaystyle\psi_{l}$ $\displaystyle\supset$
$\displaystyle(1,2,\frac{-1}{2})\oplus(1,1,\frac{-1}{2}-\frac{\sqrt{3}\beta}{2})\oplus(1,1,\frac{-2\gamma}{\sqrt{6}}-\frac{\beta}{2\sqrt{3}}-\frac{1}{2}),$
(2) $\displaystyle Q_{i}$ $\displaystyle\supset$
$\displaystyle(3,2,\frac{1}{6})\oplus(3,1,\frac{1}{6}+\frac{\sqrt{3}\beta}{2})\oplus(3,1,\frac{1}{6}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}}),$
(3) $\displaystyle Q_{3}$ $\displaystyle\supset$
$\displaystyle(3,2,\frac{1}{6})\oplus(3,1,\frac{1}{6}-\frac{\sqrt{3}\beta}{2})\oplus(3,1,\frac{1}{6}-\frac{\beta}{2\sqrt{3}}-\frac{2\gamma}{\sqrt{6}}),$
(4) $\displaystyle\phi_{1}$ $\displaystyle\supset$
$\displaystyle(1,2,\frac{-1}{2})\oplus(1,1,\frac{-1}{2}+\frac{\sqrt{3}\beta}{2})\oplus(1,1,\frac{-\beta}{\sqrt{3}}+\frac{2\gamma}{\sqrt{6}})$
(5) $\displaystyle\phi_{2}$ $\displaystyle\supset$
$\displaystyle(1,2,\frac{-\sqrt{3}\beta}{2})\oplus(1,1,0)\oplus(1,1,\frac{-1}{2}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}})$
(6) $\displaystyle\phi_{3}$ $\displaystyle\supset$
$\displaystyle(1,2,\frac{1}{2})\oplus(1,1,\frac{1}{2}+\frac{\sqrt{3}\beta}{2})\oplus(1,1,\frac{1}{2}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}})$
(7) $\displaystyle\phi_{4}$ $\displaystyle\supset$
$\displaystyle(1,2,\frac{1}{2})\oplus(1,1,\frac{1}{2}+\frac{\sqrt{3}\beta}{2})\oplus(1,1,0).$
(8)

Name | 341 representation | 331 representation | Components | F
---|---|---|---|---
$\psi_{l}$ |  (1,4,$\frac{-1}{2}-\frac{\beta}{2\sqrt{3}}-\frac{\gamma}{2\sqrt{6}})$ |  (1,3,$\frac{-1}{2}-\frac{\beta}{2\sqrt{3}}$)$\oplus$(1,1,$\frac{-2\gamma}{\sqrt{6}}-\frac{\beta}{2\sqrt{3}}-\frac{1}{2}$) |  ($\nu_{l},l,E_{l}^{q_{1}},N_{l}^{q_{2}}$) | 3
$l^{c}$ | (1,1,1) | (1,1,1) | $l^{c}$ | 3
$l_{X}^{c}$ |  (1,1,$\frac{1}{2}+\frac{\sqrt{3}}{2}\beta$) |  (1,1,$\frac{1}{2}+\frac{\sqrt{3}}{2}\beta$) | $l_{X}^{c}$ | 3
$Q_{i}(i=1,2)$ |  (3,$\bar{4}$,$\frac{1}{6}+\frac{\beta}{2\sqrt{3}}+\frac{\gamma}{2\sqrt{6}})$ |  (3,$\bar{3}$,$\frac{1}{6}+\frac{\beta}{2\sqrt{3}}$)$\oplus$(1,1,$\frac{1}{6}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}}$) | $(u,d,D_{i}^{q_{1}},U_{i}^{q_{2}})$ | 2
$Q_{3}$ |  (3,4,$\frac{1}{6}-\frac{\beta}{2\sqrt{3}}-\frac{\gamma}{2\sqrt{6}})$ |  (3,3,$\frac{1}{6}-\frac{\beta}{2\sqrt{3}}$)$\oplus$(1,1,$\frac{1}{6}-\frac{\beta}{2\sqrt{3}}-\frac{2\gamma}{\sqrt{6}}$) | $(d,u,U_{3}^{q_{1}},D_{3}^{q_{2}})$ | 1
$u^{c}$ |  ($\bar{3}$,1,$\frac{-2}{3}$) |  ($\bar{3}$,1,$\frac{-2}{3}$) | $u^{c}$ | 3
$d^{c}$ |  ($\bar{3}$,1,$\frac{1}{3}$) |  ($\bar{3}$,1,$\frac{1}{3}$) | $d^{c}$ | 3
$U^{c}_{1,2}$ |  ($\bar{3}$,1,$\frac{-1}{6}-\frac{\sqrt{3}\beta}{2}$) |  ($\overline{3}$,1,$\frac{-1}{6}-\frac{\sqrt{3}\beta}{2}$) | $U^{c}_{1,2}$ | 2
$U^{c}_{3}$ |  ($\bar{3}$,1,$\frac{-1}{6}+\frac{\sqrt{3}\beta}{2}$) |  ($\overline{3}$,1,$\frac{-1}{6}+\frac{\sqrt{3}\beta}{2}$) | $U^{c}_{3}$ | 1
$D^{c}_{1,2}$ |  ($\bar{3}$,1,$\frac{-1}{6}-\frac{\beta}{2\sqrt{3}}-\frac{2\gamma}{\sqrt{6}}$) |  ($\bar{3}$,1,$\frac{-1}{6}-\frac{\beta}{2\sqrt{3}}-\frac{2\gamma}{\sqrt{6}}$) | $D^{c}_{1,2}$ | 2
$D^{c}_{3}$ |  ($\bar{3}$,1,$\frac{-1}{6}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}}$) |  ($\bar{3}$,1,$\frac{-1}{6}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}}$) | $D^{c}_{1,2}$ | 1
$\phi_{1}$ |  (1,$\bar{4}$,$\frac{-\beta}{\sqrt{3}}+\frac{\gamma}{2\sqrt{6}})$ |  (1,$\bar{3}$,-$\frac{1}{2}+\frac{\beta}{2\sqrt{3}}$)$\oplus$(1,1,$-\frac{\beta}{\sqrt{3}}+\frac{2\gamma}{\sqrt{6}}$) |  ($\phi_{1}^{q_{1}}$,$\phi_{1}^{q_{2}}$,$\phi_{1}^{q_{3}}$,$\phi_{1}^{q_{4}})$ | 1
$\phi_{2}$ |  (1,$\bar{4}$,$\frac{-1}{2}+\frac{\beta}{2\sqrt{3}}+\frac{\gamma}{2\sqrt{6}})$ |  (1,$\bar{3}$,$-\frac{\beta}{\sqrt{3}}$)$\oplus$(1,1,$-\frac{1}{2}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}}$) |  ($\phi_{2}^{q_{1}}$,$\phi_{2}^{q_{2}}$,$\phi_{2}^{q_{3}}$,$\phi_{2}^{q_{4}})$ | 1
$\phi_{3}$ |  (1,$\bar{4}$,$\frac{1}{2}+\frac{\beta}{2\sqrt{3}}+\frac{\gamma}{2\sqrt{6}})$ |  (1,$\bar{3}$,$\frac{1}{2}+\frac{\beta}{2\sqrt{3}}$)$\oplus$(1,1,$\frac{1}{2}+\frac{\beta}{2\sqrt{3}}+\frac{2\gamma}{\sqrt{6}}$) |  ($\phi_{3}^{q_{1}}$,$\phi_{3}^{q_{2}}$,$\phi_{3}^{q_{3}}$,$\phi_{3}^{q_{4}})$ | 1
$\phi_{4}$ |  (1,$\bar{4}$,$\frac{-3\gamma}{2\sqrt{6}})$ |  (1,$\bar{3}$,$\frac{1}{2}$+$\frac{\beta}{2\sqrt{3}}$)$\oplus$(1,1,0) |  ($\phi_{4}^{q_{1}}$,$\phi_{4}^{q_{2}}$,$\phi_{4}^{q_{3}}$,$\phi_{4}^{q_{4}})$ | 1
Table 1:  Particles content of the 341 models for generic $\beta$ and $\gamma$
parameters where F represents the number of flavors.

It turns out that this quarks and leptons replication is not the only way to
have a model free from the $SU(4)_{L}\otimes U(1)_{X}$ gauge anomalies. In
this work and following ref ref8 , we build a new unique gauge anomaly free
model without exotic electric charges baptized flipped 341 as an extension of
where the previous scheme of construction is reversed that is all the quarks
generations transform under the same representation while leptons are not.
Thus, a flavor changing neutral current (FCNC)is expected at the tree level in
the lepton sector through the exchange of new neutral gauge bosons
$Z^{\prime}$ and $Z^{\prime\prime}$ of the model like in the rare leptonic
decays of the form $l_{i}\longrightarrow l_{i}l_{k}\overline{l}_{k}$ and
$l_{i}\longrightarrow l_{k}\gamma$ $(i\neq k)$.
This paper is organized as follows: in section 2, we present the basic
ingredients of the construction and features of the new flipped 341 model. In
section 3, the anomaly cancellation is discussed. In section 4, fermions
masses are emphasized at the tree level. In section 5, the FCNC and mixing in
the leptonic sector are explored via the study of the rare processes
$\mu\longrightarrow ee\overline{e}$ and $\mu\longrightarrow e\gamma$ and
stringent inequalities and bounds on the $Z^{\prime}$ and $Z^{\prime\prime}$
masses are obtained. Finally, in section 6, we draw our conclusions.

##  2 A flipped 341 model

In this model, the three quarks generations transform under the same conjugate
fundamental representations while the leptons generations lie in different
representations. It is important to mention that contrary to the ordinary 341
models, a fourth lepton generation is required in order to cancel the gauge
anomaly $[SU(4)]^{3}_{L}$ where the number of fermions quadruplets must be
equal to the number of fermions anti quadruplets. Table 2 summarizes the
particles content in this model where there are no fractionally electric
charged different from $\mp\frac{2}{3}$ and $\mp\frac{1}{3}$ for exotic quarks
and there are no integer electric charges different from 0 and $\mp 1$ for
leptons and gauge bosons (model without exotic electric charges). Moreover,
this model is unique, in the sense that the decomposition of the lepton
10-plet $\psi_{e}$ contains as subgroups
(1,6,X)$\oplus$(1,3,X’)$\oplus$(1,1,Q). To get a real (1,6,X) representation,
the value of $\beta$ should be equal to $\frac{1}{\sqrt{3}}$ ref8 . The
restriction to non exotic electric charges only allows for eight different
anomaly free 341 models ref12  and for a giving value of $\beta$ one has tow
models with $\gamma=\frac{1}{\sqrt{6}}$ and $\gamma=\frac{-2}{\sqrt{6}}$. In
the flipped version of the first model one has fractionally charged leptons
which are forbidden.111Fractionally charged leptons, produced abundantly in
the early universe, would lead to an unacceptable cosmology, because there are
no decay modes for the lightest of these exotic states ref8 .. Therefore, the
only acceptable flipped 341 is the one with $\beta=\frac{1}{\sqrt{3}}$,
$\gamma=\frac{-2}{\sqrt{6}}$.
The 10-plet $\psi_{e}$ and 15-plet S contain $\Sigma^{\mp}$, $\Sigma^{0}$,
$\triangle^{\mp}$ and $\triangle^{0}$ which stands for triplets sub-
representations, $(\nu_{e},e)$, $(H^{+}_{S},H^{{}^{\prime}0}_{S})$ and
$(H^{{}^{\prime\prime}0}_{S},H^{-}_{2S})$ form $SU(2)_{L}$ doublets,
($\beta^{+}$, $\beta^{0}$, $N_{e}^{0}$,
$(H^{0}_{S}$,$H^{+}_{1S}$,$H_{1S}^{{}^{\prime}+})$ and
$(H^{{}^{\prime}0}_{1S}$,$H^{{}^{\prime}-}_{2S}$,$H_{S}^{{}^{\prime}-})$
compose the triplet denoted by (1,3,$\frac{1}{3}$),
(1,$\overline{3}$,$\frac{2}{3}$) and (1,3,$\frac{-2}{3}$) respectively, while,
$E_{e}^{-}$, $\sigma^{+}$, $\sigma_{S}^{0}$ and $H^{0}_{1S}$ are singlets.
In our flipped 341 model and in order to generate masses for particles one has
to have four scalar fields $\phi_{1}$, $\phi_{2}$, $\phi_{3}$, $\phi_{4}$ and
a new scalar 15-plet S which transforms under the adjoint representation 15 of
the group SU(4).

Name | 341 representation | 331 representation | components | F
---|---|---|---|---
$\psi_{e}$ | (1,10,0) |  (1,6,$\frac{-1}{3}$)+(1,3,$\frac{1}{3}$)+(1,1,1) | $\left(\begin{array}[]{cccc}\Sigma^{+}&\frac{\Sigma^{0}}{\sqrt{2}}&\frac{\nu_{e}}{\sqrt{2}}&\frac{\beta^{+}}{\sqrt{2}}\\\ \frac{\Sigma^{0}}{\sqrt{2}}&\Sigma^{-}&\frac{e^{-}}{\sqrt{2}}&\frac{\beta^{0}}{\sqrt{2}}\\\ \frac{\nu_{e}}{\sqrt{2}}&\frac{e^{-}}{\sqrt{2}}&E_{e}^{-}&\frac{N_{e}^{0}}{\sqrt{2}}\\\ \frac{\beta^{+}}{\sqrt{2}}&\frac{\beta^{0}}{\sqrt{2}}&\frac{N_{e}^{0}}{\sqrt{2}}&\sigma^{+}\\\ \end{array}\right)$ | 1
$\psi_{\alpha}(\alpha=\mu,\tau)$ |  (1,4,$\frac{-1}{2}$) |  (1,3,$\frac{-2}{3}$)+(1,1,0) |  ($\nu_{\alpha},l_{\alpha},E^{-}_{\alpha},N_{\alpha}^{0}$) | 2
$\widetilde{\psi}$ |  (1,$\overline{4}$,$\frac{-1}{2}$) |  (1,3,$\frac{-1}{3}$)+(1,1,-1) |  ($\widetilde{e}^{-},\widetilde{\nu},\widetilde{N}^{0},\widetilde{E}^{-}$) | 1
$l_{\alpha}^{c}$ | (1,1,1) | (1,1,1) | $l^{c}_{\alpha}$ | 6
$Q_{\alpha}$ |  (3,$\overline{4}$,$\frac{1}{6})$ |  (3,$\bar{3}$,$\frac{1}{3}$)+(3,1,$\frac{-1}{3}$) | $(d_{\alpha},u_{\alpha},U_{\alpha},D_{\alpha})$ | 3
$u^{c}$ |  ($\overline{3}$,1,$\frac{-2}{3}$) |  ($\overline{3}$,1,$\frac{-2}{3}$) | $u^{c}$ | 6
$d^{c}$ |  ($\overline{3}$,1,$\frac{1}{3}$) |  ($\overline{3}$,1,$\frac{1}{3}$) | $d^{c}$ | 6
$\phi_{1}$ |  (1,$\overline{4}$,$\frac{-1}{2})$ |  (1,$\bar{3}$, $\frac{-1}{3}$)+(1,1,-1) |  ($\phi_{1}^{-}$,$\phi_{1}^{0}$,$\phi_{1}^{{}^{\prime}0}$,$\phi_{1}^{{}^{\prime}-})$ | 1
$\phi_{2}$ |  (1,$\overline{4}$,$\frac{-1}{2})$ |  (1,$\bar{3}$,$\frac{-1}{3}$)+(1,1,-1) |  ($\phi_{2}^{-}$,$\phi_{2}^{0}$,$\phi_{2}^{{}^{\prime}0}$,$\phi_{2}^{{}^{\prime}-})$ | 1
$\phi_{3}$ |  (1,$\overline{4}$,$\frac{1}{2})$ |  (1,$\bar{3}$,$\frac{2}{3}$)+(1,1,0) |  ($\phi_{3}^{0}$,$\phi_{3}^{+}$,$\phi_{3}^{{}^{\prime}+}$,$\phi_{3}^{{}^{\prime}0})$ | 1
$\phi_{4}$ |  (1,$\overline{4}$,$\frac{1}{2})$ |  (1,$\bar{3}$,$\frac{2}{3}$)+(1,1,0) |  ($\phi_{4}^{0}$,$\phi_{4}^{+}$,$\phi_{4}^{{}^{\prime}+}$,$\phi_{4}^{{}^{\prime}0})$ | 1
S | (1,15,0) |  (1,8,0)+(1,3,$\frac{-2}{3}$)+(1,$\bar{3}$,$\frac{2}{3}$)+(1,1,0) | $\left(\begin{array}[]{cccc}\triangle^{+}&\frac{\triangle^{0}}{\sqrt{2}}&\frac{H^{+}_{S}}{\sqrt{2}}&\frac{H^{0}_{S}}{\sqrt{2}}\\\ \frac{\triangle^{0}}{\sqrt{2}}&\triangle^{-}&\frac{H^{\prime 0}_{S}}{\sqrt{2}}&\frac{H^{+}_{1S}}{\sqrt{2}}\\\ \frac{H^{{}^{\prime\prime}0}_{S}}{\sqrt{2}}&\frac{H^{-}_{2S}}{\sqrt{2}}&H_{1S}^{0}&\frac{H_{1S}^{{}^{\prime}+}}{\sqrt{2}}\\\ \frac{H^{{}^{\prime}0}_{1S}}{\sqrt{2}}&\frac{H_{2S}^{{}^{\prime}-}}{\sqrt{2}}&\frac{H^{{}^{\prime}-}_{S}}{\sqrt{2}}&\sigma_{S}^{0}\\\ \end{array}\right)$ | 1
Table 2:  The complete anomaly free fermions content and scalar sectors in the
flipped 341 model with their flavors (F), both
$(1,3,0)+(1,2,\frac{-1}{2})+(1,1,-1)$ and
$(1,3,0)+(1,2,\frac{-1}{2})+(1,2,\frac{1}{2})+(1,1,0)$ are the decomposition
of $(1,6,\frac{-1}{3})$ ref8  and $(1,8,0)$ respectively.

##  3 Anomalies cancellation

Based on to the fact that $\mathrm{Tr}[T^{a}]=0$ and $\mathrm{Tr}[\tau^{a}]=0$
where $T^{a}$ and $\tau^{a}$ are the generators of the Lie gauge groups SU(3)
and SU(4) respectively, all the triangle gauge anomalies are automatically
canceled expect those of the non-trivial ones: $[SU(4)_{L}]^{3}$,
$[SU(3)_{C}]^{2}\otimes U(1)_{X}$, $[SU(4)_{L}]^{2}\otimes U(1)_{X}$,
$[SU(3)_{C}]^{3}$, $[U(1)_{X}]^{3}$ and $[Grav]^{2}\otimes U(1)_{X}$. From the
particles content shown in table 2, we notice that our flipped 341 model is
free from all the gauge triangle anomalies. In fact, the total contribution of
the $[SU(4)_{L}]^{3}$ anomaly comes from ref99 :

$\mathcal{A}^{abc}(4_{L})\bigg{(}\sum_{Q_{mL},f_{iL}}4_{L}-\sum_{Q_{nL}}\bar{4}_{L}\bigg{)}=\mathcal{A}^{abc}(4_{L})(n_{4L}-n_{\bar{4}_{L}})$
(9)

where the anomaly coefficient
$\mathcal{A}^{abc}(4_{L})=\mathrm{Tr}(T^{a}_{L},\\{T^{b}_{L},T^{c}_{L}\\})$
with $\mathcal{A}^{abc}(4_{L})=-\mathcal{A}^{abc}(\bar{4}_{L})$ and $4_{L}$
respectively $\bar{4}_{L}$ are SU(4) quadruplet and anti-quadruplets
fundamental representations. Here $n_{4L}$ and $n_{\bar{4}_{L}}$ are the
number of left-handed fermions quadruplets and anti-quadruplets respectively.
This anomaly cancels only if the number of the quadruplets in the fundamental
representation $4_{L}$ equals to the number of the quadruplets in the
conjugate fundamental representation $\bar{4}_{L}$ as it is the case in our
model (see Table 2). We remind that the 10-plet $\psi_{e}$ contributes as much
as eight quadruplets in the group SU(4) ref9
$(\mathcal{A}(10)=8\mathcal{A}(4))$, together with the remaining lepton
generations $\psi_{\mu}$ and $\psi_{\tau}$, the contribution of all lepton
generations equals to $10\mathcal{A}(4)$. While, the arrangement of the three
quarks families in the fundamental conjugate representation $\overline{4}$
makes their contribution equal to $9\mathcal{A}(4)$. Thus, to ensure the
cancellation of the $[SU(4)_{L}]^{3}$ anomaly, a new exotic lepton
$\widetilde{\psi}$ lies in the conjugate representation $\overline{4}$ must be
introduced.
The cancellation of the $(SU(3)_{C})^{3}$ anomaly requires the introduction of
the charge conjugate of each quark field as an $SU(4)_{L}\otimes U(1)_{X}$
singlet for which the quantum number X coincides with the electric charge Q
ref12 .
The flipped 341 model is free from the $SU(4)_{L}\otimes U(1)_{X}$,
$SU(3)_{C}\otimes U(1)_{X}$, $(Grav)^{2}\otimes U(1)_{X}$ and $(U(1)_{X})^{3}$
anomalies only if its particles content satisfied the following conditions
respectively:

$\displaystyle\frac{1}{2}\bigg{(}\sum X_{l}^{L}$ $\displaystyle+$
$\displaystyle 3\sum X_{q}^{L}\bigg{)}=0,$ (10)
$\displaystyle\frac{1}{2}\bigg{(}4\sum X_{q}^{L}$ $\displaystyle-$
$\displaystyle\sum_{Sing}X_{q}^{R}\bigg{)}=0,$ (11) $\displaystyle 4\sum
X_{l}^{L}+12\sum X_{q}^{L}$ $\displaystyle-$ $\displaystyle 3\sum
X_{q}^{R}-\sum X^{R}_{l}=0,$ (12) $\displaystyle
4\sum(X_{l}^{L})^{3}+12\sum(X_{q}^{L})^{3}$ $\displaystyle-$ $\displaystyle
3\sum_{Sing}(X^{R}_{q})^{3}-\sum_{Sing}(X^{R}_{l})^{3}=0.$ (13)

where $X_{l(q)}^{L}$ and $X_{l(q)}^{R}$ are the quantum numbers associated to
the $U(1)_{X}$ group of the left $(L)$ and right $(R)$ handed leptons $(l)$
and quarks $(q)$. The factors 4 appear to take into account all the components
of the quadruplets and anti-quadruplets ref99 , the factors 3 appear because
the quarks are in $SU(3)_{C}$ triplets ref99  and the factors 12 represent the
multiplication of $3\times 4$, the quark generations are quadruplets and
triplets under the group SU(4) and SU(3) respectively.
The special content of the flipped 341 model ensures the cancellation of the
$(Grav)^{2}\otimes U(1)_{X}$ anomaly which requires that the sum of all the
$U(1)_{X}$ charges yields to zero as it mentioned in equation (13). Table 3
shows the cancellation of the gauge anomalies in the flipped 341 model. Notice
that the sum of all rows of each column in the table 3 multiplied by the
number of flavors (F) vanishes. Therefore, the model is gauge anomalies free.

Field | $[SU(3)_{C}]^{3}$ | $[SU(4)_{L}]^{3}$ | $[SU(4)_{L}]^{2}U(1)_{X}$ | $[SU(3)_{C}]^{2}U(1)_{X}$ | $[U(1)_{X}]^{3}$ | $(Grav)^{2}U(1)_{X}$ | F
---|---|---|---|---|---|---|---
$\psi_{e}$ | 0 | 4 | 0 | 0 | 0 | 0 | 1
$\psi_{\alpha}$ | 0 | $\frac{1}{2}$ | $\frac{-1}{4}$ | 0 | $\frac{-1}{2}$ | -2 | 2
$\widetilde{\psi}$ | 0 | $\frac{-1}{2}$ | $\frac{-1}{4}$ | 0 | $\frac{-1}{2}$ | -2 | 1
$l_{\alpha}^{c}$ | 0 | 0 | 0 | 0 | 1 | 1 | 6
$Q_{\alpha}$ | 2 | $\frac{-3}{2}$ | $\frac{1}{4}$ | $\frac{1}{3}$ | $\frac{1}{18}$ | 2 | 3
$u^{c}_{\alpha}$ | $\frac{-1}{2}$ | 0 | 0 | $\frac{-1}{3}$ | $\frac{-8}{9}$ | -2 | 6
$d^{c}_{\alpha}$ | $\frac{-1}{2}$ | 0 | 0 | $\frac{1}{6}$ | $\frac{1}{9}$ | 1 | 6
Table 3:  Gauge anomalies fields contributions in the flipped 341 model.

##  4 Fermions masses

Following ref ref8 , we assume that there exists a stable, charge-preserving
vacuum state. In this spirit, we allow all neutral scalar quadruplets
components in Table 2 to have a non-zero vacuum expectation value (VEV) ref8 :

$\langle\phi_{1}\rangle=\left(\begin{array}[]{ccc}0\\\ k_{1}\\\ n_{1}\\\
0\end{array}\right)\langle\phi_{2}\rangle=\left(\begin{array}[]{ccc}0\\\
k_{2}\\\ n_{2}\\\
0\end{array}\right)\langle\phi_{3}\rangle=\left(\begin{array}[]{ccc}k_{3}\\\
0\\\ 0\\\
n_{3}\end{array}\right)\langle\phi_{4}\rangle=\left(\begin{array}[]{ccc}k_{4}\\\
0\\\ 0\\\ n_{4}\end{array}\right)$

Note that in order to eliminate the unwanted Yukawa Lagrangian interactions
(nonphysical terms) such as terms like $\nu_{e}e$, $\Sigma^{0}e$ and many
others, only the elements $S_{33}$ and $S_{44}$ develop vacuum expectation
values (VEVs). Therefore, the VEV of the scalar 15-plet S is:

$\langle S\rangle=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\
0&0&n_{1S}&0\\\ 0&0&0&n_{2S}\end{array}\right)$

The gauge symmetry in the flipped 341 model is broken to the Standard Model
one via the following steps:

$\displaystyle SU(3)_{c}\otimes S$ $\displaystyle U(4)_{L}\otimes U(1)_{X}$
(14) $\displaystyle\Downarrow n_{3},n_{4},n_{1S},n_{2S}$ $\displaystyle
SU(3)_{c}\otimes S$ $\displaystyle U(3)_{L}\otimes U(1)_{X^{\prime}}$
$\displaystyle\Downarrow n_{1},n_{2}$ $\displaystyle SU(3)_{c}\otimes S$
$\displaystyle U(2)_{L}\otimes U(1)_{Y}$ $\displaystyle\Downarrow
k_{1},k_{2},k_{3},k_{4}$ $\displaystyle SU(3)_{c}$ $\displaystyle\otimes
U_{em}$

Where $k_{i}(i=1,2,3,4)\ll n_{1},n_{2}\ll n_{i}(i=3,4,1S,2S)$. Regarding the
Yukawa Lagrangian $\mathcal{L}_{Y}$ it has the form:

$\mathcal{L}_{Y}=\mathcal{L}_{Leptons}+\mathcal{L}_{Quarks}$ (15)

The lepton sector allows the combination of $SU(4)_{L}$ quadruples and anti
quadruplets as:

$\displaystyle\mathcal{L}_{Leptons}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{2}y_{\alpha\beta}^{l(i)}\psi_{\alpha}l^{c}_{\beta}\phi_{i}+h.c$
(16)

and for quarks:

$\mathcal{L}_{Quarks}=\sum_{i=1}^{2}y_{\alpha\beta}^{u(i)}Q_{\alpha}u^{c}_{\beta}\phi^{*}_{i}+\sum_{i=3}^{4}y^{d(i)}_{\alpha\beta}Q_{\alpha}d^{c}_{\beta}\phi^{*}_{i}+h.c$
(17)

where $\beta$ and $\alpha$ stand for the flavors of the fermions fields.
Notice that an additional term in the Lagrangian (16) is introduced in order
to generate a mass to the exotic leptons $\Sigma$:

$\mathcal{L}_{Leptons}=\frac{y^{\prime}}{\Lambda}\psi_{e}^{aa^{\prime}}(\psi_{e})^{bb^{\prime}}(S)^{cc^{\prime}}S^{dd^{\prime}}\epsilon_{abcd}\epsilon_{a^{\prime}b^{\prime}c^{\prime}d^{\prime}},$
(18)

Where ${y^{\prime}}$ is the Yukawa coupling, $\Lambda$ represents the cutoff
scale where new physics is expected and $\epsilon_{abcd}$ is the Levi Civita
tensor.
Regarding the masses of quarks, substituting the VEVs into the Yukawa
Lagrangian (17), we extract the tree level masses of the up and down quarks
respectively:

$m^{u}=\left(\begin{array}[]{cccc}y_{\alpha\beta}^{u(1)}k_{1}+y_{\alpha\beta}^{u(2)}k_{2}\\\
y_{\alpha\beta}^{u(1)}n_{1}+y_{\alpha\beta}^{u(2)}n_{2}\\\ \end{array}\right)$
$m^{d}=\left(\begin{array}[]{cccc}y_{\alpha\beta}^{d(3)}k_{3}+y_{\alpha\beta}^{d(4)}k_{4}\\\
y_{\alpha\beta}^{d(3)}n_{3}+y_{\alpha\beta}^{d(4)}n_{4}\\\ \end{array}\right)$

where we have written $m^{u}$ (resp. $m^{d}$) in the basis
($u_{\alpha},U_{\alpha}$) and ($u_{\beta}^{c}$)
(resp.($d_{\alpha},D_{\alpha}$) and ($d_{\beta}^{c}$)),
$y_{\alpha\beta}^{u(i)}$(i=1,2), $y_{\alpha\beta}^{d(j)}$ (j=2,3) represent
the Yukawa couplings.
At the tree level, the light charged lepton mass matrix $m^{l}$ has been
extracted from the Yukawa Lagrangian (16) and it is written in the basis
($l_{\alpha},E_{\alpha}$) and ($l^{c}_{\beta}$))

$m^{l}=\left(\begin{array}[]{ccccccccccc}y_{\alpha\beta}^{l(1)}k_{1}+y_{\alpha\beta}^{l(2)}k_{2}\\\
y_{\alpha\beta}^{l(1)}n_{1}+y_{\alpha\beta}^{l(2)}n_{2}\\\ \end{array}\right)$

The mass matrix $m^{l}$ shows that the charged leptons $\mu$ and $\tau$
acquire their masses at the tree level except the electron, it remains
massless. Moreover, notice also that in contrast with the flipped 331 model
ref8 , the following Lagrangian:

$\mathcal{L}=y_{\alpha}^{l}\psi_{e}l^{c}_{\alpha}S$ (19)

is not allowed since it is not gauge invariant under both $U(1)_{X}$ and
$SU(4)_{L}$ symmetries. Therefore, the exotic lepton $E_{e}$ and the electron
e remain massless. Thus, in order to generate masses to both $E_{e}$ and e,
one has to introduce the following effective dimension-5 operator:

$\mathcal{L}_{eff}=\frac{\lambda^{(ij)}_{\alpha
l}}{\Lambda}(\psi_{e})^{aa^{\prime}}l^{c}_{\alpha}(\phi_{i})_{b}(\phi_{j})_{b^{\prime}}\delta_{ab}\delta_{a^{\prime}b^{\prime}}.$
(20)

In this case:

$\displaystyle m_{E_{e}}$ $\displaystyle=$
$\displaystyle\frac{\lambda^{(ij)}_{E_{e}}}{\Lambda}n_{i}n_{j},$ (21)
$\displaystyle m_{e}$ $\displaystyle=$
$\displaystyle\frac{\lambda^{(ij)}_{e}}{\sqrt{2}\Lambda}(k_{i}n_{j}+n_{i}k_{j}).$
(22)

The VEVs $n_{i},n_{j}$ (i,j= 1,2), $\Lambda$ and $k_{i}$ are of $\mathcal{O}$
(TeV) and $\mathcal{O}$ (GeV) respectively. Taking $n_{j}$=$\Lambda$, then the
electron and $E_{e}$ masses are proportional to the SM symmetry breaking
$k_{i}$ and heavy scales respectively.
Concerning the exotic leptons $\Sigma^{\mp}$ and $\Sigma^{0}$, their masses
can be obtained from the Lagrangian (18) to get:

$\displaystyle m_{\Sigma^{\mp}}$ $\displaystyle=$
$\displaystyle\frac{4y^{\prime}}{\Lambda}n_{1S}n_{2S},$ $\displaystyle
m_{\Sigma^{0}}$ $\displaystyle=$
$\displaystyle\frac{-2y^{\prime}}{\Lambda}n_{1S}n_{2S}$ (23)

Where the VEVs $n_{iS}$ (i=1,2) are of $\mathcal{O}$ (TeV). In the case of our
interest the exotic leptons $\beta^{+}$, $\beta^{0}$, $N_{e}^{0}$ and
$\sigma^{+}$ receive their masses from the following Lagrangian:

$\mathcal{L}_{eff}=Y^{L}(\psi_{e})^{aa^{\prime}}(\widetilde{\psi})_{b}(\phi_{j})_{c}\delta_{ab}\delta_{a^{\prime}c}.$
(24)

where j=3,4. Notice that we have identified the charged and neutral elements
of $\widetilde{\psi}$ namely, $\widetilde{e}$, $\widetilde{\nu}$,
$\widetilde{N}^{0}$ and $\widetilde{E}^{-}$ with the charge conjugated right-
handed components of leptons already introduced in the electron generation
$\psi_{e}$ which are $\beta^{+}$, $\beta^{0}$, $N_{e}^{0}$ and $\sigma^{+}$
respectively.

$\widetilde{\psi}_{L}=\left(\begin{array}[]{ccc}\widetilde{e}_{L}^{-}\\\
\widetilde{\nu}_{L}\\\ \widetilde{N}^{0}_{L}\\\ \widetilde{E}^{-}_{L}\\\
\end{array}\right)\equiv\left(\begin{array}[]{ccc}(\beta^{+}_{R})^{c}\\\
(\beta^{0}_{R})^{c}\\\ (N^{0}_{eR})^{c}\\\ (\sigma^{+}_{R})^{c}\\\
\end{array}\right)$

We have used this identification to avoid the presence of charged exotic
lepton with masses of the order of the electroweak scale ref99 . Therefore,
the exotic leptons $\beta^{+}$, $\beta^{0}$, $N_{e}^{0}$ and $\sigma^{+}$
receive their masses denoted by $m_{L}$ and $m_{\sigma^{+}}$ respectively from
the Lagrangian (24) to get:

$\displaystyle m_{L}$ $\displaystyle=$ $\displaystyle
Y^{L}\frac{n_{j}}{\sqrt{2}},$ $\displaystyle m_{\sigma^{+}}$ $\displaystyle=$
$\displaystyle Y^{\sigma^{+}}n_{j}.$ (25)

With $Y^{L}$ and $Y^{\sigma^{+}}$ are the Yukawa coupling and L stands for the
exotic leptons $\beta^{+}$, $\beta^{0}$ and $N_{e}^{0}$.

##  5 FCNC

The ordinary versions of the 341 models predict the existence of new heavy
neutral gauge bosons $Z^{\prime}$ and $Z^{\prime\prime}$ which have universal
couplings with leptons, while, the quarks couplings with the new gauge bosons
are non-universal. The existence of a non diagonal matrix when we rotate the
flavor basis into the mass eigenstates leading to the occurrence of FCNC in
the quark sector. We suppose that the major sources contributing to the lepton
flavor violating processes come from the direct gauge interactions between the
charged leptons and new massive gauge bosons $Z^{{}^{\prime}}$ and
$Z^{{}^{\prime\prime}}$, while the others such as scalar contributions are
considered to be small and neglected ref13 .
This scheme is reversed in the flipped 341 model, since two of leptons
families are arranged differently from the remaining generations whereas the
quarks families lie in the same representation. Therefore, the FCNC occurs
only in the lepton sector at the tree level through the exchange of the gauge
bosons $Z^{\prime}$ and $Z^{\prime\prime}$. In the context of our model, we
study the lepton flavor violation processes namely $\mu\longrightarrow
ee\overline{e}$ and $\mu\longrightarrow e\gamma$. In fact, the lepton neutral
currents of $Z^{\prime}$ and $Z^{\prime\prime}$ are described by the following
Lagrangian:

$\displaystyle\mathcal{L}_{NC}\supset-g\overline{F}\gamma^{\mu}\bigg{(}T_{3\mu}A_{3}+T_{8}A_{8\mu}+T_{15}A_{15\mu}+Xg_{X}B_{\mu}\bigg{)}F$
(26)

where $T_{i}=\lambda_{i}/2$ and $\lambda_{i}$ (i=3, 8, 15) are the diagonal
Gell-Mann matrices in the group SU(4), g and $g_{X}$ represent the gauge
couplings of the $SU(4)_{L}$ and $U(1)_{X}$ respectively, X is the charge
associated to the group $U(1)_{X}$ and F runs over all fermions multiplets.
Here $T_{i}\psi_{\alpha}=\frac{1}{2}\lambda_{i}\psi_{\alpha}$ (i=3, 8, 15) and
$T_{i}\psi_{e}=\frac{1}{2}(\lambda_{i}\psi_{e}+\psi_{e}\lambda_{i})$(i=3, 8,
15) ref13 , the right-handed leptons $e_{aR}$ and $E_{aR}$ do not participate
in FCNC ref13 . Then we obtain:

$\displaystyle\mathcal{L}_{NC}$ $\displaystyle\supset$
$\displaystyle\frac{g}{2}\bigg{(}\frac{1-t_{W}^{2}}{\sqrt{3-t_{W}^{2}}}\bigg{)}\bigg{(}\overline{\nu}_{L}\gamma^{\mu}T_{\nu}\nu_{L}+\overline{l}_{L}\gamma^{\mu}T_{l}l_{L}+\overline{E}_{L}\gamma^{\mu}T_{E}E_{L}+\overline{N}_{L}\gamma^{\mu}T_{\nu}N_{L}\bigg{)}Z^{\prime}_{\mu}$
(27) $\displaystyle+$
$\displaystyle\frac{g}{2}\bigg{(}\frac{-2}{\sqrt{6+4t^{2}_{X}}}\bigg{)}\bigg{(}\overline{\nu}_{L}\gamma^{\mu}T^{\prime}_{\nu}\nu_{L}+\overline{l}_{L}\gamma^{\mu}T_{l}l^{\prime}_{L}+\overline{E}_{L}\gamma^{\mu}T^{\prime}_{E}E_{L}+\overline{N}_{L}\gamma^{\mu}T^{\prime}_{\nu}N_{L}\bigg{)}Z^{\prime\prime}_{\mu}.$

Where $t_{X}^{2}=s_{W}^{2}/(1-2s_{W}^{2})$ ref133 , $s_{W}$, $c_{W}$ and
$t_{W}$ are the sine, cosine and tangent of the electroweak mixing angle.
Whereas $T_{i}$ and $T^{\prime}_{i}$ are:

$\displaystyle T_{\nu}=T_{l}$ $\displaystyle=$
$\displaystyle\text{diag}\bigg{(}1,(-1-t_{W}^{2})/(1-t_{W}^{2}),(-1-t_{W}^{2})/(1-t_{W}^{2})\bigg{)},$
(28) $\displaystyle T_{E}$ $\displaystyle=$ $\displaystyle\text{
diag}\bigg{(}(4-2t_{W}^{2})/(1-t_{W}^{2}),(2-2t_{W}^{2})/(1-t_{W}^{2}),(2-2t_{W}^{2})/(1-t_{W}^{2})\bigg{)},$
(29) $\displaystyle T_{N}$ $\displaystyle=$
$\displaystyle\text{diag}\bigg{(}2/(1-t_{W}^{2}),0,0\bigg{)},$ (30)
$\displaystyle T^{\prime}_{\nu}$ $\displaystyle=$ $\displaystyle
T^{\prime}_{l}=T^{\prime}_{E}=\text{diag}\bigg{(}1,(1-2t_{X}^{2})/2,(1-2t_{X}^{2})/2)\bigg{)}$
(31) and (32) $\displaystyle T^{\prime}_{N}$ $\displaystyle=$
$\displaystyle\text{diag}\bigg{(}-1,-(3+2t_{X}^{2})/2,-(3+2t_{X}^{2})/2\bigg{)}.$
(33)

Changing from the flavor basis into the mass basis
$l_{L}=V_{lL}l^{\prime}_{L}$, we get the following Lagrangian:

$\mathcal{L}_{NC}\supset\frac{g}{2}\bigg{(}\frac{1-t_{W}^{2}}{\sqrt{3-t_{W}^{2}}}\bigg{)}\overline{l^{\prime}}_{L}\gamma^{\mu}\bigg{(}V^{\dagger}_{lL}T_{l}V_{lL}\bigg{)}l^{\prime}_{L}Z^{\prime}_{\mu}-\frac{g}{2}\bigg{(}\frac{-2}{\sqrt{6+4t^{2}_{X}}}\bigg{)}\overline{l^{\prime}}_{L}\gamma^{\mu}\bigg{(}V^{\dagger}_{lL}T^{\prime}_{l}V_{lL}\bigg{)}l^{\prime}_{L}Z^{\prime\prime}_{\mu}.$
(34)

Thus,

$\mathcal{L}_{NC}\supset\frac{g}{2}\bigg{(}\frac{1-t_{W}^{2}}{\sqrt{3-t_{W}^{2}}}\bigg{)}\overline{l^{\prime}}_{L}\gamma^{\mu}(V^{*}_{lL})_{\alpha
i}(V_{lL})_{\beta
j}l^{\prime}_{L}Z^{\prime}_{\mu}+\frac{g}{2}\bigg{(}\frac{-2}{\sqrt{6+4t^{2}_{X}}}\bigg{)}\overline{l^{\prime}}_{L}\gamma^{\mu}(V^{*}_{lL})_{\alpha
i}(V_{lL})_{\beta j}l^{\prime}_{L}Z^{\prime\prime}_{\mu}.$ (35)

Here $l^{\prime}$ can be $e,\nu,E,N$ and $i\neq j$ for flavor changing.
Regarding the decays $\mu\longrightarrow ee\overline{e}$ and
$\mu\longrightarrow e\gamma$ are among the rare decay processes which are used
to search for the charged lepton flavor violation . The corresponding Feynman
diagrams are shown in Figure 1. The first vertex in the left diagram shows the
tree level FCNC coupling of $Z^{\prime}$($Z^{\prime\prime}$) boson (cLF
changing), whereas the right diagram represents the one loop level process
where $l_{i}$ can be any lepton. Here we will consider the internal fermions
line to be either $\mu$ or e, so that we will have only one FCNC
$Z^{\prime}\mu e$ ($Z^{\prime\prime}\mu e$) vertex.

Figure 1:  The decays $\mu\longrightarrow ee\overline{e}$ and
$\mu\longrightarrow e\gamma$ via the neutral gauge bosons
$Z^{\prime}$($Z^{\prime\prime}$).

The transferred momentum, whose maximal value is about the muon mass, is much
smaller than $M_{Z^{\prime}(Z^{\prime\prime})}$ ref13 , therefore, the
Branching ratio of the processes $\mu\longrightarrow ee\overline{e}$ and
$\mu\longrightarrow e\gamma$ through the exchange of $Z^{\prime}$ and
$Z^{\prime\prime}$ are found to be respectively:

$\displaystyle Br(\mu\longrightarrow ee\overline{e})$ $\displaystyle=$
$\displaystyle
M^{4}_{W}|U_{e\mu}|^{2}\bigg{(}\frac{(1-t_{W}^{2})^{2}}{2M_{Z^{\prime}}^{2}(3-t_{W}^{2})}+\frac{2}{M_{Z^{\prime\prime}}^{2}(6+4t^{2}_{X})}\bigg{)}^{2},$
(36) $\displaystyle Br(\mu\longrightarrow e\gamma)$ $\displaystyle=$
$\displaystyle\frac{24\alpha}{16\pi}M^{4}_{W}|U_{e\mu}|^{2}\bigg{(}\frac{(1-t_{W}^{2})^{2}}{M_{Z^{\prime}}^{2}(3-t_{W}^{2})}+\frac{4}{M^{2}_{Z^{\prime\prime}}(6+4t^{2}_{X})}\bigg{)}^{2}$
(37)

where $U_{e\mu}$ is the mixing matrix for the lepton. We have neglected the
contribution coming from the electron internal line as it is proportional to
$m_{e}/m_{\mu}$ ref16  and considered the electrons as massless particles.
Using the experimental upper limit of the branching ratio
$\text{Br}(\mu\longrightarrow ee\overline{e})\leq 10^{-12}$ ref15  together
with the current experimental limit $\text{Br}(\mu\longrightarrow e\gamma)\leq
10^{-13}$ ref14 , we obtain, an upper limit on the lepton flavor violating
matrix element $|U_{\mu e}|\leq 1.66\times
10^{-3}\frac{M^{2}_{Z^{\prime}}}{TeV^{2}}$ and a stringent bound on the gauge
bosons masses $M_{Z^{\prime}}\leq 0.597M_{Z^{\prime\prime}}$. Now if
$M_{Z^{\prime}}$$\approx$ $\mathcal{O}$(1 TeV), one has $|U_{\mu
e}|$$\approx$$1.66\times 10^{-3}$ and
$M_{Z^{\prime\prime}}$$\approx$$\mathcal{O}$(1.68 TeV). Figure 2 shows the
variation of $M_{Z^{\prime}}$ and $M_{Z^{\prime\prime}}$ as a function of
$|U_{\mu e}|$ and the dashed areas represent the allowed region where the
constraints are verified.

Figure 2:  The variation of $M_{Z^{\prime}}$ and $M_{Z^{\prime\prime}}$ as a
function of $|U_{e\mu}|$.

##  6 Conclusion

In this work, we have constructed a new unique model based on the gauge group
$SU(3)\otimes SU(4)\otimes U(1)$ called the flipped 341 model where all the
three quarks families arrange in the same representation whereas leptons
generations are not. The anomalies cancellation requires the introduction of
new extra exotic leptons a 10 plet $\psi_{e}$ and a quadruplet
$\widetilde{\psi}$.
In the flipped 341 model, a dangerous flavor changing neutral current can
appear in the lepton sector through the interactions with the neutral gauge
bosons $Z^{\prime}$ and $Z^{\prime\prime}$, whereas, it is absent in the quark
sector which makes this model unable to explain the anomalies of both kaon and
B decays ref8  in contrast to the ordinary versions of the 341 models. The
comparison of our results with the experimental data concerning the branching
ratios of the two processes $\mu\longrightarrow ee\overline{e}$ and
$\mu\longrightarrow e\gamma$ lead to a stringent constraint on the matrix
element and on the masses of the gauge bosons.

## Acknowledgments

We would like to thank Prof. Renato M. Fonseca for fruitful discussions. We
also thank Uladzimir Khasianevich for useful and clarifying comments. We are
very grateful to the Algerian ministry of higher education and scientific
research and DGRSDT for the financial support.

## References

  * (1) A. G. Dis, P. R. D. Pinheirob, C. A. de S. Pires and P. S. Rodrigues da Silva, _A compact 341 model at TeV scale_ , Ann.Phys.349, 232 (2014), arXiv:1309.6644 [hep-ph].
  * (2) N.Mebarki, M.Djouala, J.Mimouni and H.Aissaoui, _A New Compact 341 Model: Higgs Decay Modes_ , FISICPAC-2018 - University of Sharjah - UAE,11-13 November 2018.
  * (3) M.Djouala and N. Mebarki, _Diphton Higgs decays in the compact 341 model_ , $2^{nd}$ JIPC conference, 14-15 December 2016, Constantine, Algeria.
  * (4) M.Djouala and N. Mebarki, _Exploring the Higgs Bosons in the Compact 341 Model_ , $1^{st}$ Mediterranean Conference on Higgs Physics, 23-26 September 2019, Tangier-Morocco.
  * (5) J.M.Cabarcas and J.A.Rodriguez, _Electric charge quantization in models_ , Mod. Phys. Lett. A29, 1450032 (2014), arXiv:1303.5332 [hep-ph].
  * (6) William A. Ponce and Luis A. Sanchez, _Systematic study of the $SU(3)_{c}\otimes SU(4)_{L}\otimes U(1)_{X}$ gauge symmetry_, Int. J. Mod. Phys. A,, Issue 06, pp. 435-447 (2007), arxiv:0607175 [hep-ph].
  * (7) F.Pisano and V.Pleitez, _model for the electroweak interactions_ , Phys.Rev.D51 , 3865, (1995), arXiv:9401272 [hep-ph].
  * (8) A.Jaramiloo and L.A.Sanchez, _FCNC, CP violation and implications for some rare decays in an extension of the standard model_ , Phys.Rev.D84.115001, (2011) arXiv:1110.3363 [hep-ph].
  * (9) R. M.Fonseca and M. Hirsch, _A flipped 331 model, JHEP08, 003 (2016)_ , arxiv:1606.01109 [hep-ph].
  * (10) N.Yamatsu, _Finite dimensional Lie algebras and their representations for unified model building_ , arXiv:1511.08771 [hep-ph].
  * (11) S. Descotes-Genon, M. Moscati and G. Ricciardi, _Non-minimal 331 model for Lepton Flavour Universality Violation in $b\longrightarrow sll$ decays_, PhysRevD.98.115030, arXiv:1711.03101 [hep-ph]
  * (12) H.N.Long, L.T.Hue and D.V.Loi, _Electroweak theory based on $SU(4)_{L}\times U(1)_{X}$ gauge group, Phys.Rev. D94 (2016) no.1, 015007_, arXiv:1605.07835[hep-ph].
  * (13) W.A.Ponce and L.A.Sanchez, _Systematic study of the gauge symmetry_ , Mod.Phys.Lett.A22, 435 (2007), arXiv:0607175 [hep-ph].
  * (14) D.T.Huong, D.N.Dinh, L.D.Thien, and P.Van Dong, _Dark matter and flavor changing in the flipped 3-3-1 model_ , JHEP08, 051 (2019), arXiv:1906.05240 [hep-ph].
  * (15) Duong Van Loi and Phung Van Dong, _Kinetic mixing effect in noncommutative $B-L$ gauge theory_, JHEP09, 054 (2019), arXiv:1906.10577 [hep-ph].
  * (16) M. Tanabashi et al. (Particle Data Group), Phys. Rev. D98, 030001 (2018).
  * (17) S. Descotes-Genon, M. Moscati and G. Ricciardi, _Non minimal 331 model for lepton flavor universality violation in decays_ , Phys.Rev.D98, 115030 (2018).
  * (18) R.Mohanta, _Effect of FCNC mediated Z boson on lepton flavor violating decays_ , Eur.Phys. J. C71, 1625 (2011), arXiv:1011.4184[hep-ph] .


# Gap resonance in the classical dynamics of the current-biased Josephson
tunnel junctions

Dmitri V. Averin  Department of Physics and Astronomy, Stony Brook University,
SUNY, Stony Brook, NY 11794-3800

###### Abstract

This work reports a novel time-domain expression for the current-response
kernels of a Josephson tunnel junction between BCS superconductors with in
general different energy gaps, and its use to simulate the classical dynamics
of such junctions. The simulations show a dynamic regime characterized by the
resonance between the Josephson oscillations and the gap oscillations in the
asymptotics of the current kernels that has not been studied previously. The
resonance manifests itself as a hysteresis in the dc current-voltage
characteristics of the current-biased junctions in the voltage range above the
energy gap, in addition to the usual “inertial” hysteresis characteristic for
tunnel junctions at voltages below the energy gap. Features of the IV curves
related to the gap resonance, including the above-the-gap hysteresis, should
manifest themselves in many structures and devices utilizing high-quality
Josephson tunnel junctions with relatively large, but achievable current
densities.

##  I Introduction

A tunnel junction: two conductors separated by a layer of insulator that is
sufficiently thin to allow for charge transfer between the conductors at a
rate that is large enough to produce a noticeable current [1, 2], has a long
and still continuing history as a source of novel phenomena in physics, and a
tool for both practical devices and scientific experiments. Probably the best
studied and widely used are the Josephson tunnel junctions between
superconductors, which exhibit Josephson effect [3, 4], a quantum-coherent
transfer of Cooper pairs, which leads to supercurrent flow across the tunnel
barrier. An important feature of this effect, which makes possible its various
applications, is that the Josephson supercurrent has the same magnitude as the
typical current associated with electron tunneling in the junction [5, 3].
This is possible because the magnitude of the supercurrent, as a quantum-
coherent process, is determined by the amplitude of tunneling, in contrast to
the probability of tunneling for electron current. The amplitude of the two-
electron transfer of a Cooper pair has the same magnitude as the probability
$D$ of one-electron transmission through a tunnel barrier, and is very small,
$D\ll 1$, for typical tunnel barriers.

Josephson tunnel junctions possess several other characteristics, starting
with the strongly nonlinear dynamics of the dynamic junction variables:
Josephson phase difference $\phi(t)$ and the voltage $V(t)$ across the
junction, that make them ideal building blocks for various devices of
superconductor electronics, including digital circuits [6, 7]. Another
features of the high-quality tunnel junctions with $D\ll 1$, is that there is
a wide range of voltages $V$ around zero, where the dissipative electron
tunneling is strongly or, ideally, completely suppressed, and the junction
dynamics is reduced to the quantum-coherent Cooper-pair tunneling. Tunnel
junctions in this regime demonstrate various effects of the macroscopic
quantum dynamics of the Josephson phase $\phi$ [8, 9, 10, 11, 12, 13, 14],
which at present play most important role in the development of
superconducting quantum computing circuits (see, e.g., [15, 16, 17, 18, 19]).
As a less daunting task, the same suppression of the dissipative electron
tunneling makes it possible to use Josephson tunnel junctions for the
development of classical circuits for thermodynamically reversible computing
[20].

Classical current response of a Josephson tunnel junction to an arbitrary
fixed time-dependent voltage $V(t)$ is well understood for a long time [21,
22] and is described in details in classic textbooks [23, 24]. Nevertheless,
although the assumption of the fixed bias voltage $V(t)$ is natural for the
theoretical calculation of the current in the junction, in practice, typical
bias conditions are quite different, with the simplest case being the current
bias. Nonlinear nature of the Josephson dynamics makes the task of analyzing
the junction transport properties in these realistic situations quite
nontrivial even if the basic voltage-biased current response is known. Since
the general nonlinear dynamics of the Josephson phase $\phi$ is described more
naturally in time domain, additional layer of difficulty is added to this
problem by the fact that in the basic theoretical description [21, 22], the
current response is more naturally calculated in the energy (i.e., frequency)
representation. Expression for the current response of a Josephson tunnel
junction directly in time domain was obtained only for junctions with equal
energy gaps of the two electrodes [25]. This work suggests a more general
expression for the time-domain current response of a Josephson tunnel junction
between two BCS superconductors with different energy gaps. This expression
makes it possible to analyze more conveniently Josephson dynamics of the
current-biased junctions and calculate the current-voltage characteristics of
these junctions. The calculations show a novel “gap resonance” regime of the
Josephson dynamics in current-biased tunnel junctions.

Next Section discusses the main result of the first part of this work,
expression for the time-domain current response of a tunnel junction, which is
derived in Appendix A. Section III describes the application of this
expression to the simulation of the Josephson dynamics of the current-biased
tunnel junctions and calculates the current-voltage characteristics of the
junction under several different conditions. Section IV estimates the junction
parameters under which the gap-resonance regime of the Josephson dynamics is
realized, and concludes. Appendix B gives some details of the numerical
procedure used in the simulations of the junction dynamics.

##  II Current kernels

Consider a Josephson tunnel junction between the two BCS superconductors with
in general different energy gaps $\Delta_{1}$ and $\Delta_{2}$. If the bias
voltage $V(t)$ between the superconductors has an arbitrary dependence on time
$t$, the tunnel current $I_{T}(t)$ that flows in the junction can be expressed
through the time-dependent Josephson phase difference $\phi(t)$ across the
junction (see, e.g., [23, 24]):

$\displaystyle
I_{T}(t)=\int_{-\infty}^{t}dt^{\prime}\Big{[}I_{p}(t-t^{\prime})\sin\frac{\phi(t)+\phi(t^{\prime})}{2}\;\;\;\;\;\;\;\;$
$\displaystyle+I_{qp}(t-t^{\prime})\sin\frac{\phi(t)-\phi(t^{\prime})}{2}\Big{]}\,,\;\;\;\;\dot{\phi}(t)=2eV(t)/\hbar\,.\;\;\;$
(1)

As shown in Appendix A, the pair current $I_{p}$ and the quasiparticle
$I_{qp}$ kernels in this equation have the following analytical form for the
tunnel junction with the normal-state resistance $R_{N}$:

$\displaystyle I_{p}(t)=-\frac{\pi\Delta_{1}\Delta_{2}}{2e\hbar
R_{N}}\Big{[}J_{0}\big{(}\frac{\Delta_{1}t}{\hbar}\big{)}A_{0}\big{(}\frac{\Delta_{2}t}{\hbar},\beta_{2}\big{)}+A_{0}\big{(}\frac{\Delta_{1}t}{\hbar},\beta_{1}\big{)}$
$\displaystyle\cdot
J_{0}\big{(}\frac{\Delta_{2}t}{\hbar}\big{)}\Big{]}\,,\;\;\;I_{qp}(t)=\frac{-\hbar}{eR_{N}}\delta^{\prime}(t-0)+\frac{\pi\Delta_{1}\Delta_{2}}{2e\hbar
R_{N}}\;\;\;\;\;$ (2)
$\displaystyle\cdot\Big{[}J_{1}\big{(}\frac{\Delta_{1}t}{\hbar}\big{)}A_{1}\big{(}\frac{\Delta_{2}t}{\hbar},\beta_{2}\big{)}+A_{1}\big{(}\frac{\Delta_{1}t}{\hbar},\beta_{1}\big{)}J_{1}\big{(}\frac{\Delta_{2}t}{\hbar}\big{)}\Big{]}\,.\;\;\;$

The functions $J_{0,1}$ here are the Bessel functions of the first kind, which
for $z>0$ can be defined by the following relations [26]:

$J_{0}(z)=\frac{2}{\pi}\int_{1}^{\infty}\frac{dx\sin
zx}{\sqrt{x^{2}-1}}\,,\;\;J_{1}(z)=-\frac{d}{dz}J_{0}(z)\,,$ (3)

while the functions $A_{0,1}$ depend on the temperature $T$ of the
superconductors through the ratios $\beta_{j}\equiv\Delta_{j}/2T$, $j=1,2$,
and are defined, also for $z>0$, as

$A_{0}(z,\beta)=-\frac{2}{\pi}\int_{1}^{\infty}\frac{dx\cos
zx}{\sqrt{x^{2}-1}}\tanh\beta x\,,$ (4)
$A_{1}(z,\beta)=-\frac{d}{dz}A_{0}(z,\beta)=-\frac{2}{\pi}\int_{1}^{\infty}\frac{dxx\sin
zx}{\sqrt{x^{2}-1}}\tanh\beta x\,.$

The first, $\delta^{\prime}$ term in the quasiparticle kernel $I_{qp}$ (2)
accounts for the normal-state Ohmic contribution $I_{N}$ to the current:
$I_{N}(t)=V(t)/R_{N}$.

At low temperatures, when $\tanh\beta x\equiv 1$, the functions $A_{j}$ are
reduced to the Bessel functions $Y_{j}(z)$, $j=0,1$, of the second kind:

$A_{j}(z,\beta\rightarrow\infty)=Y_{j}(z)\,.$

As a result, in the low-temperature limit, the current kernels (2) are
expressed through the well-studied Bessel functions only:

$\displaystyle I_{p}(t)=-\frac{\pi\Delta_{1}\Delta_{2}}{2e\hbar
R_{N}}\Big{[}J_{0}\big{(}\frac{\Delta_{1}t}{\hbar}\big{)}Y_{0}\big{(}\frac{\Delta_{2}t}{\hbar}\big{)}+Y_{0}\big{(}\frac{\Delta_{1}t}{\hbar}\big{)}J_{0}\big{(}\frac{\Delta_{2}t}{\hbar}\big{)}\Big{]},$
$\displaystyle
I_{qp}(t)=\frac{-\hbar}{eR_{N}}\delta^{\prime}(t-0)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$
(5) $+\frac{\pi\Delta_{1}\Delta_{2}}{2e\hbar
R_{N}}\Big{[}J_{1}\big{(}\frac{\Delta_{1}t}{\hbar}\big{)}Y_{1}\big{(}\frac{\Delta_{2}t}{\hbar}\big{)}+Y_{1}\big{(}\frac{\Delta_{1}t}{\hbar}\big{)}J_{1}\big{(}\frac{\Delta_{2}t}{\hbar}\big{)}\Big{]}\,.$

This result generalizes to different energy gaps $\Delta_{1}\neq\Delta_{2}$
the equations obtained by Harris [25] for equal gaps. Also, one can check that
the pair and the quasiparticle kernels (5) reproduce the frequency-domain
expressions for these kernels at $T=0$ in terms of the elliptic integrals [21,
22] in the cases for which the corresponding integrals of the Bessel functions
[27] are available (see, e.g., Eq. (13) below).

As a brief aside, one can note that the generalization to different energy
gaps implies that Eqs. (1) to (5) can be used to describe the time-dependent
current response of a normal metal/superconductor (NS) junction in the tunnel
limit. Indeed, for $\Delta_{2}=0$, the pair current $I_{p}$ vanishes, while in
the limit $\Delta_{2}\rightarrow 0$, the quasiparticle kernel is obtained from
$A_{1}(z,\beta_{2})$ evaluated for $z,\beta_{2}\rightarrow 0$, with
$z/\beta_{2}=2Tt/\hbar$. For $z\rightarrow 0$ (but $z\neq 0$), the integral
$A_{1}(z,\beta_{2})$ is determined by the range of large values of the
integration variable, where it is reduced to

$\displaystyle\Delta_{2}A_{1}(z,\beta_{2})=-\frac{2\Delta_{2}}{\pi}\int_{0}^{\infty}dx\sin
zx\tanh\beta_{2}xe^{-\epsilon x}\Big{|}_{\epsilon\rightarrow 0}$
$\displaystyle=-\frac{\Delta_{2}}{\pi\beta_{2}}\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}z/2\beta_{2}}{n^{2}+(z/2\beta_{2})^{2}}=-\frac{2T}{\sinh(\pi
Tt/\hbar)}\,.\;\;\;$ (6)

This gives for the current-response kernel of the NS tunnel junction:

$I_{qp}^{(NS)}(t)=\frac{-\hbar}{eR_{N}}\delta^{\prime}(t-0)-\frac{\pi
T\Delta}{e\hbar R_{N}}\frac{J_{1}(t\Delta/\hbar)}{\sinh(\pi Tt/\hbar)}\,.$ (7)

This expression coincides with the $D\rightarrow 0$ limit of the general
expression for the time-dependent response of an NS junction with an arbitrary
electron transparency $D$ [28].

The next Section describes the application of the expressions for the current
kernels presented in this Section to the simulations of the classical
Josephson dynamics of the current-biased junctions. An outline of the
derivation of these expressions is discussed in Appendix A.

##  III Gap resonance in the dynamics of the current-biased junctions

As an application of the current kernels obtained in the previous Section, I
now discuss the simulations of the classical dynamics of a Josephson phase
$\phi$ across a current-biased Josephson tunnel junction formed by the BCS
superconductors with in general different gap energies $\Delta_{1}$ and
$\Delta_{2}$. In this work, the simulations are limited to the low-temperature
regime, when the current kernels are given by Eq. (5) containing only the
Bessel functions. The “brute-force” approach to the simulation of the junction
dynamics employed here consists of the direct numerical evaluation of the
current kernels (as described in the Appendix B), their direct numerical
integration at each step of the time evolution, and subsequent solution of the
integro-differential equations for the junction dynamics that are described
below. This approach differs from the existing approaches to the simulation of
the Josephson dynamics in the tunnel-junction model [30, 31, 32, 33], which
use the scheme that avoids the need to solve the integro-differential
equations, and in particular, to evaluate the integrals (1) for the tunnel
current at each time step of the Josephson evolution. The scheme is based on
the approximate representation of the current kernels by the sums of the
exponentials, which effectively makes it possible to take the integrals
analytically. An important advantage of this computation scheme is that it is
very efficient and enables one to treat practical present-day circuits of
superconductor electronics with the number of junctions that can reach into
hundreds of thousands. On the other hand, because of the singular nature of
the current kernels which are characterized by the very slow and oscillatory
decay with time, this approach can create some uncertainty as to whether the
features observed in the numerical simulations are real or the consequences of
the approximations used for the current kernels [34]. While the brute-force
scheme adopted in this work is much slower, it does not create this
uncertainty. The main new result of the simulations presented below is the
demonstration of the “gap resonance” regime of the junction dynamics which
produces the hysteresis in the current-voltage characteristics of the current-
biased junction appearing in the range where the dc component of the voltage
$V$ across the junction is above the gap value
$V_{g}=(\Delta_{1}+\Delta_{2})/e$. In what follows, I present the results of
the detailed numerical simulations of the dynamics of the current-biased
Josephson tunnel junctions that show the gap resonance, and provide the semi-
quantitative analytical explanation of the above-the-gap hysteresis.

Classical Josephson dynamics of a Josephson tunnel junction shunted by an
Ohmic resistor $R_{S}$ and biased by in general time-dependent external
current $I(t)$ (see the equivalent circuit in the inset in Fig. 1) is governed
by the standard set of coupled dynamic equations for the phase $\phi$ across
the junction and the instantaneous voltage $V$ on the junction capacitance
$C$:

$\dot{\phi}=\frac{2eV}{\hbar}\,,\;\;\;C\dot{V}=I(t)-\frac{V}{R_{S}}-I_{T}(t)\,.$
(8)

Here, the tunnel current $I_{T}$ is given by Eq. (1) with the current kernels
(2). Since the calculation of tunnel current $I_{T}$ (1) requires integration
over time, the resulting evolution equations (8) for the phase $\phi$ are the
integro-differential equations. The first step towards their numerical
solution is to express them in the dimensionless form. To do this, it is
natural to normalize time $t$ to the frequency $\Omega$ of the gap
oscillations, $\Omega=(\Delta_{1}+\Delta_{2})/\hbar$, and the voltage $V$ – to
half of the gap voltage $V_{t}=(\Delta_{1}+\Delta_{2})/e$. Then, it is
convenient to introduce the auxiliary dimensionless parameters
$\alpha=1+R_{N}/R_{S}$ characterising the degree to which the junction is
shunted, with $\alpha=1$ describing the unshunted junction, and
$\lambda=(\Delta_{1}+\Delta_{2})^{2}/(\pi\Delta_{1}\Delta_{2})$ characterising
the magnitude of the difference between the two energy gaps of the junction
electrodes. With this definition, $\lambda=4/\pi$ corresponds to the junction
with equal energy gaps.

The main physical parameter that controls the nature of the phase dynamics in
the tunnel junction is the ratio of the time constant $R_{N}C$ of the junction
capacitance to the gap time:

$\beta=\frac{\Delta_{1}+\Delta_{2}}{\hbar}R_{N}C\,.$ (9)

For the purpose of the discussion of the junction dynamics below, it will be
convenient to view $\beta$ simply as the parameter characterizing the
magnitude of the junction capacitance, ranging from the negligibly small
capacitance at $\beta=0$ to large capacitance for $\beta\gg 1$. It should be
noted, however, that this view is misleading in one respect. For typical
tunnel junctions, the time constant $R_{N}C$ (and therefore $\beta$) is
independent of the junction area, in contrast to the absolute magnitude of the
junction capacitance, because the area dependence of the junction resistance
$R_{N}$ compensates that of the capacitance. Both the time constant $R_{N}C$
and $\beta$ depend only on the effective electron transparency of the tunnel
barrier that is characterized usually through the “critical current density”,
decreasing with increasing current density as the inverse of it.

With this set of the dimensionless parameters and variables, the dynamic
equations (8) with the appropriate current kernels take the form:

$\displaystyle\beta\lambda\dot{v}=i(t)-\alpha\lambda
v+\int_{0}^{\infty}d\tau\Big{[}i_{p}(\tau)\sin\frac{\phi(t)+\phi(t-\tau)}{2}$
$\displaystyle-
i_{qp}(\tau)\sin\frac{\phi(t)-\phi(t-\tau)}{2}\Big{]}\,,\;\;\;\;\;\;\dot{\phi}=v\,,\;\;\;\;\;\;\;$
(10)

where both times $t$ and $\tau$ are normalized to the gap frequency, and the
current kernels are expressed through the Bessel functions:

$\displaystyle
i_{p}(\tau)=J_{0}(r_{1}\tau)Y_{0}(r_{2}\tau)+J_{0}(r_{2}\tau)Y_{0}(r_{1}\tau)\,,$
$\displaystyle
i_{qp}(\tau)=J_{1}(r_{1}\tau)Y_{1}(r_{2}\tau)+J_{1}(r_{2}\tau)Y_{1}(r_{1}\tau)\,.$
(11)

Here $r_{j}$ are the partial gaps: $r_{j}=\Delta_{j}/(\Delta_{1}+\Delta_{2})$,
$j=1,2$, i.e., $r_{1}+r_{2}=1$.

Conventions adopted above dictate that the current in the junction is
normalized to

$I_{N}=\frac{\pi\Delta_{1}\Delta_{2}}{2e(\Delta_{1}+\Delta_{2})R_{N}}\,,$ (12)

e.g., $i=I/I_{N}$. The normalization current $I_{N}$ is different from the
junction critical current $I_{C}$, which is expressed in terms of the complete
elliptic integral of the first kind $K(m)$ as (see, e.g., [23, 24])

$I_{C}=\int_{0}^{\infty}dtI_{p}(t)=\frac{2\Delta_{1}\Delta_{2}}{e(\Delta_{1}+\Delta_{2})R_{N}}K\Big{(}\frac{|\Delta_{1}-\Delta_{2}|}{\Delta_{1}+\Delta_{2}}\Big{)}\,.$
(13)

For completeness, one should note that the elliptic integral in Eq. (13) is
defined as in Ref. [27], not as in Ref. [26], in which case the argument of
$K$ should be squared. In the case of equal gaps, the relation between the two
currents simplifies to $I_{N}=I_{C}/2$.

One of the main qualitative features of the current in a Josephson tunnel
junction are the singularities (peaks and jumps) of the current components
resulting from the singularity of the density of state in the superconducting
electrodes at the edge of the energy gap. In the time domain, this gap-edge
singularity manifests itself in a very slow decay of the quasiparticle and the
pair current kernels. Using the standard expressions for the leading
asymptotic behavior of the Bessel functions at large arguments, we see that at
large time, $\tau\gg 1$, the current kernels (11) simplify to:

$-i_{p}(\tau)=i_{qp}(\tau)=\frac{2}{\pi\sqrt{r_{1}r_{2}}}\frac{\cos\tau}{\tau}\,.$
(14)

Although this expression is derived for the current kernels at vanishing
temperature, the finite temperature $T$ does not change the slow-decaying
nature of this asymptotic behavior. Using Eq. (4) for the finite-temperature
regime, one can show that the amplitude of the asymptotic (14) decreases as
$\tanh\Delta/T$ as a function of temperature $T$, while its time dependence
remain the same. This means that the problems of divergence resulting from the
slow decay of the current kernels are not cured by finite temperatures. This
conclusion is qualitatively consistent with the fact that the long-time
asymptotic behavior (14) of the current kernels is produced by the singularity
of the superconducting density of states at energy $\Delta$ which is not
affected by the temperature. Not too close to the superconducting transition,
the temperature changes only the occupation of the states at the gap edge,
which leads to the decrease of the overall magnitude of the asymptotics (14)
without affecting its dependence on time.

Expression (14) for the long-time asymptotics of the current kernels shows
that the main integral (1) for the tunnel current is not absolutely
convergent. When the evolution of the Josephson phase is such that the
function $\cos\tau$ in the asymptotics (1) taken together with the
corresponding functions of the Josephson phase in Eq. (1) has a non-vanishing
average over the period of the gap oscillations, the tunnel current diverges
logarithmically because of the slow, $1/\tau$, decay of the current kernels
(14) with time. At fixed bias voltage, this leads to the logarithmic “Riedel”
singularity [35] of the supercurrent amplitude at the gap voltage. Of course,
in real Josephson junctions, this singularity is smeared by several mechanisms
(see, e.g., [36]), the two most basic one being small static variations of the
energy gap over the area of the junction, and electron-phonon relaxation
producing finite lifetime of the quasiparticle energy states in the
superconducting electrodes of the junction. In junctions with finite electron
transparency $D$ (assumed to be vanishingly small in tunnel junctions), finite
$D$ itself directly smears the Riedel singularity [37]. Smearing of any kind
turns Riedel singularity into a finite “Riedel peak” of the supercurrent,
which was actively studied since the early days of investigations of the
Josephson effect – see, e.g., [38, 39, 40, 41, 42, 43]. These investigations
showed that the typical magnitude of the energy gap smearing in tunnel
junctions is quite small, on the order of 0.01 of the energy gap itself.
Different gap smearing mechanisms produce different specific time dependence
of the decoherence factors for the gap oscillations. For instance, static
variations of the gap over the junction area results in the Gaussian decay of
the gap oscillations, while the finite lifetime of the quasiparticle states,
in the most basic approximation, leads to the exponential decay. For realistic
weak smearing of the gap, however, precise nature of the smearing is not
important. In time domain, weak smearing affects only the long-time asymptotic
behavior (14) of the current kernels. It suppresses the logarithmic divergence
of the current by effectively limiting the integration range in Eq. (1) for
the tunnel current to some finite time $t_{0}$. As a result, only this cut-off
time, and not the precise time profile of the cut-off is important in the
regime of weak gap smearing, when $t_{0}$ is much larger than the gap
oscillation period $2\pi/\Omega$.

Figure 1: Current-voltage characteristic of a current-biased Josephson tunnel
junction between two BCS superconductors. Parameters $\beta$ and $w$ are
defined by Eqs. (9) and (15), respectively.

This means that in order to obtain a well-defined tunnel current $I_{T}(t)$,
and from this, a consistent description of the Josephson dynamics of a tunnel
junction based on Eqs. (8), (10, and (11), one needs to include in the model a
certain smearing mechanism of the energy gap which, in time domain, ensures a
faster decay of the current kernels than in the ideal BCS regime described by
Eq. (14). Numerical simulations discussed below assume a model of small
Gaussian fluctuations of the gap across the junction area, which produces an
extra Gaussian factor multiplying the current kernels (2) and improving
convergence of the current integrals (1). When the kernels are expressed
through the dimensionless time $\tau$, as in Eq. (11), the pair kernel is
transformed then like this:

$i_{p}(\tau)\;\;\;\rightarrow\;\;\;e^{-w^{2}\tau^{2}}i_{p}(\tau),$ (15)

with the similar transformation of the quasiparticle current kernel
$i_{qp}(\tau)$. Parameter $w$ in Eq. (15) characterizes the gap smearing, and
effectively cuts off the current kernels at time $t_{0}\simeq 1/(\Omega w)$.
As follows from the experimental results [38, 39, 40, 41, 42, 43], one can use
an estimate $w\simeq 0.01$, although the gap smearing strength varies for
different junction materials and structures. Such a Gaussian suppression
factor for the gap oscillations in time-domain implies that the magnitude of
the total energy gap of the junction exhibits small Gaussian fluctuations
around the mean value $\Delta_{1}+\Delta_{2}$ of relative magnitude $w$.

Figure 2: Current-voltage characteristic of a current-biased Josephson tunnel
junction with negligible junction capacitance and equal energy gaps of the
electrodes.

I now discuss the results for the dc current-voltage characteristics (“$IV$
curves”) of a current-biased Josephson tunnel junction obtained by the direct
numerical solution of the integro-differential equations (10) with the current
kernels (11) regularized according to Eq. (15). Although finite shunt
resistance $R_{S}$ does not make the simulation of the junction dynamics any
more difficult, to avoid extra parameters that are not relevant to our main
purpose, only unshunted junctions are considered below. Figure 1 shows an
example of an $IV$ curve calculated for the typical value of the broadening
parameter, $w=0.01$, and a moderate junction capacitance $\beta=0.3$. The main
feature of the $IV$ curve is that it exhibits two hysteresis loops: one at
voltages below the gap voltage $V_{g}=(\Delta_{1}+\Delta_{2})/e$, and one –
above the gap. The hysteresis at $V<V_{g}$ is well-understood in terms of the
inertia introduced into the phase dynamics by the junction electric
capacitance $C$, and is viewed as the characteristic feature of the tunnel-
junction behavior in the Josephson dynamics. Indeed, for large capacitance and
weak dissipation characteristic for tunnel junctions, the bias current needs
to be reduced well below the critical current in order to stop the time
evolution of the phase. In contrast to the inertial hysteresis, to the best of
author’s knowledge, the above-the-gap hysteresis has not been analyzed or
discussed previously, despite some indirect evidence in the literature [44,
45].

Figure 3: Current-voltage characteristics of a current-biased Josephson tunnel
junction: (a) with suppressed supercurrent: $i_{p}(\tau)\rightarrow
0.1i_{p}(\tau)$, and (b) with increased broadening of the gap singularity.

Figure 1 presents an example of the $IV$ curve of the “asymmetric” junction
with different energy gaps, $\Delta_{1}=0.5\Delta_{2}$. By making similar
calculations for other gap ratios one can see that the difference of the
energy gaps of the two junction electrodes does not have qualitative effects
on the junction current-voltage characteristics. To simplify the discussion of
the results, in what follows, I treat numerically only the case of equal gaps,
$\Delta_{1}=\Delta_{2}\equiv\Delta$. Consistent with its inertial mechanism,
the hysteresis at $V<V_{g}$ in the $IV$ curves becomes larger with increasing
junction capacitance, i.e., increasing $\beta$. As can be seen from Fig. 2,
which shows the $IV$ curve for vanishing junction capacitance, $\beta=0$, the
hysteresis at $V\gtrsim V_{g}$ becomes, by contrast, more pronounced with
decreasing capacitance $C$. With increasing $\beta$, this hysteresis shrinks,
and one can see that for equal gaps and $w=0.01$, it disappears for
$\beta>0.6$. Inertial hysteresis at $V<V_{g}$ is indeed smaller in Fig. 2, but
does not vanishes completely, as might have been expected for vanishing
capacitance. The most probable reason for this is the reactive component of
the tunnel current that results in a complex time evolution of the phase,
which gives a non-vanishing effective correction to the junction capacitance.
To understand the origin of the hysteresis at $V\gtrsim V_{g}$, one can note
from Fig. 2 that one, “ vertical”, side of the hysteresis is very close the
gap voltage, $V\simeq V_{g}$. This suggests that the hysteresis is related to
the Riedel singularity of the supercurrent in a tunnel junction. This
consideration is supported by the $IV$ curves calculated for the junction
dynamics with reduced pair component of the tunnel current (something that can
be achieved in principle, e.g., with the help of a magnetic field) or
increased broadening parameter $w$. As one can see from Fig. 3, in both cases,
above-the-gap hysteresis is indeed suppressed. Reduction of the pair component
of the current leads also to a smaller inertial hysteresis because of the
reduced junction critical current. All this makes the junction $IV$ curve
shown in Fig. 3a very close to just the dependence of the quasiparticle
component of the tunnel current on voltage in the regime of the fixed dc bias
voltage $V$. Extra (but still small) broadening of the gap singularity
illustrated in Fig. 3b does not affect the inertial hysteresis while reducing
the hysteresis at $V\gtrsim V_{g}$. This confirms that the latter hysteresis
has the same origin as the Riedel peak of the supercurrent at fixed bias
voltage.

To develop a more detailed understanding of this hysteresis, it is helpful to
look at the time evolution of the Josephson phase $\phi(t)$ on the vertical
branch of the hysteresis, $V\simeq V_{g}$. An example of this evolution is
shown in Fig. 4 which plots $\phi(t)$ for the junction with vanishing
capacitance, $\beta=0$, biased on the vertical branch of the hysteresis:
$I=2.0I_{C}$, $V\simeq 2\Delta/e$. One sees from this plot that the phase
evolves in a step-like manner, with the junction spending large fraction of
the Josephson oscillation period with $\phi(t)\simeq(\pi/2)\mbox{mod}(2\pi)$.
Since $V\simeq V_{g}$ in this regime, the Josephson frequency
$\omega_{J}=2eV/\hbar$ is simply related to the gap oscillations frequency,
$\omega_{J}=2\Omega$, and as shown below, the integral of the asymptotics of
the pair current indeed gives a non-vanishing contribution to the dc tunnel
current.

Figure 4: Typical time evolution of the Josephson phase across the current-
biased Josephson tunnel junction on the vertical branch of the above-the-gap
hysteresis. The plotted phase evolution corresponds to the point $I=2.0I_{C}$,
$V\simeq 2\Delta/e$ on the $IV$ curve shown in Fig. 2. Note that the time $t$
in this plot is normalized to the total gap frequency, i.e., the period of the
gap oscillations is $2\pi$. The period of the Josephson oscillations, which
have the frequency that is twice large than the gap frequency, is $\pi$,
consistent with the plot.

In more details, we assume the simplest model of the phase evolution
qualitatively consistent with the actual staircase-like time dependence shown
in Fig. 4:

$\phi(t)=\frac{\pi}{2}+2\pi n(t),$ (16)

where $n(t)$ is an integer that increases by 1 with each period
$2\pi/\omega_{J}$ of the Josephson oscillations. For such time dependence of
the phase, $\sin\phi(t)\equiv 1$, while
$\sin[\phi(t)-\phi(t^{\prime})]/2\equiv 0$. This means that besides the Ohmic
component, only the pair current gives nonvanishing contribution to the total
tunnel current, and

$\sin\frac{\phi(t)+\phi(t^{\prime})}{2}=\sin\phi(t)\cos\frac{\phi(t)-\phi(t^{\prime})}{2}$

$=\cos\pi[n(t)+n(t^{\prime})]=(-1)^{n(t)+n(t^{\prime})}\,.$

Under the assumption that the Riedel singularity is smeared only weakly, the
pair component of the tunnel current is dominated in time domain by its large-
time asymptotics (14). In this approximation, expression for the pair current
in dimensional form is:

$I(t)=\frac{\sqrt{\Delta_{1}\Delta_{2}}}{eR_{N}}\int\frac{d\eta}{\eta}\cos\Omega\eta\sin\frac{\phi(t)+\phi(t-\eta)}{2}\,,$
(17)

where the lower limit of the integral is set by condition of validity of the
long-time asymptotics of the Bessel functions to $\eta\simeq\Omega^{-1}$,
while the upper limit is determined by the width of the Riedel singularity,
e.g., by Eq. (15): $\eta\simeq(w\Omega)^{-1}$. For small width of the Riedel
singularity, $w\ll 1$, the integral in Eq. (17) can be found in two steps.
First, one averages out the rapidly-varying terms of frequency $\Omega$ in Eq.
(17) by integrating over one period of the gap oscillations neglecting the
slow-varying factors, including the factor $1/\eta$. Then, one integrates the
slow-varying terms with frequencies on the order of $w\Omega$, with high-
frequency terms already averaged out.

On the vertical branch of the above-the-gap hysteresis, the deviation $\delta
V$ of the dc voltage $V$ across the junction from the gap value is small. As
will be seen soon, it is proportional to the width of the Riedel singularity,
$\delta V\equiv V-V_{g}\propto w$. This means that the deviation of the
Josephson frequency from $2\Omega$ is also small, and can be neglected when
evaluating the integral (17) for the pair current over one period of the gap
oscillations:

$I^{\prime}(t)=\frac{\sqrt{\Delta_{1}\Delta_{2}}}{eR_{N}}\,\frac{\Omega}{2\pi}\int_{0}^{2\pi/\Omega}d\eta\cos\Omega\eta\,(-1)^{n(t)+n(t-\eta)}$

$=\frac{2}{\pi}\frac{\sqrt{\Delta_{1}\Delta_{2}}}{eR_{N}}|\sin\Omega t|\;.$
(18)

In this expression, it is assumed that the time $t$ is counted starting at one
of the jumps of the staircase-like evolution of the Josephson phase (Fig. 4).
To get the full magnitude of the pair current for the phase (16) we now need
to take the second step and integrate the current (18) over the full range of
the gap oscillations in Eq. (17). In the approximation discussed above, the
corresponding integral is:

$\int_{\Omega^{-1}}^{\infty}\frac{d\eta}{\eta}\cos[(\Omega-\omega_{J}/2)\eta]e^{-w\Omega\eta}=\ln\big{[}w^{2}+(\delta
V/V_{g})^{2}\big{]}^{-1/2}\,.$

Combining these two results, we obtain an expression for the pair current on
the vertical side of the above-the-gap hysteresis in the approximation (16)
for $\phi(t)$:

$I^{\prime}(t)=\frac{2}{\pi}\frac{\sqrt{\Delta_{1}\Delta_{2}}}{eR_{N}}\ln\big{[}w^{2}+(\delta
V/V_{g})^{2}\big{]}^{-1/2}|\sin\Omega t|\;.$ (19)

Although, obviously, the assumption (16) of the sharp steps in the evolution
of the phase $\phi(t)$ and the resulting time dependence (19) of the pair
current do not provide an exact self-consistent solution of the evolution
equations (8), they describe the main qualitative feature of the smoother
self-consistent evolution obtained by direct numerical solution of Eqs. (8)
and shown in Fig. 4. The staircase-like Josephson phase in Fig. 4 produces a
finite dc component of the supercurrent in the junction, an estimate of which
is obtained from Eq. (19):

$I=\frac{8}{\pi^{2}}I_{R}\,,\;\;I_{R}\equiv\frac{\sqrt{\Delta_{1}\Delta_{2}}}{2eR_{N}}\,\ln[w^{2}+(\delta
V/V_{g})^{2}]^{-1/2}.$ (20)

By contrast, for fixed bias voltage $V$ across the junction, the regime that
approximately describes the higher-voltage branch of the above-the gap
hysteresis, the pair current does not have a dc component. If the phase
evolves simply as $\phi(t)=\phi_{0}+\omega_{J}t$, the non-vanishing
contribution to the pair current Eq. (17) comes from the term:

$\sin\frac{\phi(t)+\phi(t^{\prime})}{2}=\sin\phi(t)\cos\frac{\omega_{J}(t^{\prime}-t)}{2}\,.$

The same calculation as done above for the phase (16) shows then that the
current (17) in this regime is:

$I^{\prime}(t)=I_{R}\sin\phi(t)\,,$ (21)

and indeed averages to zero over time $t$. This also shows that the current
$I_{R}$ introduced in Eq. (20) has the meaning of the supercurrent amplitude
in the junction voltage-biased in the vicinity of the Riedel singularity.

In summary, the difference between the vanishing and finite dc components of
the pair current for the “fixed-voltage” linear evolution of the phase on one
branch of the hysteresis, and the staircase-like evolution (Fig. 4) on the
other branch, explains the mechanism of the above-the-gap hysteresis in the
current-bias Josephson tunnel junctions. This means that the hysteresis is
closely related to the logarithmic gap singularity of the supercurrent in
tunnel junctions. Since the magnitude of the dc component (20) of the
supercurrent varies with the dc voltage $V$ across the junction on the scale
$wV_{g}$, for small smearing $w$ of the gap singularity, the side of the
hysteresis near the gap voltage looks nearly vertical (Fig. 2), with the
voltage $V$ pinned down to the gap voltage $V_{g}$. Also, from the point of
view of this explanation, the fact that the hysteresis is suppressed with
increasing junction capacitance is natural, since qualitatively, the
capacitance makes the phase evolution more linear, suppressing the steps
responsible for the dc part of the supercurrent.

##  IV Conclusion

In summary, a novel expression for the time-domain current kernels of a
Josephson tunnel junctions between two ideal BCS superconductors with in
general different energy gaps is derived and used to simulate the dynamics of
the current-biased junctions. The simulations make it possible to study in
detail the current-voltage characteristics of such junctions and reveal their
qualitatively new features, the main one – an additional hysteresis at
voltages above the gap voltage of the junction. The hysteresis is related to
the Riedel singularity of the junction supercurrent and is explained in terms
of the step-like time dependence of the Josephson oscillations at the gap
voltage, that can be viewed as characteristic feature of the resonance between
the Josephson oscillations and the gap oscillations with frequency
$\Omega=(\Delta_{1}+\Delta_{2})/\hbar$ in the pair current kernel of the
junction, where $\Delta$’s are the energy gaps of the two junction electrodes.
This regime of the gap resonance exists in the junctions with the relatively
short $RC$ time constant, $\beta\equiv\Omega R_{N}C\lesssim 0.6$, where
$R_{N}$ is the normal-state resistance, and $C$ – capacitance of the junction.
In typical Josephson tunnel junctions, $RC$ time constant is nearly
independent of the junction area, and characterizes the “critical current
density” $j$. Since $\beta$ is proportional to the gap frequency, the
condition that $\beta$ is not too large is more difficult to fulfill in
junctions employing superconductors with large energy gaps, like niobium that
is used frequently to produce the large-scale circuits of superconductor
electronics. Aluminum junctions, which are also quite standard, have
significantly smaller gap $\Delta\simeq 0.2$ meV, so that the gap frequency is
$\Omega\simeq 0.6$ THz and reaching the regime of the gap resonance,
$\beta\lesssim 0.6$ requires junctions with the time constant $R_{N}C$ about 1
ps. High critical current density can be reached in tunnel junctions with
insulator layer thickness below 1 nm. Combined with the relative dielectric
constant $\epsilon\simeq 8$ of the standard aluminum oxide tunnel barrier,
this gives the electric capacitance per unit area $C\simeq 7\cdot 10^{-14}$
F/$\mu$m2, and $R_{N}C\simeq 1$ ps for the resistance of unit area $15$
Ohm$\cdot\mu$m2. This translates into a critical current density $j\simeq 2$
kA/cm2 that should be achievable in junctions which still preserve both the
high-quality tunnel barriers and superconducting electrodes with sharp edges
of the energy gap.

The author would like to thank Alexander Kirichenko and Pavel Shevchenko for
sharing the unpublished results. This work was supported by the IARPA
Supertools program through the Synopsys and Hypres subcontracts. The author
would like to acknowledge Timur Fillipov, Deep Gupta, Amol Inamdar, Dmitri
Kirichenko, Michail Kupriyuanov, Konstantin Likharev, Vasili Semenov, and
Stephen Whiteley for the stimulating discussions of the results.

## Appendix A

This appendix presents the derivation of the expressions for the current
kernels discussed in Sec. I. Since the calculation of the tunnel current in a
Josephson junction is described in details in the literature, only the steps
leading specifically to Eqs. (2) are outlined.

Using the standard method of the tunnel Hamiltonian, the current $I_{T}(t)$ in
a Josephson tunnel junction can be expressed as in Eq. (1) with the current
kernels determined as follows. First, the standard evaluation of the pair
creation/annihalation terms in the two junction electrodes (see, e.g., [23])
defines the pair-current kernel:

$\displaystyle I_{p}(t)=\frac{2\Delta_{1}\Delta_{2}}{\pi e\hbar
R_{N}}\mbox{Im}\big{[}F_{1}(t)F_{2}(t)\big{]}\,,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$
$\displaystyle F_{j}(t)=\\!\int\frac{d\xi}{\xi}f(\xi)\rho_{j}(\xi)e^{-i\xi
t/\hbar}\\!,\rho_{j}(\xi)=\frac{|\xi|\theta(\xi^{2}-\Delta_{j}^{2})}{\sqrt{\xi^{2}-\Delta_{j}^{2}}}.\;\;\;\;\;\;\;\;$
(22)

Here $\xi$ is the energy of the quasiparticle states of the BCS Hamiltonian,
$f(\xi)$ is the Fermi distribution function at temperature $T$, and
$\rho_{j}(\xi)$ is the superconducting density of states in the $j$th
electrode. Taking into account in Eq. (22) that

$\mbox{Im}\big{[}e^{-i\xi t/\hbar}e^{-i\zeta t/\hbar}\big{]}=-\sin\frac{\xi
t}{\hbar}\cos\frac{\zeta t}{\hbar}-\cos\frac{\xi t}{\hbar}\sin\frac{\zeta
t}{\hbar}\,,$

and

$f(\xi)+f(-\xi)=1\,,\;\;\;f(\xi)-f(-\xi)=-\tanh(\xi/2T)\,,$

one directly obtains Eqs. (2) and (4) for $I_{p}(t)$.

Similarly to Eq. (22) for the pair kernel, equation for the quasiparticle
current kernel $I_{qp}(t)$ is:

$\displaystyle I_{qp}(t)=\frac{2}{\pi e\hbar
R_{N}}\mbox{Im}\big{[}G_{1}(t)G_{2}(t)\big{]}\,,$ $\displaystyle G_{j}(t)=\int
d\xi f(\xi)\rho_{j}(\xi)e^{-i\xi t/\hbar}.$ (23)

The same transformations of Eq. (23) as for the Eq. (22) would give:

$\displaystyle I_{qp}(t)=\frac{2}{\pi e\hbar R_{N}}\int_{0}^{\infty}d\xi
d\zeta\rho_{1}(\xi)\rho_{2}(\zeta)\;\;\;\;$
$\displaystyle\cdot\Big{(}\sin\frac{\xi t}{\hbar}\cos\frac{\zeta
t}{\hbar}\tanh\frac{\xi}{2T}+\cos\frac{\xi t}{\hbar}\sin\frac{\zeta
t}{\hbar}\tanh\frac{\zeta}{2T}\Big{)},\;\;\;\;\;$ (24)

and reproduce directly the integral part of the quasiparticle kernel (2) in
terms of the function $J_{1}(z)$ and $A_{1}(z,\beta)$. This conclusion,
however, misses the normal-state Ohmic contribution to the current determined
by the singular part of the integrals for $J_{1}(z)$ and $A_{1}(z,\beta)$
which are not properly convergent. A more accurate transformation of the
expression for $I_{qp}$ requires to explicitly separate the singular part of
these integrals, something that can be achieved by subtracting the normal
density of states [equal to 1 in the conventions of Eqs. (22) and (23)] from
the superconducting ones:

$\rho_{1}\rho_{2}=(\rho_{1}-1)(\rho_{2}-1)+(\rho_{1}-1)+(\rho_{2}-1)+1$ (25)

in Eq. (23).

Then, one needs to do the calculation similar to that for the pair current
kernel (22) individually for each term in the expansion (25). For the last,
normal-state term, repeating the steps similar to those in Eq. (6), one
obtains:

$G_{N}(t)\equiv\int d\xi f(\xi)e^{-i\xi t/\hbar}=\frac{i\pi T}{\sinh[(\pi
Tt/\hbar)+i0]}\,.$ (26)

With this result, Eq. (23) shows that the normal-state contribution to the
quasiparticle current kernel is:

$\displaystyle I_{N}(t)=\frac{2}{\pi e\hbar
R_{N}}\mbox{Im}\big{[}G_{N}(t)\big{]}^{2}$ $\displaystyle=-\frac{2\pi
T^{2}}{e\hbar R_{N}}\mbox{Im}\frac{1}{\sinh^{2}[(\pi Tt/\hbar)+i0]}\,,$
$\displaystyle=\frac{2T}{eR_{N}}\frac{d}{dt}\mbox{Im}\coth[(\pi
Tt/\hbar)+i0]=-\frac{2\hbar}{eR_{N}}\delta^{\prime}(t).\;\;\;\;$ (27)

This expression can be transformed into the corresponding part of Eq. (2). To
provide the context for this transformation, it should be mentioned that all
current kernels discussed in this work are obtained by the perturbation theory
in tunneling which leads to Eq. (1) for the current. Although the structure of
the perturbation theory itself ensures through Eq. (1) that the current
response to the bias voltage has the causal structure, the perturbative
current kernels do not have this property. As can be seen from the integral
form of the current kernel components (3) and (4) (not the final form in terms
of the Bessel functions) the kernels $I_{p}(t)$, $I_{qp}(t)$ have the property

$I_{p}(-t)=-I_{p}(t)\,,\;\;\;\;\;I_{qp}(-t)=-I_{qp}(t)\,.$

They do not vanish for $t<0$ as would be natural for the causal response, and
as indeed happens in the non-perturbative calculations (see, e.g., [28]). For
the regularized “superconducting” part of the current kernels this is a formal
property, since only the $t>0$ parts enter expression for the current. The
singular nature of the normal part (27), however, makes it necessary to
resolve the singularity using the relation

$2\int^{\infty}_{0}dt\delta(t)=1=\int^{\infty}_{0}dt\delta(t-0)\,.$

Equation (2) is written in the form that avoids this discussion.

The next term in Eq. (25) one should deal with is $\rho-1$, i.e.

$\displaystyle\int d\xi f(\xi)[\rho(\xi)-1]e^{-i\xi
t/\hbar}=\int_{0}^{\infty}d\xi[\rho(\xi)-1](\cos\frac{\xi t}{\hbar}$
$\displaystyle+i\sin\frac{\xi t}{\hbar}\tanh\frac{\xi}{2T})\,.$

Written in dimensionless form ($x=\xi/\Delta$, $z=\Delta t/\hbar$), its real
part is:

$\Delta\int_{0}^{\infty}dx\Big{[}\frac{x\theta(x-1)}{\sqrt{x^{2}-1}}-1\Big{]}\cos
zx\,,$ (28)

and is absolutely convergent. The second term in the brackets in this
expression affects the value of the integral only at $z=0$:

$\int_{0}^{\infty}dx\cos zxe^{-\epsilon x}\Big{|}_{\epsilon\rightarrow
0}=\pi\delta(z)\,,$

while the first term is the derivative of the integral (3) for $J_{0}(z)$.
This gives us the following expression for $J_{1}(z)$ for $z>0$:

$J_{1}(z)=-\frac{2}{\pi}\int_{1}^{\infty}dx\frac{x\cos zx}{\sqrt{x^{2}-1}}\,.$

For $z=0$, the value of $J_{1}(0)=0$ is determined by the full integral (28)
with removed divergence.

Making use of Eq. (6), one can express the dimensionless form of the imaginary
part of the $\rho-1$ term as

$\int_{0}^{\infty}dx\Big{[}\frac{x\theta(x-1)}{\sqrt{x^{2}-1}}-1\Big{]}\sin
zx\tanh\beta x=-\frac{\pi}{2}(A_{1}+S)\,,$ (29)

where $A_{1}$ is defined in Eq. (4) and

$S=\frac{1}{\beta\sinh(\pi Tt/\hbar)}\,.$

Since the total integral (29) is absolutely convergent, we know that
$(A_{1}+S)\Big{|}_{z\rightarrow 0}=0$, i.e. $A_{1}$ diverges in the same way
as $S$ for small $t$.

Combining the results for the real and imaginary parts of the $\rho-1$ term we
get:

$\int d\xi f(\xi)[\rho(\xi)-1]e^{-i\xi
t/\hbar}=-\frac{\pi\Delta}{2}[J_{1}(z)+i(A_{1}+S)]\,.$ (30)

Also, Eq. (26) for $G_{N}(t)$ can be written as

$G_{N}(t)=\frac{\pi\Delta}{2}\Big{[}\frac{2\hbar}{\Delta}\delta(t)+iS(t)\Big{]}\,.$
(31)

Using Eqs. (30) and (31) to evaluate the contributions to the quasiparticle
current kernel (23) of all three terms (besides the normal-state term) in the
expansion (25), one gets the following total:

$\displaystyle\frac{\pi^{2}\Delta_{1}\Delta_{2}}{4}\mbox{Im}\Big{\\{}[J_{1}+i(A_{1}+S)]\big{|}_{\Delta_{1}}[J_{1}+i(A_{1}+S)]\big{|}_{\Delta_{2}}$
$\displaystyle-[J_{1}+i(A_{1}+S)]\big{|}_{\Delta_{1}}[\frac{2\hbar}{\Delta}\delta(t)+iS]\big{|}_{\Delta_{2}}$
$\displaystyle-[J_{1}+i(A_{1}+S)]\big{|}_{\Delta_{2}}[\frac{2\hbar}{\Delta_{1}}\delta(t)+iS]\big{|}_{\Delta_{1}}\Big{\\}}$
$=\frac{\pi^{2}\Delta_{1}\Delta_{2}}{4}\Big{\\{}J_{1}\big{(}\frac{\Delta_{1}t}{\hbar}\big{)}A_{1}\big{(}\frac{\Delta_{2}t}{\hbar},\beta_{2}\big{)}+J_{1}\big{(}\frac{\Delta_{2}t}{\hbar}\big{)}A_{1}\big{(}\frac{\Delta_{1}t}{\hbar},\beta_{1}\big{)}\Big{\\}},$

where in the last line we took into account that $(A_{1}+S)=0$ for $t=0$
ensuring that the products $(A_{1}+S)\delta(t)$ vanish. Together with the
normal-state contribution (27) this equation gives the quasiparticle current
kernel (2).

## Appendix B

The main element of the numerical procedure used to simulate the Josephson
dynamics of the current-biased tunnel junctions in Sec. III is the calculation
of the Bessel functions that determine the current kernels (5), and through
them, the time-dependent current response of the junction. The four Bessel
functions that enter the expressions (5) for the current kernels were
calculated from the standard expressions. For small values of the
dimensionless argument $x$, one can use the power-series expansions (see,
e.g., [26])

$J_{n}(x)=\sum_{m=0}^{\infty}\frac{(-1)^{m}}{m!(n+m)!}\Big{(}\frac{x}{2}\Big{)}^{2m+n}\,,$
(32)

The current kernels (5), include only the first two functions, $n=0,1$. In
practice, expansion (32) gives converging results for $x\lesssim 30$. Keeping
up to $m=100$ terms in the series one can obtain $J_{0,1}$ with accuracy
better that $10^{-4}$ in this range of the argument, as one can check either
by comparison to a different, asymptotic expansion of these functions
discussed below, or, by comparison with the numerical tables in [26] for the
values of $x$ for which the tables are available.

Figure 5: Two Bessel functions of the first kind calculated from Eqs. (32) for
$x\leq 30$ and from Eqs. (34) and (35) for $x>30$. The inserts show that even
on a fine scale (on the order of $10^{-2}$), the two expansions, power series
at small arguments and the asymptotic expansions at large arguments, agree
very well, without a noticeable transition.

For the Bessel functions $Y_{n}(x)$ of the second kind, the power-series
expansion similar to Eq. (32) is

$\displaystyle
Y_{n}(x)=\frac{2}{\pi}J_{n}(x)\Big{(}\ln\frac{x}{2}+\gamma\Big{)}-\frac{1}{\pi}\sum_{m=0}^{n-1}\frac{(n-1-m)!}{m!}\;\;\;\;\;\;\;$
$\displaystyle\cdot\Big{(}\frac{x}{2}\Big{)}^{2m-n}\\!\\!-\frac{1}{\pi}\sum_{m=0}^{\infty}\frac{(-1)^{m}(x/2)^{2m+n}}{m!(n+m)!}\Big{(}\sum_{k=1}^{n+m}\\!\frac{1}{k}+\sum_{k=1}^{m}\\!\frac{1}{k}\Big{)},\;\;\;\;\;\;\;\;$
(33)

where $\gamma=0.5772156649$ is Euler’s constant. The discussion of Eq. (32)
above applies to Eq. (33) as well. The power series of up to $m=100$ terms
converges for $x\lesssim 30$ and gives $Y_{0,1}$ with accuracy better that
$10^{-4}$ in this range, as one can check by comparison to the asymptotic
expansion of these functions.

Figure 6: Two Bessel functions of the second kind calculated from Eqs. (33)
for $x\leq 30$ and from Eqs. (34) and (35) for $x>30$. Similarly to Fig. 5,
the inserts show that even on a fine scale on the order of $10^{-2}$, the two
methods agree very well, without a noticeable transition.

Very slow convergence of the Bessel functions, which decay at large argument
$x$ only as $x^{-1/2}$, implies that the integrals in Eq. (1) that determines
the tunnel current do not converge by themselves under some conditions. For
the realistic strength of the gap-smearing mechanisms which ensure convergence
of these integrals, one needs to integrate the current kernels over the time
range that extend to $x\sim 100$ in dimensionless units. This means that the
power-series expansions for the Bessel functions discussed above are not
sufficient, and one needs to use also the asymptotic expansions for large $x$.
Moreover, to ensure that the asymptotic expansions match with sufficient
accuracy the power-series expansions in the range $10\lesssim x\lesssim 30$,
where they are both applicable, one needs to go beyond the leading terms in
the expansions. Keeping three terms should give the asymptotics with the
accuracy on the order of $10^{-5}$ in the relevant range of the argument $x$.
Then, for the four functions that enter the current kernels (5), the
appropriate asymptotic expansions are [29]:

$\displaystyle J_{0}(x)=\Big{(}\frac{2}{\pi
x}\Big{)}^{1/2}\big{\\{}P\cos\big{(}x-\frac{\pi}{4}\big{)}+Q\sin\big{(}x-\frac{\pi}{4}\big{)}\big{\\}},$
$\displaystyle Y_{0}(x)=\Big{(}\frac{2}{\pi
x}\Big{)}^{1/2}\big{\\{}P\sin\big{(}x-\frac{\pi}{4}\big{)}-Q\cos\big{(}x-\frac{\pi}{4}\big{)}\big{\\}},\;\;\;\;\;\;$
(34) $\displaystyle
P=1-\frac{9}{128x^{2}}+\frac{3675}{2^{15}x^{4}}\,,\;\;\;\;Q=\frac{1}{8x}-\frac{75}{1024x^{3}}\,,\;\;\;\;\;\;\;\;$

and

$\displaystyle J_{1}(x)=\Big{(}\frac{2}{\pi
x}\Big{)}^{1/2}\big{\\{}P^{\prime}\cos\big{(}x-\frac{3\pi}{4}\big{)}+Q^{\prime}\sin\big{(}x-\frac{3\pi}{4}\big{)}\big{\\}},\;\;\;\;\;\;$
$\displaystyle Y_{1}(x)=\Big{(}\frac{2}{\pi
x}\Big{)}^{1/2}\big{\\{}P^{\prime}\sin\big{(}x-\frac{3\pi}{4}\big{)}-Q^{\prime}\cos\big{(}x-\frac{3\pi}{4}\big{)}\big{\\}},\;\;\;\;\;\;\;$
(35) $\displaystyle
P^{\prime}=1+\frac{15}{128x^{2}}-\frac{4725}{2^{15}x^{4}}\,,\;\;\;\;Q^{\prime}=-\frac{3}{8x}+\frac{105}{1024x^{3}}\,.\;\;\;\;\;\;\;\;$

Four Bessel functions that determine the current response (1) of a Josephson
tunnel junction between the two ideal BCS superconductor with in general
different energy gaps at zero temperature are shown in Figs. 5 and 6. The
functions are calculated from the power series (32)and (33) for small
arguments, $x\leq 30$, and from the asymptotic expansions (34) and (35) for
large arguments, $x>30$. The main point illustrated by these plots, explicitly
by the insets, is the agreement between the small-argument and large-argument
expansions with very high accuracy. This accuracy is better than $10^{-4}$, as
can roughly be estimated from these plots.

With the Bessel functions calculated accurately and in the arbitrary large
range of the arguments as described above, Eqs. (10) and (11) that govern the
Josephson dynamics of the current-biased junctions can be solved directly by
the simplest numerical procedures, integration of the current kernels at each
time step to find the tunnel current, and then propagation of the values of
the Josephson phase and the voltage actoss the junction according to the
evolution equations. The results of these simulations are presented in Sec.
III.

## References

  * [1] I. Giaever, “Energy gap in superconductors measured by electron tunneling”, Phys. Rev. Lett. 5, 147 (1960); “Electron tunneling between two superconductors”, ibid. 5, 464 (1960).
  * [2] M.H. Cohen, L.M. Falicov, and J.C. Phillips, “Superconductive tunneling”, Phys. Rev. Lett. 8, 316 (1962).
  * [3] B.D. Josephson, “Possible new effects in superconductive tunnelling”, Phys. Lett. 1 251 (1962).
  * [4] P.W. Anderson and J.M. Rowell, “Probable observation of the Josephson tunnel effect”, Phys. Rev. Lett. 10, 230 (1963).
  * [5] J. Nicol, S. Shapiro, and P.H. Smith“Direct measurement of the superconducting energy gap”, Phys. Rev. Lett. 5, 461 (1960).
  * [6] K.K. Likharev and V.K. Semenov, “RSFQ logic/memory family: a new Josephson junction technology for sub-terahertz-clock-frequency digital systems”, IEEE Trans. Appl. Supercond. 1, 3 (1991).
  * [7] S.K. Tolpygo, “Superconductor digital electronics: scalability and energy efficiency issues”, Low Temperature Physics 42, 361 (2016).
  * [8] A.J. Leggett and A. Garg, “Quantum mechanics versus macroscopic realism: is the flux there when nobody looks?”, Phys. Rev. Lett. 54, 857 (1985).
  * [9] D.V. Averin, A.B. Zorin, and K.K. Likharev, “Bloch oscillations in small Josephson junctions”, Sov. Phys. JETP 61, 407 (1985).
  * [10] Y. Nakamura, Yu.A. Pashkin, and J.S. Tsai, “Coherent control of macroscopic quantum states in a single-Cooper-pair box”, Nature 398, 786 (1999).
  * [11] J.R. Friedman,V. Patel,W. Chen, S.K. Tolpygo,and J.E. Lukens, “Quantum superposition of distinct macroscopic states”, Nature 406, 43 (2000).
  * [12] V.E. Manucharyan, J. Koch, L.I. Glazman and M.H. Devoret, “Fluxonium: Single Cooper-pair circuit free of charge offsets”, Science 326, 113 (2009).
  * [13] M.T. Bell, I.A. Sadovskyy, L.B. Ioffe, A.Yu. Kitaev, and M.E. Gershenson, “Quantum superinductor with tunable nonlinearity”, Phys. Rev. Lett. 109, 137003 (2012).
  * [14] I.M. Pop, B. Doucot, L. Ioffe, I. Protopopov, F. Lecocq, I. Matei, O. Buisson, and W. Guichard, “Experimental demonstration of Aharonov-Casher interference in a Josephson junction circuit”, Phys. Rev. B 85, 094503 (2012).
  * [15] J. Clarke and F.K. Wilhelm, “Superconducting quantum bits”, Nature 453, 1031 (2008).
  * [16] M.W. Johnson, M.H.S. Amin, S. Gildert, et al., “Quantum annealing with manufactured spins”, Nature 473, 194 (2011).
  * [17] M.H. Devoret and R.J. Schoelkopf, “Superconducting circuits for quantum information: an outlook”, Science 339, 1169 (2013).
  * [18] R. Barends, J. Kelly, A. Megrant, et al., “Superconducting quantum circuits at the surface code threshold for fault tolerance”, Nature 508, 500 (2014).
  * [19] F. Arute, K. Arya, R. Babbush, et al., “Quantum supremacy using a programmable superconducting processor”, Nature 473, 194 (2019).
  * [20] V.K. Semenov, G.V. Danilov, and D.V. Averin, “Negative-inductance SQUID as the basic element of reversible Josephson-junction circuits”, IEEE Trans. Appl. Supercond. 13, 938 (2003); Jie Ren, V.K. Semenov, Yu.V. Polyakov, D.V. Averin, and J.S. Tsai, “Progress towards reversible computing with nSQUID arrays”, ibid. 19, 961 (2009).
  * [21] N.R. Werthamer, “Nonlinear self-coupling of Josephson radiation in superconducting tunnel junctions”, Phys. Rev. 147, 255 (1966).
  * [22] A.I. Larkin and Yu. N. Ovchinnikov, “Tunnel effect between superconductors in an alternating field”, Zh. Eksp. Teor. Fiz. 51, 1535 (1966) [Sov. Phys. JETP 24, 1035 (1967)].
  * [23] A. Barone and G. Paterno, “Physics and applications of the Josephson effect” (John Wiley, 1982).
  * [24] K.K. Likharev, “Dynamics of Josephson junctions and circuits”, (Gordon and Breach, 1986).
  * [25] R.E. Harris, “Intrinsic response time of a Josephson tunnel junction”, Phys. Rev. B 13, 3818 (1976).
  * [26] M. Abramowitz and I.A.Stegun, Eds., “Handbook of mathematical functions”, (National Bureau of Standards, 1964), Ch. 9.
  * [27] A.P. Prudnikov, Yu.A. Brychkov, and O.I. Marichev, Integrals and series, vol. 2. Special functions, (Gordon and Breach, 1992) Sec. 2.13.19.
  * [28] D.V. Averin, G. Wang, and A.S. Vasenko, “Time-dependent Andreev reflection”, Phys. Rev. B 102, 144516 (2020).
  * [29] G.E. Andrews, R. Askey, and R. Roy, “Special functions” (Cambridge Univ. Press, 1999) Ch. 4.
  * [30] A.A. Odintsov, V.K. Semenov, and A.B. Zorin, “Specific problems of numerical analysis of the Josephson junction circuits”, IEEE Trans. Mag. 23, 763 (1987).
  * [31] D.R. Gulevich, V.P. Koshelets, and F.V. Kusmartsev, “Josephson flux-flow oscillator: The microscopic tunneling approach”, Phys. Rev. B 96, 024515 (2017).
  * [32] http://pscan2sim.org
  * [33] D.R. Gulevich, “MiTMoJCo: Microscopic tunneling model for Josephson contacts”, Comp. Phys. Commun. 251, 107091 (2020).
  * [34] A.F. Kirichenko and P.S. Shevchenko, unpublished (2018).
  * [35] E. Riedel, “Zum Tunneleffekt bei Supraleitem im Mikrowellenfeld”, Z. Naturforschg. 19A, 1634 (1964).
  * [36] A.B. Zorin, I.O. Kulik, K.K. Likharev, and J.R. Schrieffer, “The sign of the interference current component in superconducting tunnel junctions”, Sov. J. Low Temp. Phys. 5, 537 (1979).
  * [37] D.V. Averin and A. Bardas, “AC Josephson effect in a single quantum channel”, Phys. Rev. Lett. 75, 1831 (1995).
  * [38] C.A. Hamilton and S. Shapiro, “Experimental demonstration of the Riedel peak”, Phys. Rev. Lett. 26, 426 (1971).
  * [39] S.A. Buckner, T.F. Finnegan, and D.N. Langenberg, “Riedel singularity in Sn–Sn-Oxide–Sn Josephson tunnel junctions”, Phys. Rev. Lett. 28, 150 (1972).
  * [40] B. Kofoed and K. Saermark, “Observation of Riedel-peak in weak links”, Phys. Lett. A 55), 72 (1975).
  * [41] G. Vernet and R. Adde, “Direct observation of the Riedel singularity in a superconducting point contact”, Appl. Phys. Lett. 28, 559 (1976).
  * [42] S. Morita, S. Imai, I. Fukushi, S. Takaki, Y. Takeuti, and N. Mikoshiba, “Damping factor of the Riedel peak in bridge-type Josephson junctions” J. Phys. Soc. Jap. 52, 617 (1983).
  * [43] D. Winkler, “Measurements of the Riedel peak in superconducting aluminum tunnel junctions”, IEEE Trans. Appl. Supercond. 3, 2234 (1993).
  * [44] D.G. McDonald, E.G. Johnson, and R.E. Harris, “Modeling Josephson junctions”, Phys. Rev. B 13, 1028 (1976).
  * [45] A.B. Zorin, K.K. Likharev, and S.I. Turovets, “Dynamics of Josephson tunnel junctions with a finite-width Riedel peak”, IEEE Trans. Mag. 19, 629 (1983).