---
author:
- |
    \
    Theoretical Physics Group, RIKEN, Wako 2-1, Saitama 351-0198, Japan\
    E-mail:
- |
    Thomas Blum\
    Physics Department, University of Connecticut, Storrs, CT 06269-3046, USA\
    RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, New York 11973, USA\
    E-mail:
- |
    Taku Izubuchi\
    Kanazawa University, Institute for Theoretical Physics, Kanazawa 920-1192, Japan\
    RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, New York 11973, USA\
    E-mail:
- |
    Norikazu Yamada\
    High Energy Accelerator Research Organization [(]{}KEK[)]{}, Tsukuba, Ibaraki 305-0801, Japan\
    The Graduate University for Advanced Studies, Tsukuba, Ibaraki 305-0801, Japan\
    E-mail:
title: 'Hadronic light-by-light scattering contribution to the muon $g-2$ from lattice QCD:Methodology'
---

The muon anomalous magnetic moment ($g-2$) has been measured very accurately at BNL [@Bennett:2002jb] $$\Delta a_\mu({\rm EXP}) = 6 \times 10^{-10}.  
  \label{eq:BNL_accuracy}$$ Theoretical interest in the muon $g-2$ stems from the fact that it can be used as a probe of unknown microscopic structures, such as supersymmetry with its breaking scale $m_S \sim {\cal O}(100)$ GeV and large $\tan \beta (\ge 10)$, and TeV scale gravity [@Arkani-Hamed:1998rs], each of which can give the additional contribution $a_\mu({\rm new\ physics}) = {\cal O}(10 \sim 100) \times 10^{-10}$. We recall that both of these two structures are the most promising possibilities to solve the hierarchy problem, i.e., the hierarchy between the electroweak scale and the GUT or Planck scale.  Thus, the currently available accuracy of $\Delta a_\mu({\rm EXP})$ can provide us with an opportunity to examine if the underlying structures ensuring such a hierarchy exist. 

Table \[table:theory\_exp\] summarizes the current status of the muon $g-2$ [@Hertzog:2005mx], showing $2.7 \sigma$ deviation of $a_\mu({\rm EXP})$ from the standard model prediction $a_\mu({\rm SM})$ [@Bennett:2002jb; @Hertzog:2005mx].  However, the uncertainty in $a_\mu({\rm SM})$, which is attributed to the QCD contribution, must be reduced to elaborate such a discrepancy. 

                                           $a_\mu \times 10^{10}$ $\Delta a_\mu \times 10^{10}$
  -------------------------------------- ------------------------ -------------------------------
  QED                                           $11\ 658\ 471.94$ $0.14$
  [QCD]{}                                          ${\red 695.4}$ ${\red 7.3}$
  Weak                                                     $15.4$ $0.22$
  Theory                                         $11\ 659\ 182.7$ ${\red 7.4}$
  Experiment                                   $11\ 659\ 208\ \ $ $6$
  $a_\mu({\rm EXP}) - a_\mu({\rm SM})$                   [25.3]{} [9.5]{}

  : Comparison of $a_\mu({\rm SM})$ with $a_\mu({\rm EXP})$[]{data-label="table:theory_exp"}

Table \[table:muong-2fromQCD\] summarizes the QCD contribution to $a_\mu({\rm SM})$ as well as its various components which are relevant in considering $\Delta a_\mu({\rm EXP})$ in Eq.  (\[eq:BNL\_accuracy\]).  The QCD contribution starts at ${\cal O}(\alpha_{\rm em}^2)$ through the hadronic contribution to the vacuum polarization of the photon, which has been computed with lattice QCD [@Blum:2002ii], and by using chiral perturbation theory and a vector meson model to fit the lattice data [@Aubin].  The improvement of the accuracy of this contribution either relies on the further precise measurement of $\sigma(e^+ e^- \rightarrow {\rm hadrons})(\sqrt{s})$ [@Hagiwara:2003da], or on the improvement of the lattice calculation with a statistical uncertainty less than 5 % [@Aubin] (but with less well-known systematic uncertainty). 

                                                                             $a_\mu \times 10^{10}$ $\Delta a_\mu \times 10^{10}$
  ------------------------------------------------------------------------ ------------------------ -------------------------------
  hadronic vacuum polarization (${\cal O}(\alpha_{\rm em}^{{\blue 2}})$)                    $693.4$ $6.4$
  hadronic vacuum polarization (${\cal O}(\alpha_{\rm em}^{{\blue 3}})$)                    $-10.0$ $0.6$
  [hadronic light-by-light]{} (${\cal O}(\alpha_{\rm em}^{{\blue 3}})$)               ${\red 12.0}$ ${\red 3.5}$
  Total QCD                                                                                 $695.4$ $7.3$

  : QCD contribution to the muon $g-2$ [@Hertzog:2005mx][]{data-label="table:muong-2fromQCD"}

