1 Introduction

The \(a_0(980)\) and \(f_0(980)\) scalar mesons were the first observed members of a family of exotic resonances in QCD which are located very close to an inelastic two-particle (or quasi two-particle) threshold (see the review [1]). The \(a_0(980)\) resonance was discovered a long time ago and seen in both the \(K\bar{K}\) and the \(\pi \eta \) channels [2, 3] but its properties are still imprecisely known. This is partly because \(\pi \eta \) production experiments using an \(\eta \) beam, analogous to those which have allowed to determine the \(\pi \pi \) or \(\pi K\) phase shifts, are not feasible. The \(a_0(980)\) properties have to be determined solely from final-state rescattering effects. As shown by Flatté [4], who proposed a simple replacement for the Breit-Wigner resonance formula accounting for the \(\pi \eta -K\bar{K}\) coupled-channel dynamics, ambiguous results for the \(a_0\) width can be obtained. In the PDG [5], indeed, its width is simply quoted as lying in a range from 50 to 100 MeV. Beyond the value of the width, one would like to determine the positions of the poles on the unphysical Riemann sheets. For resonance states close to an inelastic threshold, there are three sheets (II, III, IV in the case of the \(a_0\)) which are physically relevant (i.e. a pole in one of these can be close to the physical region). It is also clear that a better determination of the resonance properties is closely tied to a better knowledge of the physical scattering T-matrix.

A theoretically motivated treatment of final-state interactions becomes prohibitively difficult for multiparticle final states. In this regard, photon–photon to meson-meson scattering amplitudes are very favourable processes.Footnote 1 They are free of initial-state interactions and satisfy dispersion relations which can be constrained, in the case of \(\pi \pi \) or \(\pi \eta \) by both soft photon [7, 8] and soft pion [9] low-energy theorems. As was illustrated in the seminal papers [10, 11] a predictive representation can be implemented at the level of the partial-waves with a simple modelling of the left-hand cut. From an experimental point of view \(\gamma \gamma \rightarrow \eta \pi ^0\) cross-sections were first measured by the Crystal Ball collaboration [12]. Measurements with much higher statistics were recently performed by the Belle collaboration [13]. There are also experimental results on the \(\eta \rightarrow \pi ^0\gamma \gamma \) decay amplitude [14, 15]. Recent experimental data exist also for \(\gamma \gamma \rightarrow {K^0}{\bar{K}^0}\) [16] and \(\gamma \gamma \rightarrow {K^+}{K^-}\) [17] which will be considered in our study.

There are two puzzling aspects in the data analysis performed by the Belle collaboration which we wish to reconsider: a) they find that the \(a_0(980)\) peak seems to be best described by an ordinary (i.e. essentially elastic) Breit-Wigner function and b) they find an excited \(a_0\), which could correspond to the \(a_0(1450)\), but has a width, \(\varGamma =65^{+2.1}_{-5.4}\) MeV, much smaller than the PDG average as well as a significantly smaller mass. These data have been re-analysed recently [18] based on a specific meson-meson T-matrix model [19] applied to the \(\pi \eta \) S-wave, from which the corresponding \(\gamma \gamma \) amplitude is deduced from a Muskhelishvili-Omnès (MO) construction [20, 21]. Including also the \(a_2(1320)\) resonance but no \(a_0(1450)\) resonance a good description of the data up to 1.1 GeV and a qualitative description up to 1.4 GeV has been achieved. The \(a_0(980)\) pole, in this model, lies on the fourth Riemann sheet. Interestingly, an analogous pole location was found in a lattice QCD calculation of the T-matrix [22] who implemented a coupled-channel generalisation of Lüscher’s single-channel method [23]. This result, however, corresponds to a pion mass \({m_{\pi }}=391\) MeV and it is not known how the pole location evolves upon varying \({m_{\pi }}\) (see, however, Ref.  [24]). Fits to the Belle data with a conventional \(a_0(1450)\) were performed in Ref.  [25] but only the integrated cross-sections were considered in that work. Detailed global descriptions of photon–photon scattering to two mesons were proposed also in Ref. [26] focusing mostly on the \(I=0,2\) channels. In the \(I=1\) sector, they considered the \(K\bar{K}\) channel but not \(\pi \eta \).

We perform here a global fit which takes into account both \(\pi \eta \) and \(K\bar{K}\) photon–photon data including all the differential cross-sections up to 1.4 GeV. In the \(I=1\) sector, we use the S-wave coupled-channel T-matrix model developed in Ref.  [27] which satisfies unitarity, proper analyticity properties and matches to the chiral expansion up to the next-to-leading order at low energy. The S-wave photon–photon amplitudes are then deduced from a general MO representation involving two subtraction constants and implementing a simple description of the left-hand cut from cross-channel vector-meson exchanges. This is quite similar to Ref.  [18], we differ mainly by using SU(2) chiral symmetry, which allows to fix one of the subtraction constants through a soft pion theorem.Footnote 2 The \(J=2\) partial-waves are described more phenomenologically as a sum of cross-channel resonance exchange and a direct \(a_2(1320)\) Breit-Wigner amplitude. The \(\gamma \gamma \rightarrow (K\bar{K})_{I=1}\) amplitudes are then combined with \(I=0\) amplitudes taken from a previous work [28] which considered \(\gamma \gamma \rightarrow \) \((\pi \pi )_{I=0,2}\), \((K\bar{K})_{I=0}\) in order to reconstruct the physical \({K^+}{K^-}\), \({K^0}{\bar{K}^0}\) amplitudes. A global fit of photon–photon to meson-meson data was performed some time ago based on unitarised chiral amplitudes [29], this model was also applied to the \(\eta \rightarrow {\pi ^0}\gamma \gamma \) decay [30, 31]. A remarkable qualitative agreement with the data available at the time was achieved using a single arbitrary parameter in the S-wave. Today, a much larger data set is available and the precision has increased significantly. We use here a model for the T-matrix involving six parameters which will be determined by performing fits to these data.

An interesting aspect of the \(\gamma \gamma \rightarrow \pi \eta \) amplitudes is that they can be probed experimentally both in the scattering regime: \(s_{\gamma \gamma } \ge (m_\eta +{m_{\pi }})^2\) and in the decay regime: \(0 \le s_{\gamma \gamma } \le (m_\eta -{m_{\pi }})^2 \). In the decay region the amplitudes are largely dominated by the light vector meson exchanges in the cross-channels: \(\gamma \pi \rightarrow \rho ,\omega \rightarrow \gamma \eta \), which were first computed in Ref.  [32]. In this low-energy region the rescattering contributions can be estimated in the SU(3) chiral expansion [33]. They proceed essentially via \(\gamma \gamma \rightarrow {K^+}{K^-}\rightarrow {\pi ^0}\eta \) but, at low energy, the \({\pi ^+}{\pi ^-}\) contribution (which is isospin violating) is not negligible and must also be included [33]. We will present a dispersive calculation of this contribution which gives rise to a cusp in the energy distribution \(d\varGamma ^{\eta \rightarrow {\pi ^0}\gamma \gamma }/ds_{\gamma \gamma }\).

The plan of the paper is as follows. After recalling some general properties of photon–photon amplitudes and introducing the notation in Sect. 2, we write the unitarity relations for the partial-waves in Sect. 3 and present the modelling of the left-hand cut of these. In Sect. 4 we recall the derivation of a MO representation for the S-waves. A dispersive representation for the isospin violating S-wave, valid at low energy, is also derived. The comparison with the experimental data on photon–photon scattering is performed in Sect. 5 and the information that can be deduced on the \(a_0\) resonances are discussed.

2 Basic ingredients

2.1 Kinematics

We consider two-photon to two-meson scattering amplitudes

$$\begin{aligned} \gamma (q_1)\gamma (q_2) \rightarrow M_1(p_1) M_2(p_2) \end{aligned}$$
(1)

with \(M_1M_2= \pi \eta \) or \(K\bar{K}\). The Mandelstam variables are defined as usual by

$$\begin{aligned} s=(q_1+q_2)^2,\quad t=(q_1-p_1)^2,\quad u=(q_1-p_2)^2\ . \end{aligned}$$
(2)

The various physical regions in the s, \(t-u\) Mandelstam plane for the \(\gamma \gamma \rightarrow \pi \eta \), \(\gamma \pi \rightarrow \gamma \eta \) and \(\eta \rightarrow \gamma \gamma \pi \) processes are shown on Fig. 1. In the centre-of-mass frame of the two photons the momenta are expressed as

$$\begin{aligned} q_1= & {} \dfrac{\sqrt{s}}{2}\begin{pmatrix} 1\\ \hat{z}\\ \end{pmatrix}, \quad q_2=\dfrac{\sqrt{s}}{2}\begin{pmatrix} 1\\ -\hat{z}\\ \end{pmatrix} \nonumber \\ p_1= & {} \dfrac{1}{2\sqrt{s}}\begin{pmatrix} s+\varDelta _{12}\\ \sqrt{\lambda _{12}(s)}\,\hat{v}\\ \end{pmatrix},\quad p_2=\dfrac{1}{2\sqrt{s}}\begin{pmatrix} s-\varDelta _{12}\\ -\sqrt{\lambda _{12}(s)}\,\hat{v}\\ \end{pmatrix}\nonumber \\ \end{aligned}$$
(3)

where

$$\begin{aligned} \varDelta _{12}= & {} m_1^2-m_2^2\nonumber \\ \lambda _{12}(s)= & {} (s-(m_1+m_2)^2)(s-(m_1-m_2)^2)\ . \end{aligned}$$
(4)

Taking \(\hat{z}\) to a unit vector along the z axis and \(\hat{v}\) to be a unit vector with polar angles \(\theta \), \(\phi \) such that \(\hat{z}\cdot \hat{v}=\cos \theta \), we can express t, u in terms of \(\cos \theta \) as,

$$\begin{aligned} t= & {} \dfrac{1}{2}\left( m_1^2+m_2^2-s + \sqrt{\lambda _{12}(s)}\cos \theta \right) \nonumber \\ u= & {} \dfrac{1}{2}\left( m_1^2+m_2^2-s - \sqrt{\lambda _{12}(s)} \cos \theta \right) \ . \end{aligned}$$
(5)
Fig. 1
figure 1

Physical regions for \(\gamma \gamma \rightarrow \pi \eta \), \(\gamma \pi \rightarrow \gamma \eta \) scattering and \(\eta \rightarrow \gamma \gamma \pi ^0\) decay. The black square indicates the soft pion point

