remarkRemark hypothesisHypothesis claimClaim Epps Effect and the Signature of
Short-Term Momentum TradersJ. Busca and L. Thomir [][nocite]ex_supplement

# Epps Effect and the Signature of Short-Term Momentum Traders††thanks:
Submitted to the editors 9/11/2023.

Jérôme Busca JB Quantitative Solutions, LLC<EMAIL_ADDRESS>Léon Thomir CentraleSupélec, work completed while a student at ESILV ()
<EMAIL_ADDRESS>
###### Abstract

It is a well-documented fact that the correlation function of the returns on
two ”related” assets is generally increasing as a function of the horizon $h$
of these returns. This phenomenon, termed the Epps Effect, holds true in a
wide variety of markets, and there is a large body of literature devoted to
its theoretical justification. Our focus here is to describe and understand a
deviation to the Epps effect, observed in the context of the foreign exchange
and cryptocurrency markets. Specifically, we document a sharp local maximum of
the cross-correlation function of returns on the Euro EUR/USD and Bitcoin
BTC/USD pairs as a function of $h$. Our claim is that this anomaly reveals the
activity of short-term momentum traders.

###### keywords:

Epps Effect, Momentum, Cryptocurrency, Forex

60G15, 91G60, 60G55, 91B28, 91B70

##  1 Introduction

In 1979, Epps [5] highlighted the presence of a significant drop in the
correlation between stocks when decreasing the time horizon of returns $h$.
This phenomenon was first observed on stocks [11, 3, 2] and then in other
markets such as foreign exchange [7]. Several theoretical justifications have
been proposed to account for the phenomenon. The Epps effect can, for
instance, be explained by a lead-lag phenomenon among specific stocks [9], or
by the asynchronous nature of ticks in liquid markets [8], although Tóth et
al. established in a subsequent paper that tick asynchrony can’t fully
explains the formation of Epps curves [10] (an observation that was confirmed
in a recent study [4]).

Our focus in this paper is the cross-correlation function of the returns on
the EUR/USD (Euro rate in US dollar) and BTC/USD (Bitcoin price) pairs. This
choice was driven by the fact that each pair is the flagship asset in its own
market — forex and crypto-currencies, respectively — and that one can strongly
suspect an interesting interplay between the traditional currency market, and
the relatively new crypto-currency market.

It is widely known that the collective actions of traders play a central role
in the dynamics of financial markets. The interactions and decisions of
individual traders, motivated by a myriad of factors, combine to form a
composite force known as the market factor. We will use this key insight to
model traders’ actions to first order as buying and selling both assets
EUR/USD and BTC/USD at the same time, in the manner of an index (this can also
be seen as ”buying or selling the dollar”).

This paper is designed as follows. In Section 2, we introduce the data set we
use and compute the experimental cross-correlation function $\rho(h)$. In
Section 3, we use a simple Gaussian model to explain the peak we observe in
$\rho$. Finally, in Section 4, we build a more realistic agent-based Monte
Carlo simulation which, once calibrated, shows good agreement with the data.

##  2 Empirical Analysis

In this section, we conduct an empirical analysis on forex exchange and
cryptocurrency markets through the leading pairs EUR/USD and BTC/USD and
highlight the presence of the Epps curve, along with a deviation to it.

###  2.1 Data

To be specific, we actually chose to conduct our analysis on the EUR/USDT and
BTC/USDT pairs. Tether (USDT) is a type of cryptocurrency referred to as a
stablecoin, designed to have a value which is pegged to the US Dollar. Our
choice was guided by the availability of high-quality data sets, as well as
the necessity to avoid data issues such as asynchrony between the FX and
crypto markets, among other factors.

We used a best bid/best offer data set from Binance (the largest crypto
exchange), provided by Tardis, over the period 27 November 2020 – 19 July
2022. As the cryptocurrency market operates continuously, unlike the FX
market, we removed weekends to align Bitcoin data with FX data. Then, we
calculated the simple-mid prices $P(t)$ for further analysis.

###  2.2 Correlations

We define the log-return of length $h$ for an asset of mid price $P(t)$ to be:

$\displaystyle r_{h}(t)=\log\left(\frac{P(t)}{P(t-h)}\right),$

and we denote by $\rho(h)$ the correlations between the log-returns on
EUR/USDT (1) and BTC/USDT (2) as a function of the length $h$:

(1) $\rho(h)=\frac{\left\langle\left(r_{h}^{1}-\langle
r_{h}^{1}\rangle\right)\left(r_{h}^{2}-\langle
r_{h}^{2}\rangle\right)\right\rangle}{\sigma^{1}\sigma^{2}},$

where the bracket $\langle.\rangle$ denotes time average, and the standard
deviations of the returns are defined as usual by:

$\displaystyle\sigma^{k}=\sqrt{\left\langle\left(r^{k}_{h}\right)^{2}\right\rangle-{\langle
r^{k}_{h}\rangle}^{2}},\ \ \ k=1,2.$

We also compute a confidence interval at the 95% level around (1) using a
standard Fisher transformation [6].

###  2.3 Empirical Results

Figure 1: Correlation between EUR/USDT and BTC/USDT returns as a function of
horizon h (1)

Figure 1 represents the cross-correlation $\rho(h)$ between EUR/USDT and
BTC/USDT returns. We observe a classic Epps effect — an overall increase in
correlation with $h$. However, if we zoom in on higher frequencies (Figure 2),
we observe an anomaly which manifests itself as a sharp peak around 60
seconds, with value 0.15, followed by a decrease to 0.12. These fluctuations
are statistically significant.

Figure 2: Correlation between EUR/USDT and BTC/USDT returns as a function of
horizon h (1)

##  3 A Gaussian Model of Momentum

In this section we describe a simple Gaussian model which explains how
momentum traders’ activity can generate a peak in the cross-correlation
function. We assume in the following there are only two assets, with prices
$s_{t}=(s^{1}_{t},s^{2}_{t})$, and that the market consists of: i) a momentum
trader; ii) a noise trader; and iii) a market-maker. The momentum trader is
assumed to trade solely in an equally-weighted index based on $s^{1}$ and
$s^{2}$ and to use a simple momentum rule with window $\tau>0$. For
computational ease, we assume the trader looks at the simple returns (as
opposed to percentage or log returns) of the two assets. With these rules, his
inventory (in both asset 1 and 2) at time $t$ is given by:

(2)
$p_{t}=\bar{p}\left(s^{1}_{t}-s^{1}_{t-\tau}+s^{2}_{t}-s^{2}_{t-\tau}\right),$

where $\bar{p}>0$ is a constant. As for the noise trader, we model his
inventory $X_{t}$ using an Ornstein-Uhlenbeck process with zero mean:

(3) $dX^{k}_{t}=-\lambda X^{k}_{t}dt+\sigma dW^{k}_{t},\ \ k=1,2,$

