[a]Alejandro Miranda

# Isospin-breaking corrections to $\tau^{-}\to\pi^{-}\pi^{0}\nu_{\tau}$ decays
and the muon $g-2$

###### Abstract

We review the isospin-breaking and electromagnetic corrections to the
$\tau^{-}\to\pi^{-}\pi^{0}\nu_{\tau}$ decays, which are used as an input to
the two-pion contributions to the hadronic vacuum polarization (at LO) of the
anomalous magnetic moment ($a_{\mu}$). We extend previous analyses by
Cirigliano et al. working with ChPT with resonances. As an outcome, we improve
the agreement between this determination and the other based on $e^{+}e^{-}$
data. The new results are in better agreement with an old estimation that uses
Vector Dominance Model (VMD), and the discrepancy between the SM prediction
and the combined results from BNL and FNAL is reduced to $2.1\,\sigma$ at
$\mathcal{O}(p^{4})$ and $2.3\sigma$ at $\mathcal{O}(p^{6})$.

##  1 Introduction

The Standard Model of particle physics [1, 2, 3] is one of the most successful
theories that describes strong, weak and electromagnetic interactions. One of
the most precisely measured quantities in particle physics is the muon
anomalous magnetic moment ($a_{\mu}$). A long-standing discrepancy between
theory and experiment about $3-4$ standard deviations has been observed.

The agreement between the latest measurement performed at Fermilab National
Accelerator Laboratory (FNAL) Muon $g-2$ Experiment [4] with the previous one
at Brookhaven [5], allows to combine them and get

$a_{\mu}^{\mathrm{Exp}}=116\,592\,061(41)\times 10^{-11}\quad(0.35\text{
ppm}).$

The latest SM prediction [6] is

$a_{\mu}^{\mathrm{SM}}=116\,591\,810(43)\times 10^{-11}\quad(0.37\text{
ppm}).$

Therefore, the difference, $\Delta
a_{\mu}=a_{\mu}^{\text{Exp}}-a_{\mu}^{\text{SM}}=(251\pm 59)\times 10^{-11}$,
increases the significance at $4.2\,\sigma$. This deviation from the Standard
Model could be a sign of New Physics around the corner.

The uncertainty of the SM prediction is dominated by the hadronic
contributions. Specifically, the HVP, LO contribution is dominated by the
$\pi\pi$ cut (about $\sim 70\%$ of the overall value). At present, the most
competitive estimation is obtained via dispersion relations together with
$e^{+}e^{-}\to\text{hadrons}$ data. Alternatively, the CVC relation between
electromagnetic and weak form factors in the isospin-limit allows using tau-
data for this purpose. When tau data is employed for $a_{\mu}^{HVP,LO}$, the
discrepancy between the SM prediction and the experimental measurement amounts
to $2.4\,\sigma$ [7]. The difference with respect to the $e^{+}e^{-}$-based
evaluation could be owing to NP effects, hinting at a lepton universality
violation in the corresponding non-standard vector and/or tensor couplings at
low-energies [8, 9, 10]. Nowadays, $e^{+}e^{-}$-based estimation has
superseded the tau-based one due to the very high precision achieved in the
$\sigma(e^{+}e^{-}\to\text{ hadrons})$ measurements.

Additionally to the data-driven approach, lattice QCD determinations of
$a_{\mu}^{\text{HVP, LO}}$ have achieved a significant improvement. Although
the lattice estimations are not yet competitive with the $e^{+}e^{-}$
evaluations [6], a very accurate computation made by the BMW coll. [11]
alleviates the tension concerning the SM prediction at one sigma level.

The paper is organized as follows. In Sec. 2 we review the main features of
the $\tau^{-}\to\pi^{-}\pi^{0}\nu_{\tau}\gamma$ decays as well as the
theoretical framework. We show some decay observables in Sec. 3. Then, in Sec.
4 we evaluate $a_{\mu}^{\text{HVP, LO}|_{\pi\pi}}$ using tau data. Finally,
our conclusions are presented in Sec. 5.

##  2 $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$ decays

###  2.1 Amplitude

For the radiative decay
$\tau^{-}\left(P\right)\to\pi^{-}\left(p_{-}\right)\pi^{0}\left(p_{0}\right)\nu_{\tau}\left(q\right)\gamma\left(k\right)$,
we can split the contribution due to the bremsstrahlung off the initial tau
lepton from the one coming from the hadronic part.

The most general structure for these decays can be written as [12, 13]

$\begin{split}T=e\,G_{F}V_{ud}^{*}\epsilon^{\mu}(k)^{*}&\left\\{F_{\nu}\,\bar{u}\left(q\right)\gamma^{\nu}\left(1-\gamma_{5}\right)\left(m_{\tau}+\not{P}-\not{k}\right)\gamma_{\mu}\,u\left(P\right)\right.\\\
&\left.+\left(V_{\mu\nu}-A_{\mu\nu}\right)\bar{u}\left(q\right)\gamma^{\nu}\left(1-\gamma_{5}\right)u\left(P\right)\right\\},\end{split}$
(1)

where $F_{\nu}\equiv\left(p_{0}-p_{-}\right)_{\nu}f_{+}\left(s\right)/2P\cdot
k$, with the charged pion vector form factor $f_{+}(s)$ defined through
$\left\langle\pi^{0}\pi^{-}|\bar{d}\gamma^{\mu}u|0\right\rangle=\sqrt{2}f_{+}(s)(p_{-}-p_{0})^{\mu}$
and $s=(p_{-}+p_{0})^{2}$. Gauge invariance
($\epsilon_{\mu}\to\epsilon_{\mu}+k_{\mu}$) implies the Ward identities

$k_{\mu}V^{\mu\nu}=\left(p_{-}-p_{0}\right)^{\nu}f_{+}\left(s\right),\quad
k_{\mu}A^{\mu\nu}=0.$ (2)

Imposing eq. (2) and Lorentz invariance, we have the following expression for
the vector structure-dependent tensor

$\small\begin{split}V^{\mu\nu}&=f_{+}\left[\left(P-q\right)^{2}\right]\frac{p_{-}^{\mu}\left(p_{-}+k-p_{0}\right)^{\nu}}{p_{-}\cdot
k}-f_{+}\left[\left(P-q\right)^{2}\right]g^{\mu\nu}\\\
&+\frac{f_{+}\left[\left(P-q\right)^{2}\right]-f_{+}\left(s\right)}{\left(p_{0}+p_{-}\right)\cdot
k}\left(p_{0}+p_{-}\right)^{\mu}\left(p_{0}-p_{-}\right)^{\nu}\\\
&+v_{1}\left(g^{\mu\nu}\,p_{-}\cdot
k-p_{-}^{\mu}k^{\nu}\right)+v_{2}\left(g^{\mu\nu}\,p_{0}\cdot
k-p_{0}^{\mu}k^{\nu}\right)\\\ &+v_{3}\left(p_{0}\cdot
k\,p_{-}^{\mu}-p_{-}\cdot
k\,p_{0}^{\mu}\right)p_{-}^{\nu}+v_{4}\left(p_{0}\cdot
k\,p_{-}^{\mu}-p_{-}\cdot
k\,p_{0}^{\mu}\right)\left(p_{0}+p_{-}+k\right)^{\nu},\end{split}$ (3)

and for the axial one

$\begin{split}A^{\mu\nu}&=ia_{1}\,\epsilon^{\mu\nu\rho\sigma}\,\left(p_{0}-p_{-}\right)_{\rho}k_{\sigma}+ia_{2}\,W^{\nu}\,\epsilon^{\mu\lambda\rho\sigma}k_{\lambda}\,p_{-\rho}\,p_{0\sigma}\\\
&+ia_{3}\,\epsilon^{\mu\nu\rho\sigma}k_{\rho}\,W_{\sigma}+ia_{4}\,\left(p_{0}+k\right)^{\nu}\,\epsilon^{\mu\lambda\rho\sigma}\,k_{\lambda}\,p_{-\rho}\,p_{0\sigma},\\\
\end{split}$ (4)

where $W\equiv P-q=p_{-}+p_{0}+k$. The structure-dependent contributions to
these tensor structures are contained in the four vector ($v_{i}$) and the
four axial-vector ($a_{i}$) form factors. For the axial structure, the
Schouten’s identity has been used.

Using $\left(P-q\right)^{2}=s+2\left(p_{0}+p_{-}\right)\cdot k$, it is easy to
show that Low’s theorem [14] is manifestly satisfied

$\begin{split}V^{\mu\nu}&=f_{+}\left(s\right)\frac{p_{-}^{\mu}}{p_{-}\cdot
k}\left(p_{-}-p_{0}\right)^{\nu}+f_{+}\left(s\right)\left(\frac{p_{-}^{\mu}k^{\nu}}{p_{-}\cdot
k}-g^{\mu\nu}\right)\\\
&+2\frac{df_{+}\left(s\right)}{d\,s}\left(\frac{p_{0}\cdot k}{p_{-}\cdot
k}p_{-}^{\mu}-p_{0}^{\mu}\right)\left(p_{-}-p_{0}\right)^{\nu}+\mathcal{O}\left(k\right).\end{split}$
(5)

###  2.2 Vector Form Factors

In $R\chi T$ [15, 16, 17, 18], the diagrams that contribute to the vector form
factors of the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$ decays are shown in
Figs. 1, 2 and 3.

$\rho^{-}$$\otimes$$\pi^{-}$$\gamma$$\pi^{0}$

$\rho^{-}$$\gamma$$\otimes$$\pi^{-}$$\pi^{0}$

$\rho^{-}$$\otimes$$\pi^{-}$$\gamma$$\pi^{0}$

$\rho^{0}$$\pi^{-}$$\pi^{0}$$\otimes$$\gamma$

$\omega$$\pi^{-}$$\otimes$$\pi^{0}$$\gamma$

$a_{1}^{-}$$\pi^{0}$$\otimes$$\pi^{-}$$\gamma$

$\pi^{-}$$\rho^{0}$$\pi^{0}$$\otimes$$\pi^{-}$$\gamma$

Figure 1: One-resonance exchange contributions from the $R\chi T$ to the
vector form factors of the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$ decays.

$\rho^{-}$$\rho^{-}$$\otimes$$\gamma$$\pi^{-}$$\pi^{0}$

$\rho^{-}$$\rho^{0}$$\otimes$$\pi^{-}$$\pi^{0}$$\gamma$

$\rho^{-}$$\omega$$\otimes$$\pi^{-}$$\pi^{0}$$\gamma$

$\rho^{0}$$\omega$$\pi^{-}$$\otimes$$\pi^{0}$$\gamma$

$\rho^{0}$$\gamma$$\rho^{-}$$\otimes$$\pi^{-}$$\pi^{0}$

$\rho^{-}$$a_{1}^{-}$$\otimes$$\pi^{0}$$\pi^{-}$$\gamma$

$\pi^{0}$$a_{1}^{-}$$\rho^{0}$$\otimes$$\pi^{-}$$\gamma$

$\otimes$$\rho^{-}$$\pi^{0}$$\pi^{-}$$\rho^{0}$$\pi^{-}$$\gamma$

Figure 2: Two-resonance exchange contributions from the $R\chi T$ to the
vector form factors of the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$ decays.
$\otimes$$\rho^{-}$$\pi^{0}$$a_{1}^{-}$$\rho^{0}$$\pi^{-}$$\gamma$ Figure 3:
Three-resonance exchange contributions from the $R\chi T$ to the vector form
factors of the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$ decays.

For the vector form factors, we have

$\displaystyle v_{1}$
$\displaystyle=v_{1}^{0}+v_{1}^{R}+v_{1}^{RR}+v_{1}^{RRR}+v_{GI1}^{R+RR},$
(6a) $\displaystyle v_{2}$
$\displaystyle=v_{2}^{0}+v_{2}^{R}+v_{2}^{RR}+v_{2}^{RRR}+v_{GI2}^{R+RR},$
(6b) $\displaystyle v_{3}$
$\displaystyle=v_{3}^{0}+v_{3}^{R}+v_{3}^{RR}+v_{3}^{RRR}+v_{GI3}^{R+RR},$
(6c) $\displaystyle v_{4}$
$\displaystyle=v_{4}^{0}+v_{4}^{R}+v_{4}^{RR}+v_{4}^{RRR}+v_{GI4}^{R+RR},$
(6d)

