\section{Introduction}
Polarizabilities of hadrons characterize their response to applied electromagnetic field which cannot be attributed to point-like particles, and are of fundamental interest for low-energy QCD. Experimental measurement of polarizabilities is challenging, so most data are available for lightest mesons $\pi^\pm,\pi^0$, but there is long-standing discrepancy in the values (see papers \cite{Ivanov:2015ntu,Moinester:2019sew,Moinester:2022tba} for review of theoretical and experimental status of the meson polarizability problem). Among the reported experimental results, only the most recent data on charged pion polarizability by COMPASS collaboration at CERN~\cite{COMPASS:2014eqi} is consistent with Chiral Perturbation Theory. The leading-order result of ChPT~\cite{Bijnens:1987dc} is equivalent to the value found in Ref.~\cite{Terentev:1972ix} based on the hypothesis of partially conserved axial-vector current.

The polarizabilities were investigated theoretically within Chiral Perturbation Theory up to two loops~\cite{Gasser:2006qa,Burgi:1996qi,Burgi:1996mm,Bellucci:1994eb,Gasser:2005ud}, within the methods of Lattice QCD~\cite{Fiebig:1988en,Lee:2005dq,Lujan:2014kia,Freeman:2014kka,Luschevskaya:2015cko,Lujan:2016ffj,Bali:2017ian,Niyazi:2021jrz,Bignell:2020dze,Ding:2020hxw,Wilcox:2021rtt} and various phenomenological models~\cite{Lvov:1980st,Volkov:1980cj,Volkov:1985nk,Bernard:1988gp,Ivanov:1991kw,Dorokhov:1997rv,Filkov:1982cx,Donoghue:1993kw}. 
Several studies~\cite{Ivanov:1991kw,Filkov:1998rwz,Dorokhov:1997rv} found that dominating part of pion polarizabilities is due to $\sigma$ meson, even though experimental status of $\sigma$ meson is  controversial.

The latter motivates investigation of polarizabilities within the Domain Model of QCD vacuum and hadronization~\cite{Kalloniatis:2003sa,Nedelko:2016gdk,Nedelko:2016vpj,Nedelko:2014sla,Nedelko:2020bba} in which meson fields appear as nonlocal collective colorless excitations of the quark, anti-quark  (heavy and light ones) and gluon fields in the presence of the confining vacuum fields that are represented by the statistical ensemble of the almost everywhere homogeneous Abelian (anti-)self-dual gluon fields. A highly nonlocal effective meson action contains information about the  strong, weak and electromagnetic interactions of mesons as well as their two-point correlation functions.    In particular, the model consistently describes various phenomena related to confinement and chiral symmetry realization, as well as specific limiting regimes like the heavy quark limit.  At high momentum transfer, quark and gluon propagators in the presence of the vacuum field reproduce parton picture, while they are strongly modified in the infrared regime.  Calculated within this approach meson masses, including Regge spectrum of excited states of mesons, their decay constants and form factors  are in good agreement with experimental values~\cite{Efimov:1995uz,Burdanov:1996uw,Kalloniatis:2003sa,Nedelko:2016gdk,Nedelko:2016vpj,Nedelko:2014sla,Nedelko:2020bba}. In the present study, polarizabilities are calculated in the same framework. Meson masses are identified as the poles of the nonlocal meson propagators. The specific feature of this approach particularly relevant to the present study is the absence of the poles at real momenta in the nonlocal propagators of the scalar meson-like composite fields. Light  scalar mesons as  quark-anti-quark states analogous to the pseudoscalar ones are absent in the physical spectrum of the stable collective excitations.  This property is due to  peculiarities of realization of chiral symmetry in the presence of Abelian (anti)-self-dual gluon mean field due to the quark quasi-zero modes. At the same time, virtual scalar meson-like fields contribute to the amplitudes of various processes.  In has to be noted that physical scalar states occur in the hyperfine splitting of the orbital excitations of the vector mesons with the masses above 1 GeV. 


In Section~\ref{section_polarizabilities} we calculate polarizabilities and find that charged pion polarizability is consistent with the results of ChPT~\cite{Gasser:2006qa} and COMPASS data~\cite{COMPASS:2014eqi}. The contribution of intermediate scalar meson fields to polarizability of pion and kaon is substantial,  especially for the neutral ones.  For charged pion and kaon polarizability  intermediate scalar mesons turn out to be less  important.



\section{Effective meson action\label{section_model}}