with $\lambda,\sigma>0$ constant (independent of $k$ for simplicity) and
$W^{1}_{t}$, $W^{2}_{t}$ independent Brownians, and with independent Gaussian
initial conditions equal to the stationary distribution :

(4) $X^{k}_{0}\,{\buildrel d\over{=}}\,{\cal
N}\left(0,\sigma^{2}/2\lambda\right),\ \ k=1,2.$

The assumption that the noise trader’s inventory should follow (3) is fairly
natural and can be shown to hold, for instance, in the classic Avellaneda-
Stoikov model [1], in the limit of large speed of trading (when Poisson
processes reach their diffusion limit).

Lastly, the market-maker is simply the counterparty (liquidity provider) to
the momentum and noise traders. If we further assume a linear market impact
for all liquidity-taking trades, with elasticity $\theta>0$, along with some
noise with volatility $\nu>0$, we can easily write down the price process
increment for the two assets:

(5) $ds^{k}_{t}=\theta\left(dp_{t}+dX^{k}_{t}\right)+\nu dZ^{k}_{t},\ \
k=1,2,$

where $Z^{k}_{t}$ are standard Brownian motions, independent of each other and
of the $W^{k}_{t}$’s. Integrating (5) and choosing initial conditions that
don’t generate extra constants for simplicity, we get:

(6)
$\left\\{\begin{split}p_{t}&=\bar{p}\left(s^{1}_{t}-s^{1}_{t-\tau}+s^{2}_{t}-s^{2}_{t-\tau}\right)\\\
s^{k}_{t}&=\theta\left(p_{t}+X^{k}_{t}\right)+\nu Z^{k}_{t},\ \
k=1,2.\end{split}\right.$

The price process therefore satisfies:

(7)
$s^{k}_{t}=\varepsilon\left(s^{1}_{t}-s^{1}_{t-\tau}+s^{2}_{t}-s^{2}_{t-\tau}\right)+\theta
X^{k}_{t}+\nu Z^{k}_{t},\ \ k=1,2,$

where $\varepsilon=\bar{p}\theta>0$ can be interpreted as a non-dimensional
coupling parameter of the system. In the following, we always assume
$\varepsilon\ll 1$ and expand all relevant quantities to first order in
$\varepsilon$. From (7), the price process can be rewritten

(8) $\left\\{\begin{split}(1-\varepsilon)s^{1}_{t}-\varepsilon
s^{2}_{t}&=-\varepsilon\left(s^{1}_{t-\tau}+s^{2}_{t-\tau}\right)+\theta
X^{1}_{t}+\nu Z^{1}_{t}\\\ \\\ -\varepsilon
s^{1}_{t}+(1-\varepsilon)s^{2}_{t}&=-\varepsilon\left(s^{1}_{t-\tau}+s^{2}_{t-\tau}\right)+\theta
X^{2}_{t}+\nu Z^{2}_{t}.\end{split}\right.$

Noting that $\frac{1-\varepsilon}{1-2\varepsilon}\simeq 1+\varepsilon$, from
(8) we have, to first order in $\varepsilon$

(9)
$\left\\{\begin{split}s^{1}_{t}&\simeq-\varepsilon\left(s^{1}_{t-\tau}+s^{2}_{t-\tau}\right)+\left(1+\varepsilon\right)\left(\theta
X^{1}_{t}+\nu Z^{1}_{t}\right)+\varepsilon\left(\theta X^{2}_{t}+\nu
Z^{2}_{t}\right)\\\ \\\
s^{2}_{t}&\simeq-\varepsilon\left(s^{1}_{t-\tau}+s^{2}_{t-\tau}\right)+\varepsilon\left(\theta
X^{1}_{t}+\nu Z^{1}_{t}\right)+\left(1+\varepsilon\right)\left(\theta
X^{2}_{t}+\nu Z^{2}_{t}\right),\\\ \end{split}\right.$

so that

(10) $\begin{split}s^{k}_{t}&\simeq\theta X^{k}_{t}+\nu
Z^{k}_{t}+\varepsilon\theta\left(X^{1}_{t}-X^{1}_{t-\tau}\right)+\varepsilon\nu\left(Z^{1}_{t}-Z^{1}_{t-\tau}\right)\\\
&+\varepsilon\theta\left(X^{2}_{t}-X^{2}_{t-\tau}\right)+\varepsilon\nu\left(Z^{2}_{t}-Z^{2}_{t-\tau}\right),\
\ k=1,2\\\ \end{split}$

and the $h$-horizon return is given by:

(11)
$\begin{split}s^{k}_{t+h}-s^{k}_{t}&=\theta\left(X^{k}_{t+h}-X^{k}_{t}\right)+\nu\left(Z^{k}_{t+h}-Z^{k}_{t}\right)\\\
&+\varepsilon\theta\left(\left(X^{1}_{t+h}-X^{1}_{t}\right)-\left(X^{1}_{t+h-\tau}-X^{1}_{t-\tau}\right)\right)\\\
&+\varepsilon\theta\left(\left(X^{2}_{t+h}-X^{2}_{t}\right)-\left(X^{2}_{t+h-\tau}-X^{2}_{t-\tau}\right)\right)\\\
&+\varepsilon\nu\left(\left(Z^{1}_{t+h}-Z^{1}_{t}\right)-\left(Z^{1}_{t+h-\tau}-Z^{1}_{t-\tau}\right)\right)\\\
&+\varepsilon\nu\left(\left(Z^{2}_{t+h}-Z^{2}_{t}\right)-\left(Z^{2}_{t+h-\tau}-Z^{2}_{t-\tau}\right)\right)\\\
\end{split}$

The quantities of interest are

(12) $c_{kl}(h):=\langle s^{k}_{t+h}-s^{k}_{t},s^{l}_{t+h}-s^{l}_{t}\rangle,\
\ k,l=1,2,$

as well as the cross-correlation function

(13) $\rho(h):=\frac{c_{12}(h)}{\sqrt{c_{11}(h)c_{22}(h)}}$

We can now state our main result, an explicit formula for the cross-
correlation function (13).

Main Result

######  Theorem 3.1.

To first order in $\varepsilon$, we have:

(14) $\frac{1}{2\varepsilon}\rho(h)=\frac{\theta^{2}\left(1-e^{-\lambda
h}-e^{-\lambda\tau}+\frac{1}{2}e^{-\lambda(h+\tau)}+\frac{1}{2}e^{-\lambda\left|h-\tau\right|}\right)+\xi
h\wedge\tau}{\theta^{2}\left(1-e^{-\lambda h}\right)+\xi h},$

where $h\wedge\tau=\mathrm{min}(h,\tau)$, and $\xi$ is the parameter
$\xi=\frac{\nu^{2}}{\sigma^{2}/\lambda}$

To establish this result, we will need the following lemma.

######  Lemma 3.2.

If $m\in\left\\{0,1\right\\}$, we have, for all $t,h>0$:

(15) $\langle
Z_{t+h}-Z_{t},Z_{t+h-m\tau}-Z_{t-m\tau}\rangle=\left(h-m\tau\right)_{+},$

where $a_{+}=\mathrm{max}(a,0)$; and

(16) $\langle
X_{t+h}-X_{t},X_{t+h-m\tau}-X_{t-m\tau}\rangle=\frac{\sigma^{2}}{2\lambda}\left(2e^{-\lambda
m\tau}-e^{-\lambda(h+m\tau)}-e^{-\lambda\left|h-m\tau\right|}\right),$

where $Z_{t}$ is either $Z^{1}_{t}$ or $Z^{2}_{t}$, and $X_{t}$ is either
$X^{1}_{t}$ or $X^{2}_{t}$.

Proof of Lemma 3.2

$Z$ is a standard Brownian motion. Thus, using ordinary stochastic calculus

$\langle Z_{t},Z_{s}\rangle=t\wedge s\hbox{\ \ for all }s,t>0.$

Therefore

$\begin{split}\langle
Z_{t+h}-Z_{t},Z_{t+h-m\tau}-Z_{t-m\tau}\rangle&=(t+h-m\tau)+(t-m\tau)\\\
&-t\wedge(t+h-m\tau)-(t-m\tau)\wedge(t+h)\\\
&=(t+h-m\tau)+(t-m\tau)+(-t)\vee(-t-h+m\tau)\\\
&+(-t+m\tau)\vee(-t-h),\end{split}$

where $a\vee b=\hbox{max}(a,b)$. Since $\mu+a\vee b=(a+\mu)\vee(b+\mu)$ for
all $\mu$, we find the above is equal to

$(h-m\tau)_{+}+(-h-m\tau)_{+}=(h-m\tau)_{+},$

since $m,h,\tau\geq 0$, which establishes (15).

$X$ is an Ornstein-Uhlenbeck process following (3) and (4). Therefore

$X_{t}=e^{-\lambda t}X_{0}+\sigma\int_{0}^{t}e^{-\lambda(t-u)}dW_{u},$

and for all $s,t>0$ we have

$\begin{split}\langle
X_{t},X_{s}\rangle&=e^{-\lambda(t+s)}\frac{\sigma^{2}}{2\lambda}+\sigma^{2}e^{-\lambda(t+s)}\int_{0}^{t\wedge
s}e^{2\lambda u}du\\\
&=e^{-\lambda(t+s)}\frac{\sigma^{2}}{2\lambda}\left(1+\left(e^{2\lambda
t\wedge s}-1\right)\right)\\\
&=\frac{\sigma^{2}}{2\lambda}e^{-\lambda\left|t-s\right|},\end{split}$

since $t+s-2t\wedge s=\left|t-s\right|$, from which (16) follows easily.

Let’s now prove the main result, Theorem 3.1. By symmetry
$c_{11}(h)=c_{22}(h)$ for all $h>0$ and, using (11) and Lemma 3.2, we have, to
first order in $\varepsilon$:

(17)
$c_{11}(h)=c_{22}(h)=\frac{\sigma^{2}\theta^{2}}{\lambda}\left(1-e^{-\lambda
h}\right)+\nu^{2}h\\\ $

and

(18)
$\frac{1}{2\varepsilon}c_{12}(h)=\frac{\sigma^{2}\theta^{2}}{\lambda}\left(1-e^{-\lambda
h}-e^{-\lambda\tau}+\frac{1}{2}e^{-\lambda(h+\tau)}+\frac{1}{2}e^{-\lambda\left|h-\tau\right|}\right)+\nu^{2}\left(h-(h-\tau)_{+}\right);$

and since $h-(h-\tau)_{+}=h\wedge\tau$, combining (17) and (18) and dividing
by $\sigma^{2}/\lambda$, we proved (14).

Now that we have an explicit formula for the cross-correlation function
$\rho(h)$ with (14), we can plot its typical shape. The plot below illustrates
$\rho(h)$ for the following values of the parameters: $\lambda=0.03162$,
$\xi=0.0001$, $\theta=0.6$, $\tau=66$, $\varepsilon=0.0505$. These values come
from the calibration of the model to our data set (see next section).

Figure 3: Cross-correlation $\rho(h)$ in (14)

We observe on Figure 3 that the cross-correlation $\rho(h)$ in (14) exhibits a
sharp peak (kink) at  $\mathbf{h=\tau}$, the trading horizon of the momentum
trader. We therefore established — in the simplified setting of our Gaussian
model — how the signature of index momentum traders can easily be detected
with the help of the cross-correlation function.

##  4 Agent-based Simulation

In this section, we build an agent-based Monte Carlo simulation which extends
our simplistic Gaussian model by incorporating more realistic features. As in
the Gaussian model, we assume in the following there are only two assets, with
prices $s_{t}=(s^{1}_{t},s^{2}_{t})$, and that the market consists of: i) a
momentum trader; ii) a noise trader; and iii) a market-maker.