where $v_{i}^{0}$ is the $\mathcal{O}(p^{4})$ contribution in Ref. [13], and
$v_{i}^{R}$, $v_{i}^{RR}$, $v_{i}^{RRR}$ and $v_{GIi}^{R+RR}$, which are the
subleading contributions up to $\mathcal{O}(p^{6})$, can be found in App. C.
in Ref. [19].

###  2.3 Axial-Vector Form Factors

The axial form factors at chiral $\mathcal{O}\left(p^{4}\right)$ [13] get
contributions from the Wess-Zumino-Witten functional [20, 21]:

$a_{1}^{0}\equiv\frac{1}{8\pi^{2}F^{2}},\qquad
a_{2}^{0}\equiv\frac{-1}{4\pi^{2}F^{2}\left[\left(P-q\right)^{2}-m_{\pi}^{2}\right]}.$
(7)

The diagrams contributing to these two expressions are shown in fig. 4.

$\pi^{-}$$\gamma$$\pi^{0}$$\otimes$

$\pi^{-}$$\pi^{-}$$\gamma$$\pi^{0}$$\otimes$

Figure 4: Anomalous diagrams contributing to the axial tensor amplitude
$A^{\mu\nu}$ at $\mathcal{O}\left(p^{4}\right)$.

$\rho^{-}$$\gamma$$\otimes$$\pi^{-}$$\pi^{0}$

$\rho^{-}$$\pi^{0}$$\otimes$$\pi^{-}$$\gamma$

$\rho^{0}$$\pi^{-}$$\otimes$$\pi^{0}$$\gamma$

$\omega$$\pi^{-}$$\pi^{0}$$\otimes$$\gamma$

$\rho^{0}$$\pi^{-}$$\otimes$$\pi^{-}$$\pi^{0}$$\gamma$

$\pi^{-}$$\rho^{-}$$\otimes$$\pi^{0}$$\pi^{-}$$\gamma$

$\pi^{-}$$\rho^{-}$$\otimes$$\gamma$$\pi^{-}$$\pi^{0}$

$\pi^{-}$$\omega$$\otimes$$\pi^{-}$$\pi^{0}$$\gamma$

$a_{1}^{-}$$\otimes$$\pi^{-}$$\gamma$$\pi^{0}$

Figure 5: One-resonance exchange contributions from the $R\chi T$ to the
axial-vector form factors of the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$
decays.

$\rho^{0}$$\omega$$\pi^{-}$$\otimes$$\pi^{0}$$\gamma$

$\rho^{-}$$\omega$$\pi^{0}$$\otimes$$\pi^{-}$$\gamma$

$\omega$$\gamma$$\rho^{-}$$\otimes$$\pi^{-}$$\pi^{0}$

$\otimes$$\pi^{-}$$\pi^{-}$$\rho^{0}$$\omega$$\pi^{0}$$\gamma$

$\otimes$$\pi^{-}$$\pi^{0}$$\rho^{-}$$\omega$$\pi^{-}$$\gamma$

$a_{1}^{-}$$\rho^{-}$$\otimes$$\pi^{0}$$\pi^{-}$$\gamma$

$a_{1}^{-}$$\rho^{0}$$\otimes$$\pi^{-}$$\pi^{0}$$\gamma$

$a_{1}^{-}$$\rho^{-}$$\otimes$$\gamma$$\pi^{-}$$\pi^{0}$

$a_{1}^{-}$$\omega$$\otimes$$\pi^{-}$$\pi^{0}$$\gamma$

Figure 6: Two-resonance exchange contributions from the $R\chi T$ to the
axial-vector form factors of the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$
decays.

$\otimes$$a_{1}^{-}$$\pi^{0}$$\rho^{-}$$\omega$$\pi^{-}$$\gamma$

$\otimes$$a_{1}^{-}$$\pi^{-}$$\rho^{0}$$\omega$$\pi^{0}$$\gamma$

Figure 7: Three-resonance exchange contributions from the $R\chi T$ to the
axial-vector form factors of the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$
decays.

For the axial form factors, we have

$\displaystyle a_{1}$
$\displaystyle=a_{1}^{0}+a_{1}^{R}+a_{1}^{RR}+a_{1}^{RRR},$ (8a)
$\displaystyle a_{2}$
$\displaystyle=a_{2}^{0}+a_{2}^{R}+a_{2}^{RR}+a_{2}^{RRR},$ (8b)
$\displaystyle a_{3}$ $\displaystyle=a_{3}^{R}+a_{3}^{RR}+a_{3}^{RRR},$ (8c)
$\displaystyle a_{4}$ $\displaystyle=a_{4}^{R}+a_{4}^{RR}+a_{4}^{RRR},$ (8d)

where $a_{i}^{R}$, $a_{i}^{RR}$ and $a_{i}^{RRR}$ include up to
$\mathcal{O}\left(p^{6}\right)$ contributions. The diagrams contributing to
the Eq. (8) are shown in figures 4-7. These expressions can be found in App.
D. in Ref. [19].

###  2.4 SD constraints

Including the complete set of operators [17, 18] that start contributing to
the $\mathcal{O}(p^{6})$ LECs, we have too many parameters which are allowed
by the discrete symmetries of QCD and chiral symmetry that prevent us from
making phenomenology predictions. For a detailed discussion you can see
Section 1.5 in Ref. [19].

Imposing the asymptotic behavior of: the pion vector form factor, the $V-A$
correlator, the scalar form factor, and the $S-P$ correlator, the following
constraints are found [15, 16, 22, 23, 24, 25, 26] 111Other important
relations are in App. A. :

$\small\begin{split}F_{V}G_{V}=F^{2},&\qquad F_{V}^{2}-F_{A}^{2}=F^{2},\\\
F_{V}^{2}M_{V}^{2}=F_{A}^{2}M_{A}^{2},&\qquad 4c_{d}c_{m}=F^{2},\\\
8\left(c_{m}^{2}-d_{m}^{2}\right)=F^{2},&\qquad
c_{m}=c_{d}=\sqrt{2}d_{m}=F/2.\end{split}$ (9)

The asymptotic behavior of the $2-$point Green function at
$\mathcal{O}(p^{4})$ predicts,

$F_{V}=\sqrt{2}F\,,\quad G_{V}=\frac{F}{\sqrt{2}}\,,\quad F_{A}=F\,,$ (10)

for the couplings of the $R\chi T$ Lagrangian [16]. Conversely, the $2$\- and
$3$-point Green function up to $\mathcal{O}(p^{6})$ [17, 18, 27] 222The
contributions from operators with more than one resonance field are taken into
account. determine

$F_{V}=\sqrt{3}F\,,\quad G_{V}=\frac{F}{\sqrt{3}}\,,\quad F_{A}=\sqrt{2}F\,.$
(11)

We will hereinafter refer to the constraints from the $2$\- and $3$-point
Green functions as ‘$F_{V}=\sqrt{2}F$’ and ‘$F_{V}=\sqrt{3}F$’, respectively.

Since the $\kappa_{i}^{V}$ couplings are related with the $\omega$ exchange
which is known to give an important contribution to the
$\tau\to\pi\pi\gamma\nu_{\tau}$ decays, we perform a global fit using the
relations for the resonance saturation of the anomalous sector LECs at NLO
[18], the eqs. in App. A and the estimation of the LECs in [28]. The fit
outcomes are in App. B.

##  3 Decay observables

The differential decay width [29] is given by

$d\Gamma=\frac{\lambda^{1/2}\left(s,m_{\pi^{0}}^{2},m_{\pi^{-}}^{2}\right)}{2\left(4\pi\right)^{6}m_{\tau}^{2}s}\overline{\left|\mathcal{M}\right|^{2}}\,dE_{\gamma}\,dx\,ds\,d\cos\theta_{-}\,d\phi_{-},$
(12)

where $\overline{\left|\mathcal{M}\right|^{2}}$ is the unpolarized spin-
averaged squared amplitude that corresponds to the
$\tau^{-}\to\pi^{-}\pi^{0}\gamma\,\nu_{\tau}$ decays, and $E_{\gamma}$ is the
photon energy in the $\tau$ rest frame. It is not worth to write down here the
full analytical expression for $\overline{\left|\mathcal{M}\right|^{2}}$. The
kinematics of this process can be found in App. C

###  3.1 Decay spectrum

Integrating the Eq. (12) over $E_{\gamma}$, $x$, $\cos\theta_{-}$ and
$\phi_{-}$, we obtain the $\pi^{-}\pi^{0}$ hadronic invariant distribution.
Since the decay spectrum is IR divergent due to soft photons, we require to
introduce a photon energy cut, $E_{\gamma}^{cut}$, which is related to the
experimental resolution.

In figure 8, we can see the prediction for the decay spectrum for
$E_{\gamma}^{cut}=300\,\text{MeV}$. The dotted line indicates the limit where
all the structure-dependent form factors vanish, i.e. $v_{i}=a_{i}=0$. The
predictions at $\mathcal{O}(p^{4})$ using $F_{V}=\sqrt{2}F$ and
$F_{V}=\sqrt{3}F$, which are discussed in Sec. 2.4, are denoted by the dashed
and solid line, respectively. The dotdashed red line corresponds to taking the
limit where all the couplings at $\mathcal{O}(p^{6})$ vanish except for those
constrained by SD and the band overestimates the corresponding uncertainties.
Including the $\mathcal{O}(p^{6})$ corrections, the decay spectrum receives a
noticeable enhancement at low $s$.

Figure 8: The $\pi^{-}\pi^{0}$ hadronic invariant mass distributions for
$E_{\gamma}^{cut}=300\,\text{MeV}$. The solid and dashed lines represent the
$\mathcal{O}\left(p^{4}\right)$ corrections using $F_{V}=\sqrt{3}F$ and
$F_{V}=\sqrt{2}F$, respectively. The dotted line represents the Bremsstrahlung
contribution (SI). The dotdashed red line corresponds to using only SD
constraints at $\mathcal{O}\left(p^{6}\right)$ and the blue shaded region
overestimates the corresponding uncertainties.

###  3.2 Branching ratio

By integrating upon the $s$ variable the decay spectrum, we obtain the
branching ratio for the $\tau\to\pi\pi\gamma\nu_{\tau}$ decays. Figure 9 shows
the prediction for the branching ratio as a function of $E_{\gamma}^{cut}$
obtained using the different orders of approximation for the structure-
dependent terms in Eqs. (3) and (4).

According to figs. 8 and 9, measurements of the $\pi\pi$ invariant mass and
the partial decay width, for a reasonable cut on $E_{\gamma}$ (at low enough
energies the inner bremmstrahlung contribution hides completely any structure-
dependent effect), could decrease substantially the uncertainty of the
$\mathcal{O}\left(p^{6}\right)$ computation 333The photon spectrum measurement
would also help to this task [19].. This was already emphasized in Ref. [13]
but remained unmeasured at BaBar and Belle. We hope these data can finally be
acquired and analyzed at Belle-II.

Figure 9: Branching ratio for the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$
decays as a function of $E^{cut}_{\gamma}$. The dotted line represents the
Bremsstrahlung contribution, the solid and dashed lines represent the
$\mathcal{O}\left(p^{4}\right)$ corrections using $F_{V}=\sqrt{3}F$ and
$F_{V}=\sqrt{2}F$, respectively. The dotdashed red line is the
$\mathcal{O}\left(p^{6}\right)$ contribution using only SD constraints and
neglecting all other couplings. The blue shaded region overestimates the
$\mathcal{O}\left(p^{6}\right)$ uncertainties.

##  4 Hadronic Vacuum Polarization

The leading contributions to the hadronic vacuum polarization (HVP) employing
a dispersion relation [30], which complies with unitarity and analyticity, are
given by

