††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Real-complex transition driven by quasiperiodicity: a new universality class
beyond $\mathcal{PT}$ symmetric one

Tong Liu  Department of Applied Physics, School of Science, Nanjing University
of Posts and Telecommunications, Nanjing 210003, China  Xu Xia  Chern
Institute of Mathematics and LPMC, Nankai University, Tianjin 300071, China

(September 3, 2024)

###### Abstract

We study a one-dimensional lattice model subject to non-Hermitian
quasiperiodic potentials. Firstly, we strictly demonstrate that there exists
an interesting dual mapping relation between $|a|<1$ and $|a|>1$ with regard
to the potential tuning parameter $a$. The localization property of $|a|<1$
can be directly mapping to that of $|a|>1$, the analytical expression of the
mobility edge of $|a|>1$ is therefore obtained through spectral properties of
$|a|<1$. More impressive, we prove rigorously that even if the phase
$\theta\neq 0$ in quasiperiodic potentials, the model becomes
non-$\mathcal{PT}$ symmetric, however, there still exists a new type of real-
complex transition driven by non-Hermitian disorder, which is a new
universality class beyond $\mathcal{PT}$ symmetric class.

###### pacs:

71.23.An, 71.23.Ft, 05.70.Jk

##  I Introduction

In quantum physics, the conservation of energy and probability demands that a
closed system exhibits real energies, which leads to all physical observables
must be represented by Hermitian operators Shankar . However, when an open
quantum system couples to a surrounding environment, the Hamiltonian of the
system becomes non-Hermitian and the physical process can be effectively
described through a complex eigenenergy Lindblad . Recently, non-Hermitian
systems have been attracted growing interest motivated by experimental
realization of photonic lattices Photonic1 ; Photonic2 ; Photonic3 ; Photonic4
and ultracold atomic gases ultracold1 ; ultracold2 ; ultracold3 ; ultracold4 .
A celebrated paradigm in the non-Hermitian physics was parity-time
($\mathcal{PT}$) symmetric class, namely, a large class of non-Hermitian
Hamiltonians, discovered by Bender and Boettcher Bender1 ; Bender2 , can
exhibit entirely real spectra as long as they commute with the $\mathcal{PT}$
operator. From an intuitive point of view, a necessary (but not sufficient)
condition for $\mathcal{PT}$ symmetry to preserve is that the involved complex
potential should satisfy $V(n)=V^{*}(-n)$ in the discrete lattice, where $n$
is the lattice site number PT_1 ; PT_2 ; PT_3 ; PT_4 ; PT_5 ; PT_6 . Although
the impact of $\mathcal{PT}$ symmetry in real quantum systems is still
debated, broad research interests have been put into the study of
$\mathcal{PT}$ symmetry-breaking transitions between real and complex
eigenenergies in various non-Hermitian systems Skin1 ; Skin2 ; Yao1 ; Yao2 .

Here a natural question arises: can other class of non-Hermitian systems host
entirely real spectra beyond $\mathcal{PT}$ symmetric class? Disordered
systems provide a potential candidate random ; quasiperiodic . A paradigm to
understand the Anderson localization in one dimension (1D) is the famous
Aubry-André (AA) model AA ; Harper , one of the most studied example
displaying a localization-delocalization transition Laurent1 ; Laurent2 ;
Laurent3 ; Laurent4 . Another interesting aspect of AA-like models is the
presence of mobility edges separating extended from localized states
mobility_1 ; mobility_2 ; mobility_3 ; Liu1 . The combination of non-Hermitian
properties and disordered systems is an interesting and active research topic
Longhi1 ; Longhi2 ; Gong ; Kawabata . A seminal work dealing with disorder in
non-Hermitian systems is the Hatano-Nelson model Hatano1 ; Hatano2 , in which
an asymmetric hopping caused by an imaginary gauge field results in a
localization-delocalization transition Hatano3 . Other non-Hermitian models
with either random or quasiperiodic disorder have been investigated, in which
the localization properties of systems are the main concerns and contents
Jiang ; LiuY ; Zeng .

With regard to the real-complex transition of spectra, in a pioneering study,
Hamazaki et al. numerically found that a real-complex transition occurred in
non-Hermitian interacting systems without $\mathcal{PT}$ symmetry Hamazaki ,
the real energy region was located in the localized phase of the model with
time-reversal symmetry. In the further research, Liu et al. numerically
demonstrated that the real energy region was located in the extended phase of
a quasiperiodic $p$-wave superconductor chain without $\mathcal{PT}$ symmetry
Liu2 . However, these available studies are limited to numerical simulations,
there are very few rigorous results to clarify the underlying physics of the
real energy. In this work, we are devoted to introduce an exactly solvable
model and analytically demonstrate that there exists a new type of real-
complex transitions driven by non-Hermitian disorder, which is different from
$\mathcal{PT}$ symmetric class.

##  II Model and The dual mapping

We consider one-dimensional non-Hermitian quasiperiodic model, described by
the following eigenvalue equation,

$E\psi_{n}=\psi_{n+1}+\psi_{n-1}+\frac{V}{1-ae^{i(2\pi\alpha
n+\theta)}}\psi_{n},$ (1)

where $V$ is the complex potential strength, $\alpha$ is irrational, and
$\psi_{n}$ is the amplitude of wave function at the $n$th lattice. We choose
to unitize the nearest-neighbor hopping amplitude and a typical choice for
parameter $\alpha$ is $\alpha=(\sqrt{5}-1)/2$. The phase $\theta$ locates in
the center of this study, which is not simply taken as zero. $|a|\neq 1$ is
the potential tuning parameter. When $|a|<1$ and $\theta=0$, Eq. (1) is
$\mathcal{PT}$ symmetric, Liu et al. Liu1  derived the analytical expression
of the energy-dependent mobility edge by means of a generalized self-duality
relation. In this work we demonstrate that the $|a|>1$ case is the dual
counterpart of the $|a|<1$ case, and the localization property of $|a|<1$ can
be directly mapping to that of $|a|>1$. Moreover, we prove rigorously that
when $\theta\neq 0$, Eq. (1) becomes non-$\mathcal{PT}$ symmetric, however,
there still exists a real-complex transition in the energy spectrum, which is
unconventional and different from $\mathcal{PT}$ symmetric class.

Now we unveil a dual mapping relation between $|a|<1$ and $|a|>1$ of Eq. (1).
Note the non-Hermitian quasiperiodic potential $\frac{V}{1-ae^{ix}}$ can be
written as Taylor’s expansion,

$\frac{V}{1-ae^{ix}}=V\sum_{m=0}^{\infty}a^{m}e^{imx},$ (2)

then,

$\frac{V}{1-ae^{ix}}=\frac{-V}{1-1/ae^{-ix}}+V,$ (3)

thus, spectral properties of $\frac{V}{1-ae^{ix}}$ ($|a|<1$) can be dual
mapping to these of the potential $\frac{-V}{1-ae^{-ix}}+V$ ($|a|>1$).
Consequently, utilizing the mobility edge energy $E_{c}=a+1/a$ ($|a|<1$), with
regard to $|a|>1$ we can directly obtain the exact explicit form of the
mobility edge:

$E_{m}=V-a-1/a.$ (4)

To support the analytical result given above, we now present detailed
numerical analysis of Eq. (1). In the disordered system, the localization
property of wave functions can be measured by the inverse participation ratio
(IPR) IPR . For any given normalized wave function, the corresponding IPR is
defined as $\text{IPR}=\sum_{n=1}^{L}\left|\psi_{n}\right|^{4},$ which
measures the inverse of the number of sites being occupied by particles. It is
well known that the IPR of an extended state scales like $L^{-1}$ which
approaches zero in the thermodynamic limit. However, for a localized state,
since only finite number of sites are occupied, the IPR is finite even in the
thermodynamic limit. In Fig. 1 we show the numerical IPR diagram in the $[{\rm
Re}(E),V]$ plane, different colours of the eigenvalue curves indicate
different magnitudes of the IPR of the corresponding wave functions. The black
eigenvalue curves denote the extended states, and the bright yellow eigenvalue
curves denote the localized states. It is clearly demonstrating a mobility
edge separating localized from extended states along the blue line defined by
Eq. (4).

Figure 1: (Color online) The real part of eigenvalues of Eq. (1) and IPR as a
function of $V$ with the parameter $a=2$. The total number of sites is set to
be $L=500$. The blue solid lines represent the boundary between spatially
localized and delocalized states, i.e., the mobility edge $E_{m}=V-a-1/a$. In
the inset, we plot the typical probability density $|\psi|$ with eigenvalues
above and below $E_{m}$ ($V=2$), as expected, they are the extended state
(black) and the localized state (red), respectively.

##  III Real-complex transition

Besides the localization, the most novel discovery of this work is that there
exists a new type of real-complex transitions in Eq. (1) driven by non-
Hermitian disorder, rather than $\mathcal{PT}$ symmetry.

For the sake of simplicity, here we focus on the mathematical proof of the
$|a|<1$ case. First let us multiply both sides of Eq. (1) by $e^{ink}$ and sum
over $n$. After setting $f(k)=\sum_{n}e^{ink}\psi_{n}$, one obtains

$a[E-2\cos(k+2\pi\alpha)]f(k+2\pi\alpha)=[E-V-2\cos(k)]f(k),$ (5)

According to Sarnak’s method Sarnak , the spectrum of Eq. (5) is govern by a
charateristic function defined as

$\displaystyle G(E)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\log\left|E-V-2\cos(k)\right|dk$
(6)
$\displaystyle-\frac{1}{2\pi}\int_{0}^{2\pi}\log\left|E-2\cos(k+2\pi\alpha)\right|dk,$

Here, we don’t show more complicated mathematical proofs, but directly quote
Sarnak’s Lemma Sarnak :

(i) There is no energy spectrum when $G(E)<\log|a|$.

(ii) When $G(E)=\log|a|$, Eq. (5) has dense spectrum, the corresponding wave
functions are localized.

(iii) When $G(E)>\log|a|$, wave functions of Eq. (5) are extended. However,
the spectrum of the system must be within the energy interval
$U_{E}=[V-2,V+2]$, which is derived by the equation $E-V-2\cos(k)=0$.

The conclusion (iii) is the central result, it demonstrates that when the
system is in the extended phase, the spectrum must be restricted to a real
interval, not a complex number, it should be noted that the above proof does
not involve $\mathcal{PT}$ symmetry.

When the spectrum $E$ is restricted in $U_{E}$, Eq. (6) can be reduced to

$G(E)=\log|\frac{E-\sqrt{E^{2}-4}}{2}|$ (7)

From the conclusion (ii), we can obtain the metal-insulator transition
boundary, i.e., the mobility edge $E_{c}=a+1/a$ derived by the equation
$G(E_{c})=\ln(a)$. Remarkably, this expression is exactly alike that derived
by the generalized self-duality relation. And the total real spectrum region
is enclosed by three lines $E_{min}=V-2$, $E_{max}=V+2$ and $E_{c}=a+1/a$.

Figure 2: (Color online) The eigenvalues of Eq. (1) with the parameter $a=0.5$
under open boundry conditions. The total number of sites is set to be $L=500$.
(a) and (c) use the same potential strength $V=2$ but different phase $\theta$
for comparison. Similarly, (b) and (d) use the same potential strength $V=3$
but different phase $\theta$ for comparison. The blue solid lines represent
the boundary between the real and complex energy spectrum, it is clearly shown
that the real spectrum region of four cases is restricted in $E_{min}=V-2$ and
$E_{c}=a+1/a$.