In this simulation, we relax the constraint that Poisson processes reach their
diffusion limit. Thus, we use the ”original” form of Avellaneda Stoikov’s
model [1].
On the other hand, to avoid unbounded inventories, we impose a cap on the
momentum trader’s positions. For the sake of realism, we also take into
account the tick sizes $\eta$.

###  4.1 Noise trader

We model the noise trader following [1]. Specifically, we assume trade
executions occur at the ask (resp. bid) price at a rate described by Poisson
processes $N_{t}^{a}$ and $N_{t}^{b}$ with respective intensity parameters
$\lambda^{a,b}_{t}(\delta^{a,b}_{t})=A^{a,b}e^{-k^{a,b}\delta^{a,b}_{t}}$,
where $\delta^{a,b}_{t}$ is the half-spread between the mid and the ask (resp.
bid) price, and $A^{a,b},k^{a,b}>0$ are constants measuring the liquidity of
the market.

The inventory $q_{t}^{n}$ is then simply given by

$\displaystyle q_{t}^{n}=\psi^{n}(N_{t}^{a}-N_{t}^{b}),$

where $\psi^{n}>0$ is a parameter (trade size).

###  4.2 Momentum trader

We assume the momentum trader buys and sells an equally-weighted index based
on $s^{1}$ and $s^{2}$, denoted by $\mathrm{index}_{t}$, and uses a simple
momentum rule with window $\tau>0$. To manage risk, the momentum trader has a
maximum position constraint. We also force a long (resp. short) position when
the index moves back above (resp. below) the moving average.

Their inventory $q_{t}^{m}$ is modeled as