$a^{HVP,LO}_{\mu}=\frac{1}{4\pi^{3}}\int^{\infty}_{s_{thr}}ds\,K(s)\sigma^{0}_{e^{-}e^{+}\to
hadrons}(s),$ (13)

where $K(s)$ is a smooth QED kernel [31] concentrated at low energies, which
increases the $E\lesssim M_{\rho}$ contribution, and $\sigma_{e^{-}e^{+}\to
hadrons}^{0}(s)$ is the bare hadronic cross section. We can relate the
hadronic spectral function from $\tau$ decays to the $e^{+}e^{-}$ hadronic
cross section by including the radiative corrections and the IB effects.

For the $\pi\pi(\gamma)$ final state, we have [32, 13]:

$\sigma_{\pi\pi(\gamma)}^{0}=\left[\frac{K_{\sigma}(s)}{K_{\Gamma}(s)}\frac{d\Gamma_{\pi\pi[\gamma]}}{ds}\right]\frac{R_{IB}(s)}{S_{EW}},$
(14)

where

$\begin{split}K_{\Gamma}(s)&=\frac{G_{F}^{2}|V_{ud}|^{2}m_{\tau}^{3}}{384\pi^{3}}\left(1-\frac{s}{m_{\tau}^{2}}\right)^{2}\left(1+\frac{2s}{m_{\tau}^{2}}\right),\quad
K_{\sigma}(s)=\frac{\pi\alpha^{2}}{3s},\end{split}$ (15)

and the IB corrections

$R_{IB}(s)=\frac{FSR(s)}{G_{EM}(s)}\frac{\beta^{3}_{\pi^{+}\pi^{-}}}{\beta^{3}_{\pi^{+}\pi^{0}}}\left|\frac{F_{V}(s)}{f_{+}(s)}\right|^{2}.$
(16)

The $S_{EW}$ term encodes the SD electroweak corrections [33, 34, 35, 36, 37,
38, 39, 40] and $FSR(s)$ accounts for the final-state radiation from pions
[41, 42]. The $\beta^{3}_{\pi^{+}\pi^{-}}/\beta^{3}_{\pi^{+}\pi^{0}}$ term is
a phase space factor and the last term in $R_{IB}(s)$ is a ratio between the
neutral ($F_{V}(s)$) and the charged ($f_{+}(s)$) pion form factor.

The $G_{EM}(s)$ function is obtained by adding up the contributions due to
virtual and real photons and integrating over the $u\equiv(P-p_{-})^{2}$
variable,

$\begin{split}\left.\frac{d\Gamma}{ds}\right|_{\pi\pi(\gamma)}=&\frac{G_{F}^{2}|V_{ud}|^{2}m_{\tau}^{3}S_{EW}}{384\pi^{3}}\left|f_{+}(s)\right|^{2}\left(1-\frac{s}{m_{\tau}^{2}}\right)^{2}\left(1-\frac{4m_{\pi}^{2}}{s}\right)^{3/2}\left(1+\frac{2s}{m_{\tau}^{2}}\right)G_{EM}(s).\end{split}$
(17)

The results are shown in figure 10 for the different approximations of the
$G_{EM}(s)$ function using the dispersive (left-hand) and the exponential
(right-handed) representation of the pion form factor. The $G_{EM}^{0}(s)$
contribution was obtained using the leading Low approximation in Eq. (5).

Figure 10: Correction function $G_{EM}^{(0)}\left(s\right)$ (blue dashed
line). The solid line shows the $G_{EM}(s)$ function neglecting the structure-
dependent part (SI), i.e. by taking $v_{i}=a_{i}=0$, the dashed and dotted
lines are the $\mathcal{O}\left(p^{4}\right)$ $G_{EM}(s)$ function (with
either $F_{V}=\sqrt{2}F$ or $F_{V}=\sqrt{3}F$ constraints). The blue shaded
region is the full $\mathcal{O}\left(p^{6}\right)$ contribution, including
(overestimated) uncertainties. The left-hand side plot corresponds to the
dispersive parametrization [43] while the right-hand side corresponds to the
Guerrero-Pich parametrization [44] of the form factor (the latter was used in
Ref. [13]).

We can estimate the effect of each IB correction through $\Delta
a_{\mu}^{HVP,LO}[\pi\pi]$ [13],

$\Delta
a_{\mu}^{HVP,LO}=\frac{1}{4\pi^{3}}\int_{s_{1}}^{s_{2}}ds\,K(s)\left[\frac{K_{\sigma}(s)}{K_{\Gamma}(s)}\frac{d\Gamma_{\pi\pi[\gamma]}}{ds}\right]\left(\frac{R_{IB}(s)}{S_{EW}}-1\right),$
(18)

which measures the departure from the isospin-limit, i.e. $R_{IB}(s)=1$ and
$S_{EW}=1$. It is challenging to evaluate the corrections owing to the ratio
of the form factors. For this enterprise, we have followed two alternatives,
dubbed FF1 and FF2:

  * •

In FF1, we use for the $\rho-\omega$ mixing parameter
$\theta_{\rho\omega}=(-3.5\pm 0.7)\times 10^{-3}\>\mathrm{GeV}^{2}$ [13] and
$\Gamma_{\rho^{0}}-\Gamma_{\rho^{+}}=0.3\pm 1.3$ MeV,
$m_{\rho^{\pm}}-m_{\rho^{0}}=0.7\pm 0.8$ MeV and $m_{\rho^{0}}=775.26\pm 0.25$
MeV from PDG [45].

  * •

In FF2, we use the same numerical input as in FF1 except by the rho width,
which is
$\Gamma_{\rho^{0}\to\pi^{+}\pi^{-}\gamma}-\Gamma_{\rho^{\pm}\to\pi^{\pm}\pi^{0}\gamma}=0.45\pm
0.45\,\text{MeV}$ [13].

The outcomes are summarized in Table 1 using DR form factor. The results
obtained for the $G_{EM}^{(0)}$ and the complete
$\mathcal{O}\left(p^{4}\right)$ contribution (with $F_{V}=\sqrt{2}F$) agree
with those in Ref. [13]. The uncertainties at $\mathcal{O}(p^{4})$ were
obtained using the dashed and dotdashed red line in fig. 10. On the other
hand, the errors at $\mathcal{O}(p^{6})$ were estimated using the blue region
in the same plot.

$\left[s_{1},s_{2}\right]$ | $\mathrm{S_{EW}}$ | $\mathrm{PS}$ | $\mathrm{FSR}$ | $\mathrm{FF1}$ | $\mathrm{FF2}$ | $\mathrm{EM}$ | $\mathrm{EM}$
---|---|---|---|---|---|---|---
|  |  |  |  |  | $\,\mathcal{O}(p^{4})$ | $\,\mathcal{O}(p^{6})$
$\left[4m_{\pi}^{2},1\,\mathrm{GeV}^{2}\right]$ | $-101.1$ | $-74.1$ | $+44.7$ | $+41.8\pm 49.0$ | $+78.4\pm 24.5$ | $-17.0^{+5.7}_{-15.4}$ | $-74.8^{+63.5}_{-44.0}$
$\left[4m_{\pi}^{2},2\,\mathrm{GeV}^{2}\right]$ | $-103.1$ | $-74.4$ | $+45.5$ | $+40.9\pm 48.9$ | $+77.6\pm 24.0$ | $-16.0^{+5.7}_{-15.9}$ | $-75.9^{+65.6}_{-45.5}$
$\left[4m_{\pi}^{2},3\,\mathrm{GeV}^{2}\right]$ | $-103.1$ | $-74.5$ | $+45.5$ | $+40.9\pm 48.9$ | $+77.6\pm 24.0$ | $-15.9^{+5.7}_{-16.0}$ | $-75.9^{+65.7}_{-44.6}$
$\left[4m_{\pi}^{2},m_{\tau}^{2}\right]$ | $-103.1$ | $-74.5$ | $+45.5$ | $+40.9\pm 48.9$ | $+77.6\pm 24.0$ | $-15.9^{+5.7}_{-16.0}$ | $-75.9^{+65.7}_{-45.6}$
Table 1: Contributions to $\Delta a_{\mu}^{HVP,LO}$ in units of $10^{-11}$
using the DR form factor as the reference one.

An important cross-check is the branching fraction
$B_{\pi\pi^{0}}=\Gamma(\tau\to\pi\pi^{0}\nu_{\tau})/\Gamma_{\tau}$ which can
be directly measured in experiments. It can also be evaluated from the
isovector component of the $e^{+}e^{-}\to\pi^{+}\pi^{-}(\gamma)$ cross section
after taking into account the IB corrections. The branching fraction is given
by

$B_{\pi\pi^{0}}^{CVC}=B_{e}\int_{4m_{\pi}^{2}}^{m_{\tau}^{2}}ds\,\sigma_{\pi^{+}\pi^{-}(\gamma)}(s)\mathcal{N}(s)\frac{S_{EW}}{R_{IB}(s)},$
(19)

where

$\mathcal{N}(s)=\frac{3\left|V_{ud}\right|^{2}}{2\pi\alpha_{0}^{2}m_{\tau}^{2}}s\left(1-\frac{s}{m_{\tau}^{2}}\right)^{2}\left(1+\frac{2s}{m_{\tau}^{2}}\right).$
(20)

Using the most recent data obtained from BaBar [46] 444We thank to Alex
Keshavarzi and Bogdan Malaescu for providing us tables with the measurement of
the $e^{+}e^{-}\to\pi^{+}\pi^{-}(\gamma)$ cross section. for the
$e^{+}e^{-}\to\pi^{+}\pi^{-}(\gamma)$ cross section, we obtain

$B_{\pi\pi^{0}}^{CVC}=(24.68\pm 0.11\pm 0.10\pm 0.01\pm 0.01\pm
0.02^{+0.03}_{-0.00})\,\%,\text{ at }\mathcal{O}(p^{4}),$ (21)

and

$B_{\pi\pi^{0}}^{CVC}=(24.70\pm 0.11\pm 0.10\pm 0.01\pm 0.01\pm
0.02^{+0.21}_{-0.01})\,\%,\text{ at }\mathcal{O}(p^{6}),$ (22)

where the first error corresponds to the statistical experimental uncertainty
on $\sigma_{\pi\pi(\gamma)}$, the second is related to uncertainty on the
$\rho^{+}-\rho^{0}$ width difference, the third to the uncertainty in the
$\rho^{+}-\rho^{0}$ mass difference, the fourth to the uncertainty of the
$\rho-\omega$ mixing and the fifth corresponds to the corrections induced by
FSR on $B_{\pi\pi^{0}}^{CVC}$, which reduces $\sim-0.20(2)\%$ the $\pi\pi$
branching fraction. The last error is related to the $G_{EM}(s)$ function.

These results are in good agreement with the value reported by the Belle [47]
collaboration, $B_{\pi\pi^{0}}^{\tau}=(25.24\pm 0.01\pm 0.39)\%$, where the
first uncertainty is statistical and the second is systematic. Nonetheless,
they are in some tension with the very precise ALEPH measurement
$B_{\pi\pi^{0}}^{\tau}=(25.471\pm 0.097\pm 0.085)\%$ [48].

We show in fig. 11 the prediction for the $e^{+}e^{-}\to\pi^{+}\pi^{-}$ cross
section using the data reported by Belle [47] (as it is the most precise
measurement of this spectrum) for the normalized spectrum
$(1/N_{\pi\pi})(dN_{\pi\pi}/ds)$ compared to the last measurements from BaBar
[46] and KLOE [49] 555We have chosen to show in the comparison these two
$e^{+}e^{-}$ data sets as the results from both Colls. are those deviating the
most, and thus mainly responsible from the tension in
$\sigma(e^{+}e^{-}\to\pi^{+}\pi^{-})$..

In fig. 11 the $\tau$-based prediction is obtained using the
$\mathcal{O}(p^{4})$ result for $G_{EM}(s)$, with the estimated uncertainty
from missing higher-order corrections given by the result at
$\mathcal{O}(p^{6})$ (employing only the SD constraints). The blue dotdashed
line shown overestimates the error at $\mathcal{O}(p^{6})$.

