# Kinematical higher-twist corrections in $\gamma^{*}\gamma\to M\bar{M}$

Cédric Lorcé<EMAIL_ADDRESS>CPHT, CNRS, Ecole Polytechnique,
Institut Polytechnique de Paris, Route de Saclay, 91128 Palaiseau, France
Bernard Pire<EMAIL_ADDRESS>CPHT, CNRS, Ecole Polytechnique,
Institut Polytechnique de Paris, Route de Saclay, 91128 Palaiseau, France
Qin-Tao Song<EMAIL_ADDRESS>CPHT, CNRS, Ecole Polytechnique, Institut
Polytechnique de Paris, Route de Saclay, 91128 Palaiseau, France  School of
Physics and Microelectronics, Zhengzhou University, Zhengzhou, Henan 450001,
China

###### Abstract

We estimate kinematical higher-twist (up to twist 4) corrections to the
$\gamma^{*}(q_{1})\gamma(q_{2})\to M(p_{1})\bar{M}(p_{2})$ amplitudes at large
$Q^{2}=-q_{1}^{2}$ and small $s=(q_{1}+q_{2})^{2}$, where $M$ is a scalar or
pseudoscalar meson. This process is known to factorize at leading twist into a
perturbatively calculable coefficient function and generalized distribution
amplitudes (GDAs). The kinematical higher-twist contributions of order
$s/Q^{2}$ and $m^{2}/Q^{2}$ turn out to be important in the cross section,
considering the kinematics accessible at Belle and Belle II. We present
numerical estimates for the cross section for
$\gamma^{*}\gamma\to\pi^{0}\pi^{0}$ with the $\pi\pi$ GDA extracted from Belle
measurements and with the asymptotic $\pi\pi$ GDA as inputs to study the
magnitude of the kinematical corrections. To see how the target mass
corrections of order $m^{2}/Q^{2}$ affect the cross section, we also perform
the calculation for $\gamma^{*}\gamma\to\eta\eta$ by using a model $\eta\eta$
GDA. In the range $s>1$ GeV2, the kinematical higher-twist corrections account
for $\sim 15\%$ of the total cross section, an effect which is not negligible.
Since $\pi\pi$ GDAs are the best way to access the pion energy-momentum tensor
(EMT), our study demonstrates that an accurate evaluation of EMT form factors
requires the inclusion of kinematical higher-twist contributions.

##  I Introduction

Generalized distribution amplitudes (GDAs) Muller:1994ses ; Diehl:1998dk ;
Polyakov:1998ze  – sometimes called two-meson distribution amplitudes – are
hadronic matrix elements closely related to generalized parton distributions
(GPDs) Diehl:2003ny ; Belitsky:2005qn ; Boffi:2007yc ; Goeke:2001tz . They
involve the same bilocal quark (or gluon) operator on the light cone and
correspond to $s$-$t$ crossed helicity matrix elements. GDAs can be accessed
in $e(k_{1})\gamma\to e^{\prime}(k_{2})M(p_{1})\bar{M}(p_{2})$ reactions in
$e^{+}e^{-}$ collisions, in the kinematical range where
$Q^{2}=-(k_{1}-k_{2})^{2}$ is large but $s=(p_{1}+p_{2})^{2}$ is much smaller
than $Q^{2}$. They have already been the subject of careful studies at Belle
Belle:2015oin  and were extracted in a leading-twist analysis in Ref.
Kumano:2017lhr . They are also important in the understanding of heavy meson
three-body decays, in particular in the quest for a precise determination of
Cabibbo-Kobayashi-Maskawa (CKM) matrix elements Chen:2002th ; Wang:2015uea ;
Li:2016tpn ; Jia:2021uhi .

As the studies of GPDs allow us to perform nucleon tomography through a
Fourier transform in the transverse coordinate space Burkardt:2000za ;
Ralston:2001xs ; Diehl:2002he , GDAs open the way to an impact-parameter
picture Pire:2002ut  of the exclusive hadronization process $q\bar{q}\to
M(p_{1})\bar{M}(p_{2})$. The GPDs and GDAs are also used to investigate the
matrix elements of the energy-momentum tensor (EMT) Ji:1996ek ; Ji:1996nm ;
Diehl:2000uv ; Polyakov:1998ze  for hadrons in the spacelike and timelike
regions, respectively. One can extract mass, pressure and shear force
distributions of hadrons with the spacelike EMT form factors Polyakov:2002yz ;
Goeke:2007fp ; Mai:2012cx ; Polyakov:2018zvc ; Burkert:2018bqq ;
Kumericki:2019ddg ; Lorce:2018egm ; Dutrieux:2021nlz . Since there is no
experimental facility where pion GPDs can be directly measured (see however
Refs. Amrath:2008vx ; Chavez:2021koz ; Chavez:2021llq ), the studies of
$\pi\pi$ GDAs are a necessary tool to access the pion EMT. The spacelike pion
EMT form factors can be obtained from the timelike ones by using dispersion
relations, and in this process the $\pi\pi$ GDAs and pion EMT form factors at
$s>1$ GeV2 should be included so as to make the integrals convergent.
Therefore, this goal necessitates to extract GDAs in a sufficiently large
$s$-range, thus demanding a control as precise as possible of kinematical
higher-twist corrections to the amplitudes, which are proportional to
$s/Q^{2}$ and $m^{2}/Q^{2}$, with $m$ the meson mass. While the contribution
of one higher-twist process has been previously discussed in Ref.
Lansberg:2006fv , it only matters in a very limited kinematical region, namely
the near forward or near backward regions. The phenomenological necessity of a
sizeable (genuine) twist-$4$ contribution to the $\gamma^{*}\gamma\to\rho\rho$
amplitude has also been pointed out Anikin:2005ur . A complete understanding
of higher-twist corrections to the process

$\displaystyle\gamma^{*}(q_{1})\gamma(q_{2})\to M(p_{1})\bar{M}(p_{2})\,,$ (1)

is however a difficult task which is far from being achieved.

Meanwhile, a separation of kinematical and dynamical contributions in the
product of two electromagnetic currents
$T\\{j_{\mu}^{\text{em}}(z_{1}x)j_{\nu}^{\text{em}}(z_{2}x)\\}$ was proven in
Refs. Braun:2011dg ; Braun:2011zr ; Braun:2011th  and applied to the deeply-
virtual Compton scattering (DVCS) reaction Braun:2012bg . The kinematical
corrections come from two types of operators, namely the subtraction of traces
in the leading-twist operators and the higher-twist operators which can be
reduced to the total derivatives of the leading-twist ones. The subtraction of
traces was applied in Ref. Nachtmann:1973mr  to the reaction of Deep Inelastic
Scattering (DIS), leading to target mass corrections. The kinematical
corrections in DVCS can be considered as a generalization of these target mass
corrections. However, higher-twist operators which can be reduced to the total
derivatives of the leading-twist ones will also contribute to the DVCS
reaction, since nonforward matrix elements are used. As pointed out in Refs.
Braun:2011dg ; Braun:2011zr ; Braun:2011th , the distinction between two types
of kinematical corrections is not Lorentz invariant and has no physical
meaning. Both contributions should therefore better always be added together.
Since the same operator governs the physics of the reaction (1), one may use
the same techniques to improve our understanding of the $s$-dependence of its
amplitude. We thus study here the kinematical higher-twist corrections to the
amplitude of the reaction (1), in the kinematical domain suitable for a
collinear QCD factorization framework where the leading-twist amplitude can be
written as the convolution of a perturbatively calculable coefficient function
and GDAs Diehl:2000uv .

In Sect. II, we describe the kinematics of the $\gamma^{*}\gamma\to M\bar{M}$
process and recall the basic properties of GDAs. In Sect. III, we recall the
results of Refs. Braun:2011dg ; Braun:2011zr ; Braun:2011th  and the
definitions of the higher-twist kinematical operators. In Sect. IV, we derive
the helicity amplitudes for the reaction (1), including the kinematical
higher-twist contributions. Sect. V shows our numerical estimates of the
kinematical higher-twist contributions to the cross section for both $\pi\pi$
and $\eta\eta$ cases. We briefly present our conclusions in Sect. VI.
Appendices A and B provide technical details for the calculation of helicity
amplitudes.

##  II Kinematics and generalized distribution amplitudes

Figure 1: Kinematics of the process $\gamma^{*}(q_{1})\gamma(q_{2})\to
M(p_{1})\bar{M}(p_{2})$ in the center of mass of the meson pair; the virtual
photon is emitted by the electron with four-momentum $q_{1}=k_{1}-k_{2}$.

To describe the process (1), we define the lightlike vectors $n$ and
$\tilde{n}$ in a convenient way so that they can be expressed by the momenta
of the spacelike virtual photon $q_{1}$ and the real photon $q_{2}$,

$\displaystyle\tilde{n}=q_{1}+(1-\tau)q_{2},\qquad n=q_{2},$ (2)

where $\tau=s/(Q^{2}+s)$, $Q^{2}=-q_{1}^{2}$, and
$s=(q_{1}+q_{2})^{2}=(p_{1}+p_{2})^{2}$. The polar angle of the meson ($M$)
momenta $\theta$ is illustrated in Fig. 1, and is defined as

$\displaystyle\cos{\theta}=\frac{2q_{1}\cdot(p_{2}-p_{1})}{\beta_{0}\,(Q^{2}+s)},\qquad\beta_{0}=\sqrt{1-\frac{4m^{2}}{s}},$
(3)

where $m$ is the meson mass. For convenience, a new variable
$\zeta_{0}=-\beta_{0}\cos{\theta}$ is introduced instead of $\cos\theta$

$\displaystyle\zeta_{0}=\frac{(p_{2}-p_{1})\cdot n}{(p_{2}+p_{1})\cdot n},$
(4)

but the final amplitudes will be expressed in terms of $\cos{\theta}$. If the
$z$-axis is chosen so that the $t$-$z$ plane contains the lightlike vectors
$n$ and $\tilde{n}$, then only $\Delta=p_{2}-p_{1}$ has a transverse momentum,
$\Delta=\zeta_{0}(\tilde{n}-\tau n)+\Delta_{T}$. Using the on-shell condition,
we obtain $\Delta_{T}^{2}=4m^{2}-(1-\zeta_{0}^{2})s$.

The amplitude for $\gamma^{*}\gamma\to M\bar{M}$ is defined as

$\displaystyle A_{\mu\nu}=i\int d^{4}x\,e^{-ir\cdot
x}\langle\bar{M}(p_{2})M(p_{1})|\,T\\{j_{\mu}^{\text{em}}(z_{1}x)j_{\nu}^{\text{em}}(z_{2}x)\\}\,|0\rangle,\\!\\!$
(5)

where $r=z_{1}q_{1}+z_{2}q_{2}$, and the constraint $z_{1}-z_{2}=1$ is imposed
for real constants $z_{1}$ and $z_{2}$. Owing to the electromagnetic gauge
invariance, one can decompose this amplitude as Braun:2012bg

$\displaystyle
A^{\mu\nu}=-A^{(0)}\,g_{\perp}^{\mu\nu}+A^{(1)}\,\frac{\Delta_{\alpha}g_{\perp}^{\alpha\nu}}{Q}\left(\tilde{n}^{\mu}+(1-\tau)n^{\mu}\right)+\frac{1}{2}\,A^{(2)}\,\Delta_{\alpha}\Delta_{\beta}(g_{\perp}^{\alpha\mu}g_{\perp}^{\beta\nu}-\epsilon_{\perp}^{\alpha\mu}\epsilon_{\perp}^{\beta\nu})+A^{(3)\mu}\,n^{\nu}$
(6)

with $g_{\perp}^{\mu\nu}$ and $\epsilon_{\perp}^{\mu\nu}$ given by

$\displaystyle
g_{\perp}^{\mu\nu}=g^{\mu\nu}-\frac{n^{\mu}\tilde{n}^{\nu}+n^{\nu}\tilde{n}^{\mu}}{n\cdot\tilde{n}},\qquad\epsilon_{\perp}^{\mu\nu}=\epsilon^{\mu\nu\alpha\beta}\,\frac{\tilde{n}_{\alpha}n_{\beta}}{n\cdot\tilde{n}}.$
(7)

The last term in Eq. (6) is of no interest since it does not contribute to any
observable, and the rest of them can be expressed in terms of the GDAs if the
factorization conditions $Q^{2}\gg s,\Lambda_{\text{QCD}}^{2}$ are satisfied.
The leading-twist amplitude was first presented in Ref. Diehl:2000uv  with the
help of a twist-2 GDA $\Phi_{q}(z,\zeta_{0},s)$ for an isoscalar meson pair,

$\displaystyle\langle\bar{M}(p_{2})M(p_{1})|\,\bar{q}(z_{1}n)\not{n}q(z_{2}n)\,|0\rangle=2P\cdot
n\int dz\,e^{2i\left[zz_{1}+(1-z)z_{2}\right]P\cdot
n}\,\Phi_{q}(z,\zeta_{0},s),$ (8)

where $P=(p_{1}+p_{2})/2$, $\Phi_{q}$ is the GDA for the quark flavor $q$,
$\bar{q}(z_{1}n)\not{n}q(z_{2}n)$ is the leading-twist vector operator (a
light-like Wilson line joining the points $z_{1}n$ and $z_{2}n$ is implied),
and $z_{1}-z_{2}=1$ is not a necessary condition. This matrix element can
alternatively be expressed in terms of double distributions (DDs) as
Teryaev:2001qm

$\displaystyle\langle\bar{M}(p_{2})M(p_{1})|\,\bar{q}(z_{1}n)\not{n}q(z_{2}n)\,|0\rangle=\int
d\beta\,d\alpha\left[f_{q}(\beta,\alpha)\,\Delta\cdot
n-g_{q}(\beta,\alpha)\,2P\cdot n\right]e^{-il_{z_{1}z_{2}}\cdot n}$ (9)

with $f_{q}$ and $g_{q}$ having support on the rhombus $|\alpha|+|\beta|\leq
1$ and assumed to vanish at the boundary, and

$\displaystyle
l_{z_{1}z_{2}}=(z_{2}-z_{1})\left[\beta\,\frac{\Delta}{2}-(\alpha+1)P\right]-2z_{1}P.$
(10)

Then, one can easily relate the GDA to double distributions

$\displaystyle\Phi_{q}(z,\zeta_{0},s)=2\int
d\beta\,d\alpha\,\delta(y+\alpha-\beta\zeta_{0})\left[f_{q}(\beta,\alpha)\,\zeta_{0}-g_{q}(\beta,\alpha)\right],$
(11)

where $y=2z-1$. Since the meson pair is produced with charge conjugation
$\mathcal{C}=+1$, one can obtain the relations

$\displaystyle f_{q}(\beta,\alpha)$
$\displaystyle=f_{q}(\beta,-\alpha),\qquad\phantom{--}g_{q}(\beta,\alpha)=-g_{q}(\beta,-\alpha),$
(12) $\displaystyle f_{q}(\beta,\alpha)$
$\displaystyle=-f_{q}(-\beta,-\alpha),\qquad
g_{q}(\beta,\alpha)=-g_{q}(-\beta,-\alpha)$

from charge conjugation invariance. Assuming that the DDs vanish at the
boundaries, Eq. (9) can be put in the form

$\langle\bar{M}(p_{2})M(p_{1})|\,\bar{q}(z_{1}n)\not{n}q(z_{2}n)\,|0\rangle=\frac{2i}{z_{12}}\int
d\beta\,d\alpha\,\phi_{q}(\beta,\alpha)\,e^{-il_{z_{1}z_{2}}\cdot n},$ (13)

where the notation $z_{12}=z_{1}-z_{2}$ is used. A new distribution

$\displaystyle\phi_{q}(\beta,\alpha)=\partial_{\beta}f_{q}(\beta,\alpha)+\partial_{\alpha}g_{q}(\beta,\alpha)$
(14)

is introduced with symmetry
$\phi_{q}(\beta,\alpha)=\phi_{q}(\beta,-\alpha)=\phi_{q}(-\beta,-\alpha)$, in
order to simplify the calculation of the amplitudes thanks to the property

$\displaystyle\int
d\beta\,d\alpha\,\phi_{q}(\beta,\alpha)\left[a+b\,\alpha^{n}+c\,\beta^{m}\right]=0,$
(15)

where $a$, $b$ and $c$ are constants which are independent of $\alpha$ and
$\beta$, and the exponents $n$ and $m$ are odd numbers. Although the
intermediate calculations involve the DD $\phi_{q}(\alpha,\beta)$, the final
results will be presented in terms of the GDA using

$\displaystyle\frac{\partial\Phi_{q}(z,\zeta_{0},s)}{\partial z}=4\int
d\beta\,d\alpha\,\delta((2z-1)+\alpha-\beta\zeta_{0})\,\phi_{q}(\beta,\alpha).$
(16)

##  III Operator product expansion and helicity amplitudes