$q_{t}^{m}=\left\\{\begin{split}&\mathrm{min}(\mathrm{max}(q_{t-1}^{m},0)+\psi^{m},\
q_{\mathrm{max}}^{m})\ \ \ \ \hbox{if}\ \ \mathrm{index}_{t-1}>\mu_{t-1}\\\
\\\ &\mathrm{max}(\mathrm{min}(q_{t-1}^{m},0)-\psi^{m},\
-q_{\mathrm{max}}^{m})\ \ \hbox{otherwise,}\end{split}\right.$

where $\mu_{t}$ is the $\tau-$moving average of $\mathrm{index}_{t}$,
$\psi^{m}>0$ is the trade size, and $q_{\mathrm{max}}^{m}>0$, is the maximum
absolute inventory.

###  4.3 Market maker

The market-maker is the liquidity provider. It is therefore modeled as the
counterparty to the noise and momentum traders. His inventory is simply
$q_{t}^{mm}=-q_{t}^{n}-q_{t}^{m}$. He adjusts the half-spread
$\delta^{a,b}_{t}$ as a function of k, q, the tick size $\eta$, the
volatility, and a risk aversion parameter $\gamma$ as described in [1].

Furthermore, the impact of liquidity-taking orders on the mid price is assumed
to be linear with elasticity $\theta>0$, on top of a Brownian noise with
volatility $\nu>0$. Before rounding to tick size $\eta$, the mid price
processes are therefore given by

$\displaystyle
s_{t}^{k}=s_{t-1}^{k}+\theta^{k}(q_{t}^{n}+q_{t}^{m}-q_{t-1}^{n}-q_{t-1}^{m})+\nu^{k}dZ^{k}_{t},\
\ k=1,2,$

where $dZ^{k}_{t}$ are Brownian increments.

###  4.4 Monte Carlo Simulation

We ran a discrete Monte-Carlo simulation of the model above, with time
resolution $\mathrm{dt}=0.5s$, over a duration $T=2\cdot 10^{7}\mathrm{dt}$,
i.e. around 115 trading days. We took realistic initial values for $s^{1}$ and
$s^{2}$, respectively $1.10$ USD and $30,000$ USD.

We then calibrated the parameters of our model to our data set. Here are the
fitted values: $\eta=[1e-4,1e-2]$, $A=[1,1]$, $\theta=[2.7e-11,2.7e-6]$,
$\nu=[3e-5,2.04]$, $\gamma=[6.46e-6,3.47e-9],k=[3466,34.66]$,
$\psi^{n}=[100000,4]$, $\psi^{m}=[3000000,120]$, $\tau=500s$,
$q_{\mathrm{rm}}^{m}=[6500000,250]$.

Figure 4 shows the empirical correlation function (blue), the agent-based
simulated correlation function (orange) and the correlation function of our
simple Gaussian model (green) as a function of horizon h, expressed in
seconds.

Figure 4: Comparison of the cross-correlation functions $\rho(h)$

We observe that the cross-correlation function of our agent-based simulation
exhibits a good fit to the data up to h = 120 seconds, before dropping and
failing to maintain a correlation level of 0.12. We suspect that the impact of
agents with longer trading horizons is the cause of a persistent correlation
for larger values of $h$.

As to the Gaussian model, we see a sharp peak around h = $\tau=66$ seconds,
followed by a slow drop to the 0.08 level. Due to its limitations, it is
impossible to create a ”smoother” or larger peak in correlation.

The main difference in the choice of parameters between the Gaussian and
agent-based models is the value $h=\tau$, where the peak of $\rho$ is reached
in the Gaussian model. Indeed, for the simulated model, it is necessary to set
a higher value for $\tau$ because of the switch to a discrete model and the
resulting effect on the noise trader’s orders. We made available the Python
code for the simulation at

https://github.com/RimohtL/EppsEffect

## References

  * [1] M. Avellaneda and S. Stoikov, High-frequency trading in a limit order book, Quantitative Finance, 8 (2008), pp. 217–224, https://doi.org/10.1080/14697680701381228.
  * [2] G. Bonanno, F. Lillo, and R. N. Mantegna, High-frequency cross-correlation in a set of stocks, Quantitative Finance 96-104, 1 (2001), pp. 96–104, https://doi.org/10.1080/713665554.
  * [3] K. Chan, Imperfect information and cross-autocorrelation among stock prices, The Journal of Finance, 48 (1993), pp. 1211–1230, https://doi.org/https://doi.org/10.1111/j.1540-6261.1993.tb04752.x.
  * [4] P. Chang, E. Pienaar, and T. Gebbie, The $\mathrm{E}$pps effect under alternative sampling schemes, Physica A: Statistical Mechanics and its Applications 126329, 583 (2021), https://doi.org/10.1016/j.physa.2021.126329.
  * [5] T. W. Epps, Comovements in stock prices in the very short run, Journal of the American Statistical Association, 74 (1979), pp. 291–298, https://doi.org/10.2307/2286325.
  * [6] R. A. Fisher, On the ”probable error” of a coefficient of correlation deduced from a small sample, Metron, 1 (1921), p. 3–32.
  * [7] M. Lundin, M. Dacorogna, and U. Müller, Correlation of high frequency financial time series, The Financial Markets Tick by Tick, Pierre Lequeux, ed., (1999), pp. 91–126.
  * [8] R. Renò, A closer look at the epps effect, International Journal of Theoretical and Applied Finance 6(1), 87-102, 6 (2003), pp. 87–102, https://doi.org/10.1142/S0219024903001839.
  * [9] B. Tóth and J. Kertész, Increasing market efficiency: Evolution of cross-correlations of stock returns, Physica A: Statistical Mechanics and its Applications 360 505-515, 360 (2006), pp. 505–515, https://doi.org/10.1016/j.physa.2005.06.058.
  * [10] B. Tóth and J. Kertész, The epps effect revisited, Quantitative Finance, 9 (2009), pp. 793–802, https://doi.org/10.1080/14697680802595668.
  * [11] A. A. Zebedee and M. Kasch-Haroutounian, A closer look at co-movements among stock returns, Journal of Economics and Business, 61 (2007), pp. 279–294, https://doi.org/10.1016/j.jeconbus.2008.11.001.


[a] Willem E. A. Verplanke

# Pseudoscalar transition form factors and the hadronic light-by-light
contribution to the muon $g-2$

Antoine Gérardin  Jana N. Guenther  Lukas Varnhorst  and  for the Budapest-
Marseille-Wuppertal Collaboration

###### Abstract

We present preliminary results from our calculation of the pseudoscalar
transition form factors of the $\eta$ and $\eta^{\prime}$ mesons using
staggered quarks on $N_{f}=2+1+1$ gauge ensembles generated by the Budapest-
Marseille-Wuppertal collaboration. These transition form factors are an
important input for the hadronic light-by-light contribution to the muon
$(g-2)$. We first elaborate on the extraction of the masses of the $\eta$ and
$\eta^{\prime}$ mesons, that mix under the dynamics of QCD. Thereafter, we
show our preliminary results for the pseudoscalar transition form factors,
focusing on the $\eta$ meson in the absence of mixing.

##  1 Introduction

Understanding the discrepancy between the theory calculation and the
experimental value of the muon anomalous magnetic moment ($a_{\mu}$) [1] is
one of the most pressing challenges in present-day particle physics. The
theoretical error is dominated by the hadronic vacuum polarization (HVP) and
the hadronic light-by-light (HLbL) contributions. A recent computation of the
leading order HVP by the Budapest-Marseille-Wuppertal (BMW) collaboration has
however shed new light on the issue [2]. In particular, it reduces the tension
with the experimental value of $a_{\mu}$.

The HLbL interaction is shown in Figure 1.

Figure 1: Hadronic light-by-light diagram and its decomposition into the
dominant pseudoscalar poles. A wobbly line indicates a photon, the straight
line the muon and a blob the non-perturbative hadronic interactions encoded in
the pseudoscalar TFF.

This contribution results from the scattering of four photons by means of a
non-perturbartive hadronic interaction. Two model-independent routes to
calculate the diagram are nowadays applied, namely the dispersive [3, 4, 5]
and lattice [6, 7] computations. Crucial input for the dispersive approach are
the transition form factors (TFFs) of the $\pi^{0}$, $\eta$ and
$\eta^{\prime}$ mesons, which are directly linked to pseudoscalar pole
contributions to the $a_{\mu}^{\mathrm{HLbL}}$. The TFF describes the
interaction of one on-shell pseudoscalar meson with two off-shell photons. The
pion TFF has already been calculated [8, 9], and we provide here our first
results towards a calculation of the $\eta$ and $\eta^{\prime}$ TFFs. Note
finally that, in experiment, the $\eta$ and $\eta^{\prime}$ TFFs are difficult
to access in the doubly virtual regime, where the photons have the same
virtuality, so a theoretical calculation is of particular importance. ETM has
already presented preliminary results of a calculation of the $\eta$,
$\eta^{\prime}$ TFFs [10]. The mixing between the two states encumbers the
calculation, but only a 20% precision on the contribution to
$a_{\mu}^{\mathrm{HLbL}}$ is necessary to match future experimental precision
[11].

In line with this interest, the second aim of this project is a calculation of
the $\eta$ and $\eta^{\prime}$ masses and mixing angles, using staggered
quarks. Recent calculations of these quantities have been done using Wilson-
Clover quarks [12] and twisted-mass quarks [13], and have shown good
correspondence with experiment. The $\eta$ and $\eta^{\prime}$ mesons are
interesting particles to simulate, as the large mass of the $\eta^{\prime}$
meson is directly related to the non-perturbative dynamics of QCD and the
anomaly in the $U(1)_{A}$ charge symmetry of the QCD Lagrangian [14]. Besides,
flavor singlet mesons require the use of disconnected correlators, which are
more challenging to compute and couple directly to the sea. The fourth-root
trick employed for the sea quarks in staggered lattice QCD simulations has
also been questioned [15, 16]. The $\eta^{\prime}$ meson is hence of interest,
and can therefore provide another hint towards the ability of staggered quarks
to reproduce the axial anomaly in the continuum limit.

##  2 Simulation details

We use $N_{f}=2+1+1$ dynamical staggered fermions with four steps of stout
smearing generated by the BMW collaboration [2]. These gauge ensembles are at
nearly physical pion and kaon mass. We plan to exploit five different lattice
spacings in the available range and consider $L=3$, $4$ and $6$ fm boxes for
finite-size effect studies. Simulations are performed in the isospin limit
where $m_{u}=m_{d}\equiv m_{\ell}$. In Table 1 we summarize the details of the
two ensembles used in this preliminary study. The three-point functions have
only been computed on the coarse ensemble.

|  $\beta$ |  a[fm] |  $L/a\times T/a$ |  # conf
---|---|---|---|---
three-point |  3.7000 |  0.1315 |  32 $\times$ 64 |  900
two-point |  3.8400 |  0.0952 |  32 $\times$ 64 |  1100
Table 1: Summary of two ensembles with lattice spacing, lattice size and
number of gauge configurations.

##  3 $\eta$ and $\eta^{\prime}$ mass

###  3.1 Mass extraction

The analysis of the $\eta-\eta^{\prime}$ system starts by the quark model. One
introduces an SU(3) octet and singlet field

$\displaystyle\eta_{8}(x)$
$\displaystyle=\frac{1}{\sqrt{6}}\left(\overline{u}\gamma_{5}u(x)+\overline{d}\gamma_{5}d(x)-2\overline{s}\gamma_{5}s(x)\right),$
(1) $\displaystyle\eta_{0}(x)$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\overline{u}\gamma_{5}u(x)+\overline{d}\gamma_{5}d(x)+\overline{s}\gamma_{5}s(x)\right).$
(2)