Our target here is another type of QCD contribution, namely [*the hadronic light-by-light scattering [(]{}h-lbl[)]{} contribution*]{} to the muon $g-2$ shown in Fig.  \[fig:hadronic\_light-by-light\], which also affects the physical interpretation of the muon $g-2$ in viewing Table \[table:muong-2fromQCD\].  The diagram in Fig.  \[fig:hadronic\_light-by-light\] arises through the elastic scattering amplitude of two (off-shell) photons [*by QCD*]{} (the blob).  To date, this contribution could be estimated by purely theoretical calculation.  So far, it has been calculated only based on the hadronic picture [@hadronic_lbl; @pseudoscalar].  Thus the first principle calculation based on lattice QCD is particularly desirable. 

(200,160)(-120,0) (-50,0)(-100,0) (0,0)(-50,0) (50,0)(0,0) (100,0)(50,0) (80,-5)\[t\][$\mu$]{} (0,80)(50,50)(45)[Gray]{}[Gray]{} (-50,0)(-40,50)[3]{}[4]{} (0,0)(0,30)[3]{}[3]{} (50,0)(40,50)[3]{}[4]{} (0,130)(0,160)[3]{}[2.5]{} (-50,0)[3]{} (-0,0)[3]{} (50,0)[3]{} (0,130)[3]{} (-40,50)[3]{} (0,30)[3]{} (40,50)[3]{} (150,100)\[t\][elastic scattering amplitude]{} (150,85)\[t\][of two photons by QCD]{} (22,15)(13,-40,220) (-25,15)(13,-40,220) (22,20)\[t\][$l_1$]{} (-25,20)\[t\][$l_2$]{}