A separation of kinematical and dynamical contributions in the time-ordered
product of two electromagnetic currents
$i\,T\\{j_{\mu}^{\text{em}}(z_{1}x)j_{\nu}^{\text{em}}(z_{2}x)\\}$ was
recently proved in Refs. Braun:2011dg ; Braun:2011zr ; Braun:2011th . The
kinematical contributions only involve the leading-twist distributions,
whereas unrelated genuine higher-twist distributions are necessary for the
dynamical contributions. One can thus improve the description of reactions
where two photons are involved by including the kinematical corrections,
without any knowledge of the higher-twist distributions. A complete
calculation of kinematical corrections was performed up to the twist-4
accuracy for DVCS with a (pseudo)scalar target in Ref. Braun:2012bg . In this
work we shall apply similar techniques to calculate the kinematical higher-
twist contributions in the $s$-$t$ crossed channel of DVCS, namely the
reaction $\gamma^{*}\gamma\to M\bar{M}$. The kinematical contributions to the
operator $i\,T\\{j_{\mu}^{\text{em}}(z_{1}x)j_{\nu}^{\text{em}}(z_{2}x)\\}$
were given to twist-4 accuracy by Braun:2012bg ; Braun:2011dg ; Braun:2011zr ,

$\displaystyle
T_{\mu\nu}=\frac{-1}{\pi^{2}x^{4}z_{12}^{3}}\left\\{x^{\alpha}\left[S_{\mu\alpha\nu\beta}\mathbb{V}^{\beta}-i\epsilon_{\mu\alpha\nu\beta}\mathbb{W}^{\beta}\right]+x^{2}\left[(x_{\mu}\partial_{\nu}+x_{\nu}\partial_{\mu})\mathbb{X}+(x_{\mu}\partial_{\nu}-x_{\nu}\partial_{\mu})\mathbb{Y}\right]\right\\},$
(17)

where the convention $\epsilon_{0123}=1$ is adopted for the antisymmetric
tensor and $S^{\mu\alpha\nu\beta}$ is defined as

$\displaystyle
S^{\mu\alpha\nu\beta}=g^{\mu\alpha}g^{\nu\beta}-g^{\mu\nu}g^{\alpha\beta}+g^{\mu\beta}g^{\nu\alpha}.$
(18)

In Eq. (17), $\mathbb{V}_{\mu}$ and $\mathbb{W}_{\mu}$ contain contributions
of twist 2, twist 3 and twist 4, whereas $\mathbb{X}$ and $\mathbb{Y}$ are
purely twist 4, see Appendix A for the detailed expressions. In practice, the
spinor formalism Braun:2008ia ; Braun:2009vc  is used to calculate the
amplitudes, since the expression of $T_{\mu\nu}$ becomes more compact and it
is easier to figure out the twist of each term in the corresponding matrix
elements.

In order to calculate the helicity amplitudes of Eq. (6), the photon
polarization vectors are required. Choosing the momentum of the virtual photon
along the z-axis, as shown in Fig. 1, its polarization vectors read
Diehl:2000uv

$\displaystyle\epsilon_{0}^{\mu}=\frac{1}{Q}(|q_{1}^{3}|,0,0,q_{1}^{0}),\quad\epsilon_{\pm}^{\mu}=\frac{1}{\sqrt{2}}(0,\mp
1,-i,0),$ (19)

where the lower indices $\pm$ and $0$ indicate the helicities of the photon.
The polarization vectors $\tilde{\epsilon}$ of the real photon only have the
transverse components, and they are related to the ones of the virtual photon
as $\tilde{\epsilon}_{\pm}=-\epsilon_{\mp}$. In the reaction
$\gamma^{*}\gamma\to M\bar{M}$, the helicity amplitudes are defined as

$\displaystyle A_{ij}=\epsilon_{i}^{\mu}\tilde{\epsilon}_{j}^{\nu}A_{\mu\nu},$
(20)

and there are only three independent helicity amplitudes owing to parity
invariance, as one can check from Eq. (6). Here we choose the independent
helicity amplitudes as $A_{++}$, $A_{0+}$ and $A_{-+}$, then one obtains

$\displaystyle A_{++}=A_{--}=A^{(0)},\quad
A_{0+}=-A^{(1)}(\Delta\cdot\epsilon_{-}),\quad
A_{-+}=-A^{(2)}(\Delta\cdot\epsilon_{-})^{2}.$ (21)

At leading twist, the operator product expansion of
$i\,T\\{j_{\mu}^{\text{em}}(z_{1}x)j_{\nu}^{\text{em}}(z_{2}x)\\}$ leads to
the nonlocal operator

$\displaystyle
O_{++}(z_{1}n,z_{2}n)=\sum_{q}e_{q}^{2}\,\bar{q}(z_{1}n)\not{n}q(z_{2}n)$ (22)

with a lightlike separation. Since we are interested in the reaction
$\gamma^{\ast}\gamma\rightarrow M\bar{M}$ with a charge conjugation even final
state, one can safely neglect the contribution of the strange quark in the
case of a $\pi$ meson pair,

$\displaystyle
O_{++}(z_{1}n,z_{2}n)=e_{u}^{2}\,\bar{u}(z_{1}n)\not{n}u(z_{2}n)+e_{d}^{2}\,\bar{d}(z_{1}n)\not{n}d(z_{2}n)=\chi\left[\bar{u}(z_{1}n)\not{n}u(z_{2}n)+\bar{d}(z_{1}n)\not{n}d(z_{2}n)\right],$
(23)

where $\chi=5e^{2}/18$ is obtained thanks to the isospin symmetry111The
amplitudes associated with $\bar{u}(z_{1}n)\not{n}u(z_{2}n)$ and
$\bar{d}(z_{1}n)\not{n}d(z_{2}n)$ are the same for an isosinglet $\pi\pi$
state.. One needs however to add $e_{s}^{2}\,\bar{s}(z_{1}n)\not{n}s(z_{2}n)$
to the operator $O_{++}(z_{1}n,z_{2}n)$ in the case of a $K$ meson pair. The
kinematical higher-twist contributions in the operator product expansion of
$i\,T\\{j_{\mu}^{\text{em}}(z_{1}x)j_{\nu}^{\text{em}}(z_{2}x)\\}$ are related
to the operator $\mathcal{O}_{++}^{t=2}(z_{1},z_{2})$, where the separation
$x$ is now not necessarily lightlike. We thus need to use the leading-twist
projector $\Pi(x,n)$ defined in Refs. Braun:2011dg ; Braun:2011zr ;
Braun:2011th ,

$\displaystyle\langle\bar{M}(p_{2})M(p_{1})|\,\mathcal{O}_{++}^{t=2}(z_{1},z_{2})\,|0\rangle=\Pi(x,n)\langle\bar{M}(p_{2})M(p_{1})|\,O_{++}(z_{1}n,z_{2}n)\,|0\rangle.$
(24)

Since the dependence on $n$ is always carried by a function of the type
$e^{-il\cdot n}$ in Eq. (13), the action of the leading-twist projector is
simply given by

$\displaystyle\left[\Pi e^{-il\cdot n}\right](x)=e^{-il\cdot
x}+\frac{x^{2}l^{2}}{4}\int_{0}^{1}dv\,v\,e^{-ivl\cdot x}+\mathcal{O}(x^{4}).$
(25)

Up to $1/Q^{2}$-accuracy, one obtains

