11institutetext: Albert Einstein Center for Fundamental Physics, Institute for
Theoretical Physics, University of Bern, Sidlerstrasse 5, 3012 Bern,
Switzerland 22institutetext: Paul Scherrer Institut, 5232 Villigen PSI,
Switzerland 33institutetext: University of Vienna, Faculty of Physics,
Boltzmanngasse 5, 1090 Vienna, Austria

#  PSI-PR-21-13, UWThPh 2021-7
Short-distance constraints for the longitudinal component
of the hadronic light-by-light amplitude: an update

G. ColangeloBern  F. HagelsteinPSI  M. HoferichterBern  L. LaubBern  P.
StofferVienna

###### Abstract

We reassess the impact of short-distance constraints for the longitudinal
component of the hadronic light-by-light amplitude on the anomalous magnetic
moment of the muon, $a_{\mu}=(g-2)_{\mu}/2$, by comparing different solutions
that have recently appeared in the literature. In particular, we analyze the
relevance of the exact axial anomaly and its impact on $a_{\mu}$ and conclude
that it remains rather limited. We show that all recently proposed solutions
agree well within uncertainties on the numerical estimate of the impact of
short-distance constraints on $a_{\mu}$, despite differences in the concrete
implementation. We also take into account the recently calculated perturbative
corrections to the massless quark loop to update our estimate and outline the
path towards future improvements.

††journal: Eur. Phys. J. C

##  1 Introduction

The recent measurement of $(g-2)_{\mu}$ by the Fermilab Muon $g-2$
collaboration Abi:2021gix ; Albahri:2021ixb ; Albahri:2021kmg ;
Albahri:2021mtf  has made the discrepancy with the Standard-Model prediction
Aoyama:2020ynm ; Aoyama:2012wk ; Aoyama:2019ryr ; Czarnecki:2002nt ;
Gnendiger:2013pva ; Davier:2017zfy ; Keshavarzi:2018mgv ; Colangelo:2018mtw ;
Hoferichter:2019gzf ; Davier:2019can ; Keshavarzi:2019abf ; Hoid:2020xjs ;
Kurz:2014wya ; Melnikov:2003xd ; Masjuan:2017tvw ; Colangelo:2017qdm ;
Colangelo:2017fiz ; Hoferichter:2018dmo ; Hoferichter:2018kwz ;
Gerardin:2019vio ; Bijnens:2019ghy ; Colangelo:2019lpu ; Colangelo:2019uex ;
Blum:2019ugy ; Colangelo:2014qya  more serious, by bringing it from the
$3.7\sigma$ to the $4.2\sigma$ level when combined with the previous
Brookhaven measurement bennett:2006fi . The concrete perspective of additional
reductions of the experimental uncertainty in the near future—mainly from
subsequent Runs at Fermilab Grange:2015fou , but also the future J-PARC
experiment Abe:2019thb  using a different technique—makes the need of further
theoretical improvements more urgent. As is well known, the two main sources
of theoretical uncertainties are both hadronic. The largest one is the
hadronic vacuum polarization (HVP) contribution, and the second hadronic
light-by-light (HLbL) scattering. A lot of work has been devoted to reducing
the uncertainty in the latter by separately analyzing each of the different
contributions to the HLbL amplitude. These involve different intermediate
states and their calculation requires a good understanding of the relevant
physics. The present status has been summarized in the White Paper (WP) on the
Standard Model prediction of $(g-2)_{\mu}$ Aoyama:2020ynm , resulting in a
phenomenological estimate in agreement with lattice QCD Blum:2019ugy ;
Chao:2021tvp . The behavior of the HLbL amplitude for asymptotic values of its
arguments is fixed by QCD and represents an important, global constraint,
which has a significant impact on the estimate of this contribution. A
detailed understanding of which intermediate states play a role in satisfying
this constraint is crucial to estimate its impact and reduce the overall
theoretical uncertainty.

There are two different regimes of asymptotic momenta and correspondingly two
different constraints. For $g-2$ kinematics, with one photon in the static
limit, the HLbL amplitudes only depend on the squared momenta of the remaining
three photons. The first regime is when two of them are much larger than the
third, whereas all of them being about equally large defines the second. They
represent two limiting cases of a possible continuum of short-distance
constraints (SDCs) and we will refer to them as SDC1 ($q_{1,2}^{2}\gg
q_{3}^{2}$, $q_{1,2}^{2}\gg\Lambda^{2}_{\text{QCD}}$) and SDC2 ($q_{1}^{2}\sim
q_{2}^{2}\sim q_{3}^{2}\gg\Lambda^{2}_{\text{QCD}}$). Melnikov and Vainshtein
(MV) were the first to derive SDC1 Melnikov:2003xd  and, in particular, to
point out that in the chiral limit and for asymptotic values of $q_{1,2}^{2}$,
the leading $1/q_{3}^{2}$ behavior of the longitudinal part receives no
corrections neither at large nor at small $q_{3}^{2}$ values: in other words,
the $1/q_{3}^{2}$ dependence is exact across the whole range. This is a direct
consequence of the axial anomaly Adler:1969gk ; Bell:1969ts ; Bardeen:1969md ;
Wess:1971yu ; Witten:1983tw , see Sect. 2. The transverse part in turn is
constrained by the celebrated non-renormalization theorems Vainshtein:2002nv ;
Knecht:2003xy  for the vector–vector–axial-vector (VVA) correlator. The SDC2
case was also discussed by MV on the basis of the quark loop, but its
derivation has only recently been put on a firm basis by using the operator
product expansion (OPE) Bijnens:2019ghy . Moreover, both non-perturbative
Bijnens:2020xnl  and perturbative corrections Bijnens:2021jqo  to the OPE have
recently been calculated, thereby reducing one source of uncertainty.

There has been much interest in finding a way to satisfy these SDCs beyond the
model solution discussed by MV Melnikov:2003xd . We proposed a Regge model of
pseudoscalar resonances Colangelo:2019lpu ; Colangelo:2019uex , whereas a
solution based on the resummation of a tower of axial-vector resonances in a
holographic model of QCD (hQCD) was put forward in two more recent papers
Leutgeb:2019gbz ; Cappiello:2019hwh . A completely different approach based on
a set of interpolants between long and short distance has been adopted by
Lüdtke and Procura (LP) Ludtke:2020moa . Some of these works either appeared
or were published after the WP, where the estimate about the impact of the
SDCs and of the axial-vector contribution is significantly lower than what was
estimated in Melnikov:2003xd . While none of the most recent papers has
criticized the estimate in the WP, there are statements in Cappiello:2019hwh
that their results also agree with those in Melnikov:2003xd , and with
Leutgeb:2019gbz , who in turn conclude that the MV model is not the correct
way to implement SDC1, which makes the whole situation rather confusing. Given
the relevance of the SDCs, which currently represent the largest contribution
to the theoretical uncertainty of $a_{\mu}^{\text{HLbL}}$ Aoyama:2020ynm , it
is important to understand the differences between these solutions, clarify to
what extent they agree and where exactly differences arise, and reassess the
current situation.

From a theoretical point of view, the solution based on hQCD is particularly
relevant and appealing as it represents the first hadronic model of QCD based
on axial vectors that exactly satisfies the axial anomaly and SDC1 in the
chiral limit. As we will discuss in Sect. 3 on the basis of general, model-
independent arguments, the solution has to arise from a resummation of an
infinite tower of axial vectors, as it does in the hQCD models and as is
expected when fulfilling SDCs with hadronic states Peris:1998nj ;
Bijnens:2003rc . The model also has the advantage that the resummation can be
performed analytically, but its simplicity comes at the price of lack of
flexibility: once a number of inputs is used to pin down the free parameters
in the model, any further quantity can be predicted and shows some
discrepancies with QCD phenomenology. In particular, in the simplest models on
which we will focus here, the fulfillment of the asymptotic constraints in
general generates tensions with phenomenological low-energy constraints
Leutgeb:2019gbz ; Cappiello:2019hwh .

The solution originally proposed by MV also exactly satisfies the axial
anomaly and SDC1 in the chiral limit, but achieved this goal by a mere
truncation: every hadronic contribution beyond the pion pole for $g-2$
kinematics was simply dropped. Such an approximation is very well justified
for the three-point function $\langle VVA\rangle$, as explicitly shown in
Cappiello:2019hwh . For what concerns the HLbL amplitude the situation is
different: the MV model extrapolates the OPE expression to low $q_{1,2}^{2}$,
where it cannot be justified. It was first pointed out in Colangelo:2019lpu ;
Colangelo:2019uex  that the largest contribution to $a_{\mu}$ in the MV model
comes from the low-energy region, where additional intermediate states would
contribute. The first main point of this paper is then to demonstrate that
this conclusion applies to all the recently proposed implementations,
summarized in Sect. 4, explaining why there is general consensus on the
numerical impact despite significant differences in the implementations
themselves, see Sect. 5.

For instance, the hQCD and Regge models differ in their use of pseudoscalar
vs. axial-vector states. In Colangelo:2019lpu ; Colangelo:2019uex  the main
motivation for adopting a Regge model of pseudoscalar resonances is related to
a peculiar property of their contribution to the HLbL tensor. In a dispersive
approach Colangelo:2017qdm ; Colangelo:2017fiz ; Hoferichter:2013ama ;
Colangelo:2014dfa ; Colangelo:2014pva ; Colangelo:2015ama  the contribution of
narrow-width resonances to the HLbL tensor is in general ambiguous as it
depends on the basis in which the calculation is performed, unless a set of
sum rules is satisfied. Only in the case of pseudoscalars are these sum rules
automatically satisfied. The drawback of using pseudoscalar resonances is that
in the chiral limit they do not couple to the axial current and therefore
cannot contribute to the anomaly. We argued in Colangelo:2019lpu ;
Colangelo:2019uex  that such a model would nonetheless represent a useful tool
to make a realistic evaluation of the impact of the SDCs in the physical
world, i.e., away from the chiral limit.

This view has been challenged by MV in Melnikov:2019xkq : they emphasized the
importance of the axial anomaly as an exact constraint in the chiral limit,
and considered its fulfillment essential in order to make a reliable estimate
of the impact of the SDCs on $a_{\mu}$. By a detailed comparison between the
MV and the hQCD model we will show, however, that the relevance of the exact
axial anomaly in determining the four-point function is limited to a kinematic
region whose impact on the calculation of $a_{\mu}^{\text{HLbL}}$ is very
small. This is one of the most important conclusions of this paper, which
extends and confirms the findings in Colangelo:2019lpu ; Colangelo:2019uex ,
in line with earlier studies of the relevant momentum regions Bijnens:2007pz .

Since each of the models discussed here cannot claim to be a faithful
representation of QCD but at best be a tool to capture the essential features
thereof in connection with a particular aspect of the $a_{\mu}$ calculation,
it is instructive to compare all three of them, even if they rely on different
degrees of freedom to fulfill the SDCs. Anticipating our conclusions, we will
find a satisfactory agreement between the hQCD and the pseudoscalar Regge
models. We will then use the latter to update our earlier estimate taking into
account the recently calculated perturbative and non-perturbative corrections
to the OPE Bijnens:2020xnl ; Bijnens:2021jqo , see Sect. 6. This serves only
to illustrate the current status, because we believe that it is possible to
incorporate the good theoretical properties of hQCD models into our dispersive
formalism, after developing a coherent formulation of axial vectors in the
narrow-width approximation. This is the direction in which future work will
evolve, as we will sketch in the outlook in Sect. 7.

##  2 The longitudinal OPE and non-renormalization theorems

We concentrate here on the OPE for the longitudinal amplitude, which concerns
only one of the functions in the HLbL tensor, namely the $\hat{\Pi}_{1}$
function introduced in Colangelo:2017fiz . Only this function contains the
contribution of the pion pole, which can be written as111We work with
$q_{4}^{2}=0$, but $q_{4}\neq 0$ and, for simplicity, only consider the
isospin-triplet component for now. The Mandelstam variables are defined as
$s=(q_{1}+q_{2})^{2}$, $t=(q_{1}+q_{3})^{2}$, $u=(q_{2}+q_{3})^{2}$.

$\displaystyle\hat{\Pi}_{1}(q_{1}^{2},q_{2}^{2},q_{3}^{2},0;s,t,u)$
$\displaystyle=\frac{F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})F_{\pi\gamma\gamma^{*}}(q_{3}^{2})}{s-M_{\pi}^{2}}$
$\displaystyle+\tilde{G}(q_{1}^{2},q_{2}^{2},q_{3}^{2},0;s,t,u)\;.$ (1)

The transition form factor (TFF) of the pion is a single function, which
appears both in the doubly-virtual and in the singly-virtual case:
$F_{\pi\gamma^{*}\gamma^{*}}(q^{2},0)=F_{\pi\gamma\gamma^{*}}(q^{2})$. The
function $\tilde{G}$ collects all additional contributions not containing any
poles at $s=M_{\pi}^{2}$.

For the muon $g-2$ calculation we need to take the limit $q_{4}\to 0$, which
changes the kinematics as follows: $s=q_{3}^{2}$, $t=q_{2}^{2}$, and
$u=q_{1}^{2}$, leading to

$\displaystyle\bar{\Pi}_{1}(q_{1}^{2},q_{2}^{2},q_{3}^{2})$
$\displaystyle\equiv\hat{\Pi}_{1}(q_{1}^{2},q_{2}^{2},q_{3}^{2},0;q_{3}^{2},q_{2}^{2},q_{1}^{2})$
$\displaystyle=\frac{F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})F_{\pi\gamma\gamma^{*}}(q_{3}^{2})}{q_{3}^{2}-M_{\pi}^{2}}$
$\displaystyle+G(q_{1}^{2},q_{2}^{2},q_{3}^{2})\;,$ (2)