The polarisation vectors of the two photons \(\epsilon _1(q_1,\lambda )\), \(\epsilon _2(q_2,\lambda ')\) in this frame read

$$\begin{aligned} \epsilon _1(q_1,\lambda )=\frac{1}{\sqrt{2}}\begin{pmatrix} 0 \\ -\lambda \\ -i \\ 0 \\ \end{pmatrix},\quad \epsilon _2(q_2,\lambda ')=\frac{1}{\sqrt{2}}\begin{pmatrix} 0 \\ \lambda ' \\ -i \\ 0 \\ \end{pmatrix} \end{aligned}$$
(6)

and they satisfy \(\epsilon _1\cdot q_1 = \epsilon _2\cdot q_2 = \epsilon _1\cdot q_2= \epsilon _2\cdot q_1=0\ . \)

2.2 Amplitudes

We denote the \(\gamma \gamma \rightarrow \eta \pi \) helicity amplitude as \(L_{\lambda \lambda '}\)

(7)

(we have factored out the electric charge e and the dependence on the azimuthal angle \(\phi \)) while for the \(\gamma \gamma \rightarrow K\bar{K}\) amplitudes we use the same notation as in previous work: \(K^c_{\lambda \lambda '}(s,t)\) for charged kaons and \(K^n_{\lambda \lambda '}(s,t)\) for neutral kaons. Helicity amplitudes are convenient for performing the partial wave expansion [34]. They can be expressed in terms of tensor amplitudes using the reduction formulas, e.g.

$$\begin{aligned}&\hbox {e}^{i(\lambda -\lambda ')\phi }L_{\lambda \lambda '}(s,t)\nonumber \\&\quad = \epsilon _1^\mu (q_1,\lambda )\epsilon _2^\nu (q_2,\lambda ') W_{\mu \nu }(q_1,q_2,p_1,p_2)\ . \end{aligned}$$
(8)

By gauge invariance, the tensor amplitude \(W_{\mu \nu }\) must satisfy the two Ward identities

$$\begin{aligned} q_1^\mu W_{\mu \nu }= q_2^\nu W_{\mu \nu }=0\ . \end{aligned}$$
(9)

One can form two independent tensors which satisfy (9) which we take as

$$\begin{aligned} T_1^{\mu \nu }= & {} \dfrac{1}{2} s g^{\mu \nu }-q_1^\nu q_2^\mu \nonumber \\ T_2^{\mu \nu }= & {} 2s\varDelta ^\mu \varDelta ^\nu +4 q_1.\varDelta \,q_2.\varDelta \, g^{\mu \nu }\nonumber \\&-4q_2.\varDelta \, q_1^\nu \varDelta ^\mu -4q_1.\varDelta \, q_2^\mu \varDelta ^\nu \end{aligned}$$
(10)

with

$$\begin{aligned} \varDelta =p_1-p_2 \ . \end{aligned}$$
(11)

The tensor amplitude \(W^{\mu \nu }\) can then be expressed in terms of two scalar amplitudes A, B

$$\begin{aligned} W^{\mu \nu }(q_1,q_2,p_1,p_2)= A(s,t,u) T_1^{\mu \nu } + B(s,t,u) T_2^{\mu \nu } \end{aligned}$$
(12)

which satisfy dispersion relations. Using Eqs. (5), (6) one can easily express the helicity amplitudes \(L_{\lambda \lambda '}\) in terms of the two scalar amplitudes,

$$\begin{aligned} L_{++}= & {} L_{--}=\dfrac{s}{2}\,A(s,t)+s\left( 2m_1^2+2m_2^2-s\right) \,B(s,t)\nonumber \\ L_{+-}= & {} L_{-+}=\sin ^2\theta \lambda _{12}(s) B(s,t)\ . \end{aligned}$$
(13)

Assuming unpolarised photon beams the differential cross-section reads,

$$\begin{aligned} \frac{d\sigma ^{\gamma \gamma \rightarrow M_1M_2}}{d\cos \theta }=\frac{\pi \alpha ^2}{4 s^2} \sqrt{\lambda _{12}(s)} \left( \vert L_{++}\vert ^2 + \vert L_{+-}\vert ^2\right) \ . \end{aligned}$$
(14)

Concerning the \(\eta \rightarrow \pi \gamma \gamma \) decay amplitude, the double differential distribution in the Dalitz plot reads,

$$\begin{aligned} \frac{d^2\varGamma ^{\eta \rightarrow \gamma \gamma \pi } }{ds dt}=\frac{\alpha ^2}{8\pi m_\eta ^3} \left( \vert L_{++}\vert ^2 + \vert L_{+-}\vert ^2\right) \ \end{aligned}$$
(15)

and the distribution as a function of s only (which is the one available experimentally) is given by

$$\begin{aligned} \dfrac{d\varGamma ^{\eta \rightarrow \gamma \gamma \pi } }{ds}= & {} \dfrac{\alpha ^2}{32\pi m_\eta ^3} \, \sqrt{\lambda _{12}(s)} \nonumber \\&\times {\displaystyle \int _{-1}^1}d\cos \theta \left( \vert L_{++}\vert ^2 + \vert L_{+-}\vert ^2\right) \ . \end{aligned}$$
(16)

2.3 Isospin

Using the following (usual) isospin assignments for the pions and the kaons

$$\begin{aligned} \begin{pmatrix} {\pi ^+}\\ {\pi ^0}\\ {\pi ^-}\\ \end{pmatrix}\sim & {} \begin{pmatrix} -{\vert 11 \rangle }\\ {\vert 10 \rangle }\\ {\vert 1,-1 \rangle }\\ \end{pmatrix}\nonumber \\ \begin{pmatrix} {K^+}\\ {K^0}\\ \end{pmatrix}\sim & {} \begin{pmatrix} {\vert {{1\over 2}{1\over 2}} \rangle }\\ {\vert {{1\over 2}{-1\over 2}} \rangle }\\ \end{pmatrix},\quad \begin{pmatrix} {\bar{K}^0}\\ {K^-}\\ \end{pmatrix}\sim \begin{pmatrix} {\vert {{1\over 2}{1\over 2}} \rangle }\\ -{\vert {{1\over 2}{-1\over 2}} \rangle }\\ \end{pmatrix} \end{aligned}$$
(17)

while \(\eta \sim {\vert 0,0 \rangle }\), the relations between the amplitudes \(\gamma \gamma \rightarrow {K^+}{K^-}\), \({K^0}{\bar{K}^0}\) and the isospin amplitudes \(\gamma \gamma \rightarrow (K \bar{K})_{I=0,1}\) read,

$$\begin{aligned} \left( \begin{array}{l} K^0_{\lambda \lambda '}\\ K^1_{\lambda \lambda '} \end{array}\right) = \left( \begin{array}{rr} -\sqrt{\frac{1}{2}} &{} -\sqrt{\frac{1}{2}}\\ -\sqrt{\frac{1}{2}} &{} \sqrt{\frac{1}{2}} \end{array}\right) \left( \begin{array}{l} K^c_ {\lambda \lambda '}\\ K^n_ {\lambda \lambda '} \end{array}\right) \ , \end{aligned}$$
(18)

and the analogous relations between the amplitudes \(\gamma \gamma \rightarrow {\pi ^+}{\pi ^-}\), \({\pi ^0}{\pi ^0}\) and the corresponding isospin \(I=0,2\) amplitudes (which will also be needed) is

$$\begin{aligned} \left( \begin{array}{l} H^0_{\lambda \lambda '}\\ H^2_{\lambda \lambda '} \end{array}\right) = \left( \begin{array}{rr} -\sqrt{\frac{2}{3}} &{} -\sqrt{\frac{1}{3}}\\ -\sqrt{\frac{1}{3}} &{} \sqrt{\frac{2}{3}} \end{array}\right) \left( \begin{array}{r} \sqrt{2}\,H^c_ {\lambda \lambda '}\\ H^n_ {\lambda \lambda '} \end{array}\right) \ . \end{aligned}$$
(19)

3 Partial waves: unitarity, analyticity

3.1 Right-hand cut and unitarity relations

It is convenient to collect the three \(I=1\) scattering amplitudes \(\pi \eta \rightarrow \pi \eta \), \(\pi \eta \rightarrow K\bar{K}\) and \(K\bar{K}\rightarrow K\bar{K}\) into a \(2\times 2\) matrix

$$\begin{aligned} \varvec{T}\equiv \begin{pmatrix} T^{\pi \eta \rightarrow \pi \eta } &{} T^{\pi \eta \rightarrow K\bar{K}}\\ T^{\pi \eta \rightarrow K\bar{K}} &{} T^{ K\bar{K}\rightarrow K\bar{K}}\\ \end{pmatrix} \end{aligned}$$
(20)

and we can define the partial wave expansions as

$$\begin{aligned} \varvec{T}(s,t,u)=16\pi \sum (2j+1)\varvec{T}_j(s) P_j(\cos \theta ) \end{aligned}$$
(21)

where \(\theta \) is the scattering angle in the centre-of-mass system. The unitarity relation for the partial waves are easily derived and reads

$$\begin{aligned} \mathrm{Im}\,[\varvec{T}_j(s)]= \varvec{T}_j(s)\varvec{\varSigma (s)} \varvec{T}^*_j(s) \end{aligned}$$
(22)

with

$$\begin{aligned} \varvec{\varSigma (s)} =\begin{pmatrix} \dfrac{ \sqrt{\lambda _{\pi \eta }(s)}}{s}\,\theta (s-({m_{\pi }}+m_\eta )^2) &{} 0\\ 0 &{} \dfrac{ \sqrt{\lambda _{KK}(s)}}{s}\,\theta (s-4m_K^2)\\ \end{pmatrix}.\nonumber \\ \end{aligned}$$
(23)

Concerning the \(I=1\) photon–photon amplitudes we define the partial-wave expansion as

$$\begin{aligned} \left( \begin{array}{l} L_{\lambda \lambda '}(s,t)\\ K^1_{\lambda \lambda '}(s,t) \end{array}\right) =\sum _j (2j+1) \left( \begin{array}{l} l_{j,\lambda \lambda '}(s)\\ k^1_{j,\lambda \lambda '}(s) \end{array}\right) \,d^j_{\lambda -\lambda ',0}(\theta )\ .\nonumber \\ \end{aligned}$$
(24)

The unitarity relations for the S-waves, as will be implemented below, read

$$\begin{aligned} \mathrm{Im}\,\begin{pmatrix} l_{0++}(s)\\ k^1_{0++}(s)\\ \end{pmatrix}= \varvec{T}^*_j(s)\varvec{\varSigma (s)} \begin{pmatrix} l_{0++}(s)\\ k^1_{0++}(s)\\ \end{pmatrix}\ , \end{aligned}$$
(25)

which also give the discontinuities of the partial-waves (extended to complex values of s) across the right-hand cut.

3.2 \(\pi \pi \) (isospin violating) contribution in S-wave unitarity

The unitarity relation written above (25) collects the contributions from the \(\pi \eta \) and the \(K\bar{K}\) states. Let us also consider the contribution from the \(\pi ^+\pi ^-\) state, which has the form

$$\begin{aligned} \begin{array}{ll} \mathrm{Im}\,[{\langle \pi ^0\eta \vert {{\mathcal {T}}} \vert \gamma \gamma \rangle }]_{\pi \pi }= &{} \dfrac{1}{2}\displaystyle \int d\varPhi (p_+,p_-)\, {\langle \pi ^0\eta \vert {{\mathcal {T}}} \vert {\pi ^+}{\pi ^-} \rangle }\\ &{} \times {\langle {\pi ^+}{\pi ^-}\vert {{\mathcal {T}}}^\dagger \vert \gamma \gamma \rangle } \end{array}\nonumber \\ \end{aligned}$$
(26)

where \(d\varPhi \) is the phase-space integration measure. This contribution, being proportional to the \(\pi \pi \rightarrow \eta \pi \) amplitude, is isospin violating but it is enhanced at low energy due the large size of the \(\gamma \gamma \rightarrow \pi ^+\pi ^-\) amplitude [33]. The ChPT evaluation performed in Ref.  [33] amounts to using the \(O(p^2)\) tree-level amplitudes for both \(\pi \pi \rightarrow \eta \pi \) and \(\gamma \gamma \rightarrow \pi ^+\pi ^-\) in Eq. (26). An evaluation which goes beyond the chiral expansion can be performed which we now discuss.

We consider only the S-wave contribution in Eq. (26) and a restricted kinematical region such that \(s,t, u < 1\) \(\hbox {GeV}^2\). In such a region, the \({\pi ^+}{\pi ^-}\rightarrow \eta {\pi ^0}\) amplitude can be approximated in terms of three one-variable functions \(M_0\), \(M_1\), \(M_2\) (see [36, 37]),

$$\begin{aligned} T^{{\pi ^+}{\pi ^-}\rightarrow \eta {\pi ^0}}(s,t,u)= & {} -\epsilon _L\Bigg [M_0(s){-}\frac{2}{3}M_2(s){+}(s{-}u)M_1(t)\nonumber \\&+(s{-}t)M_1(u){+}M_2(t){+}M_2(u)\Bigg ]\ .\nonumber \\ \end{aligned}$$
(27)

An isospin violating parameter \(\epsilon _L\) has been factorised which may be taken as [38]

$$\begin{aligned} \epsilon _L=\dfrac{ \left( {m}^2_{K^0}- {m}^2_{K^+}\right) _{QCD}}{3\sqrt{3} F_\pi ^2}\ . \end{aligned}$$
(28)

The three \(M_I\) functions obey a set of coupled Khuri-Treiman integral equations, see Ref.  [38] for a complete review of work on this subject. We will use here the evaluation of the \({K^0}-{K^+}\) QCD mass difference from Ref.  [39] (updated in [38]) based on experimental data on \(\eta \rightarrow 3\pi \) decays: \(\left( {m}^2_{K^0}-{m}^2_{K^+}\right) _{QCD}=(6.24\pm 0.38)\times 10^{-3}\) \(\hbox {GeV}^2\), which gives

$$\begin{aligned} \epsilon _L=0.141\pm 0.009\ . \end{aligned}$$
(29)

The amplitudes corresponding to a given \(\pi \pi \) isospin state \(I,I_z\)

$$\begin{aligned} {{\mathcal {M}}}^{I\,I_z}\equiv {\langle \eta \pi \vert T\vert \pi \pi ;I\,I_z \rangle } \end{aligned}$$
(30)

are easily expressed using crossing symmetry and the Wigner-Eckart theorem. In the unitarity relation (26) the amplitudes with \(I=0,2\), \({{\mathcal {M}}}^{00}\) and \({{\mathcal {M}}}^{20}\), are needed which have the following expressions

(31)

in terms of the \(M_I\) functions. We denote the angular integrals of these amplitudes as

$$\begin{aligned}&\frac{1}{2}\displaystyle \int _{-1}^1 dz\, {{\mathcal {M}}}^{00}(s,t,u) \equiv \sqrt{3}\,\epsilon _L \left( M_0(s)+ \hat{M}_0(s)\right) \nonumber \\&\frac{1}{2}\displaystyle \int _{-1}^1 dz\, {{\mathcal {M}}}^{20}(s,t,u) \equiv -\frac{2\sqrt{6}}{3}\,\epsilon _L \left( M_2(s)+\hat{M}_2(s)\right) \ .\nonumber \\ \end{aligned}$$
(32)

With this notation, the \(\pi \pi \) contribution to the unitarity relation is finally expressed as follows,

$$\begin{aligned}&\mathrm{disc}\,[{l}_{0,++}(s)]_{\pi \pi } =\epsilon _L \dfrac{\sqrt{3}}{32\pi }\sqrt{\dfrac{s-4m_\pi ^2}{s}}\nonumber \\&\qquad \times \bigg [ \left( h^{0}_{0,++}(s)\right) ^*(M_0(s)+\hat{M}_0(s))\nonumber \\&\qquad -\dfrac{2\sqrt{2}}{3} \left( h^{2}_{0,++}(s)\right) ^*(M_2(s)+\hat{M}_2(s)) \bigg ] \end{aligned}$$
(33)

where

$$\begin{aligned} \mathrm{disc}\,[{l}_{0,++}(s)]_{\pi \pi }\equiv \frac{{l}_{0,++}(s+i\epsilon ) -{l}_{0,++}(s-i\epsilon )}{2i}\ , \end{aligned}$$
(34)

and \(h^I_{0,++}(s)\) are the two S-wave \(\gamma \gamma \rightarrow (\pi \pi )^I\) amplitudes with \(I=0,2\). This \(\pi \pi \) discontinuity of \({l}_{0++}\) can be estimated from Eq. (33) using inputs from Ref.  [28] for \(\gamma \gamma \rightarrow \pi \pi \) and from [27] for \(\pi \pi \rightarrow \eta \pi \). The result of this estimate is illustrated on Fig. 2 and compared with the chiral calculation at NLO. The dispersive evaluation displays a square-root singularity at \(s=(m_\eta -{m_{\pi }})^2\) induced by the endpoint of the left-hand cut in the functions \(\hat{M}_0\), \(\hat{M}_2\) which overlaps with the right-hand cut as a result of the instability of the \(\eta \). As a further consequence, the phase of the partial-waves \(M_I+\hat{M}_I\) violate Watson’s theorem and do not cancel with the phases of \(h_{0,++}^I\) such that the discontinuity has both a real and an imaginary part.

Fig. 2
figure 2

Discontinuity of the \(\gamma \gamma \rightarrow \pi \eta \) S-wave amplitude (real part) across the \({\pi ^+}{\pi ^-}\) cut computed using dispersive results for the \(\gamma \gamma \rightarrow \pi \pi \) and \(\pi \pi \rightarrow \eta \pi \) amplitudes, compared to the chiral \(O(p^4)\) result

3.3 Left-hand cut: born amplitudes

A left hand cut in the \(\gamma \gamma \) partial-waves is generated by singularities in the cross-channels \(\gamma P_1 \rightarrow \gamma P_2\). In the case when \(P_1=P_2={K^+}\) or \({\pi ^+}\) the leading singularity is the kaon or pion pole and the corresponding (so-called Born) amplitudes read

$$\begin{aligned} A_P^{Born}(s,t,u)= & {} \dfrac{s}{\left( t-m^2_P\right) \left( u-m^2_P\right) }\nonumber \\ B_P^{Born}(s,t,u)= & {} \dfrac{1}{2\left( t-m^2_P\right) \left( u-m^2_P\right) } \end{aligned}$$
(35)

with \(P={K^+},\ {\pi ^+}\). We will need the \(I=1\) component of the \(K^+\) Born amplitude projected on the S-wave which reads,

$$\begin{aligned} k^{1,Born}_{0,++}(s)= & {} -\dfrac{2\sqrt{2}\, m^2_{{K^+}}}{s}\, L_{K^+}(s)\ ,\nonumber \\ L_{K^+}(s)= & {} \dfrac{1}{\beta _{{K^+}}(s)}\log \dfrac{1+\beta _{{K^+}}(s)}{1-\beta _{{K^+}}(s)} \end{aligned}$$
(36)

with \(\beta _P(s)=\sqrt{1-4m^2_P/s}\). We also recall the expressions of the \(J=2\) Born helicity amplitudes

$$\begin{aligned} k^{1,Born}_{2++}(s)&= -\dfrac{2m^2_{K^+}}{s\,\beta ^2_{K^+}(s)}\left[ \left( \beta ^2_{K^+}(s)-3\right) \,L_{K^+}(s)+6 \right] \nonumber \\ k^{1,Born}_{2+-}(s)&= \dfrac{\sqrt{6}}{4s \beta ^2_{K^+}(s)}\bigg [ \left( 1-\beta ^2_{K^+}(s)\right) ^2\,L_{K^+}(s)\nonumber \\&\quad +\dfrac{10}{3}\beta ^2_{K^+}(s)-2\bigg ]\ . \end{aligned}$$
(37)

3.4 Left-hand cut: vector meson exchanges

Leading contributions to the left-hand cut in the \({\pi ^0}\eta \) amplitudes are modelled from the \(\rho \), \(\omega \), \(\phi \) vector meson exchanges. We define the \(VP\gamma \) coupling constants \(G_{VP}\) through simple Lagrangians

$$\begin{aligned} {{\mathcal {L}}}_{VP\gamma }= e G_{VP}\, \epsilon ^{\mu \nu \alpha \beta } F_{\mu \nu }\partial _\alpha {P}V_\beta \end{aligned}$$
(38)

from which one can relate their values to the decay widths of the vector mesons via

$$\begin{aligned} \varGamma _{V\rightarrow P\gamma }= \alpha C^2_{VP} \frac{ \left( m_V^2-m_P^2\right) ^3}{6 m_V^3}\ . \end{aligned}$$
(39)

The vector-exchange contributions to the \(\gamma \gamma \rightarrow P_1 P_2\) amplitudes are easily computed,

$$\begin{aligned} A_V^{P_1P_2}(s,t,u)= & {} G_{VP_1} G_{VP_2} \dfrac{-4t-2m_1^2-2m_2^2+s}{2\left( m_V^2-t\right) }\nonumber \\&+(t\leftrightarrow u)\nonumber \\ B_V^{P_1P_2}(s,t,u)= & {} G_{VP_1} G_{VP_2} \dfrac{1}{4\left( m_V^2-t\right) } +(t\leftrightarrow u)\ . \end{aligned}$$
(40)

They give rise to poles in the zero-width approximation, which is adequate in the kinematical regions of interest here. The helicity amplitudes are

$$\begin{aligned}&L_{V,++}^{P_1P_2}(s,\theta )=C_{VP_1}C_{VP_2}\dfrac{-s t}{m_V^2-t} +(t\leftrightarrow u)\nonumber \\&L_{V,+-}^{P_1P_2}(s,\theta )=C_{VP_1}C_{VP_2} \dfrac{\sin ^2\theta \lambda _{12}(s)}{4(m_V^2-t)} +(t\leftrightarrow u) \end{aligned}$$
(41)

and the corresponding partial-waves with \(J=0,2\) have the following form

(42)

with

$$\begin{aligned} X_V(s)=\frac{s-m_1^2-m_2^2+2m_V^2}{\sqrt{\lambda _{12}(s)}}\ \end{aligned}$$
(43)

which is the cosine of the scattering angle when \(t=m_V^2\). The function \(L_V(s)\) is given by the angular integral

$$\begin{aligned} L_V(s)=\int _{-1}^1 dz\, \frac{s+2m_V^2- m_1^2-m_2^2}{\lambda _{12}(s)(1-z^2) +4m_V^2(s-s_V)} \end{aligned}$$
(44)

with

$$\begin{aligned} s_V=-\frac{\left( m_V^2-m_1^2\right) \left( m_V^2-m_2^2\right) }{m_V^2}\ , \end{aligned}$$
(45)

it can be expressed in terms of \(X_V\) as

$$\begin{aligned} L_V(s)=\frac{\log (X_V(s)+1)-\log (X_V(s)-1)}{\sqrt{\lambda _{12}(s)}}\ . \end{aligned}$$
(46)

We note that the partial-wave \(l^V_{0++}(s)\) has a soft pion Adler zero at \(s=s_A\) which can be approximated as

$$\begin{aligned} s_A= & {} m_\eta ^2+m_\pi ^2\left( 1 +\dfrac{m_\eta ^2}{m_V^2}\left( -\dfrac{2}{3} +\dfrac{4}{135}\dfrac{{m_{\pi }}^2m_\eta ^2}{m_V^4}\right. \right. \nonumber \\&\left. \left. -\dfrac{8}{8505}\dfrac{{m_{\pi }}^4m_\eta ^4}{m_V^8}+\cdots \right) \right) \ . \end{aligned}$$
(47)

From the integral representation (44) one sees that the function \(L_V\) has endpoint singularities at \(z=\pm 1\) when \(s=s_V\), which thus corresponds to a branch point of \(L_V(s)\). Another endpoint singularity occurs when \(s=\infty \). When \(s > s_V\), the denominator remains strictly positive. Therefore, \(L_V\) is an analytic function of s with a cut on the negative real axis: \( -\infty< s < s_V\). The discontinuity of \(L_V\) along the cut is easily determined

$$\begin{aligned} \mathrm{Im}\,[L_V(s+i\epsilon )]=-\frac{\pi }{\sqrt{\lambda _{12}(s)}}\theta (s_V-s) \end{aligned}$$
(48)

from which one deduces the left-cut discontinuities of the vector-exchange partial-waves

$$\begin{aligned} \dfrac{1}{\pi }\mathrm{Im}\,\left[ l^V_{0,++}(s)\right]= & {} 2C_{V\pi }C_{V\eta }\, \dfrac{s\,m_V^2}{\sqrt{\lambda _{12}(s)}}\, \theta (s_V-s) \nonumber \\ \dfrac{1}{\pi }\mathrm{Im}\,\left[ l^V_{2,++}(s)\right]= & {} C_{V\pi }C_{V\eta }\, \dfrac{s\, m_V^2}{\sqrt{\lambda _{12}(s)}} \left( 3X_V^2(s)-1\right) \nonumber \\&\times \theta (s_V-s)\nonumber \\ \dfrac{1}{\pi }\mathrm{Im}\,\left[ l^V_{2,+-}(s)\right]= & {} -\dfrac{\sqrt{6}}{8} C_{V\pi }C_{V\eta }\, \sqrt{\lambda _{12}(s)}\, {\left( 1-X_V^2(s)\right) ^2} \nonumber \\&\times \theta (s_V-s) \ . \end{aligned}$$
(49)

We will use these discontinuities in the Muskhelishvili-Omnès representations below. We must consider also the vector exchange contributions for the \(\gamma \gamma \rightarrow (K\bar{K})_{I=1}\) amplitudes, we denote the relevant combination of coupling constants as

$$\begin{aligned} \tilde{C}_{K^*}^{(1)}\equiv \frac{1}{\sqrt{2}}\left( -C^2_{K^*K^+}+ C^2_{K^*K^0}\right) \ . \end{aligned}$$
(50)

The imaginary parts of the \(J=0,2\) partial-wave amplitudes along the left-hand cut read

$$\begin{aligned} \dfrac{1}{\pi }\mathrm{Im}\,\left[ k^V_{0,++}(s)\right]= & {} 2 \tilde{C}_{K^*}^{(1)}\, \dfrac{s\, m_{K^*}^2}{\sqrt{\lambda _{KK}(s)}} \theta (s_{K^*}-s) \nonumber \\ \dfrac{1}{\pi }\mathrm{Im}\,\left[ k^V_{2,++}(s)\right]= & {} \,\tilde{C}_{K^*}^{(1)} \dfrac{s\, m_{K^*}^2}{\sqrt{\lambda _{KK}(s)}} \left( 3X_{K^*}^2(s)-1\right) \nonumber \\&\times \theta (s_{K^*}-s) \nonumber \\ \dfrac{1}{\pi }\mathrm{Im}\,\left[ k^V_{2,+-}(s)\right]= & {} -\dfrac{\sqrt{6}}{8} \tilde{C}_{K^*}^{(1)}\, \sqrt{\lambda _{KK}(s)}\, {\left( 1-X_{K^*}^2(s)\right) ^2}\nonumber \\&\times \theta (s_{K^*}-s) \ . \end{aligned}$$
(51)

The updated values of the couplings \(C_{VP}\) are collected in Table 1 below. (a)

Table 1 Radiative widths of vector mesons and corresponding coupling constants. The relative signs of the couplings are determined assuming flavour symmetry

4 Representations of the \(J=0,2\) partial-waves

4.1 Muskhelishvili-Omnès representations for the S-waves

In order to write a dispersive representation for \(l_{0++}\) some knowledge concerning its asymptotic behaviour is needed. Let us then consider the angular integral,

$$\begin{aligned} l_{0++}(s)\equiv \int _0^1 dz\, L_{++}(s,t)\ \end{aligned}$$
(52)

in the \(s\rightarrow \infty \) limit. There are two regions of the angular variable z for which the behaviour of the integrand is known: (a) when z is close to 0, then \(t\sim u\sim -s\sim -\infty \). In this regime, \(L_{++}(s,t)\) can be estimated from QCD-based methods [40, 41] according to which one has

$$\begin{aligned} L_{++}(s,t) \ll L_{+-}(s,t) \sim \frac{\alpha _s(s)}{s}\ . \end{aligned}$$
(53)

(b)  When z is close to 1, then \(|t|<< s\) which is the region where Regge theory applies, the leading contribution from the vector meson trajectory gives

$$\begin{aligned} L_{++}(s,t) \sim \beta _V(t) (\alpha ' s)^{\alpha _V+\alpha ' t} \end{aligned}$$
(54)

with \(\alpha _V\simeq 0.5\), \(\alpha '\simeq 0.9\) \(\hbox {GeV}^{-1}\) and \(\beta _V(t)\) is a smooth function when \(t<0\). Assuming only that the integrand evolves smoothly between these two regimes when \(0\le z\le 1\), one deduces that the \(l_{0++}(s)\) should not grow faster than \(\sqrt{s}\) when \(s\rightarrow \infty \).

Furthermore, the \(J=0\) partial-waves obey soft-photon theorems [7, 8] which imply that the ratios \(l_{0++}(s)/s\) and \((k^1_{0++}(s)-k^{1,Born}_{0++}(s))/s\) remain finite when \(s\rightarrow 0\). Therefore, they can be expressed as unsubtracted dispersion relations in terms of the left-hand and right-hand cuts discontinuities,

$$\begin{aligned} l_{0++}(s)= & {} s\left[ \dfrac{1}{\pi }{\displaystyle \int _{-\infty }^{{{s}_V}} }ds'\, \dfrac{\mathrm{Im}\,[l_{0++}(s')]}{s' (s'-s)} \right. \nonumber \\&\left. +\dfrac{1}{\pi }{\displaystyle \int _{m_+^2}^\infty }ds'\, \dfrac{\mathrm{Im}\,[l_{0++}(s')]}{s' (s'-s)} \right] \nonumber \\ k^1_{0++}(s)= & {} k_{0++}^{1,Born}(s)+s\left[ \dfrac{1}{\pi }{\displaystyle \int _{-\infty }^{{s_{K^*}}} }ds'\, \dfrac{\mathrm{Im}\,\left[ k^1_{0++}(s')\right] }{s' (s'-s)} \right. \nonumber \\&\left. +\dfrac{1}{\pi }{\displaystyle \int _{m_+^2}^\infty }ds'\, \dfrac{\mathrm{Im}\,\left[ k^1_{0++}(s')\right] }{s' (s'-s)}\right] \end{aligned}$$
(55)