$\displaystyle\langle\bar{M}(p_{2})M(p_{1})|\,\mathcal{O}_{++}^{t=2}(z_{1},z_{2})\,|0\rangle=\chi\,\frac{2i}{z_{12}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left[e^{-il_{z_{1}z_{2}}\cdot
x}+\frac{x^{2}l_{z_{1}z_{2}}^{2}}{4}\int_{0}^{1}dv\,v\,e^{-ivl_{z_{1}z_{2}}\cdot
x}\right],$ (26)

where $\phi=\phi_{u}+\phi_{d}$ and the second term provides a twist-4
contribution. In addition to the leading-twist operator
$\mathcal{O}_{++}^{t=2}(z_{1}x,z_{2}x)$, there are also the higher-twist
operators

$\displaystyle\mathcal{O}_{1}(z_{1},z_{2})$
$\displaystyle=\left[i\mathbf{P}^{\mu},\,\left[i\mathbf{P}_{\mu},\,\mathcal{O}_{++}^{t=2}(z_{1},z_{2})\right]\right],$
$\displaystyle\mathcal{O}_{2}(z_{1},z_{2})$
$\displaystyle=\left[i\mathbf{P}^{\mu},\,\frac{\partial}{\partial
x^{\mu}}\mathcal{O}_{++}^{t=2}(z_{1},z_{2})\right],$ (27)

which contribute to kinematical higher-twist corrections. Using Eq. (26), the
matrix elements of $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ can be expressed up
to $1/Q^{2}$-accuracy as

$\displaystyle\langle\bar{M}(p_{2})M(p_{1})|\,\mathcal{O}_{1}(z_{1},z_{2})\,|0\rangle$
$\displaystyle=-\chi\,\frac{2i}{z_{12}}\,s\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,e^{-il_{z_{1}z_{2}}\cdot x},$
$\displaystyle\langle\bar{M}(p_{2})M(p_{1})|\,\mathcal{O}_{2}(z_{1},z_{2})\,|0\rangle$
$\displaystyle=\chi\,\frac{2i}{z_{12}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left[2P\cdot
l_{z_{1}z_{2}}\,e^{-il_{z_{1}z_{2}}\cdot x}+iP\cdot
x\,l_{z_{1}z_{2}}^{2}\int_{0}^{1}dv\,v\,e^{-ivl_{z_{1}z_{2}}\cdot x}\right].$
(28)

Since the operators $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ contain total
derivatives, their matrix elements vanish in the forward limit and need not be
considered in DIS. They provide however corrections of order $m^{2}/Q^{2}$ and
$s/Q^{2}$ in the reaction $\gamma^{*}\gamma\to M\bar{M}$.

##  IV helicity amplitudes in terms of GDAs

In the following we will calculate the helicity amplitudes of
$\gamma^{*}\gamma\to M\bar{M}$, adopting similar techniques to the ones used
for DVCS in Ref. Braun:2012bg . There are three independent helicity
amplitudes, which can be expressed in terms of DDs,

$\displaystyle A_{0+}$
$\displaystyle=2\chi\,\frac{\Delta\cdot\epsilon_{-}}{Q}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta\,\frac{\ln(F)}{F-1},$
$\displaystyle A_{-+}$
$\displaystyle=-\chi\,\frac{(\Delta\cdot\epsilon_{-})^{2}}{2n\cdot\tilde{n}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta^{2}\partial_{F}\\!\left[\frac{1-2F}{F-1}\ln(F)\right],$
$\displaystyle A_{++}$ $\displaystyle=\chi\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left\\{2\ln(F)+\left[\frac{s}{n\cdot\tilde{n}}\,(F-\alpha)+\frac{\beta^{2}\Delta_{T}^{2}}{4n\cdot\tilde{n}}\,\partial_{F}\right]\frac{1}{F-1}\left[\frac{\ln(F)}{2}-\text{Li}_{2}(1)+\text{Li}_{2}(F)\right]\right\\},$
(29)

where $\Delta_{T}^{2}=g_{\perp}^{\mu\nu}\Delta_{\mu}\Delta_{\nu}$ and
$i\epsilon$ is omitted in the functions of $\ln$ and $\text{Li}_{2}$ since it
will not contribute to the amplitudes. Details of the calculations can be
found in Appendix B. $A_{0+}$ and $A_{-+}$ are proportional to
$\Delta\cdot\epsilon_{-}$ and $(\Delta\cdot\epsilon_{-})^{2}$ as indicated by
Eq. (21), respectively, and the amplitudes do not depend on $z_{1}$ and
$z_{2}$ which indicates that the translation invariance is recovered in the
physical amplitudes. The function $F(\alpha,\beta)$ is defined as

$\displaystyle F(\alpha,\beta)=\frac{\alpha-\beta\zeta_{0}+1}{2},$ (30)

where $F=0$ and $F=1$ correspond to the quark momentum fractions $z=1$ and
$z=0$ of the GDAs, respectively.

We notice that there are three types of integrals expressed by DDs in the
obtained amplitudes, namely

$\displaystyle I_{1}$ $\displaystyle=\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,Y(F),$ $\displaystyle I_{2}$
$\displaystyle=\int d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta Y(F),$
$\displaystyle I_{3}$ $\displaystyle=\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta^{2}\partial_{F}Y(F),$ (31)

where $Y(F)$ is some function of $F$. Inserting the identity $\int
dy\,\delta(\beta\zeta_{0}-y-\alpha)=1$ into the integrals above, one can
reexpress the integrals in terms of GDAs by using Eq. (16),

$\displaystyle I_{1}$
$\displaystyle=-\frac{1}{2}\int_{0}^{1}dz\,\Phi(z,\zeta_{0},s)\,\partial_{z}Y(1-z),$
$\displaystyle I_{2}$
$\displaystyle=-\frac{\partial}{\partial\zeta_{0}}\int_{0}^{1}dz\,\Phi(z,\zeta_{0},s)\,Y(1-z),$
$\displaystyle I_{3}$
$\displaystyle=2\,\frac{\partial^{2}}{\partial\zeta_{0}^{2}}\int_{0}^{1}dz\,\Phi(z,\zeta_{0},s)\,Y(1-z),$
(32)

where $y=2z-1$ and $\zeta_{0}=-\beta_{0}\cos{\theta}$ as defined in Eq. (3).
Therefore, we can write the helicity amplitudes as

$\displaystyle A^{(0)}$
$\displaystyle=\chi\left\\{\left(1-\frac{s}{2Q^{2}}\right)\int_{0}^{1}dz\,\frac{\Phi(z,\eta,s)}{1-z}-\frac{s}{Q^{2}}\int_{0}^{1}dz\,\frac{\Phi(z,\eta,s)}{z}\,\ln(1-z)\right.$
$\displaystyle\qquad-\left.\left(\frac{2s}{Q^{2}}\,\eta+\frac{\Delta_{T}^{2}}{\beta_{0}^{2}Q^{2}}\frac{\partial}{\partial\eta}\right)\frac{\partial}{\partial\eta}\int_{0}^{1}dz\,\frac{\Phi(z,\eta,s)}{z}\left[\frac{\ln(1-z)}{2}+\text{Li}_{2}(1-z)-\text{Li}_{2}(1)\right]\right\\},$
$\displaystyle A^{(1)}$
$\displaystyle=\frac{2\chi}{\beta_{0}Q}\frac{\partial}{\partial\eta}\int_{0}^{1}dz\,\Phi(z,\eta,s)\,\frac{\ln(1-z)}{z},$
$\displaystyle A^{(2)}$
$\displaystyle=-\frac{2\chi}{\beta_{0}^{2}Q^{2}}\frac{\partial^{2}}{\partial\eta^{2}}\int_{0}^{1}dz\,\Phi(z,\eta,s)\,\frac{2z-1}{z}\,\ln(1-z),$
(33)

where $\eta=\cos\theta$ and $\Phi=\Phi_{u}+\Phi_{d}$. The GDA for $s$ quarks
is also required in some reactions such as $\gamma^{\ast}\gamma\rightarrow
K^{0}\bar{K}^{0}$ with a charge conjugation-even $K$ meson pair, and we just
need to replace $\chi\Phi$ with
$e_{u}^{2}\Phi_{u}+e_{d}^{2}\Phi_{d}+e_{s}^{2}\Phi_{s}$ in the above
amplitudes. One can clearly see the $\mathcal{O}(s/Q^{2})$ corrections in the
amplitudes, and the target mass correction of order $\mathcal{O}(m^{2}/Q^{2})$
is implicit since it appears in the term $\Delta_{T}^{2}/Q^{2}$ by considering
$\Delta_{T}^{2}=4m^{2}-(1-\zeta_{0}^{2})s$. In general, charge conjugation-
even GDAs can be expanded as Diehl:2000uv

$\displaystyle\Phi(z,\cos\theta,s)=6\,z(1-z)\sum_{\begin{subarray}{c}n=1\\\
n\,\text{odd}\end{subarray}}^{\infty}\sum_{\begin{subarray}{c}l=0\\\
l\,\text{even}\end{subarray}}^{n+1}\tilde{B}_{nl}(s)\,C_{n}^{(3/2)}(2z-1)\,P_{l}(\cos\theta),$
(34)

where $C_{n}^{(3/2)}(x)$ are Gegenbauer polynomials and $P_{l}(x)$ are
Legendre polynomials. Due to this general expression for GDAs, the
singularities of $\frac{1}{z}$, $\frac{1}{1-z}$ and $\ln{(1-z)}$ in the
helicity amplitudes will be compensated by the GDA when $z\rightarrow 0$ and
$z\rightarrow 1$. As a consequence, the amplitudes have no end-point
singularities. In the asymptotic limit ($Q^{2}\rightarrow\infty$), only the
terms with $n=1$ survive,

$\displaystyle\Phi(z,\cos\theta,s)=18\,z(1-z)(2z-1)\left[\tilde{B}_{10}(s)+\tilde{B}_{12}(s)P_{2}(\cos\theta)\right],$
(35)

where the first and second terms correspond to the S-wave and D-wave
production of a meson pair, respectively. The nonvanishing helicity-flip
amplitudes $A_{-+}(A^{(2)})$ and $A_{0+}(A^{(1)})$ indicate the existence of a
D-wave GDA.

##  V Numerical estimates of the higher-twist kinematical contributions

###  V.1 $\pi\pi$ GDA extracted from Belle measurements

Figure 2:  Differential cross section for $e\gamma\to e\pi^{0}\pi^{0}$
calculated with $\pi\pi$ GDA extracted from Belle measurements through a
leading-twist analysis Kumano:2017lhr . Dashed curves show the twist-2
results, while solid curves include the kinematical higher-twist
contributions. The selected values are $s_{e\gamma}=30$ GeV2,
$Q^{2}=9~{}(16,25)$ GeV2 and $\cos\theta=0.2~{}(0.4,0.6,0.8)$ as indicated on
the different panels. Figure 3:  Ratio $d\sigma(2+3+4)/d\sigma(2)$ with the
$\pi\pi$ GDA extracted from Belle measurements, same conventions as in Fig. 2.

The process $\gamma^{*}\gamma\to M\bar{M}$ can be measured in $e^{+}e^{-}$
collisions, which are accessible at KEKB and SuperKEKB. In Ref. Diehl:2000uv ,
the differential cross section for $e\gamma\to eM\bar{M}$ is expressed as

$\displaystyle\frac{d\sigma}{dQ^{2}\,ds\,d(\cos\theta)\,d\varphi}$
$\displaystyle=\frac{\alpha_{\text{em}}^{3}\beta_{0}}{16\pi
s_{e\gamma}^{2}}\,\frac{1}{Q^{2}(1-\epsilon)}\,\Big{[}|A_{++}|^{2}+|A_{-+}|^{2}+2\epsilon\,|A_{0+}|^{2}-2\sqrt{\epsilon(1+\epsilon)}\cos\varphi\,\text{Re}(A_{++}^{\ast}A_{0+}-A_{-+}^{\ast}A_{0+})$
$\displaystyle\quad-2\epsilon\cos(2\varphi)\,\text{Re}(A_{++}^{\ast}A_{-+})\Big{]},$
(36)

where $\varphi$ is the azimuthal angle of the meson pair as illustrated in
Fig. 1 and $s_{e\gamma}$ is the center-of-mass squared energy of $e\gamma$.
$\epsilon$ is defined as usual by

$\displaystyle\epsilon=\frac{1-y}{1-y+y^{2}/2}\quad\text{with}\quad
y=\frac{Q^{2}+s}{s_{e\gamma}}.$ (37)

In the reaction $e\gamma\rightarrow e\pi\pi$, there are two types of
contributions to the cross section. The final $\pi^{+}\pi^{-}$ with negative
charge conjugation couples to a virtual photon, and its contribution is
expressed in terms of the pion electromagnetic form factor. When the charge
conjugation of $\pi\pi$ is positive, the pion pair can be $\pi^{0}\pi^{0}$ or
$\pi^{+}\pi^{-}$. Using factorization, this type of contribution is determined
by GDAs, which we are interested in. The $\pi^{+}\pi^{-}$ GDA is equal to the
one of $\pi^{0}\pi^{0}$ due to the isospin symmetry. However, since
$\pi^{0}\pi^{0}$ are identical bosons, $\cos\theta$ will be restricted to
$0\leq\cos\theta\leq 1$ in Eq. (36). After integration over $\theta$ and
$\varphi$, the cross section for $C$-even $\pi^{+}\pi^{-}$ production is then
twice as large as the one for $\pi^{0}\pi^{0}$ .

In 2016, the Belle Collaboration released the measurements of differential
cross section for $\gamma^{*}+\gamma\rightarrow\pi^{0}+\pi^{0}$ Belle:2015oin
. Since the final state is $\pi^{0}\pi^{0}$, there is no contribution from the
pion electromagnetic form factor. The twist-2 $\pi\pi$ GDA was extracted by
using the leading-twist amplitude Kumano:2017lhr . We use this pion GDA to
estimate the cross section for $e\gamma\to e\pi\pi$ where the integral over
$\varphi$ is performed in Eq. (36),

$\displaystyle\frac{d\sigma}{dQ^{2}\,ds\,d(\cos\theta)}=\frac{\alpha_{\text{em}}^{3}\beta_{0}}{8s_{e\gamma}^{2}}\frac{1}{Q^{2}(1-\epsilon)}\left[|A_{++}|^{2}+|A_{-+}|^{2}+2\epsilon\,|A_{0+}|^{2}\right].$
(38)

In order to show the size of the higher-twist kinematical contributions, Eqs.
(33) and (38) are used to calculate the cross section, and the results are
depicted as the solid lines in Fig. 2. The dashed lines represent the leading-
twist cross sections. Considering the kinematics of Belle measurements, we
choose the values $Q^{2}=9,16,25$ GeV2, $s\in(0.25,4)$ GeV2, and we set
$s_{e\gamma}=30$ GeV2 which is the typical value at Belle. In Fig. 2, black
lines denote $\cos\theta=0.2$ and orange lines correspond to $\cos\theta=0.4$,
while $\cos\theta=0.6$ and $\cos\theta=0.8$ are depicted as red and blue,
respectively. As $Q^{2}$ increases, kinematical contributions become less
important, which is consistent with the fact that the kinematical
contributions are suppressed by $1/Q$ or $1/Q^{2}$. The kinematical
contributions cannot be neglected in the region where $\sqrt{s}\geq 1$ GeV.
The helicity-flip amplitudes $A_{-+}$ and $A_{0+}$ receive only contributions
from the D-wave GDA, and a large difference between two types of cross
sections is displayed around the D-wave resonance region of $f_{2}(1270)$ in
Fig. 2. Hence, the study of the amplitudes $A_{-+}$ and $A_{0+}$ will be
important for the investigation of this resonance region. The kinematical
higher-twist corrections contribute $\sim 15\%$ to the cross section on
average, if one restricts the process $e\gamma\to e\pi\pi$ to the kinematics
of Belle measurements.

In Fig. 3, we also present the ratio $d\sigma(2+3+4)/d\sigma(2)$ where
$d\sigma(i)$ ($i=2,3,4$) is the twist-$i$ contribution to the cross section,
and the colors of the lines indicate different values of $\cos\theta$ as in
Fig. 2. In this figure, the contributions of the kinematical higher-twist
corrections are quite clear, so we can infer that the kinematical corrections
cannot be neglected when $\sqrt{s}>1$ GeV. Around $\sqrt{s}\sim 1.5$ GeV, the
kinematical corrections are dominant in the cross section with
$\cos\theta=0.8$; this appears because the twist-2 cross section is quite tiny
when calculated with the GDA extracted from Belle measurements; this GDA may
however not be accurate in this region since the uncertainties of Belle
measurements are quite large there; this ratio may thus not reflect the real
physics around $\sqrt{s}\sim 1.5$ GeV.

As we have seen, the kinematical corrections are not negligible in the region
$\sqrt{s}>1$ GeV, which turns out to be important for the studies of the pion
EMT form factors. Indeed, since pion GPDs cannot easily be measured in
experiments, GDAs offer a way to investigate the timelike EMT form factors of
pions. The spacelike EMT form factors can then be obtained from the timelike
ones by using dispersion relations, in which case the timelike EMT form
factors of $\sqrt{s}>1$ GeV are needed to be included numerically. As a
consequence, it is important to use the most accurate description of the cross
section with the inclusion of kinematical contributions.

As pointed out above, the uncertainties of Belle measurements Belle:2015oin
are quite large, and the statistical errors are dominant. However, this
situation will be improved substantially soon, since the Belle II
collaboration just started taking data at the SuperKEKB with a much higher
luminosity. Precise measurements of $\gamma^{*}+\gamma\rightarrow M+\bar{M}$
are expected in the near future, and an accurate description of the amplitudes
for the study of GDAs requires the inclusion of kinematical contributions up
to twist 4.

###  V.2 Asymptotic pion GDA

Figure 4:  Differential cross section for $e\gamma\to e\pi^{0}\pi^{0}$ with
the asymptotic $\pi\pi$ GDA described in the text, same conventions as in Fig.
2. Figure 5:  Ratio $d\sigma(2+3+4)/d\sigma(2)$ with the asymptotic
$\pi^{0}\pi^{0}$ GDA described in the text, same conventions as in Fig. 2.

The asymptotic pion GDA used in our calculation is taken from Eq. (68) of Ref.
Diehl:1998dk ,

$\displaystyle\Phi(z,\cos\theta,s)=20\,z(1-z)(2z-1)R_{\pi}\left[\frac{-3+\beta_{0}^{2}}{2}\,e^{i\delta_{0}}+\beta_{0}^{2}e^{i\delta_{2}}P_{l}(\cos\theta)\right],$
(39)

where $\delta_{0}$ and $\delta_{2}$ are $\pi\pi$ elastic scattering phase
shifts in the isospin 0 channel Bydzovsky:2016vdx ; Bydzovsky:2014cda ;
Surovtsev:2010cjf . $R_{\pi}=0.5$ represents the momentum fraction carried by
quarks in the pion meson. In this asymptotic GDA, we do not include the
contribution of the $f_{2}$ resonance. However, we believe it is reasonable to
use this GDA here, since our purpose is not to predict the cross section for
$e\gamma\to e\pi^{0}\pi^{0}$ precisely, but to estimate the magnitude of
kinematical higher-twist contributions and determine whether one can neglect
them or not in the cross section.

In Fig. 4, we show the cross section for $e\gamma\to e\pi^{0}\pi^{0}$ with
fixed $Q^{2}$ and $\cos\theta$, the dashed lines are the twist-2 cross
sections, while the solid ones indicate the cross sections with kinematical
contributions included. The colors of the lines denote different values of
$\cos\theta$ as indicated on the different panels of the figure. Similarly to
the case of the extracted $\pi\pi$ GDA, the kinematical corrections are
important to describe the cross section in the region of $\sqrt{s}>1$ GeV. As
$Q^{2}$ increases, the kinematical contributions become less important.
Compared with Fig. 2, the magnitude of the cross sections are similar at
different $Q^{2}$, even though the asymptotic GDA is very different from the
extracted GDA from Belle measurements. We also present the ratio of
$d\sigma(2+3+4)/d\sigma(2)$ in Fig. 5, where the colors of the lines indicate
different values of $\cos\theta$ as in Fig. 4. In this figure, we can see that
the kinematical contributions can account up to about $40\%$ of the cross
section, which is of course not negligible. Compared with Fig. 3, the
magnitude of the ratios are slightly smaller than the ones obtained from
extracted GDA from Belle measurements.

###  V.3 Model for $\eta\eta$ GDA

Figure 6:  Differential cross section for $e\gamma\to e\eta\eta$ with the
model $\eta\eta$ GDA described in the text, same conventions as in Fig. 2.

As mentioned in the introduction, the kinematical contributions are expected
to be proportional to $m^{2}/Q^{2}$ and $s/Q^{2}$. Since the mass of $\pi$ is
quite small compared to the $Q^{2}$ values of Belle measurements, it is
interesting to check the kinematical contributions for the production of a
pair of slightly heavier mesons, such as $K$ or $\eta$ mesons. On the one
hand, this helps understand how the target mass corrections of order
$\mathcal{O}(m^{2}/Q^{2})$ affect the cross section. On the other hand, the
$K\bar{K}$ and $\eta\eta$ GDAs can also be measured by Belle and Belle II
experiments; the Belle Collaboration indeed released the cross section for
$\gamma^{*}+\gamma\rightarrow K_{0}+\bar{K}_{0}$ in 2018 Belle:2017xsz . It
will therefore be necessary to investigate the kinematical corrections for the
production of $K$ and $\eta$ meson pairs. Unfortunately, there is almost no
information on these GDAs at present. Here, we just simply replace the mass of
$\pi$ with the one of $\eta$ in Eq. (39) and keep other parameters unchanged,
then use this GDA to estimate the cross section and the ratio of various twist
contributions for $e\gamma\to e\eta\eta$ with $s_{e\gamma}=30$ GeV2. The
differential cross section for $e\gamma\to e\eta\eta$ is shown in Fig. 6, and
the ratio of $d\sigma(2+3+4)/d\sigma(2)$ is presented in Fig. 7. The values of
$Q^{2}$ are chosen as $Q^{2}=$ 9, 16 and 25 GeV2 together with 1.2 GeV
$\leq\sqrt{s}\leq$ 2.2 GeV, and the black (orange, red, blue) lines in Fig. 6
and Fig. 7. represent $\cos\theta=0.2$ (0.4, 0.6, 0.8). The kinematical
contributions account for somewhat less than $40\%$ of the cross section,
which cannot be neglected either. The kinematical higher-twist contributions
have a significant impact on the cross section even in the region
($\sqrt{s}\sim 1.2$ GeV) which is close to the $\eta\eta$ threshold. Compared
with the ratios of Figs. (3) and (5), the kinematical higher-twist
contributions are always negative, and those negative contributions can only
come from the amplitude of $A_{++}$, since $A_{+-}$ and $A_{0-}$ always
contribute to the cross section positively as indicated by Eqs. (33) and (38).
The importance of the kinematical contributions does diminish as $m$ increases
from the pion mass to the $\eta$ mass, simply because the negative kinematical
higher-twist contributions from $A_{++}$ are compensated by the positive ones
from $A_{0-}$ and $A_{+-}$.

Figure 7: Ratio $d\sigma(2+3+4)/d\sigma(2)$ with the model $\eta\eta$ GDA
described in the text, same conventions as in Fig. 2.

##  VI Summary

In this paper we presented a complete calculation of kinematical higher-twist
corrections for the helicity amplitudes of the reaction
$\gamma^{*}+\gamma\rightarrow M+\bar{M}$ up to twist 4, where only the
leading-twist GDA is involved in the description of the cross section. In case
of $\pi\pi$ production, we use two types of GDAs to estimate the kinematical
higer-twist contributions in the cross section, namely the $\pi\pi$ GDA
extracted from the Belle measurements and the asymptotic $\pi\pi$ GDA. Even
though those two GDAs are very different, both of them lead to kinematical
corrections which cannot be neglected for $\gamma^{*}+\gamma\rightarrow\pi\pi$
considering the kinematics accessible at Belle and Belle II. Moreover, the
relative magnitude of the higher-twist corrections is also comparable for the
two types of $\pi\pi$ GDAs at different values of $Q^{2}$, as seen from Figs.
(3) and (5).

Due to the small pion mass, only kinematical corrections of the type
$\mathcal{O}(s/Q^{2})$ contribute to the cross section for
$\gamma^{\ast}+\gamma\rightarrow\pi\pi$. The production of a pair of slightly
heavier mesons is needed to check how the target mass corrections
$m^{2}/Q^{2}$ affect the cross section. Since the $\eta\eta$ GDA is an unknown
quantity at present, we calculate this effect with a model GDA identical to
the asymptotic $\pi\pi$ GDA except that the mass of $\eta$ is used. This
calculation indicates that kinematical corrections are also not negligible in
this case. Furthermore, the negative kinematical corrections from the
amplitude $A_{++}$ are dominant over the higher-twist contributions from
$A_{+-}$ and $A_{0-}$, and the kinematical corrections are always negative in
the cross section, which is different from the $\pi\pi$ case where kinematical
corrections can go both ways.

In conclusion, let us stress that while the uncertainties of present Belle
measurements Belle:2015oin  are too large for our study to invalidate the
conclusions of Kumano:2017lhr , the situation should change substantially in a
near future since Belle II collaboration just started taking data at the
SuperKEKB with a much higher luminosity. An accurate description of the
amplitudes in terms of GDAs will therefore require the inclusion of
kinematical contributions up to twist 4. Note also that precise measurements
of $\gamma^{*}+\gamma\rightarrow M+\bar{M}$ for various mesons will be of
utmost importance to address the questions of the pion EMT form factors and of
the impact-parameter representation of GDAs Pire:2002ut .

In the future, this work can be extended to the production of other meson
pairs, for example the $\gamma^{\ast}\gamma\rightarrow\pi\eta$ channel which
should help unraveling the quark and gluon structure of hybrid meson
($J^{PC}=1^{-+}$) Anikin:2006du . The scattering amplitude is likely to be
sensitive to sizeable kinematical higher twist contributions. This channel has
recently been advocated Teryaev:2022pke  to be related to shear viscosity of
quarks in hadronic matter. The production of a pair of vector mesons should
also be discussed, opening the way to a meaningful extraction of the EMT form
factors for $\rho$ or $\omega$ mesons.

Similar relations apply to the timelike process amplitude $\gamma^{*}\to
M\bar{M}\gamma$ which opens another access to GDAs Lu:2006ut  through the
interference with the initial-state radiation amplitude in the process
$e^{+}e^{-}\to M\bar{M}\gamma$, as experimentally proven by the BABAR
collaboration BaBar:2015onb . The extension of our work to the process
$\gamma^{*}+\gamma\rightarrow N+\bar{N}$ will also be needed if Belle II
detector is able to detect this channel.

##  VII Acknowledgements

We acknowledge useful discussions with Lech Szymanowski, Pieter Taels, Oleg
Teryaev, Sadaharu Uehara and Jakub Wagner. Qin-Tao Song is supported by the
National Natural Science Foundation of China under Grant Number 12005191 and
the China Scholarship Council for visiting Ecole Polytechnique.

##  Appendix A Expressions for the T-product of two electromagnetic currents

Here we give the detailed expressions of $\mathbb{V}_{\mu}$,
$\mathbb{W}_{\mu}$, $\mathbb{X}$ and $\mathbb{Y}$ Braun:2011dg ; Braun:2011th
, which are used in the T-product of two electromagnetic currents in Eq. (17).

$\displaystyle\mathbb{V}_{\mu}(z_{1},z_{2})$
$\displaystyle=\mathfrak{B}_{\mu}(z_{1},z_{2})-\mathfrak{B}_{\mu}(z_{2},z_{1})+x_{\mu}\Delta\mathbb{A}(z_{1},z_{2}),$
$\displaystyle\mathbb{W}_{\mu}(z_{1},z_{2})$
$\displaystyle=-\mathfrak{B}_{\mu}(z_{1},z_{2})-\mathfrak{B}_{\mu}(z_{2},z_{1}),$
$\displaystyle\mathbb{X}(z_{1},z_{2})$
$\displaystyle=\mathbb{C}(z_{1},z_{2})-\mathbb{C}(z_{2},z_{1}),$
$\displaystyle\mathbb{Y}(z_{1},z_{2})$
$\displaystyle=-\mathbb{C}(z_{1},z_{2})-\mathbb{C}(z_{2},z_{1}),$ (40)

where
$\Delta\mathbb{A}(z_{1},z_{2})=\mathbb{A}(z_{1},z_{2})-\mathbb{A}(z_{2},z_{1})$
is a pure twist-4 operator,

$\displaystyle\mathbb{A}(z_{1},z_{2})$
$\displaystyle=\frac{1}{4}\int_{0}^{1}du\left\\{u^{2}\ln{u}\,z_{1}z_{2}\,\mathcal{O}_{1}(z_{1}u,z_{2}u)+\left[\left(z_{2}\partial_{z_{2}}-\frac{z_{1}}{z_{12}}-\ln{u}\,z_{2}\partial_{z_{2}}^{2}z_{12}\right)\mathcal{R}(uz_{1},uz_{2})-(1\leftrightarrow
2)\right]\right\\}.$ (41)

The function $\mathcal{R}(z_{1},z_{2})$ is related to the total derivative
operators $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ through

$\displaystyle\mathcal{R}(z_{1},z_{2})=z_{12}\int_{z_{2}}^{z_{1}}\frac{dw}{z_{12}}\int_{z_{2}}^{w}\frac{dw_{1}}{z_{12}}\,\frac{w_{1}-z_{2}}{z_{1}-w_{1}}\left[\frac{1}{2}\,S_{+}\mathcal{O}_{1}(w,w_{1})-(S_{0}-1)\mathcal{O}_{2}(w,w_{1})\right],$
(42)

where $S_{+}$ and $S_{0}$ are differential operators of $w$ and $w_{1}$,

$\displaystyle S_{+}$
$\displaystyle=w^{2}\partial_{w}+2w+w_{1}^{2}\partial_{w_{1}}+2w_{1},$
$\displaystyle S_{0}$ $\displaystyle=w\partial_{w}+w_{1}\partial_{w_{1}}+2.$
(43)

The operator $\mathfrak{B}_{\mu}(z_{1},z_{2})$ contains all twists starting
from twist 2,

$\displaystyle\mathfrak{B}_{\mu}(z_{1},z_{2})=\mathfrak{B}_{\mu}^{t=2}(z_{1},z_{2})+\mathfrak{B}_{\mu}^{t=3}(z_{1},z_{2})+\mathfrak{B}_{\mu}^{t=4}(z_{1},z_{2}).$
(44)

The twist-2 and twist-3 parts are defined as

$\displaystyle\mathfrak{B}_{\mu}^{t=2}(z_{1},z_{2})=$
$\displaystyle\frac{1}{2}\,\partial_{\mu}\int_{0}^{1}du\,\mathcal{O}_{++}^{t=2}(uz_{1},uz_{2}),$
$\displaystyle\mathfrak{B}_{\mu}^{t=3}(z_{1},z_{2})=$
$\displaystyle\frac{1}{4}\int_{0}^{1}du\,u\int_{z_{2}}^{z_{1}}\frac{dw}{z_{12}}\left\\{\left[i\mathbf{P}^{\nu},K_{\alpha\nu\beta\mu}(x^{\alpha}\partial^{\beta})z_{1}\mathcal{O}_{++}^{t=2}(z_{1}u,wu)+K_{\mu\alpha\nu\beta}(x^{\alpha}\partial^{\beta})z_{2}\mathcal{O}_{++}^{t=2}(wu,z_{2}u)\right]\right.$
$\displaystyle\left.+\ln(u)\,\partial_{\mu}x^{2}\partial_{\nu}\left[i\mathbf{P}^{\nu},z_{1}\mathcal{O}_{++}^{t=2}(z_{1}u,\,wu)+z_{2}\mathcal{O}_{++}^{t=2}(wu,z_{2}u)\right]\right\\},$
(45)

where

$\displaystyle
K^{\mu\alpha\nu\beta}=(g^{\mu\alpha}g^{\nu\beta}-g^{\mu\nu}g^{\alpha\beta}+g^{\mu\beta}g^{\nu\alpha})-i\epsilon^{\mu\alpha\nu\beta}.$
(46)

The twist-4 part $\mathfrak{B}_{\mu}^{t=4}(z_{1},z_{2})$ is more complex than
the twist-2 and twist-3 ones,

$\displaystyle\mathfrak{B}_{\mu}^{t=4}(z_{1},z_{2})$
$\displaystyle=\frac{x^{2}}{8}\,\partial_{\mu}\int_{0}^{1}\frac{du}{u^{2}}\bigg{\\{}u^{2}(1-u^{2}+u^{2}\ln{u})\,z_{1}z_{2}\mathcal{O}_{1}(z_{1}u,\,z_{2}u)$
$\displaystyle\quad-\left[\left((1-u^{2})\left(z_{2}\partial_{z_{2}}-\frac{z_{1}}{z_{12}}\right)+(1-u^{2}+u^{2}\ln{u})\,z_{2}\partial_{z_{2}}^{2}z_{12}\right)\mathcal{R}(uz_{1},uz_{2})-(1\leftrightarrow
2)\right]\bigg{\\}}.$ (47)

$\mathbb{C}(z_{1},z_{2})$ is pure twist 4 and can be expressed as

$\displaystyle\mathbb{C}(z_{1},z_{2})=$
$\displaystyle-\frac{1}{4}\int_{0}^{1}\frac{du}{u^{2}}\,\mathcal{R}(uz_{1},uz_{2}).$
(48)

##  Appendix B Calculation techniques for helicity amplitudes

There are two helicity-flip amplitudes $A_{0+}$ and $A_{-+}$ in Eq. (21).
Angular momentum conservation implies that they are proportional to a given
power of the transverse momentum transfer as
$A_{0+}\propto\Delta\cdot\epsilon_{-}$ and
$A_{-+}\propto(\Delta\cdot\epsilon_{-})^{2}$. Therefore, the twist-4 part of
Eq. (17) will be beyond the accuracy of this work due to the additional factor
of $\Delta\cdot\epsilon_{-}$,

$\displaystyle A_{0+}=A_{0+}^{t=2}+A_{0+}^{t=3},\qquad
A_{-+}=A_{-+}^{t=2}+A_{-+}^{t=3},$ (49)

We consider the matrix element of the twist-2 part of
$T_{0+}=T_{\mu\nu}\epsilon_{0}^{\mu}\tilde{\epsilon}_{+}^{\nu}$ and substitute
it into Eq. (5),

$\displaystyle A_{0+}^{t=2}$ $\displaystyle=\int d^{4}x\,e^{-ir\cdot
x}\,\langle\bar{M}(p_{2})M(p_{1})|\,T_{0+}^{t=2}\,|0\rangle$
$\displaystyle=-\chi\,\frac{\Delta\cdot\epsilon_{-}}{Q}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta\int_{0}^{1}du\,\frac{4n\cdot\tilde{n}}{(r+ul_{z_{1}z_{2}})^{2}+i\epsilon}$
$\displaystyle=2\chi\,\frac{\Delta\cdot\epsilon_{-}}{Q}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta\,\frac{\ln(F-i\epsilon)-\ln(z_{1})}{F-z_{1}},$
(50)

where $(r+ul_{z_{1}z_{2}})^{2}=-2\,n\cdot\tilde{n}\,(uF+(1-u)z_{1})$ is
obtained by neglecting the terms of order $\mathcal{O}(s,m^{2})$ since they
will not contribute at the $1/Q^{2}$ accuracy. Let us mention that $i\epsilon$
can be omitted in $\ln(F-i\epsilon)$ since $F$ is always positive ($0\leq
F\leq 1$) and there is no branch cut. Indeed, $\ln(F)$ is divergent when $F=0$
$(z=1)$, but the GDA vanishes at $z=1$. Similarly, the twist-3 part is
obtained,

$\displaystyle A_{0+}^{t=3}$ $\displaystyle=\int d^{4}x\,e^{-ir\cdot
x}\,\langle\bar{M}(p_{2})M(p_{1})|\,T_{0+}^{t=3}\,|0\rangle$
$\displaystyle=2\chi\,\frac{\Delta\cdot\epsilon_{-}}{Q}\,z_{2}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta\int_{0}^{1}du\,u\int_{z_{2}}^{z_{1}}dw\,\frac{(2n\cdot\tilde{n})^{2}}{\left[(r+ul_{wz_{2}})^{2}+i\epsilon\right]^{2}}$
$\displaystyle=2\chi\,\frac{\Delta\cdot\epsilon_{-}}{Q}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta\,\frac{-z_{2}}{F-z_{1}}\left[\frac{\ln(F-i\epsilon)}{F-1}-\frac{\ln(z_{1})}{z_{2}}\right]$
(51)

with

$\displaystyle(r+ul_{wz_{2}})^{2}=-2n\cdot\tilde{n}\left[z_{1}-uw+(w-z_{2})F\right].$
(52)

We note that although the twist-2 and twist-3 amplitudes depend on $z_{1}$ and
$z_{2}$, this dependence disappears in their sum,

$\displaystyle A_{0+}=2\chi\,\frac{\Delta\cdot\epsilon_{-}}{Q}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta\,\frac{\ln(F)}{F-1}\,,$ (53)

which indicates that translation invariance is recovered in the physical
amplitudes. Once again $i\epsilon$ is omitted in $\ln(F-i\epsilon)$ as well.

The calculation of $A_{-+}$ is quite similar to the one of $A_{0+}$,

$\displaystyle A_{-+}^{t=2}$
$\displaystyle=-2\chi\,(\Delta\cdot\epsilon_{-})^{2}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta^{2}\int_{0}^{1}du\,u\,\frac{1}{(r+ul_{z_{1}z_{2}})^{2}+i\epsilon}$
$\displaystyle=\chi\,\frac{(\Delta\cdot\epsilon_{-})^{2}}{n\cdot\tilde{n}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta^{2}\,\partial_{F}\left[\frac{F\ln(F-i\epsilon)}{F-z_{1}}-\frac{z_{1}\ln(z_{1})}{F-z_{1}}\right],$
$\displaystyle A_{-+}^{t=3}$
$\displaystyle=-2\chi\,(\Delta\cdot\epsilon_{-})^{2}\,n\cdot\tilde{n}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta^{2}\int_{0}^{1}du\,u^{2}\int_{z_{2}}^{z_{1}}dw\left[\frac{z_{1}(w-z_{1})}{\left[(r+ul_{z_{1}w})^{2}+i\epsilon\right]^{2}}+\frac{z_{2}(z_{2}-w)}{\left[(r+ul_{wz_{2}})^{2}+i\epsilon\right]^{2}}\right]$
$\displaystyle=-\chi\,\frac{(\Delta\cdot\epsilon_{-})^{2}}{2n\cdot\tilde{n}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta^{2}\,\partial_{F}\left[\frac{-\ln(F-i\epsilon)}{F-1}+\frac{2z_{1}\ln(F-i\epsilon)}{F-z_{1}}-\frac{2z_{1}\ln(z_{1})}{F-z_{1}}\right].$
(54)

For the sum, one gets

$\displaystyle
A_{-+}=-\chi\,\frac{(\Delta\cdot\epsilon_{-})^{2}}{2n\cdot\tilde{n}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta^{2}\,\partial_{F}\left[\frac{1-2F}{F-1}\,\ln(F)\right],$
(55)

which does not depend on $z_{1}$ and $z_{2}$, and which is proportional to
$(\Delta\cdot\epsilon_{-})^{2}$ as shown in Eq. (21).

The calculation of $A_{++}$ is more lengthy than the ones of helicity-flip
amplitudes, since it contains the contributions of twist 2, twist 3 and twist
4.

$\displaystyle A_{++}=A_{++}^{t=2}+A_{++}^{t=3}+A_{++}^{t=4},\qquad
A_{++}=\epsilon_{+}^{\mu}\tilde{\epsilon}_{+}^{\nu}A_{\mu\nu}.$ (56)

Taking the trace of Eq. (6), one obtains Braun:2012bg

$\displaystyle A_{++}=-\frac{1}{2}A_{\,\
\mu}^{\mu}+A_{\mu\nu}\,\frac{(n^{\mu}\tilde{n}^{\nu}-\tilde{n}^{\mu}n^{\nu})}{2n\cdot\tilde{n}}.$
(57)

We simplify the operator expansion of $T_{++}$ in Eq. (17) by using Eq. (57),
and we take the matrix elements of $T_{++}$ to obtain the contribution of each
twist to $A_{++}$ Braun:2012bg ,

$\displaystyle A_{++}^{t=2}$
$\displaystyle=-\int\frac{d^{4}x}{\pi^{2}}\,\frac{e^{-ir\cdot
x}}{x^{4}}\,\langle\bar{M}(p_{2})M(p_{1})|\,\mathcal{O}_{++}^{t=2}(z_{1},\,z_{2})\,|0\rangle,$
$\displaystyle A_{++}^{t=3}$
$\displaystyle=2\int\frac{d^{4}x}{\pi^{2}}\,\frac{e^{-ir\cdot
x}}{x^{4}}\,(x\cdot\epsilon_{+})\,\epsilon_{-}^{\mu}\langle\bar{M}(p_{2})M(p_{1})|\,\mathfrak{B}_{\mu}^{t=3}(z_{1},z_{2})+\mathfrak{B}_{\mu}^{t=3}(z_{2},z_{1})\,|0\rangle,$
$\displaystyle A_{++}^{t=4}$
$\displaystyle=-\frac{1}{4}\int\frac{d^{4}x}{\pi^{2}}\,\frac{e^{-ir\cdot
x}}{x^{2}}\int_{0}^{1}du\,\langle\bar{M}(p_{2})M(p_{1})|\left\\{z_{1}z_{2}u^{2}\mathcal{O}_{1}(z_{1}u,z_{2}u)\right.$
$\displaystyle\quad\left.-[z_{2}\partial_{2}^{2}z_{12}\mathcal{R}(z_{1}u,z_{2}u)-(1\leftrightarrow
2)]-2\,\mathcal{R}(z_{2},z_{1})\right\\}|0\rangle,$ (58)

where the expressions for $\mathfrak{B}_{\mu}^{t=3}(z_{1},z_{2})$ and
$\mathcal{R}(z_{2},z_{1})$ can be found in Appendix A.

The calculation of $A_{++}^{t=2}$ is rather straightforward with the help of
Eq. (26),

$\displaystyle A_{++}^{t=2}=$ $\displaystyle 2\chi\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left[\ln((r+ul_{z_{1}z_{2}})^{2}+i\epsilon)-l_{z_{1}z_{2}}^{2}\int_{0}^{1}du\,u\,\frac{1}{(r+ul_{z_{1}z_{2}})^{2}+i\epsilon}\right]$
$\displaystyle=$ $\displaystyle 2\chi\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left\\{\ln(F-i\epsilon)+\frac{\beta^{2}\Delta_{T}^{2}}{8n\cdot\tilde{n}}\,\partial_{F}\left[z_{1}\,\frac{\ln(F-i\epsilon)-\ln(z_{1})}{F-z_{1}}\right]\right.$
$\displaystyle+\left.\frac{s}{2n\cdot\tilde{n}}\left[z_{1}z_{2}+(F-z_{1})\alpha+(1-F)F\right]\partial_{F}\left[\frac{F\ln(F-i\epsilon)}{F-z_{1}}-\frac{z_{1}\ln(z_{1})}{F-z_{1}}\right]\right\\},$
(59)

where Eq. (15) is used to eliminate the irrelevant terms, and the first term
in Eq. (59) is actually the twist-2 amplitude given by Ref. Diehl:1998dk . The
twist-3 contribution is expressed as

$\displaystyle A_{++}^{t=3}=$ $\displaystyle\frac{\chi}{n\cdot\tilde{n}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\beta\left[\zeta_{0}s-\frac{\Delta_{T}^{2}}{4}\,\beta\,\partial_{F}\right]\frac{\ln(F-i\epsilon)}{F-1},$
(60)

and is independent of $z_{1}$ and $z_{2}$. There are several terms
contributing to the amplitude of $A_{++}^{t=4}$ in Eq. (58), which we denote
as

$\displaystyle A_{++}^{t=4}=A_{(1)}^{t=4}+A_{(2+3)}^{t=4}+A_{(4)}^{t=4}.$ (61)

$A_{(1)}^{t=4}$ is calculated by using Eq. (28),

$\displaystyle A_{(1)}^{t=4}=-\chi\,\frac{z_{1}z_{2}s}{n\cdot\tilde{n}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\,\partial_{F}\left[\frac{F\ln(F-i\epsilon)}{F-z_{1}}-\frac{z_{1}\ln(z_{1})}{F-z_{1}}\right].$
(62)

When calculating the other terms in $A_{++}^{t=4}$, one finds divergences in
$A_{(2)}^{t=4}$ and $A_{(3)}^{t=4}$ associated with the real (onshell) photon.
Here, the photon is set as an offshell one so as to regularize the
divergences,

$\displaystyle q_{2}$ $\displaystyle\rightarrow q_{2}+(\xi-
z_{1})(q_{1}+q_{2}),$ $\displaystyle r=z_{1}q_{1}+z_{2}q_{2}$
$\displaystyle\rightarrow r=-q_{2}+\xi(q_{1}+q_{2}),$ (63)

and we take $\xi=z_{1}$ at the end of the calculations to obtain the final
results. As shown in Eq. (58), the operator $\mathcal{R}(z_{1},z_{2})$ is
involved in the remaining terms in $A^{++}_{t=4}$. First, we simplify the
matrix element of $\mathcal{R}(z_{1},z_{2})$,

$\displaystyle\langle\bar{M}(p_{2})M(p_{1})|\,\mathcal{R}(z_{1},z_{2})\,|0\rangle$
$\displaystyle=-2i\chi\int_{z_{2}}^{z_{1}}dw\,\int_{z_{2}}^{w}\,\frac{dw_{1}}{z_{12}}\frac{w_{1}-z_{2}}{z_{1}-w_{1}}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left[\frac{s}{2}\,S_{+}+2(S_{0}-1)P\cdot
l_{ww_{1}}+iP\cdot x\,l_{ww_{1}}^{2}\right]\frac{e^{-il_{ww_{1}}\cdot
x}}{w-w_{1}},$ (64)

where $z_{12}$ is kept so that $z_{12}=1$ is not a necessary condition. One
can make the replacements $z_{1}\rightarrow uz_{1}$ and $z_{2}\rightarrow
uz_{2}$ to obtain $\mathcal{R}(z_{1}u,z_{2}u)$. We substitute the matrix
elements $\mathcal{R}(z_{1}u,z_{2}u)$ and $\mathcal{R}(z_{2}u,z_{1}u)$ into
Eq. (58) and we take $\xi=z_{1}$,

$\displaystyle A_{(2+3)}^{t=4}$ $\displaystyle=\frac{\chi}{2}\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left\\{\frac{s}{n\cdot\tilde{n}}\left[\alpha+(F-1)F\,\partial_{F}\right]-\frac{\beta^{2}\Delta_{T}^{2}}{4n\cdot\tilde{n}}\,\partial_{F}\right\\}\frac{1}{F-1}$
$\displaystyle\quad\times\left\\{z_{2}\left[\frac{F\ln(F-i\epsilon)}{F-z_{1}}-\frac{z_{1}}{z_{2}}\,\frac{F-1}{F-z_{1}}\,\ln(z_{1})\right]+z_{1}\left[\frac{(F-1)\ln(1-F-i\epsilon)}{F+z_{2}}+\frac{1-F}{F+z_{2}}\,\ln(z_{1})\right]\right\\},$
(65)

Similarly, the last contribution to $A_{++}^{t=4}$ is given by

$\displaystyle A_{(4)}^{t=4}=$ $\displaystyle\chi\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left\\{\frac{s}{n\cdot\tilde{n}}\left[\alpha+(F-1)F\,\partial_{F}\right]-\frac{\beta^{2}\Delta_{T}^{2}}{4n\cdot\tilde{n}}\,\partial_{F}\right\\}\frac{1}{1-F}\left[\ln(F-i\epsilon)-\text{Li}_{2}(1)+\text{Li}_{2}(F+i\epsilon)\right],$
(66)

where the replacements $\alpha\rightarrow-\alpha$, $\beta\rightarrow-\beta$
and $F\rightarrow 1-F$ are used to simplify the amplitude.

Summing over all contributions leads finally to

$\displaystyle A_{++}$
$\displaystyle=A_{++}^{t=2}+A_{++}^{t=3}+A_{(1)}^{t=4}+A_{(2+3)}^{t=4}+A_{(4)}^{t=4}$
$\displaystyle=\chi\int
d\beta\,d\alpha\,\phi(\beta,\alpha)\left\\{2\ln(F)+\left[\frac{s}{n\cdot\tilde{n}}\,(F-\alpha)+\frac{\beta^{2}\Delta_{T}^{2}}{4n\cdot\tilde{n}}\,\partial_{F}\right]\frac{1}{F-1}\left[\frac{\ln(F)}{2}-\text{Li}_{2}(1)+\text{Li}_{2}(F)\right]\right\\},$
(67)

where the $z_{1}$ and $z_{2}$ dependences disappear as expected, and where
$i\epsilon$ is omitted in the functions of $\ln(F-i\epsilon)$ and
$\text{Li}_{2}(F-i\epsilon)$.

## References

  * (1) D. Müller, D. Robaschik, B. Geyer, F. M. Dittes and J. Hořejši, Fortsch. Phys. 42 (1994), 101-141.
  * (2) M. Diehl, T. Gousset, B. Pire and O. Teryaev, Phys. Rev. Lett. 81 (1998), 1782-1785.
  * (3) M. V. Polyakov, Nucl. Phys. B 555 (1999), 231.
  * (4) M. Diehl, Phys. Rept. 388 (2003), 41-277.
  * (5) A. V. Belitsky and A. V. Radyushkin, Phys. Rept. 418 (2005), 1-387.
  * (6) S. Boffi and B. Pasquini, Riv. Nuovo Cim. 30 (2007) no.9, 387-448.
  * (7) K. Goeke, M. V. Polyakov and M. Vanderhaeghen, Prog. Part. Nucl. Phys. 47 (2001), 401-515.
  * (8) M. Masuda et al. [Belle], Phys. Rev. D 93 (2016) no.3, 032003.
  * (9) S. Kumano, Q. T. Song and O. V. Teryaev, Phys. Rev. D 97 (2018) no.1, 014020.
  * (10) C. H. Chen and H. N. Li, Phys. Lett. B 561 (2003), 258-265.
  * (11) W. F. Wang, H. N. Li, W. Wang and C. D. Lü, Phys. Rev. D 91 (2015) no.9, 094024.
  * (12) Y. Li, A. J. Ma, W. F. Wang and Z. J. Xiao, Phys. Rev. D 95 (2017) no.5, 056008.
  * (13) M. K. Jia, C. Q. Zhang, J. M. Li and Z. Rui, Phys. Rev. D 104 (2021) no.7, 073001.
  * (14) M. Burkardt, Phys. Rev. D 62 (2000), 071503 [erratum: Phys. Rev. D 66 (2002), 119903].
  * (15) J. P. Ralston and B. Pire, Phys. Rev. D 66 (2002), 111501.
  * (16) M. Diehl, Eur. Phys. J. C 25 (2002), 223-232 [erratum: Eur. Phys. J. C 31 (2003), 277-278].
  * (17) B. Pire and L. Szymanowski, Phys. Lett. B 556 (2003), 129-134.
  * (18) X. D. Ji, Phys. Rev. Lett. 78 (1997), 610-613.
  * (19) X. D. Ji, Phys. Rev. D 55 (1997), 7114-7125.
  * (20) M. Diehl, T. Gousset and B. Pire, Phys. Rev. D 62 (2000), 073014.
  * (21) M. V. Polyakov, Phys. Lett. B 555 (2003), 57-62.
  * (22) K. Goeke, J. Grabis, J. Ossmann, M. V. Polyakov, P. Schweitzer, A. Silva and D. Urbano, Phys. Rev. D 75 (2007), 094021.
  * (23) M. Mai and P. Schweitzer, Phys. Rev. D 86 (2012), 096002.
  * (24) M. V. Polyakov and P. Schweitzer, Int. J. Mod. Phys. A 33 (2018) no.26, 1830025.
  * (25) C. Lorcé, H. Moutarde and A. P. Trawiński, Eur. Phys. J. C 79 (2019) no.1, 89.
  * (26) V. D. Burkert, L. Elouadrhiri and F. X. Girod, Nature 557 (2018) no.7705, 396-399.
  * (27) K. Kumerički, Nature 570 (2019) no.7759, E1-E2.
  * (28) H. Dutrieux, C. Lorcé, H. Moutarde, P. Sznajder, A. Trawiński and J. Wagner, Eur. Phys. J. C 81 (2021) no.4, 300.
  * (29) D. Amrath, M. Diehl and J. P. Lansberg, Eur. Phys. J. C 58 (2008), 179-192.
  * (30) J. M. M. Chávez, V. Bertone, F. De Soto Borrero, M. Defurne, C. Mezrag, H. Moutarde, J. Rodríguez-Quintero and J. Segovia, Phys. Rev. Lett. 128 (2022) no.20, 202501.
  * (31) J. M. M. Chavez, V. Bertone, F. De Soto Borrero, M. Defurne, C. Mezrag, H. Moutarde, J. Rodríguez-Quintero and J. Segovia, Phys. Rev. D 105 (2022) no.9, 094012.
  * (32) J. P. Lansberg, B. Pire and L. Szymanowski, Phys. Rev. D 73 (2006), 074014.
  * (33) I. V. Anikin, B. Pire and O. V. Teryaev, Phys. Lett. B 626 (2005), 86-94.
  * (34) V. M. Braun and A. N. Manashov, Phys. Rev. Lett. 107 (2011), 202001.
  * (35) V. M. Braun and A. N. Manashov, JHEP 01 (2012), 085.
  * (36) V. M. Braun and A. N. Manashov, Prog. Part. Nucl. Phys. 67 (2012), 162-167.
  * (37) V. M. Braun, A. N. Manashov and B. Pirnay, Phys. Rev. D 86 (2012), 014003.
  * (38) O. Nachtmann, Nucl. Phys. B 63 (1973), 237-247.
  * (39) O. V. Teryaev, Phys. Lett. B 510 (2001), 125-132.
  * (40) V. M. Braun, A. N. Manashov and J. Rohrwild, Nucl. Phys. B 807 (2009), 89-137.
  * (41) V. M. Braun, A. N. Manashov and J. Rohrwild, Nucl. Phys. B 826 (2010), 235-293.
  * (42) P. Bydžovský, R. Kamiński and V. Nazari, Phys. Rev. D 94 (2016) 11, 116013.
  * (43) P. Bydžovský, R. Kamiński and V. Nazari, Phys. Rev. D 90 (2014)11, 116005.
  * (44) Y. S. Surovtsev, P. Bydzovsky, R. Kaminski and M. Nagy, Phys. Rev. D 81 (2010), 016001.
  * (45) M. Masuda et al. [Belle], Phys. Rev. D 97 (2018) no.5, 052003.
  * (46) I. V. Anikin, B. Pire, L. Szymanowski, O. V. Teryaev and S. Wallon, Eur. Phys. J. C 47 (2006), 71-79.
  * (47) O. Teryaev, [arXiv:2204.09742 [hep-ph]].
  * (48) Z. Lu and I. Schmidt, Phys. Rev. D 73 (2006), 094021 [erratum: Phys. Rev. D 75 (2007), 099902].
  * (49) J. P. Lees et al. [BaBar], Phys. Rev. D 92 (2015) no.7, 072015.


# Coupled Cavity Model: Correctness and Limitations

M.I.Ayzatsky
National Science Center
Kharkov Institute of Physics and Technology (NSC KIPT),
610108, Kharkov, Ukraine<EMAIL_ADDRESS>
###### Abstract

Results of analysis of correctness and limitations of the classical Coupled
Cavity Model are presented in the paper. It is shown that in the case of an
infinite chain of resonators, there are spurious solutions of the
characteristic equation. These spurious solutions do not violate the
correctness of direct numerical calculations, but their existence makes it
difficult (or even impossible) to use approximate WKB methods for analysing
chains with slowly varying parameters.

##  1 Introduction

The classical Coupled Cavity Model (CCM) is widely used for design of
RF/microwave devices, in particular,the accelerating structures (see, for
example, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 8, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20] ), the travelling wave tubes (see, for example,[21]), the narrow-band
bandpass filters (see, for example, [22, 23, 24]) and so on. It is also found
wide application in optic and metamaterial structures (see, for example, [25,
26]).

These devises are based on chains of coupled resonators. Such chains belong to
the class of closed structured waveguides - waveguides that consist of
similar, but not always identical, cells that couples through openings in
dividing walls. Each cell couples only with two neighbouring cells. In opening
structured waveguides (optic and metamaterial structures [25, 26]) each cell
couples with all cells of the chain.

Despite the differences, the general approach was the same - it is necessary
to find eigen modes of noncoupling cells and construct the coupling
coefficients of intercoupled resonators and the external quality factors of
the input and output resonators. The general coupling matrix is of importance
for representing a wide range of coupled resonators topologies.

On the base of simple model of infinite chain of cylindrical resonators that
couple through circular openings in thing dividing walls we show that the
classical CCM has problems in describing field distributions in closed
structured waveguides.

##  2 Finite Chain of Resonators. Main equations

Consider a chain of $N_{R}$ cylindrical resonators with annular discs of zero
thickness. The first and last resonators are connected through cylindrical
openings to semi-infinite cylindrical waveguides. We will consider only
axially symmetric TM fields with $E_{z},E_{r},H_{\varphi}$ components. Time
dependence is $\exp(-i\omega{t})$.

In each resonator we expand the electromagnetic field with the short-circuit
resonant cavity modes

${\vec{E}^{(k)}}=\sum\limits_{q}{e_{q}^{(k)}\,\vec{E}_{q}^{(k)}(\vec{r})},$
(1)
${\vec{H}^{(k)}}=i\sum\limits_{q}{h_{q}^{(k)}\,\vec{H}_{q}^{(k)}(\vec{r})},$
(2)

where $q=\left\\{{0,\,m,\,n}\right\\}$.

$\,\vec{E}_{q}^{(k)},\,\,\vec{H}_{q}^{(k)}\,$ are the solutions of homogenous
Maxwell equations

$\begin{gathered}rot\,\vec{E}_{q}^{(k)}=\,i\,\omega_{q}^{(k)}{\mu_{0}}\vec{H}_{q}^{(k)}\,\,,\hfill\\\
rot\,\vec{H}_{q}^{(k)}=-i\,\omega_{q}^{(k)}{\varepsilon_{0}}\vec{E}_{q}^{(k)}\,\hfill\\\
\end{gathered}$ (3)

with boundary condition $\,{\vec{E}_{\tau}}=0$ on the metal surface

$E_{m,\,n,z}^{(k)}={J_{0}}\left({\frac{{{\lambda_{m}}}}{{{b_{k}}}}r}\right)\cos\left({\frac{\pi}{{{d_{k}}}}n\left({\,z-{z_{k}}}\right)}\right),$
(4)
$H_{m,\,n,\,\varphi}^{(k)}=-i\,\omega_{n,\,m}^{(k)}\frac{{{\varepsilon\varepsilon_{0}}\,{b_{k}}}}{{{\lambda_{m}}}}{J_{1}}\left({\frac{{{\lambda_{m}}}}{{{b_{k}}}}r}\right)\cos\left({\frac{\pi}{{{d_{k}}}}n\,\left({\,z-{z_{k}}}\right)}\right),$
(5)
$E_{m,\,n,\,r}^{(k)}=\frac{{{b_{k}}}}{{{\lambda_{m}}}}\frac{{\pi\,n}}{{{d_{k}}}}{J_{1}}\left({\frac{{{\lambda_{m}}}}{{{b_{k}}}}r}\right)\sin\left({\frac{\pi}{{{d_{k}}}}n\,\left({\,z-{z_{k}}}\right)}\right),$
(6)
$\omega_{0,\,m,\,n}^{(k)2}=\frac{c^{2}}{\varepsilon}\left({{\left({\frac{{{\lambda_{m}}}}{{{b_{k}}}}}\right)}^{2}}+{{\left({\frac{{\pi\,n}}{{{d_{k}}}}}\right)}^{2}}\right),$
(7)

where ${J_{0}}({\lambda_{m}})=0$, $b_{k}$, $d_{k}$ \- k-th resonator radius
and length, $a_{k}$ \- opening radius between k-th and (k-1)-th resonators.

Amplitudes $e_{q}^{\left(k\right)}$ can be found if the tangential electric
fields on the openings are known

$\left({\omega_{q}^{(k)2}-{\omega^{2}}}\right)e_{q}^{(k)}\,\,=\frac{{i\,\omega_{q}^{(k)}}}{{N_{q}^{(k)}}}\left({\oint\limits_{{S_{k}}}{[\vec{E}_{c}^{(k)}\vec{H}_{q}^{(k)*}]\,d\vec{S}+\oint\limits_{{S_{k+1}}}{[\vec{E}_{c}^{(k+1)}\vec{H}_{q}^{(k)*}]\,d\vec{S}}}}\right),$
(8)
$\frac{{\omega_{0,m,n}^{(k)2}}}{{N_{n,\,m}^{(k)}}}=\frac{{2{c^{2}}\,\lambda_{m}^{2}}}{{{\varepsilon_{0}}\,\left|\varepsilon\right|\varepsilon{\sigma_{n}}b_{k}^{4}\pi\,{d_{k}}J_{1}^{2}({\lambda_{m}})\,}}.$

In the semi-infinite waveguides the electromagnetic field can be expanded in
terms of the TM eigenmodes
$\vec{\mathcal{E}}_{s}^{(w,p)},\,\,\vec{\mathcal{H}}_{s}^{(w,p)}$ of a
circular waveguide ($p=1,2$ )

$\begin{gathered}{{\vec{H}}^{(w,p)}}=\sum\limits_{s}{\left({G_{s}^{(p)}\vec{\mathcal{H}}_{s}^{(w,p)}+G_{-s}^{(k)}\vec{\mathcal{H}}_{-s}^{(w,p)}}\right)}\hfill\\\
{{\vec{E}}^{(w,p)}}=\sum\limits_{s}{\left({G_{s}^{(p)}\vec{\mathcal{E}}_{s}^{(w,p)}+G_{s}^{(-p)}\vec{\mathcal{E}}_{-s}^{(w,p)}}\right)}\hfill\\\
\end{gathered}$ (9)

where ${z_{w,1}}={z_{1}},\,\,\,{z_{w,2}}={z_{{N_{R}}+1}}$,

$\mathcal{E}_{s,z}^{(w,p)}={J_{0}}\left({\frac{{{\lambda_{s}}}}{{{b_{w,p}}}}r}\right)\exp\left\\{{\gamma_{s}^{(w,p)}(z-{z_{w,p}})}\right\\},\hfill\\\
$ (10)
$\mathcal{E}_{s,\,r}^{(w,p)}=-\frac{{{b_{w,p}}}}{{{\lambda_{s}}}}\gamma_{s}^{(w,p)}{J_{1}}\left({\frac{{{\lambda_{s}}}}{{{b_{w,p}}}}r}\right)\exp\left\\{{\gamma_{s}^{(w,p)}(z-{z_{w,p}})}\right\\},\hfill\\\
$ (11)
$\mathcal{H}_{s,\varphi}^{(w,p)}=-i\,\omega\frac{{{\varepsilon_{0}}{b_{w,p}}}}{{{\lambda_{s}}}}{J_{1}}\left({\frac{{{\lambda_{s}}}}{{{b_{w,p}}}}r}\right)\exp\left\\{{\gamma_{s}^{(w,p)}(z-{z_{w,p}})}\right\\},\hfill\\\
$ (12)
$\gamma_{s}^{(w,k)2}=\frac{1}{{b_{w,k}^{2}}}\left({\lambda_{s}^{2}-\frac{{b_{w,k}^{2}{\omega^{2}}}}{{{c^{2}}}}}\right).$
(13)

Further consideration will be based on the Moment Method. We shall use the
Bessel functions as the testing functions
(${\psi_{s}}(x)={J_{1}}({\lambda_{s}}x),\,\,x\in[0,1]$) and the complete set
of functions that fulfil the edge condition on the diaphragm rims as the basis
functions (the Meixner basis). We use such Meixner basis [27]

${\varphi_{s}}(r)=2\sqrt{\pi}\frac{{\Gamma(s+1)}}{{\Gamma(s-0.5)}}\frac{1}{{\sqrt{1-{r^{2}}}}}P_{2s-1}^{-1}\left({\sqrt{1-{r^{2}}}}\right),$
(14)

where $P_{n}^{m}\left(x\right)$ are Legendre functions (or spherical
functions) of the first kind,

${\varphi_{s}}\left({\frac{r}{{{a_{k}}}}}\right)\xrightarrow[r\to{a_{k}}]{}\frac{{{C_{s}}}}{{\sqrt{1-{{\left({\frac{r}{{{a}}}}\right)}^{2}}}}}$

Tangential electric fields on the openings we expand into such series

$E_{r}^{(k)}=\sum\limits_{s}{\mathcal{E}_{s}^{(k)}}{\varphi_{s}}\left({\frac{r}{{{a}}}}\right).$
(15)

Using a procedure proposed in [28, 29, 30, 31, 32], we get such matrix
equations

$\left({{T^{(1,1)}}-{\varepsilon^{-1}}{W^{(1)}}}\right){C^{(1)}}-{T^{(1,3)}}{C^{(2)}}=-{R^{(1,1)}}e_{010}^{(1)}+{\varepsilon^{-1}}V$
(16)

$\begin{gathered}\left({{T^{(k,1)}}+{T^{(k,2)}}}\right){C^{(k)}}-{T^{(k,3)}}{C^{(k+1)}}-{T^{(k,4)}}{C^{(k-1)}}=-{R^{(k,1)}}e_{010}^{(k)}+{R^{(k,2)}}e_{010}^{(k-1)},\,\hfill\\\
\,\,\,\,k=2,3,...,{N_{R}}\hfill\\\ \hfill\\\ \end{gathered}$ (17)

$-{T^{({N_{R}}+1,4)}}{C^{({N_{R}})}}+\left({{T^{({N_{R}}+1,2)}}-{\varepsilon^{-1}}{W^{(2)}}}\right){C^{({N_{R}}+1)}}={R^{(1{N_{R}}+1,2)}}e_{010}^{({N_{R}}))},$
(18)

$\begin{gathered}\left({\omega_{010}^{(k)2}-{\omega^{2}}}\right)e_{010}^{(k)}=\hfill\\\
-2\frac{{{b_{k}}}}{{{d_{k}}}}\frac{{\omega_{010}^{(k)2}}}{{\varepsilon
J_{1}^{2}({\lambda_{1}}){\lambda_{1}}}}\left({-\frac{{a_{k}^{2}}}{{b_{k}^{2}}}\sum\limits_{s}{C_{s}^{(k)}{j_{2s-1}}\left({\frac{{{\lambda_{1}}{a_{k}}}}{{{b_{k}}}}}\right)}+\frac{{a_{k+1}^{2}}}{{b_{k}^{2}}}\sum\limits_{s}{C_{s}^{(k+1)}{j_{2s-1}}\left({\frac{{{\lambda_{1}}{a_{k+1}}}}{{{b_{k}}}}}\right)}}\right),\hfill\\\
k=1,2,...,{N_{R}}.\hfill\\\ \end{gathered}$ (19)

Matrices $T,\,W$ and vector $R$ are defined in Appendix 1. We introduce the
fundamental solution of system (16) - (18) as the solution of such difference
matrix equations

$\begin{gathered}\left({{T^{(1,1)}}-{\varepsilon^{-1}}{W^{(1)}}}\right){Y^{(1,{k_{0}})}}-{T^{(1,3)}}{Y^{(2,{k_{0}})}}=-{R^{(1,1)}}{\delta_{1,{k_{0}}}},\hfill\\\
\left({{T^{(k,1)}}+{T^{(k,2)}}}\right){Y^{(k,{k_{0}})}}-{T^{(k,3)}}{Y^{(k+1,{k_{0}})}}-{T^{(k,4)}}{Y^{(k-1,{k_{0}})}}=\hfill\\\
-{R^{(k,1)}}{\delta_{k,{k_{0}}}}+{R^{(k,2)}}{\delta_{k-1,{k_{0}}}},k=2,3,...,{N_{R}},\hfill\\\
-{T^{({N_{R}}+1,4)}}{Y^{({N_{R}},{k_{0}})}}+\left({{T^{({N_{R}}+1,2)}}-{\varepsilon^{-1}}{W^{(2)}}}\right){Y^{({N_{R}}+1,{k_{0}})}}={R^{({N_{R}}+1,2)}}{\delta_{{N_{R}},{k_{0}}}},\hfill\\\
\end{gathered}$ (20)

$\begin{gathered}\left({{T^{(1,1)}}-{\varepsilon^{-1}}{W^{(1)}}}\right){Y^{(1)}}-{T^{(1,3)}}{Y^{(2)}}=={\varepsilon^{-1}}V,\hfill\\\
\left({{T^{(k,1)}}+{T^{(k,2)}}}\right){Y^{(k)}}-{T^{(k,3)}}{Y^{(k+1)}}-{T^{(k,4)}}{Y^{(k-1)}}=0,\hfill\\\
k=2,3,...,{N_{R}},\hfill\\\
-{T^{({N_{R}}+1,4)}}{Y^{({N_{R}})}}+\left({{T^{({N_{R}}+1,2)}}-{\varepsilon^{-1}}{W^{(2)}}}\right){Y^{({N_{R}}+1)}}=0\hfill\\\
\end{gathered}$ (21)

Then the solution of (16) - (18) is written as

$C_{s}^{(k)}=\sum\limits_{{k_{0}}=1}^{{N_{R}}}{Y_{s}^{(k,{k_{0}})}e_{010}^{({k_{0}})}}+Y_{s}^{(k)}$
(22)

Its substitution in (19) gives the basic CCM equation

$\left({\omega_{010}^{(k)2}-{\omega^{2}}}\right)e_{010}^{(k)}\,\,=\omega_{010}^{(k)2}\sum\limits_{{k_{0}}=1}^{{N_{R}}}{{\alpha^{(k,{k_{0}})}}e_{010}^{({k_{0}})}}+\omega_{010}^{(k)2}{F^{(k)}},$
(23)

where

${\alpha^{(k,{k_{0}})}}=-2\frac{{{b_{k}}}}{{{d_{k}}J_{1}^{2}({\lambda_{1}}){\lambda_{1}}}}\left({-\frac{{a_{k}^{2}}}{{b_{k}^{2}}}\sum\limits_{s}{Y_{s}^{(k,{k_{0}})}{j_{2s-1}}\left({\frac{{{\lambda_{1}}{a_{k}}}}{{{b_{k}}}}}\right)}+\frac{{a_{k+1}^{2}}}{{b_{k}^{2}}}\sum\limits_{s}{Y_{s}^{(k+1,{k_{0}})}{j_{2s-1}}\left({\frac{{{\lambda_{1}}{a_{k+1}}}}{{{b_{k}}}}}\right)}}\right)$
(24)
${F^{(k)}}=-2\frac{{{b_{k}}}}{{{d_{k}}J_{1}^{2}({\lambda_{1}}){\lambda_{1}}}}\left({-\frac{{a_{k}^{2}}}{{b_{k}^{2}}}\sum\limits_{s}{Y_{s}^{(k)}{j_{2s-1}}\left({\frac{{{\lambda_{1}}{a_{k}}}}{{{b_{k}}}}}\right)}+\frac{{a_{k+1}^{2}}}{{b_{k}^{2}}}\sum\limits_{s}{Y_{s}^{(k+1)}{j_{2s-1}}\left({\frac{{{\lambda_{1}}{a_{k+1}}}}{{{b_{k}}}}}\right)}}\right)$
(25)

Mode amplitudes in cylindrical waveguides can be found using the following
formulas

$\begin{gathered}G_{-1}^{(1)}=1+2\frac{{a_{1}^{2}{\lambda_{1}}}}{{J_{1}^{2}\left({{\lambda_{s}}}\right)b_{w,1}^{2}\gamma_{1}^{(1)}{b_{w,1}}}}\sum\limits_{s^{\prime}}{C_{s^{\prime}}^{(1)}{j_{2s^{\prime}-1}}\left({\frac{{{\lambda_{1}}{a_{1}}}}{{{b_{w,1}}}}}\right)}\hfill\\\
G_{-s}^{(1)}=2\frac{{a_{1}^{2}{\lambda_{s}}}}{{J_{1}^{2}\left({{\lambda_{s}}}\right)b_{w,1}^{2}\gamma_{s}^{(1)}{b_{w,1}}}}\sum\limits_{s^{\prime}}{C_{s^{\prime}}^{(1)}{j_{2s^{\prime}-1}}\left({\frac{{{\lambda_{s}}{a_{1}}}}{{{b_{w,1}}}}}\right)},\,\,\,\,s=2,3,...\hfill\\\
\end{gathered}$ (26)

$G_{s}^{(2)}=-2\frac{{a_{{N_{R}}+1}^{2}{\lambda_{s}}}}{{J_{1}^{2}\left({{\lambda_{s}}}\right)b_{w,2}^{2}\gamma_{s}^{(2)}{b_{w,2}}}}\sum\limits_{s^{\prime}}{C_{s^{\prime}}^{({N_{R}}+1)}{j_{2s^{\prime}-1}}\left({\frac{{{\lambda_{s}}{a_{{N_{R}}+1}}}}{{{b_{w,2}}}}}\right)},\,\,s=1,2,...$
(27)

The model described above has no limitations. We get the exact coupling matrix
(24) and the vector of distributed external sources (25).

##  3 Infinitive Chain of Resonators

Consider an infinitive (${N_{R}}\rightarrow\infty$) homogeneous
($a_{k}=a,\,b_{k}=b,\,d_{k}=d$) chain. Equations (17),(19) we write as

$e_{010}^{(k)}\,\,=-q\,\sum\limits_{s=1}^{\infty}{{R_{s}}\left({\mathcal{E}_{s}^{(k)}-\mathcal{E}_{s}^{(k+1)}}\right)},$
(28)
$2{T^{(1)}}{\mathcal{E}^{(k)}}-{T^{(2)}}\,{\mathcal{E}^{(k+1)}}-{T^{(2)}}{\mathcal{E}^{(k-1)}}=R\left({e_{010}^{(k)}-e_{010}^{(k-1)}}\right),$
(29)

where

$q=\frac{{\omega_{010}^{2}}}{{\left({\omega_{010}^{2}-{\omega^{2}}}\right)}}\frac{{\,{a^{3}}}}{{{b^{2}}d\,\,J_{1}^{2}({\lambda_{1}})}},$
(30)
${R_{s}}=\frac{1}{2}\frac{{{J_{0}}\left({\frac{{a{\lambda_{1}}}}{b}}\right)}}{{{{\left({\frac{{a{\lambda_{1}}}}{b}}\right)}^{2}}-{{\left({{\lambda_{s}}}\right)}^{2}}}}$
(31)

###  3.1  Systems of matrix equations

Substitution (28) into (29) gives

$2{\bar{T}^{(1)}}{\mathcal{E}^{(k)}}-{\bar{T}^{(3)}}\,{\mathcal{E}^{(k+1)}}-{\bar{T}^{(3)}}{\mathcal{E}^{(k-1)}}=0,$
(32)

where

$\bar{T}_{s,s^{\prime}}^{(1,2)}=\frac{{a^{2}}}{{b^{2}}}\sum\limits_{m}{\bar{\Lambda}_{m}^{(1,2)}}\frac{{{\lambda_{m}}}}{{J_{1}^{2}({\lambda_{m}})}}\frac{{{J_{0}}\left({\frac{{{a}{\lambda_{m}}}}{{{b}}}}\right){j_{2s-1}}\left({\frac{{{\lambda_{m}}{a}}}{{{b}}}}\right)}}{{{{\left({\frac{{{a}{\lambda_{m}}}}{{{b}}}}\right)}^{2}}-\lambda_{s^{\prime}}^{2}}},$
(33)
$\begin{gathered}\bar{\Lambda}_{m}^{(1)}=\frac{{d\,\cosh(d{h_{m}})}}{{bd{h_{m}}\,\sinh(d{h_{m}})}},\hfill\\\
\bar{\Lambda}_{m}^{(2)}=\frac{d}{{bd{h_{m}}\,\sinh(d{h_{m}})}}.\hfill\\\
\end{gathered}$ (34)

It can be shown that the matrix $T^{(2)}$ is invertible and the equation (32)
can be rewritten as

$\bar{T}{\mathcal{E}^{(k)}}=\,{\mathcal{E}^{(k+1)}}+{\mathcal{E}^{(k-1)}},$
(35)

where

$\bar{T}=2{\bar{T}^{(2)-1}}{\bar{T}^{(1)}}.$ (36)

The general solution of the difference equation (32) is [31, 32]

${\mathcal{E}^{(k)}}=\sum\limits_{s=1}^{\infty}{\bar{C}_{s}^{(1)}\bar{\lambda}_{s}^{(1)k}{{\bar{U}}_{s}}}+\sum\limits_{s=1}^{\infty}{\bar{C}_{s}^{(2)}\bar{\lambda}_{s}^{(2)k}{{\bar{U}}_{s}}}$
(37)

where

$\bar{T}{U_{s}}={\bar{\mu}_{s}}{\bar{U}_{s}},$ (38)
$\begin{gathered}\bar{\lambda}_{s}^{2}-{{\bar{\mu}}_{s}}{{\bar{\lambda}}_{s}}+1=0,\hfill\\\
\bar{\lambda}_{s}^{(1,2)}=\frac{{{{\bar{\mu}}_{s}}}}{2}\pm\frac{1}{2}\sqrt{\bar{\mu}_{s}^{2}-4},\hfill\\\
\end{gathered}$ (39)
$\begin{gathered}{{\bar{\mu}}_{s}}=\bar{\lambda}_{s}^{(1)}+\bar{\lambda}_{s}^{(2)},\hfill\\\
\bar{\lambda}_{s}^{(1)}\bar{\lambda}_{s}^{(2)}=1.\hfill\\\ \end{gathered}$
(40)

We used EVCRG programm from IMSL Fortran Numerical Library to calculate
$\bar{\lambda}_{s}$. Results are presented in Tab.1 ( first and second
columns, angles are given in degrees). It can be seen that chosen frequency
f=2.856 GHz lay in the first passband ($E_{01}$) and others modes
($E_{02},E_{03}...$) are evanescent with zero angles. These results are fitted
well the existing theories of periodic waveguides. For the frequency that lay
in the first propagation zone we have two propagating waves and the infinitive
number of evanescent waves with characteristic multipliers with zero phases.
Between the first and the second propagation zones there are only evanescent
oscillations with the phase shift per cell equals $\pi$. The other evanescent
oscillations have zero phase shift. At the end of the second stop band
evanescent $E_{01}$ wave with $\pi$-shift transforms into propagating $E_{02}$
wave with $\pi$-shift too [33]. Therefore, in the second passband $E_{02}$
wave has negative dispersion.

> Table 1 Calculated values of $\bar{\lambda}_{s}$
> ---
> f=2.856 GHz, a=0.99 cm, b=4.08896 cm, d=3.4989 cm
> Matrix $4\times 4$ | Matrix $3\times 3$ | CCM (9) | CCM(7)
> 3.85E+08$\angle 0^{\circ}$ |  | 194.06$\angle 0^{\circ}$ |
> 0.00E+00$\angle 0^{\circ}$ |  |  5.15E-03$\angle 0^{\circ}$ |
> 1.62E+06$\angle 0^{\circ}$ | 1.64E+06$\angle 0^{\circ}$ | 193.67$\angle-89.96^{\circ}$ | 194.06$\angle-59.96^{\circ}$
> 6.18E-07$\angle 0^{\circ}$ | 6.18E-07$\angle 0^{\circ}$ | 5.16E-03$\angle 89.96^{\circ}$ | 5.15E-03$\angle 59.96^{\circ}$
> 6.09E+03 $\angle 0^{\circ}$ | 6.09E+03 $\angle 0^{\circ}$ | 193.67$\angle-89.96^{\circ}$ | 194.06$\angle-59.96^{\circ}$
> 1.64E-04$\angle 0^{\circ}$ | 1.64E-04$\angle 0^{\circ}$ | 5.16E-03$\angle 89.96^{\circ}$ | 5.16E-03$\angle 59.96^{\circ}$
> 1.0$\angle-119.994^{\circ}$ | 1.0$\angle-119.994^{\circ}$ | 1.0$\angle-119.997^{\circ}$ | 1.0$\angle-119.997^{\circ}$
> 1.0$\angle 119.994^{\circ}$ | 1.0$\angle 119.994^{\circ}$ | 1.0$\angle 199.997^{\circ}$ | 1.0$\angle 119.997^{\circ}$

###  3.2  Coupled Cavity Model

The equation (29) we rewrite as

$2{T^{(1)}}{\mathcal{E}^{(k)}}-{T^{(3)}}\,{\mathcal{E}^{(k+1)}}-{T^{(3)}}{\mathcal{E}^{(k-1)}}=R\sum\limits_{j=-\infty}^{\infty}{\left({e_{010}^{(j)}-e_{010}^{(j-1)}}\right){\delta_{k,j}}}=R\sum\limits_{j=-\infty}^{\infty}{e_{010}^{(j)}\left({{\delta_{k,j}}-{\delta_{k,j+1}}}\right)}.$
(41)

Let’s introduce the fundamental solution as the solution of such difference
equation ( $i$ \- the number of diaphragm, $j$ \- the number of cell with a
source)

$2{T^{(1)}}{Y^{(i,j)}}-{T^{(3)}}{Y^{(i+1,j)}}-{T^{(3)}}{Y^{(i-1,j)}}=R{\delta_{i,j}}$
(42)

We will suppose that eigenvalues ${\mu_{s}}$ of the matrix
$T=2{T^{(2)-1}}{T^{(1)}}$ ($T{U_{l}}={\mu_{l}}{U_{l}}$) take such values that
the quantities

$\rho_{s}^{(1,2)}=\frac{{{\mu_{s}}}}{2}\pm\frac{1}{2}\sqrt{\mu_{s}^{2}-4}$
(43)

are the real numbers. In this case, we can impose the condition at infinity

$\left|{{Y^{(i,j)}}}\right|\xrightarrow[i\to\pm\infty]{}0$ (44)

and the solution of (26) we can write as

${Y^{(g,0)}}=\sum\limits_{l=1}^{\infty}{A_{l}}{U_{l}}\left\\{\begin{gathered}\lambda_{l}^{(1)g},\,\,\,g<0,\hfill\\\
\lambda_{l}^{(2)g},\,\,\,g\geqslant 0,\hfill\\\ \end{gathered}\right.$ (45)

where

$\begin{gathered}{A_{l}}=\frac{{{r_{l}}}}{{\lambda_{l}^{(2)}-\lambda_{l}^{(1)}}},\hfill\\\
R=\sum\limits_{l=1}^{\infty}{{r_{l}}{U_{l}}}.\hfill\\\ \end{gathered}$ (46)

The solution of equation (41) can be written as

$\mathcal{E}_{s}^{(k)}=\sum\limits_{j=-\infty}^{j=-\infty}{\left({Y_{s}^{(k,j)}-Y_{s}^{(k-1,j)}}\right)e_{010}^{(j)}}$
(47)

and the equation (28) takes the form

$\begin{gathered}e_{010}^{(k)}\,\,=-q\,\sum\limits_{j=-\infty}^{\infty}{e_{010}^{(j)}\sum\limits_{s=1}^{\infty}{{R_{s}}\left({2Y_{s}^{(k,j)}-Y_{s}^{(k-1,j)}-Y_{s}^{(k+1,j)}}\right)}}\;=\hfill\\\
=-q\,\sum\limits_{j=-\infty}^{\infty}{e_{010}^{(k+j)}\sum\limits_{s=1}^{\infty}{{R_{s}}\left({2Y_{s}^{(k,j+k)}-Y_{s}^{(k-1,j+k)}-Y_{s}^{(k+1,j+k)}}\right)}}.\hfill\\\
\end{gathered}$ (48)

From (45) it follows that ${Y^{(i,j+g)}}={Y^{(i-g,j)}}$,
${Y^{(i,0)}}={Y^{(-i,0)}}$ and we get the final equation of CCM

$\left({\omega_{010}^{2}-{\omega^{2}}}\right)e_{010}^{(k)}=-\omega_{010}^{2}\sum\limits_{k_{0}=-\infty}^{\infty}{e_{010}^{(k+k_{0})}{\alpha^{\left(k_{0}\right)}}}$
(49)

where

$\alpha^{\left(k_{0}\right)}\left(\omega\right)=\frac{{{a^{3}}}}{{{b^{2}}d\,J_{1}^{2}({\lambda_{1}})}}\sum\limits_{s=1}^{\infty}{{R_{s}}\left({2Y_{s}^{(k_{0},0)}-Y_{s}^{(k_{0}-1,0)}-Y_{s}^{(k_{0}+1,0)}}.\right)}$
(50)

It is often assumed that we can neglect ”long couplings”

$\left({\omega_{010}^{2}-{\omega^{2}}}\right)e_{010}^{(k)}=-\omega_{010}^{2}\sum\limits_{k_{0}=-N}^{N}{e_{010}^{(k+k_{0})}{\alpha^{\left(k_{0}\right)}}},$
(51)

It was also assumed [18, 19] that the sum in (50) can be also truncated

${\alpha^{\left(k_{0}\right)}}=\frac{{{a^{3}}}}{{{b^{2}}d\,J_{1}^{2}({\lambda_{1}})}}\sum\limits_{s=1}^{S}{{R_{s}}\left({2Y_{s}^{k_{0}j,0)}-Y_{s}^{(k_{0}-1,0)}-Y_{s}^{(k_{0}+1,0)}}\right)},$
(52)
$\begin{gathered}2T_{s^{\prime},s}^{(1)}Y_{s}^{(i,0)}-T_{s^{\prime},s}^{(2)}Y_{s}^{(i+1,0)}-T_{s^{\prime},s}^{(2)}Y_{s}^{(i-1,0)}={R_{s^{\prime}}}{\delta_{i,0}},\,\,\hfill\\\
i=-I+1,...,0,...I-1;\,\,\,I>N,\,\,\,\,s,s^{\prime}=1,...,S,\hfill\\\
2T_{s^{\prime},s}^{(1)}Y_{s}^{(-I,0)}-T_{s^{\prime},s}^{(2)}Y_{s}^{(-I+1,0)}=0,\hfill\\\
2T_{s^{\prime},s}^{(1)}Y_{s}^{(I,0)}-T_{s^{\prime},s}^{(2)}Y_{s}^{(I-1,0)}=0.\hfill\\\
\end{gathered}$ (53)

Equation (51) is widely used for description of different objects. We made
assumption that $\rho_{s}^{(1,2)}$ are the real numbers. Tab.2 shows that it
is true.

> Table 2 Calculated values of $\mu_{s}$ and $\rho_{s}$
> ---
> f=2.856 GHz, a=0.99 cm, b=4.08896 cm, d=3.4989 cm
> I=6, N=5, S=5
> $\mu_{s}$ | $\rho_{s}^{(2)}$ | $\rho_{s}^{(1)}$
> (-1.93E+02,0.0E+0) | 5.17E-03 | 1.93E+02
> (2.16E+04,0.0E+0) | 4.61E-05 | 2.16E+04
> (2.02E+07,0.0E+0) | 4.93E-08 | 2.02E+07
> (1.26E+10,0.0E+0) | 7.91E-11 | 1.26E+10
> (5.32E+12,0.0E+0) | 1.88E-13 | 5.32E+12

If we seek the solution of (51) as $e_{010}^{(k)}\sim\lambda^{k}$, we get the
characteristic equation

$\sum\limits_{{k_{0}}=-N}^{N}{\lambda^{k_{0}}{\alpha^{\left(k_{0}\right)}\left(\omega\right)}}+\frac{{\left({\omega_{010}^{2}-{\omega^{2}}}\right)}}{{\omega_{010}^{2}}}=0.$
(54)

Solutions of this equation are given in Table 1 (third and fourth columns) for
different N. We see that propagating $E_{01}$ wave has the same phase shift
for different approaches, but other solutions are very different. We obtain
solutions that represent ”non-physical” evanescent waves with phase shifts
that differ from 0 and $\pi$. Amplitudes of these solutions are nearly the
same, but phases strongly depend on the number of interacting oscillators.

To understand the reasons for such solutions, we consider the simplest case
when the tangential electric field is described by one coefficient

$E_{r}^{(k)}={\mathcal{E}_{1}^{(k)}}\varphi_{1}\left(r/a\right).$ (55)

Then (28),(29) transform into

$e_{010}^{(k)}\,\,=-q{R_{1}}\left({\mathcal{E}_{1}^{(k)}-\mathcal{E}_{1}^{(k+1)}},\right)\
$ (56)
${T^{(E)}}\mathcal{E}_{1}^{(k)}={R_{1}}\left({e_{010}^{(k)}-e_{010}^{(k-1)}}\right),$
(57)

where operator
${T^{(E)}}=2T_{11}^{(1)}-T_{11}^{(3)}{\sigma^{(+)}}\,-T_{11}^{(3)}{\sigma^{(-)}}$,
${\sigma^{(+)}},\,\,{\sigma^{(-)}}$ \- shift operators
(${\sigma^{(+)}}{x^{(k)}}={x^{(k+1)}}$,
${\sigma^{(-)}}{x^{(k)}}={x^{(k-1)}}$).

Substitution (56) into (57) gives following equation

$2\left({T_{11}^{(1)}+qR_{1}^{2}}\right)\mathcal{E}_{1}^{(k)}-\left({T_{11}^{(3)}+qR_{1}^{2}}\right)\mathcal{E}_{1}^{(k+1)}\,-\left({T_{11}^{(3)}+qR_{1}^{2}}\right)\mathcal{E}_{1}^{(k-1)}=0.$
(58)

The solution of this equation is

$\mathcal{E}_{1}^{(k)}={C_{1}}\lambda_{1}^{k}+{C_{2}}\lambda_{2}^{k},$ (59)

where $\lambda_{1}$ and $\lambda_{2}$ are the solutions of the characteristic
equation

${\lambda^{2}}-2\frac{{\left({T_{11}^{(1)}+qR_{1}^{2}}\right)}}{{\left({T_{11}^{(3)}+qR_{1}^{2}}\right)}}\lambda+1=0.$
(60)

Solutions to this equation are (f=2.856 GHz, a=0.99 cm, b=4.08896 m, d=3.4989
cm): $\lambda_{1}=1{\angle}119.41^{\circ}$,
$\lambda_{2}=1{\angle}-119.41^{\circ}$ ( compare with results in Table 1).

We can also apply the operator $T^{(E)}$ to the right and left sides of
equation (56) and get the equation like (58)

$2\left({T_{11}^{(1)}+qR_{1}^{2}}\right)e_{010}^{(k)}-\left({T_{11}^{(3)}+qR_{1}^{2}}\right)e_{010}^{(k+1)}\,-\left({T_{11}^{(3)}+qR_{1}^{2}}\right)e_{010}^{(k-1)}=0.$
(61)

The solution of equation (57) can be also written using the solution
$Y_{1}^{(k,j)}$ of the fundamental equation

$2T_{11}^{(1)}Y_{1}^{(k,j)}-T_{11}^{(3)}Y_{1}^{(k+1,j)}\,-T_{11}^{(3)}Y_{1}^{(k-1,j)}={R_{1}}{\delta_{k,j}}.$
(62)

If $\left|{{\chi_{2}}}\right|<1$, where

${\chi_{1,2}}=\frac{{T_{11}^{(1)}}}{{T_{11}^{(3)}}}\pm\sqrt{{{\left({\frac{{T_{11}^{(1)}}}{{T_{11}^{(3)}}}}\right)}^{2}}-1},$
(63)

then the solution of (62) is

$Y_{1}^{(k,j)}=\left\\{\begin{gathered}A\chi_{1}^{k-j},\,k-j<0\,\,,\hfill\\\
A\chi_{2}^{k-j},\,k-j\geqslant 0,\hfill\\\ \end{gathered}\right.$ (64)
$A=\frac{{{R_{1}}}}{{2\left({T_{11}^{(1)}-{\chi_{2}}T_{11}^{(3)}}\right)}}.$
(65)

The solution of (57) is

$\mathcal{E}_{1}^{(k)}=\sum\limits_{j}{Y_{1}^{(k,j)}\left({e_{010}^{(j)}-e_{010}^{(j-1)}}\right)}.$
(66)

Substitution (66) into (56) gives the final equation

$\begin{gathered}e_{010}^{(k)}\,\,=-q{R_{1}}\sum\limits_{j=-\infty}^{\infty}{e_{010}^{(k+j)}\left({2Y_{1}^{(k-j,k)}-Y_{1}^{(k-j-1,k)}-Y_{1}^{(k-j+1,k)}}\right)}=\hfill\\\
=-q{R_{1}}Ae_{010}^{(k)}2\left({1-{\chi_{2}}}\right)-q{R_{1}}A\sum\limits_{j=1}^{\infty}{\left({e_{010}^{(k+j)}+e_{010}^{(k-j)}}\right)\left({2\chi_{2}^{j}-\chi_{2}^{j+1}-\chi_{2}^{j-1}}\right)}=\hfill\\\
=-q{R_{1}}2Ae_{010}^{(k)}\left({1-{\chi_{2}}}\right)-q{R_{1}}A\left({2-{\chi_{2}}-{\chi_{1}}}\right)\sum\limits_{j=1}^{\infty}{\chi_{2}^{j}\left({e_{010}^{(k+j)}+e_{010}^{(k-j)}}\right)}.\hfill\\\
\end{gathered}$ (67)

The characteristic equation of the difference equation (67) coincide with the
equation

$1+q{R_{1}}2A\left({1-{\chi_{2}}}\right)-q{R_{1}}A\frac{{{{\left({1-{\chi_{2}}}\right)}^{2}}}}{{{\chi_{2}}}}\sum\limits_{j=1}^{N}{\left({{\lambda^{j}}+{\lambda^{-j}}}\right)\chi_{2}^{j}}=0,$
(68)

when $N\rightarrow\infty$.

The equation (68) we can rewrite as

$\begin{gathered}{\lambda^{2}}-2\frac{{\left[{T_{11}^{(1)}+qR_{1}^{2}}\right]}}{{\left({T_{11}^{(3)}+qR_{1}^{2}}\right)}}\lambda+1=-qR_{1}^{2}\frac{{\left({1-{\chi_{2}}}\right)}}{{\left({1+{\chi_{2}}}\right)}}\times\hfill\\\
\frac{{\left[{\left({{\chi_{2}}-\lambda}\right)\lambda{{\left({\lambda{\chi_{2}}}\right)}^{N}}+\left({\lambda{\chi_{2}}-1}\right){{\left({{\chi_{2}}/\lambda}\right)}^{N}}}\right]}}{{\left({T_{11}^{(3)}+qR_{1}^{2}}\right)}}.\hfill\\\
\end{gathered}$ (69)

If we suppose that

$\left|{\lambda{\chi_{2}}}\right|<1,\,\,\left|{{\chi_{2}}/\lambda}\right|<1$
(70)

and the number N tends to infinity, the right hand side of the equation (69)
tends to zero and this 2(N+1) order characteristic equation transforms into
the two order equation that coincide with the equation (60).

The equation (69) under the conditions (70) has two solutions very close to
the solutions of the equation (60) , and 2N additional solutions.

For considered above parameter values f=2.856 GHz, a=0.99 cm, b=4.08896 cm,
d=3.4989 cm ${\chi_{2}}$ = 5.03E-003. Therefore, conditions (70) are fulfilled
in passband ($\left|\lambda\right|=1)$ and even in some part of stopband.

A distinctive feature of the difference equation (67) is that despite the
presence of an infinite sum (infinite order), its characteristic equation is
of the second order (at least in some frequency domain).

The above consideration clarifies the reason for the appearance of the
spurious solutions in the main equation of the Coupled Cavity Model (51) -
this is a consequence of the truncation of the infinite sum.

It is well known that without the evanescent eigenfields we cannot get the
correct characteristics of inhomogeneities in waveguides. Since we cannot
operate numerically with infinite sums, using this approach for study non-
homogeneous infinite chains may become problematic because the CCM does not
correctly describe the evanescent eigenfields.

To clarify the correctness of using CCM to study infinite inhomogeneous
resonator chains it is necessary to compare the values of the coupling
coefficients calculated on the basis of the approximate equations (52)-(53)
with the values obtained using the exact equations (16) - (24). The number of
resonators in the finite chain must be sufficient to exclude the influence of
boundaries on the coupling coefficients of cells in the middle of the chain.
We have made a series of calculations for different inhomogeneities inside the
finite chain and the same non-uniform insertions in the infinitive chain.
Their results show that there is such number $N$ in (53) that the coupling
coefficients calculated on the basis of approaches (52) and (24) are almost
the same.

Therefore, the use of the CCM for the numerical study of the characteristics
of infinitive resonator chains gives correct results. The model of an infinite
chain of resonators plays an important role in calculating the electrodynamic
characteristics of inhomogeneities in structured waveguides. Indeed, based on
the CCM of an infinitive chain, one can study the processes of wave
propagation in inhomogeneous resonator chains. This can be done if we assume
that there are homogeneous fragments before and after inhomogeneities. In this
case, in homogeneous fragments at a sufficient distance from the interfaces
(when all evanescent waves are damped), one can look for amplitudes in the
form

$e_{010}^{(k)}\,=\left\\{{\begin{array}[]{*{20}{c}}{\exp\left\\{{i{\varphi_{1}}\left({k-{k_{1}}}\right)}\right\\}+R\exp\left\\{{-i{\varphi_{1}}\left({k-{k_{1}}}\right)}\right\\},\,\,k<{k_{1}},}\\\
{T\exp\left\\{{i{\varphi_{2}}\left({k-{k_{2}}}\right)}\right\\}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,k>{k_{2}},}\end{array}}\right.$
(71)

where $R$ and $T$ are the reflection and transition coefficients. This
approach avoids the use of an ideal load or tuned couplers. This is especially
useful when we are studying frequency dependencies.

But existence of spurious solutions make it difficult (or even impossible) to
use the approximate methods for analysing the chains with slow varying
parameters (WKB approach). As shown by [34], in order to obtain the WKB
equations from the original difference equation, it is necessary to know the
local characteristic multipliers. It is unknown how the spurious solutions
will influence on WKB equations, since their number can be great. But from the
physical point of view, including the spurious solutions in consideration will
be wrong.

It is necessary to use other approaches with correct eigenmode basis for
analysing the chains with slow varying parameters [32, 35, 36].

##  4 Conclusions

Summarizing the above, we can draw the following conclusions:

1\. We can use the CCM to numerical study the characteristics of infinitive
inhomogeneous resonator chains.

2\. We cannot use the CCM to obtain approximate equations (the WKB approach)
for the analysis of chains with slowly varying parameters.

##  5 Appendix 1

$T_{s,s^{\prime}}^{(k,1)}=\frac{{a_{k}^{2}}}{{b_{k}^{2}\,}}\sum\limits_{m}{\Lambda_{m}^{(1,k)}}\frac{{{\lambda_{m}}}}{{J_{1}^{2}({\lambda_{m}})}}\,\frac{{{j_{2s^{\prime}-1}}\left({\frac{{{\lambda_{m}}{a_{k}}}}{{{b_{k}}}}}\right){J_{0}}\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k}}}}}\right)}}{{{{\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k}}}}}\right)}^{2}}-{{\left({{\lambda_{s}}}\right)}^{2}}}}$