The above derivation ignores the phase $\theta$. When $\theta=0$, Eq. (1)
still preserve the $\mathcal{PT}$ symmetry, this may be the reason why most
works LiuY ; Zeng ; PTPT1 ; PTPT2  classified this type of real-complex
transition as the $\mathcal{PT}$ symmetry. To clarify the difference, we
further focus on the $\theta\neq 0$ case, now Eq. (1) obviously does not
satisfy the $\mathcal{PT}$ symmetry. Concerning the phase $\theta$ in
quasiperiodic potentials, Simon uses the mathematical ergodic theory to prove
that the eigenvalues are almost independent of the phase Simon . At the same
time, he also pointed out that the unique ergodic is completely independent of
the phase, but he did not give the details of the proof. We add here the proof
that eigenvalues of our model is completely phase-independent. We mark the
tridiagonal matrix of the Hamiltonian as $H_{\theta_{1}}$ and $H_{\theta_{2}}$
when two phases are $\theta_{1}$ and $\theta_{2}$, respectively. All
eigenvalues of $H_{\theta_{1}}$ and $H_{\theta_{2}}$ are denoted as
$\Sigma_{1}$ and $\Sigma_{2}$, consequently.

Then for any integer $n$, there exists the following estimate

$\displaystyle\|H_{\theta_{1}}-H_{\theta_{2}}\|$
$\displaystyle\leqslant\frac{2V}{(1-a)^{2}}|\sin(\theta_{1}-\theta_{2}-nw){\pi}|$
(8) $\displaystyle\leqslant\frac{2\pi
V}{(1-a)^{2}}{\left|\theta_{1}-\theta_{2}-n\omega\right|_{\mathbb{R}/\mathbb{Z}}},$

thus, we obtain

$h(\Sigma_{1},\Sigma_{2})\leq\frac{2\pi
V}{(1-a)^{2}}{\left|\theta_{1}-\theta_{2}-n\omega\right|_{\mathbb{R}/\mathbb{Z}}}.$
(9)

where $h$ represents the minimum distance of two sets $\Sigma_{1}$ and
$\Sigma_{2}$. Because the frequency $\alpha$ is an irrational number, we can
take a list of $n_{k}$ to make
$\left|\theta_{1}-\theta_{2}-n_{k}\omega\right|_{\mathbb{R}/\mathbb{Z}}$
approach zero, which means $h=0$. Consequently, we strictly prove
$\Sigma_{1}=\Sigma_{2}$, irrelevant to $\theta_{1}$ and $\theta_{2}$. Thus
even if $\theta\neq 0$, i.e., it doesn’t exist the $\mathcal{PT}$ symmetry,
the real spectrum of the system still stays the same, these results strongly
demonstrate that the real-complex transition driven by $\exp(ix)$-like non-
Hermitian disorder belongs to a unique class, rather than the $\mathcal{PT}$
symmetric class.

Figure 3: (Color online) The eigenvalues of Eq. (1) with the parameter $a=0.5$
under open boundry conditions, $\alpha$ is set to a rational number $1/4$. The
total number of sites is set to be $L=500$. In this situation, the $\theta\neq
0$ case has different spectrum with the $\theta=0$ case, and the real-complex
transition is within the logical framework of the $\mathcal{PT}$ symmetry
broken.

In order to verify the theoretical analysis, we implement numerical
calculations of eigenvalues. As shown in Fig. 2, even if the phase $\theta$ is
different, the eigenvalues of the system with the same $V$ are same, and the
mobility edge $E_{c}=a+1/a$ exactly separates the real from complex spectrum,
thus numerical results are in excellent agreement with the analytical results.
We have also checked other combinations of parameters and get the same results
as expected. In addition, there exist several deviated points in Fig. 2, this
is because, not like finite periodic systems, disordered systems are
statistical and based on large samples. The proportion of deviated points will
decrease sharply and be almost negligible with the increase of sample size.

A possible reason of being classified as the $\mathcal{PT}$ symmetric class in
previous works is to ignore the difference of the rational and irrational of
the parameter $\alpha$. For $\alpha$ being rational, the system is perfectly
periodic, in this situation the real-complex transition is indeed driven by
the $\mathcal{PT}$ symmetry, see Fig. 3. When the phase $\theta=0$, the system
preserves the $\mathcal{PT}$ symmetry, the eigenvalues are real, while the
phase $\theta\neq 0$, the system doesn’t preserve the $\mathcal{PT}$ symmetry,
the eigenvalues are complex. For $\alpha$ being irrational, as we proved
above, this real-complex transition violates the logical framework of the
$\mathcal{PT}$ symmetry broken, hence can’t be classified as the
$\mathcal{PT}$ symmetric class.

##  IV Summary

In summary, we have shown an interesting dual mapping between $|a|<1$ and
$|a|>1$ in a non-Hermitian quasiperiodic model, the spectral properties of the
two can be mapping to each other. More impressively, we mathematically prove
that the eigenenergies in the extended phase are fully real, and an
unconventional real-complex transition occurs accompanied by the metal-
insulator transition. The meaning of unconventionality is that even if our
model doesn’t preserve the conventional $\mathcal{PT}$ symmetry, the real-
complex transition in the spectrum also appears exactly as our theoretical
predictions. Our findings act in cooperation with the previous numerical
results of non-Hermitian disorder, and establish that the real-complex
transition driven by non-Hermitian disorder is a new universality class, which
is different from $\mathcal{PT}$ symmetric class.

###### Acknowledgements.

T. L. acknowledges the Natural Science Foundation of Jiangsu Province (Grant
No. BK20200737) and NUPTSF (Grant No. NY220090 and No. NY220208). X. X. is
supported by Nankai Zhide Foundation.

## References

  * (1) R. Shankar, Principles of Quantum Mechanics (Plenum Press, 1994).
  * (2) G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * (3) L. Feng, R. El-Ganainy, and L. Ge, Non-Hermitian photonics based on parity time symmetry, Nat. Photonics 11, 752 (2017).
  * (4) A. Cerjan, S. Huang, M. Wang, K. P. Chen, Y. Chong, and M. C. Rechtsman, Experimental realization of a Weyl exceptional ring, Nat. Photonics 13, 623 (2019).
  * (5) L. Xiao, X. Zhan, Z. H. Bian, K. K. Wang, X. Zhang, X. P.Wang, J. Li, K. Mochizuki, D. Kim, N. Kawakami, W. Yi, H. Obuse, B. C. Sanders, and P. Xue, Observation of topological edge states in parity time symmetric quantum walks, Nat. Phys. 13, 1117 (2017).
  * (6) J. M. Zeuner, M. C. Rechtsman, Y. Plotnik, Y. Lumer, S. Nolte, M. S. Rudner, M. Segev, and A. Szameit, Observation of a Topological Transition in the Bulk of a Non-Hermitian System, Phys. Rev. Lett. 115, 040402 (2015).
  * (7) I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (8) J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Observation of parity-time symmetry breaking transitions in a dissipative Floquet system of ultracold atoms, Nat. Commun. 10, 855 (2019).
  * (9) W. Gou, T. Chen, D. Xie, T. Xiao, T.-S. Deng, B. Gadway, W. Yi, and B. Yan, Tunable Nonreciprocal Quantum Transport through a Dissipative Aharonov-Bohm Ring in Ultracold Atoms, Phys. Rev. Lett. 124, 070402 (2020).
  * (10) F. Alex An, E. J. Meier, and B. Gadway, Engineering a Flux-Dependent Mobility Edge in Disordered Zigzag Chains, Phys. Rev. X 8, 031045 (2018).
  * (11) C. M. Bender and S. Boettcher, Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * (12) C. M. Bender, Making sense of non-Hermitian Hamiltonians, Rep. Prog. Phys. 70, 947 (2007).
  * (13) C. Yuce, Topological phase in a non-Hermitian PT symmetric system, Phys. Lett. A 379, 1213 (2015).
  * (14) D.-W. Zhang, L.-Z. Tang, L.-J. Lang, H. Yan, and S.-L. Zhu, Non-Hermitian topological Anderson insulators, Sci. China-Phys. Mech. Astron. 63, 267062 (2020).
  * (15) D. S. Borgnia, A. J. Kruchkov, and R.-J. Slager, Non-Hermitian Boundary Modes and Topology, Phys. Rev. Lett. 124, 056802 (2020).
  * (16) K. Yamamoto, M. Nakagawa, K. Adachi, K. Takasan, M. Ueda, and N. Kawakami, Theory of Non-Hermitian Fermionic Superfluidity with a Complex-Valued Interaction, Phys. Rev. Lett. 123, 123601 (2019).
  * (17) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-Hermitian physics and PT symmetry, Nat. Phys. 14, 11 (2018).
  * (18) H. Shen, B. Zhen, and L. Fu, Topological Band Theory for Non-Hermitian Hamiltonians, Phys. Rev. Lett. 120, 146402 (2018).
  * (19) N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Topological Origin of Non-Hermitian Skin Effects, Phys. Rev. Lett. 124, 086801 (2020).
  * (20) K. Zhang, Z. Yang, and C. Fang, Correspondence between Winding Numbers and Skin Modes in Non-Hermitian Systems, Phys. Rev. Lett. 125, 126402 (2020).
  * (21) S. Yao and Z. Wang, Edge states and topological invariants of non-Hermitian systems, Phys. Rev. Lett. 121, 086803 (2018).
  * (22) S. Yao, F. Song, and Z. Wang, Non-Hermitian chern bands, Phys. Rev. Lett. 121, 136802 (2018).
  * (23) J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect, Direct observation of Anderson localization of matter waves in a controlled disorder, Nature (London) 453, 891 (2008).
  * (24) G. Roati, C. D. Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Nature (London) 453, 895 (2008).
  * (25) S. Aubry and G. André, Analyticity breaking and Anderson localization in incommensurate lattices, Ann. Isr. Phys. Soc. 3, 133 (1980).
  * (26) P. G. Harper, The General Motion of Conduction Electrons in a Uniform Magnetic Field, with Application to the Diamagnetism of Metals, Proc. Phys. Soc. London Sect. A 68, 874 (1955).
  * (27) F. Jendrzejewski, A. Bernard, K. Müller, P. Cheinet, V.Josse, M. Piraud, L. Pezzé, L. Sanchez-Palencia, A. Aspect, and P. Bouyer, Three-dimensional localization of ultracold atoms in an optical disordered potential, Nat. Phys. 8, 398 (2012).
  * (28) H. Yao, H. Khouldi, L. Bresque, and L. Sanchez-Palencia, Critical behavior and fractality in shallow one-dimensional quasiperiodic potentials, Phys. Rev. Lett. 123, 070405 (2019).
  * (29) H. Yao, T. Giamarchi, and L. Sanchez-Palencia, Lieb-Liniger bosons in a shallow quasiperiodic potential: Bose glass phase and fractal Mott Lobes, Phys. Rev. Lett. 125, 060401 (2020).
  * (30) L. Sanchez-Palencia, Ultracold gases: At the edge of mobility, Nat. Phys. 11, 525 (2015).
  * (31) J. Biddle and S. Das. Sarma, Predicted mobility edges in one-dimensional incommensurate optical lattices: An exactly solvable model of Anderson localization, Phys. Rev. Lett. 104, 070601 (2010).
  * (32) S. Ganeshan, J. H. Pixley, and S. Das Sarma, Nearest neighbor tight binding models with an exact mobility edge in one dimension, Phys. Rev. Lett. 114, 146601 (2015).
  * (33) Y. Wang, X. Xia, L. Zhang, H. Yao, S. Chen, J. You, Q. Zhou, and X.-J. Liu, One dimensional quasiperiodic mosaic lattice with exact mobility edges, Phys. Rev. Lett. 125, 196604 (2020).
  * (34) T. Liu, H. Guo, Y. Pu, and S. Longhi, Generalized Aubry-André self-duality and Mobility edges in non-Hermitian quasi-periodic lattices, Phys. Rev. B 102, 024205 (2020).
  * (35) S. Longhi, Topological phase transition in non-Hermitian quasicrystals, Phys. Rev. Lett. 122, 237601 (2019).
  * (36) S. Longhi, Metal-insulator phase transition in a non-Hermitian Aubry-André-Harper Model, Phys. Rev. B 100, 125157 (2019).
  * (37) Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Topological Phases of Non-Hermitian Systems, Phys. Rev. X 8, 031079 (2018).
  * (38) K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-Hermitian physics, Phys. Rev. X 9, 041015 (2019).
  * (39) N. Hatano and D. R. Nelson, Localization Transitions in Non-Hermitian Quantum Mechanics, Phys. Rev. Lett. 77, 570 (1996).
  * (40) N. Hatano and D. R. Nelson, Vortex pinning and non-Hermitian quantum mechanics, Phys. Rev. B 56, 8651 (1997).
  * (41) N. Hatano and D. R. Nelson, Non-Hermitian Delocalization and Eigenfunctions, Phys. Rev. B 58, 8384 (1998).
  * (42) H. Jiang, L.-J. Lang, C. Yang, S.-L. Zhu, and S. Chen, Interplay of non-Hermitian skin effects and Anderson localization in non-reciprocal quasi-periodic lattices, Phys. Rev. B 100, 054301 (2019).
  * (43) Y. Liu, X.-P. Jiang, J. Cao, and S. Chen, Non-Hermitian mobility edges in one-dimensional quasicrystals with parity-time symmetry, Phys. Rev. B 101, 174205 (2020).
  * (44) Q.-B. Zeng, Y. B. Yang, and Y. Xu, Topological phases in non-Hermitian Aubry-André-Harper models, Phys. Rev. B 101, 020201 (2020).
  * (45) R. Hamazaki, K. Kawabata, and M. Ueda, Non-Hermitian Many-Body Localization, Phys. Rev. Lett. 123, 090603 (2019).
  * (46) T. Liu, S. Cheng, H. Guo, and G. Xianlong, Fate of Majorana zero modes, exact location of critical states, and unconventional real-complex transition in non-Hermitian quasiperiodic lattices, Phys. Rev. B 103, 104203 (2021).
  * (47) D. J. Thouless, A relation between the density of states and range of localization for one dimensional random systems, J. Phys. C: Solid State Phys. 5, 77 (1973).
  * (48) P. Sarnak, Spectral behavior of quasi periodic potentials, Commun. Math. Phys. 84, 377 (1982).
  * (49) S. Longhi, Phase transitions in a non-Hermitian Aubry-André-Harper model, Phys. Rev. B 103 054203 (2021).
  * (50) Y. Liu, Y. Wang, X.-J. Liu, Q. Zhou, and S. Chen, Exact mobility edges, PT-symmetry breaking and skin effect in one-dimensional non-Hermitian quasicrystals, Phys. Rev. B 103, 014203 (2021).
  * (51) H. L. Cycon, R. G. Froese, W. Kirsch, and B. Simon, Schrödinger operators: With application to quantum mechanics and global geometry[M]. Springer, 2009.
  * (52) A. Schreiber, K. N. Cassemiro, V. Potocek, A. Gabris, I. Jex, and Ch. Silberhorn, Decoherence and disorder in quantumwalks: From ballistic spread to localization, Phys. Rev. Lett. 106, 180403 (2011).
  * (53) I. D. Vatnik, A. Tikan, G. Onishchukov, D. V. Churkin, and A. A. Sukhorukov, Anderson localization in synthetic photonic lattices, Sci. Rep. 7, 4301 (2017).
  * (54) S. Derevyanko, Disorder-aided pulse stabilization in dissipative synthetic photonic lattices, Sci. Rep. 9, 12883 (2019).