The $\eta_{8}(x)$ field is a Goldstone realization of spontaneously broken
symmetry group $G_{F}=U(1)_{V}\otimes SU(3)_{L}\otimes
SU(3)_{R}\overset{SSB}{\longrightarrow}U(1)_{V}\otimes SU(3)_{V}$ of the
classical QCD Lagrangian. On the other hand, $\eta_{0}$ is the would-be
Goldstone boson realization of the $U(1)_{A}$ charge symmetry, that is broken
at the quantum level. Both fields have quantum numbers $J^{PC}=0^{-+}$, like
the $\pi^{0}$. However, no mixing with the pion occurs in our lattice
simulations since we work in the isospin limit.

We consider the matrix of correlators built from the lattice operators
$O_{8}$, $O_{0}$ that interpolate the SU(3) fields

$\displaystyle\mathcal{C}(t)$ $\displaystyle\equiv\begin{pmatrix}\langle
O_{8}(t)O_{8}^{\dagger}(0)\rangle&\langle O_{8}(t)O_{0}^{\dagger}(0)\rangle\\\
\langle O_{0}(t)O_{8}^{\dagger}(0)\rangle&\langle
O_{0}(t)O_{0}^{\dagger}(0)\rangle\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\frac{1}{3}\left(C_{\ell}+2C_{s}-4D_{\ell
s}+2D_{\ell\ell}+2D_{ss}\right)&\frac{\sqrt{2}}{3}\left(C_{\ell}-D_{\ell
s}-D_{ss}-C_{s}+2D_{\ell\ell}\right)\\\
\frac{\sqrt{2}}{3}\left(C_{\ell}-D_{\ell
s}-D_{ss}-C_{s}+2D_{\ell\ell}\right)&\frac{1}{3}\left(2C_{\ell}+C_{s}+4D_{\ell\ell}+4D_{\ell
s}+D_{ss}\right)\end{pmatrix},$ (3)