Figure 11: Comparison between the different data sets from BaBar (above) and
KLOE (below) with $\Delta\Gamma_{\pi\pi\gamma}=1.5\,\mathrm{MeV}$ (left-hand)
and $\Delta\Gamma_{\pi\pi\gamma}=0.45\,\mathrm{MeV}$ (right-hand) for FF1 and
FF2, respectively. The blue region corresponds to the experimental error on
$\sigma_{\pi\pi(\gamma)}$. The solid and dashed lines represent the
contributions with $F_{V}=\sqrt{3}F$ and $F_{V}=\sqrt{2}F$ at
$\mathcal{O}(p^{4})$, respectively. The dotted line is the SI contribution.
The red line depicts the envelope of $G_{EM}(s)$ at $\mathcal{O}(p^{6})$, that
overestimates the uncertainty at this order. The blue dotdashed line is the
$\mathcal{O}(p^{6})$ contribution using only SD constraints.

From fig. 11, we observe good agreement between the BaBar data and the $\tau$
decays prediction (slightly better for FF1). The previous comparisons make us
consider our evaluation with FF1 the reference one (so that its difference
with FF2 will assess the size of the error induced by IB among the
$\rho\to\pi\pi\gamma$ decay channels) 666We, nevertheless, recall that recent
BESIII data [50] and evaluations within the Hidden Local Symmetry model [51,
52, 53, 54] agree better with the KLOE data than with BaBar’s..

Taking into account all di-pion tau decay data from the ALEPH [48], Belle
[47], CLEO [55] and OPAL [56] Colls., we get the combined tau-data
contribution

$10^{10}\cdot
a_{\mu}^{HVP,LO|_{\pi\pi,\tau\;\mathrm{data}}}\,=\,519.6\pm{2.8_{\mathrm{spectra+BRs}}}{{}^{+1.9}_{-2.1}}_{\mathrm{IB}}\,,$
(23)

at $\mathcal{O}(p^{4})$ and

$10^{10}\cdot
a_{\mu}^{HVP,LO|_{\pi\pi,\tau\;\mathrm{data}}}\,=\,514.6\pm{2.8_{\mathrm{spectra+BRs}}}{{}^{+5.0}_{-3.9}}_{\mathrm{IB}}\,,$
(24)

at $\mathcal{O}(p^{6})$.

When eqs. (23) and (24) are supplemented with the four-pion tau decays
measurements (up to $1.5$ GeV) and with $e^{+}e^{-}$ data at larger energies
in these modes (and with $e^{+}e^{-}$ data in all other channels making up the
hadronic cross section), we get [57, 7]

$10^{10}\,\cdot\,a_{\mu}^{HVP,LO|_{\tau\;\mathrm{data}}}\,=\,705.7\pm{2.8_{\mathrm{spectra+BRs}}}{{}^{+1.9}_{-2.1}}_{\mathrm{IB}}\pm
2.0_{\mathrm{e^{+}e^{-}}}\pm 0.1_{\mathrm{narrow\,res}}\pm
0.7_{\mathrm{QCD}}\,,$ (25)

at $\mathcal{O}(p^{4})$, and

$10^{10}\,\cdot\,a_{\mu}^{HVP,LO|_{\tau\;\mathrm{data}}}\,=\,700.7\pm{2.8_{\mathrm{spectra+BRs}}}{{}^{+5.0}_{-3.9}}_{\mathrm{IB}}\pm
2.0_{\mathrm{e^{+}e^{-}}}\pm 0.1_{\mathrm{narrow\,res}}\pm
0.7_{\mathrm{QCD}}\,,$ (26)

at $\mathcal{O}(p^{6})$ and we have also included the uncertainties
corresponding to using $e^{+}e^{-}$ data for those contributions not covered
by tau decay measurements and to the inclusion of narrow resonances and the
perturbative QCD part.

When all other (QED, EW and subleading hadronic) contributions are added to
eqs. (25) and (26) according to Ref. [6], the $4.2\sigma$ [6] deficit of the
SM prediction with respect to the experimental average (FNAL+BNL) [4, 5] is
reduced to

$\Delta a_{\mu}\equiv a_{\mu}^{exp}-a_{\mu}^{SM}=(12.5\pm 6.0)\cdot
10^{-10}\,,$ (27)

at $\mathcal{O}(p^{4})$, and

$\Delta a_{\mu}\equiv a_{\mu}^{exp}-a_{\mu}^{SM}=(17.5^{+6.8}_{-7.5})\cdot
10^{-10}\,,$ (28)

at $\mathcal{O}(p^{6})$, which are $2.1$ and $2.3\,\sigma$, respectively.

In figure 12 we show a comparison between our $\mathcal{O}(p^{4})$ and
$\mathcal{O}(p^{6})$ calculation with respect to the estimation based in the
$e^{+}e^{-}$ data driven [6] and the lattice results from the BMW
collaboration [11].

Figure 12: Comparison between the experimental values of $a_{\mu}$ from BNL
[5] and FNAL [4] with respect to the Muon $g-2$ Theory Initiative recommended
value [6], the lattice QCD calculation from the BMW collaboration [11] and our
results [19].

##  5 Conclusions

There is a global effort in improving the hadronic contributions to $a_{\mu}$.
Specifically, dedicated studies to improve the HVP part from lattice,
dispersion relations and improved $e^{+}e^{-}$ data and Monte Carlos are being
undertaken.

The observables for the $\tau\to\pi\pi\gamma\nu_{\tau}$ decays have the
potential to reduce drastically the errors in our estimation.

Our IB corrections improve the agreement between $e^{+}e^{-}$ and tau data, on
the spectrum and the branching ratio.

Evaluating the HVP, LO contributions from tau data, we get
$a_{\mu}^{HVP,LO|_{\tau\;\mathrm{data}}}\,=\,(705.7^{+4.0}_{-4.1})\cdot
10^{-10}$ at $\mathcal{O}(p^{4})$, and
$a_{\mu}^{HVP,LO|_{\tau\;\mathrm{data}}}\,=\,(700.7^{+6.1}_{-5.2})\cdot
10^{-10}$ at $\mathcal{O}(p^{6})$. This reduces the anomaly $\Delta
a_{\mu}\equiv a_{\mu}^{exp}-a_{\mu}^{SM}$ to $2.1$ and $2.3\,\sigma$,
respectively.

##  Appendix A Short-distance constraints

For the parameters contributing to the leading-order chiral LECs [16, 22, 23,
24, 25, 26]:

$\begin{split}F_{V}G_{V}=F^{2},&\qquad F_{V}^{2}-F_{A}^{2}=F^{2},\\\
F_{V}^{2}M_{V}^{2}=F_{A}^{2}M_{A}^{2},&\qquad 4c_{d}c_{m}=F^{2},\\\
8\left(c_{m}^{2}-d_{m}^{2}\right)=F^{2},&\qquad
c_{m}=c_{d}=\sqrt{2}d_{m}=F/2.\end{split}$ (29)

For the even-intrinsic parity sector [17, 58]:

$\begin{split}\lambda_{13}^{P}&=0,\quad\lambda_{17}^{S}=\lambda_{18}^{S}=0,\\\
\lambda_{17}^{A}&=0,\quad\lambda_{21}^{V}=\lambda_{22}^{V}=0.\end{split}$ (30)

The analysis of the $\left\langle VAS\right\rangle$ Green function yields
[18]:

$\begin{split}\kappa_{2}^{S}=\kappa_{14}^{A}=0,\quad\kappa_{4}^{V}&=2\kappa_{15}^{V},\quad\kappa_{6}^{VA}=\frac{F^{2}}{32F_{A}F_{V}},\\\
F_{V}\left(2\kappa_{1}^{SV}+\kappa_{2}^{SV}\right)&=2F_{A}\kappa_{1}^{SA}=\frac{F^{2}}{16\sqrt{2}c_{m}}.\end{split}$
(31)

The study of the $\left\langle VAP\right\rangle$ and $\left\langle
SPP\right\rangle$ Green functions yield the following restrictions on the
resonance couplings [17, 59, 60]:

$\begin{split}\sqrt{2}\lambda_{0}=-4\lambda_{1}^{VA}-\lambda_{2}^{VA}-\frac{\lambda_{4}^{VA}}{2}-\lambda_{5}^{VA}&=\frac{1}{2\sqrt{2}}\left(\lambda^{\prime}+\lambda^{\prime\prime}\right),\\\
\sqrt{2}\lambda^{\prime}=\lambda_{2}^{VA}-\lambda_{3}^{VA}+\frac{\lambda_{4}^{VA}}{2}+\lambda_{5}^{VA}&=\frac{M_{A}}{2M_{V}},\\\
\sqrt{2}\lambda^{\prime\prime}=\lambda_{2}^{VA}-\frac{\lambda_{4}^{VA}}{2}-\lambda_{5}^{VA}&=\frac{M_{A}^{2}-2M_{V}^{2}}{2M_{V}M_{A}},\\\
\lambda_{1}^{PV}=-4\lambda_{2}^{PV}=-\frac{F\sqrt{M_{A}^{2}-M_{V}^{2}}}{4\sqrt{2}d_{m}M_{A}},\quad\lambda_{1}^{PA}&=\frac{F\sqrt{M_{A}^{2}-M_{V}^{2}}}{16\sqrt{2}d_{m}M_{V}}.\end{split}$
(32)

For the odd-intrinsic parity sector [18]:

$\begin{split}\kappa_{14}^{V}&=\frac{N_{C}}{256\sqrt{2}\pi^{2}F_{V}},\quad
2\kappa_{12}^{V}+\kappa_{16}^{V}=-\frac{N_{C}}{32\sqrt{2}\pi^{2}F_{V}},\quad\kappa_{17}^{V}=-\frac{N_{C}}{64\sqrt{2}\pi^{2}F_{V}},\quad\kappa_{5}^{P}=0,\\\
\kappa^{VV}_{2}&=\frac{F^{2}+16\sqrt{2}d_{m}F_{V}\kappa_{3}^{PV}}{32F_{V}^{2}}-\frac{N_{C}M_{V}^{2}}{512\pi^{2}F_{V}^{2}},\quad
8\kappa_{2}^{VV}-\kappa_{3}^{VV}=\frac{F^{2}}{8F_{V}^{2}}.\end{split}$ (33)

##  Appendix B Fit results

Neglecting all the other contributions, we find

$\displaystyle\kappa^{V}_{1}$ $\displaystyle=(-2.1\pm 0.7)\cdot 10^{-2}\text{
GeV}^{-1},$ (34a) $\displaystyle\kappa^{V}_{2}$ $\displaystyle=(-8.8\pm
9.1)\cdot 10^{-3}\text{ GeV}^{-1},$ (34b) $\displaystyle\kappa^{V}_{3}$
$\displaystyle=(2.2\pm 5.8)\cdot 10^{-3}\text{ GeV}^{-1},$ (34c)
$\displaystyle\kappa^{V}_{6}$ $\displaystyle=(-2.1\pm 0.3)\cdot 10^{-2}\text{
GeV}^{-1},$ (34d) $\displaystyle\kappa^{V}_{7}$ $\displaystyle=(1.2\pm
0.5)\cdot 10^{-2}\text{ GeV}^{-1},$ (34e) $\displaystyle\kappa^{V}_{8}$
$\displaystyle=(3.1\pm 0.9)\cdot 10^{-2}\text{ GeV}^{-1},$ (34f)
$\displaystyle\kappa^{V}_{9}$ $\displaystyle=(-0.1\pm 5.9)\cdot 10^{-3}\text{
GeV}^{-1},$ (34g) $\displaystyle\kappa^{V}_{10}$ $\displaystyle=(-5.9\pm
9.6)\cdot 10^{-3}\text{ GeV}^{-1},$ (34h) $\displaystyle\kappa^{V}_{11}$
$\displaystyle=(-3.0\pm 0.6)\cdot 10^{-2}\text{ GeV}^{-1},$ (34i)
$\displaystyle\kappa^{V}_{12}$ $\displaystyle=(1.0\pm 0.8)\cdot 10^{-2}\text{
GeV}^{-1},$ (34j) $\displaystyle\kappa^{V}_{13}$ $\displaystyle=(-5.3\pm
1.1)\cdot 10^{-3}\text{ GeV}^{-1},$ (34k) $\displaystyle\kappa^{V}_{18}$
$\displaystyle=(4.7\pm 0.8)\cdot 10^{-3}\text{ GeV}^{-1}.$ (34l)