$T_{s,s^{\prime}}^{(k,2)}=\frac{{a_{k}^{2}}}{{b_{k-1}^{2}}}\sum\limits_{m}{\Lambda_{m}^{(1,k-1)}}\frac{{{\lambda_{m}}}}{{J_{1}^{2}({\lambda_{m}})}}\,\,\frac{{{j_{2s^{\prime}-1}}\left({\frac{{{\lambda_{m}}{a_{k}}}}{{{b_{k-1}}}}}\right){J_{0}}\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k-1}}}}}\right)}}{{{{\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k-1}}}}}\right)}^{2}}-{{\left({{\lambda_{s}}}\right)}^{2}}}}$

$T_{s,s^{\prime}}^{(k,3)}=\frac{{a_{k+1}^{2}}}{{b_{k}^{2}}}\sum\limits_{m}{\Lambda_{m}^{(2,k)}}\frac{{{\lambda_{m}}}}{{J_{1}^{2}({\lambda_{m}})}}\,\,\frac{{{j_{2s^{\prime}-1}}\left({\frac{{{\lambda_{m}}{a_{k+1}}}}{{{b_{k}}}}}\right){J_{0}}\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k}}}}}\right)}}{{{{\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k}}}}}\right)}^{2}}-{{\left({{\lambda_{s}}}\right)}^{2}}}}$