where
$G(q_{1}^{2},q_{2}^{2},q_{3}^{2})=\tilde{G}(q_{1}^{2},q_{2}^{2},q_{3}^{2},0;q_{3}^{2},q_{2}^{2},q_{1}^{2})$.
We stress that taking the limit $q_{4}\to 0$ starting from the representation
in (2) unambiguously leads to (2). The splitting between the first and second
term is inherited from the splitting between pole term and the rest for
general kinematics, but is nonetheless unique. If one wanted to identify the
pole term directly for $g-2$ kinematics, the first term in (2) should only
have its residue as numerator, thereby separating any additional $q_{3}^{2}$
dependence carried by the TFFs: in other words, separating the pseudoscalar
poles from the vector-meson ones in the TFFs (which correspond to cuts from
$2\pi$, $3\pi$, etc. intermediate states, if one does not take the narrow-
width approximation). Both definitions of the pion pole (for general or $g-2$
kinematics) are possible and the relation between the two is completely
understood. The connection of this aspect with the SDCs has been discussed in
detail in Knecht:2020xyr . Here we adopt the splitting between pion pole and
the function $G$ given in (2) and concentrate our discussion on the latter.

In the MV limit, $\hat{q}^{2}\equiv q_{1}^{2}=q_{2}^{2}\gg q_{3}^{2}$,
$\hat{q}^{2}\gg\Lambda_{\text{QCD}}^{2}$, with no constraints on $q_{3}^{2}$,
the function $\bar{\Pi}_{1}$ reads:

$\displaystyle\bar{\Pi}_{1}(\hat{q}^{2},\hat{q}^{2},q_{3}^{2})$
$\displaystyle=\frac{F_{\pi\gamma^{*}\gamma^{*}}(\hat{q}^{2},\hat{q}^{2})F_{\pi\gamma\gamma^{*}}(q_{3}^{2})}{q_{3}^{2}-M_{\pi}^{2}}$
$\displaystyle+G(\hat{q}^{2},\hat{q}^{2},q_{3}^{2})\;,$ (3)

which can be further simplified taking into account the leading-order OPE for
the pion TFF Lepage:1979zb ; Lepage:1980fj :

$F_{\pi\gamma^{*}\gamma^{*}}(\hat{q}^{2},\hat{q}^{2})=-\frac{2F_{\pi}}{3\hat{q}^{2}}+\mathcal{O}(\hat{q}^{-3})\;.$
(4)

We now carry out the separation between the pion pole in $g-2$ kinematics (in
the chiral limit) from the rest, and rewrite the expression for
$\bar{\Pi}_{1}$ as follows:

$\displaystyle\bar{\Pi}_{1}(\hat{q}^{2},\hat{q}^{2},q_{3}^{2})$
$\displaystyle=-\frac{2F_{\pi}}{3\hat{q}^{2}}\bigg{[}\frac{F_{\pi\gamma\gamma}}{q_{3}^{2}}$
$\displaystyle+\frac{F_{\pi\gamma\gamma^{*}}(q_{3}^{2})-F_{\pi\gamma\gamma}}{q_{3}^{2}}+\mathcal{O}(M_{\pi}^{2})\bigg{]}$
$\displaystyle+G(\hat{q}^{2},\hat{q}^{2},q_{3}^{2})+\mathcal{O}(\hat{q}^{-3})\;,$
(5)

where

$F_{\pi\gamma\gamma}:=F_{\pi\gamma\gamma^{*}}(0)=\frac{1}{4\pi^{2}F_{\pi}}\;.$
(6)

Since we know how the amplitude has to behave in the chiral limit
Melnikov:2003xd :

$\displaystyle\bar{\Pi}_{1}(\hat{q}^{2},\hat{q}^{2},q_{3}^{2})\Big{|}_{m_{q}=0}=-\frac{1}{6\pi^{2}}\frac{1}{\hat{q}^{2}q_{3}^{2}}+\mathcal{O}(\hat{q}^{-3})\;,$
(7)

we have to conclude that Colangelo:2019uex

$\displaystyle G(\hat{q}^{2},\hat{q}^{2},q_{3}^{2})\Big{|}_{m_{q}=0}$
$\displaystyle=\frac{2F_{\pi}}{3\hat{q}^{2}}\frac{F_{\pi\gamma\gamma^{*}}(q_{3}^{2})-F_{\pi\gamma\gamma}}{q_{3}^{2}}\bigg{|}_{m_{q}=0}$
$\displaystyle+\mathcal{O}(\hat{q}^{-3})\;.$ (8)

This remarkable result is actually a consequence of the non-renormalization of
the axial anomaly Adler:1969gk ; Bell:1969ts ; Bardeen:1969md ; Wess:1971yu ;
Witten:1983tw , as first discussed in Melnikov:2003xd  (see Vainshtein:2002nv
; Knecht:2003xy  for a full account of non-renormalization theorems for the
$VVA$ correlator):

$\displaystyle
w_{L}(q_{1}^{2},q_{2}^{2},(q_{1}+q_{2})^{2})=\frac{2N_{C}}{(q_{1}+q_{2})^{2}}+\mathcal{O}(M_{\pi}^{2})\;.$
(9)

The expression on the right-hand side looks like the pion-pole contribution in
the chiral limit, though without the (properly normalized) TFF in the
numerator. As discussed in detail in Masjuan:2020jsf ; Knecht:2020xyr  this
implies a constraint between the contribution of the pion and that of
transverse degrees of freedom, such that the only effect of their contribution
in the chiral limit is to replace the pion TFF by its value at
$q_{3}^{2}=q_{4}^{2}=0$. A crucial point is that the constraint (8) applies
for arbitrary values of $q_{3}^{2}$, but only for large values of
$\hat{q}^{2}\gg\Lambda_{\text{QCD}}^{2}$ (and $\hat{q}^{2}\gg q_{3}^{2}$),
where the OPE in the MV limit is valid: for non-asymptotic values of
$q_{1,2}^{2}$ the connection between the four- and the three-point function
gets lost and nothing can be inferred on the behavior of
$G(q_{1}^{2},q_{2}^{2},q_{3}^{2})$. The key point in assessing the relevance
of the non-renormalization theorem for $\langle VVA\rangle$ for the numerical
evaluation of $a_{\mu}$ thus concerns the weight of the integration region in
which the constraint applies, as we will discuss in detail in Sect. 5.

##  3 On axial-vector contributions to the longitudinal function in the
dispersive approach

Since the non-renormalization theorems on the $\langle VVA\rangle$ function
interrelate transverse and longitudinal degrees of freedom, it is clear that
axial-vector states play a role in fulfilling (8). Within the dispersive
framework for HLbL, the inclusion of axial-vector mesons suffers from two
closely related difficulties: on the one hand, the contribution of narrow
states depends on the choice of basis—these ambiguities apply to all narrow
resonances beyond the pseudoscalar ones and have been recently discussed for
scalar contributions Danilkin:2021icn . As only the full HLbL scattering
amplitude needs to be basis independent, a phenomenological evaluation of
axial-vector effects thus must proceed in accordance with the other
contributions. The axial-vector exchanges discussed in the context of SDCs,
both in hQCD Leutgeb:2019gbz ; Cappiello:2019hwh  and in other implementations
Melnikov:2003xd ; Melnikov:2019xkq ; Roig:2019reh ; Masjuan:2020jsf ;
Pauk:2014rta ; jegerlehner:2017gek , typically refer to a Lagrangian model,
which can differ by non-pole pieces from a dispersive definition, depending on
the choice of basis.222Note that, in addition, all such estimates assume the
validity of a narrow-width approximation. Since the main branching fractions
proceed into three- or higher-multiplicity final states, a full dispersive
treatment of axial-vector intermediate states is difficult, but for $S$\- and
$D$-wave resonances that decay predominantly into two-meson states a
comparison to implementations in terms of $\gamma^{*}\gamma^{*}$ amplitudes
GarciaMartin:2010cw ; Hoferichter:2011wk ; Moussallam:2013una ;
Danilkin:2018qfn ; Hoferichter:2019nlq ; Danilkin:2019opj ; Danilkin:2017lyn ;
Lu:2020qeo  is possible, see Danilkin:2021icn  for the scalar case.

The second difficulty in the dispersive approach concerns kinematic
singularities: while the basis of Colangelo:2017fiz  is free from kinematic
singularities in the dispersed Mandelstam variable, it still contains
singularities in the photon virtualities, with residues that vanish for the
entire HLbL contribution due to the presence of sum rules. As narrow
resonances do not fulfill the sum rules individually, a further ambiguity in
their contribution is introduced by the subtraction scheme of the singular
parts, which again disappears only for the entire HLbL contribution. This
second complication does not affect pseudoscalar or scalar contributions, but
appears for axial and higher-spin resonances in the basis of Colangelo:2017fiz
. By employing the sum rules, we have now constructed a new basis that
explicitly removes all kinematic singularities from axial-vector
contributions, while leaving pseudoscalar and scalar contributions unaltered,
thereby solving this second issue in the case of axial-vector contributions.
In this basis the contribution of a single axial-vector meson (with mass
$M_{A}$) to the function $G$ takes the form:

$\displaystyle
G(q_{1}^{2},q_{2}^{2},q_{3}^{2})=\frac{G_{2}(q_{1}^{2},q_{2}^{2})G_{1}(q_{3}^{2})}{M_{A}^{6}}\;,$
(10)

where

$\displaystyle G_{2}(q_{1}^{2},q_{2}^{2})$
$\displaystyle=(q_{1}^{2}-q_{2}^{2})\mathcal{F}_{1}(q_{1}^{2},q_{2}^{2})$
$\displaystyle+q_{1}^{2}\mathcal{F}_{2}(q_{1}^{2},q_{2}^{2})+q_{2}^{2}\mathcal{F}_{2}(q_{2}^{2},q_{1}^{2})\;,$
$\displaystyle G_{1}(q^{2})$
$\displaystyle=\mathcal{F}_{1}(q^{2},0)+\mathcal{F}_{2}(q^{2},0)=\frac{G_{2}(q^{2},0)}{q^{2}}\;,$
(11)

and $\mathcal{F}_{1,2}(q_{1}^{2},q_{2}^{2})$ are two of the three TFFs of an
axial-vector meson, see Hoferichter:2020lap  for the precise definitions. The
third one, $\mathcal{F}_{3}$, does not appear in the expression above but is
related to $\mathcal{F}_{2}$ by the symmetry properties of the TFFs:

$\displaystyle\mathcal{F}_{1}(q_{1}^{2},q_{2}^{2})=-\mathcal{F}_{1}(q_{2}^{2},q_{1}^{2})\;,$
$\displaystyle\mathcal{F}_{2}(q_{2}^{2},q_{1}^{2})=-\mathcal{F}_{3}(q_{1}^{2},q_{2}^{2})\;.$
(12)

The expression (3) shows that the dispersive contribution of axial-vector
mesons to the function $G$ has the form of non-pole terms, but does not
vanish. Our new basis avoids any kinematic singularities in (10) and makes the
dependence on the virtualities unambiguous for basis changes that preserve
this property, up to terms that are subleading for $q_{i}^{2}\gg M_{A}^{2}$.

As the remaining ambiguities become irrelevant for asymptotic virtualities,
(10) leads to an interesting model-independent conclusion. The light-cone
expansion determines the asymptotic behavior of
$\mathcal{F}_{1}=\mathcal{O}(1/q_{i}^{6})$,
$\mathcal{F}_{2}=\mathcal{O}(1/q_{i}^{4})$, with coefficients determined via
decay constants in analogy to (4), see Hoferichter:2020lap . This implies
that, asymptotically,

$G_{2}(\hat{q}^{2},\hat{q}^{2})=\mathcal{O}(1/\hat{q}^{2})\;,\quad
G_{1}(q_{3}^{2})=\mathcal{O}(1/q_{3}^{4})\;.$ (13)

Figure 1: Singly-virtual TFF of the ground-state axial-vector meson
$f_{1}(1285)$: comparison of the dipole ansatz used to fit the L3 data (black
curve and gray band) Achard:2001uu , the HW2 hQCD model representations (blue
dashed and turquoise dot-dashed curves) Leutgeb:2019gbz , and the one obtained
from (14) using as input the $\pi^{0}$, $\eta$, and $\eta^{\prime}$ TFFs
(green, yellow, and red curves) from Colangelo:2019uex .

Since (10) factorizes into parts dependent on $q_{1,2}^{2}$ and $q_{3}^{2}$,
respectively, we find that (8) decomposes into two equations that need to be
fulfilled separately:

$\displaystyle\lim_{\hat{q}^{2}\to\infty}x\frac{G_{2}(\hat{q}^{2},\hat{q}^{2})}{M_{A}^{4}}$
$\displaystyle=-\frac{2}{3\hat{q}^{2}}+\mathcal{O}(\hat{q}^{-3})\;,$ (14)
$\displaystyle\frac{G_{1}(q_{3}^{2})}{xM_{A}^{2}}$
$\displaystyle=-F_{\pi}\frac{F_{\pi\gamma\gamma^{*}}(q_{3}^{2})-F_{\pi\gamma\gamma}}{q_{3}^{2}}\;,$