These values are in good agreement with our earlier estimation,
$|\kappa_{i}^{V}|\lesssim 0.025\,\mathrm{GeV}^{-1}$ [19].

##  Appendix C Kinematics

For the $\tau^{-}\to\pi^{-}\pi^{0}\gamma\nu_{\tau}$ decays, we have the
following integration region

$\small\begin{split}\mathcal{D}=&\left\\{E_{\gamma}^{min}\leq E_{\gamma}\leq
E_{\gamma}^{max},\,x_{min}\leq x\leq x_{max},\,s_{min}\leq s\leq
s_{max},-1\leq\cos\theta_{-}\leq+1,\,0\leq\phi_{-}\leq
2\pi\right\\},\end{split}$ (35)

with boundaries

$\begin{array}[]{rcl}\frac{(m_{\tau}^{2}-s+x)}{4m_{\tau}^{2}}-\frac{\lambda^{1/2}\left(s,x,m_{\tau}^{2}\right)}{4m_{\tau}}\leq&E_{\gamma}\left(s,x\right)&\leq\frac{(m_{\tau}^{2}-s+x)}{4m_{\tau}}+\frac{\lambda^{1/2}\left(s,x,m_{\tau}^{2}\right)}{4m_{\tau}},\\\
4m_{\pi}^{2}\leq&s\left(x\right)&\leq\left(m_{\tau}-\sqrt{x}\right)^{2},\\\
0\leq&x&\leq\left(m_{\tau}-2m_{\pi}\right)^{2},\\\ \end{array}$ (36)

or interchanging the last two limits,

$\begin{array}[]{rcl}0\leq&x\left(s\right)&\leq\left(m_{\tau}-\sqrt{s}\right)^{2},\\\
4m_{\pi}^{2}\leq&s&\leq m_{\tau}^{2}.\\\ \end{array}$ (37)

There are other ways to write these,

$\begin{array}[]{rcl}4m_{\pi}^{2}\leq&s\left(x,E_{\gamma}\right)&\leq\frac{(m_{\tau}-2E_{\gamma})(2m_{\tau}E_{\gamma}-x)}{2E_{\gamma}}\\\
0\leq&x\left(E_{\gamma}\right)&\leq\frac{2E_{\gamma}(m_{\tau}^{2}-4m_{\pi}^{2}-2m_{\tau}E_{\gamma})}{m_{\tau}-2E_{\gamma}},\\\
E_{\gamma}^{cut}\leq&E_{\gamma}&\leq\frac{m_{\tau}^{2}-4m_{\pi}^{2}}{2m_{\tau}},\\\
\end{array}$ (38)

or exchanging $x\leftrightarrow E_{\gamma}$,

$\begin{array}[]{rcl}\frac{(m_{\tau}^{2}+x-4m_{\pi}^{2})}{4m_{\tau}}-\frac{\lambda^{1/2}\left(x,m_{\tau}^{2},4m_{\pi}^{2}\right)}{4m_{\tau}}\leq&E_{\gamma}\left(s\right)&\leq\frac{(m_{\tau}^{2}+x-4m_{\pi}^{2})}{4m_{\tau}}+\frac{\lambda^{1/2}\left(x,m_{\tau}^{2},4m_{\pi}^{2}\right)}{4m_{\tau}},\\\
0\leq&x&\leq(m_{\tau}-2m_{\pi})^{2},\\\ \end{array}$ (39)

and

$\begin{array}[]{rcl}0\leq&x\left(s,E_{\gamma}\right)&\leq\frac{2E_{\gamma}(m_{\tau}^{2}-s-2E_{\gamma}m_{\tau})}{m_{\tau}-2E_{\gamma}}\\\
4m_{\pi}^{2}\leq&s\left(E_{\gamma}\right)&\leq
m_{\tau}(m_{\tau}-2E_{\gamma}),\\\
E_{\gamma}^{cut}\leq&E_{\gamma}&\leq\frac{m_{\tau}^{2}-4m_{\pi}^{2}}{2m_{\tau}}.\\\
\end{array}$ (40)

Further, interchanging $s\leftrightarrow E_{\gamma}$, we get

$\begin{array}[]{rcl}E_{\gamma}^{cut}\leq&E_{\gamma}\left(s\right)&\leq\frac{m_{\tau}^{2}-s}{2m_{\tau}},\\\
4m_{\pi}^{2}\leq&s&\leq m_{\tau}(m_{\tau}-2E_{\gamma}^{cut}).\\\ \end{array}$
(41)

## References

  * [1] S.L. Glashow, _The renormalizability of vector meson interactions_ , _Nucl. Phys._ 10 (1959) 107.
  * [2] A. Salam and J.C. Ward, _Weak and electromagnetic interactions_ , _Nuovo Cim._ 11 (1959) 568.
  * [3] S. Weinberg, _A Model of Leptons_ , _Phys. Rev. Lett._ 19 (1967) 1264.
  * [4] Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * [5] Muon g-2 collaboration, _Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL_ , _Phys. Rev._ D73 (2006) 072003 [hep-ex/0602035].
  * [6] T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , _Phys. Rept._ 887 (2020) 1 [2006.04822].
  * [7] M. Davier, A. Höcker, B. Malaescu, C.-Z. Yuan and Z. Zhang, _Update of the ALEPH non-strange spectral functions from hadronic $\tau$ decays_, _Eur. Phys. J. C_ 74 (2014) 2803 [1312.1501].
  * [8] J. Miranda and P. Roig, _Effective-field theory analysis of the $\tau^{-}\to\pi^{-}\pi^{0}\nu_{\tau}$ decays_, _JHEP_ 11 (2018) 038 [1806.09547].
  * [9] V. Cirigliano, A. Falkowski, M. González-Alonso and A. Rodríguez-Sánchez, _Hadronic $\tau$ Decays as New Physics Probes in the LHC Era_, _Phys. Rev. Lett._ 122 (2019) 221801 [1809.01161].
  * [10] S. Gonzàlez-Solís, A. Miranda, J. Rendón and P. Roig, _Exclusive hadronic tau decays as probes of non-SM interactions_ , _Phys. Lett. B_ 804 (2020) 135371 [1912.08725].
  * [11] S. Borsanyi et al., _Leading hadronic contribution to the muon magnetic moment from lattice QCD_ , _Nature_ 593 (2021) 51 [2002.12347].
  * [12] J. Bijnens, G. Ecker and J. Gasser, _Radiative semileptonic kaon decays_ , _Nucl. Phys._ B396 (1993) 81 [hep-ph/9209261].
  * [13] V. Cirigliano, G. Ecker and H. Neufeld, _Radiative tau decay and the magnetic moment of the muon_ , _JHEP_ 08 (2002) 002 [hep-ph/0207310].
  * [14] F.E. Low, _Bremsstrahlung of very low-energy quanta in elementary particle collisions_ , _Phys. Rev._ 110 (1958) 974.
  * [15] G. Ecker, J. Gasser, A. Pich and E. de Rafael, _The Role of Resonances in Chiral Perturbation Theory_ , _Nucl. Phys._ B321 (1989) 311.
  * [16] G. Ecker, J. Gasser, H. Leutwyler, A. Pich and E. de Rafael, _Chiral Lagrangians for Massive Spin 1 Fields_ , _Phys. Lett._ B223 (1989) 425.
  * [17] V. Cirigliano, G. Ecker, M. Eidemuller, R. Kaiser, A. Pich and J. Portolés, _Towards a consistent estimate of the chiral low-energy constants_ , _Nucl. Phys._ B753 (2006) 139 [hep-ph/0603205].
  * [18] K. Kampf and J. Novotný, _Resonance saturation in the odd-intrinsic parity sector of low-energy QCD_ , _Phys. Rev._ D84 (2011) 014036 [1104.3137].
  * [19] J.A. Miranda and P. Roig, _New $\tau$-based evaluation of the hadronic contribution to the vacuum polarization piece of the muon anomalous magnetic moment_, _Phys. Rev. D_ 102 (2020) 114017 [2007.11019].
  * [20] J. Wess and B. Zumino, _Consequences of anomalous Ward identities_ , _Phys. Lett._ 37B (1971) 95.
  * [21] E. Witten, _Global Aspects of Current Algebra_ , _Nucl. Phys._ B223 (1983) 422.
  * [22] A. Pich, _Colorless mesons in a polychromatic world_ , in _Phenomenology of large N(c) QCD. Proceedings, Tempe, USA, January 9-11, 2002_ , pp. 239–258, 2002, DOI [hep-ph/0205030].
  * [23] S. Weinberg, _Precise relations between the spectra of vector and axial vector mesons_ , _Phys. Rev. Lett._ 18 (1967) 507.
  * [24] M.F.L. Golterman and S. Peris, _The 7/11 Rule: An Estimate of m(rho) / f(pi)_ , _Phys. Rev._ D61 (2000) 034018 [hep-ph/9908252].
  * [25] M. Jamin, J.A. Oller and A. Pich, _S wave K pi scattering in chiral perturbation theory with resonances_ , _Nucl. Phys._ B587 (2000) 331 [hep-ph/0006045].
  * [26] M. Jamin, J.A. Oller and A. Pich, _Strangeness changing scalar form-factors_ , _Nucl. Phys._ B622 (2002) 279 [hep-ph/0110193].
  * [27] P. Roig and J.J. Sanz Cillero, _Consistent high-energy constraints in the anomalous QCD sector_ , _Phys. Lett._ B733 (2014) 158 [1312.6206].
  * [28] S.-Z. Jiang, Z.-L. Wei, Q.-S. Chen and Q. Wang, _Computation of the $O(p^{6})$ order low-energy constants: An update_, _Phys. Rev. D_ 92 (2015) 025014 [1502.05087].
  * [29] A. Flores-Tlalpa, _Modelo de dominancia de mesones para decaimientos semileptónicos de sabores pesados_ , Ph.D. thesis, CINVESTAV, IPN, 2008.
  * [30] M. Gourdin and E. De Rafael, _Hadronic contributions to the muon g-factor_ , _Nucl. Phys._ B10 (1969) 667.
  * [31] S.J. Brodsky and E. De Rafael, _SUGGESTED BOSON - LEPTON PAIR COUPLINGS AND THE ANOMALOUS MAGNETIC MOMENT OF THE MUON_ , _Phys. Rev._ 168 (1968) 1620.
  * [32] V. Cirigliano, G. Ecker and H. Neufeld, _Isospin violation and the magnetic moment of the muon_ , _Phys. Lett._ B513 (2001) 361 [hep-ph/0104267].
  * [33] A. Sirlin, _Radiative corrections to g(v)/g(mu) in simple extensions of the su(2) x u(1) gauge model_ , _Nucl. Phys. B_ 71 (1974) 29.
  * [34] A. Sirlin, _Current Algebra Formulation of Radiative Corrections in Gauge Theories and the Universality of the Weak Interactions_ , _Rev. Mod. Phys._ 50 (1978) 573.
  * [35] A. Sirlin, _Large m(W), m(Z) Behavior of the O(alpha) Corrections to Semileptonic Processes Mediated by W_ , _Nucl. Phys. B_ 196 (1982) 83.
  * [36] W. Marciano and A. Sirlin, _Radiative Corrections to beta Decay and the Possibility of a Fourth Generation_ , _Phys. Rev. Lett._ 56 (1986) 22.
  * [37] W. Marciano and A. Sirlin, _Electroweak Radiative Corrections to tau Decay_ , _Phys. Rev. Lett._ 61 (1988) 1815.
  * [38] E. Braaten and C.-S. Li, _Electroweak radiative corrections to the semihadronic decay rate of the tau lepton_ , _Phys. Rev. D_ 42 (1990) 3888.
  * [39] W.J. Marciano and A. Sirlin, _Radiative corrections to pi(lepton 2) decays_ , _Phys. Rev. Lett._ 71 (1993) 3629.
  * [40] J. Erler, _Electroweak radiative corrections to semileptonic tau decays_ , _Rev. Mex. Fis._ 50 (2004) 200 [hep-ph/0211345].
  * [41] J.S. Schwinger, _PARTICLES, SOURCES, AND FIELDS. VOL. 3_ (1989).
  * [42] M. Drees and K.-i. Hikasa, _Scalar top production in e+ e- annihilation_ , _Phys. Lett. B_ 252 (1990) 127.
  * [43] D. Gómez Dumm and P. Roig, _Dispersive representation of the pion vector form factor in $\tau\to\pi\pi\nu_{\tau}$ decays_, _Eur. Phys. J. C_ 73 (2013) 2528 [1301.6973].
  * [44] F. Guerrero and A. Pich, _Effective field theory description of the pion form-factor_ , _Phys. Lett. B_ 412 (1997) 382 [hep-ph/9707347].
  * [45] Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * [46] BABAR collaboration, _Precise Measurement of the $e^{+}e^{-}\to\pi^{+}\pi^{-}(\gamma)$ Cross Section with the Initial-State Radiation Method at BABAR_, _Phys. Rev._ D86 (2012) 032013 [1205.2228].
  * [47] Belle collaboration, _High-Statistics Study of the tau- — > pi- pi0 nu(tau) Decay_, _Phys. Rev. D_ 78 (2008) 072006 [0805.3773].
  * [48] ALEPH collaboration, _Branching ratios and spectral functions of tau decays: Final ALEPH measurements and physics implications_ , _Phys. Rept._ 421 (2005) 191 [hep-ex/0506072].
  * [49] KLOE collaboration, _Precision measurement of $\sigma(e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}\gamma)/\sigma(e^{+}e^{-}\rightarrow\mu^{+}\mu^{-}\gamma)$ and determination of the $\pi^{+}\pi^{-}$ contribution to the muon anomaly with the KLOE detector_, _Phys. Lett._ B720 (2013) 336 [1212.4524].
  * [50] BESIII collaboration, _Measurement of the $e^{+}e^{-}\to\pi^{+}\pi^{-}$ cross section between 600 and 900 MeV using initial state radiation_, _Phys. Lett. B_ 753 (2016) 629 [1507.08188].
  * [51] M. Benayoun, P. David, L. DelBuono and F. Jegerlehner, _Upgraded Breaking Of The HLS Model: A Full Solution to the $\tau^{-}e^{+}e^{-}$ and $\phi$ Decay Issues And Its Consequences On g-2 VMD Estimates_, _Eur. Phys. J. C_ 72 (2012) 1848 [1106.1315].
  * [52] M. Benayoun, P. David, L. DelBuono and F. Jegerlehner, _An Update of the HLS Estimate of the Muon g-2_ , _Eur. Phys. J. C_ 73 (2013) 2453 [1210.7184].
  * [53] M. Benayoun, P. David, L. DelBuono and F. Jegerlehner, _Muon $g-2$ estimates: can one trust effective Lagrangians and global fits?_, _Eur. Phys. J. C_ 75 (2015) 613 [1507.02943].
  * [54] M. Benayoun, L. Delbuono and F. Jegerlehner, _BHLS 2, a New Breaking of the HLS Model and its Phenomenology_, _Eur. Phys. J. C_ 80 (2020) 81 [1903.11034].
  * [55] CLEO collaboration, _Hadronic structure in the decay tau- — > pi- pi0 neutrino(tau)_, _Phys. Rev. D_ 61 (2000) 112002 [hep-ex/9910046].
  * [56] OPAL collaboration, _Measurement of the strong coupling constant alpha(s) and the vector and axial vector spectral functions in hadronic tau decays_ , _Eur. Phys. J. C_ 7 (1999) 571 [hep-ex/9808019].
  * [57] M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, _A new evaluation of the hadronic vacuum polarisation contributions to the muon anomalous magnetic moment and to $\mathbf{\boldsymbol{\alpha}(m_{Z}^{2})}$_, _Eur. Phys. J._ C80 (2020) 241 [1908.00921].
  * [58] A. Guevara, P. Roig and J.J. Sanz-Cillero, _Pseudoscalar pole light-by-light contributions to the muon $(g-2)$ in Resonance Chiral Theory_, _JHEP_ 06 (2018) 160 [1803.08099].
  * [59] V. Cirigliano, G. Ecker, M. Eidemuller, A. Pich and J. Portolés, _The < VAP > Green function in the resonance region_, _Phys. Lett._ B596 (2004) 96 [hep-ph/0404004].
  * [60] V. Cirigliano, G. Ecker, M. Eidemuller, R. Kaiser, A. Pich and J. Portolés, _The < SPP > Green function and SU(3) breaking in K(l3) decays_, _JHEP_ 04 (2005) 006 [hep-ph/0503108].


