On axials and pseudoscalars in the hadronic light-by-light contribution to the muon $(g-2)$

Despite recent developments, there are a number of conceptual issues on the hadronic light-by-light (HLbL) contribution to the muon $(g-2)$ which remain unresolved. One of the most controversial ones is the precise way in which short-distance constraints get saturated by resonance exchange, particularly in the so-called Melnikov-Vainshtein (MV) limit. In this paper we address this and related issues from a novel perspective, employing a warped five-dimensional model as a tool to generate a consistent realization of QCD in the large$-N_c$ limit. This approach differs from previous ones in that we can work at the level of an effective action, which guarantees that unitarity is preserved and the chiral anomaly is consistently implemented at the hadronic level. We use the model to evaluate the inclusive contribution of Goldstone modes and axial mesons to the HLbL. We find that both anomaly matching and the MV constraint cannot be fulfilled with a finite number of resonances (including the pion) and instead require an infinite number of axial states. Our numbers for the HLbL point at a non-negligible role of the axial mesons, which is closely linked to a correct implementation of QCD short-distance constraints.


I. INTRODUCTION
The anomalous magnetic moment of the muon is one of the most precise tests of the Standard Model dynamics. Besides the dominant electromagnetic contribution, one is also testing the weak and strong interactions. The present experimental value is given by a exp µ = 116592091(54)(33) × 10 −11 [1], where statistical errors are the largest source of uncertainty. The upcoming experiments at FNAL [2] and J-PARC [3] are expected to reduced the experimental error by a factor four down to 1.6 × 10 −10 , much smaller than the current theoretical uncertainty.
If the E821 experimental number is confirmed, the discrepancy between the experimental value and the theoretical prediction, currently at about 3.5σ, would rise up to a 7σ effect with the projected new precision. It is therefore essential to have good control over the theoretical estimate.
The theoretical prediction for a µ is overwhelmingly dominated by electromagnetic and, to a much lesser extend, weak effects (see [4] for a review). Hadronic effects have a very modest contribution but are extremely difficult to evaluate. The present theoretical number, a SM µ = 116591823(1)(34)(26) × 10 −11 [1], is dominated by the hadronic uncertainties (second and third error sources). The largest hadronic contribution comes from the hadronic vacuum polarization, which can be rather cleanly connected to existing data on e + e − scattering [5][6][7]. In contrast, the (subleading) hadronic light-by-light contribution is more remote from experiment.
The physics involved in the HLbL contribution is sen-sitive to nonperturbative hadronic dynamics and cannot be calculated from first principles, except in some particular kinematical limits. One is therefore bound to use nonperturbative techniques. General arguments, based on chiral symmetry and the large-N c limit, can be used to assess the relevance of the different contributions [8].
There is general consensus that the neutral pion exchange provides the largest effect, and all determinations agree on its value. However, the status of the axial meson contribution and the subleading 1/N c effects is not as satisfactory. According to estimates, these contributions could be responsible for at least half the value of the HLbL, so a better understanding of them is definitely needed.
Progress on the light-by-light front is nowadays pursued along three main avenues: hadronic models, lattice simulations and dispersion relation approaches. Hadronic models provide, by far, the largest pool of HLbL determinations. In some cases the models are rather broad in scope [9][10][11][12], while in some others cases [13][14][15] the focus is on specific contributions. The main strategy behind these approaches is to come up with hadronic form factors able to successfully interpolate between low energies, where experimental data is available, and high energies, where the OPE of QCD is valid. The different models can then be understood as different ways to build interpolating functions between these limiting cases. In principle, the more constraints that a model satisfies, the more reliable their predictions should be. This has motivated a lot of work to increase the number of known short-distance constraints and test their impact on a µ [15,16].
Dispersion relation techniques [17][18][19][20] are an interesting alternative, where the goal is to bring the HLbL determination as close as possible to experimental data. However, at least at this point in time, it is not clear how these methods can implement the short-distance constraints from perturbative QCD. In turn, lattice simulations are rapidly becoming competitive [21] and should eventually give us the most precise determination.
However, the effort to bring the HLbL under better theoretical control also requires to address some conceptual issues that are still open. The core of the problem is to understand how short-distances are saturated by the different hadronic states. In the HLbL, understanding how this duality relations work turns out to be a highly nontrivial task.
One notorious example is the short-distance constraint discussed in [15], which was claimed to increase substantially the value for the HLbL contribution. Attempts to incorporate the constraint into form factor models have led to a number of proposals [4,22]. However, they all turn out to be problematic. Without a better understanding of how this constraint happens to be fulfilled, the state of affairs with the HLbL cannot be considered satisfactory.
In order to address the previous point, it is clear that one has to go beyond form factor parametrizations and be able to compute in terms of hadronic states at the level of correlators. This can be done if one borrows techniques from QFT in extra dimensions. It is well-known that, starting from a five-dimensional theory, the compactification to four dimensions gives rise to an infinite number of modes, which can be interpreted as mesons [23][24][25][26]. These constructions can be taylored to break chiral symmetry spontaneously, such that in the infrared limit one recovers chiral perturbation theory. If, additionally, the five-dimensional theory lives in an Anti-de-Sitter (AdS) gravitational background, the breaking of the associated conformal symmetry mimics the almost conformal behavior of QCD at large momenta. As a result, these theories approximate remarkably well both the long and short-distance behavior of QCD correlators. Finally, if the five-dimensional model is endowed with a Chern-Simons form, one obtains a four-dimensional theory with the chiral anomaly consistently implemented at all energy scales (see, e.g. [27]). Following the AdS/CFT prescription [28][29][30], there is a well-defined procedure to express the effective action in terms of external sources, out of which correlators can be computed with functional differentiation. The five-dimensional model is therefore used here simply as a technical device to end up with a consistent four-dimensional theory of hadrons.
In this paper we employ the simplest model implementing all of the abovementioned features to evaluate the Goldstone and the axial contributions to the HLbL correlator. As opposed to previous approaches, this determination is inclusive, in the sense that it contains all the axial excitations and unitarity is thus automatically built-in.
Our analysis clarifies a number of points. First, it shows that the Melnikov-Vainshtein condition can only be fulfilled by a collective effect of the axial resonances and not by a finite number of form factor. This explains, in particular, why the attempts to fulfill the condition with form factors were problematic, no matter their degree of sophistication. Second, this collective axial effect is intimately connected with having anomaly matching at all energies. With a finite number of axial mesons, anomaly matching simply fails. This shows that axial mesons have a more prominent role than previously assumed: in particular. Our number for the joint pseudoscalar and axial contribution to the HLbL is a µ = 127(17)·10 −11 , largely compatible with the findings of other hadronic models that incorporate the Melnikov-Vainshtein constraint. However, we argue that part of the axial mesons contribution has been misidentified as coming from pseudoscalars. As a result, the axial meson contribution to the HLbL is larger than previously estimated.
The structure of this paper is as follows. In Sect.II we introduce our toy model and bring it to the form of a four-dimensional effective action. We highlight how the pion and the (infinite towers of) vector and-axial vector resonances arise. In Sect.III we derive the expression for the HLbL electromagnetic tensor as predicted by the model in a closed form. The HLbL tensor is naturally split in two pieces, which collect the longitudinal and transverse polarizations of the resonances exchanged. The structure of the longitudinal piece is explored in detail in Sect.IV, and in Sect.V we show its relation with the triangle anomaly by analyzing the V V A correlator. The Melnikov-Vainshtein condition is discussed in Sect.VI. We show explicitly how the model fulfills it and make the connection with anomaly matching. In Sect.VII we give our numbers for the pion and axial meson contributions to the muon anomalous magnetic moment and compare them with previous estimates. Concluding remarks are given in Sect.VIII, while technical aspects are collected in three Appendices.