The diagram in Fig.  \[fig:hadronic\_light-by-light\] evokes the following naive approach; we calculate [*repeatedly*]{} the correlation function of [*four*]{} hadronic electromagnetic currents by lattice QCD [*with respect to two independent four-momenta*]{} $l_1,\,l_2$ of off-shell photons, and integrate it over $l_1,\,l_2$. Such a task is too difficult to accomplish with use of supercomputers available in the foreseeable future.  Here we propose a practical method to calculate the h-lbl contribution by using the lattice (QCD + [*QED*]{}) simulation; we compute $$\begin{aligned}
 && 
 \left< 
  \begin{picture}(70,40)(-5,20) 
       \SetWidth{1.2} 
    \ArrowArc(30,40)(15,90,180) 
    \ArrowArc(30,40)(15,180,270) 
    \ArrowArc(30,40)(15,270,360) 
    \ArrowArc(30,40)(15,0,90) 
    \Text(30,43)[t]{\scriptsize quark} 
       \SetWidth{1.2} 
    \SetColor{Red} 
    \Photon(30,55)(30,65){1.5}{2.5} 
    \Photon(30,25)(30,5){1.5}{4} 
       \SetColor{Green} 
    \Vertex(30,55){2} 
    \Vertex(30,25){2} 
       \SetWidth{1.2} 
    \SetColor{Blue} 
    \ArrowLine(60,5)(30,5) 
    \ArrowLine(30,5)(0,5) 
       \SetColor{Green} 
    \Vertex(30,5){2}
  \end{picture} 
 \right>_{{\rm QCD} + {\rm quenched\ QED}_A} 
  - 
 \left< 
  \begin{picture}(70,20)(20,0) 
      \SetWidth{1.2} 
   \ArrowArc(55,0)(15,90,180) 
   \ArrowArc(55,0)(15,180,270) 
   \ArrowArc(55,0)(15,270,360) 
   \ArrowArc(55,0)(15,0,90) 
   \Text(55,3)[t]{\scriptsize quark} 
      \SetWidth{1.2} 
   \SetColor{Red} 
   \Photon(55,15)(55,25){1.5}{2.5} 
   \Photon(55,-15)(55,-35){1.5}{4} 
      \SetColor{Green} 
   \Vertex(55,15){2} 
   \Vertex(55,-15){2} 
  \end{picture} 
 \right>_{{\rm QCD} + {\rm quenched\ QED}_B} 
  \nonumber \\ 
 && 
 \qquad 
 \qquad \qquad \qquad \qquad \qquad \qquad 
 \qquad \,  
  \left< 
   \begin{picture}(70,20)(-5,0) 
        \SetWidth{1.2} 
    \SetColor{Blue} 
    \ArrowLine(60,5)(30,5) 
    \ArrowLine(30,5)(0,5) 
        \SetColor{Green} 
    \Vertex(30,5){2} 
   \end{picture} 
  \right>_{{\rm quenched\ QED}_A} \, , 
   \label{eq:method_figure} \end{aligned}$$ amputate the external muon lines, and project the magnetic form factor, and divide by the factor $3$. In Eq.  (\[eq:method\_figure\]) the red line denotes the free photon propagator $D_{\mu\nu}(x,\,y)$ in the non-compact lattice QED solved in an appropriate gauge fixing condition.  The black line denotes the full quark propagator $S_f(x,\,y;U,\,u)$ for a given set of $SU(3)_C$ gauge configuration $\left\{U_{x,\,\mu}\right\}$ and $U(1)_{\rm em}$ gauge configuration $\left\{u_{x,\,\mu}\right\}$, where the sum over relevant flavors $f$ is implicitly assumed.  The blue line represents the full muon propagator $s(x,\,y;\,u)$. The average $\left<,\,\right>$ above means the one over the unquenched $SU(3)_C$ gauge configurations and/or the quenched $U(1)_{\rm em}$ gauge configurations [^1] as specified by the subscript attached to it.  Since two statistically independent averages over $U(1)_{\rm em}$ gauge configurations appear in the second term, they are distinguished by the labels $A$, $B$.