with $x$ an unknown, but constant factor that depends on the axial-vector
TFFs. While the asymptotic form of $G_{2}$ thus matches, the axial-vector
contribution to $G_{1}(q_{3}^{2})$ decreases too fast, mirroring the need for
an infinite tower of axial-vector states in the hQCD models. This mismatch can
also be illustrated by comparing $Q^{2}G_{1}(-Q^{2})/G_{1}(0)$ evaluated from
the axial-vector TFFs with (14), see Fig. 1, where we concentrated on the
contribution from $\mathcal{F}_{2}$, given that $\mathcal{F}_{1}$ is
suppressed for several reasons: for small virtualities its antisymmetry
implies $\mathcal{F}_{1}(-Q^{2},0)\sim Q^{2}$, for large virtualities due to
the asymptotic behavior, and phenomenologically due to small couplings
Zanke:2021wiq . The comparison curves for $\eta$ and $\eta^{\prime}$ show that
this qualitative behavior does not depend on the isospin channel, reinforcing
that a single state is not sufficient to implement the SDCs.

##  4 Three approaches to satisfy short-distance constraints

In this section we compare the different solutions to the SDCs that have been
proposed so far in the literature Melnikov:2003xd ; Colangelo:2019lpu ;
Colangelo:2019uex ; Leutgeb:2019gbz ; Cappiello:2019hwh ; Melnikov:2019xkq ,
in terms of the different representations of the functions $w_{L}(q^{2})\equiv
w_{L}(q^{2},0,q^{2})$ and $G(q_{1}^{2},q_{2}^{2},q_{3}^{2})$, as constrained
by the non-renormalization theorem (9) and the asymptotic behavior (8) in the
chiral limit. We consider the models as they are, in other words, we take each
one of them as an approximation to the total contribution to $w_{L}$ and the
longitudinal SDC. The question of building a better model, possibly by
combining features or degrees of freedom of the present ones will be touched
upon in Sect. 7. The analysis based on interpolants Ludtke:2020moa  will be
included in the numerical comparison in the following section.

###  4.1 The Melnikov–Vainshtein model

After deriving the SDC for the HLbL tensor, MV go beyond the asymptotic limit
and propose a model that by construction satisfies (7):

$\displaystyle\bar{\Pi}_{1}^{\mathrm{MV}}(q_{1}^{2},q_{2}^{2},q_{3}^{2})$
$\displaystyle=\frac{F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})F_{\pi\gamma\gamma}}{q_{3}^{2}-M_{\pi}^{2}}\;,$
(15)

with the shift in the pole position from zero to $M_{\pi}^{2}$ as the only
effect of the light quark masses considered. This implies that, even though no
additional contributions beyond the pion pole are introduced explicitly, such
additional contributions are implicitly assumed to be completely determined
not only in the asymptotic region, as implied by (8), but everywhere:

$\displaystyle G^{\text{MV}}(q_{1}^{2},q_{2}^{2},q_{3}^{2})=$ $\displaystyle-
F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})$
$\displaystyle\frac{F_{\pi\gamma\gamma^{*}}(q_{3}^{2})-F_{\pi\gamma\gamma}}{q_{3}^{2}}+\mathcal{O}(M_{\pi}^{2})\;.$
(16)

Equation (16) is a very strong assumption, with no apparent physical
justification: it extrapolates a constraint only valid at asymptotically high
energies (8) to all possible values of $q_{1,2}^{2}$, all the way down to
$q_{1}^{2}=q_{2}^{2}=0$. As such, it has the potential to significantly affect
the value of the HLbL contribution to $a_{\mu}$—a quantity most sensitive to
low $q_{i}^{2}$. That this indeed happens has already been shown explicitly in
Colangelo:2019uex , and will be discussed in more detail below.

Regarding the three-point function, the MV model reads

$\displaystyle w_{L}^{\text{MV}}(q^{2})=\frac{2N_{C}}{q^{2}-M_{\pi}^{2}}\;,$
(17)

which again amounts to including as only chiral correction the one that shifts
the pole in the pion propagator. Since chiral corrections become negligible at
large $q^{2}$, (17) is a well-motivated model at all $q^{2}$, with small
deviations from the truth expected only at intermediate energies (chiral
corrections may become more sizable for the $\eta$ and $\eta^{\prime}$
channels).

###  4.2 The Leutgeb–Rebhan and Cappiello–Catà–D’Ambrosio–Greynat–Iyer models

In two recent papers, Leutgeb–Rebhan Leutgeb:2019gbz  (LR) and Cappiello et
al. Cappiello:2019hwh  (CCDGI) have proposed models based on hQCD to satisfy
the SDCs. As the discussion in Sect. 3 shows, a solution in terms of a single
axial-vector meson (per isospin channel) is essentially excluded, and indeed
in these models the solution emerges from a resummation of an infinite tower
of axial-vector mesons. For simplicity, we concentrate here on the model
presented by CCDGI in Cappiello:2019hwh , which is equivalent to the HW2 model
in Leutgeb:2019gbz , although the two groups make different choices for the
parameters.333We will not consider the HW1 model discussed in Leutgeb:2019gbz
simply because it does not offer a simple analytic representation like HW2.
The representation of this model for the function $G$ reads (in the notation
of CCDGI)

$\displaystyle G^{\text{HW2}}(q_{1}^{2},q_{2}^{2},q_{3}^{2})$
$\displaystyle=-\frac{F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})F_{\pi\gamma\gamma^{*}}(q_{3}^{2})}{q_{3}^{2}}$
(18)
$\displaystyle-\frac{F_{\pi\gamma\gamma}^{2}}{q_{3}^{2}}\int_{0}^{z_{0}}dz\,\alpha^{\prime}(z)v_{1}(z)v_{2}(z)v_{3}(z)\;,$

where $\alpha(z)=1-z^{2}/z_{0}^{2}$, with $z_{0}=(\\!\sqrt{2}\pi
F_{\pi})^{-1}$ and

$v_{i}(z)=zQ_{i}\left[K_{1}(zQ_{i})+\frac{K_{0}(z_{0}Q_{i})}{I_{0}(z_{0}Q_{i})}I_{1}(zQ_{i})\right]\;,$
(19)

with $K_{n}(x)$ and $I_{n}(x)$ modified Bessel functions Cappiello:2019hwh .
The same function $v_{i}(z)$ also determines the pion TFF

$F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})=-F_{\pi\gamma\gamma}\int_{0}^{z_{0}}dz\,\alpha^{\prime}(z)v_{1}(z)v_{2}(z)\;.$
(20)

This representation correctly reproduces the high-$q^{2}$ limit of the TFF
shown in (4), the Brodsky–Lepage limit of the singly-virtual pion TFF
Lepage:1979zb ; Lepage:1980fj :

$F_{\pi\gamma\gamma^{*}}(\hat{q}^{2})=-\frac{2F_{\pi}}{\hat{q}^{2}}+\mathcal{O}(\hat{q}^{-3}),$
(21)

as well as, by construction, the normalization at $q_{1}^{2}=q_{2}^{2}=0$. A
convenient rewriting for $G$ is

$\displaystyle G^{\text{HW2}}(q_{1}^{2},q_{2}^{2},q_{3}^{2})$
$\displaystyle=-F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})\frac{F_{\pi\gamma\gamma^{*}}(q_{3}^{2})-F_{\pi\gamma\gamma}}{q_{3}^{2}}$
$\displaystyle-\frac{F_{\pi\gamma\gamma}^{2}}{q_{3}^{2}}\int_{0}^{z_{0}}dz\,\alpha^{\prime}(z)v_{1}(z)v_{2}(z)\bar{v}_{3}(z),$
(22)

where $\bar{v}_{3}(z)=v_{3}(z)-1$, as it shows that there is no divergence at
$q_{3}^{2}=0$ (the integral vanishes for $q_{3}^{2}\to 0$). Note that the
first term in (22) coincides, up to chiral corrections, with the MV model,
which can therefore be viewed as a truncated version of the hQCD model. In the
HW2 model, however, the first and the second non-factorizable term always come
together as they have the same physical origin: both arise from the
resummation of the whole tower of axial-vector mesons. In the numerical
analysis below we will see that, while the first term is dominant for
asymptotic values of $q_{1}^{2}\sim q_{2}^{2}$, for low momenta, they are
equally important and in fact cancel each other.

We also observe that the HW2 model offers a compact and convenient
representation for the function $\bar{\Pi}_{1}$:

$\displaystyle\bar{\Pi}_{1}^{\text{HW2}}$
$\displaystyle=F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})\frac{F_{\pi\gamma\gamma}}{q_{3}^{2}-M_{\pi}^{2}}$
$\displaystyle\times\left[1+\frac{M_{\pi}^{2}(F_{\pi\gamma\gamma^{*}}(q_{3}^{2})-F_{\pi\gamma\gamma})}{q_{3}^{2}F_{\pi\gamma\gamma}}\right]$
$\displaystyle-\frac{F_{\pi\gamma\gamma}^{2}}{q_{3}^{2}}\int_{0}^{z_{0}}dz\,\alpha^{\prime}(z)v_{1}(z)v_{2}(z)\bar{v}_{3}(z),$
(23)

where one can clearly see that the corrections to the $1/q_{3}^{2}$ behavior,
i.e., the pion pole in $g-2$ kinematics, vanish in the chiral limit—the
integral behaves as $\mathcal{O}(\hat{q}^{-4})$.

Finally, for the $VVA$ correlation function the HW2 model gives

$\displaystyle
w_{L}^{\text{HW2}}(q^{2})=\frac{2N_{C}}{q^{2}-M_{\pi}^{2}}\left[1+\frac{M_{\pi}^{2}(F_{\pi\gamma\gamma^{*}}(q^{2})-F_{\pi\gamma\gamma})}{q^{2}F_{\pi\gamma\gamma}}\right]\;,$
(24)

which shows that the first term again corresponds to the MV model. In this
case it is evident that the corrections to the MV model are of
$\mathcal{O}(M_{\pi}^{2})$ for any value of $q^{2}$, and therefore expected to
be small everywhere. This expression also shows that the first term in the
representation of the four-point function (23) takes the form
$F_{\pi\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})F_{\pi\gamma\gamma}w_{L}(q_{3}^{2})/(2N_{C})$.
While there is no harm in approximating the $w_{L}$ function with (17), it is
dropping the non-factorizable second term in (23) that amounts to an
uncontrolled approximation. Its numerical impact will be shown in Sect. 5.

###  4.3 Regge model of excited pseudoscalars

In the model we presented in Colangelo:2019lpu ; Colangelo:2019uex , we
considered only the contribution of excited pseudoscalars to the function $G$:

$\displaystyle G^{\text{eP}}(q_{1}^{2},q_{2}^{2},q_{3}^{2})$
$\displaystyle=\sum_{i=1}^{\infty}\frac{F_{P_{i}\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})F_{P_{i}\gamma\gamma^{*}}(q_{3}^{2})}{q_{3}^{2}-M_{P_{i}}^{2}}\;.$
(25)

Clearly, by dropping axial-vector intermediate states, which contribute to
this function according to (10), we are transferring their unique role in the
chiral limit to the pseudoscalars, which amounts to effectively changing their
chiral behavior, in particular the coupling of the excited pseudoscalars to
the axial-vector current, which has to vanish in the chiral limit. This
procedure cannot be strictly justified, but is similar in spirit to models
that use constituent quark masses. In order to remove some of the model
dependence, after matching to the behavior dictated by the OPE we are
replacing our model in the asymptotic region with the perturbative QCD quark
loop.

We have imposed as constraint to our model that it satisfy (8) only for
$q_{3}^{2}\gg\Lambda_{\text{QCD}}$, which is a less ambitious goal than the
one reached by both models described above:

$\displaystyle\lim_{\hat{q}^{2}\to\infty}\hat{q}^{2}G^{\text{eP}}(\hat{q}^{2},\hat{q}^{2},q_{3}^{2})$
$\displaystyle=-\frac{1}{6\pi^{2}q_{3}^{2}}+\mathcal{O}(q_{3}^{-3})\;.$ (26)

By construction, our model takes into account singularities that are known to
be present in the spectrum of QCD (the low-lying pseudoscalar excitations).
The resummation of all higher excitations is used essentially to achieve the
matching to the asymptotic behavior, but its precise form is inessential.

The resulting representation for the longitudinal component of the $VVA$
correlator becomes444Note that here we are using a peculiarity of our model
that
$F_{P_{i}\gamma^{*}\gamma^{*}}(\hat{q}^{2},\hat{q}^{2})=-\frac{1}{6\pi^{2}F_{\pi\gamma\gamma}\hat{q}^{2}}+\mathcal{O}(\hat{q}^{-3})\;$
for all $i$.

$\displaystyle
w_{L}^{\text{eP}}(q^{2})=\frac{2N_{C}}{F_{\pi\gamma\gamma}}\left[\frac{F_{\pi\gamma\gamma^{*}}(q^{2})}{q^{2}-M_{\pi}^{2}}+\sum_{i=1}^{\infty}\frac{F_{P_{i}\gamma\gamma^{*}}(q^{2})}{q^{2}-M_{P_{i}}^{2}}\right]\;.$
(27)

We stress that this model was not conceived to approximate this function other
than for asymptotic values of its argument, and its use in Colangelo:2019lpu ;
Colangelo:2019uex  was limited to the four-point function. However, we find it
instructive to provide this expression and compare it numerically to the other
models.

##  5 Numerical comparison of the three models

In this section we compare the three models numerically, first for the
$w_{L}(q^{2})$ function, then the $G$ function and its contribution to
$a_{\mu}$.

Figure 2:  Comparison of the MV, the HW2, and our model (solid) for the
isovector component of $-Q^{2}w_{L}(-Q^{2})$ as a function of $Q$. Our model
is also broken down into pion and sum of excited pseudoscalar contributions.

###  5.1 The $w_{L}$ function