11institutetext: Colorado State University, Fort Collins CO 80523, USA
11email<EMAIL_ADDRESS>22institutetext: University of Tennessee,
Chattanooga TN 37403 USA Z
22email<EMAIL_ADDRESS>33institutetext: Colorado State University,
Fort Collins CO 80523, USA
33email<EMAIL_ADDRESS>
# A Multi-Domain Multi-Task Approach for Feature Selection from Bulk RNA
Datasets

Karim Salta 11 0009-0002-4589-3335 Tomojit Ghosh 22 0000-0001-6875-3840
Michael Kirby 33 0000-0001-9802-9263

###### Abstract

In this paper a multi-domain multi-task algorithm for feature selection in
bulk RNAseq data is proposed. Two datasets are investigated arising from mouse
host immune response to Salmonella infection. Data is collected from several
strains of collaborative cross mice. Samples from the spleen and liver serve
as the two domains. Several machine learning experiments are conducted and the
small subset of discriminative across domains features have been extracted in
each case. The algorithm proves viable and underlines the benefits of across
domain feature selection by extracting new subset of discriminative features
which couldn’t be extracted only by one-domain approach.

###### Keywords:

Sparse Feature Selection Multi-Domain Multi-Task Learning Bulk RNA VAE HPC.

##  1 Introduction

In the field of bioinformatics, researchers often use microarray or next-
generation sequencing techniques to study the expression levels of genes, with
each sample typically having tens of thousands of features. The large number
of features often necessitates the use of feature selection algorithms to
improve the performance of machine learning tasks such as classification given
that many of the observed features may be unrelated to the biological
phenomenon of interest. In this way, feature selection algorithms can be used
to determine the processes related to a biological mechanism, e.g., the host
immune response to infection. Other downstream benefits of feature selection
include data visualization and understanding, reduced storage requirements,
and faster computations.

Multi-domain feature extraction addresses the problem of leveraging data from
disparate sources and is related to the more general problem of multi-domain
learning (MDL) [27]. In this paper we address a special case of MDL that
involves the classification of data related to the host immune response to
infection. The host consists of multiple lines of the Collaborative Cross
mouse, the pathogen under consideration is Salmonella. The two domains consist
of gene expression data collected from liver and spleen tissues. The proposed
machine learning approach has the potential to identify novel biomarkers whose
signals are are too weak to be captured by analyzing domains individually. The
methodology will be demonstrated by selecting potentially important biomarkers
that appear to be amplified in strength by the multi-domain data synthesis for
characterizing biological processes that exist simultaneously in different
tissues.

In a variety of existing methods the designs vary from shallow to deep, with
the networks optimized for regression, classification, dimensionality
reduction or a combination of multiple tasks, i.e. multi-task learning (MTL)
[26]. Most of the feature selection methods fall into the category of MLT
methods, with the the objective function considering the combination of
different goals resulting in better generalization of results. However, the
majority of them focus on a single domain feature selection. In this this
paper we suggest a new MDL method with a multi-task objective (MDL/MTL)
function, or, in terms of [28], a multi-domain multi-task (MDMT) method. The
MDMT methods have been used widely in Natural Language Processing applications
but much less for the analysis of biological data sets.

The feature selection task is frequently performed by the introduction of
$l_{p}$-norms with $p=0,1,2$, or a combination of norms of weights at
intermediate layers, possibly stacked closer to the input. The $l_{0}$ is of
course appealing but the most expensive to compute [7, 8]. However, in [24] it
was shown that $l_{1}$-norm based methods can be quite competitive when
applied to biological data, while not having the complications of dealing with
$l_{0}$-norm. When applied to biological task, the $l_{0}$-norm based methods
may require a priori knowledge of the size of a subset of biologically
significant features, which is certainly not the case for many explorative
tasks. Moreover, when the domains are very different the subsets of important
features can be very different as well, hence it becomes even harder to
approximate the number of selected features. At the same time, $l_{0}$-norm
based methods haven’t been battle-tested across domain as much as $l_{1}$-norm
based methods, intuitively, the discrete distribution can negatively effect
the alignment of disparate domains. Considering all the pros and cons, in this
paper we opt to employ $\ell_{1}$-norm based sparsity promotion.

The contributions of this paper include the following: we propose a new
sparsity promoting MDMT architecture for feature selection; this approach uses
a new masking term that restricts the features that contribute to the cost
function; we demonstrate the utility of developed algorithm on gene expression
dataset including liver and spleen domains; we conclude that our MDMT approach
allows to find new features that are significantly discriminative only across
two domains, i.e., are identified when the data is restricted to a single
domain. The promise of this approach is an enriched picture of the host immune
response that has the potential to lead to a better understanding of the
biological process across tissues.

The organization of this paper is as follows: In Related Work, we present a
review of the related articles, situating our study within the broader context
of the field and highlighting key contributions from prior research.
Methodology details the method we employed, outlining the design of a neural
network behind our chosen approach and the techniques utilized for design. In
Data, we describe the data used in our investigation, elaborating on pre-
processing steps and labelled groups. Experiment delves into the computational
experimental setup, discussing the training process, and criteria for rough-
tuning. In Results, we present the results of our experiments, providing an
interpretation and discussion of our findings. Finally, in Conclusions, we
summarize the highlights of the research and contributions, the implications
of our results, and potential avenues for future work in this domain.

##  2 Related Work

The majority of the feature selection methods study one domain. In what
follows we survey a variety of different techniques and algorithms, including
linear, non-linear methods and the methods exploiting neural networks. In one
of the most influential papers [14], feature selection is cast as a regression
task with $l_{1}$ regularization of the norm of discriminative vector. This
has become a common approach, see also, [17, 18, 16, 15]. The Lasso-type
methods fail to capture nonlinear interactions between the features. The non-
linear methods developed as the kernelized modifications of Lasso method
showed decent efficiency when applied to biological data [19, 20, 21]. Note
that there are also some other methods aiming to sparsify a signal in latent
dimensions [22, 23].

Deep Feature Selection DFS [7] is one of the first deep neural network
algorithms designed specifically for feature selection. DFS employs a one-to-
one sparsity layer at the input. The weights on these single connections are
penalized with a $\ell_{1}$-norm minimization of norm of weights of this layer
in a spirit of [15]; the resulting non-zero weights in the sparsity layer
correspond to the selected features. In [3] autoencoders are proposed for
feature selection with a sparsity layer used in a fashion similar to DFS. Now
the $\ell_{1}$-norm of the weights of sparsity layer is minimized jointly with
the reconstruction error. Concrete Autoencoders (CAE) [4] use a concrete
selector layer as the first layer in autoencoder setting based on continuous
relaxations of concrete random variables suggested in [5]. A supervised CAE
method reported in [6] was apparently susceptible to overfitting with limited
data. The FsNet paper [6] addressed this problem by introducing small weight-
predictor networks. In terms of design [6] is one of the closest to the design
developed in this paper, i.e., our method is based on two neural networks:
autoencoder and classifier in latent space, but the approach to sparsification
is different, and, most importantly, the method in [6] is designed for one
domain and the implication of extracted features is quite different.