To explain the mechanism underlying our method, let us look at the first term of Eq.  (\[eq:method\_figure\]) [*perturbatively with respect to QED*]{}.  Its magnetic components up to ${\cal O}(\alpha_{\rm em}^3)$ consist of $$\begin{aligned}
  && {\cal O}(\alpha_{\rm em}^{2}) 
    + {\cal O}(\alpha_{\rm em}^{3}) \qquad \qquad 
 \left< 
  \begin{picture}(80,20)(20,0) 
      \SetWidth{1.2} 
   \ArrowArc(60,0)(15,90,180) 
   \ArrowArc(60,0)(15,180,270) 
   \ArrowArc(60,0)(15,270,360) 
   \ArrowArc(60,0)(15,0,90) 
   \Text(60,3)[t]{\scriptsize quark}
      \SetColor{Red} 
   \SetWidth{1.2} 
   \Photon(60,15)(60,25){1.5}{2.5} 
   \Photon(60,-15)(60,-35){1.5}{4} 
      \SetColor{Green} 
   \Vertex(60,15){2} 
   \Vertex(60,-15){2} 
  \end{picture} 
 \right>_{\rm QCD} 
  \nonumber \\ 
 && \qquad \qquad \qquad \qquad \qquad \quad \ \  
  \left< 
  \begin{picture}(80,20)(0,0) 
      \SetWidth{1.2} 
   \SetColor{Blue} 
   \ArrowLine(70,3)(50,3) 
   \Line(50,3)(40,3) 
   \Line(40,3)(30,3) 
   \ArrowLine(30,3)(10,3) 
      \SetWidth{1.2} 
   \SetColor{Red} 
   \PhotonArc(40,3)(10,180,360){1.5}{6} 
      \SetColor{Green} 
   \Vertex(50,3){2} 
   \Vertex(40,3){2} 
   \Vertex(30,3){2} 
  \end{picture} 
 \right>_{\rm 1 + {\cal O}(\alpha_{\rm em})}
 \nonumber \\ 
 \nonumber \\ 
 \nonumber \\ 
 && {\cal O}(\alpha_{\rm em}^{3}) \quad 
   \left< 
   \begin{picture}(80,20)(20,0) 
        \SetWidth{1.2} 
    \ArrowArc(60,0)(15,90,270) 
    \ArrowArc(60,0)(15,270,450) 
    \Text(60,3)[t]{\scriptsize quark} 
        \SetColor{Red} 
    \SetWidth{1.2} 
    \Photon(60,15)(60,25){1.5}{2.5} 
    \Photon(50,-10)(43,-35){1.5}{5} 
    \Photon(60,-15)(60,-35){1.5}{4} 
    \Photon(70,-10)(77,-35){1.5}{5} 
        \SetColor{Green} 
    \Vertex(60,15){2} 
    \Vertex(50,-10){2} 
    \Vertex(60,-15){2} 
    \Vertex(70,-10){2}
   \end{picture} 
  \right>_{\rm QCD} 
  \quad  
    \left< 
   \begin{picture}(80,20)(20,0) 
       \SetWidth{1.2} 
    \ArrowArc(60,0)(15,90,180) 
    \ArrowArc(60,0)(15,180,270) 
    \ArrowArc(60,0)(15,270,360) 
    \ArrowArc(60,0)(15,0,90) 
    \Text(60,3)[t]{\scriptsize quark}
       \SetColor{Red} 
    \SetWidth{1.2} 
    \Photon(60,15)(60,25){1.5}{2.5} 
    \Photon(60,-15)(60,-35){1.5}{4} 
       \SetColor{Green} 
    \Vertex(60,15){2} 
    \Vertex(60,-15){2} 
   \end{picture} 
  \right>_{{\rm QCD} \times {\cal O}(\alpha_{\rm em})}
 \nonumber \\ 
 && \qquad \qquad  
   \begin{picture}(80,20)(20,0) 
      \SetWidth{1.2} 
   \SetColor{Blue} 
   \ArrowLine(107,3)(89,3) 
   \Line(89,3)(72,3) 
   \Line(72,3)(55,3) 
   \ArrowLine(55,3)(37,3) 
      \SetColor{Green} 
   \Vertex(89,3){2} 
   \Vertex(72,3){2} 
   \Vertex(55,3){2} 
  \end{picture} 
  \qquad \qquad \ \   
    \begin{picture}(80,20)(20,0) 
    \SetWidth{1.2} 
  \SetColor{Blue} 
  \ArrowLine(100,3)(80,3) 
  \Line(80,3)(70,3) 
  \Line(70,3)(60,3) 
  \ArrowLine(60,3)(40,3) 
    \SetWidth{1.2} 
  \SetColor{Red} 
  \PhotonArc(70,3)(10,180,360){1.5}{6} 
    \SetColor{Green} 
  \Vertex(80,3){2} 
  \Vertex(70,3){2} 
  \Vertex(60,3){2} 
  \end{picture} 
   \, .  
   \label{eq:O(alpha)-expansion_first-term} 