with \(m_+=m_\pi +m_\eta \). The isospin symmetry limit has been assumed here, implying that the \(\eta \) meson is stable and the discontinuities are real. At low energies the left-cut discontinuities are dominated by the vector meson exchanges. We will introduce subtractions, in order to reduce the influence of the higher energy regions where the discontinuities are not known. The right-cut discontinuities are given by the unitarity relations. Unitarity is known to be saturated with two channels to a very good approximation in the region of the \(a_0(980)\) resonance. We will assume that elastic unitarity holds below the \(K\bar{K}\) threshold and that two-channel unitarity remains a sufficiently good approximation up to \(\sqrt{s}=1.4\) GeV. Under the assumption of two-channel unitarity the dispersion relations (55) form a set of coupled inhomogeneous Muskhelishvili equations. They can be solved in terms of a two-channel MO matrix which satisfies a homogeneous set of coupled equations in terms of the T matrixFootnote 3

$$\begin{aligned} \varvec{\varOmega }_0(s)=\dfrac{1}{\pi }\int _{m_+^2}^\infty \frac{ds'}{s'-s} \varvec{T}(s)\varvec{\varSigma }(s)\varvec{\varOmega }^*_0(s)\ . \end{aligned}$$
(56)

Asymptotic conditions on the phase-shifts are imposed which ensure that the set of equations (56) has a unique solution once initial conditions at \(s=0\) are given (see Sect. 3.2 in Ref.  [27] and references therein). At \(s=0\) one can take