A numerical comparison between the hQCD model and the MV model for the $w_{L}$
function was provided and discussed in Cappiello:2019hwh  and clearly showed
that the difference between the two models amounts to chiral corrections
neglected by MV, which the hQCD model estimates to be numerically very small.
The picture that emerges from the comparison is that the $w_{L}$ function is
essentially determined by its low-energy (fixed by the pion pole) and its
high-energy behavior, with no room for any structure in between. In Fig. 2 we
repeat the comparison for the isovector channel and show in addition the
contribution of the pion if one includes its transition form factor in the
numerator—in other words, according to the dispersive definition of the pion
contribution for general kinematics. The difference between the $\pi^{0}$ and
the CCDGI/HW2 curves is the contribution of the axials, but its main effect is
to remove the TFF from the numerator, as the minute difference to the MV model
shows. The hQCD models confirm that the MV model appears to be an excellent
approximation to the true $w_{L}$ function in QCD. It is instructive to see
algebraically how the contribution of the axial-vector states manages to
remove the TFF from the pion-pole contribution and also to be able to estimate
the additional corrections to it, but for all practical purposes, and unless
the highest precision is required, the MV model seems to provide an excellent
description of $w_{L}$. In the same plot we also show our model of excited
pseudoscalars. This is designed to agree with the other two for asymptotic
values of $q^{2}$, as indeed it does. At low energy, where the pion
contribution dominates, it also agrees with the other two, but the transition
region is not as smooth and shows some structure, with up to 30% discrepancy
with the other two models. There is nothing to be read into this discrepancy
other than the fact that the model was never designed to provide a good
description of $w_{L}$: the structure it shows in the intermediate region just
reflects the fact that it was not required to fulfill any constraints here. If
required, the model could be refined to improve the transition between low-
and high-energy constraints. Whether this is of relevance in the context of
the four-point function will be discussed in the following.

Figure 3: Comparison between the non-factorizable (NF) and the MV term (second
and first term in (22), respectively, for $q_{1}^{2}=q_{2}^{2}=-Q^{2}$ and
$q_{3}^{2}=-Q_{3}^{2}$) in the CCDGI model (set 2).

###  5.2 The $G$ function and the role of $w_{L}$

Such a good understanding of the function $w_{L}$ and the simple and accurate
description provided by the MV model raises the question which role the three-
point function plays in determining the four-point function: how strongly does
the accurate knowledge of $w_{L}$ constrain $\bar{\Pi}_{1}$ or $G$? MV have
shown that for asymptotic values of $q_{1}^{2}\sim q_{2}^{2}$, the leading
behavior of $G$ is completely fixed by $w_{L}$, but how relevant is the
asymptotic region in determining the contribution of $G$ to $a_{\mu}$? In
Melnikov:2019xkq  it is argued that the kinematic region $q_{1}^{2}\sim
q_{2}^{2}\gg q_{3}^{2}$ “provides the largest contribution to
$a_{\mu}^{\text{HLbL}}$,” but we are not aware of any quantitative basis for
such a statement. When we proposed an alternative way to fulfill the SDCs and
compared to MV Colangelo:2019lpu ; Colangelo:2019uex , we showed that the
large difference between ours and the MV model arose precisely in the
low-$q_{i}^{2}$ region and, moreover, that the largest contribution to
$a_{\mu}^{\text{HLbL}}$ in the MV model itself came from the same region.

Figure 4:  Comparison of the MV (dot-dashed), the HW2 (dashed), and our model
(solid) for $-Q_{3}^{2}G(-Q^{2},-Q^{2},-Q_{3}^{2})/F_{\pi\gamma\gamma}^{2}$ as
a function of $Q$ for different values of $Q_{3}$.

The hQCD models, which satisfy exactly the anomaly and the MV constraints,
allow us to test the approximation made in the MV model in a more quantitative
way. As discussed above, there are two approximations made by MV in their
representation of $\bar{\Pi}_{1}$: the first is to neglect chiral correction
in $w_{L}(q_{3}^{2})$, which is a very good one as we have just seen, but the
second one is to neglect non-factorizable corrections, which in the hQCD
models are given by the integral term in (23). To establish the relevance of
the function $w_{L}$ in the calculation of the four-point function and its
contribution to $a_{\mu}$ we can therefore compare the non-factorizable and
the MV term: the region where the latter dominates is the region where the MV
limit matters, i.e., where the $w_{L}$ function plays an important role, and
the MV model is a good approximation. This comparison is shown in Fig. 3: even
for the modest requirement that the non-factorizable term amounts to at most
30% of the MV term, the minimum value of $Q$ for which this is satisfied is
above $2\,\text{GeV}$. At the matching point $Q_{\text{match}}=1.7$ GeV,
adopted in Colangelo:2019lpu ; Colangelo:2019uex  for the transition between
the hadronic and the pQCD description, and which will again be used below, the
non-factorizable term is at least a 50% correction to the MV term.

| MV | CCDGI | LR | LP | PS Regge
---|---|---|---|---|---
| set 1 | set 2 | HW2 |  HW2${}_{\text{UV-fit}}$
| $\Delta a_{\mu}^{\pi/a_{1}}\times 10^{11}$
$Q^{2}_{i}>Q^{2}_{\text{match}}\quad\forall i$ | $1.4$ | $0.5$ | $0.8$ | $0.6$ | $0.8$ | $0.9$ | $0.7$
$Q^{2}_{1,2}>Q^{2}_{\text{match}}>Q^{2}_{3}$ | $1.4$ | $0.8$ | $1.0$ | $0.8$ | $1.0$ | $0.3$ | $0.4$
$Q^{2}_{i,3}>Q^{2}_{\text{match}}>Q^{2}_{j}\quad i\neq j\neq 3$ | $0.8$ | $0.2$ | $0.3$ | $0.2$ | $0.3$ | $0.4$ | $0.4$
$Q^{2}_{i}>Q^{2}_{\text{match}}>Q^{2}_{j,k}\quad i\neq j\neq k$ | $0.8$ | $0.3$ | $0.4$ | $0.3$ | $0.4$ | $0.3$ | $0.2$
$Q^{2}_{\text{match}}>Q^{2}_{i}\quad\forall i$ | $11.8$ | $2.2$ | $1.7$ | $2.3$ | $1.8$ | $0.7$ | $1.0$
Total | $16.2$ | $4.0$ | $4.2$ | $4.2$ | $4.3$ | $2.6$ | $2.7$
| $\Delta a_{\mu}^{\eta/f_{1}+\eta^{\prime}/f^{\prime}_{1}}\times 10^{11}$
$Q^{2}_{i}>Q^{2}_{\text{match}}\quad\forall i$ | $3.4$ | $1.4$ | $1.7$ | $1.7$ | $2.5$ | $2.5$ | $3.1$
$Q^{2}_{1,2}>Q^{2}_{\text{match}}>Q^{2}_{3}$ | $2.1$ | $2.1$ | $2.3$ | $2.5$ | $3.0$ | $0.6$ | $1.1$
$Q^{2}_{i,3}>Q^{2}_{\text{match}}>Q^{2}_{j}\quad i\neq j\neq 3$ | $1.9$ | $0.6$ | $0.7$ | $0.6$ | $0.9$ | $1.2$ | $1.6$
$Q^{2}_{i}>Q^{2}_{\text{match}}>Q^{2}_{j,k}\quad i\neq j\neq k$ | $1.7$ | $0.8$ | $0.9$ | $0.9$ | $1.1$ | $0.7$ | $0.9$
$Q^{2}_{\text{match}}>Q^{2}_{i}\quad\forall i$ | $12.9$ | $5.6$ | $5.1$ | $6.8$ | $5.4$ | $1.5$ | $3.1$
Total | $22.1$ | $10.4$ | $10.7$ | $12.5$ | $12.8$ | $6.5$ | $9.9$
Grand total ($\pi/a_{1}+\eta/f_{1}+\eta^{\prime}/f^{\prime}_{1}$) | $38.3$ | $14.3$ | $14.9$ | $16.7$ | $17.1$ | $9.1$ | $12.6$
Table 1:  Contribution of $G$ to $a_{\mu}$ (referred to as the longitudinal SD
contribution in Colangelo:2019lpu ; Colangelo:2019uex  and the longitudinal
axial-vector contribution in Leutgeb:2019gbz ; Cappiello:2019hwh ) from the
isovector and isoscalar plus isosinglet channels broken down in different
integration regions ($Q_{\text{match}}=1.7\,\text{GeV}$). The notation for the
mixed regions includes the respective crossed versions, e.g., the second line
gives the contribution from $\bar{\Pi}_{1}$ in the region
$Q_{1,2}^{2}>Q_{\text{match}}^{2}>Q_{3}^{2}$ and from $\bar{\Pi}_{2}$ in the
region $Q_{1,3}^{2}>Q_{\text{match}}^{2}>Q_{2}^{2}$, in such a way that the
region in which the SDC1 applies is contained in this (and partly the first)
row, while the third row has a scaling $1/Q^{4}$ in the hard momenta. Due to
different mixing patterns the $\eta/f_{1}$ and $\eta^{\prime}/f_{1}^{\prime}$
contributions cannot be compared separately. Note that the Regge-model
contribution to the asymptotic region is not yet replaced by the OPE result.
The numbers for LP refer to the “reference interpolant” of Ludtke:2020moa .
The HW1 model, which we have not considered here, gives a higher contribution
$\Delta a_{\mu}=23.2\times 10^{-11}$ Leutgeb:2019gbz . All entries are
understood to be accurate at the level of $\pm 0.1$ due to the applied
numerical integration methods, other (model-dependent) errors are not shown.

This suggests that if we compare the MV and the hQCD models at the level of
the function $G$, the agreement is going to be much worse than for the
function $w_{L}$. To verify this expectation, we plot the isovector component
of the function $-Q_{3}^{2}G(-Q^{2},-Q^{2},-Q_{3}^{2})$ as a function of
$Q^{2}$, see Fig. 4. The dashed curves show the HW2 model, whereas the dot-
dashed ones only show the first term of (22), which corresponds to the MV
model. The plot shows very clearly that the latter two versions tend to the
same asymptotic limit, even for low values of $Q_{3}$, but that they differ
significantly at low $Q^{2}$. While at large $Q^{2}$ the second term in (22)
is subdominant and can be neglected, it compensates exactly the first one for
low $Q^{2}$, so that their sum vanishes. This is expected for axial-vector
mesons, and is already seen in Fig. 4. If one keeps just the first term, i.e.,
the MV model, this grows at low $Q^{2}$, reaching a finite limit for $Q^{2}=0$
(not visible in the plot range, because it is quite large). Fig. 4 also
displays our model (solid curves), which has non-vanishing limits for
$Q^{2}=0$, too. This is well understood, however, because excited
pseudoscalars do couple to two real photons (even in the chiral limit). As
detailed in Colangelo:2019uex , these couplings are compatible with the
available phenomenological information, which still suffers from large
uncertainties. By construction, our model also agrees with the other two for
large $Q^{2}$ and $Q_{3}^{2}$, whereas for low $Q_{3}^{2}$ it does not agree
well with the other two even as $Q^{2}$ grows.

###  5.3 Contribution to $a_{\mu}$ of the function $G$

We can now address the question of how these differences are reflected in the
calculation of the contribution to $a_{\mu}$. We do so by breaking down the
contributions from different kinematic regions and separating the isovector
channel from the isoscalar and isosinglet ones. Identifying the isoscalar and
isosinglet pieces with the physical states ignores mixing effects, which
implies that the $\eta/f_{1}$ and $\eta^{\prime}/f_{1}^{\prime}$ cannot be
compared separately. In general, the correct implementation of mixing effects
requires two mixing angles (see Gan:2020aco  for a review), but the
differences can be illustrated based on the simple U(3) formula

$\frac{\Gamma(P\to\gamma\gamma)}{\Gamma(P^{\prime}\to\gamma\gamma)}=\frac{M_{P}}{M_{P^{\prime}}}\cot^{2}\big{(}\theta_{A}-\theta_{0}),\,\theta_{0}=\arcsin\frac{1}{3},$
(28)

which for $P=f_{1}$ gives $\theta_{A}=62(5)^{\circ}$ Achard:2001uu ;
Achard:2007hm , but $\theta_{A}=84.8(6)^{\circ}$ for $P=\eta$. For this
reason, we only compare the sum of $\eta$ and $\eta^{\prime}$ with the sum of
$f_{1}$ and $f_{1}^{\prime}$ contributions, which are not affected by this
ambiguity. Table 1 shows that, although the CCGDI/HW2 and our model completely
differ in the degrees of freedom that are used to satisfy the relevant SDCs,
they give similar numerical contributions to $a_{\mu}$. The MV model, which
satisfies (8) exactly, much like CCGDI/HW2, but by neglecting any contribution
to $\bar{\Pi}_{1}$ beyond the pion pole in $g-2$ kinematics, gives instead a
much larger contribution.

The breakdown of the contribution to $a_{\mu}$ in different integration
regions shows that there is in general a rather good agreement (with a few
exceptions) between the Regge and the CCDGI/HW2 models. In particular in the
pion/$a_{1}$ channel, the agreement is very good in the “asymptotic” region,
the first row in the table. At low $q^{2}$ there are differences, but these
are expected, because the two models describe different degrees of freedom
there. The situation is similar in the
$\eta/f_{1}+\eta^{\prime}/f_{1}^{\prime}$ channels, where again the largest
differences occur in the low-energy region. However, there are also some non-
negligible differences even in the large-$Q_{i}^{2}$ region, which might be
related to the fact that the HW2 models do not fully saturate the SDC2
Cappiello:2019hwh . Overall, the PS Regge and also the hQCD models are largely
compatible with the LP interpolants, which are independent of the choice of
degrees of freedom. As far as the MV model is concerned, all regions where at
least one of the $Q_{i}$ is large are in reasonable agreement with the other
two models, but it is the region where all $Q_{i}$ are small where the MV
model estimates significantly larger effects; as expected, since in this
region the truncation of the non-factorizable contributions, see Sect. 5.2,
cannot be justified. The table also shows that the kinematic region
$Q_{1}^{2}\sim Q_{2}^{2}\gg Q_{3}^{2}$, $Q_{1,2}^{2}>Q_{\text{match}}$, all
contained in part of the first row and in the second row, provides a small
contribution to the total. This is particularly true for the MV model.