$$ The diagram in the first line gives the ${\cal O}(\alpha_{\rm em}^2)$-contribution.  The ${\cal O}(\alpha_{\rm em})$-corrections to its muon part and to its quark part induce the ${\cal O}(\alpha_{\rm em}^3)$-contributions shown in the first line and in the right diagram on the second line respectively.  We recall that the QED gauge configurations in the first term of Eq.  (\[eq:method\_figure\]) are commonly shared by the quark part and the muon part.  Hence, the photons can be exchanged between the two parts.  As a consequence, the left diagram in the second line of Eq.  (\[eq:O(alpha)-expansion\_first-term\]) is induced at ${\cal O}(\alpha_{\rm em}^3)$, which takes the form of our target, Fig.  \[fig:hadronic\_light-by-light\].  Alternatively, the quark and muon parts in the first and third diagrams in Eq.  (\[eq:O(alpha)-expansion\_first-term\]) are connected only by a single photon attached a priori.  The second term in Eq.  (\[eq:method\_figure\]) also contains those extra diagrams.  Thus, by subtracting the second term from the first term, we may extract the h-lbl contribution.  The quantities evaluated in our method (\[eq:method\_figure\]) are constructed from [*two*]{} currents for both terms, which are surely less noisy than the case of [*four*]{} currents encountered in the naive approach.  Amazingly, the only difference between the first and second terms of Eq.  (\[eq:method\_figure\]) is the ways averaging over the $U(1)_{\rm em}$ gauge configurations.  The h-lbl contribution should thus emerge as such a subtle difference in averaging procedure. 

For our method to work efficiently, it is important to implement the cancellation of ${\cal O}(\alpha_{\rm em}^2)$-term in Eq.  (\[eq:O(alpha)-expansion\_first-term\]) without too much statistics.  Sharing the same ensemble of $U(1)_{\rm em}$ gauge configurations between the first term and the ${\rm QED}_A$ part of the second term is one such method since the ${\cal O}(\alpha_{\rm em}^2)$-correction to the muon vertex in both terms will be highly correlated.  Furthermore, we organize Eq.  (\[eq:method\_figure\]) in the manner; $$\begin{aligned}
 && 
 H_\mu(t_F,\,x_c,\,t_I) 
  \nonumber \\ 
 && \  
 \equiv 
 \frac{Z_\mu}{3} 
 \left< 
  \left< 
  \sum_{x,\,y} 
   \sum_{\vec{x}_F} 
    e^{-i \vec{p}_F \cdot \vec{x}_F}\, 
    s(x_F,\,x;\,u_A)\, 
      e \gamma_\rho    \sum_{\vec{x}_I} 
    e^{i \vec{p}_I \cdot \vec{x}_I}\, 
    s(x,\,x_I;\,u_A) 
  \right.  
 \right.  
  \nonumber \\ 
 && \qquad \ 
 \times 
 D_{\rho \lambda}(x,\,y) 
  \nonumber \\ 
 && \qquad  
 \times 
  \sum_f Q_f^2 
  \left( 
   {\rm tr} 
   \left( 
    \gamma_\mu     S_f(x_c,\,y;\,U,\,u_A)\, 
    e \gamma_\lambda     S_f(y,\,x_c;\,U,\,u_A) 
   \right) 
  \right.  
  \nonumber \\ 
 && \quad \qquad \qquad \ \ 
 \left.  
 \left.  
  \left.  
   - 
   \left< 
    {\rm tr} 
    \left( 
     \gamma_\mu      S_f(x_c,\,y;\,U,\,u_B) 
     e \gamma_\mu      S_f(y,\,x_c;\,U,\,u_B) 
    \right) 
   \right>_{u_B} 
  \right) 
 \right>_{u_A} 
\right>_U 
  \, , 
  \label{eq:method_equation}\end{aligned}$$ where the local form for the electromagnetic currents is assumed to avoid the complicate expression, and the wave function renormalization constant $Z_\mu$ of the muon is multiplied to remove the UV singularities in the vertex function of the muon.  The ordering of the averages in Eq.  (\[eq:method\_equation\]) is organized in such a way that the cancellation of ${\cal O}(\alpha_{\rm em}^2)$-”components” can be realized as accurately as possible at the level of each configuration $(U,\,u_A)$. The nested sum for $x,\,y$ could be efficiently carried out by working in the momentum space, where the Fourier transform of $D(x,\,y)$ is known exactly. 