The effective meson action deduced via bosonization of the one-gluon exchange interaction of the quark currents in the presence of  the the random Abelian (anti-)self-dual vacuum gluon fields has the form~\cite{Efimov:1995uz,Kalloniatis:2003sa,Nedelko:2016gdk,Nedelko:2016vpj,Nedelko:2014sla,Nedelko:2020bba}:
\begin{align}
\nonumber
Z&={\cal N}
\int D\phi_{\cal Q}
\exp\left\{-\frac{\Lambda^2}{2}\frac{h^2_{\cal Q}}{g^2 C^2_\mathcal{Q}}\int d^4x 
\phi^2_{\cal Q}(x)
-\sum\limits_{k=2}^\infty\frac{1}{k}W_k[\phi]\right\},
\\
\label{effective_meson_action}
W_k[\phi]&=
\sum\limits_{{\cal Q}_1\dots{\cal Q}_k}h_{{\cal Q}_1}\dots h_{{\cal Q}_k}
\int d^4x_1\dots\int d^4x_k
\Phi_{{\cal Q}_1}(x_1)\dots \Phi_{{\cal Q}_k}(x_k)
\Gamma^{(k)}_{{\cal Q}_1\dots{\cal Q}_k}(x_1,\dots,x_k),
\\
\label{physical_meson_fields}
\Phi_{{\cal Q}}(x)&=\int \frac{d^4p}{(2\pi)^4}e^{ipx}{\mathcal O}_{{\mathcal Q}{\mathcal Q}'}(p)\tilde\phi_{{\mathcal Q}'}(p),\quad C_\mathcal{Q}=C_J,\quad C^2_{S/P}=2C^2_{V/A}=\frac{1}{9}.
\end{align}
The condensed index $\mathcal{Q}\equiv\{aJLn\}$ includes all quantum numbers of a meson, $\Lambda$ is a scale  related to the strength of the vacuum gluon field, and  finally  to the value of gluon condensate $\langle g^2 F^2 \rangle$.  The physical color neutral meson fields $\phi_{{\mathcal Q}}$ are obtained by means of  orthogonal transformation ${\mathcal O}_{{\mathcal Q}{\mathcal Q}'}$ of fields $\Phi_{{\mathcal Q}}$.   The quadratic part of the action for $\phi_{{\mathcal Q}}$ is diagonal with respect to all quantum numbers. The masses of mesons correspond to the poles of nonlocal propagators
\begin{equation*}
D_\mathcal{Q}(p^2)=h_\mathcal{Q}^{-2} \left(\frac{\Lambda^2}{g^2 C^2_\mathcal{Q}}+\tilde\Gamma^{(2)}_{\mathcal{Q}}(p^2)\right)^{-1}
\end{equation*}
and can be found as zeroes of inverse propagator from equation
\begin{align}
\label{mass-eq}
0&=
\frac{\Lambda^2}{g^2 C^2_\mathcal{Q}}+\tilde\Gamma^{(2)}_{\cal Q}(-M^2_{\cal Q}),
\end{align}
where $\tilde\Gamma^{(2)}_{\cal Q}$ is two-point correlation function diagonalized with respect to all quantum numbers. 
Constants $h_\mathcal{Q}$ are defined by equation
\begin{align}
\nonumber
1&=h^2_{\cal Q}
\left.\frac{d}{dp^2}\tilde\Gamma^{(2)}_{\cal Q}(p^2)\right|_{p^2=-M^2_{\cal Q}}
\end{align}
which ensure that residue at the pole of the propagator is equal to unity. The results of calculation of the masses of various mesons as well as analytical expressions of $\tilde\Gamma^{(2)}_{\cal Q}(p^2)$ can be found in Ref.~\cite{Nedelko:2016gdk}. 

Correlation functions include ``connected'' and ``disconnected'' contributions of quark loops in the background field. For example, the two-point nonlocal vertex function $\tilde\Gamma^{(2)}_{\cal QQ'}(p)$ is given by
\begin{equation}
\label{Gammak}
\Gamma^{(2)}_{{\cal Q}_1{\cal Q}_2}=
\overline{G^{(2)}_{{\cal Q}_1{\cal Q}_2}(x_1,x_2)}-
\Xi_2(x_1-x_2)\overline{G^{(1)}_{{\cal Q}_1}G^{(1)}_{{\cal Q}_2}},
\end{equation}
where $\Xi$ is correlation function of the background field which belongs to the statistical ensemble of the almost everywhere homogeneous Abelian (anti-)self-dual fields.
Quark loops $G^{(k)}_{{\cal Q}_1\dots{\cal Q}_k}$ are averaged over the ensemble of background field configurations with measure $d\sigma_B$:
\begin{eqnarray}
\label{barG}
&&\overline{G^{(k)}_{{\cal Q}_1\dots{\cal Q}_k}(x_1,\dots,x_k)}
=\int d \sigma_B
{\rm Tr}V_{{\cal Q}_1}\left(x_1\right)S\left(x_1,x_2\right)\dots
V_{{\cal Q}_k}\left(x_k\right)S\left(x_k,x_1\right),
\\
&&\overline{G^{(l)}_{{\cal Q}_1\dots{\cal Q}_l}(x_1,\dots,x_l)
G^{(k)}_{{\cal Q}_{l+1}\dots{\cal Q}_k}(x_{l+1},\dots,x_k)}
=
\nonumber\\
\nonumber
&&\int d \sigma_B
{\rm Tr}\left\{
V_{{\cal Q}_1}\left(x_1\right)S\left(x_1,x_2\right)\dots
V_{{\cal Q}_k}\left(x_l\right)S\left(x_l,x_1\right)
\right\}\times
\\
&&
{\rm Tr}\left\{
V_{{\cal Q}_{l+1}}\left(x_{l+1}\right)S\left(x_{l+1},x_{l+2}\right)\dots
V_{{\cal Q}_k}\left(x_k\right)S\left(x_k,x_{l+1}\right)
\right\}.
\nonumber
\end{eqnarray}
Here $S(x,y)$ is the quark propagator and $V_{{\cal Q}}$ are  nonlocal  meson-quark-antiquark vertices. 

The quark propagator and meson-quark vertices in the presence of the almost everywhere homogeneous fields are approximated by those in the homogeneous (anti-)self-dual Abelian background field. The averaging over mean field ensemble is achieved by averaging the quark loops over configurations of the homogeneous  background fields, supplemented by taking into account  $n$-point correlators of the mean fields $\Xi_n$. The averaging is performed over self-dual and anti-self-dual Abelian (anti-)self-dual configurations and their  directions in Euclidean and color spaces. Averaging over spatial directions in $R^4$ is performed with the help of generating formula
\begin{equation}
\label{averaging_over_vacuum_field}
\langle\exp(if_{\mu\nu}J_{\mu\nu})\rangle=\frac{\sin\sqrt{2\left(J_{\mu\nu}J_{\mu\nu}\pm J_{\mu\nu}\widetilde{J}_{\mu\nu}\right)}}{\sqrt{2\left(J_{\mu\nu}J_{\mu\nu}\pm J_{\mu\nu}\widetilde{J}_{\mu\nu}\right)}},
\end{equation}
where $J_{\mu\nu}$ is an arbitrary antisymmetric tensor.
Tensor $f_{\mu\nu}$ is an appropriately normalized Abelian (anti-)self-dual background field with strength $\Lambda$:
\begin{gather}
\nonumber
\hat B_\mu=-\frac{1}{2}\hat n B_{\mu\nu}x_\nu, \ \hat n = t^3\cos\xi+t^8\sin\xi,
\\
\label{b-field}
\tilde{B}_{\mu\nu}=\frac12\epsilon_{\mu\nu\alpha\beta}B_{\alpha\beta}=\pm B_{\mu\nu}, \  \hat{B}_{\rho\mu}\hat{B}_{\rho\nu}=4\upsilon^2\Lambda^4\delta_{\mu\nu},\\
\nonumber
f_{\alpha\beta}=\frac{\hat{n}}{2\upsilon\Lambda^2}B_{\alpha\beta}, \  \upsilon=\mathrm{diag}\left(\frac16,\frac16,\frac13\right), \ f_{\mu\alpha}f_{\nu\alpha}=\delta_{\mu\nu},
\end{gather}
where the upper sign in ``$\pm$'' should be taken for self-dual field, and the lower for anti-self-dual field.
Nonlocal vertices $V^{aJln}_{\mu_1\dots\mu_l}$ are given by formulas
\begin{gather}
V^{aJln}_{\mu_1\dots\mu_l}= {\cal C}_{ln}\mathcal{M}^a\Gamma^J F_{nl}\left(\frac{\stackrel{\leftrightarrow}{\cal D}^2\!\!\!
(x)}{\Lambda^2}\right)T^{(l)}_{\mu_1\dots\mu_l}\left(\frac{1}{i}\frac{\stackrel{\leftrightarrow}{\cal D}\!(x)}{\Lambda}\right),
\label{qmvert}\\
{\cal C}^2_{ln}=\frac{l+1}{2^ln!(n+l)!},\quad F_{nl}(s)=s^n\int_0^1 dt t^{n+l} \exp(st),
\nonumber\\
{\stackrel{\leftrightarrow}{\mathcal{D}}}\vphantom{D}^{ff'}_{\mu}=\xi_f\stackrel{\leftarrow}{\mathcal{D}}_{\mu}-\ \xi_{f'}\stackrel{\rightarrow}{\mathcal{D}}_{\mu}, 
\ \
\stackrel{\leftarrow}{\mathcal{D}}_{\mu}\hspace*{-0.3em}(x)=\stackrel{\leftarrow}{\partial}_\mu+\ i\hat B_\mu(x),  \ \ 
\stackrel{\rightarrow}{\mathcal{D}}_{\mu}\hspace*{-0.3em}(x)=\stackrel{\rightarrow}{\partial}_\mu-\ i\hat B_\mu(x), 
\nonumber\\
\xi_f=\frac{m_{f'}}{m_f+m_{f'}},\ \xi_{f'}=\frac{m_{f}}{m_f+m_{f'}}.
\nonumber
\end{gather}
Here $\mathcal{M}^a$ and $\Gamma^J$ are flavor and Dirac matrices corresponding to a given meson field, $\xi_f,\xi_{f'}$ provide that $x$ is the center of mass of a meson, $n,l$ are radial and orbital quantum numbers, respectively.
Radial part  $F_{nl}$ is defined by the propagator of the gluon fluctuations  charged with respect to the Abelian background, $T^{(l)}$ are irreducible tensors of four-dimensional rotation group.
Propagator of the quark with mass $m_f$ in the presence of the homogeneous Abelian (anti-)self-dual field  has the form 
\begin{align}
\label{quark_propagator}
S_f(x,y)&=\exp\left(-\frac{i}{2}\hat n x_\mu  B_{\mu\nu}y_\nu\right)H_f(x-y),
\\
\tilde H_f(p)&=\frac{1}{2\upsilon \Lambda^2} \int_0^1 ds e^{(-p^2/2\upsilon \Lambda^2)s}\left(\frac{1-s}{1+s}\right)^{m_f^2/4\upsilon \Lambda^2}
\nonumber\\
&\quad\times \left[\vphantom{\frac{s}{1-s^2}}p_\alpha\gamma_\alpha\pm is\gamma_5\gamma_\alpha f_{\alpha\beta} p_\beta
+m_f\left(P_\pm+P_\mp\frac{1+s^2}{1-s^2}-\frac{i}{2}\gamma_\alpha f_{\alpha\beta}\gamma_\beta\frac{s}{1-s^2}\right)\right],
\nonumber
\end{align}
where anti-Hermitian representation of Dirac matrices is used, and ``$\pm$'' signs are arranged in accordance with formula~\eqref{b-field}. The translation-invariant part $H_f$ of the propagator is an analytical function in the finite complex momentum plane and matches the behavior of free Dirac propagator at large Euclidean momentum. The analyticity of quark propagator is interpreted as confinement of dynamical quarks.

Overall, the mass spectrum of the  ground state and excited mesons composed of light and heavy quarks is described rather accurately, in complete agreement with
expectations based on confinement (Regge mass spectrum of radially and orbitally excited states) and chiral symmetry breaking (light pseudoscalar and heavy vector nonets, etc.) as well as asymptotic heavy-quark relations.  

A peculiar property of the quark propagator in the mean gluon field under consideration is that, unlike the case of the pseudoscalar and vector ground-state mesons, there are no real solutions to equation~\eqref{mass-eq} for their parity partners, the  ground-state scalar and axial mesons. Axial and scalar mesons appear with a mass above $1$ GeV in the hyperfine splitting of the orbital excitations. The inverse propagators of pseudoscalar and scalar meson fields are shown in Fig.~\ref{figure_pi_propagator} which manifestly illustrates the absence of a pole for the scalar field propagator.  This feature
 is particularly relevant to the present study. Contributions of intermediate scalar meson-like fields to various processes are available and can be computed  but  a controversial  issue of  existence  of the light scalar mesons does not occur.
\begin{figure}
\includegraphics[scale=1]{pi_sigma_propagator}
\caption{Left-hand side of Eq.~\ref{mass-eq} (which is proportional to inverse propagator) for pion and ground-state scalar quark-antiquark field with respect to Euclidean momenta $p^2$. The zero of inverse pion propagator (correspondingly, the pole of propagator $D_\pi(p^2)$) is located at $-M_\pi^2$, while inverse propagator $D_S^{-1}(p^2)$ has no zeroes at real momenta.\label{figure_pi_propagator}}
\end{figure}

Though the meson-quark coupling constant  $h_\mathcal{Q}$ is obviously undefined for meson-like composite fields if corresponding Eq.~\eqref{mass-eq} has no solutions, it is convenient to keep universal notation. Such fields can only be virtual, and $h_\mathcal{Q}$ cancels out in final expressions ($h^2_\mathcal{Q}$ for two vertices cancel $h_\mathcal{Q}^{-2}$ in propagator $D_\mathcal{Q}$).

Electromagnetic interactions are included in gauge-invariant way with the help of expansions (see~\cite{Nedelko:2016gdk,Nedelko:2021dsh})
\begin{align}
\nonumber
S_f(x,y|A)&=S_f(x,y)+\sum_{n=1}^\infty \left(Q_fe\right)^n \int dz_1\cdots \int dz_n S_f(x_1,z_1)\gamma_{\mu_1}A_{\mu_1}(z_1)\cdots S_f(z_{i-1},z_i)\gamma_{\mu_i}A_{\mu_i}(z_i)\cdots S_f(z_n,y),\\
\nonumber
V_\mathcal{Q}(x|A)&=V_\mathcal{Q}(x)+\sum_{n=1}^\infty e^n \int dz_1\cdots \int dz_n V_{\mathcal{Q};\mu_1\dots\mu_n}(x,z_1,\dots,z_n)A_{\mu_1}(z_1)\cdots A_{\mu_n}(z_n),
\end{align}
where $Q_f$ is electric  charge of a quark with flavor $f$ in units of electron charge $e$, and meson-photon vertices appear due to nonlocality of meson-quark interactions. 
One-photon and two-photon meson vertices are given by
\begin{gather}
\label{1-photon-vertex}
V_{\mathcal{Q};\mu}(x;q)= \int_0^1 d\tau \frac{1}{\tau} \frac{\partial}{\partial q_\mu} \left\{-Q_f V_{\mathcal{Q}}\left( \stackrel{\leftrightarrow}{\mathcal{D}}(x) -iq\tau\xi\right)+ Q_{f'} V_{\mathcal{Q}}\left(\stackrel{\leftrightarrow}{\mathcal{D}}(x) + iq\tau\xi'\right)\right\},\\
\label{2-photon-vertex}
\begin{split}
V_{\mathcal{Q};\mu\nu}(x;q_1,q_2)=\frac{1}{2}\int_0^1 d\tau_1 \int_0^1 d\tau_2 \frac{1}{\tau_1 \tau_2}\frac{\partial}{\partial q_{1\mu}} \frac{\partial}{\partial q_{2\nu}}  
& \left[Q_f Q_{f'}V_\mathcal{Q}\left( \stackrel{\leftrightarrow}{\mathcal{D}}(x) -iq_1\tau_1\xi -iq_2\tau_2\xi\right)\right.
\\
&-Q_f Q_{f'} V_\mathcal{Q}\left( \stackrel{\leftrightarrow}{\mathcal{D}}(x) -iq_1\tau_1\xi +iq_2\tau_2\xi\right)\\
&-Q_{f'} Q_f V_\mathcal{Q}\left( \stackrel{\leftrightarrow}{\mathcal{D}}(x) +iq_1\tau_1\xi -iq_2\tau_2\xi\right)\\
&\left.+Q_{f'} Q_{f'} V_\mathcal{Q}\left( \stackrel{\leftrightarrow}{\mathcal{D}}(x) +iq_1\tau_1\xi +iq_2\tau_2\xi\right)\right].
\end{split}
\end{gather}


Effective meson action takes the form
\begin{align}
\nonumber
Z&={\cal N}
\int D\phi_{\cal Q}\int DA_\mu
\exp\left\{-\frac{1}{4}\int d^4x\ F_{\mu\nu}F_{\mu\nu}-\frac{\Lambda^2}{2}\frac{h^2_{\cal Q}}{g^2 C^2_\mathcal{Q}}\int d^4x 
\phi^2_{\cal Q}(x)
-\Gamma(A)
-\sum\limits_{k=2}^\infty\frac{1}{k}W_k[\phi|A]\right\},
\\
\nonumber
W_k[\phi|A]&=
\sum\limits_{{\cal Q}_1\dots{\cal Q}_k}h_{{\cal Q}_1}\dots h_{{\cal Q}_k}
\int d^4x_1\dots\int d^4x_k
\Phi_{{\cal Q}_1}(x_1)\dots \Phi_{{\cal Q}_k}(x_k)
\Gamma^{(k)}_{{\cal Q}_1\dots{\cal Q}_k}(x_1,\dots,x_k|A),
\end{align}
where
\begin{equation*}
\Gamma(A)=\int d\sigma_B \mathrm{Tr}\log\left[1+ Qe\gamma_\mu A_\mu(x)S(x,y)\right],
\end{equation*}
and
$\Gamma^{(k)}_{{\cal Q}_1\dots{\cal Q}_k}(x_1,\dots,x_k|A)$ are obtained from $\Gamma^{(k)}_{{\cal Q}_1\dots{\cal Q}_k}(x_1,\dots,x_k)$ by substitutions
\begin{equation*}
S_f(x,y)\to S_f(x,y|A),\quad
V_\mathcal{Q}(x)\to V_\mathcal{Q}(x|A).
\end{equation*}

The free  parameters of the model are scale $\Lambda$ (related to the scalar gluon condensate $\langle g^2F^2\rangle$), infrared limits of dynamical quark masses and strong coupling $\alpha_s$ which has been determined by fitting to the masses of mesons $\pi,\rho,K,K^*,J/\psi,\Upsilon,\eta'$ (for details see~\cite{Nedelko:2016gdk}). Yet one more parameter is the mean domain size related to the topological susceptibility of QCD vacuum, but it is irrelevant to the present study, and we skip it.

As it has been mentioned,  diagonalization of the quadratic part of the effective action~\eqref{effective_meson_action} with respect to the radial quantum number is a part of the calculation procedure. In practice some finite number of excited states can be taken into account.  
As it has been analyzed in~\cite{Nedelko:2016gdk}, though typically about five lowest radial states have to be taken into account for robust stability of the computation, a consistent overall description of the mass spectrum of mesons  is achieved irrespective to a number of accounted  radial excitations. Just the values of the free parameters  have to be adjusted when the number of accounted radial excitation changes. 
 
\section{Evaluation of polarizabilities\label{section_polarizabilities}}

The polarizabilities are defined by the Compton scattering amplitude
\begin{equation*}
P(p)+\gamma(q,\varepsilon)\to P(p')+\gamma(q',\varepsilon')
\end{equation*}
of a pseudoscalar meson $P$
\begin{equation*}
{}_\text{out}\langle P(p')\gamma(q',\varepsilon')|P(p)\gamma(q,\varepsilon) \rangle_\text{in}=i(2\pi)^4\delta^{(4)}(p'+q'-p-q)\varepsilon^\mu(q)\varepsilon^{*\nu}(q')M_{\mu\nu}.
\end{equation*}
We concentrate on the electric $\alpha_\text{E}$ and magnetic $\beta_\text{M}$ dipole polarizabilities which appear in expansion of the amplitude in small photon momenta as
\begin{equation*}
\varepsilon^\mu(q)\varepsilon^{*\nu}(q')M_{\mu\nu}=-2e^2\vec{\varepsilon}\cdot \vec{\varepsilon} {\,}'^{*}+8\pi M(\alpha_\text{E}\omega\omega'\,\vec{\varepsilon}\cdot \vec{\varepsilon} {\,}'^{*}+\beta_\text{M}\left(\vec{\varepsilon}\times \vec{q}\,\right)\cdot\left(\vec{\varepsilon}{\,}'^{*}\times \vec{q}{\,}'\right))+\dots,
\end{equation*}
where $M$ is the mass of a pseudoscalar meson.
The tensor $M_{\mu\nu}$ can be separated into two parts
\begin{equation}
\label{born_and_nonborn}
M_{\mu\nu}=M_{\mu\nu}^\text{Born}+M_{\mu\nu}^\text{NB},
\end{equation}
where the part $M_{\mu\nu}^\text{NB}$ describes the response of a meson as a composite system to the applied electromagnetic field. The term $M_{\mu\nu}^\text{Born}$ given by
\begin{equation}
\label{born_scattering}
M_{\mu\nu}^\text{Born}=e^2\left[2g_{\mu\nu}-\frac{(2p_\mu+q_\mu)(2p'_\nu+q'_\nu)}{(p+q)^2-M^2}-\frac{(2p_\nu-q'_\nu)(2p'_\mu-q_\mu)}{(p-q')^2-M^2}\right]
\end{equation}
describes real Compton scattering of a structureless pseudoscalar particle.

In the case of real Compton scattering ($q^2={q'}^2=0$, $\varepsilon^\mu(q) q_\mu=\varepsilon^{*\mu}(q')q'_\mu=0$), the tensor $\mathcal{M}_{\mu\nu}$ contains only two independent tensor structures~\cite{Tarrach:1975tu,Bardeen:1968ebo,Lvov:2001zdg}
\begin{align}
\label{compton_tensor_parametrization}
\varepsilon^\alpha(q)\varepsilon^{*\beta}(q')M_{\alpha\beta}&=\varepsilon^\alpha(q)\varepsilon^{*\beta}(q')\left(A T_{1\alpha\beta}+B T_{2\alpha\beta}\right),\\
\nonumber
T_{1\alpha\beta}&=-\frac{t}{2}g_{\alpha\beta}-q_{\beta}q'_{\alpha},\\
\nonumber
T_{2\alpha\beta}&=-\frac{t}{2}P_\alpha P_\beta+\nu^2 g_{\alpha\beta}-\nu(P_\alpha q_\beta+P_\beta q'_\alpha),  
\end{align}
where
\begin{equation*}
P=\frac{p+p'}{2},\quad s=(p+q)^2,\quad t=(q-q')^2,\quad u=(p-q')^2,\quad \nu=\frac{1}{4}(s-u).
\end{equation*}
In accordance with Eqs.~\eqref{born_and_nonborn} and~\eqref{born_scattering}, amplitudes $A$ and $B$ can be split in two parts
\begin{align*}
A(\nu^2,t)&=A^\text{Born}(\nu^2,t)+A^\text{NB}(\nu^2,t),\\
B(\nu^2,t)&=B^\text{Born}(\nu^2,t)+B^\text{NB}(\nu^2,t),
\end{align*}
where $A^\text{Born},B^\text{Born}$ are given by
\begin{equation*}
A^\text{Born}(\nu^2,t)=-\frac{e^2t}{(s-M_\pi^2)(u-M_\pi^2)},\quad B^\text{Born}(\nu^2,t)=-\frac{8e^2}{(s-M_\pi^2)(u-M_\pi^2)}.
\end{equation*}
The electric $\alpha_\text{E}$ and magnetic $\beta_\text{M}$ polarizabilities are related to $A^\text{NB},B^\text{NB}$ by means of equations (see \cite{Guiasu:1978dz,Guiasu:1979sz})
\begin{align}
\label{alpha_plus_beta_definition}
\alpha_\text{E}+\beta_\text{M}&=-\frac{M}{8\pi}B^\text{NB}(0,0),\\
\label{alpha_minus_beta_definition}
\alpha_\text{E}-\beta_\text{M}&=-\frac{1}{8\pi M}\left(2A^\text{NB}(0,0)+M^2B^\text{NB}(0,0)\right).
\end{align}

The diagrams that contribute to Compton tensor are shown in Fig.~\ref{compton_scattering_diagrams}.
\begin{figure}
\begin{tabular}{c@{\hspace*{1em}}c@{\hspace*{1em}}c@{\hspace*{1em}}c}
\includegraphics[scale=1]{compton_scattering1}&
\includegraphics[scale=1]{compton_scattering2}&
\includegraphics[scale=1]{compton_scattering3}&
\includegraphics[scale=1]{compton_scattering4}\\
(a)&(b)&(c)&(d)
\end{tabular}
\caption{
Diagrams contributing to Compton scattering tensor $M_{\mu\nu}$. The grey circles denote all possible one-particle irreducible contributions.
\label{compton_scattering_diagrams}}
\end{figure}
Consider contribution of diagrams (a),(b),(c) in Fig.~\ref{compton_scattering_diagrams} which form gauge-invariant combination.
\begin{figure}
\begin{tabular}{c@{\hspace*{1em}}c@{\hspace*{1em}}c@{\hspace*{1em}}c@{\hspace*{1em}}c}
\includegraphics[scale=.75]{2p2m_1}&
\includegraphics[scale=.75]{2p2m_2}&
\includegraphics[scale=.75]{2p2m_3}&
\includegraphics[scale=.75]{2p2m_4}&
\includegraphics[scale=.75]{2p2m_5}\\
(a)&(b)&(c)&(d)&(e)
\\[1ex]
\includegraphics[scale=.75]{2p2m_6}&
\includegraphics[scale=.75]{2p2m_7}&
\includegraphics[scale=.75]{2p2m_8}&
\includegraphics[scale=.75]{2p2m_9}&
\includegraphics[scale=.75]{2p2m_10}\\
(f)&(g)&(h)&(i)&(j)
\end{tabular}
\caption{The diagrams contributing to $\Gamma_{\mu\nu}$. The diagrams related by crossings are not shown.\label{figure_2p2m}}
\end{figure}
The diagrams (b) and (c) contain kinematic singularities that are canceled by corresponding Born terms in Eq.~\eqref{born_scattering}.
One can notice that among these diagrams, only diagram (a) contains tensors proportional to $g_{\alpha\beta}$ and can be parametrized as
\begin{align*}
M^{(\text{a})}_{\alpha\beta}&=\left(-\frac{t}{2}A(\nu^2,t)+\nu^2B(\nu^2,t)\right) g_{\alpha\beta}+ \text{other tensor structures},\\
M^{(\text{a})\text{NB}}_{\alpha\beta}&=\left(-\frac{t}{2}A(\nu^2,t)+\nu^2B(\nu^2,t)-2\right) g_{\alpha\beta}+ \text{other tensor structures},
\end{align*}
according to formulas~\eqref{born_scattering} and~\eqref{compton_tensor_parametrization}. The amplitudes $A,B$ that appear in definition of polarizabilities~\eqref{alpha_plus_beta_definition},\eqref{alpha_minus_beta_definition} can be extracted from coefficient of $g_{\alpha\beta}$ in $M_{\alpha\beta}^\text{NB}$ with the help of formulas
\begin{align*}
A^\text{NB}(0,0)&=-2\frac{\partial}{\partial t}\left.\left(-\frac{t}{2}A(\nu^2,t)+\nu^2B(\nu^2,t)-2\right)\right|_{t=0,\nu=0},\\
B^\text{NB}(0,0)&=\frac{1}{2}\frac{\partial^2}{\partial^2 \nu}\left.\left(-\frac{t}{2}A(\nu^2,t)+\nu^2B(\nu^2,t)-2\right)\right|_{t=0,\nu=0}.
\end{align*}
It is therefore sufficient to calculate only diagrams (a) of gauge-invariant combination of (a),(b) and (c) in order to extract electric and magnetic dipole polarizabilities. This is more straightforward because diagram (a) does not contain kinematic poles.
One-loop contributions of this type are shown in Fig.~\ref{figure_2p2m} (the Feynman rules in Euclidean space are given by formulas~\eqref{effective_meson_action},\eqref{Gammak} and \eqref{barG} for loops, Eqs.~\eqref{qmvert},\eqref{1-photon-vertex},\eqref{2-photon-vertex} describe non-local vertices which should be multiplied by corresponding $h_\mathcal{Q}$, the local vertices are the same as in QED, the quark propagators are defined by Eq.\eqref{quark_propagator}). For example, diagram (e) in Fig.\ref{figure_2p2m} corresponds to
\begin{equation*}
\begin{split}
e^2h_P^2 \int\! d\sigma_B\int\! d^4x \int\! d^4y \int\! d^4z &\exp(ipx+iqx-ip'y-iq'z)\\
&\times(-1) \text{Tr}V_{;\mu}^{aP00}(x;-q)S(x,y)V^{bP00}(y)S(y,z)Q\gamma_\nu S(z,y)\\
&+\text{crossed term},
\end{split}
\end{equation*}
for two external pseudoscalar ground-state mesons ($\delta$-function resulting from conservation of total momentum is included). The loop integrals are finite due to non-local meson vertices, so no regularization is needed. With meson vertices and quark propagators given by formulas~\eqref{qmvert} and~\eqref{quark_propagator}, the space and momentum integrals are Gaussian and can be computed analytically. The averaging over background field is performed with the help of formula~\eqref{averaging_over_vacuum_field}. The remaining integrals over proper times that originate from vertices and propagators are computed numerically. 

The contribution to the amplitude corresponding to diagram (d) in Fig.~\ref{compton_scattering_diagrams} with intermediate scalar fields is separately gauge-invariant and can be represented as
\begin{equation}
\label{intermediate_scalar_meson}
\varepsilon^\alpha(q)\varepsilon^{*\beta}(q')M^{(d)}_{\alpha\beta}=\varepsilon^\alpha(q)\varepsilon^{*\beta}(q')\sum_S \Gamma_{\alpha\beta}^{S\gamma\gamma} D_S h_P^2\Gamma_{SPP},
\end{equation}
where $\Gamma_{S\gamma\gamma}$ parametrizes $S\to \gamma\gamma$ subprocess, $D_S$ is a scalar meson field propagator,
$h_P^2\Gamma_{SPP}(t)$ describes transition of scalar field to a couple of pseudoscalar mesons (corresponding one-loop diagrams are shown in Fig.~\ref{figure_intermediate_scalar}). Diagram $M^{(d)}$ does not contain tensor $T_{2\alpha\beta}$ (see formula~\eqref{compton_tensor_parametrization}) and hence contributes only to amplitude $A$.
\begin{figure}
\begin{tabular}{c@{\hspace*{1em}}c@{\hspace*{1em}}c}
\includegraphics[scale=.75]{meson_v1p1_2photons}&
\includegraphics[scale=.75]{meson_v0p2_2photons}&
\includegraphics[scale=.75]{meson_2mesons}\\[1ex]
(a)&(b)&(c)
\end{tabular}
\caption{
One-loop diagrams contributing to $M_{\alpha\beta}^S$ in formula~\eqref{intermediate_scalar_meson} (crossed diagrams are not shown). Diagrams (a),(b) are related to $\Gamma_{S\gamma\gamma}$, diagram (c) results in $h_P^2\Gamma_{SPP}$.
\label{figure_intermediate_scalar}}
\end{figure}

The computational complexity of the amplitudes which we need to evaluate in order to extract polarizabilities quickly grows with increasing the radial number $n$
taken into account for diagonalization of the quadratic part of the action. In the present paper, for numerical calculation   only $n=0$ states are taken into account, such that the  matrix $\mathcal{O}$ in Eq.~\eqref{physical_meson_fields} is reduced to the  unity matrix, and the higher radial states are neglected.
The values of parameters in this lowest approximation with respect to the ``radial excitation mixing'' are  given in Table~\ref{values_of_parameters} (for detailed discussion  see~\cite{Nedelko:2016gdk}). The values of polarizabilities found in the present work are presented in Table~\ref{polarizabilities_table}.

\begin{table}
\begin{tabular}{|*{4}{@{\hspace*{0.5em}}c@{\hspace*{0.5em}}|}}
\hline
$m_{u/d}$(MeV)&$m_s$(MeV)&$\Lambda$(MeV)&$\alpha_s$\\
\hline
$174.8$&$393.3$&$439.7$&$6.23$\\
\hline
\end{tabular}
\caption{Values of parameters used for calculations of polarizabilities.\label{values_of_parameters}}
\end{table}


\begin{table}
\begin{tabular}{|c@{\hspace*{1em}}c|c|c|c||c|c|}
\hline
&&Compton scattering&intermediate scalar fields&total&experiment&ChPT\\
\hline
\multirow{2}{*}{$\pi^\pm$}&$\alpha_\text{E}+\beta_\text{M}$&$0.13$&0&$0.13$&$0.5\pm 0.5_\text{stat}$~\cite{COMPASS:2014eqi}&$0.16$~\cite{Gasser:2006qa}\\
&$\alpha_\text{E}-\beta_\text{M}$&$4.82$&$1.48$&$6.3$&$4.0\pm 1.2_\text{stat}\pm 1.4_\text{syst}$~\cite{COMPASS:2014eqi}&$5.7$~\cite{Gasser:2006qa}\\
\hline
\multirow{2}{*}{$\pi^0$}&$\alpha_\text{E}+\beta_\text{M}$&$0.71$&0&$0.71$&$0.98\pm 0.03$~\cite{Filkov:1998rwz}&$1.15$~\cite{Gasser:2005ud}\\
&$\alpha_\text{E}-\beta_\text{M}$&$-0.26$&$1.48$&$1.22$&$-1.6\pm 2.2$~\cite{Filkov:1998rwz}&$-1.9$~\cite{Gasser:2005ud}\\
\hline
\multirow{2}{*}{$K^\pm$}&$\alpha_\text{E}+\beta_\text{M}$&$0.41$&0&$0.41$&&\\
&$\alpha_\text{E}-\beta_\text{M}$&$1.22$&$0.47$&$1.69$&&\\
\hline
\multirow{2}{*}{$K^0,\bar{K}^0$}&$\alpha_\text{E}+\beta_\text{M}$&$0.62$&0&$0.62$&&\\
&$\alpha_\text{E}-\beta_\text{M}$&$-0.29$&$0.17$&$-0.12$&&\\
\hline
\end{tabular}
\caption{Numerical results for polarizabilities of pseudoscalar mesons in Gaussian units of $10^{-4}\ \text{fm}^3$. ``Compton scattering'' corresponds to contribution of diagrams (a),(b),(c) in Fig.~\ref{compton_scattering_diagrams} with corresponding intermediate pseudoscalar mesons. ``Intermediate scalar fields'' labels contribution diagram (d) in Fig.~\ref{compton_scattering_diagrams} with intermediate scalar quark-antiquark fields.\label{polarizabilities_table}}
\end{table}

\section{Summary}
We investigated dipole polarizabilities of the light pseudoscalar mesons in the framework of the nonlocal effective meson action obtained within the mean field approach to QCD vacuum. The specific feature of the approach is the possibility to evaluate the contribution of intermediate scalar meson fields composed of light quarks to pion and kaon polarizabilities avoiding at the same time the issue of the existence of light scalar mesons as physical meson states. 

The values of charged pion polarizabilities calculated in the present study are in agreement with COMPASS data and most recent two-loop ChPT calculation~\cite{Gasser:2006qa}. In contrast with Refs.~\cite{Ivanov:1991kw,Filkov:1998rwz,Dorokhov:1997rv}, the main contribution comes from one-loop diagrams, while contribution of intermediate scalar field is less important. Prediction of Lattice QCD for polarizabilities depends on parameters such as the lattice volume, lattice spacings, quark masses. The value of dipole magnetic polarizability of charged pion found in paper~\cite{Luschevskaya:2015cko} with the finest lattice is $\beta_{\pi^\pm}=-2.06\pm 0.76\times 10^{-4}\ \text{fm}^3$ supports data of COMPASS collaboration~\cite{COMPASS:2014eqi}, ChPT~\cite{Gasser:2006qa} and findings of this paper. 

The computation has been performed in the lowest approximation with respect to the mixing of radially excited states in the effective meson action~\eqref{effective_meson_action}, and it would be interesting and important to check the stability of obtained results in this respect by accounting the higher radial excitations,
which will also allow one to estimate the polarizabilities of  radially   excited pion and kaon states.  However, the latter has mostly  purely theoretical importance as experimental measurement seems to be hardly achievable.


\section{Acknowledgments}
We are grateful for resources provided by ``Govorun'' supercomputer at the Joint Institute for Nuclear Research.