Figure 5: Contribution to $a_{\mu}$ for $Q_{i}\geq Q_{\text{match}}$: the
longitudinal part of the massless perturbative QCD quark loop (dotted red),
the MV model (dot-dashed green), the CCDGI model (dashed and dotted-dashed
magenta), the LR HW2 model (dashed and dotted-dashed turquoise), the LP model
(solid black), and our model (solid blue). The blue point indicates the final
value in (30).

Another way to visualize the impact on $a_{\mu}$ of the different kinematic
regions is to plot the contribution to $a_{\mu}$ as a function of a lower
cutoff $Q_{\text{min}}<Q_{i}$, as shown in Fig. 5. The plot shows again that
although there are differences in the $Q_{i}$ dependence between the Regge and
the hQCD models, these are not so significant with respect to $a_{\mu}$ and
lead to a similar final number. The MV model, on the other hand, only comes
close to the other two for large values of $Q_{\text{min}}$, whereas it
estimates a much larger effect in the region of low $Q_{i}$. Finally, we have
also shown the result obtained with the interpolants by LP, which is
compatible with both the Regge as well as the hQCD models, even though
somewhat lower for low-$Q^{2}$. This may have to do with the fact that it does
not include any explicit resonances and lacks the corresponding low-$Q^{2}$
enhancements. However, as explained in Ludtke:2020moa , the method of
interpolants can be generalized to explicitly include resonance contributions,
once their model-independent description becomes available (and might offer a
valuable alternative to the resummation of a tower of states). For axial-
vector states this is not yet the case, however: a phenomenologically driven
evaluation seems within reach at least for the $f_{1}$ contribution
Zanke:2021wiq , but it will require a detailed understanding of sum-rule
ambiguities.

The present numerical comparison seems to be at odds with the conclusions
drawn by CCDGI Cappiello:2019hwh , who claim to be in agreement with the MV
estimate. They reach this conclusion on the basis of two comparisons: a
detailed one at the level of the $\langle VVA\rangle$ correlation function and
one at the level of the total contribution to $a_{\mu}$. The first one has
been discussed above and indeed shows that the two models agree very well.
However, the comparison of the contribution to $a_{\mu}$ at the level of the
total without separation of the poles of the ground-state pseudoscalars risks
to be misleading: the total number is dominated by the poles due to the
Goldstone bosons and even small differences in the evaluation of the latter
(necessary because of our improved understanding of their TFF) may obscure the
comparison for the remainder. LR Leutgeb:2019gbz , whose model coincides
algebraically with that of CCDGI and numerically differs very little, make the
comparison after first subtracting the Goldstone-boson poles and come to the
same conclusion we reached here.

##  6 Impact of the perturbative corrections to the OPE

We can now evaluate the effect of the recently calculated gluonic corrections
to the OPE Bijnens:2021jqo , as illustrated in Fig. 6. The NLO corrections
lead to a reduction of the massless quark loop that for integrated quantities
tends to evaluate around $1-\alpha_{s}/\pi$, e.g., one has at the symmetric
point $Q_{1}=Q_{2}=Q_{3}=Q$ Bijnens:2021jqo

$\displaystyle\frac{\bar{\Pi}_{1}\big{|}_{\text{NLO}}}{\bar{\Pi}_{1}\big{|}_{\text{LO}}}$
$\displaystyle=C_{1}\frac{\alpha_{s}}{\pi}\;,$ $\displaystyle C_{1}$
$\displaystyle=-\frac{75\Delta^{(1)}-2\Delta^{(3)}-360\zeta_{3}}{54}\approx-0.86\;,$
(29)

where $\Delta^{(n)}=\psi^{(n)}(1/3)-\psi^{(n)}(2/3)$ in terms of the polygamma
function $\psi^{(n)}$ and $\zeta_{3}\approx 1.202$ (the coefficient of the
$\alpha_{s}$ corrections becomes exactly $-1$ in the MV limit Ludtke:2020moa
). In the following, we use the full corrections from Bijnens:2021jqo .

Figure 6: Matching between the NLO OPE and the Regge model for pseudoscalars
(red curve). The blue curve shows the contribution of only the excited
pseudoscalars, excluding the $\pi^{0}$, $\eta$, and $\eta^{\prime}$. The gray
band refers to the uncertainty in setting the $\alpha_{s}$ input, see
Bijnens:2021jqo  for more details.

As the hQCD models and the Regge model agree reasonably well at the numerical
level, we will rely only on the latter in the following. We essentially repeat
the matching procedure described in Colangelo:2019uex , replacing the plain
massless quark loop with the one containing $\mathcal{O}(\alpha_{s})$
corrections. The resulting shift, even down to momentum cutoffs of about
$1\,\text{GeV}$ is small, as can be seen in Fig. 6, but the most relevant
improvement is the reduction of the uncertainties, which had been estimated to
be $\mathcal{O}(20\%)$ of the massless quark loop in Colangelo:2019uex  and is
reduced to a few percent after the NLO calculation Bijnens:2021jqo , much
smaller than the uncertainties on the hadronic side. As a consequence, the
procedure to determine the matching point by minimizing the total
uncertainties would not work anymore. Instead, we keep it fixed at
$1.7(5)\,\text{GeV}$ as we did in Colangelo:2019uex . With these changes, our
updated estimate of the impact of longitudinal SDCs on $a_{\mu}$ reads

$\displaystyle\Delta a_{\mu}^{\text{LSDC}}$
$\displaystyle=\left[8.7(5.3)_{\text{PS-
poles}}+4.2(1)_{q\text{-loop}}\right]\times 10^{-11}$
$\displaystyle=13(5)\times 10^{-11}\;,$ (30)

where the first number in brackets is the contribution from the region below
the matching momentum of $1.7\,\text{GeV}$, evaluated as resummation of
excited pseudoscalars, and the second from the region above $1.7\,\text{GeV}$,
evaluated with the NLO quark loop. A welcome feature of the perturbative
corrections is that they push the OPE curve down, thereby improving the
matching with the hadronic model, which now seems to work perfectly already
around $2\,\text{GeV}$.

If one takes into account the large uncertainties on the hadronic curve (not
shown in Fig. 6) and the rather small estimated perturbative and non-
perturbative uncertainties, one would be led to push the matching point
towards $1\,\text{GeV}$: this would reduce the importance of both hadronic and
model uncertainties and lead to smaller total uncertainties. To illustrate the
point we mention the number we get for a matching point at the lower end of
the range we considered, for $1.2\,\text{GeV}$: $\Delta
a_{\mu}^{\text{LSDC}}=\left[5.7(2.8)_{\text{PS-poles,
par}}+8.1(5)_{q\text{-loop}}\right]\times 10^{-11}=14(3)\times 10^{-11}$,
where the error in the hadronic model only refers to the parametric
uncertainty—reduced from the $3.6\times 10^{-11}$ it contributes to (30)—while
the remainder of the error estimate, especially the variation of the matching
point, does not adapt in a straightforward way to the lower scale. Most
notably, the central value only changes slightly, well within the
uncertainties of the hadronic model, which shows that the information coming
from the perturbative side agrees with the hadronic estimate. Of course, it is
not obvious that at such low energies power corrections beyond the ones
calculated in Bijnens:2020xnl  remain irrelevant, and to be on the safe side
one may consider increasing a bit the uncertainties on the OPE side (also
perturbative corrections, estimated in Bijnens:2021jqo  via scale variation in
$\alpha_{s}$, could play a role at such low energies). In any case, this brief
discussion is just meant to underline the relevance of the calculation of the
corrections to the massless quark loop and their possible impact in reducing
the uncertainties of this contribution: a full implementation of this is left
for future work.

##  7 Conclusions and outlook

In this paper we have discussed our current understanding of the role and
impact of longitudinal SDCs on the HLbL contribution to $(g-2)_{\mu}$ and
updated it to take into account the recent calculation of the NLO perturbative
corrections to the OPE Bijnens:2021jqo . On the low-energy, hadronic side
different solutions for the matching to the SDCs have been proposed, sometimes
accompanied by contradicting statements. To clarify the situation we have
compared these models both at the level of the longitudinal component $w_{L}$
of three-point function $\langle VVA\rangle$ and in terms of the function $G$,
defined in (2), which collects all contributions beyond the pion pole to the
$\bar{\Pi}_{1}$ function of the HLbL tensor—the only one relevant for the
longitudinal SDCs. In this way, the core assumptions and features of each
implementation become most transparent, facilitating the comparison of the
different proposed solutions.

Our conclusions can be summarized as follows:

  1. 1.

Both the original MV model and the recent hQCD models satisfy the axial
anomaly in the chiral limit exactly. When compared at the level of the three-
point function and the longitudinal component $w_{L}(q^{2})$, they agree very
well, supporting that the MV model is an excellent approximation to QCD for
this particular quantity. We have also compared our Regge model for excited
pseudoscalars Colangelo:2019uex  and showed that, as expected, it satisfies
the axial anomaly only asymptotically and at low $q^{2}$.

  2. 2.

A comparison between the MV and the hQCD models for the four-point function,
and in particular the function $G$, can be done analytically and is very
transparent: the MV model can be viewed as a truncation of the hQCD models and
amounts to dropping all contributions beyond the pion pole for $g-2$
kinematics. In the hQCD models these additional contributions are expressed in
terms of a single integral over Bessel functions, which cannot be factorized
into a function of $q_{1,2}^{2}$ and one of $q_{3}^{2}$. We have analyzed the
relative importance of the non-factorizable and the MV term and shown that the
latter is dominant only for rather large values of $q_{1,2}^{2}$: neglecting
the former term leads to a significant overestimate of the low-$q_{1,2}^{2}$
contribution. The hQCD and the MV model, which agree almost exactly on the
axial anomaly, thus differ substantially in their estimate of
$a_{\mu}^{\text{HLbL}}$.

  3. 3.

The two approaches that achieve a matching to the OPE by resumming a tower of
hadronic states provide very similar estimates of the impact on $a_{\mu}$,
even though one is based on excited pseudoscalars in a Regge model
Colangelo:2019uex  and the other on axial-vector mesons in a hQCD model
Leutgeb:2019gbz ; Cappiello:2019hwh , with the aforementioned differences in
$w_{L}(q^{2})$ in the transition region between low and high momenta. This
again shows that the role of the axial anomaly in determining the HLbL
amplitude and its contribution to $a_{\mu}$ is rather limited.

  4. 4.

The hQCD models Leutgeb:2019gbz ; Cappiello:2019hwh  provide an explicit,
analytic solution of the SDC1 in terms of a tower of axial-vector resonances,
which offers useful insights in the mechanism by which SDC1 is fulfilled.
Simple versions of these models, such as HW2, depend on very few parameters,
which can be pinned down by imposing a number of phenomenological constraints,
but once this is done further comparisons to phenomenology show discrepancies.
This can be improved by considering more complicated versions of these models,
such as HW1.

  5. 5.

In the chiral limit the axial-vector mesons have to play an important role in
satisfying the SDCs, and the hQCD models provide a concrete realization of the
underlying mechanism. In the future it will be critical to achieve a full,
model-independent understanding of how axial-vector resonances contribute to
HLbL (in analogy to scalar states Danilkin:2021icn ), at least in the narrow-
width approximation. Otherwise a combination with other contributions to HLbL
scattering evaluated within a dispersive approach would not be justified.
Here, we have presented the expression for the dispersive axial-vector
contribution to $\bar{\Pi}_{1}$ in a particular choice of basis that is
compatible with all contributions evaluated dispersively so far, but sum-rule
ambiguities that are reflected in a basis dependence still need to be
addressed together with a numerical analysis based on any TFF input.

  6. 6.

The final estimates of $\Delta a_{\mu}^{\text{LSDC}}$ obtained with the hQCD
and our Regge model agree quite well with each other as well as with a
solution of the SDCs based on interpolants Ludtke:2020moa . On this basis, we
have updated the final result given in Colangelo:2019uex  to incorporate the
perturbative corrections to the OPE calculated in Bijnens:2021jqo :

$\Delta a_{\mu}^{\text{LSDC}}=13(5)\times 10^{-11}\,.$ (31)

Even with the slight reduction of the total uncertainty, this covers all
realistic estimates of the impact of the longitudinal SDCs on
$a_{\mu}^{\text{HLbL}}$ present in the literature.

As we argued above, further reductions of the uncertainties in the HLbL
contribution due to the fulfillment of the SDCs are possible, also in view of
the smallness of the perturbative corrections to the OPE and their
uncertainties down to $\sim 1\,\text{GeV}$ Bijnens:2021jqo . This will require
an improved and less model-dependent description on the hadronic side before
trying to optimize the matching and exploiting at best the result of the
perturbative calculation. Some of the recent developments discussed here have
paved the way to this goal. Future steps in this direction include:

  * i)

fully clarifying how to evaluate the contribution of axial-vector resonances
to $a_{\mu}^{\text{HLbL}}$ in an unambiguous way;

  * ii)

understanding how to incorporate the solution of the SDCs in the chiral limit
provided by the hQCD models in a more general, dispersively motivated
framework based on axial-vector mesons;

  * iii)