Most of the feature selection methods can be grouped into three broad
categories: filter, wrapper, and embedded methods. In filter methods the
features are typically scored, ranked and thresholded with respect to some
classification task using different measures such as correlation and mutual
information [1]. Filter methods can be very fast, but the quality of extracted
features is poor in terms of robustness and adaptability to different
datasets. The wrapper methods [2] are universal methods used on top of any
learning algorithm based on practical heuristic search of a subset of $d$
features in $2^{d}$ space providing the better performance for the underlying
algorithm. They are universal and capable of obtaining great results given the
large number of samples. For small datasets in high-dimensional space they
tend to overfit, and the NP-hardness of the problem makes the computations
prohibitively expensive. In embedded methods the feature selection process is
typically performed concurrently with some learning algorithm. For example,
Iterative Feature Removal (IFR) uses the absolute weights of a sparse SVM
model as a criterion for selecting features from a high-dimensional biological
data set [25]. Our paper along with the most related works falls into the
embedded method category.

##  3 Methodology

The methods described above only address data residing in one domain. The
major question that motivated this study was what are the biological features
in datasets sampled from different domains that appear to be related to the
host immune response to infection only when studied across domains, naturally
leading the consideration of the domain alignment task along with the feature
selection. With that said and with a general design of the network in mind
we’ve been looking for the most capable in terms of domain alignment method in
application to RNA data. This led us to [30], which is an MDMT method based on
a pair of domain specific variational autoencoders (VAE’s) [31] generating
aligned embeddings for datasets of very different modalities (singe-cell RNA
and Chromatin images), and we adapt this method now enhanced with sparsity
promoting optimization constraints for feature selection. In order to find a
universal representation across tasks, the MTL methods in deep neural networks
[29] either improve the architecture of neural networks, or try to find a
balance between concurrently trained objectives. This paper benefits from both
since our network has shared subnets and at the same time we roughly fine-tune
the coefficients used in [30] along with the contribution of sparsity
promoting loss function. Utilizing both methods is also justified by the
results of generalization of unbalanced optimization methods, e.g. [32], [33],
[34], indicating that overall they don’t outperform the naive approach when
all loss functions are weighted with constant scalars.

Our proposed method is based on the network shown in Figure 1, implemented in
PyTorch and trained with the AdamW optimizer. The cost function is a weighted
combination of objective functions associated with the reconstruction,
classification and sparsification tasks. In our settings we can observe that
during the early stages of training the algorithm is learning the shared
representations of different domains such that similar samples group together
regardless of domain of origin. At the later stages, the sparsification goal
is becoming more important with a relatively higher contribution to the total
objective and this behaviour continues until the conflicting tasks reach the
balance and no further sparsification is possible without a significant loss
in classification. We run the training process multiple times. In a spirit of
embedded methods, we treat the resulting magnitude of the weights of sparse
layer as indicating the importance of different features. However, for the
post-processing we employ the frequency of selected features across all runs
as it appears to be a more robust metric for feature importance.

It was mentioned before that our design is developed based on the network
suggested in [30] with modifications. The classication task is performed in
latent space as before, but the inputs of domain-specific VAE’s are sparsified
by the shared Sparsification Layer (SL) and, naturally, the VAE’s are trained
to reconstruct only these sparsified inputs. Note that in [30], a primary goal
is the alignment of the domains in the latent space coupled with the
reconstruction of hyper-dimensional RNA data. In contrast to [30], our main
goal is the across-domains classification with the sparsification of inputs.
The choice of variational modification of autoencoders was dictated by the
distribution of latent space provided by this particular modification,
allowing further indirect "easy" and relaxed alignment through the shared
classifier. Hence, we not only train VAE’s solely to reconstruct the
sparsified signal, but also we omit the loss function minimizing the KL-
divergence across domains from the original design. The resulting network
consist of 4 subnets: shared between domains Sparse Layer (SL) and Classifier
subnets, and two domain specific Variational AutoEncoders (VAE1 and VAE2), as
depicted in Figure 1.

The SL is a one-to-one mapping: $x\rightarrow\mathbf{W}\odot x$, with the
$\ell_{1}$-norm of weights $\left\|\mathbf{W}\right\|_{1}$ penalty used to
promote sparsity.

The VAE1 and VAE2 subnets are deep fully-connected networks with the following
specifications:

  * •

Encoder: 2 linear layers with 1024 nodes with batch normalization and Relu-
activations, followed by 1 linear mapping to $\mu$ and $\sigma$ living in
128-dimensional space

  * •

Decoder: 2 linear layers with 1024 nodes with batch normalization and Relu-
activations, followed by 1 layer mapping to input space

Each Encoder performs the mapping of the sparsified input to
$\mathbf{R}^{128}$: $\mathbf{W}\odot x\rightarrow\mu,\,\sigma$, while the
Decoder maps distribution in the latent space into the input space:
$\mu+\sigma\rightarrow\tilde{x}\in\mathbf{R}^{34861}$, and the output is
further masked by the frozen weights of the sparse layer
$\mathbf{W}^{*}\odot\tilde{x}$ and fed to the MSE loss function of respective
VAE.

The Classifier is the subnet consisting of 5 linear layers with 1024 nodes
with Relu-activations, followed by 1 layer mapping to 2-dimensional space and
1 layer mapping to 1-dimensional space followed by a standard sigmoidal
activation function. It is trained to classify embeddings of inputs from both
domains in their respective latent spaces.

Figure 1: Network design

The overall objective function is a weighted sum of objective functions for
different tasks including reconstruction, normalization, classification and
sparsification, namely,

$\displaystyle\operatorname{Loss}\left(x_{1},x_{2}\right)=$
$\displaystyle\alpha\cdot\left(\operatorname{Loss}_{\text{rec
}}\left(x_{1}\right)+\operatorname{Loss}_{\text{rec
}}\left(x_{2}\right)\right)+$
$\displaystyle\beta\cdot\left(\operatorname{Loss}_{\text{var
}}\left(x_{1}\right)+\operatorname{Loss}_{\text{var
}}\left(x_{2}\right)\right)+$
$\displaystyle\gamma\cdot\left(\operatorname{Loss}_{\text{class
}}\left(x_{1}\right)+\operatorname{Loss}_{\text{class
}}\left(x_{2}\right)\right)+$
$\displaystyle\,\theta\cdot\operatorname{Loss}_{\text{sparse }}$ (1)

where

$\displaystyle\operatorname{Loss}_{\text{rec }}(x)$
$\displaystyle=MSE(\mathbf{W}\odot x,\mathbf{W}^{*}\odot\tilde{x})$
$\displaystyle\operatorname{Loss}_{var}(x)$
$\displaystyle=D_{KL}\left(N\left(\mu,\sigma^{2}I\right),N(0,I)\right)$
$\displaystyle\operatorname{Loss}_{\text{class }}(x)$
$\displaystyle=\log\operatorname{Loss}\left(p(x),\tilde{p}(x)\right)$
$\displaystyle\operatorname{Loss}_{\text{sparse }}$
$\displaystyle=\left\|\mathbf{W}\right\|_{1}$ (2)

The block-scheme of one training process is given by Algorithm 1.

Input : $\theta,\alpha,\beta,\gamma$ | weight of loss functions
---|---
$argsOpt1$, $argsOpt2$ | parameters of Adam optimizers for VAE’s
$argsOptCls$ | parameters of Adam optimizers for Classifier
$argsOptSL$ | parameters of Adam optimizers for SL

Output : $\mathbf{W}$ \- weights of SL
---

Data: $[x_{1},x_{2}]$ | list of pairs of equally seized batches sampled from both datasets
---|---
| 20 000 elements, i.e epochs of training

Initialize : SL, VAE1, VAE2, Classifier (initialize subnets)
---
OptSL = AdamW($argsSL$)
Opt1 = AdamW($argsOpt1$)
Opt2 = AdamW($argsOpt2$)
OptCls = AdamW($argsCls$)

1

2for  _$x_{1},x_{2}\in[x_{1},x_{2}]$_ do

3 $sp(x_{1}),sp(x_{2})=\mathbf{W}\odot x_{1},\mathbf{W}\odot x_{2}$

4 $\tilde{x}_{1},\mu_{1}+\epsilon\sigma_{1}=\text{VAE}1(sp(x_{1}))$

5 $\tilde{x}_{2},\mu_{2}+\epsilon\sigma_{2}=\text{VAE}2(sp(x_{2}))$

6 $\tilde{p}(x_{1})=\text{Classifier}(\mu_{1}+\epsilon\sigma_{1})$

7 $\tilde{p}(x_{2})=\text{Classifier}(\mu_{2}+\epsilon\sigma_{2})$

8 Calculate Loss as in (1)

9 Backpropagate Loss

10 Step all optimizers

11

Algorithm 1 Training

After the rough-tuning of the hyper-parameters from [30] along with the
sparsity contribution and parameters for optimizers, we set the following
hyper-parameters:

  * •

$\alpha,\,\beta,\,\gamma,\,\theta$ = 10, $10^{-4}$, 1, $10^{-4}$

  * •

LR’s for AdamW optimizer for VAE’s, SL and Classifier = $10^{-4}$

  * •

all other parameters for AdamW optimizer are set to Pytorch default

##  4 Data

The data consists of mice bulk RNA sequences extracted from two different
tissues: spleen and liver, used as two different domains for the purposes of
this paper. The mice that were exposed to Salmonella infection were monitored
and categorized by health status as tolerant, resistant, susceptible, or
delayed susceptible, the latter two being related to strains unifying the mice
who died within 1 or 3 weeks respectively. In all our experiments we combine
these latter two groups into one susceptible group. With this new labelling we
have 31 and 9 tolerant samples, 27 and 7 resistant samples and 90 and 53
susceptible samples for spleen and liver domains respectively for all infected
mice. Also, the data includes control samples representing the mice who had
never been exposed to infection labeled as "never infected". This group
accounts for 104 samples, with 93 samples from spleen and 11 samples from
liver. Phenotypes of these samples are determined based on their genetic
strains. Initial bulk RNA dataset was TMM-normalized, the outliers and
duplicates have been detected and dropped out. Finally, the domain-specific
data, i.e. combined RNA data for samples from spleen and combined RNA data for
samples from liver, have been z-scored for each domain separately and filtered
for common across tissues genes in all feature selection algorithms resulting
in data samples consisting of 34,861 genes, i.e., the dimension of the input
space.

##  5 Experiment

We consider three distinct types of experiment. In the first type the goal is
to extract a small subset of features that discriminate among the phenotypes
susceptible and tolerant. The second type extracts features discriminating
among the phenotypes susceptible and resistant. The third type extracts the
features discriminating between infected and never infected mice. The exact
Python code with the training models and post-processing utilities is
available at [35]. With a use of Ray package [36] the experiment was run on 16
V100 GPU’s in a multiprocessing mode for 10 different random samplings of
85$\%$ of data, with 90 different weights initializations for each, summing up
to 900 runs. You can see the typical evolution of training process including
60,000 epochs for all three experiments in Figures 2, 3, 4.

Figure 2: Weighted components of overall losses for phenotypes tolerant versus
susceptible across domain experiment: (a) - reconstruction errors, (b) -
classification errors, (c) - sparsity loss, (d) - overall loss.

Figure 3: Weighted components of overall losses for phenotypes resistant
versus susceptible across domain experiments: (a) reconstruction error, (b)
classification error, (c) sparsity loss, (d) total loss.

Figure 4: Weighted components of overall losses for infected mice versus never
infected mice across domain experiment: (a) reconstruction error, (b)
classification error, (c) sparsity loss, (d) total loss.

Figure 5: Latent space PCA’s for all three across domain experiments: (a)
phenotypes tolerant (TOL) versus susceptible (SUS), (b) phenotypes resistant
(RES) versus susceptible (SUS), (c) infected (INF) mice versus never infected
mice (NOT INF).