II. THE MODEL
In order to have a consistent realization of hadronic physics in the large-N c limit, we will adopt an extension of the five-dimensional model introduced in [25], which is a particular application of the AdS/CFT correspondence [28] to hadronic physics. The spirit of the model is to have a minimal setup able to capture the relevant features of QCD for the HLbL, namely a framework with conformality at very high energies, chiral symmetry breaking at low energies, and the chiral anomaly built in. Having the theory rooted in an underlying Lagrangian ensures that unitarity is satisfied and n-point functions are consistently interrelated. The AdS/CFT correspondence is useful to the extend that it provides the necessary toolbox to establish the connection between the five-dimensional formulation and the four-dimensional predictions.
The model is a five-dimensional U (3) L × U (3) R Yang-Mills-Chern-Simons theory, where the wedge product of forms is implicitly understood. Similar considerations apply to the right-handed sector. t a are the eight Gell-Mann matrices extended with t 0 = 1 3 / √ 6, normalized such that tr (t a t b ) = 1 2 δ ab . A point in the five-dimensional space has coordinates (x, z). The background metric will be chosen to be exactly Anti-de-Sitter, where µ, ν = (0, 1, 2, 3), M, N = (0, 1, 2, 3, z) and η µν has a mostly-negative signature. The fifth dimension is assumed to be compact, i.e. four-dimensional boundary branes exist at (x, 0) and (x, z 0 ), the so-called UV and IR boundary branes, respectively. In order to make contact with the hadronic states, it is convenient to trade the left and right-handed fields L µ (x, z) and R µ (x, z) for axial and vector ones through the usual relations L µ = V µ − A µ and R µ = V µ + A µ . These (massless) fields admit Kaluza-Klein decompositions and generate two infinite towers of four-dimensional modes, which become massive by absorbing the scalar modes V through higgsing. The resonance poles are determined from the solutions for ϕ V,A n (z). Working in four-dimensional momentum space, they are normalizable only for discrete values of the fourmomentum q, namely at where γ k,n is the n th root of the Bessel function J k (x). The previous equation shows that the size of the fifth dimension is an infrared quantity that sets the confinement scale.
Spontaneous chiral symmetry breaking is implemented in this model through boundary conditions on the IR brane, where low-energy physics takes place. The choice ensures that on the infrared brane only the subgroup SU (3) V is preserved. The pattern of breaking is therefore SU (3) L × SU (3) R → SU (3) V and an octet of Goldstone bosons is generated. The infrared boundary conditions make sure that all the zero modes cancel except A (0) 5 , which encodes the Goldstone degrees of freedom. In order to have a more conventional representation of the pions, it is convenient to trade the A such that the IR boundary conditions are respected, and redefines the fields as These chirally-dressed combinations make sure that the physical degrees of freedom are L ξ µ (x, z), R ξ µ (x, z) while their fifth components identically vanish.
The Goldstone degrees of freedom are captured by the combination of Wilson lines (ξ A (x, 0) ≡ ξ A (x)): which transforms as U (x) → g L (x)U (x)g † R (x). The field redefinitions in (8) are thus the way to move from a linear to a nonlinear representation of chiral symmetry breaking. With SU (3) L × SU (3) R , one can always choose ξ L (x) = ξ † R (x) ≡ u(x), such that U (x) = u 2 (x) [31]. As a result, the expressions for the chirally-dressed UV sources are The Yang-Mills action in eq. (1) is invariant under this field redefinition. In contrast, the Chern-Simons form gets shifted to where Σ L = dξ L ξ † L , and similarly for the right-handed fields. The second term can be shown to reproduce the ungauged Wess-Zumino-Witten Lagrangian, while the function is a pure boundary term in five dimensions. The connection with the associated effective fourdimensional theory is done with the AdS/CFT correspondence prescription [29,30], according to which the value of the five-dimensional fields on the UV brane are the sources of the four-dimensional operators. It is therefore convenient to split the fields as wherev(x) andâ(x) are identified with the classical sources associated with the chiral currents The functions a(x, z),ā(x, z), α(z), v(x, z) andv(x, z) can be found by solving the (linearized) five-dimensional equation of motions subject to the appropriate boundary conditions (see Sect. III).
In order to obtain the four dimensional effective action, the solutions (13) are substituted into the five dimensional action and the dependence on the fifth dimension is integrated out. The end result is a four-dimensional generating functional, out of which the correlators of the theory can be computed.
The model presented here contains an infinite tower of vector and axial-vector resonances together with the pion multiplet and is therefore a good toy model to evaluate the interplay between pions and axial-vectors in the HLbL. All observables are expressed in terms of three parameters, λ, c and z 0 . The former is normally fixed by matching the coefficient of the parton logarithm in the axial two-point function (see, e.g. [24] and Appendix B). This gives By requiring the right normalization of the chiral anomaly, one finds The parameter z 0 is a characteristic infrared scale. Actually, the simplicity of the model means that all infrared quantities depend on z 0 . For instance, the pion decay constant is given by and the resonance masses are given in eq. (5). Which parameter is chosen to fix z 0 depends on the application at hand. We will come back to this issue in Sect.VII. For the time being, we simply observe that the model predicts extremely well the splitting between the lowestlying vector and axial states, i.e. m ρ m a1 = γ 0,1 γ 1,1 ∼ 0.63 (18) but cannot account for satisfactory values for f π and m ρ simultaneously. This shortcomings can be circumvented but at the price of sophisticating the model in a rather ad hoc way. In this paper we will stick to the minimal model. As we will argue in the next sections, since not all the hadronic information is equally relevant for the HLbL, the minimal model is able to provide interesting quantitative estimates.
Another limitation of the model is that the Goldstone modes are strictly massless, i.e. there is no explicit chiral symmetry breaking. In order to be realistic, we should depart from the chiral limit and give the pion multiplet a mass. We will account for explicit chiral symmetry breaking simply by adding the pion mass in the pion propagator whenever necessary. This modification can be understood as introducing a deformation operator on the UV boundary and therefore does not jeopardize the consistency of the model. For more details, see for instance [32] and references therein.
This correlator satisfies the Ward identities: which reduce the number of independent kinematic invariants down to 43 gauge invariant tensor structures [19]. In our model, all quartic terms in vector fields are antisymmetric in flavor indices and therefore cancel due to Bose symmetry. The leading contributions are driven by the cubic interactions in the Chern-Simons term, corresponding to pion and axial-vector exchanges, which are the leading effects in the 1/N c expansion.
The corresponding diagrams are listed in fig. 2. The vertices can be extracted from the effective action, which is obtained by solving the equations of motion for axial and vector fields and plugging the solutions back into the five-dimensional action. As usual, exact solutions cannot be found and one has to resort to perturbation The HLbL diagram. The blob represent the HLbL tensor. In our conventions, momenta are pointing inwards.
theory, with the quadratic part of the Yang-Mills piece as the leading effect and the Chern-Simons term as a perturbation.
The Yang-Mills piece of the action can be written in components as From the Chern-Simons term we need to keep interactions linear in the axial field. They come from the first terms in eqs. (2) and (12). In components, one finds whered abc = tr t a {t b , t c } and the last piece comes from the boundary term in (12). Plugging the expressions for the axial and vector fields given in (13), the previous expression can be split into a transverse, longitudinal and the Goldstone contribution. In the following we will concentrate on the contributions of the neutral pion and the a 1 (1260) tower, and we will select accordingly the flavor structured 3γγ = 1 3 . For the numerical estimates in Sect. V, we will consider the contribution of the η, η ′ and the axial isosinglet towers. In order to built the diagrams of fig. 2 from the action one needs the first-order solutions for vector fields and up to second-order solutions for the axial ones. The first-order solutions are the so-called bulk-to-boundary propagators. Working in four-dimensional (Euclidean) momentum space, they read where K j (x) and I j (x) are modified Bessel functions. The axial zero-mode α(z) instead simplifies to The second-order solution for the axial field can be expressed in terms of the first-order ones as is the axial Green function (see Appendix A). The expression for the transverse and longitudinal components can be easily found by projecting out the corresponding components of the Green function, defined as Plugging the previous solutions back into the eq. (1), the relevant terms of the effective action to compute the electromagnetic tensor are: The first piece takes into account the contribution of axial vectors while the second line contains the form factor of the pion coupled to two photons, which is defined The pion propagator follows directly from the term which is the effective action for axial fields coming from eq. (21) once the expression in eq. (13) is used. Notice that the boundary term in (22) affecting the pion is no longer present. Its cancellation, leaving the expression above for the pion form factor, is a consistency check that the anomaly in our model is correctly implemented. The presence of a boundary term would be in conflict, e.g. with the OPE of this form factor at large photon virtualities. By matching (29) to the holographic expression above one finds The expression above depends on the parameters of the model, which can be traded for c, z 0 and f π . However, as opposed to other hadronic models, in our approach the functional form of F πγγ is completely fixed, so eq. (31) is actually a consistency check of the model. In the zero-momentum limit, F πγγ is determined by the chiral anomaly. Using that v(z, 0) = 1, the integral above is given by the boundary values for α(z). From eq. (16) one then obtains the well-known result which confirms that the model has the chiral anomaly correctly implemented. At very high energies, for large and equal photon momenta, one can expand (31) which matches the OPE prediction. When one photon is on-shell and the other far off-shell, one expects the Brodsky-Lepage Q −2 scaling. We indeed find Eq. (31) therefore has the correct high-and low-energy behaviour. Notice that only the leading term in the OPE is correctly reproduced by the model, with all subleading pieces identically vanishing. This is a consequence of the conformal symmetry of the AdS metric. The electromagnetic four-point function defined in (19) can now be obtained by taking the variation of (28) with respect to the boundary values of the (transverse) vector fields. Using (31), one obtains where s = (q 1 +q 2 ) 2 , t = (q 1 −q 3 ) 2 and u = (q 1 −q 4 ) 2 and we have used the shorthand notation The tensors T µ ij are defined as The closed expression of eq. (35) for the hadronic lightby-light tensor is all that is needed for the evaluation of the contribution to the anomalous magnetic moment. However, one of the virtues of having a consistent model with analytical control is that a number of issues can be examined in detail. This we will do in the following sections.