# Soft gluon fields and anomalous magnetic moment of muon

Sergei Nedelko<EMAIL_ADDRESS>Aleksei Nikolskii<EMAIL_ADDRESS>Vladimir Voronin<EMAIL_ADDRESS>Bogoliubov Laboratory of Theoretical
Physics, JINR, 141980 Dubna, Russia

###### Abstract

An impact of nonperturbatively treated soft gluon modes on the value of
anomalous magnetic moment of muon $a_{\mu}$ is studied within the mean-field
approach to QCD vacuum and hadronization. It is shown that radial excitations
of vector mesons strongly enhance contribution of hadronic vacuum polarization
to $a_{\mu}$, doubling the contribution of one-meson processes compared to the
result for ground state mesons. The mean field also strongly influences the
hadronic light-by-light scattering contribution due to the Wilson line in
quark propagators.

##  I Introduction

The measurement of the anomalous magnetic moment of muon
$a_{\mu}=(g_{\mu}-2)/2$ conducted by Muon $g-2$ Collaboration at FNAL [1]

$a_{\mu}^{\text{FNAL}}=116592040(54)\times 10^{-11}$ (1)

retained the earlier discrepancy with the Standard Model prediction (see
reviews [2, 3])

$a_{\mu}^{\text{SM}}=116591810(43)\times 10^{-11}.$ (2)

Significant fraction of the uncertainty in the SM value comes from the part
described by QCD and hadrons. Since robust analytical first-principle QCD
evaluations at the whole range of energies are not feasible at the moment,
hadronic contributions to the “consensus value” (2) are essentially extracted
from experimental data via dispersion relations supplemented by perturbative
QCD where large momentum transfer is involved. In order to gain a clearer
insight into relevant low-energy dynamics, hadronic contributions to $a_{\mu}$
were also extensively studied with the methods of lattice QCD, Dyson-Schwinger
equations, within the models of QCD vacuum and hadronization, and in the
models of physics beyond the Standard Model.

In this paper, we address the specific role which soft gluon fields may play
in the hadronic vacuum polarization (HVP) and light-by-light scattering (HLbL)
contributions to $a_{\mu}$ using the mean-field approach to physical QCD
vacuum. Specifically, we will touch this issue by calculating the lowest order
in the number of intermediate mesons contributions (quark loop and one
intermediate meson) within the domain model of QCD vacuum and hadronization
developed in Refs. [4, 5, 6, 7, 8, 9, 10]. The driving feature of the model is
the vacuum mean field represented by domain-structured configurations of
almost everywhere homogeneous Abelian (anti-)self-dual gluon field which is
treated nonperturbatively. These field configurations are characterized by
various vacuum expectation values and can be interpreted as “soft vacuum
fields” that reside in condensates of local operators in the context of
operator product expansion. Such a mean-field approach has simultaneously
provided for static and dynamic confinement – area law for Wilson loop and
absence of poles in the propagators of color charged fields in complex
momentum plane, respectively, chiral symmetry breaking and resolution of
$U_{A}(1)$ problem. Upon hadronization, the spectrum of light, heavy-light
mesons and heavy quarkonia, including excited states, decay constants and form
factors calculated with a minimal set of parameters systematically agree very
well with experimental data. The detailed description of various features of
the approach can be found in above-mentioned articles. The mean field does not
change short-distance behavior of quark, gluon and ghost propagators, but
otherwise strongly modifies propagators, in a way consistent with the results
of functional renormalization group and Lattice QCD [7].

For the present considerations, it is important that effective meson action
describes explicitly the mesons with all possible quantum numbers and their
interactions. In particular, it contains not only ground state mesons but also
their excited states. The absence of poles in the momentum representation of
the quark and gluon propagators in the finite complex momentum plane, i.e.
confinement of dynamical color, leads to the Regge spectrum of light meson
masses. This allows one to compute contribution of excited meson states to the
HVP and HLbL straightforwardly.

Another important feature of propagators in the presence of the mean gluon
field is the Schwinger gauge phase factor also known as Wilson line

$\mathcal{P}\,\exp\left[ig\int_{x}^{y}dz_{\alpha}\hat{B}_{\alpha}(z)\right],$
(3)

which provides explicit invariance of the effective meson action with respect
to the background gauge transformations. Calculation of pion transition form
factor $F_{\pi\gamma^{*}\gamma}$ within the model has demonstrated [8] that
the Wilson lines are responsible for deviation of $F_{\pi\gamma^{*}\gamma}$
from Brodsky-Lepage limit [11] at large photon virtuality. The Wilson lines
were also critically important for the description of decays of vector mesons
into a couple of pseudoscalar mesons simultaneously with their masses. It is
interesting to study influence of the mean field (via confinement and Wilson
lines) on the anomalous magnetic moment of muon $a_{\mu}$. It has to be
stressed that the same values of parameters that were used for calculations of
the masses and decay constants – are employed in the present calculations. We
consider contributions to $a_{\mu}$ schematically shown in Fig. 1, and
concentrate on one-loop quark and one-meson diagrams where the role of soft
mean field is seen manifestly.

There are two main observations. Radial excitations of vector mesons strongly
enhance the contribution of hadronic vacuum polarization to $a_{\mu}$. In one-
meson processes the total contribution of exited states is as big as the
ground-state meson contribution. The mean field strongly influences the
hadronic light-by-light scattering contribution due to the Wilson line in
quark propagators.

Figure 1: Hadronic contributions to the muon $g-2$. The first is known as
hadronic vacuum polarization (HVP), the second is hadronic light-by-light
scattering (HLbL). Shaded regions denote all possible contributions of strong
interactions.

The paper is organized as follows. Section II contains a brief review of the
model employed in this study. Section III is devoted to HVP contribution to
the anomalous magnetic moment of muon $a_{\mu}$. One-loop quark HLbL-type
contribution to $a_{\mu}$ is considered in Section IV. The results are
discussed in Section V.

##  II Description of the model

Motivation of the model employed in the present paper and technical details of
its derivation can be found in Refs. [4, 5, 6, 7, 8, 9, 10]. The defining
property of the model is domain-structured Abelian (anti-)self-dual vacuum
gluon fields which are treated nonperturbatively. The hadronization procedure
of bilocal quark currents in the presence of such vacuum fields leads to the
following effective nonlocal meson action (in Euclidean space):

$\displaystyle Z$ $\displaystyle={\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^{4}x\phi^{2}_{\cal
Q}(x)-\sum\limits_{k=2}^{\infty}\frac{1}{k}W_{k}[\phi]\right\\},$ (4)
$\displaystyle W_{k}[\phi]$ $\displaystyle=\sum\limits_{{\cal Q}_{1}\dots\cal
Q_{k}}h_{{\cal Q}_{1}}\dots h_{{\cal Q}_{k}}\int d^{4}x_{1}\dots\int
d^{4}x_{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}),$
(5) $\displaystyle\Phi_{{\cal Q}}(x)$
$\displaystyle=\int\frac{d^{4}p}{(2\pi)^{4}}e^{ipx}{\mathcal{O}}_{{\mathcal{Q}}{\mathcal{Q}}^{\prime}}(p)\tilde{\phi}_{{\mathcal{Q}}^{\prime}}(p),\quad
C_{\mathcal{Q}}=C_{J},\quad C^{2}_{S/P}=2C^{2}_{V/A}=\frac{1}{9}.$