$T_{s,s^{\prime}}^{(k,4)}=\frac{{a_{k-1}^{2}}}{{b_{k-1}^{2}}}\sum\limits_{m}{\Lambda_{m}^{(2,k-1)}}\frac{{{\lambda_{m}}}}{{J_{1}^{2}({\lambda_{m}})}}\,\frac{{{j_{2s^{\prime}-1}}\left({\frac{{{\lambda_{m}}{a_{k-1}}}}{{{b_{k-1}}}}}\right){J_{0}}\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k-1}}}}}\right)}}{{{{\left({\frac{{{a_{k}}{\lambda_{m}}}}{{{b_{k-1}}}}}\right)}^{2}}-{{\left({{\lambda_{s}}}\right)}^{2}}}}$

$\Lambda_{m}^{(1,k)}=\left\\{{\begin{array}[]{*{20}{c}}{\frac{{{d_{k}}\cosh({d_{k}}h_{m}^{(k)})}}{{{b_{k}}{d_{k}}h_{m}^{(k)}\sinh({d_{k}}h_{m}^{(k)})}}-\frac{{{d_{k}}}}{{{b_{k}}}}\frac{1}{{\left({{\pi^{2}}+d_{k}^{2}h_{\,m}^{(k)2}}\right)}},\,\,\,m=1,}\\\
{\frac{{{d_{k}}\cosh({d_{k}}h_{m}^{(k)})}}{{{b_{k}}{d_{k}}h_{m}^{(k)}\sinh({d_{k}}h_{m}^{(k)})}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,m\neq
1,}\end{array}}\right.$

$\Lambda_{m}^{(2,k)}=\left\\{{\begin{array}[]{*{20}{c}}{\frac{{{d_{k}}}}{{{b_{k}}{d_{k}}h_{m}^{(k)}\sinh({d_{k}}h_{m}^{(k)})}}-\frac{{{d_{k}}}}{{{b_{k}}}}\frac{1}{{\left({{\pi^{2}}+d_{k}^{2}h_{\,m}^{(k)2}}\right)}},\,\,\,m=1,}\\\
{\frac{{{d_{k}}}}{{{b_{k}}{d_{k}}h_{m}^{(k)}\sinh({d_{k}}h_{m}^{(k)})}},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,m\neq
1.}\end{array}}\right.$

$R_{s}^{(k,1)}=\frac{1}{2}\frac{{{J_{0}}\left({\frac{{{a_{k}}{\lambda_{1}}}}{{{b_{k}}}}}\right)}}{{{{\left({\frac{{{a_{k}}{\lambda_{1}}}}{{{b_{k}}}}}\right)}^{2}}-{{\left({{\lambda_{s}}}\right)}^{2}}}},\,\,\,\,R_{s}^{(k,2)}=\frac{1}{2}\frac{{{J_{0}}\left({\frac{{{a_{k}}{\lambda_{1}}}}{{{b_{k-1}}}}}\right)}}{{{{\left({\frac{{{a_{k}}{\lambda_{1}}}}{{{b_{k-1}}}}}\right)}^{2}}-{{\left({{\lambda_{s}}}\right)}^{2}}}},$

$W_{s^{\prime},s}^{(1)}=\frac{{a_{1}^{2}}}{{b_{w,1}^{2}}}\sum\limits_{s^{\prime\prime}}{\,\frac{{{\lambda_{s^{\prime\prime}}}}}{{J_{1}^{2}\left({{\lambda_{s^{\prime\prime}}}}\right)\gamma_{s^{\prime\prime}}^{(1)}{b_{w,1}}}}\frac{{{J_{0}}\left({\frac{{{a_{1}}{\lambda_{s^{\prime\prime}}}}}{{{b_{w,1}}}}}\right){j_{2s-1}}\left({\frac{{{\lambda_{s^{\prime\prime}}}{a_{1}}}}{{{b_{w,1}}}}}\right)}}{{{{\left({\frac{{{a_{1}}{\lambda_{ss^{\prime\prime}}}}}{{{b_{w,1}}}}}\right)}^{2}}-\lambda_{s^{\prime}}^{2}}}},\,\,\,\,\,{V_{s^{\prime}}}=\frac{{{J_{0}}\left({\frac{{{a_{1}}{\lambda_{1}}}}{{{b_{w,1}}}}}\right)}}{{{{\left({\frac{{{a_{1}}{\lambda_{1}}}}{{{b_{w,1}}}}}\right)}^{2}}-\lambda_{s^{\prime}}^{2}}}$

$W_{s^{\prime},s}^{(2)}=\frac{{a_{{N_{R}}+1}^{2}}}{{b_{w,2}^{2}}}\sum\limits_{s^{\prime\prime}}{\,\frac{{{\lambda_{s^{\prime\prime}}}}}{{J_{1}^{2}\left({{\lambda_{s^{\prime\prime}}}}\right)\gamma_{s^{\prime\prime}}^{(2)}{b_{w,2}}}}\frac{{{J_{0}}\left({\frac{{{a_{{N_{R}}+1}}{\lambda_{s^{\prime\prime}}}}}{{{b_{w,2}}}}}\right){j_{2s-1}}\left({\frac{{{\lambda_{s^{\prime\prime}}}{a_{{N_{R}}+1}}}}{{{b_{w,2}}}}}\right)}}{{{{\left({\frac{{{a_{{N_{R}}+1}}{\lambda_{s^{\prime\prime}}}}}{{{b_{w,2}}}}}\right)}^{2}}-\lambda_{s^{\prime}}^{2}}}}$

## References

  * [1] H.A. Bathe.  Theory of diffraction by small holes.  Phys. Rev., 66(7):163–182, 1944.
  * [2] V.V. Vladimirsky.  Coupling hollow electromagnetic resonators through a small hole.  ZhTF, pages 1277–1282, 1947.
  * [3] A.I. Akhiezer and Ya.B. Fainberg.  Slow waves.  UFN, 44(3):321–368, 1951.
  * [4] M.A. Allen and G.S. Kino.  On the theory of strongly coupled cavity chains.  lRE Transactions on Microwave Theory and Techniques, 8(3):362, 1960\.
  * [5] R.M. Bevensee.  Electromagnetic Slow Wave Systems.  John Wiley and Sons, 1964.
  * [6] T. Nishikawa.  Normal mode analysis of standing wave linacs.  IEEE Trans. Nucl. Sci., NS-12(3):630, 1965.
  * [7] T. Nishikawa, S. Giordano, and D. Carter.  Dispersion relation and frequency characteristics of alternating periodic structure for linear accelerators.  Rev. Sci. Instrum., 37(5):652–661, 1966.
  * [8] G.R. Swain.  Circuit analog techniques for analysis of resonantly-coupled linear accelerator structures.  Proceedings of the LINAC1966, pages 125–131, 1966.
  * [9] R. Helm.  Computer study of wave propagation, beam loading and beam blowup in the slac accelerator.  SLAC-PUB-218, 1966.
  * [10] D.E. Nagle, E.A. Knapp, and B.C Knapp.  Coupled resonator model for standing wave accelerator tanks.  Rev. Sci. Instrum., 38(11):1583–1587, 1967.
  * [11] Y. Yamazaki.  Stability of the standing-wave accelerating structure studied with a coupled resonator model.  Particle Accelerators, 32:39–44, 1990.
  * [12] W.E. Jule and R.L. Gluckstern.  Equivalent circuit for a perturbed multi-cell structure.  Proceedings of the PAC73, page 943, 1973.
  * [13] Y. Yamazaki.  Analysis of a standing-wave accelerating structure based on a coupled-resonator model.  KEK Report 99-2, 1999.
  * [14] T. Shintake.  Analysis of the transient respnse in periodic structure based on a coupled-resonator model.  Frontiers of Accelerator Technology, pages 435–454, 1999.
  * [15] J. Gao.  Analytical formulas for the resonant frequency changes due to opening apertures on cavity walls.  Nuclear Instruments and Methods in Physics Research, A311:437–443, 1992.
  * [16] P. Sobenin and B.V. Zverev.  Electrodynamic Characteristics of Accelerating Cavities.  Foundation for International Scientific and Education Cooperation, 1999\.
  * [17] V.V. Paramonov.  Equivalent circuit for a perturbed multi-cell structure.  Proceedings of the LINAC2002, page 413, 2002.
  * [18] M.I. Ayzatsky and V.V. Mytrochenko.  Coupled cavity model based on the mode matching technique.  https://arxiv.org/abs/1505.03223, 2015.
  * [19] M.I. Ayzatsky and V.V. Mytrochenko.  Coupled cavity model for disc-loaded waveguides.  https://arxiv.org/abs/1511.03093, 2015.
  * [20] M.D. Kelisani, S. Doebert, and M. Aslaninejad.  An analytical approach for beam loading compensation and excitation of maximum cavity field gradient in a coupled cavity-waveguide system.  Nuclear Instruments and Methods in Physics Research, A828:132––144, 2016.
  * [21] A.S. Gilmour.  Traveling Wave Tubes, Magnetrons, Crossed-Field Amplifiers, and Gyrotrons.  ARTECH HOUSE, 2011.
  * [22] I.C. Hunter.  Theory and design of microwave filters.  The Institution of Engineering and Technology, 2006.
  * [23] R.J. Cameron.  General coupling matrix synthesis methods for chebyshev filtering functions.  IEEE Transactions on Microwave Theory and Techniques, 47(4):433–442, 1999.
  * [24] Jia-Sheng Hong.  Microstrip Filters for RF/Microwave Applications.  John Wiley and Sons, Inc, 2011.
  * [25] M.Y. Nada, M.A.K. Othman, and F. Capolino.  Theory of coupled resonator optical waveguides exhibiting high-order exceptional points of degeneracy.  Physical Review B, 96:184304, 2013.
  * [26] Hui Liu and Shining Zhu.  Coupled magnetic resonator optical waveguides.  Laser Photonics Rev., 7(6):882–900, 2013.
  * [27] S.I. Eminov.  Meixner basis in the disk diffraction problem.  ZhTF, (3):96–99, 2007.
  * [28] S. Amari, J. Bornemann, and R. Vahldieck.  Accurate analysis of scattering from multiple waveguide discontinuities using the coupled integral equation technique.  J. Electromag. Waves Applicat., 10:1623–1644, 1996.
  * [29] S. Amari, R. Vahldieck, and J. Bornemann.  Analysis of propagation in periodically loaded circular waveguides.  IEE Proceedings Microwaves Antennas and Propagation, 146(1):50–54, 1999.
  * [30] J. Bornemann, S. Amari, and R. Vahldieck.  Analysis of waveguide discontinuities by the coupled-integral equations technique.  Recent Res. Devel. Microwave Theory and Tech., 1:25–35, 1999.
  * [31] M.I. Ayzatsky.  Model of finite inhomogeneous cavity chain and approximate methods of its analysis.  Problems of Atomic Science and Technology, (3):28–37, 2021.
  * [32] M.I. Ayzatsky.  Modelling of inhomogeneous disk-loaded waveguides: matrix difference equations and WKB approximation.  https://arxiv.org/abs/2010.10349, 2020.
  * [33] M.I. Ayzatsky.  Evanescent oscillations in accelerating structures.  https://arxiv.org/abs/physics/0011025, 2000.
  * [34] M.I. Ayzatsky.  Transformation of the linear difference equation into a system of the first order difference equations.  https://arxiv.org/abs/1806.04378, 2018.
  * [35] M.I. Ayzatsky.  Inhomogeneous travelling–wave accelerating sections and WKB approach.  http://arxiv.org/abs/2103.10664, 2021.
  * [36] M.I. Ayzatsky.  Modification of coupled integral equations method for calculation the accelerating structure characteristics.  https://arxiv.org/abs/2203.03518, 2022.