IV. LONGITUDINAL PIECE AND PION-EXCHANGE DOMINANCE
Based on large-N c arguments and dimensional power counting, there is agreement that the pion exchange contribution is the dominant piece of the HLbL.
Equation (35) contains, as expected, the pion contribution to the hadronic light-by-light tensor as the product of the πγγ form factors connected by a pion propagator. In turn, the first term on each line accounts for the contribution of the full tower of axial-vector states.
In order to understand better the structure of the HLbL tensor, it is convenient to project out its longitudinal and transverse parts, which can be done using eq. (26).
The longitudinal component, which also contains the pion contribution, can be expressed in terms of three tensorial structures, and the tensors T µνλσ and T µνλσ are the crossedsymmetric ones. Defining, as in [19], the crossing opera- The longitudinal form factors can be shown to take the simplified form The three one-pion exchange HLbL diagrams. The dashed line denotes the pion propagator and the dots the Fπγγ form factors. Photon momenta assignments are the same as in fig.1, i.e., they point inwards.
with similar expressions for W 13;24 and W 14;23 . The second line above can be easily obtained from the first term in eq. (35) by using the antisymmetry of the Levi-Civita tensors and integration by parts. Using the expression (see Appendix A) the longitudinal form factor takes the form where we have used that z 2 0 = 8λf −2 π and the definition of the pion transition form factor in eq. (31).
The previous expression shows explicitly that the contribution of the whole tower of axial states consists of a factorizable and a nonfactorizable piece in five dimensions. In four dimensions, they correspond to a propagating piece and a contact term, respectively. In the chiral limit, one can easily check that the axial propagating piece and the pion contribution cancel each other and one is left with the contact term. Explicitly, The structure of this expression, and in particular the presence the contact term, is mostly dictated by the chiral anomaly, as we will show more explicitly in the next section. At very low energies (still in the chiral limit) the integral appearing in the contact term can be easily evaluated. Using that v(z, 0) = 1 and the explicit expression for α(z), one finds The result is actually the same that one would have obtained by dropping the axial tower and considering only the pion exchange contribution. This is of course not a coincidence. At very low energies only the pion is a dynamical degree of freedom and it is entirely responsible to fulfill the chiral anomaly. This is the content of the Wess-Zumino-Witten term in chiral perturbation theory, which our model also contains, and actually fixes the value of F πγγ (0, 0), as we have already shown. It is clear that at higher energies anomaly matching requires the participation of resonance states other than the pion. However, since no first-principle description of the strong interactions exists in the intermediate energy regime, it is not known how this is implemented in detail. The result in eq. (43) is precisely the way the model implements the chiral anomaly in a consistent way at all energy scales. The expression for the resumed axial contributions and the precise cancellation of the pion contribution in the chiral limit can be therefore seen in this light as a sort of sum rule to enforce anomaly matching at all energies. This interpretation will be reinforced in the following section, where we will evaluate the effects of the chiral anomaly in a more explicit fashion.