The condensed index $\mathcal{Q}\equiv\\{aJLn\\}$ denotes all quantum numbers
of a meson (isospin, spin-parity, orbital and radial numbers, space-time
indices for the fields with nonzero total angular momentum), and $\Lambda$ is
related to the strength of the background gluon field. The orthogonal
transformation ${\mathcal{O}}_{{\mathcal{Q}}{\mathcal{Q}}^{\prime}}$
diagonalizes quadratic part of the action and thus relates physical fields
$\phi_{{\mathcal{Q}}}$ to auxiliary fields $\Phi_{{\mathcal{Q}}}$ arising by
virtue of the hadronization procedure. Meson masses are determined by the
poles of the meson propagator defined from the equation

$\displaystyle 0$
$\displaystyle=\frac{\Lambda^{2}}{g^{2}C^{2}_{\mathcal{Q}}}+\tilde{\Gamma}^{(2)}_{\cal
Q}(-M^{2}_{\cal Q}),$ (6)

solved with respect to the mass $M^{2}_{\cal Q}$ (location of the pole).
Constant $h_{\mathcal{Q}}$ is defined by equation

$\displaystyle 1$ $\displaystyle=h^{2}_{\cal
Q}\left.\frac{d}{dp^{2}}\tilde{\Gamma}^{(2)}_{\cal
Q}(p^{2})\right|_{p^{2}=-M^{2}_{\cal Q}},$ (7)

solved with respect to constant $h_{\mathcal{Q}}$ to make the residue at the
pole of meson propagator to be equal to unity. The function
$\tilde{\Gamma}^{(2)}_{\cal Q}$ is two-point correlation function
$\tilde{\Gamma}^{(2)}_{\cal QQ^{\prime}}(p)$ diagonalized with respect to all
quantum numbers:

$\tilde{\phi}_{\mathcal{Q}}(-p)\left[\mathcal{O}^{T}(p)\tilde{\Gamma}^{(2)}(p^{2})\mathcal{O}(p)\right]_{\cal
QQ^{\prime}}\tilde{\phi}_{{\mathcal{Q}}^{\prime}}(p)=\tilde{\Gamma}^{(2)}_{\cal
Q}(p^{2})\tilde{\phi}_{\mathcal{Q}}(-p)\tilde{\phi}_{\mathcal{Q}}(p).$ (8)

Two-point nonlocal vertex function $\tilde{\Gamma}^{(2)}_{\cal
QQ^{\prime}}(p)$ is given by

$\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}}},$ (9)

and analogous formulas can be written for $k$-point vertices
$\Gamma^{(k)}_{\mathcal{Q}_{1}\dots\mathcal{Q}_{k}}$. Quark loops
$G^{(k)}_{{\cal Q}_{1}\dots\cal Q_{k}}$ are averaged over the background field
with measure $d\sigma_{B}$:

$\displaystyle\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),$ (10)
$\displaystyle\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})}=$ $\displaystyle\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$
$\displaystyle{\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\\}.$ (11)

Here $S(x,y)$ is the quark propagator and $V_{{\cal Q}}$ are nonlocal meson-
quark-antiquark vertices with ${\cal Q}$ \- meson quantum numbers (see [8] for
details).

In order to take into account the basic properties of the mean field
corresponding to almost everywhere homogeneous (anti-)self-dual Abelian gauge
field represented by the statistical ensemble of domain wall networks and
still be able to perform calculations analytically, we make the following
simplifications. The quark propagator and meson-quark vertices are
approximated by those in the homogeneous (anti-)self-dual Abelian background
field, and the statistical properties of the ensemble of domain networks are
taken into account by means of averaging of the quark loops over all
configurations of the homogeneous background fields as well as via correlators
$\Xi_{n}$ in the background field ensemble.

Specifically, the quark loops $G^{(k)}$ are averaged over parity-conjugated
self-dual and anti-self-dual Abelian (anti-)self-dual background field and
direction (in coordinate and color spaces). Averaging over spatial directions
is performed with the help of generating formula

$\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)}},$ (12)

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$:

$\displaystyle\hat{B}_{\mu}=-\frac{1}{2}\hat{n}B_{\mu\nu}x_{\nu},\
\hat{n}=t^{3}\cos\xi+t^{8}\sin\xi,$
$\displaystyle\tilde{B}_{\mu\nu}=\frac{1}{2}\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},$
(13) $\displaystyle
f_{\alpha\beta}=\frac{\hat{n}}{2\upsilon\Lambda^{2}}B_{\alpha\beta},\
\upsilon=\mathrm{diag}\left(\frac{1}{6},\frac{1}{6},\frac{1}{3}\right),\
f_{\mu\alpha}f_{\nu\alpha}=\delta_{\mu\nu},$

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

$\displaystyle V^{aJln}_{\mu_{1}\dots\mu_{l}}={\cal
C}_{ln}\mathcal{M}^{a}\Gamma^{J}F_{nl}\left(\frac{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal
D}}^{2}\\!\\!\\!(x)}{\Lambda^{2}}\right)T^{(l)}_{\mu_{1}\dots\mu_{l}}\left(\frac{1}{i}\frac{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal
D}}\\!(x)}{\Lambda}\right),$ (14) $\displaystyle{\cal
C}^{2}_{ln}=\frac{l+1}{2^{l}n!(n+l)!},\quad
F_{nl}(s)=s^{n}\int_{0}^{1}dtt^{n+l}\exp(st),$
$\displaystyle{\stackrel{{\scriptstyle\leftrightarrow}}{{\mathcal{D}}}}\vphantom{D}^{ff^{\prime}}_{\mu}=\xi_{f}\stackrel{{\scriptstyle\leftarrow}}{{\mathcal{D}}}_{\mu}-\
\xi_{f^{\prime}}\stackrel{{\scriptstyle\rightarrow}}{{\mathcal{D}}}_{\mu},\ \
\stackrel{{\scriptstyle\leftarrow}}{{\mathcal{D}}}_{\mu}(x)=\stackrel{{\scriptstyle\leftarrow}}{{\partial}}_{\mu}+\
i\hat{B}_{\mu}(x),\ \
\stackrel{{\scriptstyle\rightarrow}}{{\mathcal{D}}}_{\mu}(x)=\stackrel{{\scriptstyle\rightarrow}}{{\partial}}_{\mu}-\
i\hat{B}_{\mu}(x),$
$\displaystyle\xi_{f}=\frac{m_{f^{\prime}}}{m_{f}+m_{f^{\prime}}},\
\xi_{f^{\prime}}=\frac{m_{f}}{m_{f}+m_{f^{\prime}}}.$

Here $\mathcal{M}^{a}$ and $\Gamma^{J}$ are flavor and Dirac matrices
corresponding to a given meson field, $\xi_{f},\xi_{f^{\prime}}$ 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

$\displaystyle S_{f}(x,y)$
$\displaystyle=\exp\left(-\frac{i}{2}\hat{n}x_{\mu}B_{\mu\nu}y_{\nu}\right)H_{f}(x-y),$
(15) $\displaystyle\tilde{H}_{f}(p)$
$\displaystyle=\frac{1}{2\upsilon\Lambda^{2}}\int_{0}^{1}dse^{(-p^{2}/2\upsilon\Lambda^{2})s}\left(\frac{1-s}{1+s}\right)^{m_{f}^{2}/4\upsilon\Lambda^{2}}$
$\displaystyle\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],$

where anti-Hermitean representation of Dirac matrices is used, and “$\pm$”
signs are arranged in accordance with formula (13). Due to the Wilson line
translation invariance is hold only up to the large gauge transformation

$\displaystyle S_{f}(x+a,y+a)=U(x)S_{f}(x,y)U^{+}(y),\ \
U(x)=\exp\left(-\frac{i}{2}\hat{n}x_{\mu}B_{\mu\nu}a_{\nu}\right).$ (16)

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 values of parameters of the model given in Table 1 were fitted to the
masses of ground-state mesons $\pi,\rho,K,K^{*},J/\psi,\Upsilon,\eta^{\prime}$
(see [7]), and the same values are employed in the present paper.

$m_{u/d}$(MeV) |  $m_{s}$(MeV) |  $m_{c}$(MeV) |  $m_{b}$(MeV) |  $\Lambda$(MeV) | $\alpha_{s}$
---|---|---|---|---|---
$145$ | $376$ | $1566$ | $4879$ | $416$ | $3.45$
Table 1: Values of parameters introduced in Ref. [7] and used for calculations
of mass spectrum, decay constants and form factors. The same values are
employed in the present paper.

Electromagnetic interactions are included in this scheme with the help of
expansions (see [12, 7])

$\displaystyle S_{f}(x,y|A)$
$\displaystyle=S_{f}(x,y)+\sum_{n=1}^{\infty}\left(Q_{f}e\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),$
(17) $\displaystyle V_{\mathcal{Q}}(x|A)$
$\displaystyle=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}}^{(n)}(x,z_{1},\dots,z_{n})A_{\mu_{1}}(z_{1})\cdots
A_{\mu_{n}}(z_{n}),$ (18)

where $Q$ is a diagonal matrix of quark charges in units of electron charge
$e$, and meson-photon vertices appear due to nonlocality. In the next section
we give explicit examples of such interactions. Effective meson action takes
the form

$\displaystyle Z$ $\displaystyle={\cal N}\int D\phi_{\cal Q}\int
DA_{\mu}\exp\left\\{-\frac{1}{4}\int d^{4}x\
F_{\mu\nu}F_{\mu\nu}-\frac{\Lambda^{2}}{2}\frac{h^{2}_{\cal
Q}}{g^{2}C^{2}_{\mathcal{Q}}}\int d^{4}x\phi^{2}_{\cal
Q}(x)-\Gamma(A)-\sum\limits_{k=2}^{\infty}\frac{1}{k}W_{k}[\phi|A]\right\\},$
(19) $\displaystyle W_{k}[\phi|A]$ $\displaystyle=\sum\limits_{{\cal
Q}_{1}\dots\cal Q_{k}}h_{{\cal Q}_{1}}\dots h_{{\cal Q}_{k}}\int
d^{4}x_{1}\dots\int d^{4}x_{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),$ (20)

where

$\Gamma(A)=\int
d\sigma_{B}\mathrm{Tr}\log\left[1+Qe\gamma_{\mu}A_{\mu}(x)S(x,y)\right],$

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})$
via substitutions

$S_{f}(x,y)\to S_{f}(x,y|A),\quad V_{\mathcal{Q}}(x)\to V_{\mathcal{Q}}(x|A).$

##  III Hadronic vacuum polarization contribution to muon $g-2$

We employ the following representation [13, 14] for the contribution to
$a_{\mu}$ from the first diagram in Fig. 1

$\displaystyle a^{\text{HVP}}_{\mu}$
$\displaystyle=4\pi^{2}\left(\frac{\alpha}{\pi}\right)^{2}\int\limits_{0}^{\infty}dp^{2}f(p^{2},m_{\mu}^{2})\tilde{\Pi}^{\rm
R}(p^{2}),$ (21) $\displaystyle f(p^{2},m_{\mu}^{2})$
$\displaystyle=\frac{m_{\mu}^{2}p^{2}Z^{3}(1-p^{2}Z)}{1+m_{\mu}^{2}Z^{2}p^{2}},$
(22) $\displaystyle Z$
$\displaystyle=-\left(p^{2}-\sqrt{p^{4}+4m_{\mu}^{2}p^{2}}\right)/2m_{\mu}^{2}p^{2},$

where $m_{\mu}$ is the muon mass, and $\Pi^{\rm R}(p^{2})$ is related to
hadronic vacuum polarization