$$\begin{aligned} \varvec{\varOmega }_0(0)=\begin{pmatrix} 1 &{} 0\\ 0 &{} 1\\ \end{pmatrix}\ . \end{aligned}$$
(57)

Multiplying the amplitudes \((l_{0++},k^1_{0++})\) by the inverse of the matrix \(\varvec{\varOmega }_0\) removes the right-hand cuts. We can use this property in writing once-subtracted dispersion relations for the two functions

$$\begin{aligned} \begin{pmatrix} \phi _1(s)\\ \phi _2(s)\\ \end{pmatrix}\equiv \varvec{\varOmega }_0^{-1}(s) \begin{pmatrix} l_{0++}(s)/s\\ (k^1_{0++}(s)- k_{0++}^{1,Born}(s))/s\\ \end{pmatrix}\ . \end{aligned}$$
(58)

This provides MO-type dispersive representations for the \(\gamma \gamma \) amplitudes \(l_{0++}(s)\), \(k^1_{0++}(s)\) in terms of their imaginary parts on the left-cut and two parameters, \(b_l\), \(b_k\)

$$\begin{aligned} \begin{pmatrix} l_{0++}(s)\\ k^1_{0++}(s)\\ \end{pmatrix}= & {} \begin{pmatrix} 0\\ k_{0++}^{1,Born}(s)\\ \end{pmatrix} \nonumber \\&+s\,\varvec{\varOmega }_0(s) \begin{pmatrix} b_l + L_1(s)+R_1(s) \\ b_k + L_2(s)+R_2(s) \\ \end{pmatrix}\ . \end{aligned}$$
(59)