where $C_{q}$ and $D_{qq^{\prime}}$ are respectively the connected and
disconnected correlator for quark flavor $q,q^{\prime}=\ell,s$ and we have
made use of $D_{ls}=D_{sl}$. Notice that the non-vanishing off-diagonal terms
explicit the mixing between $\eta_{8}$ and $\eta_{0}$; furthermore, even in
the $SU(3)_{F}$ limit, the disconnected contributions to $\langle
O_{0}(t)O_{0}^{\dagger}(0)\rangle$ do not vanish. Lastly, in this chiral
limit, the physical $\eta$ and $\eta^{\prime}$ state correspond exactly to the
$\eta_{8}$ and $\eta_{0}$ state.

Masses of the $\eta$ and $\eta^{\prime}$ mesons can then be obtained by
solving a Generalized Eigenvalue Problem (GEVP) [17]

$\mathcal{C}(t)v_{n}(t,t_{0})=\lambda_{n}(t,t_{0})\mathcal{C}(t_{0})v_{n}(t,t_{0}).$
(4)

Here $\mathcal{C}(t)$ is an $N\times N$ matrix, $\lambda_{n}$ are the
eigenvalues, $v_{n}$ are the eigenvectors, $t_{0}/a$ is a free parameter that
is fixed to $1$ in this analysis, and $n\in\\{1,\dots,N\\}$. From the
eigenvalues one can define an effective mass that converges to the meson mass
at large time

$E_{n}^{\textrm{eff}}(t)=\log\left(\frac{\lambda_{n}(t,t_{0})}{\lambda_{n}(t+1,t_{0})}\right)\overset{t\to\infty}{\longrightarrow}m_{n},$
(5)

in the region where the backward propagating quarks can be neglected. In the
next section we present the lattice operators that interpolate the $\eta$,
$\eta^{\prime}$ states in the staggered quark formalism.

###  3.2 Staggered mesonic operators

Staggered mesonic operators have been classified by Golterman in [18].
Following the notation of [19], there are two taste-singlet operators that
couple directly to the $\eta$, $\eta^{\prime}$ mesons

$\displaystyle\mathcal{O}_{3}(x)$
$\displaystyle=\frac{1}{6}\sum_{ijk}\epsilon_{ijk}\overline{\chi}(x)\left[\eta_{i}\Delta_{i}\left[\eta_{j}\Delta_{j}\left[\eta_{k}\Delta_{k}\right]\right]\right]\chi(x)\equiv\overline{\chi}(x)\hat{O}_{3}\chi(x),$
(6) $\displaystyle\mathcal{O}_{4}(x)$
$\displaystyle=\frac{1}{2}\eta_{4}(x)\left[\overline{\chi}(x)\hat{O}_{3}\chi_{+}(x)+\overline{\chi}_{+}(x)\hat{O}_{3}\chi(x)\right].$
(7)

Here $\chi(x)$ is a fermionic field and $\epsilon_{ijk}$ is the Levi-Civita
tensor; the symmetric shift $\Delta_{\mu}$, $\chi_{+}(x)$ and the staggered
phase factors $\eta_{\mu}(n)$ in the convention $n=(x,y,z,t)$ are

$\displaystyle\Delta_{\mu}\chi(x)$
$\displaystyle=\frac{1}{2}\left[U_{\mu}(x)\chi(x+\hat{\mu})+U_{\mu}^{\dagger}(x-\hat{\mu})\chi(x-\hat{\mu})\right],$
(8) $\displaystyle\chi_{+}(x)$
$\displaystyle=U_{4}(x)\chi(x+\hat{t}),\quad\overline{\chi}_{+}(x)=\overline{\chi}(x+\hat{t})U_{4}^{\dagger}(x),$
(9) $\displaystyle\eta_{\mu}(n)$ $\displaystyle=(-1)^{\sum_{\nu<\mu}n_{\mu}}.$
(10)

In the ubiquitous staggered operator notation $\Gamma_{D}\otimes\Gamma_{F}$,
where $\Gamma_{D}$ denotes the Dirac structure, and $\Gamma_{F}$ the taste
structure, the 3-link operator $\mathcal{O}_{3}(x)$ couples to an axial vector
$\gamma_{4}\gamma_{5}\otimes 1$ and a parity partner scalar meson
$1\otimes\gamma_{4}\gamma_{5}$. The parity partner for the 3-link operator
creates an oscillatory behavior in the effective energy for the taste-singlet
pion. This makes a trustworthy extraction of the pseudoscalar ground state
energies cumbersome. The 4-link operator $\mathcal{O}_{4}(x)$, which is non-
local in time, couples to a pseudoscalar $\gamma_{5}\otimes 1$ and an exotic
state $\gamma_{4}\otimes\gamma_{4}\gamma_{5}$. As noted in [19, 20], the
coupling to the parity partner state for such an operator is strongly
suppressed. Hence, it is our choice of mesonic operator throughout the
analysis, also because the separation by an even amount of links allows for a
computational trick, as is discussed next.

###  3.3 Pseudoscalar two-point functions

We apply several noise-reduction tricks to compute the correlation functions
effectively. Firstly, we use low-mode averaging (LMA) [21, 22] with $n=300$
modes in our $4$ fm box (and scale modes with the volume of the box).
Secondly, we apply all-mode averaging (AMA) [23, 24] for the stochastic part
of our estimator. And lastly, we employ a Venkataraman-Kilcup reduction trick
[25] for our one-point functions. In this context, it uses the staggered Dirac
operator that connects only even and odd sites, together with the 4-link
operator which only couples even (odd) and even (odd) sites. This allows one
to construct a different pseudoscalar loop estimator that has a decreased
variance. As a consequence of all these tricks, we have managed to reach the
gauge noise for the two-point disconnected correlation functions. In Figure 2,
we show the result for all our correlation functions. Note that we find a good
signal for $C_{\ell}$, $C_{s}$ and $D_{\ell\ell}$, while the signal for the
$D_{\ell s}$ and especially the $D_{ss}$ are lost in the noise quickly. This
is however not an issue, since the effect of the $D_{ss}$ correlator on the
GEVP is extremely small, and to a lesser extent this holds for the $D_{\ell
s}$ as well.

Figure 2: Different connected $C_{q}$ and disconnected $D_{qq^{\prime}}$
correlation functions contributing to the $\eta$, $\eta^{\prime}$
spectroscopic analysis.

###  3.4 Analysis techniques

Extracting the mass of the $\eta$ and $\eta^{\prime}$ mesons is hindered by
the exponential growth of the noise/signal ratio. In fact, by only using a
GEVP, the signal of the two mesons is lost before reaching a plateau. Thus, to
to circumvent this issue, we apply two analysis techniques. Firstly, we employ
a fitting trick to remove excited state contributions in the connected
correlator, which was first introduced in this context by ETM [26]. The
underlying assumption is that most excited state contamination is found in the
connected correlation functions. Since these have a relatively good signal
over noise ratio (however not constant, as we work with taste-singlet
quantities), it is possible to reliably determine the ground state of these
correlators. Connected mesonic correlation functions calculated on a lattice
with periodic boundary conditions, assume a spectral decomposition of the form