The observation done around Eq.  (\[eq:O(alpha)-expansion\_first-term\]) relies on the QED perturbation, while the simulation treats QED nonperturbatively.  Thus, for our method to work, it is necessary that (a) the simulation can be done at the small coupling constant $\alpha_{\rm em}$, possibly at its physical value $\alpha_{\rm em} = \frac{1}{137}$, and that (b) the results for the weak coupling allow perturbative interpretation.  As for (a), the [*quenched non-compact*]{} QED allows us to get the uncorrelated $U(1)_{\rm em}$ gauge configurations efficiently [@Duncan:1996xy].  As for (b), the lattice (QCD + QED) study of the electromagnetic splitting of the pion masses demonstrates that the splitting is perturbative [@Yamada], confidence that our method for the h-lbl contribution will work. 

So far we have focused on the diagram consisting of a single quark loop with electromagnetic currents inserted.  In general, Fig.  \[fig:hadronic\_light-by-light\] contains two, three or four quark loops, where each quark loop contains at least one current.  In the quenched approximation to QED, they must be computed separately.  To illustrate implication of multi-quark loop diagrams, let us consider the charged pseudoscalar meson contribution [@hadronic_lbl] but with all electromagnetic currents put on its valence quarks, for simplicity.  Three types of two quark loop diagrams provide the meson diagrams where the currents lie on both of two valence quarks, while the single quark loop diagram can provide the diagram where all the currents lie on one of two valence quarks.  As is well-known, the calculation of multi-quark loop contributions (disconnected diagrams) is difficult in lattice calculations, and will be surely harder than the single quark loop diagram.  Actually, the multi-quark loop diagrams vanish in the exact flavor $SU(3)$ limit, except for a two quark loop diagram where four currents distribute over two quark loops evenly.  It is an important issue to examine if they are subdominant and safely neglected. 

We thank RIKEN, Brookhaven National Laboratory and the U. S. Department of Energy for providing the facilities and hospitality where this work was done.  M. H. ’s work is supported in part by Grant-in-Aid for Scientific Research (15740173) in Japan. 