The functions \(L_i(s)\) are dispersive integrals over the left-hand cuts. We express them in a way which allows to easily implement the presence of an Adler zero at \(s=s_A\) in the amplitude \(l_{0++}\) (see “Appendix B”)

$$\begin{aligned}&L_i(s) {=}\dfrac{s-s_A}{\pi }\left\{ {\displaystyle \int _{-\infty }^{{s}_V}}\dfrac{ds'}{s'(s'-s_A)(s'-s)} \,D_{i1}(s')\,\mathrm{Im}\,\left[ l^V_{0++}(s')\right] \right. \nonumber \\&\qquad \left. + {\displaystyle \int _{-\infty }^{s_{K^*}}}\dfrac{ds'}{s'(s'-s_A)(s'-s)} \,D_{i2}(s')\,\mathrm{Im}\,\left[ k^V_{0++}(s')\right] \right\} \end{aligned}$$
(60)

where the functions \(D_{ij}(s)\) are the matrix elements of the inverse of the Omnès matrix,

$$\begin{aligned} \varvec{D}(s)\equiv \varvec{\varOmega }_0^{-1}(s)\ . \end{aligned}$$
(61)

The functions \(R_i(s)\), secondly, are the dispersive integrals over the right-hand cut,

$$\begin{aligned} R_i(s)= & {} -\dfrac{s-s_A}{\pi }\nonumber \\&\times {\displaystyle \int _{4m_K^2}^\infty } \dfrac{ds'}{s'(s'-s_A)(s'-s)}\, \nonumber \\&\quad \times \mathrm{Im}\,[D_{i2}(s')]\, k_{0++}^{1,Born}(s')\ . \end{aligned}$$
(62)

A relation between the parameters \(b_l\) and \(b_k\) can be derived from imposing that the amplitude \(l_{0++}\) has an Adler zero at \(s=s_A\),

$$\begin{aligned} b_l= -b_k \varOmega _{12}(s_A)/\varOmega _{11}(s_A)\ . \end{aligned}$$
(63)

The parameter \(b_k\) will eventually be fitted to the experimental data but we can estimate its order of magnitude by matching the amplitude \(k^1_{0++}\) with the SU(3) chiral expansion. Including the order \(p^4\) contributions (see “Appendix B”) and an estimate of the order \(p^6\) from the vector-meson exchange amplitudes, we obtain

$$\begin{aligned} b_k +L_2(0)+R_2(0)\simeq & {} -\dfrac{2\sqrt{2}}{F_\pi ^2}(L_9+L_{10})\nonumber \\&-\sqrt{2}\left( G^2_{K^{*0}K}-G^2_{K^{*+}K}\right) \dfrac{m_K^2}{m_{K^*}^2}\nonumber \\\simeq & {} -(0.57\pm 0.03)\ \hbox {GeV}^{-2} \end{aligned}$$
(64)

using the determination \(L_9+L_{10}= (1.44\pm 0.08)\cdot 10^{-3}\) taken from Ref. [42]. Below, the value of \(b_k\) will be fitted to the experimental data, the resulting combination \(b_k +L_2(0)+R_2(0)\) turns out to have a sign compatible with (64) and a magnitude smaller by a factor of two. The result of the dispersive construction of the S-wave amplitude \(l_{0++}\) is illustrated in Fig. 3. The corresponding result for the \(K\bar{K}\) amplitude \(k^1_{0++}\) is shown in “Appendix C” (see Fig. 13).

Fig. 3
figure 3

The S-wave amplitude \(l_{0++}\) from the dispersive construction using the central values of the fitted parameters

4.2 Dispersive construction of the isospin-violating S-wave

In Sect. 3.2 we have considered the unitarity contribution to the \(\gamma \gamma \rightarrow \pi \eta \) S-wave amplitude induced by the \({\pi ^+}{\pi ^-}\) intermediate state, which is isospin violating. Let us call \(\tilde{L}_{++}\) the isospin-violating part of the \(\gamma \gamma \rightarrow \pi \eta \) helicity amplitude and \(\tilde{l}_{0++}\) the corresponding S-wave. \(\tilde{L}_{++}\) can be defined as a matrix element,

$$\begin{aligned} \tilde{L}_{++}\equiv & {} \dfrac{1}{2}(m_d-m_u)\nonumber \\&\times {\langle {\pi ^0}(p_1)\eta (p_2)\vert (\bar{u}u{-}\bar{d}{d}) \vert \gamma (q_1,+)\gamma (q_2,+) \rangle }\ . \end{aligned}$$
(65)

We will attempt here to estimate the amplitude \(\tilde{l}_{0++}\) at low energy only and we write a dispersive representation keeping only the contribution from the \(\pi \pi \) cut,

(66)

We have used two subtractions in Eq. (66) in order to strongly reduce the influence of the energy region above 1 GeV. One subtraction constant is fixed from imposing the soft photon zero. We can then estimate the parameter \(\tilde{\lambda }\) in Eq. (66) by using the soft pion limit. This limit provides a relation between the amplitude \(\tilde{l}_{0++}(s=s_A)\) and a matrix element of the pseudo-scalar operator \(p_0=i(\bar{u}\gamma ^5u+\bar{d}\gamma ^5d)\) which, in turn, can be estimated using ChPT. This is detailed in “Appendix B.3”. Using Eqs. (B.22),  (B.23) from this appendix we obtain the following result for \(\tilde{\lambda }\)

$$\begin{aligned} \tilde{\lambda }= \dfrac{3\epsilon _L}{8\pi ^2 s_A}\,\left( 1-\dfrac{m_\pi ^2}{3F_\pi ^2}\right) \,\left( G_\pi (s_A)-\frac{1}{2} G_K(s_A)\right) \ \end{aligned}$$
(67)

where the loop functions \(G_\pi \), \(G_K\) are given in Eq. (B.1). The discontinuity \(\mathrm{disc}\,[{l}_{0++}(s)]_{\pi \pi }\) (given in Eq. (33)) has a singularity at the pseudo-threshold \(s=(m_\eta -{m_{\pi }})^2\) induced by the \(\pi \pi \rightarrow \eta \pi \) partial-wave (see Fig. 2). The integral in Eq. (66), however, is finite. It is defined by using the \(m_\eta ^2+i\epsilon \) limiting prescription exactly in the same way as those which appear in the Khuri-Treiman equations for the \(\eta \rightarrow 3\pi \) amplitude (see e.g. [36, 43]).

The result, in the low energy region relevant for \(\eta \rightarrow \pi \gamma \gamma \) is shown in Fig. 4 and compared to the corresponding chiral \(O(p^4)\) result from Eq. (B.6). The chiral and dispersive real parts agree at \(s=m_\eta ^2\), as a result of the soft pion relation, but the two amplitudes differ substantially at lower energy:

  • The cusp at the \(\pi \pi \) threshold is much more pronounced in the dispersive amplitude which is approximately five times larger in magnitude than the \(p^4\) amplitude at \(s=4m_\pi ^2\).

  • The \(p^4\) amplitude has a zero at \(s=4/3m_\pi ^2\) in contrast to the dispersive amplitude which has no zero in this region. This is because this zero is unrelated to a soft photon constraint and is an accidental feature of the \(p^4\) amplitude.

Fig. 4
figure 4

Dispersive calculation of the isospin-violating component of the \(\eta \rightarrow \pi \gamma \gamma \) S-wave amplitude (divided by s), compared to the chiral result at order \(p^4\). The shaded area shows the physical region

4.3 D-wave amplitudes modelling

For the \(J=2\) partial-waves one can write unsubtracted dispersion representations analogous to those for \(J=0\), e.g.,

$$\begin{aligned} l_{2\lambda \lambda '}(s)= & {} s^{|\lambda +\lambda '|/2}\lambda _{12}(s)\nonumber \\&\times \Bigg [\dfrac{1}{\pi }{\displaystyle \int _{-\infty }^{{{s}_V}}} ds'\, \dfrac{\mathrm{Im}\,[l_{2\lambda \lambda '}(s')]}{(s')^{|\lambda +\lambda '|/2}\lambda _{12}(s')(s'-s)}\nonumber \\&+\dfrac{1}{\pi }{\displaystyle \int _{m_+^2}^\infty } ds'\, \dfrac{\mathrm{Im}\,[l_{2\lambda \lambda '}(s')]}{(s')^{|\lambda +\lambda '|/2}\lambda _{12}(s')(s'-s)} \Bigg ] \end{aligned}$$
(68)

displaying the kinematical zeros at \(s=m_\pm ^2\) and \(s=0\). In the case of \(J=2\), however, it seems difficult to derive useful constraints from unitarity as in the case of \(J=0\) because in the important energy region of the \(a_2(1320)\) resonance there are too many contributing channels. We will therefore content with a simple Breit–Wigner estimate of the right-hand cut integral in Eq. (68).

We parametrise the coupling of the \(a_2\) resonance to the \(\eta \pi \) channel by a constant \(C^{a_2}_{\eta \pi }\) defined from the Lagrangian

$$\begin{aligned} {{\mathcal {L}}}_{a_2\eta \pi }= C^{a_2}_{\eta \pi }\, T_{\mu \nu }(x) \partial ^\mu \eta (x) \partial ^\nu \pi ^0(x)\ . \end{aligned}$$
(69)

The couplings to the \(\gamma \gamma \) channel involve two constants

(70)

The first term in (70) contributes to the \(+-\) helicity state, which is expected to be dominant, and the second term to the \(++\) helicity state. The expressions for the decay widths read,

$$\begin{aligned} \varGamma [a_2 \rightarrow \eta \pi ]= & {} \dfrac{ \left( C^{a_2}_{\eta \pi }\right) ^2}{60\pi } \dfrac{ \left( q_{\eta \pi }(m_T^2)\right) ^5}{m^2_T} \nonumber \\ \varGamma [a_2 \rightarrow \gamma \gamma ]= & {} \dfrac{e^4 m_T^3}{80\pi }\left( \left( C^{a_2}_{\gamma \gamma }\right) ^2 +\dfrac{1}{6}(D^{a_2}_{\gamma \gamma })^2 \right) \ . \end{aligned}$$
(71)

where \(q_{\eta \pi }(s)=\sqrt{\lambda _{\eta \pi }(s)/4s}\). The experimental values of the branching fractions of the main hadronic decay modes are [5]

$$\begin{aligned} B_{\eta \pi }= & {} (14.5\pm 1.2)\%,\quad B_{K\bar{K}}=(4.9\pm 0.8)\%,\nonumber \\ B_{3\pi }= & {} (70.1\pm 2.7)\%,\quad B_{\omega \pi \pi }=(10.6\pm 3.2)\% \end{aligned}$$
(72)

and the \(2\gamma \) width given by the PDG is

$$\begin{aligned} \varGamma ^{a_2}_{\gamma \gamma }= 1.00\pm 0.06\ \hbox {keV}\ . \end{aligned}$$
(73)

From this, one obtains the following values for the coupling constants

$$\begin{aligned} C^{a_2}_{\eta \pi }= & {} (10.8\pm 0.5)\ \hbox {GeV}^{-1}\nonumber \\&\times \sqrt{\left( C^{a_2}_{\gamma \gamma }\right) ^2 +\dfrac{1}{6}\left( D^{a_2}_{\gamma \gamma }\right) ^2}= (0.115\pm 0.005)\ \hbox {GeV}^{-1}\nonumber \\ \end{aligned}$$
(74)

choosing \(C^{a_2}_{\eta \pi }\) to have a positive sign. Upon performing fits to the differential cross-sections (see below) the two couplings \(C^{a_2}_{\gamma \gamma }\), \(D^{a_2}_{\gamma \gamma }\) will get separately determined. Defining a coupling constant \(C^{a_2}_{K\bar{K}}\) in the same way as \(C^{a_2}_{\eta \pi }\) we get,

$$\begin{aligned} C^{a_2}_{KK}=-(10.5\pm 0.9)\ \hbox {GeV}\ \end{aligned}$$
(75)

where the negative sign derives from flavour symmetry. The Breit–Wigner model for the \(a_2(1320)\) contributions is then taken as

$$\begin{aligned} l_{2++}^{BW}(s')= & {} \dfrac{D^{a_2}_{\gamma \gamma } C^{a_2}_{\eta \pi }}{60m_T^2} \sqrt{\dfrac{W_2(q_{\eta \pi }(m_T^2)R)}{W_2(q_{\eta \pi }(s')R)}}\nonumber \\&\times \dfrac{s'\lambda _{\eta \pi }(s')}{m_T^2-s'-im_T \varGamma _T(s)}\nonumber \\ l_{2+-}^{BW}(s')= & {} \dfrac{\sqrt{6}\,C^{a_2}_{\gamma \gamma } C^{a_2}_{\eta \pi }}{60} \sqrt{\dfrac{W_2(q_{\eta \pi }(m_T^2)R)}{W_2(q_{\eta \pi }(s')R)}}\nonumber \\&\times \dfrac{\lambda _{\eta \pi }(s')}{m_T^2-s'-im_T \varGamma _T(s)}\ . \end{aligned}$$
(76)

This form is obtained by first computing the amplitudes from the Lagrangians (69) and (70) and then modifying them by introducing a width function \(\varGamma _T(s)\) in the denominator, for which we use the same modellingFootnote 4 as in Ref.  [13]

$$\begin{aligned} \varGamma _T(s)=\varGamma _{\eta \pi }(s)+\varGamma _{K\bar{K}}(s) +\varGamma _{\rho \pi }(s)+\varGamma _{\omega 2\pi }(s) \end{aligned}$$
(77)

and a Blatt-Weisskopf related function (with \(W_2(x)=9+3x^2+x^4\)). Finally, the \(J=2\) amplitudes \(l_{2\lambda \lambda '}\) are approximated by adding the Breit-Wigner \(a_2(1320)\) contribution in the s-channel to the vector-exchange contributions in the t, u channels,

$$\begin{aligned} l_{2\lambda \lambda '}(s)= l^V_{2\lambda \lambda '}(s) + l^{BW}_{2\lambda \lambda '}(s)\ . \end{aligned}$$
(78)

The corresponding \(J=2\) \((K\bar{K})_{I=1}\) amplitudes are similarly described by a sum of three terms,

$$\begin{aligned} k^1_{2\lambda \lambda '}(s)= -\frac{1}{\sqrt{2}}k^{Born}_{2\lambda \lambda '}(s) +k^{K^*}_{2\lambda \lambda '}(s) + k^{BW}_{2\lambda \lambda '}(s)\ \end{aligned}$$
(79)

where the Breit–Wigner amplitudes are given by

$$\begin{aligned} k_{2++}^{BW}(s')= & {} \dfrac{D^{a_2}_{\gamma \gamma } C^{a_2}_{KK}}{60m_T^2} \sqrt{\dfrac{W_2(q_{KK}(m_T^2)R)}{W_2(q_{KK}(s')R)}}\nonumber \\&\times \dfrac{s\lambda _{KK}(s')}{m_T^2-s'-im_T \varGamma _T(s)}\nonumber \\ k_{2+-}^{BW}(s')= & {} \dfrac{\sqrt{6}\,C^{a_2}_{\gamma \gamma } C^{a_2}_{KK}}{60} \sqrt{\dfrac{W_2(q_{KK}(m_T^2)R)}{W_2(q_{KK}(s')R)}}\nonumber \\&\times \dfrac{\lambda _{KK}(s')}{m_T^2-s'-im_T \varGamma _T(s)}\ . \end{aligned}$$
(80)

It will be necessary to consider also the \((K\bar{K})_{I=0}\) amplitudes in order to be able to construct the \({K^+}{K^-}\) and \({K^0}{\bar{K}^0}\) amplitudes separately and compare with the experimental results. From the value of the \(K\bar{K}\) branching fraction of the \(f_2(1270)\) resonance [5]:

$$\begin{aligned} B^{f_2}_{K\bar{K}}=(4.6\pm 0.5)\%\ , \end{aligned}$$
(81)

one derives the following values for the corresponding coupling constants

$$\begin{aligned} C^{f_2}_{K\bar{K}}= & {} -(15.9\pm 0.9) \ \hbox {GeV}^{-1}\nonumber \\&\times \sqrt{ \left( C^{f_2}_{2\gamma }\right) ^2 {+} \frac{1}{6}\left( D^{f_2}_{2\gamma }\right) ^2}=(0.19\pm 0.02) \ \hbox {GeV}^{-1}\ .\nonumber \\ \end{aligned}$$
(82)

Based on nonet symmetry we have taken \(C^{f_2}_{K\bar{K}}\) and \(C^{a_2}_{K\bar{K}}\) to have the same sign.

5 Comparison with experiment

5.1 Experimental inputs

We will compare our model for the \(\gamma \gamma \) amplitudes with precise experimental data on \(\gamma \gamma \rightarrow \pi \eta \) from the Belle collaboration [13], as was done recently in Ref.  [18]. In addition, we consider also here \(\gamma \gamma \rightarrow K\bar{K}\) data in order to provide further constraints on the coupled-channel dynamics which is believed to be important for the \(a_0(980)\) resonance. Recently, high statistics experimental data have been obtained by the Belle collaboration for the \(K_S K_S\) channel [16]. Experimental data for the charged kaons channel \(K^+K^-\) are also available in this low energy range [44] but they are older and have much less statistics. We will restrict ourselves to the energy range \(E\le 1.4\) GeV: we can use 448 differential cross-section points for \(\pi \eta \) (\(0.85 \le E\le 1.39\) GeV) and 240 differential cross-section points for \(K_SK_S\) (\(1.105\le E\le 1.395\) GeV).

Table 2 Parameters of the two-channel T-matrix in the two fits
Fig. 5
figure 5

The phase-shifts and the inelasticity from the two-channel T-matrix model using the two sets of parameters corresponding to the two \(\chi ^2\) minimums

5.2 Parameters of the T-matrix

Concerning the S-wave, firstly, we employ the T-matrix model of Ref.  [27] which involves six parameters. It uses a chiral K-matrix type representation, which ensures one-channel unitarity below the \(K\bar{K}\) threshold and two-channel unitarity above, together with the chiral expansion. This kind of approach was initiated in [45], see Ref.  [46] for a review. The T-matrix is written as

$$\begin{aligned} \varvec{T}(s)=\left( \varvec{1}- \varvec{K}(s)\varvec{\varPhi }(s) \right) ^{-1} \varvec{K}(s)\ \end{aligned}$$
(83)

where the matrix \(\varvec{\varPhi }\) reads

$$\begin{aligned} \varvec{\varPhi }(s)=\begin{pmatrix} \alpha _1 +\beta _1 s +16\pi \bar{J}_{\eta \pi }(s) &{} 0 \\ 0 &{} \alpha _2 +\beta _2 s +16\pi \bar{J}_{KK}(s) \end{pmatrix}\nonumber \\ \end{aligned}$$
(84)

it involves four phenomenological polynomial parameters \(\alpha _i\), \(\beta _i\) and the one-loop functions \(\bar{J}_{P_1P_2}\) are given in (D.8). The K-matrix has the following form

$$\begin{aligned} \varvec{K}(s)= \varvec{K}_{(2)}(s) + \varvec{K}_{(4)}(s) + \varvec{K}_{(6)}(s) \end{aligned}$$
(85)

in which the subscript refers to the chiral order. The first two terms are to be computed from the chiral expansion of the scattering amplitudes involving the \(\eta \pi \) and \(K\bar{K}\) channels at order \(p^2\) and \(p^4\) respectively. The last term in Eq. (85) allows for a pole in s and involves two phenomenological parameters \(m_8\) and \(\lambda \),

$$\begin{aligned} \left[ \varvec{K}_{(6)}(s)\right] _{ij}=\lambda \frac{g_i g_j}{16\pi }\left( \frac{1}{m_8^2-s} - \frac{1}{m_8^2} \right) \ , \end{aligned}$$
(86)

The form of \(g_1\), \(g_2\) is derived from a resonance chiral Lagrangian [47]

$$\begin{aligned} g_1= & {} \dfrac{\sqrt{6}}{3F_\pi ^2}\left( c'_d\left( s-m_\eta ^2-m_\pi ^2\right) +2c'_m m_\pi ^2\right) \nonumber \\ g_2= & {} \dfrac{1}{F_\pi ^2}\left( c'_d\left( s-2m_K^2\right) +2c'_mm_K^2\right) \ \end{aligned}$$
(87)

such that \(\varvec{K}_{(6)}\) has chiral order \(p^6\) provided \(\lambda \) is O(1). In this model, the T-matrix has good analyticity properties and coincides with the chiral expansion at low energy up to \(O(p^4)\) provided that the parameters \(\alpha _i\), \(\beta _i\), \(\lambda \) are of chiral order O(1). In addition to these six phenomenological parameters, the T-matrix depends on the values of the \(O(p^4)\) chiral parameters \(L_i\) [48] and on the ratio \(c'_m/c'_d\). We will use here the set of \(L_i\) values from the \(p^6\) fit of Ref.  [49] (labelled as BE14 in that reference). The ratio \(c'_m/c'_d\) is expected to be of order \(1-2\), we will use \(c'_m/c'_d=2\) as central value and include the variation as a source of error.

Fig. 6
figure 6

Experimental \(\gamma \gamma \rightarrow \pi \eta \) differential cross-sections compared with the two fits results

Obviously, such a model which implements two-channel unitarity is mostly justified in the \(a_0(980)\) region and below. We will assume that it remains qualitatively acceptable up to \(E\simeq 1.4\) GeV. The T-matrix is computed from Eq. (83) for \(E \le E_1\), \(E_1=1.5\) GeV. In the higher energy region \(E>E_1\), the T-matrix is described through a simple interpolation of the phase-shifts and the inelasticity such that \(\delta _{11}(\infty )=2\pi \), \(\delta _{22}(\infty )=0\) and \(\eta (\infty )=1\). These conditions introduce a smooth cutoff in the integral equations satisfied by the matrix elements of the MO matrix and ensure the existence of a unique solution.

5.3 Fits results

In addition to the six S-wave T-matrix parameters listed above, further parameters must be introduced which describe couplings to the \(\gamma \gamma \) channel. In the S-wave, two parameters \(b_l\), \(b_k\) were introduced as subtraction constants. Implementing the Adler zero condition we keep only \(b_k\) as an independent parameter. In the D-wave sector, we include the values of the tensor resonance couplings \(C^{a_2}_{2\gamma }\), \(D^{a_2}_{2\gamma }\), \(C^{f_2}_{2\gamma }\), \(D^{f_2}_{2\gamma }\) in the fitting as well as the mass and width of the \(a_2(1320)\) resonance: \(m_{a_2}\), \(\varGamma _{a_2}\). In total, we thus have 6+7 parameters to be fitted.

At first, we have kept the T-matrix parameters fixed to one of the sets of values determined previously in Ref.  [27] (which, in particular, use assumed values for the pole positions of the two \(a_0\) resonances). It was not possible to obtain a good fit of the \(\gamma \gamma \) data in this manner: using these sets of parameter values one finds that the \(\pi \eta \) cross-section at the \(a_0(980)\) peak tends to be too large and the energy of the peak tends to be somewhat displaced as compared to experiment. Relaxing the T-matrix parameters, reasonably good fits become possible and we actually found two distinct minimums of the total \(\chi ^2\) combining the \(\pi \eta \) and the \(K_SK_S\) data,

$$\begin{aligned} \chi ^2\vert _{tot}= & {} (119+76)\vert _{\pi \eta } +233\vert _{K_SK_S} \qquad (\hbox {fit I})\nonumber \\= & {} (93+117)\vert _{\pi \eta } +229\vert _{K_SK_S} \qquad (\hbox {fit II})\ . \end{aligned}$$
(88)

In the case of \(\pi \eta \) the first number corresponds to the region \(E<1.1\) GeV. In this region fit II is better than fit I while the overall \(\chi ^2\) is slightly smaller in fit I (\(\chi ^2=428\)) than in fit II (\(\chi ^2=439\)). The \(\chi ^2\) is defined in a simple and naive way: the correlation matrix is assumed to be diagonal and the statistical and systematic errors provided by the Belle collaboration are added in quadrature. The searches for minimums were performed with the help of the computer code MINUIT [50].

The numerical values of the set of T-matrix parameters resulting from these two fits are shown in Table 2. One notices, in particular, that the value of the pole parameter \(m_8\) differs significantly in the two fits. Figure 5 shows the behaviour of the two phase-shifts \(\delta _{11}\), \(\delta _{22}\) and of the inelasticity parameters \(\eta \) as a function of energy, corresponding to the two fits. At low energies, the \(\pi \eta \) phase-shift changes sign and becomes negative. This low-energy behaviour is not anticipated in simple hadronic models of the \(\pi \eta \) amplitude (e.g. [51]) but it was also observed to emerge from fitting the \(\gamma \gamma \) data in Ref.  [25].

In fit II the \(\pi \eta \) scattering length (defined as in [52]) is found to have the following value

$$\begin{aligned} {m_{\pi }}a_0^{\pi \eta }= -\left( 9.5^{+0.3}_{-6.3}\right) \times 10^{-3}\ . \end{aligned}$$
(89)

For comparison, the first two terms in the chiral expansion of the scattering length give

$$\begin{aligned} {m_{\pi }}a_0^{\pi \eta }= & {} \left. 6.17\times 10^{-3}\right| _{p^2},\nonumber \\= & {} \left. 2.43\times 10^{-3}\right| _{p^2+p^4}. \end{aligned}$$
(90)

using the same set of chiral parameters as in the unitary T-matrix.

Fig. 7
figure 7

Cross-sections for \(\gamma \gamma \rightarrow \pi \eta , K_SK_S, {K^+}{K^-}\) integrated in the range \(|\cos \theta |<0.8\). The data are from Refs. [13, 16, 44], they are compared with the two fits results

At higher energies a clear difference between the two fits is the sharp increase of the \(\pi \eta \) phase-shift in fit I around \(E\simeq 1.32\) GeV, typical of a narrow resonance, which we will call \(a'_0\). Indeed, one finds a resonance pole in the T-matrix located on the third Riemann sheet with the value,Footnote 5

$$\begin{aligned} \sqrt{s_{a'_0}}=1315(4)-i\,24(3)\ \hbox {MeV}\qquad (\hbox {fit I}). \end{aligned}$$
(91)

This resonance is lighter and narrower than the standard \(a_0(1450)\). Since the phase-shift increases from \(\pi /2\) to \(\pi \) (approximately) it gives rise to a sharp dip (instead of a peak) in the S-wave cross-section. Clearly then, our fit I is quite analogous to the best fit by the Belle collaboration [13] which displays a resonance (called \(a_0(Y)\) in that reference) which has very similar features while using a parametrisation rather different from ours.Footnote 6 The S-wave amplitude in fit II also has an \(a'_0\) resonance pole but it is heavier and much broader than that of fit I,

$$\begin{aligned} \sqrt{s_{a'_0}}=1421(5)-i\,175(4)\ \hbox {MeV}\qquad (\hbox {fit II})\ . \end{aligned}$$
(92)

In this case, the width of the resonance is larger than the experimental one. One eventually does not expect a very accurate determination since the resonance lies close to the cutoff of the model. Finally, the effects of the two \(a_0\) resonances on the \(\pi \eta \) and \(K\bar{K}\) scalar form factors are shown in Fig. 11 in “Appendix A”.

Figures 6 and 8 show a sample of \(\gamma \gamma \rightarrow \pi \eta , K_SK_S\) differential cross sections comparing the experimental results with those from the two fits. The difference between the two fits is remarkably small, which is somewhat puzzling: why is the narrow \(a'_0\) resonance present in fit I not seen much more clearly? The reason for this arises from a specific interference effect between the \(J^{\lambda \lambda '}=0^{++}\) and the \(2^{++}\) amplitudes. We have seen already that the fast energy variation induced by the narrow \(a'_0\) coincides with that induced by the \(a_2(1320)\). More specifically, when \(\sqrt{s}\simeq 1.32\) GeV, the following relations hold approximately between the \(l_{0++}\) and the \(l_{2++}\) amplitudes in the case of fit I

$$\begin{aligned} \mathrm{Re}\,[l_{0++}(s)]\simeq & {} -5\mathrm{Re}\,[l_{2++}(s)]\nonumber \\ \mathrm{Im}\,[l_{0++}(s)]\simeq & {} -5\mathrm{Im}\,[l_{2++}(s)]+0.7\ . \end{aligned}$$
(93)

The sum of the \(0^{++}\) and the \(2^{++}\) partial-wave amplitudes can effectively be absorbed into the \(2^{+-}\) resonance amplitude (since the angular functions satisfy the relation \(d^2_{00}(\theta )-1 = -\sqrt{6}\, d^2_{20}(\theta )\)) and no specific \(0^{++}\) resonance effect remains. We also note that the \(a'_0\) in fit I essentially decouples from the \(K\bar{K}\) channel such that no \(0^{++}\) resonance effect is seen in this channel either. Because of these fine-tuned relations between the scalar and the tensor resonances the fit I solution is most likely to be unphysical.

One can see from Fig. 6 that the \(\pi \eta \) differential cross-sections are well described except, however, in one energy bin: \(E=1.01\) MeV. At this energy, the cross-sections are underestimated by our model, see also Fig. 7 which shows the cross-sections integrated over \(\theta \). This discrepancy could be caused by an isospin breaking effect close to the resonance peak. Our model assumes isospin symmetry and uses \(m_K\equiv m_{K^+}\) (in order to correctly implement the kaon Born amplitudes) and the peak of the cross-section occurs exactly at \(E=2m_{K^+}\). Physically, however, the \(K^+\) and the \(K^0\) have slightly different masses which should lead to two cusps in the shape of the integrated cross-section. On the experimental side, the energy resolution should be smaller than \(2m_{K^+}-2m_{K^0}\simeq 8\) MeV in order to clearly observe this effect.

Fig. 8
figure 8

Experimental \(\gamma \gamma \rightarrow K_SK_S\) differential cross-sections compared with the two fits results

The values of the remaining seven parameters included in the fit (subtraction parameter \(b_k\), tensor mesons to \(2\gamma \) coupling constants, \(a_2\) mass and width) are collected in Table 3 below. The couplings \(D^T_{2\gamma }\) are found, as expected, to be smaller in magnitude than the \(C^T_{2\gamma }\) although this suppression is only by a factor of two in the case of fit I. The values are in qualitative agreement with the PDG expectations except, however, for the \(a_2\) mass which is shiftedFootnote 7 by approximately 10 MeV. This shift is easily seen to be caused by the presence of non-resonant contributions to the \(J=2\) amplitudes modelled here by the vector-meson exchanges, \(l^V_{2,\lambda \lambda '}\). The presence of this term is essential for obtaining a correct description of the amplitude in the \(\eta \) decay region \(s < (m_\eta -{m_{\pi }})^2\). In the 1 GeV energy region, one expects some modifications induced by higher mass exchanges, but we see no reason that it should be completely cancelled. We will consider a variation of this term as a source of error below.

Table 3 Values of the coupling constants, of the subtraction constant \(b_k\) and of the mass and width of the \(a_2(1320)\) resonance resulting from the two fits

5.4 Properties of the \(a_0\) resonances

In this section we consider in more detail the properties of the two \(a_0\) resonances which can be deduced from our analysis of the photon–photon data. We will focus on the results from fit II since fit I was argued not to be physically relevant. The formulas needed to define the T-matrix elements on the unphysical Riemann sheets are given in “Appendix D”. One may define coupling constants from the residues of the resonance poles (e.g. [54] in the context of photon–photon amplitudes), the couplings to the \(\pi \eta \) and \(K\bar{K}\) channels are thus defined as

$$\begin{aligned} \left. 16\pi T_{11}^{(II)}(z)\right| _{pole}= & {} \dfrac{g^2_{a_0\pi \eta }}{z_{a_0}-z}\nonumber \\ \left. 16\pi T_{12}^{(II)}(z)\right| _{pole}= & {} \dfrac{g_{a_0\pi \eta }g_{a_0K\bar{K}}}{z_{a_0}-z} \end{aligned}$$
(94)

and similarly for the third Riemann sheet. The coupling to the \(\gamma \gamma \) channel and the associated width are defined as

$$\begin{aligned} \left. e^2l_{0++}^{(II)}(z)\right| _{pole}= \frac{g_{a_0\gamma \gamma }g_{a_0\pi \eta }}{z_{a_0}-z},\quad \varGamma _{a_0\rightarrow \gamma \gamma }= \frac{\vert g_{a_0\gamma \gamma } \vert ^2}{16\pi m_{a_0}}\ . \end{aligned}$$
(95)

The following numerical results are found on sheet II for the position of the \(a_0(980)\) pole and the corresponding coupling constants

$$\begin{aligned} \sqrt{s_{a_0}}= & {} 1000.7(7)-i\,36.6(1.3) \ (\hbox {MeV})\nonumber \\ |g_{a_0\pi \eta }|= & {} 2.17(2) \ (\hbox {GeV})\nonumber \\ |g_{a_0K\bar{K}}|= & {} 4.03(2)\ (\hbox {GeV})\nonumber \\ \varGamma _{a_0\gamma \gamma }= & {} 0.52(1) \ (\hbox {keV})\ . \end{aligned}$$
(96)

From the sheet III resonance, the position of the pole was given in Eq. (92) and the corresponding coupling constants have the following values

$$\begin{aligned} |g_{a'_0\pi \eta }|= & {} 3.15(4) \ (\hbox {GeV})\nonumber \\ |g_{a'_0K\bar{K}} |= & {} 1.89(4) \ (\hbox {GeV})\nonumber \\ \varGamma _{a'_0\gamma \gamma }= & {} 1.05(5) \ (\hbox {keV})\ . \end{aligned}$$
(97)
Table 4 Errors generated by varying the fixed parameters of the S-wave amplitudes, columns 5 and 6 refer to the left-cut (see text)

The errors quoted in the above formulas are those arising from the fitted parameters (for which MINUIT provides a correlation matrix). In addition to those, one must consider errors associated with further parameters on which the T-matrix depends, which we have assumed to be fixed up to now:

  1. 1.

    Adler zero: the value of Adler zero was fixed to \(s_A=m_\eta ^2\) but its exact value is not known, we will vary it in the range \(s_A=m_\eta ^2\pm 3m_\pi ^2\).

  2. 2.

    Set of the \(O(p^4)\) parameters \(L_i\): in the BE14 determination [49], their values are given with errors but there are strong correlations. In order to estimate an error from this source we used two different sets of \(L_i's\) which correspond to two best fits made with different assumptions (second and third column in Table 3 of Ref.  [49]).

  3. 3.

    The ratio of the scalar couplings \(c'_m/c'_d\) (see Eqs. (86), (87)) was varied between 1 and 3 to estimate an uncertainty from this source.

  4. 4.

    Left-hand cut: (a) we added a contribution from a set of axial-vector mesonsFootnote 8 and (b) we varied the products of the couplings of the vector mesons by \(\pm 20\%\).

Table 4 gives the detailed list of the errors generated from these variations on the properties of the two \(a_0\) resonances.

5.5 The \(\eta \rightarrow {\pi ^0}\gamma \gamma \) decay amplitude

We consider now the modelled amplitudes in the kinematical region relevant for the decay \(\eta \rightarrow {\pi ^0}\gamma \gamma \). Including the \(J=0\) and \(J=2\) partial-waves,Footnote 9 the distribution as a function of the \(\gamma \gamma \) invariant mass squared, s, reads (Fig. 8)

$$\begin{aligned} \dfrac{d\varGamma ^{\eta \rightarrow \gamma \gamma \pi } }{ds}= & {} \dfrac{\alpha ^2}{16\pi m_\eta ^3} \sqrt{\lambda _{12}(s)} \Bigg \{ \vert l_{0,++}(s)+\tilde{l}_{0,++}(s)\vert ^2\nonumber \\&+5\,\left| l^V_{2,++}(s) + l^{BW}_{2,++}(s)\right| ^2 \ \nonumber \\&+5\,\left| l^V_{2,+-}(s) + l^{BW}_{2,+-}(s)\right| ^2 \Bigg \} \end{aligned}$$
(98)

also accounting for the \(J=0\) isospin violating amplitude \(\tilde{l}_{0,++}\).

Fig. 9
figure 9

Contributions to the \(\eta \rightarrow \pi ^0\gamma \gamma \) energy distribution from the amplitudes fitted in the scattering region. Dotted line (red): S-wave (isospin conserving) and \(2^{++}\) D-wave, dash-dotted line (magenta): \(2^{+-}\) D-wave, dashed line: sum of the S and D waves, solid line: sum including the isospin violating amplitude \(\tilde{l}_{0++}\)

Fig. 10
figure 10

Experimental data on the \(\eta \rightarrow \pi ^0\gamma \gamma \) energy distribution [14, 15] compared with predictions from our amplitudes showing the influence of the position of the Adler zero

Figure 9 illustrates the contributions of the various amplitudes corresponding to central values of the parameters fitted in the scattering region. The \(2^{+-}\) partial-wave dominates over the other ones near \(s=0\) because the \(J^{++}\) amplitudes are suppressed by the soft-photon zero. The relative role of the S-wave increases with the energy and starts to be dominating above the \(\pi ^+\pi ^-\) threshold. The figure also shows that the isospin-violating S-wave generates a visible cusp at this threshold. The central value of the decay width generated by our amplitudes is \(\varGamma =0.237\) eV which is on the low side of the most recent experimental determinations \(\varGamma _{exp}=0.285\pm 0.031\pm 0.061\) eV (Crystal Ball at the AGS [14]), \(\varGamma _{exp}=0.33\pm 0.03\) eV (MAMI [15]). A smaller value was reported by the KLOE-2 collaboration [56] but it has not been confirmed and the data is currently being reanalysed. The amplitudes in the \(\eta \) decay region are very sensitive to the precise position of the Adler zero. This is illustrated in Fig. 10 showing the effect of varying \(s_A\) in the range \([m_\eta ^2-3m_\pi ^2,m_\eta ^2+3m_\pi ^2]\) and comparing with the experimental results. Given somewhat more precise data the value of \(s_A\) could be included in the fitting. The amplitudes in the decay region are also sensitive to the vector meson coupling constants, the value of which dominate the energy distribution near \(s=0\). Accounting for these main errors we would predict

$$\begin{aligned} \varGamma ^{\eta \rightarrow {\pi ^0}\gamma \gamma }= 0.237^{+0.060}_{-0.043}\ \hbox {eV}\ . \end{aligned}$$
(99)

6 Conclusions

In this work we have reconsidered the properties of the light isovector scalar resonances as can be determined from photon–photon scattering experimental results. For this purpose, we have implemented a standard Muskhelishvili–Omnès integral representation for the \(J=0\) amplitude in which the left-cut is modelled from light vector meson exchanges. The underlying T-matrix satisfies unitarity with two channels (\(\pi \eta ,K\bar{K}\)) and involves six phenomenological parameters. In the case of the \(J=2\) amplitudes the constraints from unitarity are more difficult to implement, a cruder description is used which consists in simply adding the cross-channel vector-exchange and the direct channel tensor resonance amplitudes. In order to constrain the free parameters as unambiguously as possible we performed fits to both \(\pi \eta \) and \(K_SK_S\) data, for which high-statistics data below 1.4 GeV are available at present, and we found two different acceptable solutions to the minimisation. Both solutions are also compatible with the available \(K^+K^-\) data from the ARGUS collaboration [44].

In one of our fits the S-wave amplitude displays a light and narrow \(a'_0\) resonance exactly similar to the one found in the Belle analysis. While this is mathematically allowed we have argued that the fit which displays a broad \(a'_0\) is likely to be more physical. Concerning the \(a_0(980)\) resonance, we find that a rather conventional picture i.e. a pole on the second sheet with a mass and width compatible with the PDG and coupling to both the \(\pi \eta \) and the \(K\bar{K}\) channels is perfectly compatible with both the \(\pi \eta \) and the \(K_SK_S\) data. This is in contrast with the Belle analysis which uses an elastic Breit-Wigner description and also with the recent analysis of Ref.  [18] in which the mass and width are found to be both significantly larger than the PDG values. Data with a better energy resolution would be useful to resolve these remaining ambiguities. The \(\gamma \gamma \rightarrow {K^+}{K^-},K_SK_S\) cross-sections close to the \(K\bar{K}\) thresholds are also very sensitive to the position of the \(a_0(980)\). Our results in this energy region are in qualitative agreement with the chiral-unitary calculations from Ref.  [29] and with the estimates made in Ref.  [57] but not with those from Ref.  [26]. Experimental data in this near-threshold region would obviously be very constraining. Finally, it will be quite interesting to see how the pole position determined in a lattice QCD simulation [22] evolves when the value of \(m_\pi \) is decreased.