V. ANOMALY MATCHING IN THE V V A CORRELATOR
The best way to uncover the role of the axial anomaly in the HLbL tensor is to look at the three-point correlator This correlator has been studied in great detail in e.g. [34,35]. In [35] it is shown that, on general grounds, Γ µνλ can be decomposed in four independent tensorial structures: one longitudinal and three transverse, which are associated with four different scalar functions. In the holographic model, the V V A correlator can be computed from the variation of the effective action (22) with respect to the sources. Rewriting the action for transverse and longitudinal axial sources, integrating by parts and using the boundary conditions for the axial and vector fields the result is (S where the first term in each line is a boundary term, while the second one has a nontrivial profile and energy  46) and (47). The diagram on the right is the pion exchange contribution.
dependence. Apart from direct vertices, the correlator also contains a pion contribution (see fig. 4). Explicitly, the pion contribution comes from where the first line comes entirely from eq. (30) by using the decomposition of eq. (13). From the previous terms in the effective action it is straightforward to obtain the expression for Γ µνλ . The longitudinal part of the correlator yields Γ µνλ (q 3 , q 4 ) = t µνλ 2c 3q 2 where t µνλ is the longitudinal tensor defined as in [35], The transverse part is less straightforward to get. The reason is that the model necessarily describes the consistent anomaly, which is the only one that can be derived from an action [36]. However, in the presence of gauge fields, only the covariant anomaly is compatible with the Ward identities. This change of prescription can be interpreted as a different definition of the chronological product T . Therefore, the result that comes out of the effective action in (22) is not (45) but Γ µνλ (q 1 , q 2 ) = i d 4 xd 4 y e i(q1·x+q2·y) whereT is associated with the consistent anomaly. In this prescription, the photon contains both transverse and longitudinal components. Taking this into account, one can check that the Ward identities are: which are indeed the ones corresponding to the consistent anomaly.
The relation between both prescriptions involves a Bardeen-Zumino polynomial of the form [36] such that in the covariant prescription one recovers the well-known Adler-Bardeen results In the particular kinematic configuration where q 4 and (q 3 +q 4 ) 2 go to zero, the four independent tensorial structures reduce to just two [35], In this limit there are two independent kinematic invariants, defined as The longitudinal function ω L is known to be fixed by the anomaly to This result is exact in the chiral limit, and gets corrections only from non-perturbative contributions. ω T instead depends on the dynamics of axial meson exchange. Both functions are however linked at very high energies, where they satisfy the well-known expression [34] Given the discussion of the previous section, it should come as no surprise that, in the chiral limit, the previous expression shows a cancellation between the pion contribution and the energy-dependent part of the longitudinal vertex. In the chiral limit, therefore, the longitudinal part is saturated by the boundary term. Notice that this cancellation between the whole tower of axial states and the pion, which is naturally implemented by the model, is the only way to have at all energies. The cancellation of the pion contribution against a collective effect of the whole axial tower (in the chiral limit) is thus a sum rule enforced by the chiral anomaly.
As constructed in (54), the Bardeen-Zumino term cancels the boundary terms in the transverse part and brings the longitudinal part to the Adler-Bardeen value. In this prescription, the predictions for ω L and ω T in Euclidean The integral in ω T can actually be solved analytically and its expression considerably simplified. In order to do so, it is convenient to rewrite [37] a(z, The first piece is a boundary term while the second one can be shown to be linear in z (see eq. (A8)). The result for both form factors can therefore be written in the rather compact form: As expected, in the chiral limit the pion contribution gets cancelled by part of the axial one such that ω L is structureless, as already noticed before. We stress that this is a consistency check that the anomaly is correctly implemented in the model. The study of the V V A correlator thus reveals that the non-factorizable piece in the HLbL tensor that we observed in the previous Section has the same origin as the contact term in the triangle. Expanding the previous expressions for large momenta, one finds which implies that This scaling is consistent with the OPE. However, as already emphasized, conformal symmetry in the fivedimensional model ensures that the OPE relations will be satisfied to leading order, but the effects of OPE condensates will in general be missed. Therefore, the expression above has the value of a consistency check rather than a prediction.
Another important observation is that the pion contribution at low energies happens to be lim q→0 ω (π) and one would be tempted to conclude that the pion saturates the anomaly at low energies. This interpretation is not wrong but it is misleading: in the m π → 0 limit, the pion and the dynamical axial contribution cancel analytically at all energies. Based on the previous derivation, one is forced to conclude that the pion does not saturate ω L , although it is fundamental to make sure that the result is consistent with the anomaly.
where a π is the slope of F πγγ at zero momentum (see eq. (90)). The predictions of the present model can be compared, for instance, with the expressions used in [15]: The results are illustrated in figs. 5 and 6, where one can see that the main difference between both models mostly affect ω T .

VI. THE MELNIKOV-VAINSHTEIN LIMIT
The closed expression that we obtained for the lightby-light tensor can be extrapolated to large Euclidean momenta and compared with the different short-distance constraints derived from the OPE of QCD. This is a necessary consistency check before the numerical analysis of the next section. The more constraints the expression satisfies, the more reliable the predictions will be. Conversely, if a constraint is not satisfied, one can estimate its impact by comparing with models that do implement it.
Relevant for the evaluation of the (g−2) µ are the limits when q 2 4 = 0. A non-trivial check is to find out how the HLbL behaves when all virtual photons have large momenta, For the longitudinal component of the HLbL tensor, the quark loop diagram in QCD gives [15] One can find the corresponding expression in the model by applying the high-energy limit to eq. (42). The first thing to note is that this limit cannot be fulfilled by the pion contribution, which falls off like Q −6 . The relevant piece comes instead from the axial tower. This is the consistent with what we observed in the previous sections. The leading term comes from the first piece in eq. (42), which gives Numerically, this amounts to 80% of the OPE coefficient. Particularly interesting are the limits where only two photons have large virtualities. Without loss of generality one can consider Q 2 1 ≃ Q 2 2 ≫ Q 2 3 ≫ Λ 2 QCD , with the remaining two possibilities generated by crossing symmetry. In the context of the (g − 2) µ , these limits were first explored in [15], where the relation with the anomaly was emphasized.
Consider the product of two electromagnetic currents: in the kinematical limit where ξ is large and all momenta are space-like. In this limit, the OPE gives 5ρ (z) =qQ 2 γ ρ γ 5 q. This limit shows that a number of the short-distance constraints relevant for the evaluation of the HLbL are actually determined by the axial anomaly.
Within the model, the previous OPE can be computed in the following way. It is rather straightforward to check that, when the limit Q 2 1 ≃ Q 2 2 ≫ Q 2 3 ≫ Λ 2 QCD is taken, the leading term comes from the first line of eq. (35), i.e. from Figs. (2.a) and (2.b). The crossed diagrams are subleading, as expected.
The pion contribution can be computed using the asymptotic expression (33), from which one concludes that it falls off with higher powers of Q 2 than the previous OPE demands. The contribution of the axial tower is the relevant piece. For its evaluation it is convenient to consider the integral In order to separately analyze the longitudinal and transverse components, we will define In the limit Q ≫ Q 3 , one can easily show that The fact that the longitudinal piece is exact to all orders in Q 3 is the manifestation of the chiral anomaly. For the softer part of the integral, using that one easily concludes that Plugging the previous expressions back into the HLbL tensor, one finds where in the second line we have identified ω L and ω T using their expressions in eqs. (60). The former displays the cancellation between the pion and axial contributions that we have observed before. We already emphasized the role of the contact term in the chiral limit. This term is actually the one that enforces the right expression for ω L and, accordingly, makes sure that the chiral anomaly is correctly implemented. Notice that the pion contribution alone is incompatible with both requirements, i.e. ω (π) This equation is clearly incompatible with ω ∼ Q −2 3 in the chiral limit. The problem is the structure of the form factor, which depends on Q 3 . For this same reason, it is clear that no single particle exchange can saturate ω L . In order to satisfy the constraint an infinite number of (axial) particles is needed. This is precisely what the contact term is indicating.
We note that the mechanism to saturate the Melnikov-Vainshtein constraint found above depends only on the correct implementation of the chiral anomaly and is therefore rather generic. In particular, this indicates that attempts to reproduce ω L with pseudoscalar form factors [4,22] are not realistic. Our result also implies that with experimental input on low-energy states alone, the MV constraint will be missed. In the next Section we will see that this has a substantial numerical impact on the HLbL, thus confirming the main conclusion of [15].

VII. NUMERICAL ANALYSIS
In eq. (35) we already wrote down the electromagnetic HLbL tensor in a closed form within our model. In order to perform our numerical analysis and be able to compare with other studies, we will employ standard modelindependent techniques.
Quite generically, the HLbL tensor can be expanded as a sum over gauge invariant Lorentz tensor structures. Out of the 138 structures [33], once Ward identities are imposed, one ends up with 43 kinematic structures. Here we will adopt the formalism introduced in [19], which builds on previous works [38,39]. The number of kinematical structures is extended such that they are free of poles and zeros. We will therefore write where the definitions of the tensors T µνλσ i can be found in [19].
Using projection techniques, the two-loop diagram of fig. 1 can be related to the muon anomalous magnetic moment as follows: where p is the muon momentum. Due to the projection above, only 19 independent linear combinations of the 54 T µνρλ i contribute to a HLbL µ [33]. Furthermore, due to the symmetries of the two-loop integral, one needs to evaluate eventually only 12 different integrals containing 12 scalar coefficientsΠ i (q 1 , q 2 , q 3 ).
Following the general analysis in [14,19], one can perform five out of the eight integrals above using Gegenbauer polynomials, regardless of the specific form ofΠ i . The resulting master formula contains then only three integrals and, in terms of Euclidean momenta, takes the form: where Q 1 and Q 2 are the radial components of the momenta. The hadronic scalar functionsΠ i are evaluated for the reduced kinematics The complete list of the integral kernelsT i (Q 1 , Q 2 , τ ) can be found in Appendix B of Ref. [40].

A. Longitudinal contributions
Given our previous discussion, it is convenient to split the contributions in longitudinal and transverse parts. One can check that the longitudinal contribution is described by the first two structures in eq. (83). They correspond to eq. (37), i.e.
and similarly for the crossed-symmetric T (2) µνλσ and T µνλσ . The scalar invariants were simplified in Sec. IV to the form where the first line accounts for the pion contribution and the remaining two are the resummation of the whole axial tower. In Sec. IV we showed that our model generates a pion transition form factor which has the correct intercept at zero momentum and the right scaling when one or the two photons have large momenta, namely However, in Sec. II we also emphasized that our model predicts which is a factor 2 bigger than the experimental number. This is a well-known shortcoming of the present model, which can be amended with more sophisticated versions of it. Such sophistications are beyond the scope of the present paper. In the following we will make two choices for the parameters c, λ and z 0 , which emphasize different energy regimes in the HLbL. This will be used as an estimate of the uncertainty in our determination of a µ .
Since the longitudinal contribution to HLbL is the dominant one, and in particular the pion form factor plays a prominent role, a reasonable criterium is to choose the parameters such that agreement with F πγγ is achieved. Based on eq. (88), fixing N c and f π to the physical values thus seems to be the reasonable choice, at the price of overshooting m ρ .
However, experimental information exists also for the slope of the form factor at low momentum a π , defined as Using the low-energy expansion for v(z, Q), one readily finds [41] a π = −m 2 (92) which is in excellent agreement with the current world average, (a π ) exp = 0.0335 (31), if one fixes z 0 to match the physical m ρ . The result is largely undershooted if one fixes z 0 with f π . In the following, we will thus consider Additionally, we will take α em = 1 137.036 ; m µ = 105.7 MeV; m π = 135 MeV (94) The first set ensures the right behaviour of F πγγ at low energies (intercept and slope). An additional perk of fixing m ρ to the physical value is that the axial multiplet masses fall into the right ballpark. At high energies, it still correctly reproduces the Q 2 fall-off behaviour, yet it fails to pin down the right coefficients. The second set matches the expected asymptotic Q 2 behaviour but gives a poor determination of m ρ and the F πγγ slope. Given the form of the kernels in the expression for a µ , which are peaked in the low GeV regime, one would be tempted to prioritize the low-energy regime. However, this depends on how fast the asymptotic regime sets in. In figs. 7 and 8 we compare the predictions for the pion form factors with both sets of parameters for one virtual photon and for both photons virtual with equal momenta, respectively. The lattice result of Ref. [42] is also included for comparison. Both parameter sets give predictions that are within the experimental errors for the singlyvirtual photon case (not shown in the plot), so there is a priori no reason to prefer one set of parameters over the other. Having the right slope is important, but since the a µ measures an integral involving the form factor, global aspects are also important. In the following, we will report our numbers for both sets. Our final number will be the average of them.
For the pion contribution we find where the left and right numbers stand for the Set 1 and Set 2 predictions. This is in excellent agreement with previous determinations based on form factor analyses, e.g.
The contributions of the isoscalar pseudoscalars will be estimated by simply using the physical values for masses and decay constants, These numbers are not the ones predicted by the model, which has an exactly massless Goldstone nonet with a common decay constant. The introduction of masses can be argued exactly as we did with the pion: their effects are limited to the pseudoscalar propagators without affecting the dynamics of the model. From the fivedimensional perspective, this can be achieved with the introduction of flavour-dependent boundary terms. For our purposes, we will simply add the different masses by hand. Breaking the degeneracy of the decay constants for the η and η ′ can also be done as long as one consistently correlates it with the corresponding isoscalar axial  channels, such that the sum rules that preserve anomaly matching remain in place. In practice, this can be done by modifying the parameter sets to for the η ′ and axial f * 1 (1420) tower, and similarly for the η and axial f 1 (1285) tower.
Our results are The transverse terms of the HLbL tensor collect the remaining axial contributions. They represent only a small correction to the HLbL value, but they are considerably harder to determine. Explicitly, they follow from with the first two scalar factors, which are the longitudinal contributions, subtracted. The expressions for the scalar functions in our model are given in Appendix C.
The final values one obtains for the whole towers of isovector and isoscalar axial states is for both sets of parameters. As discussed in the previous Section, care has to be taken to use the right input parameters for each flavour.

C. General discussion
Averaging out the results from the two sets of parameters and using the spread as an estimate of the uncertainty, our final number for the contribution of Goldstone modes and axial states is a HLbL µ = 12.7(1.7) · 10 −10 (107) The breakdown of the different contributions to this number is collected in Table I. This number can be compared with previous literature on the same contributions, e.g.
showing excellent agreement. However, care has to be exercised when analyzing the relative weight of pseudoscalars and axials in the numbers above. As we already discussed, our model shows that the pseudoscalars do not contribute to the Melnikov-Vainshtein limit. A more meaningful comparison would be to compare the longitudinal and transverse contributions. One would then find a L µ = 11.3(1.7) · 10 −10 a HLbL However, we note that when one splits the particle contributions apart from the numbers above we find a PS µ = 9.6(1.6) · 10 −10 ; a A µ = 3.0(0.2) · 10 −10 (113) In other words, the relative weight of axials is substantially bigger than previously claimed in the literature. This does not mean that previous analyses should increase their estimate for axials, as long as the Melnikov-Vainshtein constraint is satisfied. However, analyses of the lowest-lying axial contributions [49,50] are grossly underestimating the axial contribution. Our analysis shows that, since the axials are mostly constrained by the anomaly, the uncertainty in their contribution turns out to be very small. The bulk of the associated uncertainty of our number comes from the pseudoscalar contribution, which can only be improved with better experimental data. However, this does not mean that our uncertainty on the axials is realistic. As already shown in [15], the splitting of axial masses inside the multiplet (something that our model cannot reproduce) can have a sizeable effect. This effect cannot be estimated in our model, but the observation recommends to increase the axial uncertainty.
We conclude that the role of axials is larger than previously claimed. However, their effects have been properly accounted for in those calculations where the Melnikov-Vainshtein limit is enforced. They have just been misrepresented as pseudoscalar contributions. In order to avoid confusion, it is more accurate to speak about longitudinal and transverse contributions to the a HLbL µ , which are well-defined even in form factor analyses.

VIII. CONCLUSIONS
The evaluation of the muon hadronic light-by-light contains a number of conceptual issues which are hard to address using the approaches used so far. Particularly challenging is the phenomenological implementation of the so-called Melnikov-Vainshtein limit. The problem is really how to match OPE constraints with resonance exchanges.
A framework suitable to study these issues should be able to (a) evaluate hadronic effects from a Lagrangian formalism while (b) being able to reproduce the right high-energy limits of QCD correlators. In other words, one would need a consistent realization of hadronic physics at the level of correlators.
This can be done if one starts from a Lagrangian formulation in five dimensions and integrates out the fifth dimension. The spectrum of the resulting effective fourdimensional action contains an infinite number of resonances, with the quantum numbers of the fields introduced in the initial Lagrangian. In this paper we have chosen a minimal version of such constructions. The resulting four-dimensional effective action has a number of interesting features: (a) it is a theory of resonances with unitarity built-in; (b) the anomaly is consistently implemented at the hadronic level, i.e. at all energies; (c) the high-energy limit of correlators matches the pQCD predictions, such that duality is correctly implemented; and (d) it generates a phenomenologically successful pion transition form factor. With this toy model we have evaluated the contributions of pseudoscalar and axial resonances in an inclusive way. We have thereby clarified why the implementation of the Melnikov-Vainshtein limit does not look natural in a form factor analysis: the limit results from a collective effect of the axial resonances, and accordingly cannot be reproduced with a finite number of states. In particular, it does not result from heavy pseudoscalar exchange, as initially hinted in [15] and recently attempted in [22].
We have shown explicitly how the sum rule that enforces the Melnikov-Vainshtein limit is the same that implements the anomaly in the V V A triangle through a nontrivial interplay between pseudoscalars and the whole tower of axials. Our analysis also clarifies that the difference between the pion-pole and the pion-exchange contributions to the HLbL is simply the longitudinal axial contribution. We believe that it would be more accurate to describe the contributions to the HLbL as longitudinal and transverse, as we have done in this paper. where the uncertainty is simply orientative. This number is in agreement with those determinations that implement Melnikov-Vainshtein constraint and confirms the relevant weight of this constraint in the HLbL.

Our final number for the HLbL is
Note added: During the completion of this work Ref. [51] appeared. While the focus of [51] is different from the one taken in this paper, there is overlap with our work as to the role of axials. Our analysis agrees with their conclusions.
In all theΠ i above, v(z, q 2 4 ) (or v(z ′ , q 2 4 ) ) always appear without derivative, and since we are taking the q 4 → 0 limit, if follows that v(z, q 2 4 ) → 1. Thus, the integral at the soft-photon vertex contains the product of the transverse axial Green function and only one derivative both depending on the same four-momentum. This leads to simplifications.