[99]{} G. W. Bennett [*et al.*]{}  \[Muon $g-2$ Collaboration\], *Measurement of the positive muon anomalous magnetic moment to 0.7-ppm*, *Phys.  Rev. Lett.*  [**89**]{} (2002) 101804, *Erratum-ibid.*  [**89**]{} (2002) 129903 \[[hep-ex/0208001]{}\]; *Measurement of the negative muon anomalous magnetic moment to 0.7-ppm*, *Phys.  Rev. Lett.*  [**92**]{} (2004) 161802 \[[hep-ex/0401008]{}\].  N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, *The hierarchy problem and new dimensions at a millimeter*, *Phys.  Lett.* [** B 429**]{} (1998) 263 \[[hep-ph/9803315]{}\].  D. W. Hertzog \[E821 Collaboration\], *Measurement of the muon anomaly to high and even higher precision*, *Nucl.  Phys.  Proc.  Suppl.*  [**144**]{} (2005) 191 \[[hep-ex/0501053]{}\].  T. Blum, *Lattice calculation of the lowest order hadronic contribution to the muon anomalous magnetic moment*, *Phys.  Rev. Lett.*  [**91**]{} (2003) 052001 \[[hep-lat/0212018]{}\];\
M. Gockeler, R. Horsley, W. Kurzinger, D. Pleiter, P. E. L. Rakow and G. Schierholz \[QCDSF Collaboration\], *Vacuum polarisation and hadronic contribution to muon $g-2$ from lattice QCD*, *Nucl.  Phys.* [** B 688**]{} (2004) 135 \[[hep-lat/0312032]{}\].  C. Aubin and T. Blum, *Lowest order hadronic contribution to the muon $g-2$*, in proceedings of *XXIIIrd International Symposium on Lattice Field Theory*.  K. Hagiwara, A. D. Martin, D. Nomura and T. Teubner, *Predictions for $g-2$ of the muon and $\alpha({\rm QED})(M_Z^2)$*, *Phys.  Rev.* [**D 69**]{} (2004) 093003 \[[hep-ph/0312250]{}\];\
M. Davier, S. Eidelman, A. Hocker and Z. Zhang, *Updated estimate of the muon magnetic moment using revised results from $e^+ e^-$ annihilation*, *Eur.  Phys.  J.* [**C 31**]{} (2003) 503 \[[hep-ph/0308213]{}\].  M. Hayakawa, T. Kinoshita and A. I. Sanda, *Hadronic light by light scattering effect on muon $g-2$*, *Phys.  Rev. Lett.*  [**75**]{} (1995) 790 \[[hep-ph/9503463]{}\]; *Hadronic Light-by-light Scattering Contribution to Muon $g - 2$*, *Phys.  Rev.* [**D 54**]{} (1996) 3137 \[[hep-ph/9601310]{}\];\
J. Bijnens, E. Pallante and J. Prades, *Hadronic light by light contributions to the muon $g-2$ in the large $N_C$ limit*, *Phys.  Rev. Lett.*  [**75**]{} (1995) 1447, *Erratum-ibid.*  [**75**]{} (1995) 3781 \[[hep-ph/9505251]{}\]; *Analysis of the Hadronic Light-by-Light Contributions to the Muon $g-2$*, *Nucl.  Phys.* [** B 474**]{} (1996) 379 \[[hep-ph/9511388]{}\].\
M. Hayakawa and T. Kinoshita, *Pseudoscalar pole terms in the hadronic light-by-light scattering contribution to muon $g-2$*, *Phys.  Rev.* [**D 57**]{} (1998) 465, *Erratum-ibid.* [** D 66**]{} (2002) 019902 \[[hep-ph/9708227]{}\];\
M. Knecht and A. Nyffeler, *Hadronic light-by-light corrections to the muon $g-2$: The pion-pole contribution*, *Phys.  Rev.* [**D 65**]{} (2002) 073034 \[[hep-ph/0111058]{}\];\
M. Knecht, A. Nyffeler, M. Perrottet and E. De Rafael, *Hadronic light-by-light scattering contribution to the muon $g-2$: An effective field theory approach*, *Phys.  Rev. Lett.*  [**88**]{} (2002) 071802 \[[hep-ph/0111059]{}\];\
M. Hayakawa and T. Kinoshita, *Comment on the sign of the pseudoscalar pole contribution to the muon $g-2$*, hep-ph/0112102.\
I. Blokland, A. Czarnecki and K. Melnikov, *Pion pole contribution to hadronic light-by-light scattering and muon anomalous magnetic moment*, *Phys.  Rev. Lett.*  [**88**]{} (2002) 071803 \[[hep-ph/0112117]{}\];\
J. Bijnens, E. Pallante and J. Prades, *Comment on the pion pole part of the light-by-light contribution to the muon $g-2$*, *Nucl.  Phys.* [** B 626**]{} (2002) 410 \[[hep-ph/0112255]{}\];\
M. Ramsey-Musolf and M. B. Wise, *Hadronic light-by-light contribution to muon $g-2$ in chiral perturbation theory*, *Phys.  Rev. Lett.*   [**89**]{} (2002) 041601 \[[hep-ph/0201297]{}\];\
K. Melnikov and A. Vainshtein, *Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment revisited*, *Phys.  Rev.* [**D 70**]{} (2004) 113006 \[[hep-ph/0312226]{}\].  A. Duncan, E. Eichten and H. Thacker, *Electromagnetic Splittings and Light Quark Masses in Lattice QCD*, *Phys.  Rev. Lett.*  [**76**]{} (1996) 3894 \[[hep-lat/9602005]{}\].  N. Yamada, T. Blum, M. Hayakawa and T. Izubuchi, *Electromagnetic properties of hadrons with two flavor dynamical domain wall fermions*, in proceedings of *XXIIIrd International Symposium on Lattice Field Theory*, PoS(LAT2005)092. 

[^1]: For the unquenched QCD plus quenched QED to respect the gauge invariance of QED, the electromagnetic charges of sea quarks are assumed to be zero.