#  $J/\Psi$ and $\chi_{cJ}(J=0,1)$ production in electron-positron
annihilation at $\sqrt{s}=10.6$ GeV in the framework of Bethe-Salpeter
equation

Shashank Bhatnagar, Teresa Aruja, Vaishali Guleria

###### Abstract

In present work we study the production of ground and excited charmonium pairs
in $e^{-}e^{+}\rightarrow\Psi(nS)+\chi_{cJ}(nP)$ for $J=0,1$ and $n=1,2$,
through leading order (LO) tree-level diagrams $\sim
O(\alpha_{em}\alpha_{s})$, which proceed through exchange of a virtual photon
and an internal gluon line connecting two quark lines (in the triangle quark
loop part of the diagram), in the framework of $4\times 4$ Bethe-Salpeter
equation, at center of mass energy, $\sqrt{s}=10.6GeV.$ We have tried to show
that the cross sections for this process calculated with use of the four
leading order $\sim O(\alpha_{em}\alpha_{s})$ diagrams using BSE approach come
close to experimental data, which might be mainly due to the consistent
treatment of motion of quarks inside the hadrons in the BSE framework.

Department of Physics, University Institute of Sciences, Chandigarh
University, Mohali-140413, India

Key words: Bethe-Salpeter equation, double charmonium production, cross
sections

##  1\. Introduction

In recent years there has been observation of double charm quark pair
production by BABAR and Belle Collaborations [1, 2, 3] in $e^{-}e^{+}$
annihilation at center of mass energies, $\sqrt{s}=10.6$ GeV. This process is
particularly interesting in testing quarkonium production mechanism. This
process proceeds through the exchange of a single photon. Charge conjugation
symmetry requires that if one of the charmonia is $1^{--}$ (such as $J/|Psi$),
then the other should have $C=+$ such as $0^{++}(\chi_{c0})$,
$1^{++}(\chi_{c1})$, $2^{++}(\chi_{c2})$, or $0^{-+}(\eta_{c})$. In this work
we will focus on the two processes: $e^{-}+e^{+}\rightarrow J/\Psi+\chi_{c0}$,
and $e^{-}+e^{+}\rightarrow J/\Psi+\chi_{c1}$.

But there is a significant discrepancy between the experimental measurements
of total cross sections of $e^{-}e^{+}\rightarrow J/\Psi\eta_{c}$ and
$e^{-}e^{+}\rightarrow J/\Psi\chi_{cJ}$ at energy, $\sqrt{s}=10.6$GeV., and
their NRQCD predictions using leading order (LO) diagrams. It is to be
mentioned that the LO diagrams alone give cross-sections which are much less
than data. Such a large discrepancy between experimental results and
theoretical predictions has been a challenge to the understanding of
charmonium production through NRQCD.

Many studies have been done to resolve this problem. For instance Ref.[4]
showed that to obtain cross-sections from NRQCD which are consistent with
data, one has to incorporate not only the NLO QCD corrections, (where the
value of total NLO contribution to cross section is nearly twice the LO
contribution), but also incorporate relativistic corrections. Incorporation of
NLO QCD corrections and relativistic corrections brings their results close to
the lower bound set by the Babar and Belle collaborations for the double-
charmonia production.

Further, it was recently pointed out that for both the processes,
$e^{-}e^{+}\rightarrow J/\Psi\eta_{c}$ and $e^{-}e^{+}\rightarrow
J/\Psi\chi_{cJ}$ in addition to the NLO QCD corrections, the interference
between the QCD and QED tree-level diagrams[6, 5], namely the
$O(\alpha_{em}^{2}\alpha_{s})$-order terms, can also provide signiﬁcant
contributions. Their studies have further shown that the colour octet
contributions to $e^{-}e^{+}\rightarrow J/\Psi\eta_{c}$, and
$e^{-}e^{+}\rightarrow J/\Psi\chi_{cJ}$ are negligibly small, with the
dominant contributions [5, 6] coming from the colour singlet channels.

Since it is seen that the perturbative QCD calculations fail at low energy
range typically below the $J/\Psi$ mass scale, thus, at this scale we can
employ the Bethe-Salpeter (BS) framework [7, 8, 9], which deals with the
relativistic motion of quarks inside the hadrons. Thus, in this work, we
calculate the cross section for the process taking the leading order (tree-
level) diagrams $\sim O(\alpha_{em}\alpha_{s})$ in QED and QCD for $J/\Psi$
and $\chi_{cJ}$, $(J=0,1)$ production in $e^{-}e^{+}$ annihilation through a
virtual photon, and a gluon being exchanged between two quark lines in the
quark-triangle diagram shown in Fig.1. In this work we focus only on the
colour singlet contributions, which is the dominant channel [6], and derive
analytical expressions for their cross sections in the framework of Bethe-
Salpeter equation, which is an equation that is firmly rooted in field theory
and is supposed to incorporate all the relativistic effects in the bound
states of $c\bar{c}$ systems.

Thus, in this work we study the production of ground and excited charmonium
meson pairs in $e^{-}e^{+}$ annihilation, such as:
$J/\Psi(1S)+~{}\chi_{c0}(1P);~{}\Psi(2S)~{}+\chi_{c0}(1P);~{}J/\Psi(1S)+~{}\chi_{c0}(2P);~{}\Psi(2S)~{}+\chi_{c0}(2P)$
along with $e^{-}+e^{+}\rightarrow
J/\Psi(1S)~{}+\chi_{c1}(1P);~{}\Psi(2S)~{}+\chi_{c1}(1P);~{}J/\Psi(1S)+~{}\chi_{c1}(2P);~{}\Psi(2S)~{}+\chi_{c1}(2P)$
through leading order (LO) tree-level diagrams $\sim
O(\alpha_{em}\alpha_{s})$, in the framework of the Bethe-Salpeter equation at
center of mass energy, $\sqrt{s}=10.6GeV.$

We wish to mention that Bethe-Salpeter Equation (BSE) [8, 10, 9, 7] is a
conventional approach in dealing with relativistic bound state problems. Due
to its firm base in quantum field theory and being a dynamical equation based
approach, it provides a realistic description for analyzing hadrons as
composite objects and can be applied to study not only the low energy hadronic
processes but also the high energy production processes involving quarkonia as
well. To simplify the calculations, we make use of the heavy quark limit on
the propagators of quark as shown in the next section.

The paper is organised as follows: Section 2 deals with the BS equation for
$Q\bar{Q}$ system, Section 3 deals with calculation of cross section for the
process, $e^{-}+e^{+}\rightarrow J/\Psi+\chi_{c0}$, while section 4 deals with
the corresponding calculation for the process, $e^{-}+e^{+}\rightarrow
J/\Psi+\chi_{c1}$. Section 5 deals with the Discussions.

##  2\. Bethe-Salpeter equation for $Q\bar{Q}$ bound state

A quark-anti quark bound state system can be described by a Bethe-Salpeter
equation (BSE),

$S_{F}^{-1}(p_{1})\Psi(P,q)S_{F}^{-1}(-p_{2})=i\int\frac{d^{4}q^{\prime}}{(2\pi)^{4}}K(q,q^{\prime})\Psi(P,q^{\prime}),$
(1)

where $p_{1}$ and $p_{2}$ are the momenta of the two particles, with the
internal momentum of the hadron being, $q$, and external hadron momentum $P$,
and mass, $M$. In Eq.(1), $K(q,q^{\prime})$ is the interaction kernel, and
$S_{F}^{-1}(\pm p_{1,2})=\pm i{\not}p_{1,2}+m_{1,2}$ are the inverse
propagators for the quark and antiquark.

To reduce the above equation to 3D form, we make use of the Covariant
Instantaneous Ansatz on the BS kernel, $K(q,q^{\prime})$, where we can write
$K(q,q^{\prime})=K(\widehat{q},\widehat{q}^{\prime})$, where the BS kernel
depends entirely on the component of internal momentum of the hadron,
$\widehat{q}_{\mu}=q_{\mu}-\frac{q.P}{P^{2}}P_{\mu}$, which is a 3D variable,
and is orthogonal to the total hadron momentum, i.e. $\widehat{q}.P=0$, while
$\sigma P_{\mu}=\frac{q.P}{P^{2}}P_{\mu}$ is the component of $q$ that is
longitudinal to $P$. And the 4-dimensional volume element is,
$d^{4}q=d^{3}\widehat{q}Md\sigma$.

It is to be observed that the longitudinal component, $M\sigma^{\prime}$ of
$q^{\prime}$ does not appear in $K(\hat{q},\hat{q}^{\prime})$. We thus carry
out integration over the longitudinal component, $Md\sigma^{\prime}$ of the
four dimensional volume element, $d^{4}q^{\prime}$ on the right side of
Eq.(1). Here, we use of the fact that

$\psi(\hat{q}^{\prime})=\frac{i}{2\pi}\int
Md\sigma^{\prime}\Psi(P,q^{\prime}),$ (2)

leading to equation,

$S_{F}^{-1}(p_{1})\Psi(P,q)S_{F}^{-1}(-p_{2})=\int\frac{d^{3}\hat{q}^{\prime}}{(2\pi)^{3}}K(\hat{q},\hat{q}^{\prime})\psi(\hat{q}^{\prime})=\Gamma(\hat{q}),$
(3)

where, $\Gamma(\hat{q})$ is the hadron-quark vertex function, and is directly
related to the 4D wave function, $\Psi(P,q)$. Multiplying the previous
equation from the left by $S_{F}(p_{1})$, and from the right by
$S_{F}(-p_{2})$, we one can express the 4D BS wave function $\Psi(P,q)$ in
terms of $\Gamma(\hat{q})$ as,

$\Psi(P,q)=S_{1}(p_{1})\Gamma(\hat{q})S_{2}(-p_{2}),$ (4)

where, the 4D hadron-quark vertex, that enters into the definition of the 4D
BS wave function in the previous equation is,

$\Gamma(\hat{q})=\int\frac{d^{3}\hat{q}^{\prime}}{(2\pi)^{3}}K(\hat{q},\hat{q}^{\prime})\psi(\hat{q}^{\prime}).$
(5)

Following a sequence of steps we get four Salpeter equations which are the
effective 3D forms of BSE (Salpeter equations)[11]:

$\displaystyle(M-\omega_{1}-\omega_{2})\psi^{++}(\hat{q})=\Lambda_{1}^{+}(\hat{q})\Gamma(\hat{q})\Lambda_{2}^{+}(\hat{q})$
$\displaystyle(M+\omega_{1}+\omega_{2})\psi^{--}(\hat{q})=-\Lambda_{1}^{-}(\hat{q})\Gamma(\hat{q})\Lambda_{2}^{-}(\hat{q})$
$\displaystyle\psi^{+-}(\hat{q})=0.$ $\displaystyle\psi^{-+}(\hat{q})=0$ (6)

These four Salpeter equations were used to derive the mass spectral equations
of $0^{++},0^{-+},1^{--},1^{++}$ and $1^{+-}$ mesons[11, 12]. The solutions of
these mass spectral equations not only lead to their mass spectrum but also
the analytical structures of their radial wave functions, that were used for
calculations of various transitions such as decay constants of these
quarkonia, as well as $M1$ and $E1$ transitions such as, $V\rightarrow
P\gamma,P\rightarrow V\gamma$, and $V\rightarrow S\gamma,S\rightarrow
V\gamma,A^{-}\rightarrow P\gamma,P\rightarrow A^{-}\gamma$ respectively [13,
14]. The analytic structures of these radial wave functions of $S,V$ and
$A^{+}$ mesons that are used for study of charmonium production processes:
$e^{-}+e^{+}\rightarrow J/\Psi+\chi_{cJ}(J=0,1)$ in the present work, are
given in Eqs.(12), and (23).

Thus, it is to be mentioned that in our framework, the component,
$\hat{q}_{\mu}$, is always orthogonal to $P_{\mu}$ and satisfies the relation,
$\hat{q}.P=0$, irrespective of whether $q.P=0$ (i.e. $\sigma=0$), or $q.P\neq
0$ (i.e. $\sigma\neq 0$). Further, due to the Lorentz-invariant nature of
$\hat{q}^{2}$, the applicability of this framework of Covariant Instantaneous
Ansatz is valid all the way from low energy spectra to high energy transition
amplitudes[11, 14].

The interaction kernel in BSE is taken as a one-gluon-exchange type as regards
the colour and spin dependence, and thus has a general structure [11],
$K(\hat{q},\hat{q}^{\prime})=(\frac{1}{2}\lambda_{1}.\frac{1}{2}\lambda_{2})\gamma_{\mu}\times\gamma_{\mu}V(\hat{q},\hat{q}^{\prime})$,
where the scalar part of the kernel is written as a sum of the one-gluon-
exchange part and the confinement part as,
$V(\hat{q},\hat{q}^{\prime})=\frac{4\pi\alpha_{s}}{(\hat{q}-\hat{q}^{\prime})^{2}}+\frac{3}{4}\omega^{2}_{q\bar{q}}\int
d^{3}r\bigg{(}\kappa
r^{2}-\frac{C_{0}}{\omega_{0}^{2}}\bigg{)}e^{i(\hat{q}-\hat{q}^{\prime}).\vec{r}}$,
and $\kappa=(1+4\hat{m}_{1}\hat{m}_{2}A_{0}M^{2}r^{2})^{-\frac{1}{2}}$
(see[11] for details).

Regarding the parameters of the model, $\omega_{0}=0.22$ GeV. is the spring
constant, $C_{0}=0.69$ is a dimensionless constant, while
$\frac{C_{0}}{\omega_{0}^{2}}$ plays the role of ground state energy,
$\Lambda=0.0250$ GeV. is the QCD length scale, with input quark masses,
$m_{u}=0.300GeV,m_{s}=0.430GeV,m_{c}=1.490GeV$, and $m_{b}=4.690$ GeV. Our
previous studies on mass spectral calculations of heavy-light quarkonia[11]
were used to fit the input parameters of our model. The cross sections are
calculated with the inverse range parameter $\beta=0.0900$GeV.

##  3\. Cross section for $e^{-}+e^{+}\rightarrow J/\Psi+\chi_{c0}$

The four colour singlet leading-order (LO) Feynman diagrams
$O(\alpha_{em}\alpha_{s})$ for the process, $e^{-}+e^{+}\rightarrow
J/\Psi+\chi_{c0}$, two of which are given in Figure 1. These diagrams along
with two other exchange diagrams form a gauge invariant subset. In all these
diagrams, the coupling of virtual photon to the $c\bar{c}$ mesons is described
by the quark loop diagrams, where the photon-quark-anti-quark vertex is given
as $ie_{Q}\gamma_{\mu}$, where, $e_{Q}=\frac{2}{3}e$ is the charge of the $c$
quark, and there is a single gluon exchange between two quark lines. Since the
amplitude for the process from all diagrams will be equal, we can express the
total amplitude for the process as four times the amplitude from the diagram
in Fig.2.

Figure 1: Two of the lowest order Feynman diagrams for the production of a
pair of charmonia in electron-positron annihilation. Other two diagrams can be
obtained by permutations Figure 2: Momentum labeling of the ﬁrst of the two
Feynman diagrams shown in Fig. 1

The invariant amplitude $M^{1}_{fi}$ for double charmonium production,
corresponding to Fig.2, is given by the one-loop momentum integral as:

$M^{I}_{fi}=ee_{Q}g^{2}F_{12}[\bar{v}^{(s2)}({\bar{p}}_{2})\gamma_{\mu}u^{(s1)}(\bar{p}_{1})]\frac{-1}{(\bar{p}_{1}+\bar{p}_{2})^{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\int\frac{d^{4}q^{\prime}}{(2\pi)^{4}}Tr[\bar{\Psi}_{V}(P,q)\gamma_{\beta}S_{F}(p_{3})\gamma_{\mu}\bar{\Psi}_{S}(P^{\prime},q^{\prime})\gamma_{\alpha}]\frac{\delta_{\alpha\beta}}{k^{2}},$
(7)

where, $P,q$ and $\epsilon^{\lambda}$ are the external momentum, internal
momentum and polarization vector of $J/\Psi$, while, $P^{\prime},q^{\prime}$
are corresponding variables of $\chi_{c0}$.
$F_{12}=\frac{1}{2}\vec{\lambda_{1}}.\frac{1}{2}\vec{\lambda_{1}}=\frac{4}{3}$
is the colour factor, and the wave functions of both incoming and outgoing
particles are normalized to one particle per unit volume. The above expression
can thus be factorized as the scalar product of the leptonic current and the
hadronic current as, $M^{I}_{fi}=J_{\mu}^{em}J_{\mu}^{hadronic}$. The total
amplitude for the process, $M_{fi}$ will be four times the amplitude,
$M^{I}_{fi}$ for diagram 1. The 4D adjoint BS wave functions of V and S mesons
can be written as,

$\displaystyle\bar{\Psi}_{V}(P,q)=S_{F}(-p_{2})\Gamma_{V}(\hat{q})S_{F}(p_{1})$
$\displaystyle\bar{\Psi}_{S}(P^{\prime},q^{\prime})=S_{F}(-p^{\prime}_{1})\Gamma_{S}(\hat{q}^{\prime})S_{F}(p^{\prime}_{2})$
(8)

where $-p_{2}=-p^{\prime}_{2}$. Here, we consider $p_{1,2}$ as the momenta of
quark/anti-quark of the vector meson, $J/\Psi$ of total momentum, $P$, while
$p^{\prime}_{1,2}$ as the momenta of the quark/anti-quark of the scalar
($\chi_{c0}$)/ axial vector ($\chi_{c1}$) meson of momentum $P^{\prime}$.
These quark momenta can be expressed in terms of the total and the internal
momenta of respective mesons as,

$p_{1}=\frac{1}{2}P+q,~{}~{}~{}~{}p_{2}=\frac{1}{2}P-q,~{}~{}~{}~{}p^{\prime}_{1}=\frac{1}{2}P^{\prime}-q^{\prime},~{}~{}~{}~{}~{}~{}~{}p^{\prime}_{2}=\frac{1}{2}P^{\prime}+q^{\prime}.$
(9)

To simplify the calculation, we assume that the propagators of the quark and
the gluon are independent of relative momenta $q$, and $q^{\prime}$ (see [8,
9, 10]). This simplification is appropriate since the masses of heavy quarks
are large compared with the relative momenta, which are of order
$\alpha_{s}m_{c}$. Thus the momenta $p_{3}$ and $k$ of the propagators are
large in comparison with the relative momenta, $q$ and $q^{\prime}$ of the
mesons. One may expect that, in the heavy quark limit, the calculation without
taking into account the relative momenta should be a good approximation.

Thus under heavy quark approximation, we can take,
$p_{1}=p_{2}\approx\frac{1}{2}P$, and
$p_{1}^{\prime}=p_{2}^{\prime}\approx\frac{1}{2}P^{\prime}$. Thus we can write
the momentum of the internal gluon propagator as,
$k=\bar{p_{1}}+\bar{p_{2}}\approx\frac{1}{2}(P+P^{\prime})$. Thus
$k^{2}\approx-\frac{s}{4}$. Similarly, we can express the momentum $p_{3}$ of
one of the internal quark lines as,
$p_{3}=p_{1}+k=\frac{1}{2}P^{\prime}+P+q\approx\frac{1}{2}P^{\prime}+P$.
Similarly, $p_{3}^{2}$ that enters into the denominator of this propagator can
be expressed as, $p_{3}^{2}\approx-\frac{s}{2}-M^{2}-\frac{M^{\prime 2}}{4}$.

Now,we make use of the fact that we decompose $q=(\hat{q},M\sigma)$, and
$q^{\prime}=(\hat{q}^{\prime},M^{\prime}\sigma^{\prime})$, while
$e_{Q}=\frac{2}{3}e$ is the charge of the charmed quark. Further we can
express the four-dimensional volume elements, $d^{4}q=d^{3}\hat{q}Md\sigma$,
and $d^{4}q^{\prime}=d^{3}\hat{q}^{\prime}M^{\prime}d\sigma^{\prime}$, and
making use of the fact that,

$\displaystyle\int\frac{Md\sigma}{2\pi}\Psi_{V}(P,q)=\Psi_{V}(\hat{q})$
$\displaystyle\int\frac{M^{\prime}d\sigma^{\prime}}{2\pi}\Psi_{S}(P^{\prime},q^{\prime})=\Psi_{S}(\hat{q}^{\prime}).$
(10)

Now, to calculate $M_{fi}$, we use the 3D decomposition of the general form of
4D BS wave functions [15] for vector and scalar mesons. Further, making use of
the most dominant Dirac structures for both vector and scalar mesons, that
contribute maximum to calculation of any meson observable[16, 17, 18], we take
the structures of 3D wave functions as [13]:

$\begin{split}\psi_{S}(\hat{q}^{\prime})=N_{S}(M^{\prime})\phi_{S}(\hat{q}^{\prime})\\\
\psi_{V}(\hat{q})=N_{V}(iM{\not}\epsilon)\phi_{V}(\hat{q}),\end{split}$ (11)

where the scalar functions, $\phi_{V}(\hat{q})$, and
$\phi_{S}(\hat{q}^{\prime})$ are obtained as analytic solutions of mass
spectral equations for vector and scalar mesons, that are in turn obtained
through solutions of their corresponding mass spectral equations[11]. These
radial wave functions for the ground and excited states of vector and scalar
mesons are[11]:

$\displaystyle\phi_{S}(1P,\hat{q})=\sqrt{\frac{2}{3}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{S}^{5/2}}\hat{q}e^{-\frac{\hat{q}^{2}}{2\beta_{S}^{2}}},$
$\displaystyle\phi_{S}(2P,\hat{q})=\sqrt{\frac{5}{3}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{S}^{5/2}}\hat{q}(1-\frac{2\hat{q}^{2}}{5\beta_{S}^{2}})e^{-\frac{\hat{q}^{2}}{2\beta_{S}^{2}}},$
$\displaystyle\phi_{V}(1S,\hat{q})=\frac{1}{\pi^{3/4}\beta_{V}^{3/2}}e^{-\frac{\hat{q}^{2}}{2\beta_{V}^{2}}},$
$\displaystyle\phi_{V}(2S,\hat{q})=\sqrt{\frac{3}{2}}\frac{1}{\pi^{3/4}\beta_{V}^{3/2}}(1-\frac{2\hat{q}^{2}}{3\beta_{V}^{2}})e^{-\frac{\hat{q}^{2}}{2\beta_{V}^{2}}},$
(12)

where $\beta_{v}$ and $\beta_{s}$ are their inverse range parameters. To
evaluate the trace over the gamma matrices, we need the expressions for the
adjoint 3D wave functions, $\bar{\Psi}_{V}(\hat{q})$ and
$\bar{\Psi}_{V}(\hat{q})$, where we express
$\bar{\Psi}=\gamma_{4}\Psi^{\dagger}\gamma_{4}$. We then write the invariant
amplitude $M_{fi}$ in the center of mass frame (where,
$\vec{\bar{p_{1}}}=-\vec{\bar{p_{2}}}$, and $E_{1}=E_{2}(=E))$, and with
$\sqrt{s}=2E$, after evaluation of trace as,

$M_{fi}=i\frac{128}{3}\frac{ee_{Q}g^{2}}{s}[\bar{v}^{(s2)}({\bar{p}}_{2}){\not}\epsilon
u^{(s1)}(\bar{p}_{1})]\frac{m_{c}N_{V}N_{S}MM^{\prime}}{k^{2}(-\frac{s}{2}-M^{2}-\frac{M^{\prime
2}}{4}+m_{c}^{2})}\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\phi_{V}(\hat{q})\int\frac{d^{3}\hat{q}^{\prime}}{(2\pi)^{3}}\phi_{S}(\hat{q}^{\prime})$
(13)

Now, to evaluate spin averaged invariant amplitude modulus squared,
$|\bar{M}_{fi}|^{2}$, we need to average over the initial spin states of
electron and positron, and sum over the polarization states of vector meson,
that is,
$|\bar{M}_{fi}|^{2}=\frac{1}{4}\sum_{s_{1},s_{2},\lambda}M^{\dagger}_{fi}M_{fi}$,
which after evaluation of trace over gamma matrices and making use of the
normalizations,
$\Sigma_{\lambda}\epsilon_{\mu}^{\lambda}\epsilon_{\nu}^{\lambda}=\frac{1}{3}(\delta_{\mu\nu}+\frac{P_{\mu}P_{\nu}}{M^{2}})$
for polarization vectors of the vector meson can be expressed as,

$\displaystyle|\bar{M}_{fi}|^{2}=\frac{2^{22}}{3^{4}}\frac{\pi^{4}\alpha_{em}^{2}\alpha_{s}^{2}}{s^{4}}\bigg{[}-\frac{1}{3}\bar{p}_{1}.\bar{p}_{2}+\frac{2}{3}\frac{(\bar{p}_{1}.P)(\bar{p}_{2}.P)}{M^{2}}+m_{e}^{2}\bigg{]}\frac{m_{c}^{2}N_{V}^{2}N_{S}^{2}M^{2}M^{\prime
2}}{(-\frac{s}{2}-M^{2}-\frac{M^{\prime
2}}{4}+m_{c}^{2})^{2}}\xi_{V}^{2}\xi_{S}^{2}$
$\displaystyle\xi_{V}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\phi_{V}(\hat{q})$
$\displaystyle\xi_{S}=\int\frac{d^{3}\hat{q^{\prime}}}{(2\pi)^{3}}\phi_{S}(\hat{q}^{\prime}),$
(14)

which has the general form,

$|\bar{M}_{fi}|^{2}=\frac{2^{22}}{3^{4}}\frac{\pi^{4}\alpha_{em}^{2}\alpha_{s}^{2}}{s^{4}}\frac{m_{c}^{2}N_{V}^{2}N_{S}^{2}M^{2}M^{\prime
2}}{(-\frac{s}{2}-M^{2}-\frac{M^{\prime
2}}{4}+m_{c}^{2})^{2}}\xi_{V}^{2}\xi_{S}^{2}L_{\mu\mu^{\prime}}H_{mu\mu^{\prime}},$
(15)

where $L_{\mu\mu^{\prime}}$ is the leptonic tensor, while
$H_{\mu\mu^{\prime}}$ is the hadronic tensor. Here electromagnetic coupling
constant, $\alpha_{em}=\frac{e^{2}}{4\pi}$, and the strong coupling constant,
$\alpha_{s}=\frac{g^{2}}{4\pi}=\frac{1}{log(m_{c}/\Lambda_{QCD})}$, $m_{e}$ is
electron mass, while $m_{c}$ is the c-quark mass, and we have made use of
$k^{2}=-\frac{s}{4}$. Now, $M_{fi}$ in previous equation involves the dot
products of momenta, which can be expressed in the center of mass frame as,
$\bar{p}_{1}.\bar{p}_{2}=\frac{1}{2}(2m_{e}^{2}-s),~{}~{}~{}~{}\bar{p}_{1}.P\approx-E^{2}(1-cos\theta),~{}~{}~{}~{}\bar{p}_{2}.P\approx-E^{2}(1+cos\theta)$,
with $\sqrt{s}=2E=E_{cm}$, being the center of mass energy. Here,
$(\bar{p}_{1}+\bar{p}_{2})=\sqrt{s}$. It can be checked that the asymptotic
behaviour of amplitude, $M_{fi}\sim 1/s^{3}$.

Further, $N_{S}$ and $N_{V}$ are the 4D BS normalizers of vector and scalar
mesons, that are obtained through the current conservation condition,

$2iP_{\mu}=\int\frac{d^{4}q}{(2\pi)^{4}}\mbox{Tr}\left\\{\overline{\Psi}(P,q)\left[\frac{\partial}{\partial
P_{\mu}}S_{F}^{-1}(p_{1})\right]\Psi(P,q)S_{F}^{-1}(-p_{2})\right\\}+(1\rightleftharpoons
2),$ (16)

and are expressed as,

$\displaystyle
N_{V}^{-2}=4\hat{m}_{1}\widehat{m}_{2}M^{2}\frac{1}{m_{c}}\int\frac{d^{3}\widehat{q}}{(2\pi)^{3}}\phi_{V}^{2}(\widehat{q}),$
$\displaystyle N_{S}^{-2}=4\hat{m}_{1}\widehat{m}_{2}M^{\prime
2}\frac{1}{m_{1}}\int\frac{d^{3}\widehat{q}^{\prime}}{(2\pi)^{3}}\phi_{S}^{2}(\widehat{q}^{\prime}),$
(17)

respectively.

We can express $d\sigma$ in terms of the transition rate per unit volume,
$W_{fi}$ as,

$d\sigma=\frac{W_{fi}(2\pi)^{4}\delta^{(4)}(p_{1}+p_{2}-P^{\prime}-P^{\prime\prime})}{F}[\frac{d^{3}\vec{P}^{\prime}}{2E^{\prime}(2\pi)^{3}}\frac{d^{3}\vec{P}^{\prime\prime}}{2E^{\prime\prime}(2\pi)^{3}}],$
(18)

where the factor in the square brackets on the extreme right is the total
phase space volume for the two outgoing particles in the center of mass frame
considered. And $F$ is the incident flux, which is expressed as,
$F=4\sqrt{(\bar{p_{1}}.\bar{p_{2}})^{2}-m_{1}^{2}m_{2}^{2}}$ (which includes
the kinematical factors of initial particles), and can in turn be expressed in
the Center of mass frame as, $F=4|\vec{p}|\sqrt{s}$, where $|\vec{p}|$ is the
magnitude of momentum of either of the incident particles in the center of
mass frame, and $\sqrt{s}=2E$ is the center of mass energy of incident
particles. Following a series of steps, we can express the total cross section
as,

$\sigma=\frac{1}{32\pi^{2}s^{3/2}}|\vec{P^{\prime}}|\int
d\Omega^{\prime}|\bar{M}_{fi}|^{2},$ (19)

with
$|\vec{P^{\prime}}|=\sqrt{\frac{1}{s}[s-(M+M^{\prime})^{2}][s-(M-M^{\prime})^{2}]}$
being the momenta of either of the outgoing particles in the center of mass
frame, and $|\bar{M}_{fi}|^{2}$, being expressed as in Eq.(7). As regards the
behaviour of various quantities, it can be checked that
$\beta_{S}\sim\beta_{V}\sim M$, The BS normalizers behave as, $N_{S}\sim
N_{V}\sim M^{-1/2}$, while the 3D integrals, $\xi_{S}\sim\xi_{V}\sim M^{3/2}$.
Thus the cross sectional formula behaves as, $\sigma\sim M^{-2}$.

Thus we have expressed the amplitude and the cross section for the process,
$e^{+}e^{-}\rightarrow J/\Psi\chi_{c0}$ in terms of the radial wave functions
of both the produced vector and scalar quarkonia, which were analytically
obtained from solutions[11] of mass spectral equations of vector and scalar
mesons. The cross sections for
$e^{-}~{}~{}e^{+}\rightarrow\Psi(nS)~{}~{}\chi_{c0}(nP)$ calculated in this
work are given in Table 1 along with experimental data and results of other
models. It is seen that the contribution of wave functions at the origin
$\xi_{V}$ and $\xi_{S}$ have a major role to play in the calculation of the
cross section. Their numerical values are: 0.0308 GeV. and 0.0336 GeV. for
$J/Psi$ and $\chi_{c0}(1P)$ respectively.

Process |  BSE-CIA |  Expt. [3] |  [6] |  [4] |  [19]
---|---|---|---|---|---
$e^{-}e^{+}\rightarrow J/\Psi(1S)\chi_{c0}(1P)$ |  6.771 |  17$\pm$ 8$\pm$7 |  6.62 |  2.40$\pm$1.02 |  5
$e^{-}e^{+}\rightarrow\Psi(2S)\chi_{c0}(1P)$ |  3.639 |  |  |  1.00$\pm$0.42 |
$e^{-}e^{+}\rightarrow\Psi(1S)\chi_{c0}(2P)$ |  2.950 |  |  |  |
$e^{-}e^{+}\rightarrow\Psi(2S)\chi_{c0}(2P)$ |  3.176 |  |  |  |
Table 1: Cross sections for process,
$e^{-}~{}~{}e^{+}\rightarrow\Psi(nS)~{}~{}\chi_{c0}(nP)$ calculated in BSE-CIA
along with experimental data and results of other models in fb

##  4\. Cross section for $e^{-}+e^{+}\rightarrow J/\Psi+\chi_{c1}$

In lowest order the process proceeds through the same Feynman diagrams given
in Fig.1, except that now we have $J/\Psi$ and $\chi_{c1}$ in the final state,
with external momentum, $P^{\prime}$, internal momentum, $q^{\prime}$, and
polarization vector, $\epsilon^{\lambda^{\prime}}$, while the corresponding
variables for $J/\Psi$ remain the same as, $P,q$ and $\epsilon^{\lambda}$ as
in Diagram 1.

Amplitude for this process can thus be written as,

$M^{I}_{fi}=ee_{Q}g^{2}F_{12}[\bar{v}^{(s2)}({\bar{p}}_{2})\gamma_{\mu}u^{(s1)}(\bar{p}_{1})]\frac{-1}{(\bar{p}_{1}+\bar{p}_{2})^{2}}\int\frac{d^{4}q}{(2\pi)^{4}}\int\frac{d^{4}q^{\prime}}{(2\pi)^{4}}Tr[\bar{\Psi}_{V}(P,q)\gamma_{\beta}S_{F}(p_{3})\gamma_{\mu}\bar{\Psi}_{A}(P^{\prime},q^{\prime})\gamma_{\alpha}]\frac{\delta_{\alpha\beta}}{k^{2}},$
(20)

where, the 4D BS adjoint BS wave functions are: The 4D adjoint BS wave
functions of V and S mesons can be written as,

$\displaystyle\bar{\Psi}_{V}(P,q)=S_{F}(-p_{2})\Gamma_{V}(\hat{q})S_{F}(p_{1})$
$\displaystyle\bar{\Psi}_{A}(P^{\prime},q^{\prime})=S_{F}(p^{\prime}_{1})\Gamma_{A}(\hat{q}^{\prime})S_{F}(-p^{\prime}_{2})$
(21)

We make use of the definitions of quark momenta in Eq.(2), and make use of the
definitions of the quark momenta mentioned in paragraph after Eq.(9). Further,
the 4D volume elements, $d^{4}q=d^{3}\hat{q}Md\sigma$, and
$d^{4}q^{\prime}=d^{3}\hat{q}^{\prime}M^{\prime}d\sigma^{\prime}$, and we make
use of the fact that,

$\displaystyle\int\frac{Md\sigma}{2\pi}\Psi_{V}(P,q)=\Psi_{V}(\hat{q})$
$\displaystyle\int\frac{M^{\prime}d\sigma^{\prime}}{2\pi}\Psi_{A}(P^{\prime},q^{\prime})=\Psi_{A}(\hat{q}^{\prime}).$
(22)

Now, to calculate $M_{fi}$, we use the 3D decomposition of the general form of
4D BS wave functions [15] for vector and scalar mesons, making use of the most
dominant Dirac structures for both vector and scalar mesons, that contribute
maximum to calculation of any meson observable, as in [13, 14], and take,

$\displaystyle\psi_{V}(\hat{q})=N_{V}iM{\not}\epsilon^{\lambda}\phi_{V}(\hat{q})$
$\displaystyle\psi_{A}(\hat{q}^{\prime})=N_{A}iM^{\prime}\gamma_{5}{\not}\epsilon^{\lambda^{\prime}}\phi_{A}(\hat{q}^{\prime}).$
(23)

where the ground and excited radial wave functions of axial mesons that are
obtained through the solutions of their mass spectral equations are:

$\displaystyle\phi_{A}(1P,\hat{q})=\sqrt{\frac{2}{3}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{A}^{5/2}}\hat{q}e^{-\frac{\hat{q}^{2}}{2\beta_{A}^{2}}},$
$\displaystyle\phi_{A}(2P,\hat{q})=\sqrt{\frac{5}{3}}\frac{1}{\pi^{3/4}}\frac{1}{\beta_{A}^{5/2}}\hat{q}(1-\frac{2\hat{q}^{2}}{5\beta_{A}^{2}})e^{-\frac{\hat{q}^{2}}{2\beta_{A}^{2}}}$
(24)

Performing integrals over the longitudinal components, $Md\sigma$ and
$M^{\prime}d\sigma^{\prime}$ in Eq.(20) through use of Eq.(22), putting the
expressions for 3D BS wave functions $\Psi_{V}(\hat{q})$ and
$\Psi_{A}(\hat{q}^{\prime})$, evaluating trace over the gamma matrices, we can
express the amplitude, $M_{fi}$ as,

$\displaystyle
M_{fi}=\frac{2^{9}ee_{Q}g^{2}}{3s^{2}(\frac{s}{2}+M^{2}+M^{\prime
2}/2-m_{c}^{2})}N_{A}N_{V}\xi_{V}\xi_{A}[\bar{v}^{(s_{2})}(p_{2})\gamma_{\mu}u^{(s_{1})}(p_{1})]\epsilon_{\mu\nu\alpha\beta}\epsilon^{\lambda}_{\nu}\epsilon^{\lambda^{\prime}}_{\alpha}(P+\frac{1}{2}P^{\prime})_{\beta};$
$\displaystyle\xi_{V}=\int\frac{d^{3}\hat{q}}{(2\pi)^{3}}\phi_{V}(\hat{q})$
$\displaystyle\xi_{A}=\int\frac{d^{3}\hat{q^{\prime}}}{(2\pi)^{3}}\phi_{A}(\hat{q}^{\prime}),$
(25)

Now, to evaluate spin averaged invariant amplitude modulus squared,
$|\bar{M}_{fi}|^{2}$, we need to average over the initial spin states of
electron and positron, and sum over the polarization states of vector meson
and axial meson, that is,
$|\bar{M}_{fi}|^{2}=\frac{1}{4}\sum_{s_{1},s_{2},\lambda,\lambda^{\prime}}M^{\dagger}_{fi}M_{fi}$,
which after evaluation of trace over gamma matrices and making use of the
normalizations,
$\Sigma_{\lambda}\epsilon_{\mu}^{\lambda}\epsilon_{\nu}^{\lambda}=\frac{1}{3}(\delta_{\mu\nu}+\frac{P_{\mu}P_{\nu}}{M^{2}})$
for polarization vectors of the vector meson, and
$\Sigma_{\lambda^{\prime}}\epsilon_{\mu}^{\lambda^{\prime}}\epsilon_{\nu}^{\lambda^{\prime}}=\frac{1}{3}(\delta_{\mu\nu}+\frac{P^{\prime}_{\mu}P^{\prime}_{\nu}}{M^{\prime
2}})$ for polarization vectors of the axial meson, can be expressed as,

$\displaystyle|\bar{M}_{fi}|^{2}=-\frac{2^{26}\alpha_{em}^{2}\alpha_{s}^{2}\pi^{4}}{81s^{4}(\frac{s}{2}+M^{2}+M^{\prime
2}/2-m_{c}^{2})^{2}}\epsilon_{\mu\nu\alpha\beta}~{}~{}\epsilon_{\mu^{\prime}\nu^{\prime}\alpha^{\prime}\beta^{\prime}}\bigg{(}p_{1\mu^{\prime}}p_{2\mu}+p_{1\mu}p_{2\mu^{\prime}}+(-p_{1}.p_{2}+m_{e}^{2})\delta_{\mu\mu^{\prime}}\bigg{)}$
$\displaystyle\bigg{(}\delta_{\nu\nu^{\prime}}\delta_{\alpha\alpha^{\prime}}+\frac{P_{\nu}P_{\nu^{\prime}}}{M^{2}}\delta_{\alpha\alpha^{\prime}}+\frac{P_{\alpha}P_{\alpha^{\prime}}}{M^{\prime
2}}\delta_{\nu\nu^{\prime}}+\frac{P_{\nu}P_{\nu^{\prime}}P^{\prime}_{\alpha}P_{\alpha^{\prime}}}{M^{2}M^{\prime
2}}\bigg{)}(P+\frac{1}{2}P^{\prime})_{\beta}(P+\frac{1}{2}P^{\prime})_{\beta}^{\prime}N_{A}^{2}N_{V}^{2}M^{2}M^{\prime
2}\xi_{v}^{2}\xi_{A}^{2}.$ (26)

Thus, we express $|\bar{M}_{fi}|^{2}$ in terms of dot products of momenta in
the center of mass frame as,
$\bar{p}_{1}.\bar{p}_{2}=\frac{1}{2}(2m_{e}^{2}-s),~{}~{}~{}~{}\bar{p}_{1}.P\approx-E^{2}(1-cos\theta),~{}~{}~{}~{}\bar{p}_{2}.P\approx-E^{2}(1+cos\theta),\bar{p}_{1}.P^{\prime}\approx-E^{2}(1+cos\theta)=\bar{p}_{2}.P$
as,

$\displaystyle|\bar{M}_{fi}|^{2}=-\frac{2^{28}\alpha_{em}^{2}\alpha_{s}^{2}\pi^{4}}{9^{3}s^{4}(\frac{s}{2}+M^{2}+M^{\prime
2}/2-m_{c}^{2})^{2}}N_{A}^{2}N_{V}^{2}M^{2}M^{\prime
2}\xi_{v}^{2}\xi_{A}^{2}|\overline{TR}|^{2};$
$\displaystyle|\overline{TR}|^{2}=\frac{1}{128M^{2}M^{\prime
2}}\bigg{[}4M^{4}\bigg{(}64m_{e}^{2}M^{\prime 2}+s(-96M^{\prime
2}+19s-7s\cos^{2}\theta)\bigg{)}+$ $\displaystyle s\bigg{(}-8M^{\prime
6}+3s^{3}(-5+\cos^{2}\theta)+2M^{\prime
2}s(8m_{e}^{2}-5s+7s\cos^{2}\theta)+M^{\prime
4}s(15+\cos\theta(3+4\cos\theta))\bigg{)}+$ $\displaystyle
2M^{2}\bigg{(}32m_{e}^{2}(M^{\prime 4}+4M^{\prime 2}s+s^{2})+s[56M^{\prime
4}+4s^{2}(3+\cos^{2}\theta)+M^{\prime
2}s(63+\cos\theta(-8+3\cos\theta))]\bigg{)}\bigg{]}$ (27)

Thus, the total cross section can then be expressed as,

$\sigma=\frac{1}{16\pi s^{3/2}}|\vec{P^{\prime}}|\int
d\cos\theta|\bar{M}_{fi}|^{2}.$ (28)

Our results of cross sections for
$e^{-}~{}~{}e^{+}\rightarrow\Psi(nS)~{}\chi_{c1}(nP)$ are listed in table 2
along with experimental data and results of other models.

Process |  BSE-CIA |  Expt. [3] |  [6] |  [4] |  [19]
---|---|---|---|---|---
$e^{-}e^{+}\rightarrow J/\Psi(1S)\chi_{c1}(1P)$ |  9.446 |  $<$ 18 |  0.53 |  0.38$\pm$0.12 |  18
$e^{-}e^{+}\rightarrow\Psi(2S)\chi_{c1}(1P)$ |  7.535 |  |  |  |
$e^{-}e^{+}\rightarrow\Psi(1S)\chi_{c1}(2P)$ |  6.577 |  |  |  |
$e^{-}e^{+}\rightarrow\Psi(2S)\chi_{c1}(2P)$ |  5.199 |  |  |  |
Table 2: Cross sections for process,
$e^{-}~{}~{}e^{+}\rightarrow\Psi(nS)~{}~{}\chi_{c1}(nP)$ calculated in BSE-CIA
along with experimental data and results of other models in fb

##  5\. Discussions

In this work we study the production of ground and excited charmonium pairs in
$e^{-}e^{+}\rightarrow\Psi(nS)+\chi_{cJ}(nP)$ for $J=0,1$ and $n=1,2)$,
through leading order (LO) tree-level diagrams $\sim
O(\alpha_{em}\alpha_{s})$, which proceed through exchange of a virtual photon
and an internal gluon line connecting two quark lines in the triangle quark
loop part of the diagram, in the framework of $4\times 4$ Bethe-Salpeter
equation under Covariant Instantaneous Ansatz, at center of mass energy,
$\sqrt{s}=10.6GeV.$ One of the biggest challenges in calculations based on
NRQCD was that the theoretical predictions of cross sections for the
processes, $e^{-}e^{+}\rightarrow J/\Psi+\chi_{cJ}(nP)$, and
$e^{-}e^{+}\rightarrow J/\Psi+\eta_{c}$ using leading order QCD diagrams is an
order of magnitude lower than the experimental measurements [1, 2, 3] of the
same. And it was only after use of NLO QCD diagrams and incorporation of
relativistic corrections that the results of cross sections came near the
lower bound set by BABAR and Belle[1, 2, 3].

In these calculations, we try to show that the cross sections for double
charmonium production calculated with use of the four leading order $\sim
O(\alpha_{em}\alpha_{s})$ diagrams using BSE approach come close to
experimental data, which might be mainly due to the consistent treatment of
motion of quarks inside the hadrons in the BSE framework, which is firmly
rooted in field theory, though to simplify evaluation of integrals in
amplitude, $M_{fi}$ we have made use of the heavy-quark approximation on the
quark propagator, $S_{F}(p_{3})$ and the gluon propagator, where the
propagators of the quark and the gluon are independent of internal momenta
$q$, and $q^{\prime}$ [8, 9, 10]. This simplification seems appropriate since
the masses of heavy quarks are large compared with the relative momenta, which
are of order $\alpha_{s}m_{c}$. Thus the momenta $p_{3}$ and $k$ of the
propagators can be considered to be large in comparison with the internal
momenta of the two mesons.

We do expect that the above heavy-quark approximation on the quark and gluon
propagators will not lead to a value of cross section that is drastically
different from the calculation where we employ the dependence of quark and
gluon momenta on the internal momenta of the two hadrons, which we intend to
do next on lines of our previous works on $M1$ and $E1$ decays [13, 14] of
quarkonia.

However what is interesting in the present calculation is the simplicity by
which one can deduce the main features of these cross section calculations for
these interesting processes using the leading order diagrams, $\sim
O(\alpha_{em}\alpha_{s})$ in QCD in the framework of BSE, namely the
consistent treatment of internal motion of quarks in hadrons in BSE framework,
which we feel is the single most important reason for obtaining theoretical
values of cross sections of the same order of magnitude as the lower bound set
by BABAR and Belle[1, 2, 3]. We wish to point out this is further validated by
a recent calculation [20] on $e^{-}e^{+}\rightarrow J/\Psi\eta_{c}$ using
Light cone, where the authors had shown that by taking intrinsic motion of
quarks inside the hadrons, one can significantly increase the value of cross
section.

This calculation involving production of double charmonia in electron-positron
annihilation can be easily extended to studies on other processes (involving
the exchange of a single virtual photon) observed at B-factories such
as,$e^{+}e^{-}\rightarrow\gamma\chi_{cJ}$, and also to processes proceeding
through exchange of two virtual photons such as $e^{+}e^{-}\rightarrow J/\Psi
J/\Psi$.

Among other similar works, in [7], the authors study in the framework of
Bethe-Salpeter equation, the production ratio of neutral to charged kaon pair
in $e^{+}e^{-}$ annihilation below the $J/\psi$ mass by employing the two
leading Dirac structures in hadronic wave functions as per the power counting
rule [16, 17] we had proposed some time ago. We wish to extend this study to
calculation of cross section for each of these processes with the involvement
of all the possible Dirac structures.

## References

  * [1] K.F.Chen et al.(Belle collaboration), Phys. Rev. D82, 091106(R) (2010).
  * [2] B.Auger et al.(BaBar collaboration), Phys. Rev. Lett. 103, 161801 (2009).
  * [3] T.V.Ulgov, Eur.Phys.J.C33, s01,s235 (2005).
  * [4] E.Braaten, J.Lee, arxiv: hep-ph/0211085.
  * [5] Z.Sun, X.G.Wu, Y. Ma, S.J.Brodsky, Phys. Rev. D98, 094001 (2018).
  * [6] Y.Z.Jiang, Z.Sun, Eur. Phys. J. C78, 892 (2018).
  * [7] J.Yi, S-Y.Li,Y-R.Liu, Z-X.Meng, Z-G.Si, T.Yao, Phys. Rev. C102, 015201 (2020).
  * [8] X-H. Guo, H-W. Ke, X-Q. Li, X-H. Wu arxiv:0804.0949[hep-ph], X. Guo, H. Jin and X. Li, Phys. Rev. D 53, 1153 (1996).
  * [9] E.Mengesha, S.Bhatnagar, Intl. J. Mod. Phys. E20, 2521 (2011).
  * [10] H.Negash, S.Bhatnagar, Adv. in HEP 2019, 4029356 (2019).
  * [11] E.Gebrehana, S.Bhatnagar, H.Negash, Phys. Rev. D100, 054034 (2019).
  * [12] V.Guleria, S.Bhatnagar, Intl. J. Theo. Phys. 60, 3143 (2021).
  * [13] S.Bhatnagar, E.Gebrehana, Phys. Rev. D102, 094024 (2020).
  * [14] V.Guleria, E.Gebrehana, S.Bhatnagar, Phys. Rev. D104, 094045 (2021).
  * [15] C. H. L. Smith, Ann. Phys. (N.Y.) 53, 521 (1969).
  * [16] S.Bhatnagar, S-Y.Li, J. Phys. G32, 949 (2006).
  * [17] S.Bhatnagar, S.Y.Li, J.Mahecha, Int. J. Mod. Phys. E 20, 1437 (2011).
  * [18] S.Bhatnagar, J.Mahecha, Y.Mengesha, Phys. Rev. D90,014034(2014).
  * [19] K.T.Chao, L.K.Hao, arxiv: hep-ph/0209189.
  * [20] V.V.Braguta, A.K.Likhoded, A.V.Luchinsky, Phys. Rev. D72, 074019 (2005).


# Laser-driven pointed acceleration of electrons with preformed plasma lens

K.A. Ivanov<EMAIL_ADDRESS>Faculty of Physics, Lomonosov Moscow State
University, 119991, Moscow, Russia  Lebedev Physical Institute of Russian
Academy of Sciences, 119991, Moscow, Russia  D.A. Gorlova  Faculty of Physics,
Lomonosov Moscow State University, 119991, Moscow, Russia  Institute for
Nuclear Research of Russian Academy of Sciences, 117312, Moscow, Russia  I.N.
Tsymbalov  Faculty of Physics, Lomonosov Moscow State University, 119991,
Moscow, Russia  Institute for Nuclear Research of Russian Academy of Sciences,
117312, Moscow, Russia  I.P. Tsygvintsev  Keldysh Institute of Applied
Mathematics of Russian Academy of Sciences, 125047 Moscow, Russia  S.A.
Shulyapov  Faculty of Physics, Lomonosov Moscow State University, 119991,
Moscow, Russia  R.V. Volkov  Faculty of Physics, Lomonosov Moscow State
University, 119991, Moscow, Russia  A. Savel’ev  Faculty of Physics, Lomonosov
Moscow State University, 119991, Moscow, Russia  Lebedev Physical Institute of
Russian Academy of Sciences, 119991, Moscow, Russia

(October 2023)

###### Abstract

The simultaneous laser-driven acceleration and angular manipulation of the
fast electron beam is experimentally demonstrated. The bunch of multi-MeV
energy charged particles is generated during the propagation of the
femtosecond laser pulse through the near-critical plasma slab accompanied by
plasma channeling. Plasma is formed by the controlled breakdown of a thin-tape
target by a powerful nanosecond prepulse. The electron beam pointing approach
is based on the refraction of a laser pulse in the presence of a strong radial
density gradient in the breakdown of the tape with a small displacement of the
femtosecond laser beam relative to the breakdown symmetry axis. A shift of
several micrometers makes it possible to achieve beam deflection by an angle
up to 10 degrees with acceptable beam charge and spectrum conservation. This
opens up opportunities for in-situ applications for scanning objects with an
electron beam and the multistage electron beam energy gain in consecutive
laser accelerators without bulk magnetic optics for particles. Experimental
findings are supported by numerical Particle-In-Cell calculations of laser-
plasma acceleration and hydrodynamic simulations.

##  I Introduction

The generation of energetic collimated electron beams with a high charge using
laser-plasma accelerators has become one of the main area of application for
terawatt (TW) and petawatt (PW) laser systems Malka (2012); Joshi (2019);
Ledingham and Galster (2010). Thanks to Laser Wakefield Acceleration (LWFA
Kurz  _et al._ (2021)) and Direct Laser Acceleration (DLA Hussein  _et al._
(2021); Arefiev  _et al._ (2016)) in relativistic plasma channel, together
with the rapid development of high repetition rate TW-class lasers the new
phenomena related to astrophysics Takabe and Kuramitsu (2021), nuclear
photonics Nedorezov  _et al._ (2021), biomedicine and diagnostics Nakajima
(2015); Labate  _et al._ (2020), etc. has become available to a wide community
of researchers literally within the laboratory. The latter circumstance
imposes fairly severe restrictions on the size of laboratory sources utilizing
compact laser systems with peak power of up to 100 TW (including new-
generation systems with peak power of a few TW with a repetition rate up to a
kHz level Antipenkov  _et al._ (2021)). Besides, the design of the electron
beam generator itself (target assembly and adjacent elements) must have
reasonable dimensions and provide simple and stable operation as well as
control over the beam. One should also consider the bunch acceleration
efficiency in terms of mean energy and charge. Instead of a single
accelerating stage by an extremely powerful laser pulse the use of a few laser
pulses can provide consecutive bunch energy gain in a number of acceleration
stages driven by less-demanding and modest laser systems Kim  _et al._ (2013);
Zhu  _et al._ (2022); Jin  _et al._ (2019).

For various application scenarios the angular targeting and pointing of the
electron beam are of key importance, whether it concerns the optimal injection
into the subsequent stage of acceleration Massimo  _et al._ (2019);
Starodubtseva  _et al._ (2023), injection schemes Khudiakov and Pukhov (2022)
relevant to the AWAKE experiment Moschuering  _et al._ (2019), scanning of a
sample, directing the bunch into different detectors, radio medicine
applications including FLASH therapy Favaudon  _et al._ (2014), etc. For this,
standard accelerator technology can be used: magnetic lenses and dipoles. But
their key drawback is high energy selectivity, whereas the laser acceleration
of particles beam with high charge leads to rather wide energy spectrum Gahn
_et al._ (1999). Also worth noting the bulkiness and slow adjustment.
Therefore, the issue of controlling the beam pointing is of great importance
not only from the point of view of compactness, but also the preservation of a
high charge. For these purposes, several purely plasma methods have been
proposed for the LWFA scheme. Among them, we can single out the waveguide
propagation of a pulse, when a laser beam is introduced onto the axis of an
electron bunch through a curved plasma channel Nakajima (2018); Luo  _et al._
(2018). It was also proposed to use a pulse with an inclined wave front, which
will turn in a low-density plasma and take the electron beam behind it in the
Wakefield acceleration mode Mittelberger  _et al._ (2019); qing Zhu  _et al._
(2020).

In this work, we experimentally demonstrate the possibility of controlling the
ejection angle of an electron bunch with particle energy exceeding a few MeV
accelerated by a 1 TW laser pulse in a plasma channel. The latter is formed in
a plasma sheet of near critical density when a thin (15 $\mu$m) mylar tape is
hole bored by an additional nanosecond prepulse (100 mJ, 8 ns, $10^{13}$ W
cm-2), coming ahead of the main accelerating femtosecond pulse (50 mJ, 50 fs,
$5\times 10^{18}$ W cm-2) by a few nanoseconds. The prepulse forms a breakdown
with strong radial density gradient. Precise adjustment of the breakdown axis
with respect to the femtosecond pulse axis makes it possible to control the
trajectory of the relativistic plasma channel due to laser beam refraction.
The accelerated electron beam preserves its characteristics (energy, spectrum,
divergence) at a high degree. The findings are confirmed by simulations and
analytical model.

##  II Main experimental results

Previously, we proposed an efficient laser-driven DLA electron beam source
using subcritical plasma and moderate peak power (1 TW) laser system Tsymbalov
_et al._ (2020) operating at 10 Hz. The electron bunch had a $\sim$50 pC
charge for the energy range above $\sim$3 MeV. To form a plasma with the
required electron density longitudinal profile, we used a nanosecond prepulse
from a Nd:YAG laser at a wavelength of 1064 nm, focused to a peak intensity
$\sim 10^{12}$ W cm-2 onto the surface of thin plastic tape target at an angle
close to the normal. After the peak of the prepulse with delay of several
nanoseconds the terawatt radiation from the TiSa laser system (805 nm, 50 fs,
$>10^{18}$ W cm-2) was focused into the preformed plasma along the same
optical axis. With such an interaction geometry, efficient generation of an
electron beam is possible under Direct Laser Acceleration in the self-focusing
induced plasma channel. The optimal delay between the peaks of the pulses
depends on the target thickness and prepulse intensity, since it determines
nanosecond plasma evolution dynamics of the tape target. In the current
experimental campaign we used a 15 $\mu$m thick mylar tape target. The
prepulse intensity was on the level of $10^{13}$ W cm-2 with a focal diameter
13 $\mu$m spot at full width at half maximum (FWHM), which resulted in the
optimal delay where the electron beam is generated equal to $\sim$2
nanoseconds. The femtosecond pulse was focused into a 4.0$\pm$0.2 $\mu$m FWHM
spot to a peak intensity $\approx 5\times 10^{18}$ W cm-2. The relative
amplitude of the amplified spontaneous emission (ASE) pedestal in the domain
10-100 picoseconds prior to main pulse was on the level of 10-8. Obviously,
the density distribution created by the nanosecond pulse is radially dependent
and the plasma acts as a lens: one can estimate the laser beam refraction
under its propagation with a small deviation D from the longitudinal symmetry
axis of the plasma hole in the tape target. For this we implemented the
updated experimental setup (Fig.1) which makes possible to shift the axis of
the heating nanosecond pulse with respect to the femtosecond pulse axis by
tilting the mirror (2) in the Nd:YAG laser path.

Figure 1:  Experimental scheme: 1 - Nd:YAG nanosecond prepulse, 2 - adjustable
mirror for Nd:YAG radiation, 3 - femtosecond pulse, 4 - tape target, 5 -
electron beam detector (Lanex scintillation screen and camera with lens), 6 -
camera with microobjective collecting the scattered light from plasma
(perpendicular to observer), D \- deviation of the Ti:Sa laser pulse axis
relative to the Nd:YAG axis (target breakdown symmetry axis). Figure 2:  Upper
row - the averaged image of the electron beam on the detector with D = 0 (a),
5 (b) and 10 (c) $\mu$m. Second row - the electron beam stamp in single shots
at D = 0 (a), 5 (b) and 10 (c) $\mu$m. Bottom row - the plasma scattered light
signal in the optical range under the same femtosecond pulse propagation
conditions: D = 0 (coaxial, g), 5 $\mu$m (h) and 10 $\mu$m (i). Red arrow
indicates the incident laser beam direction. The red lines mark the tape
initial position.

The alignment of the axes of nano- and femtosecond pulses was controlled by
transferring images of the waist of the off-axis parabolic mirror by a
microobjective to a CCD camera with high optical resolution. The discreteness
of the D was determined by a stepper actuator on the mount of the Nd:YAG laser
mirror and was 2.5$\pm$0.5 $\mu$m. The electron beam was detected on a Lanex
scintillating screen in the direction of laser radiation. A metal filter (300
$\mu$m of tungsten) directly in front of the scintillator cut off electrons
with energies below $\sim$3 MeV without strong additional beam divergence
after scattering in the filter. Scattered light from the plasma in the visible
range was also detected in the perpendicular to laser propagation plane as the
femtosecond pulse passed through the tape. For this, the region of interest
was transferred to a CCD camera using a microobjective.

Fig.2 (top row) illustrates averaged over 100 shots images of the electron
beam on the Lanex screen for coaxial propagation of femtosecond and nanosecond
pulses (a), and for D = 5 and 10 $\mu$m, (b) and (c) correspondingly. Of
course, from shot to shot, the position of the beam varies somewhat across the
detector, which is due to the instabilities of the transverse intensity
distribution across the focal spot of the nanosecond laser, the temporal
jitter of the Q-switch, which ultimately affect the breakdown parameters in
the tape. In the selected single shots (see second row in Fig.2), the beam
divergence does not exceed 0.1 rad. The maximum absolute electron beam
displacement on the Lanex corresponds to an angle of deflection from the
original axis by $\sim$10 degrees, Fig.2(c).

No evident decrement of average scintillation screen luminosity was observed
experimentally, which suggests that the beam charge did not vary
significantly. The image of the plasma channel scattered light follows the
electron beam as seen in the third row in Fig.2. It is worth noting that the
absolute position of the tape relative to the waist of the femtosecond beam
did not require adjustment in the case of a transverse displacement of the
breakdown axis. This greatly simplifies the handling of electron beam for
applications.

##  III Numerical simulations and discussion

The deflection angle of the laser beam $\theta$ during propagation along the
Oz axis is determined by the change in the refractive index $n$ along the
radial coordinate: $\dfrac{d\theta}{dz}=n^{-1}$$\dfrac{dn}{dr}$. Whereas,
according to the Drude model, the refractive index depends on the electron
density $n_{e}$ in plasma. The preformed plasma with radial density gradient
acts as a lens, leading to the refraction of the femtosecond pulse. Additional
modulation into the local refractive index (which causes, in particular, the
formation of a plasma channel) is introduced by field ionization and the
ponderomotive action of light during the relativistic laser pulse propagation
through the plasma. We made hydrodynamic simulations of the tape ablation by
the nanosecond laser pulse. Then PIC simulations of the electron acceleration
process using the simulated plasma profile were performed. PIC simulations
were made for different deviations D from the plasma longitudinal axis. The
target ablation and boring was modeled in hydrodynamic software package 3DLINE
Krukovskiy  _et al._ (2017). The one-temperature one-fluid model of a quasi-
neutral plasma was used. The local ionization multiplicity was determined from
the stationary distribution obtained in the framework of collisional-radiative
equilibrium using the THERMOS code Vichev  _et al._ (2019). This approximation
works well on the considered scale of characteristic times ($\sim$ns) and
laser intensities ($10^{12}-10^{13}$ W cm-2), and detailed consideration of
the ionization kinetics is not required. The equation of state was calculated
using the Thomas-Fermi model with semi-empirical corrections in the phase
transition region using the FEOS code Faik  _et al._ (2018). The target
material was PET (initial density 1.38 g/cm3, temperature 270 K), thickness 15
$\mu$m. The laser parameters where similar to the experiment: 1064 nm, 8 ns
pulse duration FWHM, peak intensity $10^{13}$ W cm-2. The simulation was
carried out in two-dimensional axially symmetric (r,z) geometry on a fixed
grid with a size of 0.5x0.945 mm. Fig.3 illustrates the resulting distribution
of ions concentration with strong radial dependence. The data for the optimal
delay was then integrated into the PIC code SMILEI Derouillat  _et al._
(2018). We first carried out 3D modeling to ensure the accurate reproduction
of the experimental conditions in the case of co-axial propagation of the two
pulses. The target initially consisted of neutral carbon. To save
computational time and resources the simulations with shifted propagation from
the symmetry axis were performed in simplified 2D3V geometry. The time step
was 0.083 fs, box size - X=200$\lambda$ (propagation axis), Y=80$\lambda$ in
2D and X=200$\lambda$, Y=Z=24$\lambda$ in 3D, laser wavelength $\lambda$=0.8
$\mu$m, spatial resolution - $\lambda$/32 in 2D and $\lambda$/32 (along X) and
$\lambda$/4 (along Y,Z) in 3D. Several transversal shifts 2.5-20 $\mu$m were
simulated. The 50 fs pulse was focused in a 4 $\mu$m spot FWHM at the front
surface of the non-perturbed tape with normalized vector potential
$a_{0}$=1.5, where $a_{0}\approx 0.85\sqrt{I\lambda^{2}}$, which corresponds
to the experimentally available peak intensity of $\sim 5\times 10^{18}$ W
cm-2.

Figure 3:  Calculated distributions (r,z) of the ions concentration in plasma
of the PET tape.

Without focusing too much attention on the mechanisms of electron acceleration
and interaction features (which can be found in Tsymbalov  _et al._ (2020)),
here we focus onto deflection of the femtosecond and electron beams while
passing through the plasma slab with radial dependence of the electron
density. It should be noted that this effect can be used not only to deflect a
laser beam (as well as an electron beam Mittelberger  _et al._ (2019)), but
also to diagnose the transverse plasma inhomogeneity itself Brockington  _et
al._ (2005); Enloe  _et al._ (1987). The principle difference between the 2D
and 3D simulations occurs in the modeling of self-focusing Pukhov and Meyer-
ter Vehn (1996) resulting in a longer lasting channel in the 3D case. Fig.4(a)
demonstrates the compared evolution of the magnetic field $B_{z}$ in units of
$a_{0}$ of the pulse in 3D/2D. The trailing part of the laser pulse diverges
stronger in the 2D model during the propagation in the channel. The particles
gain a finite energy in the channel within a few betatron oscillations and
then go out of synchronism, still being trapped by the channel walls though.
This affects the accelerated electron beam divergence: $\sim$0.5 rad in the 2D
and $\sim$0.1 rad in the 3D, evaluated as the FWHM of the angular distribution
of number of particles $\dfrac{dN}{d\theta}$, see Fig.4(b). Nevertheless, the
spectra $\dfrac{dN}{dE}$ in these two cases, Fig.4(c), reasonably match in
terms of shape (above 1 MeV) and maximal energy. Estimations of beam
divergence in 3D correlates well with experimental measurements of electron
beam.

Fig.5 represents a wider image of the beam propagation, ionization and
particles motion in the simulation box (2D) with the axial and shifted (D = 10
$\mu$m) geometry. The concentration gradient promotes the beam deflection. As
a result, the plasma channel also deviates from the original axis with
conservation of its properties. It is important to note that the laser beam
deflects mainly in the region of higher plasma density near the initial tape
position, where the electron concentration $n_{e}$ after ionization lies in
the range 0.1-0.2 units of critical density $n_{cr}$ (for 800 nm radiation
$n_{cr}=1.7\times 10^{21}$ cm-3). This is supported by the ray-tracing through
graded-index medium by solving the eikonal equation proposed by A.Sharma et
al. Sharma  _et al._ (1982):
$\dfrac{d}{ds}n(\vec{r})\dfrac{d}{ds}=\nabla{n(\vec{r})}$, where $\vec{r}$ is
the position vector of the ray relative to origin, $ds$ is the element of the
beam path, $n(\vec{r})$ represents the refractive index distribution. The
solution by 4th order Runge-Kutta method is shown by the dashed lines in
Fig.5. It does not consider self-focusing but clearly depicts the refraction
in the preformed plasma lens. The resulting effect in the bended channel is
the simultaneous deflection of the accelerated electron bunch.

Figure 4:  Comparison of magnetic field of the laser pulse in 3D and 2D at
different simulations periods (a). Normalized angular distribution of
particles (b) and energy spectra (c). Figure 5:  Results of numerical 2D PIC
calculations: normalized magnetic field in the box at the background of
normalized electron density during propagation of a femtosecond pulse with
deviation D = 0 (a) and 10 $\mu$m (b). Green dots mark electrons with energies
above 3 MeV. The blue dashed lines represent the results of ray tracing of a
laser beam through an inhomogeneous plasma.

The basic beam quality parameters are summarized in Fig.6 in dependence on the
deviation D of the driver pulse for two values of $a_{0}$. The laser pulse
begins to pass over the denser plasma region at D$\textgreater$10 $\mu$m,
where due to the modest peak intensity ($a_{0}$=1.5) the acceleration regime
related to channeling and DLA struggles to evolve efficiently, which is seen
as decrease of electrons energy, Fig.6(a). Our simulations with increased peak
intensity ($a_{0}$=2.2, Fig.6(b)) demonstrate that this effect is less
pronounced, ensuring experimental applicability of the proposed approach in a
wide range of laser intensity and final electrons energies. One can also
notice that charge of the bunch is almost the same independently of the
deviation D with $a_{0}$=2.2, Fig.6(c). As mentioned above the divergence in
the 2D simulations is higher as compared to the 3D and the experiment. However
the angular spread does not tend to widen in the investigated range of D being
in agreement with the experimental data. A similar effect, considered earlier
in the context of wakefield acceleration, led to a transverse inhomogeneity of
the plasma wave, the accelerating longitudinal and focusing transverse fields
of which deflected the accelerated electron beam following the laser driver
Nakajima (2018); Mittelberger  _et al._ (2019). The deviation of the
components of the laser field from the initial direction itself determines the
change in the trajectory of the accelerated electron. The effect of electron
beam small angular deflection in the wake was also discussed in terms of
transverse density gradient in gas jet outflow Doss  _et al._ (2023). Recently
the acceleration of electrons in the curved plasma channel was demonstrated
experimentally Zhu  _et al._ (2023). Within the framework of the plasma model
of an ablated and bored thin tape, one can conclude that key bunch parameters
(charge, spectrum, divergence) are retained acceptable and comparable to the
coaxial propagation in a wide range of deviations from the central axis
ensuring beam deflection up to 10 degrees. When further increasing D the laser
pulse hits denser layers, and the acceleration regime becomes affected.
Nevertheless, already small displacements evidently change the electron beam
propagation axis by several degrees, which is more than enough for the beam
deflection from its initial axis in whole range of particles energies without
the use of energy selective magnetic optics.

Figure 6:  Energy spectra of the electron beam at varied driver laser beam
deviation D from symmetry axis for the vector potential $a_{0}$=1.5 (a) and
$a_{0}$=2.2 (b). The relative charge of the electron beam for particles above
3 MeV in dependence on the D (c).

The summary of our studies with experimental data on the pulse and electron
beam deflection (evaluated from scattered light signal and Lanex screen
correspondingly) and comparison with numerical results are presented in Fig.7,
where a good correlation is observed. This indicates that the numerical model
of the tape target boring by the nanosecond prepulse is consistent with the
experiment. Moreover, one can rely on the good evaluation of the electron
bunch acceleration mechanism established numerically.

Figure 7:  Numerical and experimental data on the deflection of the electron
beam (a) and the laser pulse (b) as a function of the femtosecond pulse
deviation D relative to the breakdown symmetry axis.

##  IV Conclusion

The angular manipulation of a laser-generated electron beam during its
acceleration in a relativistic plasma channel in the presence of a radial
electron density gradient is experimentally shown for the first time. An
inhomogeneous transverse density distribution was formed under the action of a
powerful nanosecond pulse ($10^{13}$ W cm-2) on a thin mylar tape (15 $\mu$m).
Heating and ablation lead to the tape boring with the formation of a plasma
with a strong radial dependence of the concentration. The propagation of a
relativistically intense femtosecond pulse (above $10^{18}$ W cm-2, 50 fs)
through a plasma with a slight deviation relative to the breakdown symmetry
axis is accompanied by refraction of the laser beam on the density gradient.
It has been numerically established that, within the framework of the
considered model of the breakdown and boring of a tape target by nanosecond
radiation, the deviation of the propagation axis of a femtosecond pulse from
the breakdown symmetry axis of up to 10 $\mu$m leads to deflection of the
generated electron beam by up to 10 degrees, while the beam parameters
maintain almost unchanged. The calculations are confirmed by experimental
data, where a similar deflection of the electron beam is observed on the
detector.

The experimental approach can be used to create schemes for multistage
electron acceleration in a sequence of laser accelerators. Of great importance
is the possibility of in-situ control of beam directionality for the
multistage accelerators, where the electron bunch needs to be injected at a
small angle to the axis Starodubtseva  _et al._ (2023). In the proposed
approach there is no angular dispersion in the energy of particles in the
beam, which preserves a high charge of the beam in a collimated form in
contrast to the case of using magnetic optics. Our representation, for
example, for the AWAKE experiment the scheme could be drawn as in Fig.8. The
driver femtosecond pulse provides seed electrons with 10-20 MeV energy. The
tape target may serve as the entrance window of the gas cell. It is also worth
noting that the proposed method can serve not only for the tasks described
above, but also for diagnosing the interaction region: an optically opaque
volume of dense plasma near the tape, the channel formed by a nanosecond
pulse, and the radial density distribution in it can be indirectly estimated
from angle of refraction of the escaping beam.

Figure 8:  Simplified model of the implementation of electron beam deflection
approach to the AWAKE experiment.

###### Acknowledgements.

The work was carried out with the support of the Russian Science Foundation
project 21-79-10207. A part of equipment was purchased with the support of the
Ministry of Education and Science of the Russian Federation within the
framework of the national project ”Science and universities”. The
applicability of our approach for multi stage electrons acceleration and AWAKE
experiment was considered in the frame of scientific program of the National
Center of Physics and Mathematics (project “Physics of high energy density.
Stage 2023-2025”).

## References

  * Malka (2012) V. Malka, Physics of Plasmas 19, 055501 (2012).
  * Joshi (2019) C. Joshi, Plasma Physics and Controlled Fusion 61, 104001 (2019).
  * Ledingham and Galster (2010) K. W. D. Ledingham and W. Galster, New Journal of Physics 12, 045005 (2010).
  * Kurz  _et al._ (2021) T. Kurz, T. Heinemann, M. F. Gilljohann, Y. Y. Chang, J. P. Couperus Cabadağ, A. Debus, O. Kononenko, R. Pausch, S. Schöbel, R. W. Assmann, M. Bussmann, H. Ding, J. Götzfried, A. Köhler, G. Raj, S. Schindler, K. Steiniger, O. Zarini, S. Corde, A. Döpp, B. Hidding, S. Karsch, U. Schramm, A. Martinez de la Ossa, and A. Irman, Nature Communications 12, 2895 (2021).
  * Hussein  _et al._ (2021) A. E. Hussein, A. V. Arefiev, T. Batson, H. Chen, R. S. Craxton, A. S. Davies, D. H. Froula, Z. Gong, D. Haberberger, Y. Ma, P. M. Nilson, W. Theobald, T. Wang, K. Weichman, G. J. Williams, and L. Willingale, New Journal of Physics 23, 023031 (2021).
  * Arefiev  _et al._ (2016) A. V. Arefiev, V. N. Khudik, A. P. L. Robinson, G. Shvets, L. Willingale, and M. Schollmeier, Physics of Plasmas 23, 056704 (2016).
  * Takabe and Kuramitsu (2021) H. Takabe and Y. Kuramitsu, High Power Laser Science and Engineering 9, e49 (2021).
  * Nedorezov  _et al._ (2021) V. G. Nedorezov, S. G. Rykovanov, and A. B. Savel’ev, Physics-Uspekhi 64, 1214 (2021).
  * Nakajima (2015) K. Nakajima, Proceedings of the Japan Academy, Series B 91, 223 (2015).
  * Labate  _et al._ (2020) L. Labate, D. Palla, D. Panetta, F. Avella, F. Baffigi, F. Brandi, F. Di Martino, L. Fulgentini, A. Giulietti, P. Köster, D. Terzani, P. Tomassini, C. Traino, and L. A. Gizzi, Scientific Reports 10, 17307 (2020).
  * Antipenkov  _et al._ (2021) R. Antipenkov, R. Boge, E. Erdman, M. Greco, J. T. Green, A. Grenfell, M. Horáček, Z. Hubka, L. Indra, K. Majer, T. Mazanec, P. Mazůrek, J. Naylon, J. Novák, V. Šobr, A. Špaček, P. Strkula, W. Szuba, B. Tykalewicz, P. Bakule, and B. Rus, in  _High Power Lasers and Applications_, Vol. 11777, edited by J. Hein, T. J. Butcher, P. Bakule, C. L. Haefner, G. Korn, and L. O. Silva, International Society for Optics and Photonics (SPIE, 2021) p. 117770E.
  * Kim  _et al._ (2013) H. T. Kim, K. H. Pae, H. J. Cha, I. J. Kim, T. J. Yu, J. H. Sung, S. K. Lee, T. M. Jeong, and J. Lee, Phys. Rev. Lett. 111, 165002 (2013).
  * Zhu  _et al._ (2022) X.-Z. Zhu, M. Chen, B.-Y. Li, F. Liu, X.-L. Ge, Z.-M. Sheng, and J. Zhang, Physics of Plasmas 29, 013101 (2022).
  * Jin  _et al._ (2019) Z. Jin, H. Nakamura, N. Pathak, Y. Sakai, A. Zhidkov, K. Sueda, R. Kodama, and T. Hosokai, Scientific Reports 9, 20045 (2019).
  * Massimo  _et al._ (2019) F. Massimo, A. Beck, J. Derouillat, M. Grech, M. Lobet, F. Pérez, I. Zemzemi, and A. Specka, Plasma Physics and Controlled Fusion 61, 124001 (2019).
  * Starodubtseva  _et al._ (2023) E. M. Starodubtseva, I. N. Tsymbalov, K. A. Ivanov, D. A. Gorlova, and A. B. Savel’ev, Bulletin of the Lebedev Physics Institute 50, S735 (2023).
  * Khudiakov and Pukhov (2022) V. Khudiakov and A. Pukhov, Phys. Rev. E 105, 035201 (2022).
  * Moschuering  _et al._ (2019) N. Moschuering, K. V. Lotov, K. Bamberg, F. Deutschmann, and H. Ruhl, Plasma Physics and Controlled Fusion 61, 104004 (2019).
  * Favaudon  _et al._ (2014) V. Favaudon, L. Caplier, V. Monceau, F. Pouzoulet, M. Sayarath, C. Fouillade, M.-F. Poupon, I. Brito, P. Hupé, J. Bourhis, J. Hall, J.-J. Fontaine, and M.-C. Vozenin, Science Translational Medicine 6, 245ra93 (2014).
  * Gahn  _et al._ (1999) C. Gahn, G. D. Tsakiris, A. Pukhov, J. Meyer-ter Vehn, G. Pretzler, P. Thirolf, D. Habs, and K. J. Witte, Phys. Rev. Lett. 83, 4772 (1999).
  * Nakajima (2018) K. Nakajima, Light: Science & Applications 7, 21 (2018).
  * Luo  _et al._ (2018) J. Luo, M. Chen, W. Y. Wu, S. M. Weng, Z. M. Sheng, C. B. Schroeder, D. A. Jaroszynski, E. Esarey, W. P. Leemans, W. B. Mori, and J. Zhang, Phys. Rev. Lett. 120, 154801 (2018).
  * Mittelberger  _et al._ (2019) D. E. Mittelberger, M. Thévenet, K. Nakamura, A. J. Gonsalves, C. Benedetti, J. Daniels, S. Steinke, R. Lehe, J.-L. Vay, C. B. Schroeder, E. Esarey, and W. P. Leemans, Phys. Rev. E 100, 063208 (2019).
  * qing Zhu  _et al._ (2020) C. qing Zhu, J. guang Wang, Y. fei Li, J. Feng, D. zhang Li, Y. hang He, J. hao Tan, J. long Ma, X. Lu, Y. tong Li, and L. ming Chen, Opt. Express 28, 11609 (2020).
  * Tsymbalov  _et al._ (2020) I. Tsymbalov, D. Gorlova, K. Ivanov, S. Shulyapov, V. Prokudin, A. Zavorotny, R. Volkov, V. Bychenkov, V. Nedorezov, and A. Savel’ev, Plasma Physics and Controlled Fusion 63, 022001 (2020).
  * Krukovskiy  _et al._ (2017) A. Y. Krukovskiy, V. G. Novikov, and I. P. Tsygvintsev, Mathematical Models and Computer Simulations 9, 48 (2017).
  * Vichev  _et al._ (2019) I. Vichev, A. Solomyannaya, A. Grushin, and D. Kim, High Energy Density Physics 33, 100713 (2019).
  * Faik  _et al._ (2018) S. Faik, A. Tauschwitz, and I. Iosilevskiy, Computer Physics Communications 227, 117 (2018).
  * Derouillat  _et al._ (2018) J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, Computer Physics Communications 222, 351 (2018).
  * Brockington  _et al._ (2005) S. J. E. Brockington, R. D. Horton, D. Q. Hwang, R. W. Evans, S. J. Howard, and Y. C. F. Thio, Review of Scientific Instruments 76, 063503 (2005).
  * Enloe  _et al._ (1987) C. L. Enloe, R. M. Gilgenbach, and J. S. Meachum, Review of Scientific Instruments 58, 1597 (1987).
  * Pukhov and Meyer-ter Vehn (1996) A. Pukhov and J. Meyer-ter Vehn, Phys. Rev. Lett. 76, 3975 (1996).
  * Sharma  _et al._ (1982) A. Sharma, D. V. Kumar, and A. K. Ghatak, Appl. Opt. 21, 984 (1982).
  * Doss  _et al._ (2023) C. E. Doss, R. Ariniello, J. R. Cary, S. Corde, H. Ekerfelt, E. Gerstmayr, S. J. Gessner, M. Gilljohann, C. Hansel, B. Hidding, M. J. Hogan, A. Knetsch, V. Lee, K. Marsh, B. O’Shea, P. San Miguel Claveria, D. Storey, A. Sutherland, C. Zhang, and M. D. Litos, Phys. Rev. Accel. Beams 26, 031302 (2023).
  * Zhu  _et al._ (2023) X. Zhu, B. Li, F. Liu, J. Li, Z. Bi, X. Ge, H. Deng, Z. Zhang, P. Cui, L. Lu, W. Yan, X. Yuan, L. Chen, Q. Cao, Z. Liu, Z. Sheng, M. Chen, and J. Zhang, Phys. Rev. Lett. 130, 215001 (2023).