$\tilde{\Pi}^{\rm
R}_{\mu\nu}(p)=(p^{2}\delta_{\mu\nu}-p_{\mu}p_{\nu})\tilde{\Pi}^{\rm
R}(p^{2})$ (23)

renormalized at $p^{2}=0$. HVP is defined as

$\tilde{\Pi}_{\mu\nu}(p)=\int d^{4}xe^{-ipx}\langle
j_{\mu}(x)j_{\nu}(0)\rangle,$

where $j_{\mu}=\bar{\psi}Q\gamma_{\mu}\psi$ are quark currents.

###  III.1 Quark loop

Figure 2: Quark loop contributions to HVP. Solid lines denote quark propagator
(15). Gradient shading indicates averaging the quark loop over configurations
of the mean field.

One-loop quark contribution to HVP corresponds to the diagram in Fig. 2 which
is given by

$\tilde{\Pi}_{\mu\nu}^{(1)}(p)=\sum_{f}Q_{f}^{2}\ (-1)\ {\rm
reg}\int\frac{d^{4}l}{(2\pi)^{4}}\mathrm{Tr}\
\gamma_{\mu}\tilde{H}_{f}(l+p/2)\gamma_{\nu}\tilde{H}_{f}(l-p/2),$ (24)

where reg denotes some UV regularization which preserves gauge invariance. In
this case with two local electromagnetic vertices, translation-noninvariant
Wilson lines in quark propagators cancel each other. Performing the loop
integration and renormalization at $p^{2}=0$ (see Appendix B), we find

$\displaystyle\tilde{\Pi}^{(1)\mathrm{R}}_{\mu\nu}(p)$
$\displaystyle=(p^{2}\delta_{\mu\nu}-p_{\mu}p_{\nu})\sum_{f}Q_{f}^{2}\tilde{\Pi}_{f}^{\rm
R}(p^{2}),$ $\displaystyle\tilde{\Pi}_{f}^{(1)\mathrm{R}}(p^{2})$
$\displaystyle=-\mathrm{Tr}\frac{1}{2\pi^{2}}\int_{0}^{1}ds_{1}\int_{0}^{1}ds_{2}\frac{s_{1}s_{2}+s_{1}^{2}s_{2}^{2}/3}{(s_{1}+s_{2})^{4}}\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]^{\frac{m_{f}^{2}}{4v\Lambda^{2}}}\left[\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}\right)-1\right].$

The above integral converges for any finite complex value of $p^{2}$, and its
derivative is continuous, hence $\tilde{\Pi}_{f}^{(1)\mathrm{R}}(p^{2})$ is an
analytical function of momentum as was mentioned in the Introduction. The
expression for $\tilde{\Pi}_{\mu\nu}^{(1)\mathrm{R}}(p^{2})$ is substituted
into Eq. (21) to obtain numbers given in Table 2.

| $u$ | $d$ | $s$ | $c$ | $b$ | total one-loop | N$\chi$QM
---|---|---|---|---|---|---|---
$a_{\mu}^{\text{HVP,one-loop}}$ | 284.7 | 71.2 | 26.4 | 7.2 | 0.2 | 389.6 | 533
Table 2: The contributions to $a_{\mu}$ in units of $10^{-10}$ from one-loop
quark diagram shown in Fig. 2. One-loop contribution in nonlocal chiral quark
model (N$\chi$QM) [15] is given for comparison.

###  III.2 Intermediate mesons

Figure 3: Contributions to HVP of the first order with respect to the number
of intermediate mesons. Notation is the same as in Fig. 2.

Expression for the diagrams in Fig. 3 is derived systematically from formula
(4) (see Appendix A) and is given by

$\tilde{\Pi}_{\mu\nu}^{(2)}(p)=\sum_{\mathcal{Q}}\tilde{\Gamma}_{\mathcal{Q};\mu}(p)\left[\frac{\Lambda^{2}h_{\mathcal{Q}}^{2}}{g^{2}C^{2}_{\mathcal{Q}}}+h_{\mathcal{Q}}^{2}\tilde{\Gamma}^{(2)}_{\cal
Q}(p^{2})\right]^{-1}\tilde{\Gamma}_{\mathcal{Q};\nu}(p),$ (25)

where the term

$\left[\frac{\Lambda^{2}h_{\mathcal{Q}}^{2}}{g^{2}C^{2}_{\mathcal{Q}}}+h_{\mathcal{Q}}^{2}\tilde{\Gamma}^{(2)}_{\cal
Q}(p^{2})\right]^{-1}$ (26)

is the propagator of a meson field (cf. Eqs. (6) and (7)). Only the meson
states with quantum numbers of a photon (ground-state and excited meson
fields) contribute to Eq. (25), while the contributions of other states are
identically zero. The meson propagator includes $\Gamma_{\mathcal{Q}}^{(2)}$
which is given by formula (8), and only connected piece contributes to
$\Gamma^{(2)}_{{\cal Q}_{1}{\cal Q}_{2}}$ in the case of intermediate vector
mesons:

$\Gamma^{(2)}_{{\cal V}_{1}{\cal V}_{2}}=\overline{G^{(2)}_{{\cal V}_{1}{\cal
V}_{2}}(x_{1},x_{2})}.$

In one-loop approximation, meson-photon transition amplitudes
$\tilde{\Gamma}_{\mathcal{Q};\mu}(p)$ are given by the sum of two diagrams
shown in Fig. 4 (first and second terms, correspondingly)

$\displaystyle\tilde{\Gamma}_{\mathcal{Q};\mu}(p)$
$\displaystyle=h_{\mathcal{Q}}\sum_{\mathcal{Q}}O_{\mathcal{Q}^{\prime}\mathcal{Q}}(p^{2})\int
d\sigma_{B}\int\frac{d^{4}p^{\prime}}{2\pi^{4}}\int d^{4}xe^{ipx}\int
d^{4}ye^{ip^{\prime}y}\mathrm{Tr}S(y,x)V_{\mathcal{Q}^{\prime}}(x)S(x,y)Q\gamma_{\mu}$
$\displaystyle\quad+h_{\mathcal{Q}}\sum_{\mathcal{Q}^{\prime}}O_{\mathcal{Q}^{\prime}\mathcal{Q}}(p^{2})\int
d\sigma_{B}\int d^{4}z\ e^{ipz}\
\mathrm{Tr}S(x,x)V_{\mathcal{Q}^{\prime};\mu}^{(1)}(x,z)\gamma_{\mu}.$

Figure 4:  One-loop diagrams that contribute to meson-photon transition
amplitude $\tilde{\Gamma}_{\mathcal{Q}\mu}(p)$. Notation is the same as in
Fig. 2.

Due to nonlocality of meson-quark vertices
$\tilde{\Gamma}_{\mathcal{Q}}^{(2)}$ and $\tilde{\Gamma}_{\mathcal{Q};\mu}(p)$
are UV-finite, and gauge invariance provides that

$p_{\mu}\tilde{\Gamma}_{\mathcal{Q};\mu}(p)=0.$

Therefore, only transversal parts of meson-photon transition amplitudes and
meson propagator contribute to Eq. (25), and

$\tilde{\Pi}_{\mu\nu}^{(2)}(p)=(p^{2}\delta_{\mu\nu}-p_{\mu}p_{\nu})\tilde{\Pi}^{(2)}(p^{2}).$
(27)

The explicit formulas for $\tilde{\Gamma}_{\mathcal{Q}}^{(2)}$ and
$\tilde{\Gamma}_{\mathcal{Q}\mu}$ can be found in Ref. [7], where they were
used to find spectrum and transition constants of various mesons. Substituting
the formulas into (25) and (21), we find the contributions to
$a_{\mu}^{\text{HVP}}$ corresponding to various vector meson fields and their
radial excitations ($n=0\dots 6$) which are given in Table 3.

| $V_{0}$ | $V_{1}$ | $V_{2}$ | $V_{3}$ | $V_{4}$ | $V_{5}$ | $V_{6}$ | total | N$\chi$QM
---|---|---|---|---|---|---|---|---|---
$\rho_{n}$ | $41.5$ | $27$ | $14.8$ | $9.5$ | $6$ | $3.3$ | $3.5$ | $105.5$ |
$\omega_{n}$ | $4.6$ | $3$ | $1.6$ | $1.1$ | $0.7$ | $0.4$ | $0.4$ | $11.7$ |
$\phi_{n}$ | $2.8$ | $2$ | $1.4$ | $1$ | $0.7$ | $0.5$ | $0.4$ | $8.8$ |
$J/\psi_{n}$ | $0.2$ | $0.19$ | $0.18$ | $0.17$ | $0.16$ | $0.15$ | $0.14$ | $1.2$ |
total | $49.1$ | $32.2$ | $17.9$ | $11.7$ | $7.6$ | $4.4$ | $4.4$ | $127.3$ | $13$
Table 3: The contributions to $a_{\mu}$ in units of $10^{-10}$ corresponding
to various vector mesons and their radial excitations $V_{n}$. Currently, we
take into account seven radial states $n$, and numerical values of
contributions $V_{6}$ are affected by this finite number of retained radial
states most notably. Contribution of ground state vector mesons in N$\chi$QM
[15] is shown for comparison. It can be seen that the total contribution of
radial excitations is as large as the ground state itself.

Total contribution to $a_{\mu}^{\text{HVP}}$ from intermediate vector mesons
is numerically smaller than the one of quark loop given in Table 2. Total
contribution of radially excited mesons is as large as the contribution of the
ground state mesons.

this paper | N$\chi$QM | ENJL | LMD | GNC | data-driven | lattice QCD
---|---|---|---|---|---|---
516.9 | 623(40) | 750 | 570(170) | 630(50) | 693.1(4.0) | 711.6(18.4) | 707.5(5.5)
Table 4:  Contributions to $a_{\mu}^{\text{HVP}}$ found in the present paper,
nonlocal chiral quark model (N$\chi$QM) [15], extended Nambu–Jona-Lasinio
model (ENJL) [16], lowest meson dominance (LMD) approximation to large-Nc QCD
[17, 18], gauged nonlocal constituent quark model (GNC) [19]. The total
leading-order HVP contribution to $a_{\mu}$ extracted from experimental data
[2, 20, 21, 22, 23, 24, 25] is shown in column “data-driven”. Average of the
lattice QCD results by many authors [26, 27, 28, 29, 30, 31, 32, 33, 34] for
leading-order HVP contribution collected in Ref. [2] is given in column
“lattice QCD” on the left. Central value of lattice QCD results given in Ref.
[35] is shown on the right.

The contribution to $a_{\mu}$ corresponding to quark loop (24) and
intermediate vector mesons (25), in the one-meson approximation, is given in
Table 4. The results of various other calculations are shown for the
reference. It has to be stressed that data-driven and LQCD results take into
account much more various contributions than the model calculations, the
(N$\chi$QM) result accounts also for meson loop with two intermediate mesons
which has been found to be three times bigger than the one-meson part [15]. It
is seen that the total contribution of the quark loop and one-meson
approximation underestimates the value of HVP in the present calculation. One
concludes that the contribution of two intermediate mesons has to be taken
into account. The largest of them are expected to correspond to two
intermediate pions and kaons. In the present approach, calculation of two-
meson diagrams is straightforward but somewhat bulky. Investigation of this
and several other potentially important contributions deserves a separate
systematic study and will be reported elsewhere.

##  IV Hadronic light-by-light contribution

We employ the well-known projection technique and angular averages [36, 37,
38]