$C_{q}(t)=A_{q}\left(\exp(-E_{q}^{(1)}t)+\exp(-E_{q}^{(1)}(T-t))\right)+\mathcal{O}(\exp(-E_{q}^{(2)}t)),\quad
q=\ell,s.$ (11)

Fitting $C_{q}$ to this function in a region where excited states are
negligible, one can safely determine the energy $E_{q}$ and overlap $A_{q}$ of
the ground state. Then one replaces the connected contribution in equation (3)
by this one-exponential fit on the whole time range, leaving the disconnected
contributions unchanged. If excited state contributions in the disconnected
correlators are indeed small, excited states should effectively be removed and
a plateau for the effective mass should be reached early in time. Results of
applying this method to our data are shown in Figure 3a. We cannot resolve
excited states at this stage, but errors are still large, especially for the
$\eta^{\prime}$ effective mass.

The second trick was first introduced in [13], and results from an analytic
understanding of correlation functions in finite volume [27]. In fact, for
$\vec{p}=\vec{0}$, disconnected correlation functions in finite volume, do not
tend to zero at large time if the topological charge $Q$ isn’t correctly
sampled. One finds that

$\displaystyle
D_{qq}(t)\overset{t\to\infty}{\sim}\frac{a^{5}}{T}\left(\chi_{t}-\frac{Q^{2}}{V}+\frac{c_{4}}{2V\chi_{t}}\right).$
(12)

Here, $\chi_{t}$ is the topological susceptibility and $c_{4}$ is the kurtosis
of the topological charge. Such a constant shift can be removed by taking a
discrete derivative of the matrix of correlators $\mathcal{C}(t)$,

$\displaystyle\mathcal{C}(t)\rightarrow\mathcal{C}^{\prime}(t)\equiv\mathcal{C}(t)-\mathcal{C}(t+\Delta
t).$ (13)

This shifted correlation matrix still satisfies a GEVP with $\Delta t/a$ as a
free parameter. Aside of removing this constant shift, it also helps reducing
correlations between time-slices and hence improves the point error. In Figure
3b we show the result of applying this method to our data, alongside with the
first trick, and see that now a plateau is found for $\eta$ and
$\eta^{\prime}$ and errors are reduced drastically. Further, in Figure 3c we
show the result for our correlation functions which carry one unit of
momentum. The errors are slightly smaller and the data are compatible with
$\vec{p}=\vec{0}$, motivating a use of both kinematic frames in future
analyses. Lastly, we observe a mild discretization effect on the
$\eta^{\prime}$ meson, while it’s of the order of 10% for the $\eta$ meson.

(a)

(b)

(c)

Figure 3: Effective mass plots for three scenarios. Left: after applying
method $1$ of section 3.4; Middle: Additionally applied method (2); Right:
Added the $|\vec{p}|=1\cdot\left(\frac{2\pi}{L}\right)$ result. Dashed lines
indicate the experimental values of the $\eta$ and $\eta^{\prime}$ masses.

##  4 Pseudoscalar transition form factors

(a)

(b)

Figure 4: Left: The integrand of the $\eta_{8}$ TFF for specified kinematic
condition; $t_{sep}$ refers to the seperation between the pseudoscalar and one
of the vector currents that is kept fixed. Right: Illustration of the
different topologies of the connected and disconnected contributions to the
three-point function.

In Minkowski space-time, the TFF for a pseudoscalar meson
$\mathcal{F}_{\mathrm{p}\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})$ is defined
by the following matrix elements $M_{\mu\nu}$

$M_{\mu\nu}(p,q_{1})=i\int d^{4}x\,e^{iq_{1}\cdot
x}\bra{\Omega}T\\{J_{\mu}(x)J_{\nu}(0)\\}\ket{P(p)}=\epsilon_{\mu\nu\alpha\beta}\,q_{1}^{\alpha}q_{2}^{\beta}\mathcal{F}_{\mathrm{p}\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2}),$
(14)

where $q_{1}$ and $q_{2}$ are the photon $4$-momenta, $J_{\mu}$ is the
hadronic component of the electromagnetic current and
$\epsilon_{\mu\nu\alpha\beta}$ is a $4$-rank Levi Civita tensor. On the
lattice, one can relate three-point functions to the Euclidean version of
these matrix elements [28]. We write

$M_{\mu\nu}^{E}=\frac{2E_{P}}{Z_{P}}\int_{-\infty}^{\infty}d\tau\,e^{\omega_{1}\tau}\tilde{A}_{\mu\nu}(\tau),$

where $E_{P}$ and $Z_{P}$ are respectively the pseudoscalar energy and overlap
of the pseudoscalar state with the interpolating operator, $\omega_{1}$ is a
free parameter $q_{1}=(\omega_{1},\vec{q}_{1})$,
$q_{2}=(E_{P}-\omega_{1},\vec{q}_{2})$, and

$\displaystyle\tilde{A}_{\mu\nu}(\tau)$
$\displaystyle\equiv\lim\limits_{t_{P}\to\infty}e^{E_{P}(t_{f}-t_{0})}C_{\mu\nu}^{(3)}(\tau,t_{P}),$
$\displaystyle C_{\mu\nu}^{(3)}(\tau,t_{P})$
$\displaystyle=a^{6}\sum_{\vec{x},\vec{z}}\langle
J_{\mu}(\vec{z},t_{i})J_{\nu}(\vec{0},t_{f})P^{\dagger}(\vec{x},t_{0})\rangle
e^{i\vec{p}\cdot\vec{x}}e^{-i\vec{q}_{1}\cdot\vec{z}}.$

Here we have defined $\tau=t_{i}-t_{f}$, the time-difference between the two
vector currents and $t_{P}=\mathrm{min}(t_{f}-t_{0},t_{i}-t_{0})$, the minimal
time separation between the pseudoscalar and the vector currents. In the
three-point function, $P$ is the interpolating operator for the pseudoscalar
meson, expressed in terms of $O_{8}$, $O_{0}$ and appropriate mixing
parameters; $J_{\mu}$ is implemented using the conserved vector current which
does not require renormalization [2]. Note also that, an accurate
determination of the mixing parameters and the masses of the $\eta$,
$\eta^{\prime}$ mesons is important, as they enter in the three-point
function.

$\tilde{A}_{\mu\nu}$ can be decomposed into two scalar functions
$\tilde{A}^{(1)}$, $\tilde{A}^{(2)}$ that form the integrand of the TFF (for
technical details we refer to [9]). In Figure 4, we show preliminary results
for one of these scalar functions in the $\vec{p}=\vec{0}$ frame, where we
focus on $\eta_{8}$, thus discarding any mixing with $\eta_{0}$. The latter is
important, and will be implemented in the future. Note that four different
topologies contribute to the integrand (the hyphen indicates a disconnection):
the two dominant PVV connected and P-VV disconnected loops; and the negligible
PV-V and P-V-V disconnected loops. The P-VV loop has a significant and
negative contribution and therefore we observe a large cancellation with the
PVV loop. As a consequence, a good control of the disconnected VV and P loops
is important.