The number of epochs was fixed at 20,000 based on the indication of the
flattening of the sparsity curves as long as it doesn’t effect the accuracy of
classification and the reconstruction loss is relatively small. Again, the
reconstruction is not in the primary focus of this method, i.e, it wasn’t the
major task to consider when deciding on the number of epochs, especially,
since the required for the across-domain classification alignment in the
latent space was typically achieved even after 20,000 epochs, see the Figure
5. The visualization of loss indicates that further sparsification slightly
decreases the classification accuracy for two out of three experiments, by
limiting the number of epochs we also prevent this long-run negative effect.
The PCA images for all figures represent the PCA of combined representations
of both domains in their bottleneck layers, and indicate a proper clustering
and separability across domains, i.e., good alignment in latent space.

##  6 Results

The post-processing of the sparsity layer weights was conducted in a same way
for all 3 experiments. Firstly, all the weights across 900 runs have been
aggregated, normalized and the Elbow Method was applied to find a threshold,
and later for each run the weights below the threshold were set to zero. At
the next step we calculated the frequencies of features appearing in subsets
of features with non-zero weights across all runs. The resulting distributions
of frequencies are shown in Figure 6, along with resulting number of features
selected in two consecutive steps by Elbow Method.

Figure 6: Results from the across domain (MDMT) experiments. The features
selected are ordered by frequency for 3 experiments: panels (a) and (d)
correspond to the phenotypes tolerant versus susceptible; panels (b) and (e)
correspond to resistant versus susceptible; panels (c) and (f) correspond to
infected versus never infected.

Importantly, in addition to the cross domain learning, we also selected
features for each domain separately for all three experiments. This allows us
to evaluate the distinct characteristics of single domain and multi-domain
alignment for feature extraction. In Figure 7 panels (a), (b), and (c), we can
see the distribution of features for all three across domain experiments
grouped by overlapping with features selected in one-domain experiments.

Figure 7: Grouped by overlapping distribution of features for all experiments.
Columns are associated with experiments, i.e. 1st - tolerant versus
susceptible, 2nd - resistant versus susceptible, 3rd - infected versus never
infected. (a),(b),(c) distributions of features extracted in across domains
experiments; (d),(e),(f) features extracted in separate spleen experiment,
features are grouped by overlapping with features extracted in respective
across domains experiment: "Both" - overlapping, "Spleen" - not overlapping;
(g),(h),(i) features extracted in separate liver experiment, features are
grouped by overlapping with features extracted in respective across domains
experiment: "Both" - overlapping, "Liver" - not overlapping.

"Both" features is a subset of these features that also appear in both the
spleen and liver domains. "None" refers to the features that are the features
that only appear in across domain results. "Spleen" are the features appearing
in both across domains and spleen domain results but not in the liver domain
results. Finally, "Liver" features are the features appearing in both across
domains and liver domain results but not in the spleen domain results. We can
see that for all the experiments apart from the features captured from domain-
specific experiments some new highly-weighted features were captured in across
domain experiments (the blue line denoting "None"). These features that are
only present in the across domain experiment reflect the new information being
captured by the proposed method. These correspond to biomarkers that we
suspect have a potentially unique role in the host response to infection.

Figures 7, panels (d),(g),(e),(h),(f), and (i) show the distributions of
features extracted in one domain experiments grouped by overlapping with
respective across domains results. Even though the majority of the most
discriminative features extracted in separate experiments are also extracted
in across domain experiments, some highly-weighted features extracted from one
domain experiments are apparently not captured in respective across domains
experiments.These results might indicate that the improvement in robustness is
needed, but at the same time they may contain biological insights about the
across domain importance of some features that couldn’t be obtained from
studying only one domain for Salmonella infection or about the difference in
manifestations of infection in different tissues. Based on the results from
additional experiments with other datasets most likely the latter is true.

Returning to Figure 5, we see the results of the embedding of the sparsified
input in a two-dimensional latent space using the features found in the across
domain architecture for all three experiments. We observe that the tolerant
and resistant versus susceptible experiments give excellent classification on
test data in the latent space. In contrast, infected versus control mice are
not as easily discriminated.

##  7 Conclusions

This paper proposes a novel architecture for multi-domain, multi-task feature
selection. This area of research has a relatively small literature, possibly
because of the complexity associated with simultaneously exploring multiple
incommensurate measurement domains. The proposed approach leverages prior art
in multi-domain learning while adding a masked feature selection approach that
serves to identify biologically relevant aspects of the host immune response
to infection. We demonstrate that the demands of the MDMT problem formulation
can be successfully addressed with the proposed architecture. Further, the
application of the approach leads to the discovery of novel biomarkers whose
signals appear to be amplified by the multi-domain approach; indeed, a
fraction of these biomarkers do not appear in either single experiment. Hence,
this approach holds the promise of generating new biological insights that
might go undetected using single domain methodologies. Additionally, we
observed the the MDMT features provide excellent classification results
between susceptible and tolerant or resistant phenotypes.

There are several possible modifications of our method. In this paper the
optimization problem descends in the direction of the gradient of weighted
objectives. There is growing evidence that isolating descent directions to
improve individual cost functions may lead to improved solutions.
Additionally, alternative data reduction and reconstruction mappings could be
explored reflecting recent developments in deep neural networks including
graph convolutional neural networks or transformers. This preliminary work is
focused on algorithm development; we propose to explore the biological
ramifications of the biomarkers in future work.

We believe further development in these directions could lead to even more
efficient methods and provide biologists with a new perspective on
understanding the evolution of infections in tissues.

Acknowledgements: We would like to thank Helene Andrews-Polymenis and David
Threadgill for providing the data for this investigation.

## References

  * [1] Blum, A. L., Langley, P.: Selection of relevant features and examples in machine learning, Artif. Intell., 97, 1, pp. 245–271 (1997)
  * [2] John, G.H., Kohavi, R., Pfleger, K.: Irrelevant features and the subset selection problem, Proc. 11th Int. Conf. Mach. Learn., pp. 121–129 (1994)
  * [3] Han, K., Wang, Y., Zhang, C., Li, Ch., Xu, C.: AutoEncoder Inspired Unsupervised Feature Selection, ArXiv, 1710.08310, https://arxiv.org/abs/1710.08310 (2017)
  * [4] Fatih Balin, M., Abid, A., Zou, J.Y.: Concrete autoencoders: Differentiable feature selection and reconstruction, ICML (2019)
  * [5] Maddison, C. J., Mnih, A., Teh, Y. W.: The Concrete Distribution: A Continuous Relaxation of Discrete Random Variables, ArXiv, 1611.00712, https://arxiv.org/abs/1611.00712 (2016)
  * [6] Singh, D., E., Yamada, M.: Fsnet: Feature selection network on high-dimensional biological data. ArXiv, 2001.08322, https://arxiv.org/abs/2001.08322, (2020)
  * [7] Li, Y., Chen, C., Wasserman, W.: Deep feature selection: theory and application to identify enhancers and promoters. Journal of Computational Biology, 23(5), pp. 322–-336 (2016)
  * [8] Feng J., Simon, N.: Sparse-Input Neural Networks for High-dimensional Nonparametric Regression and Classification, ArXiv, 1711.07592, https://arxiv.org/abs/1711.07592 (2017)
  * [9] Hinton, G., Salakhutdinov, R.: Reducing the dimensionality of data with neural networks. Science, 313, pp. 504–507 (2006)
  * [10] Salah, R., Pascal, V., Xavier, M., Xavier, G, Yoshua, B.: Contractive Auto-Encoders: Explicit Invariance During Feature Extraction, Proceedings of the 28th International Conference on Machine Learning, ICML ’11, pp. 833–840, New York (2011)
  * [11] Hinton, G., Osindero, S., Teh, Y.: A fast learning algorithm for deep belief nets. Neural Computation, 18, pp. 1527–1554 (2006)
  * [12] Yamada, Y., Lindenbaum, O., Negahban, S., Kluger., Y.: Feature selection using stochastic gates. International Conference on Machine Learning, PMLR, 119, pp. 10648–10659, (2020)
  * [13] Lemhadri, I., Ruan, F., Abraham, L., Tibshirani, R.: Lassonet: A neural network with feature sparsity. Journal of Machine Learning Research, 22(127), pp. 1–29, (2021)
  * [14] Tibshirani, R.: Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, 58, pp. 267–288, (1996)
  * [15] Zou, H., Hastie, T.: Regularization and variable selection via the elastic net., J. R. Stat. Soc. Series B Stat. Methodol. 67(2), pp. 301–320 (2005)
  * [16] Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., Knight, K.: Sparsity and smoothness via the fused lasso. J. Roy. Stat. Soc. B (Stat.Methodol.), 67(1), pp. 91–108 (2005)
  * [17] Zou, H.: The adaptive lasso and its oracle properties. J. Amer. Stat. Assoc., 101(476), pp. 1418–-1429 (2006)
  * [18] Meinshausen, N.: Relaxed lasso. Comput. Stat. Data Anal., 52(1), pp. 374–-393 (2007)
  * [19] Yamada, M., Jitkrittum, W., Sigal, L., Xing, E.P., Sugiyama, M.: High-dimensional feature selection by feature-wise kernelized lasso, Neural computation, 26(1), pp. 185-–207, (2014)
  * [20] Liu, H., Wasserman, L., Lafferty, J. D.: Nonparametric regression and classification with joint sparsity constraints, Proc. Adv. Neural Inf. Process. Syst., pp. 969–976 (2009)
  * [21] Shevade S. K., Keerthi, S. S.: A simple and efficient algorithm for gene selection using sparse logistic regression, Bioinformatics, 19, 17, pp. 2246–2253 (2003)
  * [22] Chan, A.B., Vasconcelos, N., Lanckriet, G.R.G.: Direct convex relaxations of sparse SVM. In: International Conference on Machine Learning, (2007)
  * [23] Gurram, P., Kwon, H.: Optimal sparse kernel learning in the empirical kernel feature space for hyperspectral classification. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, (2014)
  * [24] Ghosh, T., Karimov, K. ,Kirby, M.: Sparse Linear Centroid-Encoder: A Biomarker Selection tool for High Dimensional Biological Data. IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pp. 3012–3019 (2023)
  * [25] O’Hara, S., Wang, K., Slayden, R.A. et all.: Iterative feature removal yields highly discriminative pathways. BMC genomics, 14(1), pp. 1–15, (2013)
  * [26] Caruana, R.: Multitask learning, Machine learning, 28(1), pp. 41–75 (1997)
  * [27] Joshi, M., Cohen, W. W., Dredze M., Rosé, C.P.: Multi-domain learning: when do domains matter?, Proceedings of the 2012 Joint Conference on Empirical Methods in Natural Language Processing and Computational Natural Language Learning, pp. 1302–1312 (2012)
  * [28] Yang, Y., Hospedales, T.M.: A unified perspective on multi-domain and multi-task learning. ArXiv, 1412.7489, https://arxiv.org/abs/1412.7489 (2014)
  * [29] Ruder, S., An Overview of Multi-Task Learning in Deep Neural Networks. ArXiv, 1706.05098, https://arxiv.org/abs/1706.05098 (2017)
  * [30] Yang, K.D., Belyaeva, A., Venkatachalapathy, S. et al.: Multi-domain translation between single-cell imaging and sequencing data using autoencoders. Nat Commun 12(31) (2021). 10.1038/s41467-020-20249-2
  * [31] Kingma, D.P., Welling, M.: Auto-Encoding Variational Bayes. ArXiv, 1312.6114, https://arxiv.org/abs/1312.6114 (2013)
  * [32] Chen, Z., Badrinarayanan, V., Lee, C.Y., Rabinovich, A.: GradNorm: Gradient Normalization for Adaptive Loss Balancing in Deep Multitask Networks. International Conference on Machine Learning (2017)
  * [33] Sener, O., Koltun, V.: Multi-Task Learning as Multi-Objective Optimization. Neural Information Processing Systems (2018)
  * [34] Guo, M., Haque, A., Huang, DA., Yeung, S., Fei-Fei, L.: Dynamic Task Prioritization for Multitask Learning. Computer Vision – ECCV 2018, Lecture Notes in Computer Science, 11220, Springer, Cham, 10.1007/978-3-030-01270-0_17 (2018)
  * [35] https://github.com/kkarimov/iccs2024
  * [36] https://github.com/ray-project/ray