while our discussion here was only concerned with the SDC for the longitudinal
amplitude, a solution in terms of axial vectors can address at the same time
both the longitudinal and the transverse SDCs;

  * iv)

once the treatment of axial-vector mesons in the general dispersive formalism
will become possible, the reasons to use a Regge model of pseudoscalars as a
tool to estimate the impact of the SDCs will cease to exist: only the few
lightest excited pseudoscalars will need to be included.

Work along these lines is ongoing.

###### Acknowledgements.

We thank Jan Lüdtke and Massimiliano Procura for discussions and for providing
numbers about their model. We further thank Antonio Rodríguez-Sánchez for help
with the numerical evaluation of the perturbative correction to the OPE, and
Luigi Cappiello, Oscar Catà, Giancarlo D’Ambrosio, Nils Hermansson-Truedsson,
Abhishek Iyer, Marc Knecht, Josef Leutgeb, Jan Lüdtke, Massimiliano Procura,
and Anton Rebhan for comments on the manuscript. Financial support by the SNSF
(Project Nos. PCEFP2_181117 and PZ00P2_193383 and Grant No. 200020_175791) is
gratefully acknowledged.

## References

  * (1) B. Abi et al. [Muon $g-2$ Collaboration], Phys. Rev. Lett. 126, 141801 (4, 2021), [arXiv:2104.03281 [hep-ex]].
  * (2) T. Albahri et al. [Muon $g-2$ Collaboration], Phys. Rev. D 103, 072002 (2021), [arXiv:2104.03247 [hep-ex]].
  * (3) T. Albahri et al. [Muon $g-2$ Collaboration], Phys. Rev. A 103, 042208 (2021), [arXiv:2104.03201 [hep-ex]].
  * (4) T. Albahri et al. [Muon $g-2$ Collaboration], Phys. Rev. Accel. Beams 24, 044002 (2021), [arXiv:2104.03240 [physics.acc-ph]].
  * (5) T. Aoyama et al., Phys. Rept. 887, 1 (2020), [arXiv:2006.04822 [hep-ph]].
  * (6) T. Aoyama, M. Hayakawa, T. Kinoshita, and M. Nio, Phys. Rev. Lett. 109, 111808 (2012), [arXiv:1205.5370 [hep-ph]].
  * (7) T. Aoyama, T. Kinoshita, and M. Nio, Atoms 7, 28 (2019).
  * (8) A. Czarnecki, W. J. Marciano, and A. Vainshtein, Phys. Rev. D 67, 073006 (2003), [arXiv:hep-ph/0212229], [Erratum: Phys. Rev. D 73, 119901 (2006)].
  * (9) C. Gnendiger, D. Stöckinger, and H. Stöckinger-Kim, Phys. Rev. D 88, 053005 (2013), [arXiv:1306.5546 [hep-ph]].
  * (10) M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Eur. Phys. J. C 77, 827 (2017), [arXiv:1706.09436 [hep-ph]].
  * (11) A. Keshavarzi, D. Nomura, and T. Teubner, Phys. Rev. D 97, 114025 (2018), [arXiv:1802.02995 [hep-ph]].
  * (12) G. Colangelo, M. Hoferichter, and P. Stoffer, JHEP 02, 006 (2019), [arXiv:1810.00007 [hep-ph]].
  * (13) M. Hoferichter, B.-L. Hoid, and B. Kubis, JHEP 08, 137 (2019), [arXiv:1907.01556 [hep-ph]].
  * (14) M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Eur. Phys. J. C 80, 241 (2020), [arXiv:1908.00921 [hep-ph]], [Erratum: Eur. Phys. J. C 80, 410 (2020)].
  * (15) A. Keshavarzi, D. Nomura, and T. Teubner, Phys. Rev. D 101, 014029 (2020), [arXiv:1911.00367 [hep-ph]].
  * (16) B.-L. Hoid, M. Hoferichter, and B. Kubis, Eur. Phys. J. C 80, 988 (2020), [arXiv:2007.12696 [hep-ph]].
  * (17) A. Kurz, T. Liu, P. Marquard, and M. Steinhauser, Phys. Lett. B 734, 144 (2014), [arXiv:1403.6400 [hep-ph]].
  * (18) K. Melnikov and A. Vainshtein, Phys. Rev. D 70, 113006 (2004), [arXiv:hep-ph/0312226].
  * (19) P. Masjuan and P. Sánchez-Puertas, Phys. Rev. D 95, 054026 (2017), [arXiv:1701.05829 [hep-ph]].
  * (20) G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, Phys. Rev. Lett. 118, 232001 (2017), [arXiv:1701.06554 [hep-ph]].
  * (21) G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, JHEP 04, 161 (2017), [arXiv:1702.07347 [hep-ph]].
  * (22) M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, Phys. Rev. Lett. 121, 112002 (2018), [arXiv:1805.01471 [hep-ph]].
  * (23) M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, JHEP 10, 141 (2018), [arXiv:1808.04823 [hep-ph]].
  * (24) A. Gérardin, H. B. Meyer, and A. Nyffeler, Phys. Rev. D 100, 034520 (2019), [arXiv:1903.09471 [hep-lat]].
  * (25) J. Bijnens, N. Hermansson-Truedsson, and A. Rodríguez-Sánchez, Phys. Lett. B 798, 134994 (2019), [arXiv:1908.03331 [hep-ph]].
  * (26) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, and P. Stoffer, Phys. Rev. D 101, 051501 (2020), [arXiv:1910.11881 [hep-ph]].
  * (27) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, and P. Stoffer, JHEP 03, 101 (2020), [arXiv:1910.13432 [hep-ph]].
  * (28) T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Phys. Rev. Lett. 124, 132002 (2020), [arXiv:1911.08123 [hep-lat]].
  * (29) G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera, and P. Stoffer, Phys. Lett. B 735, 90 (2014), [arXiv:1403.7512 [hep-ph]].
  * (30) G. W. Bennett et al. [Muon $g-2$ Collaboration], Phys. Rev. D 73, 072003 (2006), [arXiv:hep-ex/0602035].
  * (31) J. Grange et al. [Muon $g-2$ Collaboration], arXiv:1501.06858 [physics.ins-det].
  * (32) M. Abe et al., PTEP 2019, 053C02 (2019), [arXiv:1901.03047 [physics.ins-det]].
  * (33) E.-H. Chao, R. J. Hudspith, A. Gérardin, J. R. Green, H. B. Meyer, and K. Ottnad, Eur. Phys. J. C 81, 651 (2021), [arXiv:2104.02632 [hep-lat]].
  * (34) S. L. Adler, Phys. Rev. 177, 2426 (1969).
  * (35) J. S. Bell and R. Jackiw, Nuovo Cim. A 60, 47 (1969).
  * (36) W. A. Bardeen, Phys. Rev. 184, 1848 (1969).
  * (37) J. Wess and B. Zumino, Phys. Lett. B 37, 95 (1971).
  * (38) E. Witten, Nucl. Phys. B 223, 422 (1983).
  * (39) A. Vainshtein, Phys. Lett. B 569, 187 (2003), [arXiv:hep-ph/0212231].
  * (40) M. Knecht, S. Peris, M. Perrottet, and E. de Rafael, JHEP 03, 035 (2004), [arXiv:hep-ph/0311100].
  * (41) J. Bijnens, N. Hermansson-Truedsson, L. Laub, and A. Rodríguez-Sánchez, JHEP 10, 203 (2020), [arXiv:2008.13487 [hep-ph]].
  * (42) J. Bijnens, N. Hermansson-Truedsson, L. Laub, and A. Rodríguez-Sánchez, JHEP 04, 240 (2021), [arXiv:2101.09169 [hep-ph]].
  * (43) J. Leutgeb and A. Rebhan, Phys. Rev. D 101, 114015 (2020), [arXiv:1912.01596 [hep-ph]].
  * (44) L. Cappiello, O. Catà, G. D’Ambrosio, D. Greynat, and A. Iyer, Phys. Rev. D 102, 016009 (2020), [arXiv:1912.02779 [hep-ph]].
  * (45) J. Lüdtke and M. Procura, Eur. Phys. J. C 80, 1108 (2020), [arXiv:2006.00007 [hep-ph]].
  * (46) S. Peris, M. Perrottet, and E. de Rafael, JHEP 05, 011 (1998), [arXiv:hep-ph/9805442].
  * (47) J. Bijnens, E. Gámiz, E. Lipartia, and J. Prades, JHEP 04, 055 (2003), [arXiv:hep-ph/0304222].
  * (48) M. Hoferichter, G. Colangelo, M. Procura, and P. Stoffer, Int. J. Mod. Phys. Conf. Ser. 35, 1460400 (2014), [arXiv:1309.6877 [hep-ph]].
  * (49) G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, JHEP 09, 091 (2014), [arXiv:1402.7081 [hep-ph]].
  * (50) G. Colangelo, M. Hoferichter, B. Kubis, M. Procura, and P. Stoffer, Phys. Lett. B 738, 6 (2014), [arXiv:1408.2517 [hep-ph]].
  * (51) G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, JHEP 09, 074 (2015), [arXiv:1506.01386 [hep-ph]].
  * (52) K. Melnikov and A. Vainshtein, arXiv:1911.05874 [hep-ph].
  * (53) J. Bijnens and J. Prades, Mod. Phys. Lett. A 22, 767 (2007), [arXiv:hep-ph/0702170].
  * (54) M. Knecht, JHEP 08, 056 (2020), [arXiv:2005.09929 [hep-ph]].
  * (55) G. P. Lepage and S. J. Brodsky, Phys. Lett. B 87, 359 (1979).
  * (56) G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22, 2157 (1980).
  * (57) P. Masjuan, P. Roig, and P. Sánchez-Puertas, arXiv:2005.11761 [hep-ph].
  * (58) I. Danilkin, M. Hoferichter, and P. Stoffer, Phys. Lett. B 820, 136502 (2021), [arXiv:2105.01666 [hep-ph]].
  * (59) P. Roig and P. Sánchez-Puertas, Phys. Rev. D 101, 074019 (2020), [arXiv:1910.02881 [hep-ph]].
  * (60) V. Pauk and M. Vanderhaeghen, Eur. Phys. J. C 74, 3008 (2014), [arXiv:1401.0832 [hep-ph]].
  * (61) F. Jegerlehner, Springer Tracts Mod. Phys. 274, 1 (2017).
  * (62) R. García-Martín and B. Moussallam, Eur. Phys. J. C 70, 155 (2010), [arXiv:1006.5373 [hep-ph]].
  * (63) M. Hoferichter, D. R. Phillips, and C. Schat, Eur. Phys. J. C 71, 1743 (2011), [arXiv:1106.4147 [hep-ph]].
  * (64) B. Moussallam, Eur. Phys. J. C 73, 2539 (2013), [arXiv:1305.3143 [hep-ph]].
  * (65) I. Danilkin and M. Vanderhaeghen, Phys. Lett. B 789, 366 (2019), [arXiv:1810.03669 [hep-ph]].
  * (66) M. Hoferichter and P. Stoffer, JHEP 07, 073 (2019), [arXiv:1905.13198 [hep-ph]].
  * (67) I. Danilkin, O. Deineka, and M. Vanderhaeghen, Phys. Rev. D 101, 054008 (2020), [arXiv:1909.04158 [hep-ph]].
  * (68) I. Danilkin, O. Deineka, and M. Vanderhaeghen, Phys. Rev. D 96, 114018 (2017), [arXiv:1709.08595 [hep-ph]].
  * (69) J. Lu and B. Moussallam, Eur. Phys. J. C 80, 436 (2020), [arXiv:2002.04441 [hep-ph]].
  * (70) M. Hoferichter and P. Stoffer, JHEP 05, 159 (2020), [arXiv:2004.06127 [hep-ph]].
  * (71) P. Achard et al. [L3 Collaboration], Phys. Lett. B 526, 269 (2002), [arXiv:hep-ex/0110073].
  * (72) M. Zanke, M. Hoferichter, and B. Kubis, JHEP 07, 106 (2021), [arXiv:2103.09829 [hep-ph]].
  * (73) L. Gan, B. Kubis, E. Passemar, and S. Tulin, arXiv:2007.00664 [hep-ph].
  * (74) P. Achard et al. [L3 Collaboration], JHEP 03, 018 (2007).