$\displaystyle
a^{\text{HLbL}}_{\mu}=-\frac{1}{48m}\mathrm{Tr}\left\\{(\not{p}-m)[\gamma_{\rho},\gamma_{\sigma}](\not{p}-m)\Gamma_{\rho\sigma}^{\text{HLbL}}\right\\},$
(28)
$\displaystyle\begin{split}\Gamma_{\rho\sigma}^{\text{HLbL}}(p)=\left.\frac{\partial}{\partial
k_{\sigma}}\Gamma_{\rho}^{\text{HLbL}}(p,p^{\prime})\right|_{k=0}&=e^{6}\int\frac{d^{4}q_{1}}{(2\pi)^{4}}\frac{d^{4}q_{2}}{(2\pi)^{4}}\frac{1}{q_{1}^{2}q_{2}^{2}(-q_{1}-q_{2})^{2}}\\\
&\quad\times\gamma_{\mu}\frac{\not{p}^{\prime}+\not{q}_{1}-m}{(p+q_{1})^{2}+m^{2}}\gamma_{\lambda}\frac{\not{p}-\not{q}_{2}-m}{(p-q_{2})^{2}+m^{2}}\gamma_{\nu}\left.\frac{\partial}{\partial
k_{\rho}}\Pi_{\mu\nu\lambda\sigma}(q_{1},q_{2},k-q_{1}-q_{2})\right|_{k=0}\end{split}$
(29)

to extract hadronic light-by-light scattering contribution to $a_{\mu}$. Here
$\Gamma_{\rho}^{\text{HLbL}}$ corresponds to the second diagram in Fig. 1
(without external legs), $k$ is momentum of outgoing photon, $p,p^{\prime}$
are momenta of muon in initial and final states. Hadronic light-by-light
tensor is defined as

$\Pi_{\mu\nu\lambda\sigma}(q_{1},q_{2},q_{3})=-\int d^{4}x\ d^{4}y\
d^{4}ze^{-iq_{1}x-iq_{2}y-iq_{3}z}\langle
0|j_{\mu}(x)j_{\nu}(y)j_{\lambda}(z)j_{\sigma}(0)|0\rangle,$ (30)

where $j_{\mu}=\bar{\psi}Q\gamma_{\mu}\psi$ are quark currents.

In the one-loop approximation HLbL tensor is given by the sum of six diagrams,
five of which differ from the one shown in Fig. 5 by permutations of photon
legs:

$\begin{split}\Pi_{\mu\nu\lambda\sigma}(q_{1},q_{2},q_{3})&=-\sum_{f}Q_{f}^{4}\int
d\sigma_{B}\left[\int d^{4}x\ d^{4}y\
d^{4}ze^{-iq_{1}x-iq_{2}y-iq_{3}z}\right.\\\
&\quad\left.\times\text{Tr}\gamma_{\mu}S_{f}(x,y)\gamma_{\nu}S_{f}(y,z)\gamma_{\lambda}S_{f}(z,0)\gamma_{\sigma}S_{f}(0,x)\vphantom{\int}\right]+\
\text{five other permutations}.\end{split}$ (31)

Figure 5: One of the possible one-loop HLbL contributions, the other five
differ by permutation of photon legs. Among these six diagrams, only three
have to be calculated, the other three differ by reversed direction of the
quark loop (see also [36]).

Wilson lines in quark propagators do not cancel in this case, and momentum
does not conserve in each vertex. The sum of six diagrams is finite, but there
is logarithmic divergence in each diagram. HLbL tensor appears in Eq. (29) in
the derivative with respect to external photon momentum

$\left.\frac{\partial}{\partial
k_{\rho}}\Pi_{\mu\nu\lambda\sigma}(q_{1},q_{2},k-q_{1}-q_{2})\right|_{k=0}$
(32)

which makes all six terms finite [36], so it can be calculated without
regularization. The external field significantly increases the size of
formulas which become too cumbersome to be presented here, but the calculation
is straightforward. After evaluation of traces, integrations and averaging
over background field with the help of Eq. (12), the derivative of HLbL tensor
(32) is substituted into Eq. (29), which in turn is substituted into formula
(28). The latter is then averaged over direction of muon momentum in Euclidean
space as suggested in Ref. [37], and analytically continued to Minkowski space
$p^{2}\to-m_{\mu}^{2}$. The resulting integral over proper times of quark
propagators (15), virtualities of loop photons and the angle between photons
is computed numerically. The numbers are presented in Table 5, column “full”.

| without BF | without WL | full | N$\chi$QM | DSE | C$\chi$QM | ENJL | HLS | data-driven | QED${}_{\text{L}}$
---|---|---|---|---|---|---|---|---|---|---
$u$ | 16.8 | 11.7 | 7.77 |  |  |  |  |  |  |
$d$ | 1 | 0.73 | 0.49 |  |  |  |  |  |  |
$s$ | 0.2 | 0.2 | 0.18 |  |  |  |  |  |  |
$c$ | 0.2 | 0.2 | 0.2 |  |  |  |  |  |  |
total | 18.24 | 12.79 | 8.63 | 11(0.9) | 10.7(0.2) | 8.2 | 2.1(0.3) | 0.97(1.11) | 9.2(1.9) | 7.87(3.06)(1.77)
Table 5: Contributions from one-loop HLbL tensor within the model to $a_{\mu}$
in units of $10^{-10}$. Note that contributions of $s$ and $c$-quarks are
comparable due to charge factors despite the large difference in masses. For
comparison, we present quark loop HLbL contributions in nonlocal quark chiral
model (N$\chi$QM) [39], Dyson-Schwinger equations (DSE) [40], constituent
chiral quark model (C$\chi$QM) [41], extended Nambu–Jona-Lasinio model (ENJL)
[42] (result depends on cutoff, “central value” is shown), hidden local
symmetry approach (HLS) [43]. Data-driven estimate of all leading-order
contributions [2] is given in column “data-driven”. Result of evaluations of
both QED and QCD on finite lattice [44] is shown in column
“QED${}_{\text{L}}$” (the first number in parentheses is statistical error,
the second is systematic).

The analogous procedure is carried out for the case when Wilson lines in quark
propagators are simply omitted in order to estimate their role in the final
result. This is achieved by substitution (cf. Eq. (15))

$S_{f}(x,y)\to H_{f}(x-y)$

in Eq. (31). The numbers are given in column “without WL” in Table 5. For the
column “without BF”, background field is set to zero, and propagators
$S_{f}(x,y)$ transform into Dirac propagators with masses given in Table 1.
The contributions to $a_{\mu}$ in this case are found with the help of
formulas given in Ref. [45]. The results show that background field and Wilson
lines affect the dynamics of light quarks most notably.

##  V Discussion

As it has already been mentioned, rather unexpected result is that radial
excitations of vector mesons turned out far from being negligible. As it is
seen from Table 3 their total contribution is even bigger than ground state
meson’s part. Qualitative reason for that can be seen in highly nontrivial
dependence of the meson propagator defined by the meson polarization operator,
quark-meson vertices Eq. (14) and $V\to\gamma$ transition amplitudes Eq. (4).
Though the masses of radial excitations of vector mesons are much larger than
the mass of ground state mesons, the propagators in Euclidean domain are very
different from the propagator of free vector field, and contributions of the
excited states turn out to be not suppressed as much as it would be naively
expected (see Fig. 6).

Figure 6: Propagators of mesons (left panel) and corresponding contributions
to $\tilde{\Pi}^{(2)}(p^{2})$ in Eq. (27) (right panel). Propagators of
ground-state $\rho$ meson and its first radial excitation $\rho(1450)$ given
by formula (26) versus Euclidean $p^{2}$ are drawn with solid and dashed
lines. The propagators approach
$\left[\Lambda^{2}h_{\mathcal{Q}}^{2}/g^{2}C_{\mathcal{Q}}^{2}\right]^{-1}$ at
large Euclidean $p^{2}$ and coincide with naive propagators
$\left(p^{2}+M_{\mathcal{Q}}^{2}\right)^{-1}$ in the vicinity of
$p^{2}=-M_{\mathcal{Q}}^{2}$. Corresponding contributions to
$\tilde{\Pi}^{(2)}(p^{2})$ decrease with radial number $n$.

It has to be stressed that all three elements—meson polarization operator,
quark-meson vertices and $V\to\gamma$ amplitudes—are completely determined by
the given mean field, and their overall relevance to various meson properties
have already been tested in computation of the masses, decay constant and form
factors.

Figure 7: The diagram for the contribution of two intermediate mesons to HVP.

Contribution coming from the lowest order in number of intermediate mesons HVP
diagrams (Figs. 2 and 3) clearly underestimates the total HVP part of
$a_{\mu}$. However, there are many contributions which have not been taken
into account yet. First of all, these are diagrams shown in Fig. 7 with two
intermediate pseudoscalar mesons including radially excited ones. There are no
ad hoc arguments that their contribution must be particularly small. As a
matter of fact, this question addresses the role of higher in
$1/N_{\mathrm{c}}$ contributions. One could also consider diagrams that would
generate imaginary part of meson propagators in Minkowski kinematics. These
diagrams are of higher order in $1/N_{\mathrm{c}}$ than the diagrams discussed
in the article. Diagrams with multiple internal meson lines could be also
examined. Calculation of these diagrams is quite complicated, the result will
be reported in due course.

The simplest quark loop contribution to HLbL part is quite close to both data-
driven and lattice values, but only if the impact of Wilson loop is properly
taken into account, see Table 5. It is clear that diagrams with intermediate
mesons, including the excited ones, have to be computed as well before coming
to some definite conclusions.

Figure 8: The quark two-loop HLbL diagram due to the mean-field correlator,
quark loops are averaged jointly, in accordance with Eq.(11).

Besides these meson-exchange diagrams, there exists a specific contribution
illustrated by the diagram in Fig. 8 which corresponds to the exchange between
virtual quarks inside loop by (infinitely) many soft gluons, represented in
the domain model by the mean-field correlator. Obviously this diagram is a
source of uncertainty unless the mean-field correlators are determined.

Systematic calculations of all these and some other contributions within the
domain model are straightforward but rather involved, and their results will
be reported over time.

##  Appendix A Contributions to HVP

In order to find HVP from formula (4) where quark fields are integrated out,
it is convenient to introduce _external_ field $A_{\mu}$ and use relation

$\tilde{\Pi}_{\mu\nu}(p)=\int d^{4}xe^{-ipx}\langle
j_{\mu}(x)j_{\nu}(0)\rangle=\int
d^{4}xe^{-ipx}\left.\frac{1}{e^{2}}\frac{\delta}{\delta
A_{\mu}(x)}\frac{\delta}{\delta A_{\nu}(0)}Z[A]\right|_{A=0},$ (33)

where $j_{\mu}=\bar{\psi}Q\gamma_{\mu}\psi$ are quark currents. Functional
$Z[A]$ is given by (cf. Eq. (19))

$Z[A]={\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^{4}x\phi^{2}_{\cal
Q}(x)-\Gamma(A)-\sum\limits_{k=2}^{\infty}\frac{1}{k}W_{k}[\phi|A]\right\\},\quad
Z[0]=1.$

Now, terms that are at most quadratic in fields $\phi_{\mathcal{Q}},A_{\mu}$
are retained in the exponential, while higher-order interactions are expanded
perturbatively:

$\displaystyle Z[A]$ $\displaystyle={\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^{4}x\phi^{2}_{\cal
Q}(x)+\frac{1}{2}e^{2}\int d^{4}x\ d^{4}y\
A_{\mu}(x)\Pi_{\mu\nu}^{(1)}(x-y)A_{\nu}(y)\right.$
$\displaystyle\quad\left.-\frac{1}{2}\int d^{4}x\ d^{4}y\ \phi_{\cal
Q}(x)h^{2}_{\cal Q}\Gamma_{\mathcal{Q}}^{(2)}(x-y)\phi_{\cal Q}(y)-e\int
d^{4}x\ d^{4}y\ \phi_{\cal Q}(x)h_{\cal
Q}\Gamma_{\mathcal{Q};\mu}(x-y)A_{\mu}(y)\right\\}$
$\displaystyle\quad\times\left(1+\text{higher-order interactions}\right).$

When higher-order terms are neglected, one performs Gaussian integration and
finds

$\displaystyle Z[A]$ $\displaystyle=\exp\left\\{\frac{1}{2}e^{2}\int d^{4}p\
\tilde{A}_{\mu}(-p)\left(\tilde{\Pi}_{\mu\nu}^{(1)}(p)+\tilde{\Gamma}_{\mathcal{Q};\mu}(p)\left[\frac{\Lambda^{2}h_{\mathcal{Q}}^{2}}{g^{2}C^{2}_{\mathcal{Q}}}+h_{\mathcal{Q}}^{2}\tilde{\Gamma}^{(2)}_{\cal
Q}(p)\right]^{-1}\tilde{\Gamma}_{\mathcal{Q};\nu}(p)\right)\tilde{A}_{\nu}(p)\right\\}$
$\displaystyle=\exp\left\\{\frac{1}{2}e^{2}\int d^{4}p\
\tilde{A}_{\mu}(-p)\left(\tilde{\Pi}_{\mu\nu}^{(1)}(p)+\tilde{\Pi}_{\mu\nu}^{(2)}(p)\right)\tilde{A}_{\nu}(p)\right\\}$
$\displaystyle=\exp\left\\{\frac{1}{2}e^{2}\int d^{4}x\ d^{4}y\
A_{\mu}(x)\left(\Pi_{\mu\nu}^{(1)}(x-y)+\Pi_{\mu\nu}^{(2)}(x-y)\right)A_{\nu}(y)\right\\}.$

Substituting this formula into Eq. (33), one recovers diagrams shown in Figs.
2 and 3.

##  Appendix B Quark loop contribution to HVP

Upon evaluating trace of Dirac matrices and performing loop momentum
integration in formula (24), one finds

$\displaystyle(-1)\ {\rm reg}\int\frac{d^{4}l}{(2\pi)^{4}}\mathrm{Tr}\
\gamma_{\mu}\tilde{H}_{f}(l+p/2)\gamma_{\nu}\tilde{H}_{f}(l-p/2)=(p^{2}\delta_{\mu\nu}-p_{\mu}p_{\nu})\tilde{\Pi}_{f1}^{(1)}(p^{2})+f_{\mu\alpha}p_{\alpha}f_{\nu\beta}p_{\beta}\tilde{\Pi}_{f2}^{(1)}(p^{2})+\delta_{\mu\nu}\Delta.$

The first two structures with

$\displaystyle\tilde{\Pi}_{f1}^{(1)}(p^{2})$
$\displaystyle=-\mathrm{Tr}\frac{1}{(4\pi)^{2}}\mathrm{reg}\int_{0}^{1}ds_{1}\int_{0}^{1}ds_{2}\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]^{\frac{m_{f}^{2}}{4v\Lambda^{2}}}\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}\right)\frac{8s_{1}s_{2}}{(s_{1}+s_{2})^{4}},$
$\displaystyle\tilde{\Pi}_{f2}^{(1)}(p^{2})$
$\displaystyle=-\mathrm{Tr}\frac{1}{(4\pi)^{2}}\mathrm{reg}\int_{0}^{1}ds_{1}\int_{0}^{1}ds_{2}\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]^{\frac{m_{f}^{2}}{4v\Lambda^{2}}}\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}\right)\frac{8s_{1}^{2}s_{2}^{2}}{(s_{1}+s_{2})^{4}},$

where the trace is over color indices, satisfy Ward identity (tensor
$f_{\mu\nu}$ is antisymmetric), and the third term is equal to zero. In order
to show this, let us represent $\Delta$ as a sum of two terms

$\displaystyle\Delta(p^{2})$
$\displaystyle=-\Delta_{1}(p^{2})-\Delta_{2}(p^{2}),$
$\displaystyle\Delta_{1}(p^{2})$
$\displaystyle=\mathrm{Tr}\frac{2v\Lambda^{2}}{(4\pi)^{2}}\mathrm{reg}\int\limits_{0}^{1}\int\limits_{0}^{1}\frac{ds_{1}ds_{2}}{(s_{1}+s_{2})^{3}}\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]^{\frac{m_{f}^{2}}{4v\Lambda^{2}}}\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}\right)$
$\displaystyle\quad\times\left[-4\frac{p^{2}}{2v\Lambda^{2}}\frac{s_{1}s_{2}}{(s_{1}+s_{2})}-4-4\frac{m_{f}^{2}}{2v\Lambda^{2}}\frac{(1-s_{1}s_{2})(s_{1}+s_{2})}{(1-s_{1}^{2})(1-s_{2}^{2})}\right],$
$\displaystyle\Delta_{2}(p^{2})$
$\displaystyle=\mathrm{Tr}\frac{2v\Lambda^{2}}{(4\pi)^{2}}\mathrm{reg}\int\limits_{0}^{1}\int\limits_{0}^{1}\frac{ds_{1}ds_{2}}{(s_{1}+s_{2})^{3}}\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]^{\frac{m_{f}^{2}}{4v\Lambda^{2}}}\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}\right)$
$\displaystyle\quad\times
s_{1}s_{2}\left[-4\frac{p^{2}}{2v\Lambda^{2}}\frac{s_{1}s_{2}}{(s_{1}+s_{2})}+4-4\frac{m_{f}^{2}}{2v\Lambda^{2}}\frac{(1-s_{1}s_{2})(s_{1}+s_{2})}{(1-s_{1}^{2})(1-s_{2}^{2})}\right].$

Then $\Delta_{1}$ can be identically rewritten as

$\displaystyle\Delta_{1}(p^{2})$ $\displaystyle=4\lim_{\rho\to
1}\mathrm{Tr}\frac{2v\Lambda^{2}}{(4\pi)^{2}}\mathrm{reg}\int\limits_{0}^{1}\int\limits_{0}^{1}\frac{ds_{1}ds_{2}}{(s_{1}+s_{2})^{3}}\rho^{-1}\frac{d}{d\rho}\rho^{-1}\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\rho\frac{p^{2}}{2v\Lambda^{2}}\right)$
$\displaystyle\quad\times\exp\left(\frac{m_{f}^{2}}{4v\Lambda^{2}}\ln\left[\frac{(1-\rho
s_{1})(1-\rho s_{2})}{(1+\rho s_{1})(1+\rho s_{2})}\right]\right)$
$\displaystyle=4\lim_{\rho\to
1}\mathrm{Tr}\frac{2v\Lambda^{2}}{(4\pi)^{2}}\frac{d}{d\rho}\mathrm{reg}\int\limits_{0}^{\rho}\int\limits_{0}^{\rho}\frac{ds_{1}ds_{2}}{(s_{1}+s_{2})^{3}}\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}+\frac{m_{f}^{2}}{4v\Lambda^{2}}\ln\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]\right)=0.$

The second term $\Delta_{2}$ is rewritten in a similar way

$\displaystyle\Delta_{2}(p^{2})$ $\displaystyle=4\lim_{\rho\to
1}\mathrm{Tr}\frac{2v\Lambda^{2}}{(4\pi)^{2}}\mathrm{reg}\int\limits_{0}^{1}\int\limits_{0}^{1}\frac{ds_{1}ds_{2}s_{1}s_{2}}{(s_{1}+s_{2})^{3}}\rho\frac{d}{d\rho}\rho\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\rho\frac{p^{2}}{2v\Lambda^{2}}\right)$
$\displaystyle\quad\times\exp\left(\frac{m_{f}^{2}}{4v\Lambda^{2}}\ln\left[\frac{(1-\rho
s_{1})(1-\rho s_{2})}{(1+\rho s_{1})(1+\rho s_{2})}\right]\right)$
$\displaystyle=4\lim_{\rho\to
1}\mathrm{Tr}\frac{2v\Lambda^{2}}{(4\pi)^{2}}\frac{d}{d\rho}\mathrm{reg}\int\limits_{0}^{\rho}\int\limits_{0}^{\rho}\frac{ds_{1}ds_{2}s_{1}s_{2}}{(s_{1}+s_{2})^{3}}\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}+\frac{m_{f}^{2}}{4v\Lambda^{2}}\ln\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]\right)=0.$

Averaging $\tilde{\Pi}_{\mu\nu}(p^{2})$ over background gluon field and
renormalizing at $p^{2}=0$, one arrives at

$\displaystyle\tilde{\Pi}^{(1)\mathrm{R}}_{\mu\nu}(p)$
$\displaystyle=(p^{2}\delta_{\mu\nu}-p_{\mu}p_{\nu})\sum_{f}Q_{f}^{2}\tilde{\Pi}_{f}^{(1)\mathrm{R}}(p^{2}),$
$\displaystyle\tilde{\Pi}_{f}^{(1)\mathrm{R}}(p^{2})$
$\displaystyle=\tilde{\Pi}^{(1)\mathrm{R}}_{f1}(p^{2})+\frac{1}{3}\tilde{\Pi}^{(1)\mathrm{R}}_{f2}(p^{2}),$
$\displaystyle\tilde{\Pi}^{(1)\mathrm{R}}_{f1}(p^{2})$
$\displaystyle=-\mathrm{Tr}\frac{1}{2\pi^{2}}\int_{0}^{1}ds_{1}\int_{0}^{1}ds_{2}\frac{s_{1}s_{2}}{(s_{1}+s_{2})^{4}}\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]^{\frac{m_{f}^{2}}{4v\Lambda^{2}}}\left[\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}\right)-1\right],$
$\displaystyle\tilde{\Pi}^{(1)\mathrm{R}}_{f2}(p^{2})$
$\displaystyle=-\mathrm{Tr}\frac{1}{2\pi^{2}}\int_{0}^{1}ds_{1}\int_{0}^{1}ds_{2}\frac{s_{1}^{2}s_{2}^{2}}{(s_{1}+s_{2})^{4}}\left[\frac{(1-s_{1})(1-s_{2})}{(1+s_{1})(1+s_{2})}\right]^{\frac{m_{f}^{2}}{4v\Lambda^{2}}}\left[\exp\left(-\frac{s_{1}s_{2}}{s_{1}+s_{2}}\frac{p^{2}}{2v\Lambda^{2}}\right)-1\right].$