##  5 Conclusions

We have shown preliminary results toward the extraction of the
$\eta,\eta^{\prime}$ mass using staggered quarks. Using the GEVP and several
noise reduction techniques, we are able to obtain a good signal for the $\eta$
and $\eta^{\prime}$ effective mass. Alongside this, we have shown promising
preliminary results for the integrand of the TFF of the $\eta_{8}$ state. In
this context we have found that the two leading contributions are the fully
connected PVV and disconnected P-VV loops, while the disconnected PV-V and
P-V-V loops are of a marginal size.

Presently, we are accumulating statistics for our P-VV contribution, to
further reduce its error. At the same time, we will explore the signal of the
$\eta$, $\eta^{\prime}$ and their TFFs, including therefore mixing between the
states. We plan to do this on several lattice spacings and volumes, improving
also the statistics on the different ensembles that we have already analyzed.
Ultimately, we will perform a continuum extrapolation of the TFFs.

##  6 Acknowledgements

This publication received funding from the Excellence Initiative of Aix-
Marseille University - A*MIDEX, a French “Investissements d’Avenir” programme,
AMX-18-ACE-005 and from the French National Research Agency under the contract
ANR-20-CE31-0016. Center de Calcul Intensif d’Aix-Marseille is acknowledged
for granting access to its high performance computing resources.

## References

  * [1] T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , _Phys. Rept._ 887 (2020) 1 [2006.04822].
  * [2] S. Borsanyi et al., _Leading hadronic contribution to the muon magnetic moment from lattice QCD_ , _Nature_ 593 (2021) 51 [2002.12347].
  * [3] G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, _Dispersive approach to hadronic light-by-light scattering_ , _JHEP_ 09 (2014) 091 [1402.7081].
  * [4] G. Colangelo, M. Hoferichter, B. Kubis, M. Procura and P. Stoffer, _Towards a data-driven analysis of hadronic light-by-light scattering_ , _Phys. Lett. B_ 738 (2014) 6 [1408.2517].
  * [5] G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, _Dispersion relation for hadronic light-by-light scattering: theoretical foundations_ , _JHEP_ 09 (2015) 074 [1506.01386].
  * [6] E.-H. Chao, R.J. Hudspith, A. Gérardin, J.R. Green, H.B. Meyer and K. Ottnad, _Hadronic light-by-light contribution to $(g-2)_{\mu}$ from lattice QCD: a complete calculation_, _Eur. Phys. J. C_ 81 (2021) 651 [2104.02632].
  * [7] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung et al., _Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment from lattice QCD_ , _Physical Review Letters_ 124 (2020) .
  * [8] A. Gérardin, H.B. Meyer and A. Nyffeler, _Lattice calculation of the pion transition form factor $\pi^{0}\to\gamma^{*}\gamma^{*}$_, _Phys. Rev. D_ 94 (2016) 074507 [1607.08174].
  * [9] A. Gérardin, H.B. Meyer and A. Nyffeler, _Lattice calculation of the pion transition form factor with $N_{f}=2+1$ Wilson quarks_, _Phys. Rev. D_ 100 (2019) 034520 [1903.09471].
  * [10] S. Burri, _Pseudoscalar-Photon Transition Form Factors from Twisted Mass Lattice QCD_ , _Asia-Pacific Symposium for Lattice Field Theory (APLAT)_ (2020) .
  * [11] A. Gérardin, _The anomalous magnetic moment of the muon: status of Lattice QCD calculations_ , _Eur. Phys. J. A_ 57 (2021) 116 [2012.03931].
  * [12] G.S. Bali, V. Braun, S. Collins, A. Schäfer and J. Simeth, _Masses and decay constants of the $\eta$ and $\eta^{\prime}$ mesons from lattice QCD_, 2021\.
  * [13] ETM collaboration, _Flavor-singlet meson decay constants from $N_{f}=2+1+1$ twisted mass lattice QCD_, _Phys. Rev. D_ 97 (2018) 054508 [1710.07986].
  * [14] G. ’t Hooft, _Computation of the Quantum Effects Due to a Four-Dimensional Pseudoparticle_ , _Phys. Rev. D_ 14 (1976) 3432.
  * [15] M. Creutz, _Why rooting fails_ , _PoS_ LATTICE2007 (2007) 007 [0708.1295].
  * [16] A.S. Kronfeld, _Lattice Gauge Theory with Staggered Fermions: How, Where, and Why (Not)_ , _PoS_ LATTICE2007 (2007) 016 [0711.0699].
  * [17] B. Blossier, M.D. Morte, G.v. Hippel, T. Mendes and R. Sommer, _On the generalized eigenvalue method for energies and matrix elements in lattice field theory_ , _Journal of High Energy Physics_ 2009 (2009) 094–094.
  * [18] M.F.L. Golterman, _STAGGERED MESONS_ , _Nucl. Phys. B_ 273 (1986) 663.
  * [19] MT(c) collaboration, _The Hadron spectrum in QCD with dynamical staggered fermions_ , _Nucl. Phys. B_ 389 (1993) 445.
  * [20] G.W. Kilcup and S.R. Sharpe, _A Tool Kit for Staggered Fermions_ , _Nucl. Phys. B_ 283 (1987) 493.
  * [21] L. Giusti, P. Hernandez, M. Laine, P. Weisz and H. Wittig, _Low-energy couplings of QCD from current correlators near the chiral limit_ , _Journal of High Energy Physics_ 2004 (2004) 013–013.
  * [22] T. DeGrand and S. Schaefer, _Improving meson two-point functions by low-mode averaging_ , _Nuclear Physics B - Proceedings Supplements_ 140 (2005) 296–298.
  * [23] G.S. Bali, S. Collins and A. Schäfer, _Effective noise reduction techniques for disconnected loops in lattice qcd_ , _Computer Physics Communications_ 181 (2010) 1570–1583.
  * [24] T. Blum, T. Izubuchi and E. Shintani, _New class of variance-reduction techniques using lattice symmetries_ , _Physical Review D_ 88 (2013) .
  * [25] L. Venkataraman and G. Kilcup, _The eta-prime meson with staggered fermions_ , hep-lat/9711006.
  * [26] C. Michael, K. Ottnad and C. Urbach, _$\eta$ and $\eta$’ mixing from lattice QCD_, _Physical Review Letters_ 111 (2013) .
  * [27] S. Aoki, H. Fukaya, S. Hashimoto and T. Onogi, _Finite volume QCD at fixed topological charge_ , _Physical Review D_ 76 (2007) .
  * [28] X. Ji and C. Jung, _Studying hadronic structure of the photon in lattice QCD_ , _Physical Review Letters_ 86 (2001) 208–211.