11institutetext: {daniel.zentai1, mihail.plesa2<EMAIL_ADDRESS>
xtendr
Budapest, Hungary

# A Multiparty Commutative Hashing Protocol based on the Discrete Logarithm
Problem

Daniel Zentai1 and Mihail Plesa2 and Robin Frot3

###### Abstract

Let $\mathcal{X}$ and $\mathcal{Y}$ be two sets and suppose that a set of
participants $P=\\{P_{1},P_{2},\dots,P_{n}\\}$ would like to calculate the
keyed hash value of some message $m\in\mathcal{X}$ known to a single
participant in $P$ called the data owner. Also, suppose that each participant
$P_{i}$ knows a secret value $x_{i}\in\mathcal{X}$. In this paper, we will
propose a protocol that enables the participants in this setup to calculate
the value $y=H(m,x_{1},x_{2},\dots,x_{n})$ of a hash function
$H:\mathcal{X}^{n+1}\rightarrow\mathcal{Y}$ such that:

  * •

The function $H$ is a one-way function.

  * •

Participants in $P\backslash\\{P_{i}\\}$ cannot obtain $x_{i}$.

  * •

Participants other than the data owner cannot obtain $m$.

  * •

The hash value $y=H(m,x_{1},x_{2},\dots,x_{n})$ remains the same regardless
the order of the secret $x_{i}$ values.

###### Keywords:

Hash functions, Discrete logarithm problem, Anonymization

##  1 Introduction

Hash functions are very common building blocks of cryptographic protocols like
digital signatures or message authentication codes. In this paper, we will
propose a protocol built around the Chaum-van Heijst-Pfitzmann hash function
[1]. The goal of our protocol is to calculate a hash value in a multiparty
setup, i.e. the calculation is made by multiple participants collaboratively.

Our motivation was the following. Suppose we have to store personal data in a
way that enables us to have access to a certain subset of non-sensitive
attributes (e.g. age or height) and also enables us to keep the
confidentiality of more sensitive attributes (e.g. passport number or name).
Anonymization and pseudonymization may look like very similar concepts, but
according to the European Union’s General Data Protection Regulation (GDPR
[8]) this is not the case. Pseudonymized data can be recovered using some
secret information (e.g. a decryption key). Anonymized data on the other hand
cannot be recovered under any circumstances, meaning we cannot use encryption
to anonymize data.

In this paper, we will propose a protocol that enables a set of participants
and a trusted server to calculate a hash value (i.e. anonymize their data)
collaboratively.

We impose the following requirements regarding the process of the calculation.
The function $H$ should be a one-way function, the participants (including the
server) should not learn either each other’s secret values or the plaintext,
and the hash value should be the same regardless of which participant is the
data owner (i.e. what is the order of the participants in the calculation).
The reason for the last requirement is that we want to be able to detect
whether two different participants have the same plaintext without knowing the
plaintext itself.

##  2 Preliminaries

In this section, we will give some basic definitions required throughout our
paper. For the interested reader, we recommend [5] for more information.

###### Definition 1

A function $\epsilon:\mathbb{N}\rightarrow\mathbb{R}$ is negligible if for all
$c\in\mathbb{N}$ there exists an $n_{0}\in\mathbb{N}$ such that for all $n\geq
n_{0}$ we have $\epsilon(n)<n^{-c}.$

###### Definition 2

A hash function $h$ is a deterministic polynomial time algorithm with
arbitrary input length and fixed output length.

###### Definition 3

A hash function $h$ is said to be one-way if it is

  * •

Easy to compute, i.e. $h$ outputs $h(x)$ in polynomial time for all $x\in
D(h)$.

  * •

Hard to invert, i.e. given a hash value $y$ no probabilistic polynomial time
algorithm can calculate an $x$ for which $y=h(x)$ with non-negligible
probability.

###### Definition 4

A hash function is said to be collision resistant if no probabilistic
polynomial time algorithm can succeed in finding two values $x\neq
x^{\prime}\in D(h)$ such that $h(x)=h(x^{\prime})$ with non-negligible
probability.

Throughout this paper, we will assume that the discrete logarithm problem
defined below is hard, i.e. it cannot be solved by any probabilistic
polynomial time algorithm.

###### Definition 5

If $G$ is a cyclic group and $a,b\in G$ are group elements then the discrete
logarithm problem includes finding an $x\in G$ such that $a^{x}=b$.

There are several cryptosystems based on the hardness of the discrete
logarithm problem. The most popular ones include the Diffie-Hellman key
exchange [3] and the ElGamal encryption [4].

##  3 Related Work

There are several anonymization techniques available, including but not
limited to $k$-anonymity, differential privacy, or synthetic data [6]. These
methods are designed either to add noise to the data, add some artificial data
to the real data set, or generalize the data such that the data subject is no
longer identifiable from it.

Our method is a bit different in a few aspects. First of all, we propose a
purely cryptographic method, the protocol requires the collaboration of
multiple participants, and the anonymized (i.e. hashed) version of a fixed
data record remains the same regardless of which participant owns the
plaintext. The main advantage of this approach is that this way we allow the
participants to calculate the intersection of their data set without revealing
sensitive identifiers.

To ensure that the hash of a fixed message $m$ remains the same regardless of
which participant owns the plaintext, we have to use commutative hashing. In
[2] the authors propose $h(x,y)=f(\min(x,y),\max(x,y))$ where $f$ is a
collision-resistant hash function. Unfortunately, a similar approach won’t
work in our case since the participants know nothing about each other’s secret
values.

###### Definition 6

A hash function $h:\mathcal{X}^{n}\rightarrow\mathcal{Y}^{n}$ is commutative,
if $h(x_{1},x_{2},\dots,x_{n})=h(\sigma(x_{1},x_{2},\dots,x_{n}))$ for all
$x_{1},x_{2},\dots,x_{n}$ and all $\sigma$ permutations.

Note that we cannot ensure collision resistance if we use a commutative hash
function, since $(x_{1},x_{2},\dots,x_{n})\neq\sigma(x_{1},x_{2},\dots,x_{n})$
(unless $\sigma$ is the identity permutation), but
$h(x_{1},x_{2},\dots,x_{n})=h(\sigma(x_{1},x_{2},\dots,x_{n}))$ by definition.
Therefore the most we can do without losing the desired functionality is to
ensure one-way property.

##  4 Hashing Protocol

As a building block of our protocol, we will use the discrete logarithm-based
collision-resistant hash function defined by Chaum, van Heijst, and Pfitzmann
[1]

###### Definition 7

Let $p$ be a prime number such that $q=\frac{p-1}{2}$ is also a prime number,
and let $a$ and $b$ be primitive elements of the multiplicative group of
$GF(p)$. Let $h:GF(q)\times GF(q)\rightarrow GF(p)\backslash\\{0\\}$ be the
following function:

$h(x,y)=a^{x}\cdot b^{y}\mod p$

The hash function $h$ defined above is collision-resistant if the discrete
logarithm problem is hard [1].

Now, suppose that we have a set of participants
$P=\\{P_{1},P_{2},\dots,P_{n}\\}$ and a trusted server represented by a
special participant $S\notin P$. Also, suppose that the data owner (i.e. the
only participant who knows $m$ is $P_{1}$, therefore $P_{1}$ initiates the
hashing protocol).

Let $h$ be the Chaum-van Heijst-Pfitzmann hash function and
$\Pi=(Gen,Enc,Dec)$ a public key encryption scheme. Suppose that each
participant $P_{i}$ has two randomly generated keys $x_{i},y_{i}\in GF(q)$.
Also, suppose that $P_{1}$ wants to hash the message $m\in GF(q)$. In
addition, we will suppose that all the messages in the protocol are sent over
a secure channel.

Our protocol proceeds as follows:

Multiparty Commutative Hashing Protocol
1. Upon receiving an upload request from $P_{1}$, $S$ generates $n$ random numbers $r_{1},r_{2},\dots,r_{n}$ and sends them to the participants in $P$. 2. Upon receiving $r_{1}$ from $S$, $P_{1}$ calculates $h_{1}=h(x_{1}+m,y_{1})$ and sends $h_{1}|Enc_{K^{pub}_{S}}(r_{1})$ to $S$, where $K^{pub}_{S}$ is the public key of $S$. 3. For all $i=2,\dots,n$, upon receiving $r_{i}$ from $S$, $P_{i}$ calculates $h_{i}=h(x_{i},y_{i})$ and sends $h_{i}|Enc_{K^{pub}_{S}}(r_{i})$ to $S$. 4. Upon receiving $h_{i}|Enc_{K^{pub}_{S}}(r_{i})$ from $P_{i}$ for all $i=1,\dots,n$, $S$ checks the random values $r_{1},r_{2},\dots,r_{n}$ and stores $\prod_{j=1}^{n}h_{j}$ if these are the same values $S$ sent out at the first step. Otherwise, $S$ outputs an error message.

Note that the purpose of the $r_{i}$ values is not authentication, these just
make it possible for the server to check that every participant took part in
the calculation and therefore $\prod_{j=1}^{n}h_{j}$ can be stored as the
output of the protocol.

###### Theorem 4.1

If the discrete logarithm problem is hard, then
$H(x_{1},\dots,x_{n},y_{1},\dots,y_{n})=\prod_{i=1}^{n}h(x_{i},y_{i})$ is a
one-way function.

Proof. Recall that if the discrete logarithm problem is hard, then $h$ is a
collision-resistant hash function. Suppose that $H$ is not a one-way function,
thus there exists a PPT algorithm $\mathcal{A}$ that can successfully find a
corresponding input to a given value of $H$. We will use $\mathcal{A}$ to find
a collision in $h$.

From the construction of $H$, it immediately follows that

$\displaystyle H(x_{1},\dots,x_{n},y_{1},\dots,y_{n})$
$\displaystyle=\prod_{i=1}^{n}h(x_{i},y_{i})$ (1)
$\displaystyle=\prod_{i=1}^{n}a^{x_{i}}\cdot b^{y_{i}}$
$\displaystyle=a^{\sum_{i=1}^{n}x_{i}}\cdot b^{\sum_{i=1}^{n}y_{i}}$
$\displaystyle=h\left(\sum_{i=1}^{n}x_{i},\sum_{i=1}^{n}y_{i}\right)$

Now let $\alpha=H(x,0,\dots,0,y,0,\dots,0)$ and run algorithm $\mathcal{A}$ on
$\alpha$. Suppose that
$\mathcal{A}(\alpha)=x_{1}^{*},\dots,x_{n}^{*},y_{1}^{*},\dots,y_{n}^{*}$.

From this, we can successfully find a preimage of $h(x,y)$ since

$\displaystyle h(x,y)$ $\displaystyle=H(x,0,\dots,0,y,0,\dots,0)$ (2)
$\displaystyle=H(x_{1}^{*},\dots,x_{n}^{*},y_{1}^{*},\dots,y_{n}^{*})$
$\displaystyle=h\left(\sum_{i=1}^{n}x_{i}^{*},\sum_{i=1}^{n}y_{i}^{*}\right)$

And that is a contradiction since $h$ is collision-resistant (therefore also a
one-way function) if the discrete logarithm problem is hard.

$\square$

Note that $H$ is still collision-resistant in a somewhat weaker sense. If we
fix the $x_{i}$ and $y_{i}$ values (e.g. we do not let the participants to
change their keys) then it is hard to find two messages $m\neq m^{\prime}$
with the same hash value.

##  5 Hashing with elliptic curves

Although the Chaum-van Heijst-Pfitzmann hash is a provable secure
construction, it is not widely used in practice because of its running time.
One way to make a discrete logarithm-based protocol more applicable in
practice, is to use elliptic curves.

###### Definition 8

Let $p$ be a prime number and $\mathcal{E}_{p}$ be an elliptic curve over the
finite field $GF(p)$. Also, let $P,Q$ be points of $\mathcal{E}_{p}$. The
elliptic curve discrete logarithm problem includes finding a $k\in\mathbb{Z}$
such that $k\cdot P=Q$.

The elliptic curve discrete logarithm problem is thought to be even harder
than the aforementioned classical discrete logarithm problem. The Chaum-van
Heijst-Pfitzmann hash function can be defined as follows using elliptic
curves.

###### Definition 9

Let $p$ be a prime number and $\mathcal{E}_{p}$ be an elliptic curve over the
finite field $GF(p)$, moreover assume that $|\mathcal{E}_{p}|$ is a prime $n$.
Let $A$ and $B$ be base points of $\mathcal{E}_{p}$ such that the discrete
logarithm of $B$ is not known. Let $h:(\mathbb{Z}\backslash
n\mathbb{Z})^{2}\rightarrow\mathcal{E}_{p}\backslash\\{\mathcal{O}\\}$ be the
following function:

$h(k,l)=k\cdot A+l\cdot B$

###### Theorem 5.1

The hash function $h$ defined above is collision-resistant if the elliptic
curve discrete logarithm problem is hard.

Proof. Suppose that we can find a collision efficiently, i.e. two pairs of
integers $(k,l)\neq(k^{\prime},l^{\prime})$ such that

$h(k,l)=k\cdot A+l\cdot B=k^{\prime}\cdot A+l^{\prime}\cdot
B=h(k^{\prime},l^{\prime}).$

We can assume that $l^{\prime}\neq l$. Otherwise, we would have
$kA=k^{\prime}A$ and thus $k=k^{\prime}$. Since $n$ is prime, $l-l^{\prime}$
is invertible modulo $n$ and this gives

$\frac{k-k^{\prime}}{l^{\prime}-l}\cdot A=B.$

We successfully solved the elliptic curve discrete logarithm problem.

$\square$

Again, suppose that we have a set of participants
$P=\\{P_{1},P_{2},\dots,P_{n}\\}$ and a trusted server represented by a
special participant $S\notin P$. Also, suppose that the data owner (i.e. the
only participant who knows $m$ is $P_{1}$, therefore $P_{1}$ initiates the
hashing protocol.

Let $h$ be the Chaum-van Heijst-Pfitzmann hash function and
$\Pi=(Gen,Enc,Dec)$ a public key encryption scheme. Suppose that each
participant $P_{i}$ has two randomly generated keys
$k_{i},l_{i}\in\mathbb{Z}$. Also, suppose that $P_{1}$ wants to hash the
message $m$ where $m$ is a point of an elliptic curve $\mathcal{E}_{p}$ over
the finite field $GF(p)$

Now our modified protocol proceeds as follows:

Multiparty Anonymization Protocol with Elliptic Curves
1. Upon receiving an upload request from $P_{1}$, $S$ generates $n$ random numbers $r_{1},r_{2},\dots,r_{n}$ and sends them to the participants in $P$. 2. Upon receiving $r_{1}$ from $S$, $P_{1}$ calculates $h_{1}=h(k_{1}+m,l_{1})$ and sends $h_{1}|Enc_{K^{pub}_{S}}(r_{1})$ to to $S$, where $K^{pub}_{S}$ is the public key of $S$. 3. For all $i=2,\dots,n$, upon receiving $r_{i}$ from $S$, $P_{i}$ calculates $h_{i}=h(k_{i},l_{i})$ and sends $h_{i}|Enc_{K^{pub}_{S}}(r_{i})$ to $S$. 4. Upon receiving $h_{i}|Enc_{K^{pub}_{S}}(r_{i})$ from $P_{i}$ for all $i=1,\dots,n$, $S$ checks the random values $r_{1},r_{2},\dots,r_{n}$ and stores $\sum_{j=1}^{n}h_{j}$ if these are the same values $S$ sent out at the first step. Otherwise, $S$ outputs an error message.

##  6 Implementation results

We instantiated the protocol both over the multiplicative group of integers
modulo $p$ and over the elliptic curve SECP256k1 [13]. The purpose of our
experiment was to analyze the running time with respect to the number of
participants. We discovered that the elliptic curve implementation is
$10\times$ faster than the implementation over the multiplicative group of
integers. For each number of participants, we run $100$ trials of the protocol
and average the running time. The results are presented in Table 1. We used
the Levenberg-Marquardt algorithm to determine the coefficients of a linear
polynomial that best fit the data [9]. Let $N$ be the number of participants.
We denote by $T_{\mathbb{Z}_{p}}\left(N\right)$ and
$T_{\mathcal{E}_{p}}\left(N\right)$ the running time in seconds of the
protocol implemented over $\mathbb{Z}_{p}$, respectively $\mathcal{E}_{p}$.
The discovered polynomials are shown in (4) and (5):

$T_{\mathcal{E}_{p}}\left(N\right)=0.008N-0.733$ (3)
$T_{\mathbb{Z}_{p}}\left(N\right)=0.13N-42.06$ (4)

Num. of participants |  Running time over $\mathcal{E}_{p}$ (s) |  Running time over $\mathbb{Z}_{p}$ (s)
---|---|---
4 |  0.028 |  0.280
8 |  0.058 |  0.561
16 |  0.111 |  1.122
32 |  0.222 |  2.245
64 |  0.442 |  4.511
128 |  0.882 |  9.065
256 |  1.772 |  18.318
512 |  3.551 |  37.178
1024 |  7.135 |  76.174
2048 |  14.431 |  160.599
4096 |  29.474 |  352.496
8192 |  61.157 |  837.490
16384 |  132.225 |  2242.677
Table 1: The running time (s) with respect to the number of participants

Figures 1 and 2 show the running time with respect to the number of
participants as well as the polynomial found by the curve fitting algorithm
when implemented over $\mathcal{E}_{p}$ and over $\mathbb{Z}_{p}$.

All experiments were performed in Python on an Apple M1 Max platform.

Figure 1: The running time of the protocol implemented over $\mathcal{E}_{p}$
Figure 2: The running time of the protocol implemented over $\mathbb{Z}_{p}$

##  7 A threshold protocol

One possible relaxation of the protocol is to ensure that any $k$-element
subset of the participants can calculate the anonymous ID for some fixed
$1<k\leq n=|P|$. As a building block, we can use Shamir’s secret sharing
scheme [10] which works as follows.
Let $\mathbb{F}$ be a finite field and let $s_{0}\in\mathbb{F}$ be the secret
we want to share with the participants. Let $a_{1},\dots,a_{k-1}$ be elements
of $\mathbb{F}$ chosen uniformly at random and
$f(x)=s_{0}+\sum_{i=1}^{k-1}a_{i}x^{i}$. Each participant $P_{i}$ is given a
point $(x_{i},f(x_{i}))$ of the polynomial $f$. With Lagrange interpolation,
any $k$ element subset of the participants can now calculate the secret
$s_{0}$ as follows.

$s_{0}=f(0)=\sum_{i=1}^{k}f(x_{i})\cdot\ell_{i}$

where

$\ell_{i}=\prod_{j=1,j\neq i}^{k}\frac{x_{j}}{x_{j}-x_{i}}$

Now we can use this secret sharing scheme to modify our anonymization protocol
in a way that any $k$ element subset of $P$ can calculate the anonymous IDs.
Let $h$ be the Chaum-van Heijst-Pfitzmann hash function and
$\Pi=(Gen,Enc,Dec)$ a Fully Homomorphic Encryption (FHE) scheme. Suppose that
each participant $P_{i}$ has a randomly generated key $x_{i}\in GF(q)$. Also,
suppose that $P_{1}$ wants to hash the message $m\in GF(q)$. In addition, we
will suppose that all the messages in the protocol are sent over a secure
channel.

Multiparty Threshold Anonymization Protocol
1. First, $S$ generates two secret $s_{0}$ and $t_{0}$ and distributes them using Shamir’s secret sharing scheme to the members in $P$. $S$ uses two different polynomials $f(x)=s_{0}+\sum_{i=1}^{k-1}a_{i}x^{i}$ and $g(x)=t_{0}+\sum_{i=1}^{k-1}b_{i}x^{i}$. Let $x_{1},\dots,x_{n}\in GF(p)$ secret random elements generated by the respective $P_{i}$. 2. Upon receiving $Enc_{K^{pub}_{i}}(x_{i})$ from $P_{i}$, $S$ sends $f(Enc_{K^{pub}_{i}}(x_{i}))=Enc_{K^{pub}_{i}}(f(x_{i}))$ and $g(Enc_{K^{pub}_{i}}(x_{i}))=Enc_{K^{pub}_{i}}(g(x_{i}))$ to $P_{i}$ where $K^{pub}_{i}$ is the public key of $P_{i}$. 3. Every pair $(P_{i},P_{i+1})$ shares with $S$ the quotient $x_{i+1}/x_{i}$ with the Multiply protocol described below. 4. Upon receiving an anonymization request from $P_{1}$, $S$ randomly chooses a $k$ element subset $Q\subseteq P$. Without loss of generality suppose now that $Q=\\{P_{1},\dots,P_{k}\\}$. Additionally, $S$ generates $k$ random values $r_{1},r_{2}\dots,r_{k}$ and sends $r_{i}$ to $P_{i}\in Q$. 5. For every participants $P_{i}\in Q$, $S$ calculates $\ell_{i}=\prod_{j=1,j\neq i}^{k}\frac{x_{j}}{x_{j}-x_{i}}$ with the stored quotients and sends $\ell_{i}$ to $P_{i}$. 6. Upon receiving $r_{1}$ and $\ell_{1}$ from $S$, $P_{1}$ calculates $h_{1}=h(m+f(x_{1})\cdot\ell_{1},g(x_{1})\cdot\ell_{1})$ and sends $h_{1}|Enc_{K_{S}^{pub}}(r_{1})$ to $S$ where $K_{S}^{pub}$ is the public key of $S$. 7. For all $i=2,\dots,k$, upon receiving $r_{i}$ from $S$, $P_{i}$ calculates $h_{i}=h(f(x_{i})\cdot\ell_{i},g(x_{i})\cdot\ell_{i})$ and sends $h_{i}|Enc_{K^{pub}_{S}}(r_{i})$ to $S$. 8. Upon receiving $h_{i}|Enc_{K^{pub}_{S}}(r_{i})$ from $P_{i}$ for all $i=1,\dots,k$, $S$ stores $\prod_{i=1}^{k}h_{i}=a^{m+\sum_{i=1}^{k}f(x_{i})\cdot\ell_{i}}\cdot b^{\sum_{i=1}^{k}g(x_{i})\cdot\ell_{i}}=a^{m+s_{0}}\cdot b^{t_{0}}$ if the random values $r_{1},r_{2},\dots,r_{k}$ are the same values $S$ sent out at the first step. Otherwise, $S$ outputs an error message.

Note that in step $5$, the server $S$ needs to calculate
$\frac{x_{i}}{x_{i}-x_{j}}$, where $x_{j}$ is not known by $P_{i}$ since it is
a secret value generated by $P_{j}$. We will use a subroutine that solves this
problem. The subroutine uses the fact that $\frac{x_{i}}{x_{i}-x_{j}}$ can be
calculated securely if we can calculate the product $x_{k}^{-1}\cdot x_{k+1}$
securely, since

$\frac{x_{i}}{x_{i}-x_{j}}=\left(\frac{x_{i}-x_{j}}{x_{i}}\right)^{-1}=\left(1-\frac{x_{j}}{x_{i}}\right)^{-1}=\left(1-x_{i}^{-1}\cdot
x_{j}\right)^{-1}$

and, if $j>i$,

$\frac{x_{j}}{x_{i}}=\prod_{k=i}^{j-1}\frac{x_{k+1}}{x_{k}}.$

Two-Party Multiplication Protocol
1. Multiply($x,y$) takes two inputs from two different parties: $x$ from $P_{1}$ and $y$ from $P_{2}$. 2. $P_{1}$ generates a random number $r_{1}$ and sends $r_{1}\cdot x$ to $P_{2}$ who then generates a random number $r_{2}$, calculates $r_{1}\cdot x\cdot r_{2}\cdot y$ and sends it to the server. 3. The server also generates a random value $r_{S}$, calculates $r_{S}\cdot r_{1}\cdot x\cdot r_{2}\cdot y$ and sends it back to $P_{1}$. 4. $P_{1}$ calculates $r^{-1}_{1}\cdot r_{S}\cdot r_{1}\cdot x\cdot r_{2}\cdot y=r_{S}\cdot x\cdot r_{2}\cdot y$ and sends it to $P_{2}$. 5. $P_{2}$ calculates $r^{-1}_{2}\cdot r_{S}\cdot x\cdot r_{2}\cdot y=r_{S}\cdot x\cdot y$ and sends it to the server. 6. The server calculates $r^{-1}_{S}\cdot r_{S}\cdot x\cdot y=x\cdot y$.

Note that in both versions of the protocol, if the server is corrupted, it can
brute-force the hash values. Indeed, if $P_{i}$ is the data owner during the
first run then $S$ learns $h(x_{j},y_{j})$ (or already knows $(s_{0},t_{0})$
in the threshold protocol) for all $j\in P\backslash\\{P_{i}\\}$. And if
$P_{u}\neq P_{i}$ is the data owner during the second run then $S$ learns
$h(x_{i},y_{j})$ (or $h(f(x_{i}),g(x_{i})$ in the threshold protocol) as well.
We can prevent this if the participants agree on a random number $R$ using
some group key exchange protocol and run the protocol on the message $R\cdot
m$ instead of $m$. Other participants besides $S$ cannot perform this attack
since every message is sent over a secure channel.

Also note that it is possible to double the size of the plaintext space if
$P_{1}$ sends $h(m_{1}+x_{1},m_{2}+y_{1})$ (or
$h(m_{1}+f(x_{1})\cdot\ell_{1},m_{2}+g(x_{1})\cdot\ell_{1})$ in the threshold
case) with $(m_{1},m_{2})\in GF(p)^{2}$.

##  8 Conclusion

We proposed a protocol that enables a set of participants to calculate a one-
way commutative hash function collaboratively in a way that all of the
participants have to take part in the computation, and this fact is verifiable
by a trusted server and no participant other than the data owner has access to
the plaintext. We also proposed a threshold protocol where any $k$-element
subset of the participants can run the protocol successfully.

A possible future direction would be to construct a decentralized version of
the protocol or make the protocol post-quantum.

Our implementation was purely experimental and was not made for efficiency.
One possible direction of research is to create an efficient implementation of
the protocol using low-level parallelizable instruction or FPGA [11][12].

## References

  * [1] D. Chaum, E. Van Heijst, and B. Pfitzmann. Cryptographically strong undeniable signatures, unconditionally secure for the signer. Advances in Cryptology-CRYPTO ’91, volume 576 of Lecture Notes in Computer Science, pages 470-484. Springer-Verlag, (1992).
  * [2] Goodrich, Michael T., and Roberto Tamassia. Efficient authenticated dictionaries with skip lists and commutative hashing. Technical Report, Johns Hopkins Information Security Institute, 2000.
  * [3] W. Diffie, M. E. Hellman. ”New Directions in Cryptography. IEEE Transactions on Information Theory. 22 (6): 644–654 (1976).
  * [4] T. ElGamal. A Public-Key Cryptosystem and a Signature Scheme Based on Discrete Logarithms. IEEE Transactions on Information Theory. 31 (4): 469–472. (1985.)
  * [5] J. Katz, Y. Lindell. Introduction to modern cryptography. CRC Press, (2007).
  * [6] Marques JF, Bernardino J. Analysis of Data Anonymization Techniques. InKEOD 2020 (pp. 235-241).
  * [7] S. Martín, P. Morillo, J.L. Villar, Computing the order of points on an elliptic curve modulo N is as difficult as factoring N, Applied Mathematics Letters, Volume 14, Issue 3, Pages 341-346, ISSN 0893-9659 (2001.)
  * [8] Regulation (EU) 2016/679 of the European Parliament and of the Council of 27 April 2016 on the protection of natural persons with regard to the processing of personal data and on the free movement of such data, and repealing Directive 95/46/EC (General Data Protection Regulation). OJ L 119/1 (2016).
  * [9] Ranganathan, A. The levenberg-marquardt algorithm. Tutorial On LM Algorithm. 11, 101-110 (2004)
  * [10] A. Shamir. How to share a secret, Communications of the ACM, (1979.)
  * [11] Brown, S., Rose, J. FPGA and CPLD architectures: A tutorial. IEEE Design & Test Of Computers. 13, 42-57 (1996)
  * [12] Peleg, A., Weiser, U. MMX technology extension to the Intel architecture. IEEE Micro. 16, 42-50 (1996)
  * [13] Chen, L., Moody, D., Randall, K., Regenscheid, A. & Robinson, A. Recommendations for discrete logarithm-based cryptography: Elliptic curve domain parameters. (Lily Chen, Dustin Moody, Karen Randall, Andrew Regenscheid, Angela Robinson, 2023)