## References

  * Albahri  _et al._ [2021] T. Albahri _et al._ (Muon g-2), Measurement of the anomalous precession frequency of the muon in the Fermilab Muon $g-2$ Experiment, Phys. Rev. D 103, 072002 (2021), arXiv:2104.03247 [hep-ex] .
  * Aoyama  _et al._ [2020] T. Aoyama _et al._ , The anomalous magnetic moment of the muon in the Standard Model, Phys. Rept. 887, 1 (2020), arXiv:2006.04822 [hep-ph] .
  * Jegerlehner and Nyffeler [2009] F. Jegerlehner and A. Nyffeler, The Muon g-2, Phys. Rept. 477, 1 (2009), arXiv:0902.3360 [hep-ph] .
  * Efimov and Nedelko [1995] G. V. Efimov and S. N. Nedelko, Nambu-Jona-Lasinio model with the homogeneous background gluon field, Phys. Rev. D51, 176 (1995).
  * Burdanov  _et al._ [1996] Ja. V. Burdanov, G. V. Efimov, S. N. Nedelko, and S. A. Solunin, Meson masses within the model of induced nonlocal quark currents, Phys. Rev. D54, 4483 (1996), arXiv:hep-ph/9601344 [hep-ph] .
  * Kalloniatis and Nedelko [2004] A. C. Kalloniatis and S. N. Nedelko, Realization of chiral symmetry in the domain model of QCD, Phys. Rev. D69, 074029 (2004), [Erratum: Phys. Rev.D70,119903(2004)], arXiv:hep-ph/0311357 [hep-ph] .
  * Nedelko and Voronin [2016] S. N. Nedelko and V. E. Voronin, Regge spectra of excited mesons, harmonic confinement and QCD vacuum structure, Phys. Rev. D93, 094010 (2016), arXiv:1603.01447 [hep-ph] .
  * Nedelko and Voronin [2017] S. N. Nedelko and V. E. Voronin, Influence of confining gluon configurations on the $P\to\gamma^{*}\gamma$ transition form factors, Phys. Rev. D95, 074038 (2017), arXiv:1612.02621 [hep-ph] .
  * Nedelko and Voronin [2015] S. N. Nedelko and V. E. Voronin, Domain wall network as QCD vacuum and the chromomagnetic trap formation under extreme conditions, Eur. Phys. J. A51, 45 (2015), arXiv:1403.0415 [hep-ph] .
  * Nedelko and Voronin [2021] S. N. Nedelko and V. E. Voronin, Energy-driven disorder in mean field QCD, Phys. Rev. D 103, 114021 (2021), arXiv:2012.07081 [hep-ph] .
  * Lepage and Brodsky [1980] G. P. Lepage and S. J. Brodsky, Exclusive Processes in Perturbative Quantum Chromodynamics, Phys. Rev. D 22, 2157 (1980).
  * Burdanov  _et al._ [1998] Y. V. Burdanov, G. V. Efimov, and S. N. Nedelko, Selfdual homogeneous gluon field and electromagnetic structure of pion, arXiv:hep-ph/9806478 (1998).
  * Lautrup  _et al._ [1972] B. E. Lautrup, A. Peterman, and E. de Rafael, Recent developments in the comparison between theory and experiments in quantum electrodynamics, Phys. Rept. 3, 193 (1972).
  * Blum [2003] T. Blum, Lattice calculation of the lowest order hadronic contribution to the muon anomalous magnetic moment, Phys. Rev. Lett. 91, 052001 (2003), arXiv:hep-lat/0212018 .
  * Dorokhov [2004] A. E. Dorokhov, Adler function and hadronic contribution to the muon g-2 in a nonlocal chiral quark model, Phys. Rev. D 70, 094011 (2004), arXiv:hep-ph/0405153 .
  * de Rafael [1994] E. de Rafael, Hadronic contributions to the muon g-2 and low-energy QCD, Phys. Lett. B 322, 239 (1994), arXiv:hep-ph/9311316 .
  * Peris  _et al._ [1998] S. Peris, M. Perrottet, and E. de Rafael, Matching long and short distances in large N(c) QCD, JHEP 05, 011, arXiv:hep-ph/9805442 .
  * Knecht [2004] M. Knecht, The Anomalous magnetic moment of the muon: A Theoretical introduction, Lect. Notes Phys. 629, 37 (2004), arXiv:hep-ph/0307239 .
  * Holdom  _et al._ [1994] B. Holdom, R. Lewis, and R. R. Mendel, Hadronic contribution to the photon vacuum polarization: A Theoretical estimate, Z. Phys. C 63, 71 (1994), arXiv:hep-ph/9304264 .
  * Davier  _et al._ [2017] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Reevaluation of the hadronic vacuum polarisation contributions to the Standard Model predictions of the muon $g-2$ and ${\alpha(m_{Z}^{2})}$ using newest hadronic cross-section data, Eur. Phys. J. C 77, 827 (2017), arXiv:1706.09436 [hep-ph] .
  * Keshavarzi  _et al._ [2018] A. Keshavarzi, D. Nomura, and T. Teubner, Muon $g-2$ and $\alpha(M_{Z}^{2})$: a new data-based analysis, Phys. Rev. D 97, 114025 (2018), arXiv:1802.02995 [hep-ph] .
  * Colangelo  _et al._ [2019] G. Colangelo, M. Hoferichter, and P. Stoffer, Two-pion contribution to hadronic vacuum polarization, JHEP 02, 006, arXiv:1810.00007 [hep-ph] .
  * Hoferichter  _et al._ [2019] M. Hoferichter, B.-L. Hoid, and B. Kubis, Three-pion contribution to hadronic vacuum polarization, JHEP 08, 137, arXiv:1907.01556 [hep-ph] .
  * Davier  _et al._ [2020] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, A new evaluation of the hadronic vacuum polarisation contributions to the muon anomalous magnetic moment and to $\mathbf{\boldsymbol{\alpha}(m_{Z}^{2})}$, Eur. Phys. J. C 80, 241 (2020), [Erratum: Eur.Phys.J.C 80, 410 (2020)], arXiv:1908.00921 [hep-ph] .
  * Keshavarzi  _et al._ [2020] A. Keshavarzi, D. Nomura, and T. Teubner, $g-2$ of charged leptons, $\alpha(M^{2}_{Z})$ , and the hyperfine splitting of muonium, Phys. Rev. D 101, 014029 (2020), arXiv:1911.00367 [hep-ph] .
  * Chakraborty  _et al._ [2018] B. Chakraborty _et al._ (Fermilab Lattice, LATTICE-HPQCD, MILC), Strong-Isospin-Breaking Correction to the Muon Anomalous Magnetic Moment from Lattice QCD at the Physical Point, Phys. Rev. Lett. 120, 152001 (2018), arXiv:1710.11212 [hep-lat] .
  * Borsanyi  _et al._ [2018] S. Borsanyi _et al._ (Budapest-Marseille-Wuppertal), Hadronic vacuum polarization contribution to the anomalous magnetic moments of leptons from first principles, Phys. Rev. Lett. 121, 022002 (2018), arXiv:1711.04980 [hep-lat] .
  * Blum  _et al._ [2018] T. Blum, P. A. Boyle, V. Gülpers, T. Izubuchi, L. Jin, C. Jung, A. Jüttner, C. Lehner, A. Portelli, and J. T. Tsang (RBC, UKQCD), Calculation of the hadronic vacuum polarization contribution to the muon anomalous magnetic moment, Phys. Rev. Lett. 121, 022003 (2018), arXiv:1801.07224 [hep-lat] .
  * Giusti  _et al._ [2019] D. Giusti, V. Lubicz, G. Martinelli, F. Sanfilippo, and S. Simula, Electromagnetic and strong isospin-breaking corrections to the muon $g-2$ from Lattice QCD+QED, Phys. Rev. D 99, 114502 (2019), arXiv:1901.10462 [hep-lat] .
  * Shintani and Kuramashi [2019] E. Shintani and Y. Kuramashi (PACS), Hadronic vacuum polarization contribution to the muon $g-2$ with 2+1 flavor lattice QCD on a larger than (10 fm$)^{4}$ lattice at the physical point, Phys. Rev. D 100, 034517 (2019), arXiv:1902.00885 [hep-lat] .
  * Davies  _et al._ [2020] C. T. H. Davies _et al._ (Fermilab Lattice, LATTICE-HPQCD, MILC), Hadronic-vacuum-polarization contribution to the muon’s anomalous magnetic moment from four-flavor lattice QCD, Phys. Rev. D 101, 034512 (2020), arXiv:1902.04223 [hep-lat] .
  * Gérardin  _et al._ [2019] A. Gérardin, M. Cè, G. von Hippel, B. Hörz, H. B. Meyer, D. Mohler, K. Ottnad, J. Wilhelm, and H. Wittig, The leading hadronic contribution to $(g-2)_{\mu}$ from lattice QCD with $N_{\rm f}=2+1$ flavours of O($a$) improved Wilson quarks, Phys. Rev. D 100, 014510 (2019), arXiv:1904.03120 [hep-lat] .
  * Aubin  _et al._ [2020] C. Aubin, T. Blum, C. Tu, M. Golterman, C. Jung, and S. Peris, Light quark vacuum polarization at the physical point and contribution to the muon $g-2$, Phys. Rev. D 101, 014503 (2020), arXiv:1905.09307 [hep-lat] .
  * Giusti and Simula [2019] D. Giusti and S. Simula, Lepton anomalous magnetic moments in Lattice QCD+QED, PoS LATTICE2019, 104 (2019), arXiv:1910.03874 [hep-lat] .
  * Borsanyi  _et al._ [2021] S. Borsanyi _et al._ , Leading hadronic contribution to the muon magnetic moment from lattice QCD, Nature 593, 51 (2021), arXiv:2002.12347 [hep-lat] .
  * Aldins  _et al._ [1970] J. Aldins, T. Kinoshita, S. J. Brodsky, and A. J. Dufner, Photon-photon scattering contribution to the sixth order magnetic moments of the muon and electron, Phys. Rev. D 1, 2378 (1970).
  * Knecht and Nyffeler [2002] M. Knecht and A. Nyffeler, Hadronic light by light corrections to the muon g-2: The Pion pole contribution, Phys. Rev. D 65, 073034 (2002), arXiv:hep-ph/0111058 .
  * Barbieri and Remiddi [1975] R. Barbieri and E. Remiddi, Electron and Muon 1/2(g-2) from Vacuum Polarization Insertions, Nucl. Phys. B 90, 233 (1975).
  * Dorokhov  _et al._ [2015] A. E. Dorokhov, A. E. Radzhabov, and A. S. Zhevlakov, Dynamical quark loop light-by-light contribution to muon g-2 within the nonlocal chiral quark model, Eur. Phys. J. C 75, 417 (2015), arXiv:1502.04487 [hep-ph] .
  * Goecke  _et al._ [2013] T. Goecke, C. S. Fischer, and R. Williams, Role of momentum dependent dressing functions and vector meson dominance in hadronic light-by-light contributions to the muon $g-2$, Phys. Rev. D 87, 034013 (2013), arXiv:1210.1759 [hep-ph] .
  * Greynat and de Rafael [2012] D. Greynat and E. de Rafael, Hadronic Contributions to the Muon Anomaly in the Constituent Chiral Quark Model, JHEP 07, 020, arXiv:1204.3029 [hep-ph] .
  * Bijnens  _et al._ [1996] J. Bijnens, E. Pallante, and J. Prades, Analysis of the hadronic light by light contributions to the muon g-2, Nucl. Phys. B 474, 379 (1996), arXiv:hep-ph/9511388 .
  * Hayakawa  _et al._ [1995] M. Hayakawa, T. Kinoshita, and A. I. Sanda, Hadronic light by light scattering effect on muon g-2, Phys. Rev. Lett. 75, 790 (1995), arXiv:hep-ph/9503463 .
  * Blum  _et al._ [2020] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Hadronic Light-by-Light Scattering Contribution to the Muon Anomalous Magnetic Moment from Lattice QCD, Phys. Rev. Lett. 124, 132002 (2020), arXiv:1911.08123 [hep-lat] .
  * Kuhn  _et al._ [2003] J. H. Kuhn, A. I. Onishchenko, A. A. Pivovarov, and O. L. Veretin, Heavy mass expansion, light by light scattering and the anomalous magnetic moment of the muon, Phys. Rev. D 68, 033018 (2003), arXiv:hep-ph/0301151 .