# Sequence-based Machine Learning Models in Jet Physics

Rafael Teixeira de Lima
SLAC National Accelerator Laboratory
<EMAIL_ADDRESS>
###### Abstract

Sequence-based modeling broadly refers to algorithms that act on data that is
represented as an ordered set of input elements. In particular, Machine
Learning algorithms with sequences as inputs have seen successfull
applications to important problems, such as Natural Language Processing (NLP)
and speech signal modeling. The usage this class of models in collider physics
leverages their ability to act on data with variable sequence lengths, such as
constituents inside a jet. In this document, we explore the application of
Recurrent Neural Networks (RNNs) and other sequence-based neural network
architectures to classify jets, regress jet-related quantities and to build a
physics-inspired jet representation, in connection to jet clustering
algorithms. In addition, alternatives to sequential data representations are
briefly discussed.

###### Contents

  1. 1 Introduction
    1. 1.1 Recurrent neural networks
    2. 1.2 Long-Term Dependencies, LSTMs and GRUs
  2. 2 Applications of RNNs to Jet Physics
    1. 2.1 Identifying Heavy Flavor Jets
    2. 2.2 Identifying Strange Jets
    3. 2.3 Identifying tau lepton jets
    4. 2.4 Identifying Top Jets
    5. 2.5 Learning a jet representation
  3. 3 Alternatives to RNNs
    1. 3.1 Recursive Neural Networks
    2. 3.2 Transformers and Deep Sets
  4. 4 Conclusion

##  1 Introduction

The area of Sequence-Based Learning in Machine Learning deals with the
concepts and algorithms used to learn from data represented as an ordered set
(sequence) of objects, each with its set of characteristics (features), in
which positional information of each object (context) is important. The idea
of contextual information as being important for the algorithms is
fundamental, since it can encode correlations between objects along the
sequence. One of the main applications of this class of models is in natural
language processing (NLP). In these cases, the sequence is often built from
words in a sentence and the algorithm must learn from it. How the learning
occurs and what is learned will depend on the application. To perform a
translation task (neural machine translation), for example, the algorithm must
output another sequence. In other instances, a summary semantic information
needs to be obtained, such as when the algorithm needs to classify a certain
sentence as positive or negative in an on-line product review, for example.

In more mathematical terms, sequence-based models aim to perform operations
$f$ on a sequence of inputs $\\{\mathbf{x}^{t}\\}$, where each entry
$\mathbf{x}^{t}$ is a vector of features, and $t$ is a position in the ordered
sequence with a length $T$, as shown in the scheme presented in Fig. 1. In
particle physics terms, the sequence can represent an ordered set of tracks
that constitutes a jet, for example, while the entry $\mathbf{x}^{t}$
represents the kinematics of the track in the position $t$ in that sequence.
An algorithm acting on the sequence $J=\\{\mathbf{x}^{t}\\}$ can then be used
to learn information about that jet, such as its flavor or charge (as will be
discussed in the following sections).

In general, different sequences being utilized for the algorithm definition
will have different lengths, just as jets can have different number of tracks.
If all sequences have fixed length, they can be collapsed into a single
feature vector and simpler algorithms, such as densely-connected NNs, can be
used. However, even for fixed-length sequence problems, sequence-based models
can outperform simpler models by exploiting the ordered nature of the input
data.

$\mathbf{x}^{t-1}$$\mathbf{x}^{t}$$\mathbf{x}^{t+1}$$f$$f$$f$$f$ Figure 1:
Scheme of a sequence-based algorithm acting on a sequence of inputs
$\mathbf{x}^{t}$.

A key feature of sequence-based models is the ability to share parameters in
different parts of the same model, i.e., the operation $f$, which is learned,
is applied to every step in the sequence. With parameter sharing, these models
are able to generalize to sequences of different lengths, using the same set
of parameters throughout the input elements. If different parameters were to
be learned individually, the desired generalizability would not be
achieavable, and the model would behave similar to a densely-connected
network.

###  1.1 Recurrent neural networks

A recurrent neural network (RNN) is a neural network implementation of the
concepts described in the previous section. Densely-connected neural networks
map an input vector of features $\mathbf{x}$ to an output vector $\mathbf{o}$.
In contrast, RNNs map a sequence of inputs ${\mathbf{x}}$ into an output,
which can be a vector or a sequence as well. This difference can be achieved
in many different ways, but RNN architectures generally present cyclical
connections between units in the same or different layers.

These interconnections between units are sequential, in the sense that each
unit’s hidden state is obtained by a combination of the previous unit’s hidden
state and the input from that step. This means that instead of a unit’s hidden
state be given by $h=f(\mathbf{x};\theta)$, it will be given by
$h^{t}=f(h^{t-1},\mathbf{x}^{t};\theta)$. On the particular case in which $f$
is given by a hyperbolic tangent function, for example, the RNN can be
represented as

$h^{t}=\tanh{\left(\mathbf{W}^{\top}\mathbf{x}^{t}+\mathbf{V}^{\top}h^{t-1}+b\right)},$
(1)

instead of a simple densely-connected unit
$h=\tanh{\left(\mathbf{W}^{\top}\mathbf{x}+b\right)}$, where $\mathbf{W}$,
$\mathbf{T}$, and $b$ represent learnable weights and biases.

Notice that in eq. 1 the weights on the operations ($\theta$, or $\mathbf{W}$,
$\mathbf{V}$, and $b$ explicitly) do not depend on the time step, or the
sequence element position, $t$, explicitly showcasing the parameter sharing
feature of the RNN. Similarly to MLPs functioning as universal approximators,
large enough RNNs have been shown to universaly approximate any measurable
sequence-to-sequence maps HAMMER2000107 .

In general, RNN architectures add other features on top of the recurrent layer
format described above. For example, in tasks where the algorithm needs to
output another sequence, each individual hidden state in the recurrent layer
might be read out into a densely connected network. In contrast, the cases
where only a single output is read at the end of the sequential layer (also
known as time-unfolded RNNs) are used to extract a summary information of the
input sequence. Even though the presence of cycles in these architectures
could potentially complicate the process of updating the network parameters
during the optimization step, backpropagation can still be applied to the
unrolled computation graph111The application of backpropagation on the time-
series unrolled RNN is known as backpropagation through time., thus no
specialized algorithms are necessary.

A simple but powerful extension of standard RNN architectures are
Bidirectional Recurrent Neural Networks (BRNNs) SchusterPaliwal1997 . For
certain sequence-based modeling applications, knowledge about backward-in-time
context might be as important as forward-in-time, only the latter of which is
exploited in standard RNNs. BRNNs manage to extend that context information by
adding a second recurrent layer to the network architecture which processes
the sequence in reversed order. Both the forward and backward RNNs are then
connected to the same output layer, providing a combined representation.

An idea related to RNN architectures which is also used in time-series and
signal processing analyses are 1D convolutional layers. In these
architectures, the convolution operation kernal acts on neighbouring time
steps, and outputs a new sequence based on its inputs. These operations act on
fixed length sequences, where empty entries can be masked - similarly to 2D
CNNs acting on sparse images with a fixed pixel grid. Parameter sharing is
also an important feature here, exploited by the use of a single convolution
operation across different time steps. One drawback of this method is its
limited sensitivity to long-term dependencies, only exploring correlations
across close neighbours, as defined by the kernel length.

###  1.2 Long-Term Dependencies, LSTMs and GRUs

When dealing with large sequences, one expected behavior of RNNs is to learn
how correlated certain entries are, regardless of how far appart they appear
in the input sequence. In practical tearms, this implies that information from
early entries in the sequence must be encoded in how the network learns about
latter entries (long-term dependencies). Unfortunately, in the learning steps,
it is common for gradients propagated through many steps to either quickly
approach zero or increase - both undesirable features in terms of
optimization. This issue is referred in the literature as the vanishing or
exploding gradient problem (e.g. in Pascanu2012 ). Vanishing gradients, for
example, can lead to long-term dependencies being given less importance
through the sequence compared to short-term ones.

This issue can be understood by imagining a simplified recurrent structure as
a linear transformation between hidden states
$h_{t}=\mathbf{V}^{\top}h_{t-1}$. This operation will be therefore performed
$T$ times when moving from the time-step 0 to the last sequence entry. This
shows that the learnable parameters in $\mathbf{V}$ will be raised to the
power of $T$, which means that weights less than 1 will tend to approximate 0
at later steps, while weights larger than 1 will quickly grow. Most modern RNN
architectures solve this problem with the usage of Long Short-Term Memory
HochreiterSchmidhuber1997  (LSTM) units or Gated Recurrent Units Cho2014
(GRUs).

LSTM units mitigate the vanishing gradient problem by introducing a new path
in the recurrent loop where the information coming from each sequence entry
can flow for long durations, possibly without interference from subsequent
hidden states. This path is dynamically gated through learnable parameters,
which means that the importance of long-term dependencies is optimized with
the rest of the network parameters. In particular, if the activation function
pertaining to a gate remains close to 0, the information from the previous
time-step will not be propagated throughout the sequence. However, the LSTM
will still use that time-step information to update its current hidden state.
This optimizable gated structure ensures that both long-term and short-term
contributions to the gradient are taken into account.

A schemematic view of three sequential LSTM units is shown in Figure 2. The
LSTM receives at the time step $t$ both the hidden state of the previous time
step ($h_{t-1}$), which is concatenated with the feature vector
$\mathbf{x}_{t}$, and an extra input (the cell state, $C_{t-1}$) which is
regulated by a forget gate ($f_{t}$). The forget gate can be a neural network
itself, with a sigmoid output which is shared across all units, ensuring that
the LSTM actively learns how much long-term correlations should be propagated
in the sequence. Next, the cell state is updated with information learned from
the previous hidden state and $\mathbf{x}_{t}$ with a dedicated neural
network, generating an intermediate cell state $\tilde{C}_{t}$. The impact of
$\tilde{C}_{t}$ in the final cell state is regulated by another neural
network, $i_{t}$. After the update, $C_{t}$ is propagated to the next time-
step. The final hidden state at this time step is derived with information
learned from the previous hidden state and $\mathbf{x}_{t}$ through a network
($o_{t}$), and also $\tilde{C}_{t}$.

Figure 2: Schematic structure of a LSTM recurrent layer. Two outputs for the
unit’s hidden state $h_{t}$ are shown, representing the case in which the LSTM
layer outputs another sequence. Image adapted from colah_lstm .

GRUs work with a similar gated structure - however, the same forget gate that
decides on the propagation of the previous cell state also determines, with an
inverse importance, how the current cell state should be updated with
$h_{t-1}$ and $\mathbf{x}_{t}$. With the LSTM nomenclature used above, the
forget gate and the cell state update gate would be given by $f_{t}=u_{t}$ and
$i_{t}=(1-u_{t})$, respectively, where $u_{t}$ is called the update gate. GRUs
comparatively uses less trainable parameters than LSTMs, which can be
benefitial for smaller datasets; its less flexible architecture can
potentially be detrimental in more complex applications - however, the studies
to be shown below that compared GRUs and LSTMs generally see similar
performances.

##  2 Applications of RNNs to Jet Physics

Representing reconstructed jets in collider experiments as images for
classification and other machine learning based tasks, has been a successful
avenue of research for years now.There are, however, a few features related to
this choice of representation that can make the process of training a computer
vision-based algorithm difficult when compared, for example, to the training
of a simple densely-connected neural network based on engineered features.

Jet images can be very sparse, i.e., containing few populated pixels222Initial
computer vision-based jet images studies have reported $5\sim 10\%$ of
activated pixels on average for 25 x 25 pixels images from signal and
background-type jets OliveiraJetImages ., making the identification of
features on individual jets a complicated task, even if identifying these
features on their averaged images might be very easy. Pre-processing steps can
be applied to reduce sparseness, for example, increasing the coarseness the
image further by combining adjacent pixels. This procedure, however,
effectively reduces the spatial information contained in the input image,
penalizing the algorithm’s performance.

Another issue arising in the pre-processing step is finding a unique
geometrical representation for these images that reflects the expected
symmetries of the problem. In general, geometrical transformations (e.g.,
rotations, translations and reflections) are used, aligning pre-defined axes
based on the jets’ spatial energy distributions - but these definitions can be
very task specific and not well generalizable. This is specially true when the
pattern of energy deposition inside these jets displays a different number of
core clusters (prongness) - for example, when comparing quark-initiated jets,
jets from a collimated hadronic W boson decay, and jets from a collimated top
quark decay.

The algorithms to be described below think about the jet instead as collection
of correlated objects, each one with its set of characteristics (such as
position and energy), and apply some of the sequence-based ML ideas discussed
above for different types of tasks. This idea is reminiscent of the actual way
jets are built in collider experiments, using sequential clustering algorithms
(for example, the widely used in LHC experiments anti-$k_{t}$ algorithm
CacciariAntiKt ), acting on low-level detector quantities, such as calorimeter
deposits or tracks.

###  2.1 Identifying Heavy Flavor Jets

The identification of jets originated from the products of the hadronization
process of heavy quarks (bottom and charm quarks), known as heavy flavor jets,
is of fundamental importance for experiments at the LHC for two main reasons.
Firstly, the Higgs boson - discovered in 2012 by ATLAS and CMS, and currently
the focus of intense experimental investigation - mainly decays to a pair of
bottom quarks, with a predicted decay branching fraction of about $56\%$
Heinemeyer:2013tqa . Secondly, the top quark, which together with the Higgs
boson can help us learn about the structure of the electroweak vacuum
Buttazzo_2013 , decays almost entirely to a bottom quark plus a W boson
Zyla:2020zbs .

Finding these Higgs and top decays is not an easy task due to the enormous
multijet background events at the LHC (events with at least two reconstructed
jets). These multijet events are produced with cross sections over 3 times
larger than events with top quarks, and 4 times larger than events with a
Higgs boson ATL-PHYS-PUB-2019-010 . Fortunately, most of these events contain
jets from light quarks (quarks other than charm, botton and top) and gluons
(collectively grouped as light jets). Therefore, learning how to separate
heavy flavor signal jets from the light jets background is necessary.

Hadrons containing bottom and charm quarks are heavy - for example, the
$B^{0}$ meson invariant mass is about 5.3 GeV/c2 Zyla:2020zbs . They decay
through the Weak force, thus having a long mean lifetime ($\tau_{B^{0}}\simeq
1.5\,\times\,10^{-12}$ s). Therefore, a $B^{0}$ with a momentum of 50 GeV/c
will travel on average over 4.5 mm before decaying. This displaced decay can
be reconstructed as a secondary vertex, i.e., a vertex that is separated from
the collision’s primary vertex (where the initial proton-proton interaction
ocurred, in the case of the LHC). The LHC experiment’s inner trackers
(reponsible for reconstructing the trajectory of charged particles) have been
built with the intent of separating these secondary vertices with good
precision, in order to identify heavy flavor jets. The decay chain of a
$B^{0}$ meson, as simulated by the ATLAS experiment, is shown in Figure 3 ATL-
PHYS-PUB-2018-025 , which also includes the tertiary vertex from the displaced
decay of the $D^{0}$ meson, containing a charm quark.

Figure 3: Simulated decay chain of a $B^{0}$ meson in the ATLAS experiment.
This event display was obtained from a simulated dataset of top quark pair
production, at center-of-mass energy of 13 TeV. It shows the displacement of
the vertices produced from $B^{0}$ (secondary vertex) and subsequent $D^{0}$
(tertiary vertex) decays, with respect to the center of the coordinate system
(primary vertex).

In general, requiring that a secondary vertex needs to be reconstructed in
order to identify a b-jet can be detrimental for a high efficiency algorithm.
However, tracks originated from this displaced location will have very
particular characteristics when compared to tracks from the primary vertex,
and will be correlated by their shared origin. Therefore, algorithms that
focus on finding correlations between tracks perform comparatively well with
respect to direct secondary vertex finding, and can provide complimentary
information for a combined heavy versus light jet discrimination.

In the ATLAS experiment, two different types of algorithms have been developed
to identify heavy flavor jets based on the likelihood of tracks being
originated from secondary vertices. Both use the tracks’ impact parameter
information, which encodes the distance of closest approach of the charged
particle’s trajectory with respect to the primary vertex. Particles originated
from the primary vertex will have small impact parameters, while particles
from secondary vertices will tend to have larger impact parameters. An
important related quantity is the impact parameter significance, in which the
distance is divided by the uncertainty in its measurement. Utilizing the
significance minimizes the impact of low-quality tracks with large mis-
measured impact parameters.

IP3D Aad_2019_IP3D , one of the first ATLAS algorithm based on tracks impact
parameter information, treats the tracks as independent entities, ignoring
possible correlations. It uses 3D histograms of transverse and longitudinal
impact parameter significances ($\mathcal{S}_{d_{0}}$ and
$\mathcal{S}_{z_{0}}$, respectively), and a track quality grade, built from
simulation and separately for b-jets, c-jets and light flavor jets. Per-flavor
conditional likelihoods are calculated from these histograms for each track in
a jet. With a naïve Bayes approach, a final jet-level likelihood is built by
multiplying the individual tracks likelihoods. However, important information
is lost with the assumption that the tracks in the jet are uncorrelated. This
can be seen in Figure 4, where a strong correlation between the
$\mathcal{S}_{d_{0}}$ distribution of the leading and subleading track in the
jet (ordered by $\mathcal{S}_{d_{0}}$) can be seen near the diagonal for
b-jets but not for light jets.

Figure 4: 2D histogram of $\mathcal{S}_{d_{0}}$ for leading (horizontal axis)
and subleading (vertical axis) tracks inside a b-jet (left) and a light flavor
jet (right). The correlation observed is an indication that the naïve Bayes
approach of the IP3D algorithm is not enough to exploit the full information
contained in the tracks’ impact parameters with respect to b-jet
identification.

More recently, algorithms exploiting RNN architectures based on LSTMs have
been proposed to perform the task described above, treating the list of tracks
within the jet as the input sequence to the algorithm. Due to the variable
number of tracks within a jet, RNNs are better suited than dense architectures
in this approach. Even though no natural track ordering is clear to the
problem, tracks with larger impact parameters are more likely to come from
heavy flavor jets, therefore, impact parameter based ordering is a good
ansatz333Empirically, the studies mentioned below have shown that certain
track orderings work better than others.. While the most basic version of this
algorithm acts on tracks only, proposals have been made to combine track and
secondary vertices information in a single LSTM-based architecture
GuestFlavorDNN .

The ATLAS implementation of the LSTM-based architecture for heavy jets
identification with tracks’ impact parameters is called RNNIP ATLAS-RNNIP . It
treats the tracks within a jet as a sequence, and uses the impact parameter
significance information and track kinematics as features. It also uses
categorical information based on the track reconstruction quality in an
embedded layer. These categories separate high quality, well measured tracks,
in which a better impact parameter resolution is expected, based on detector-
level information such as the number of hits in the innermost tracker layer.

The tracks are ordered by $\mathcal{S}_{d_{0}}$, although other orderings
(such as by track $p_{T}$) have shown similar performance. The algorithm is
trained in a simulated sample of top pairs, which provide a dataset enriched
of both heavy and light quarks, and outputs a probability of a given jet to be
a bottom jet ($b$-jet), a charm jet ($c$-jet), or a light flavor jet. The
three probabilities are then combined into a likelihood that is used for
discrimination.

A comparison between the performance of the naïve Bayes algorithm described
above (IP3D) and RNNIP is shown in Figure 5. The efficiency of identifying
$b$-jets is plotted on the horizontal axis, while one over the probability of
identifying light flavor jets (misidentification probability) as $b$-jets is
plotted on the vertical axis. The RNNIP algorithm displays a better light
flavor jet discrimination for every value of $b$-jet efficiency, and is
comparable to a boosted decision tree that combines IP3D with the secondary
vertex-based ATLAS algorithms (MV2c10 Aad_2019_IP3D ), even though it does not
explicitly reconstruct secondary vertices. Performances were measured for jets
clustered with the anti-$k_{T}$ algorithm CacciariAntiKt  with R = 0.4, with a
transverse momentum above 20 GeV, in a simulated dataset of top quark pairs,
at center-of-mass energy of 13 TeV.

Figure 5: Performance of heavy flavor identification algorithms in the ATLAS
experiment ATLAS-RNNIP . The horizontal axis shows the efficiency of correctly
identifying $b$-jets, while the vertical axis shows one over the efficiency of
incorrectly identifying light flavor jets as $b$-jets. The red dashed curve
shows the performance of the RNN-based algorithm (RNNIP), while the dashed
blue curve shows the performance of a naïve Bayes algorithm acting on similar
inputs (IP3D).

Figure 6 shows the Pearson’s correlation coefficient $\rho$ between the RNNIP
likelihood, and $\mathcal{S}_{d_{0}}$ and $\mathcal{S}_{z_{0}}$ for each track
in the sequence. It is interesting to note that stronger correlations in
b-jets are seen for impact parameter significances of the first $\sim 8$
tracks, which may be related to the expected charged particle multiplicity of
b-hadron decays. This shows that the network architecture is able to learn
contextual information from the given sequence ordering.

Figure 6: Pearson’s correlation factor between RNNIP likelihood, and
transverse and longitudinal impact parameter significances. Correlations are
shown separately for b-jets (orange), c-jets (red), light flavor jets (blue).

Heavy flavor jets identification in the CMS experiment shares many
similarities with the strategies employed by ATLAS. In particular, their final
discriminant is also a combination of information pertaining to secondary
vertexing and the set of tracks inside the jet. Two sets of algorithms have
been developed with this intent: one set combining engineered features
extracted from the jet, and one that directly uses reconstructed objects
information into a neural network architecture which includes LSTM layers. The
CMS DeepCSV algorithm Sirunyan_2018  exemplifies the first strategy, similarly
to the ATLAS MV2c10 boosted decision tree. It is based on a densely-connected
neural network with 4 hidden layers, with inputs that are defined by other
algorithms which act directly on tracks’ impact parameter information and
reconstructed secondary vertices.

Significant improvement has been observed by CMS by moving to the DeepFlavour
algorithm CMS-DP-2018-033 , which acts directly on these low-level
observables. The DeepFlavour network receives three sequences as inputs: a
sequence of charged particles (reconstructed from tracks and calorimeter
clusters), a sequence of neutral particles (calorimeter clusters with no
associated tracks), and a sequence of reconstructed secondary vertices. Each
sequence is processed by a 1D convolutional layer, which learns a shared
representation that is specific for each type of sequence. The convolutional
layer outputs are then fed to three different LSTM layers, which summarized
the sequences information into three fixed length feature vector. These
features are combined with additional jet-level information in a densely
connected network, which outputs the jet flavor probabilities.

Figure 7 summarizes the CMS $b$-jet identification performance. The red and
blue lines represent the performance of the DeepCSV and DeepFlavour algorithms
respectively, with the CMS detector conditions present during the 2017 LHC
data taking period (Phase 1), while the green line shows the DeepCSV
performance with the CMS detector conditions in 2016 (Phase 0). Between 2016
and 2017, CMS inner tracking detector was upgraded to deal with the harsher
radiation conditions at the LHC later Run 2 years. This upgrade also provided
the CMS experiment a better impact parameter resolution, directly improving
their heavy flavor identification performance. The $b$-jet efficiency is shown
as a function of the light flavor jets misidentification probability (full
lines), and as a function of the $c$-jets misidentification probability
(dashed lines). While a large improvement is seen in the performance of the
DeepCSV algorithm with the CMS Phase 1 inner tracker upgrade, an improvement
just as large is seen with the usage of DeepFlavour in terms of light flavor
jet discrimination, with an even large gain in terms of $c$-jet rejection.
Performances were measured for jets clustered with the anti-$k_{T}$ algorithm,
with a transverse momentum above 30 GeV, in a simulated dataset of top quark
pairs, at center-of-mass energy of 13 TeV.

Figure 7: Performance of heavy flavor identification algorithms in the CMS
experiment CMS-DP-2018-033 . The horizontal axis shows the efficiency of
correctly identifying $b$-jets, while the vertical axis shows the efficiency
of identifying light flavor jets (full lines) or $c$-jets (dashed lines) as
$b$-jets. The red and blue curves represent the DeepCSV and DeepFlavor
algorithms performances with the 2017 CMS detector conditions, while the green
curve represents the DeepCSV algorithm with the 2016 CMS detector conditions.

###  2.2 Identifying Strange Jets

Jets from strange quarks are grouped within the light flavor jets category for
the algorithms described above. However, for certain physics applications,
such as the direct measurement of the $|V_{ts}|$ element of the CKM matrix
through the search of the rare decay $t\rightarrow W^{+}s$ (or
$\bar{t}\rightarrow W^{-}\bar{s}$) Ali_2010_vts , discriminating strange jets
from first-generation jets is a necessity.

A LSTM-based algorithm has recently been proposed to tackle this problem
ErdmannStrangeTagging . The study is performed with a simplified detector
description model (Delphes delphes ) based on a CMS-like detector. The
algorithm is trained to discriminate between strange jets and jets from the
hadronization of up and down quarks, produced by proton-proton QCD
interactions. Similarly to the ATLAS RNNIP, the proposed algorithm acts on
sequences of tracks, with features based on their impact parameters and
kinematics with respect to the jet.

Secondary vertices are expected to be present in strange jets through the
decay of strange kaons into $\pi\pi$ and lambda baryons into $p\pi$.
Therefore, all possible secondary vertices in the jet are reconstructed by
pairing tracks with small distances of closest approach to each other. These
vertices are then used to define a track ordering based on the parameter $R$
assigned to each track. This parameter is defined either by the transverse
distance between the primary vertex and the secondary vertex to which that
track belongs, or, in the cases where the track is not associated to a
secondary vertex, the innermost tracker hit belonging to that track. If
multiple tracks receive the same $R$ (e.g., two tracks from the same secondary
vertex), their ordering is performed by $p_{T}$. This ordering ensures that
adjacent tracks belong to the same secondary vertex, which the network will
use to learn about displaced decays.

Figure 8 compares the performance of the LSTM-based strange jet discriminator
to the performance of simpler methods. In particular, the performance of using
the transverse momentum fraction $x_{K}$ and $x_{\Lambda}$ of identified kaons
and lambda baryons is also investigated. To calculate these quantities, a
selection on the invariant mass of the reconstructed secondary vertices,
consistent with the $K$ and $\Lambda$ masses, is applied. The highest $p_{T}$
$K$ and $\Lambda$ candidates of the remaining vertices are chosen, and used to
calculate $x_{K}$ and $x_{\Lambda}$ to their parent jet. Two different setups
of the LSTM are compared: one including all selected jets, and one only
including jets that contain at least one track with large transverse impact
parameter $|d_{0}|>1$ mm.

Figure 8: Performance of strange jet identification algorithms based on LSTM
architectures and strange hadron reconstruction ErdmannStrangeTagging . The
vertical axis shows the efficiency of correctly identifying strange jets,
while the horizontal axis shows the efficiency of incorrectly identifying up
and down quark jets as strange jets. The blue and orange lines show the
performance of the LSTM-based algorithm using all jets and jets with at least
one track with transverse impact parameter $|d_{0}|>1$ mm, respectively. The
red dashed and full black lines show the performance of selecting on the
transverse momentum fractions $x_{K}$ and $x_{\Lambda}$, and only on $x_{K}$,
respectively.

The overall performance of this algorithm is unfortunately limited by the
similarity between strange and first-generation quark jets - achieving a
background efficiency of $21\%$ ($63\%$) for a signal efficiency of $30\%$
($70\%$). The secondary vertices from kaon decays are not as displaced as
vertices from $b$-jets or $c$-jets and can easily be misidentified from
vertices produced by material interaction or decays of particles originated in
the hadronization process. One important improvement with respect to using
$x_{K}$ and $x_{\Lambda}$, however, is the ability to achieve higher signal
efficiencies, as shown in Figure 8.

###  2.3 Identifying tau lepton jets

Similarly to heavy flavor jet identification, identification of tau leptons is
particularly interesting due to its ties to Higgs physics. The coupling
between the Higgs and the tau lepton is the largest Higgs couplings to leptons
in the Standard Model. It is therefore an opportunity to directly measure the
structure of the Higgs Yukawa couplings to that sector of the Standard Model.

Taus decay either leptonically ($\tau\rightarrow\nu_{\tau}+\ell\nu_{\ell}$, in
which $\ell$ is an electron or a muon), or hadronically
($\tau\rightarrow\nu_{\tau}+$ hadrons). Leptonic tau decays are roughly
indistinguishable from isolated leptons in hadron collider experiments.
Therefore, tau identification focuses on hadronic taus, which represent a
branching fraction of approximately $65\%$ Zyla:2020zbs . Hadronic tau decays
usually include one or three charged pions and one or more neutral pions.
Therefore, these decays are seen in the detector as narrow jets with one or
more tracks.

Since neural pions do not leave signals on the detectors’ inner tracks, their
trajectories cannot be reconstructed as tracks. This means that if only track
information were used, a large portion of information for tau identification
would be missing. Therefore, an optimal strategy should aim to combine the
tracking and calorimetry information.

The ATLAS experiment state-of-the-art tau identification algorithm is based on
a double LSTM architecture that combines track sequences and calorimeter
deposits (clusters) sequences ATLAS-RNNTau . The algorithm has three sets of
inputs: a track sequence, a cluster sequence, and a set of high-level
variables connected to a dense layer. The track and cluster features
considered refer to their kinematics and detector-level properties, while the
the high-level ones are related to the jet itself, or engineered features
based on the collection of jet constituents.

Tracks are individually fed through dense layers with shared weights, so that
an embedded representation can be learned. The same procedure is applied to
calorimeter clusters separately. The two processed sequences, one of track
embeddings and one of cluster embeddings, are ordered by decreasing $p_{T}$ of
the original objects and used as inputs to two separate LSTM blocks. These
blocks includes two layers of LSTM units - the first one maps the sequence in
the learned representation into another sequence of the same length. The
second LSTM layer only outputs the information at the last time-step, thus
providing a summary of the input sequence. The LSTM blocks outputs are fed to
a densely-connected block, which also receives information from a densely-
connected block encoding high-level observables of the tau jet.

The training and evaluation of the architecture defined above is performed
separately for the cases in which the tau decay includes one or three tracks.
The tau decays (signal) are provided by simulation of
$\gamma^{*}\rightarrow\tau\tau$ events, while background jets are selected
from a simulation of dijet events.

Figure 9: Performance of hadronic tau jets algorithms in the ATLAS experiment
ATLAS-RNNTau . The horizontal axis shows the efficiency of correctly
identifying hadronic tau jets, while the vertical axis shows the inverse of
the efficiency of incorrectly identifying quark-initiated jets as hadronic tau
jets. The red curves show the performance exclusively on tau jets with a
single track (1-prong), while the blue curves represent tau jets with three
tracks (3-prong). The RNN-based model’s performance is shown in the full
lines, with the dashed lines representing an algorithm with similar inputs but
based on a boosted decision tree.

The performance obtained in Figure 9 compares the LSTM architecture (RNN)
optimized for tau decays with one (1-prong) and three (3-prong) tracks. It
also compares to the previous algorithm used in the ATLAS experiment, based on
a boosted decision tree (dashed lines). The LSTM-based architecture
outperforms the previous baseline for all hadronic tau efficiencies. These
improvements have been shown to be significant enough that the new
architecture has actually been used for identifying tau candidates at the
ATLAS High-Level Trigger (HLT) in 2018.

###  2.4 Identifying Top Jets

When produced at large momenta, top quarks’ decay products will start to
merge, making it more difficult to resolve them spatially in the detector. In
this regime, the entire top quark decay can be clustered into a single jet.
Usually these jets have larger radius parameters: CMS cms_largeR  utilizes R =
0.8 anti-$k_{t}$ jets and R = 1.5 Cambridge-Aachen jets cajets , while ATLAS
atlas_top  focuses on R = 1.0 anti-$k_{t}$ jets. Identifying these boosted top
objects, from a background of jets from the hadronization of lighter quarks
and gluons, is particularly interesting when searching for Beyond the Standard
Model physics which predict TeV-range resonances decaying to top pairs.

Several interesting features which are present in a boosted top jet can be
used to discriminate against a high momentum jet produced by QCD interactions.
In particular, hadronic top decays will tend to be three-pronged, with each
prong corresponding to a final state particle in the $t\rightarrow
bW\rightarrow bqq^{\prime}$ decay chain. Two important details on this chain
is that one of these prongs will be consistent with a heavy flavor jet, and
the other two will be consistent with a W boson decay.

LSTM-based architectures have been proposed for identifying these boosted top
jets EganTopTagging . The study is performed with a Delphes-based detector
simulation delphes  with a particle flow type of particle reconstruction,
overlaying minimum bias events to emulate the LHC 2016 collisions conditions,
averaging of 23 proton-proton interaction per event. Jets are clustered
following the ATLAS strategy, with the anti-$k_{T}$ algorithm and R = 1.0. The
signal top jets are obtained from simulating a beyond the Standard Model
process in which a $Z^{\prime}$ boson with masses ranging from 1400-6360 GeV
decays to hadronically decaying $t\bar{t}$ pairs. Background jets come from
the simulation pure QCD hard scattering processes (QCD jets).

Similar to tau identification, the sequence for the recurrent model is built
of calorimeter clusters. The input ordering is defined based on going through
the clustering history of the jet, starting from the final reconstructed jet,
and adding constituents to the sequence as they appear in each step. The
decision on which path to follow in each time two parent nodes merge depends
on the anti-$k_{T}$ distance metric $d_{ij}$. If the two parent nodes are
present in the list of jet constituents, they are added to the sequence
ordered by $p_{T}$. A scheme presenting how the clustering tree is used to
build the sequence ordering is shown in Figure 10. Another strategy is based
on reconstructing R = 0.2 anti-$k_{T}$ subjets, order them by descending
$p_{T}$, then adding constituents to the sequence depending to which subjet
they belong - constituents on the same subjet are also ordered by descending
$p_{T}$. These schemes are compared with purely odering on the jet
constituents $p_{T}$.

Figure 10:  Scheme of jet clustering history used to define input sequence
ordering in LSTM-based top identification algorithm EganTopTagging . The
anti-$k_{T}$ distance $d_{ij}$ is used to choose which path to follow in the
tree.

As seen in Figure 11 (left), the LSTM-based architecture with substructure
oredering outperforms the previous baseline studied by the same group top_dnn
, based on a fully connected dense neural network. The results are presented
depending on the ordering scheme as well whether the jets are trimmed or not.
Trimming trimming  is a technique that ensures robustness of the jet
kinematics, especially the jet mass, with respect to pile-up contamination by
removing R = 0.2 subjet constituents with a certain energy fraction; in this
case, subjets are removed if their energy fraction is lower than $5\%$. Figure
11 (right) shows that the impact of the ordering strategy is limited. In fact,
it appears to be smaller than the impact of applying trimming, which removes
subjet information that could be important for top identification. For more
information on top jet identification with jet images, see section XYZ.

Figure 11:  Performance of top identification with LSTM-based techniques
EganTopTagging . The horizontal axis shows the efficiency of correctly
identifying top jets, while the vertical axis shows one over the efficiency of
incorrectly identifying QCD jets as top jets. The left plot compares the LSTM
algorithm with a strategy using a dense neural network. The right plot
compares different training strategies for the LSTM, including different
sequence oderings (substructure, subjet or purely constituent $p_{T}$ based)
and with or without trimming.

###  2.5 Learning a jet representation

In general, a jet can be understood as a collection of final state particles
that are results of the hadronization and fragmentation processes of a quark
or a gluon. However, the actual procedure of reverse-engineering these
processes, which are quantum mechanical by nature, can be complicated,
especially when dealing with the busy environment of a hadron collider event.
Jet clustering algorithms aim to reduce this complexity and make the
connection between observable final state particles and the phenomenological
predictions based on partons.

While many types of jet clustering algorithms have been proposed in the past,
most recent applications utilize strategies based on sequential clustering,
such as the ones already discussed here ($k_{T}$, anti-$k_{T}$ and Cambridge-
Aachen being three well-known examples). These algorithms take advante from
the fact that the processes for the particle shower development can be
approximated well, due to factorization theorems, by a sequence of
$2\rightarrow 1$ splittings in creating a hierarchical representation of the
jet. When reverse-engineered, this sequence of splittings becomes a sequence
of mergings of jet constituents, forming a binary tree. This binary tree
represents the jet clustering history, and encodes important information about
the nature of that jet.

As presented in the previous section, the usage of RNN architectures for
learning jet labels has been a successful avenue in collider experiments,
greatly improving on previous, often already Machine Learning-based,
strategies. Treating the jet constituents simply as an input sequence,
however, des not fully exploit the full information contained in the jet
clustering history as represented by the clustering tree.

It has been suggested that through understanding the interplay between the
clustering history and neural network architectures could lead to a more
natural jet representation AndreassenJUNIPR . The study presents a framework
in which a neural network can be built based on the clustering tree obtained
from the jet clustering algorithm history. It tries to exploit the entire
physical knowledge contained within the jet, particularly with respect to
hierarchical correlations between individual constituents. This structure is
then used to learn an overall probabilistic model of the jet conditioned on
its constituents, through an unsupervised training task.

Building a probabilistic model for a class of jets presents the possibility of
having a tractable and differentiable distribution function which can inform
different aspects of jet physics. In particular, models trained for different
classes of jets, such as b-jets and light flavor jets, can be used to compute
likelihood ratios for discrimination. Models can also be sampled in order to
generate new jets, a process which can take a significant amount of time when
large datasets are required.

The framework, named JUNIPR ("Jets from UNsupervised Interpretable
PRobabilistic models"), is based on the factorization of the jet probabilistic
model into a product of probabilities given by each step of the clustering
tree. For a set of jet constituents denoted by their 4-momenta
$p_{1},\ldots,p_{n}$, it computes the probability density
$P_{\text{jet}}(\\{p_{1},\ldots,p_{n}\\})$ of this set to have arisen from the
specified model. This probability factorizes based on the $2\rightarrow 1$
clustering tree, so that
$P_{\text{jet}}(\\{p_{1},\ldots,p_{n}\\})=\prod_{t}^{t=n}P_{t}$, where $P_{t}$
represents the probability model of the branching step $t$. This factorization
is schematically shown in Figure 12.

Figure 12: JUNIPR calculation of jet-level probability model based on its
factorization to nodes in binary tree of the jet clustering history
AndreassenJUNIPR .

The study further models $P_{t}$ as the product of three independent
probabilities: $P_{\text{end}}$, probability over binary values predicting if
the tree stops after this split; $P_{\text{mother}}$, probability over
integers of how likely is that the mother in step $t$ indeed participates on
the splitting $1\rightarrow 2$; and $P_{\text{branch}}$, probability over
kinematic configurations of the three states involved in the branching step
$t$. Each one of these three probabilities are modeled with densely-connected
neural networks. They use as inputs the hidden state $h_{t}$ at branching step
$t$, calculated based on the recurrent relation below:

$h_{t}=\tanh{\left(V\cdot\left(k_{d_{1}}^{t},k_{d_{2}}^{t}\right)+W\cdot
h_{t-1}+b\right)},$ (2)

where $k_{d_{1}}^{t},k_{d_{2}}$ represent the 4-momenta of the two daughters
involved in the splitting step $t$, and $V$, $W$, and $b$ represent learnable
weights and biases. This equation has the same form as the one presented in
section 1.1 detailing the functioning of RNNs. The study also tried to replace
this simple RNN strategy in the equation above with more complex solutions,
such as LSTMs and GRUs. No significant improvement in performance was
observed; it was argued that this was due to the simplicity of the task,
acting on sequences with only two elements.

The algorithm is trained based on the full jet probabilistic model, maximizing
the log-likelihood over the examples in the training dataset, with stochastic
gradient descent. Datasets are generated by $e^{+}e^{-}$ collisions, without
detector simulation, at center-of-mass energy of 1 TeV. The two jets in the
final state are separated into hemispheres by the exclusive $k_{T}$ algorithm
exkt , and their clustering trees are obtained with the Cambridge-Aachen
algorithm.

Results of using the learned jet probability applied to jet classification,
through likelihood ratio computation, are shown in Figure 13. The blue line
represents the discrimination power between jets from quarks and jets from a
boosted hadronic Z decay by calculating the likelihood ratio between the jets
probabilities based on each hypotheses $P_{Z}(\text{jet})/P_{q}(\text{jet})$.
$P_{Z}(\text{jet})$ and $P_{q}(\text{jet})$ represent the JUNIPR probabilities
trained individually on Z jets and quark jets, respectively, and evaluated on
a given jet. The JUNIPR likelihood ratio greatly outperforms the strategy
based on engineered features representing the jet substructure.

Figure 13: Performance of JUNIPR used for binary classification, based on the
calculation of likelihood ratios AndreassenJUNIPR . The horizontal axis shows
the efficiency of identifying Z jets, while the vertical axis shows one minus
the efficiency of incorrectly identifying quark jets as Z jets. The image also
shows the discrimination performance by simply using a substructure variable
$\tau_{21}$ and constituents multiplicity information.

The framework has also been used for direct binary classification tasks
AndreassenBinaryJUNIPR . In this case, two JUNIPR networks are built based on
two different types of jets (quark jets and gluon jets, for example). The two
networks are trained with a cross entropy objective function, where the
individual jet probabilities from each jet type are defined by the JUNIPR
networks. The quark versus gluon discrimination achieved by this method was
seen to outperform standard approaches, such as CNNs on jet images. It also
significantly outperforms the strategy above of individually training JUNIPR
models and calculating their likelihood ratio.

##  3 Alternatives to RNNs

The applications of RNN architectures in collider experiments and
phenomenological studies has been an active area of recent developments,
successfully avoiding issues previously seen with CNN architectures, while
very often improving on their performances. Some features of these
architectures are, however, less desirable for certain problems. Certain the
physics problems utilizing RNNs do not have a well-defined, natural ordering
of the sequence elements. Choosing a specific order becomes a non-trivial step
of the data formatting, one that could lead to undesirable performance losses.
On the other hand, other physics problems might benefit from topological
structures that encode more complex relations than a sequence.

###  3.1 Recursive Neural Networks

Recursive Neural Networks (RecNNs) have been proposed in the literature as
possible generalizations of RNNs, in which the sequential computational graph
is replaced by a tree structure recursive_bottou . This means that a node
hidden state still depends on the previous step in the computational graph,
similarly to the RNN, but the step itself is defined by a binary tree instead
of a sequence. This feature opens up the possibility of adding extra domain
knowledge information into the network architecture itself when building the
tree. A scheme of a recursive structure based on a binary tree is shown in
Figure 14.

$\mathbf{x}^{0}$$\mathbf{x}^{1}$$\mathbf{x}^{2}$$\mathbf{x}^{3}$$h_{0}$$h_{1}$$h_{2}$$\mathbf{x}^{4}$$o$
Figure 14: Scheme of a recursive binary tree structure algorithm acting on a
sequence of inputs
$\mathbf{x}^{0},\mathbf{x}^{1},\mathbf{x}^{2},\mathbf{x}^{3},\mathbf{x}^{4}$.
The $h_{i}$ nodes correspond to hidden states performing the combination of
two other nodes, while the $o$ node represents the output.

When building a RecNN, entries in a sequence are combined via a learned
function, with a predetermined binary tree structure. Therefore, a hidden
state combining entries $\mathbf{x}^{i}$ and $\mathbf{x}^{j}$ is represented
by $h=f(\mathbf{x}^{i},\mathbf{x}^{j},\mathbf{\theta})$, where $\theta$
represents the learnable parameters. This function $f(.,.,\theta)$ is then
used in all binary combinations, which allows the architecture to act on
variable-length sequences444This is another instance of weights sharing in
neural networks, a necessity when dealing with variable length inputs..
Another advantage of RecNNs over RNNs is its lower complexity, since the
computations are not performed per sequence entry anymore.

Studies utilizing RecNNs with trees reproducing the jet clustering history as
a basis for jet representation have been performed in the context of jet
classification LouppeQCDAwareRecursiveNN ; ChengRecNNforQG ; FraserJetCharge
and will be detailed below. They show that RecNNs are able to learn a fixed-
length jet embedding from a variable length tree structure built from the jet
constituents. This embedding can be further used for different tasks, such as
classification and regression.

#### Applications to Jet Physics

Initial studies used the RecNN architecture for jet discrimination
LouppeQCDAwareRecursiveNN , using simulated events processed through a
simplified detector simulation (Delphes). The signal jets are comprised of
hadronically decaying W bosons reconstructed as a single R = 1.0 anti-$k_{T}$
jet (boosted jet), while the background is taken from purely QCD hard
scattering proton-proton collisions. The study performs a comparison between
RecNN architectures based on different binary trees, corresponding to
different jet clustering algorithms histories. The result of this comparison
is shown in Figure 15. It’s interesting to note that, even though the jets
have been initially clustered with the anti-$k_{T}$ algorithm, other
clustering histories, such as $k_{T}$, perform better. This is consistent with
previous observations that $k_{T}$ outperforms anti-$k_{T}$ in terms of
identifying substructures in jets. In general, this variation in performance
is an evidence for the strong dependence of RecNNs architectures on the choice
of binary tree topology.

The performance is also studied when adding gating to the RecNN nodes.
Similarly to LSTMs and GRUs, gated structures are used to regulate how much
information is passed through the binary tree. An improvement in performance
is observed with gating, but the importance of topology is still dominant in
the results.

This study also shows that the learned RecNN jet embedding can be used for
event-wide discrimination, identifying beyond the Standard Model processes
involving boosted W and Z jets from purely QCD processes. In this case, these
RecNN embeddings of the individual jets in an event are passed through as a
$p_{T}$-ordered sequence to a GRU-based architecture which performs the
classification task. It was observed that the usage of RecNN-based embeddings
improves significantly the event classification performance with respect to
using a GRU-only architecture acting on sequences of jets.

Figure 15: Performance of RecNN-based algorithm for identifying boosted W jets
with respect to jets produced via purely QCD interactions
LouppeQCDAwareRecursiveNN . The horizontal axis shows the efficiency of
correctly identifying W jets, while the vertical axis shows one over the
efficiency of incorrectly identifying QCD jets as W jets. The different line
colors correspond to different jet clustering algorithms used to build the
binary tree topology in the RecNN.

A similar strategy has also been employed for quark versus gluon jet
discrimination ChengRecNNforQG , showing a slight improvement on a baseline
boosted decision tree-based algorithm. Quark versus gluon discrimination is an
important avenue of work at the LHC, due to the enormous gluon background from
soft hadronic interactions, which have a strong impact on analyses with light
quarks in the final state. One important application of this class of
algorithms is identifying the hard scatter final state light jets involved in
production of the Higgs boson via vector boson fusion.

As in the previous study, simulated events are processed with Delphes and jets
are clustered with the anti-$k_{T}$ algorithm with R = 0.7 for high jets with
$p_{T}>1$ TeV and R = 0.4 for jets with lower $p_{T}$. Purely QCD events with
two jets from the hard scatter (dijet events) are produced separately for when
the hard scatter partons are gluons (background), or up, down or strange
quarks (signal). The discrimination achieved with the RecNN under different
jet reconstruction strategies is shown in Figure 16, together with a baseline
approach based on a BDT with engineered features (jet shape and kinematics).
Three types of jet reconstructions are compared: using calorimeter towers only
("nopflow"); using particle identification for neutral hadrons, photons, and
positively and negatively charged particles, encoded in one-hot vectors for
each jet constituent in the tree ("one-hot"); using a $p_{T}$ weighted jet
charge defined by the clustering tree ("ptwcharge"). Although little
difference is observed with respect to the investigated jet reconstruction
schemes, a significant improvement is obtained over the BDT baseline.

Figure 16: Performance of RecNN-based algorithm for discriminating quark and
gluon jets ChengRecNNforQG . The horizontal axis shows the efficiency of
correctly identifying quark initiated jets, while the vertical axis shows one
minus the efficiency of incorrectly identifying gluon jets as quark jets. The
red lines correspond to the RecNN models trained with different jet
reconstruction techniques. The blue line shows the performance of a baseline
BDT based on engineered features related to jet shape and kinematics.

Recent studies have also compared RecNNs to RNNs and CNNs when tasked to
estimate jet charges FraserJetCharge . In the same spirit of identifying jet
flavor, identifying jet charges can help with mitigating the enormous multijet
background present in hadronic LHC searches. Requiring that two jets forming a
neutral resonance have opposite charges could potentially reach that goal,
assuming a good charge reconstruction resolution is achieved.

For these studies, jets from up quarks were used as proxies for positively
charged jets, while jets from down quarks were used for negatively charged
jets. Jets were simulated from QCD hard scattering processes in proton-proton
interactions, and clustered from final state particles with the anti-$k_{T}$
algorithm with R = 0.4. To test the performance of CNN’s on jet images, the
jets were formated into $\delta\phi\times\delta\eta=33\times 33$ pixel box,
with pixel intensities refering to the transverse momenta going into that
specific pixel, and to the sum of track charges weighted by their momenta over
the tracks corresponding to that pixel.

Results are shown comparing different ML-based models acting on jet
constituents: standard CNNs, residual CNNs, RNNs and RecNNs. The RNN
architecture implemented is based on GRUs, with LSTMs performing similarly;
the jet constituents are ordered in the input sequence by their $p_{T}$, with
the ordering based on distance to the jet axis performing equally. These
different algorithms’ performances are presented in Figure 17 in terms of
Significance Improvement Curves (SIC), defined by
$\epsilon_{s}/\sqrt{\epsilon_{b}}$, where $\epsilon_{s}$ is the efficiency of
correctly identifying down quark initiated jets, and $\epsilon_{b}$ is the
efficiency of incorrectly identifying up quark jets as down jets. They are
compared to more standard strategies based on engineered features. Overall,
the algorithms acting on jet constituents consistently outperform the baseline
approaches.

Figure 17: Performance of ML-based algorithms acting on jet constituents to
discriminate up and down jets based on their electric charges FraserJetCharge
. The horizontal axis shows the efficiency of correctly identifying down quark
initiated jets ($\epsilon_{s}$), while the vertical axis shows the ratio
$\epsilon_{s}/\sqrt{\epsilon_{b}}$, where $\epsilon_{b}$ is the efficiency of
incorrectly identifying up quark jets as down jets.

###  3.2 Transformers and Deep Sets

As shown by some of the studies described above, the ordering choice can be
detrimental to the RNNs performance. One recent strategy proposed to overcome
this ordering dependency is attention mechanisms bahdanau2016neural . The main
working point of attention mechanisms can be understood with a simple NLP
example. Within translation problems, two languages often display different
semantic structures for the same sentence: “a yellow cat" in English becomes
“um gato amarelo" in Portuguese, with the words for “cat" and “yellow"
switching positions. This relationship is difficult to be learned through
standard RNNs in which the input and output sequences have a defined order. To
solve this issue, network architectures with attention will directly learn
correlations between the entries in the input sequence and the entries in the
output sequence, which might not be encoded in the sequences ordering.

Different strategies and architectures involving ideas related to neural
network attention have been proposed to mitigate the ordering issue. In
particular, Transformer networks vaswani2017attention , which employ a self-
attention technique, learning two-by-two correlations between the sequence
inputs themselves through attention weights. Transformers have become common
tools in NLP, with pretrained models such as BERT (Bidirectional Encoder
Representations from Transformers) devlin2019bert  being adopted by Google in
its search engine for better understanding search queries.

Another method which aims to exploit correlations on a variable-length input
structure is the Deep Sets architecture zaheer2018deep . Deep Sets are
particularly suited for situations in which the ordering is not well defined,
as it treats the variable-length input sequence as a permutation invariant
set. A similar Deep Sets architecture was used Komiske_2019_PFN  to learn
representations for events in collider experiments (Particle Flow Networks).
These Particle Flow Networks act on lists of particles, such as jet
constituents, learns a combined representation through dense layers with
shared weights, and sums them into a single, object-wide representation.

Figure 18: Performance of the RNNIP ATLAS-RNNIP  and DIPS ATL-PHYS-
PUB-2020-014  heavy flavor identification algorithms in the ATLAS experiment.
The horizontal axis shows the efficiency of correctly identifying $b$-jets,
while the vertical axis shows the inverse of the efficiency of incorrectly
identifying light flavor jets as $b$-jets. The violet band shows the
performance the DIPS performance, while the green band represents RNNIP. The
width and central value of the curves represent the standard deviation and
mean of the light flavor jets rejection for a given $b$-jet efficiency for
five different network trainings. Performances were measured for jets with a
transverse momentum above 20 GeV, in a simulated dataset of top quark pairs,
at center-of-mass energy of 13 TeV.

A similar architecture to the Particle Flow Networks was used in the ATLAS
experiment as an alternative to the RNN-based model in the context of heavy
flavor identification ATL-PHYS-PUB-2020-014 . The Deep Sets-based ATLAS heavy
flavor discriminant (DIPS) outperforms its RNN analogous (RNNIP, described in
section 2.1) version of the RNNIP ATLAS heavy flavor discriminant by up to
$20\%$ for $b$-jet efficiencies larger than $60\%$, while using the same
inputs, as seen in Figure 18. DIPS has also been found to reduce significantly
the training and evaluation time with respect to RNNIP, due to its
paralelizability. The ability of parallelizing computations on each sequence
element is an important feature of this model, particularly for applications
in which the network evaluation time is limited.

##  4 Conclusion

Sequence-based Machine Learning algorithms have a long and rich history in the
context of natural language processing. While the idea of representing a jet
as a sequence of its constituents is not new in particle physics, as evidenced
by jet clustering algorithms, the usage of ML concepts exploting this
representation is relatively recent when compared to computer vision
algorithms. Even so, the application of RNNs to jet physics in particular has
been a fruitful avenue of research in the past few years. Special attention
was given to the predictive power of these models, with the successful
application of LSTM-based neural network architectures to different types of
jet classification tasks. More generally, recurrent structures were shown to
be well-suited to describe jet clustering histories, leading to a full
probabilistic model of a jet given its constituents.

Recent work has also been focused on expanding the basic RNN ideas of
hierarchical context learning to more physics inspired architectures, such as
recursive neural networks. However, while these models achieve high precision
when the correct choice of input structure is used - either the binary tree in
a RecNN or the ordered sequence itself for RNNs - their performances can be
significantly degraded given a wrong structure choice. This is particularly
difficult to deal with when the input has variable sized length but the data
structure is not obvious. For example, how one chooses to order the set of
tracks inside the jet will depend on the task to be performed: ordering in
impact parameter significance can be suited for b-jet identification, but not
for quark versus gluon discrimination. With that in mind, algorithims in which
the data structure itself is either learned, such as transformers or graphs,
or invariant under certain problem transformations, such as permutation
invariant sets, have shown a great potential for future studies.

## Acknowledgements

RTdL would like to thank Michael Kagan for the help in reviewing the document;
the book editors, Paolo Calafiura, David Rousseau and Kazuhiro Terao; and the
other reviewers involved. This work was supported by the US Department of
Energy (DOE) under grant DE-AC02-76SF00515, and by the SLAC Panofsky
Fellowship.

## References

  * (1) B. Hammer, “On the approximation capability of recurrent neural networks,” Neurocomputing, vol. 31, no. 1, pp. 107 – 123, 2000.
  * (2) M. Schuster and K. Paliwal, “Bidirectional recurrent neural networks,” Signal Processing, IEEE Transactions on, vol. 45, pp. 2673 – 2681, 12 1997.
  * (3) R. Pascanu, T. Mikolov, and Y. Bengio, “Understanding the exploding gradient problem,” CoRR, vol. abs/1211.5063, 2012.
  * (4) S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation, vol. 9, pp. 1735–80, 12 1997.
  * (5) K. Cho, B. van Merrienboer, Ç. Gülçehre, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder-decoder for statistical machine translation,” CoRR, vol. abs/1406.1078, 2014.
  * (6) C. Olah, “Understanding LSTM Networks,” 2015.  https://colah.github.io/posts/2015-08-Understanding-LSTMs/.
  * (7) L. de Oliveira, M. Kagan, L. Mackey, B. Nachman, and A. Schwartzman, “Jet-images — deep learning edition,” JHEP, vol. 07, p. 069, 2016\.
  * (8) M. Cacciari, G. P. Salam, and G. Soyez, “The anti-$k_{t}$ jet clustering algorithm,” JHEP, vol. 04, p. 063, 2008.
  * (9) LHC Higgs Cross Section Working Group, S. Heinemeyer, C. Mariotti, G. Passarino, and R. Tanaka (Eds.), “Handbook of LHC Higgs Cross Sections: 3\. Higgs Properties,” CERN-2013-004, CERN, Geneva, 2013.
  * (10) D. Buttazzo, G. Degrassi, P. P. Giardino, G. F. Giudice, F. Sala, A. Salvio, and A. Strumia, “Investigating the near-criticality of the higgs boson,” Journal of High Energy Physics, vol. 2013, Dec 2013.
  * (11) P. Zyla et al., “Review of Particle Physics,” PTEP, vol. 2020, no. 8, p. 083C01, 2020.
  * (12) “Standard Model Summary Plots Spring 2019,” Tech. Rep. ATL-PHYS-PUB-2019-010, CERN, Geneva, Mar 2019.
  * (13) “Topological $b$-hadron decay reconstruction and identification of $b$-jets with the JetFitter package in the ATLAS experiment at the LHC,” Tech. Rep. ATL-PHYS-PUB-2018-025, CERN, Geneva, Oct 2018.
  * (14) G. Aad, B. Abbott, D. C. Abbott, A. A. Abud, K. Abeling, D. K. Abhayasinghe, S. H. Abidi, O. S. AbouZeid, N. L. Abraham, and et al., “ATLAS b-jet identification performance and efficiency measurement with $t\bar{t}$ events in pp collisions at $\sqrt{s}=13$ TeV,” The European Physical Journal C, vol. 79, Nov 2019.
  * (15) D. Guest, J. Collado, P. Baldi, S.-C. Hsu, G. Urban, and D. Whiteson, “Jet Flavor Classification in High-Energy Physics with Deep Neural Networks,” Phys. Rev., vol. D94, no. 11, p. 112002, 2016.
  * (16) ATLAS Collaboration, “Identification of Jets Containing $b$-Hadrons with Recurrent Neural Networks at the ATLAS Experiment,” Tech. Rep. ATL-PHYS-PUB-2017-003, CERN, Geneva, Mar 2017.
  * (17) A. Sirunyan, A. Tumasyan, W. Adam, F. Ambrogi, E. Asilar, T. Bergauer, J. Brandstetter, E. Brondolin, M. Dragicevic, J. Erö, and et al., “Identification of heavy-flavour jets with the cms detector in pp collisions at 13 tev,” Journal of Instrumentation, vol. 13, p. P05011–P05011, May 2018.
  * (18) “Performance of b tagging algorithms in proton-proton collisions at 13 TeV with Phase 1 CMS detector,” Jun 2018.
  * (19) A. Ali, F. Barreiro, and T. Lagouri, “Prospects of measuring the CKM matrix element $|V_{ts}|$ at the LHC,” Physics Letters B, vol. 693, p. 44–51, Sep 2010.
  * (20) J. Erdmann, “A tagger for strange jets based on tracking information using long short-term memory,” JINST, vol. 15, no. 01, p. P01021, 2020.
  * (21) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi, “Delphes 3: a modular framework for fast simulation of a generic collider experiment,” Journal of High Energy Physics, vol. 2014, Feb 2014.
  * (22) ATLAS Collaboration, “Identification of hadronic tau lepton decays using neural networks in the ATLAS experiment,” Tech. Rep. ATL-PHYS-PUB-2019-033, CERN, Geneva, Aug 2019.
  * (23) A. Sirunyan, A. Tumasyan, W. Adam, F. Ambrogi, T. Bergauer, M. Dragicevic, J. Erö, A. E. D. Valle, M. Flechl, R. Frühwirth, and et al., “Identification of heavy, energetic, hadronically decaying particles using machine-learning techniques,” Journal of Instrumentation, vol. 15, p. P06005–P06005, Jun 2020.
  * (24) M. Wobisch and T. Wengler, “Hadronization corrections to jet cross-sections in deep inelastic scattering,” in Workshop on Monte Carlo Generators for HERA Physics (Plenary Starting Meeting), pp. 270–279, 4 1998.
  * (25) M. Aaboud et al., “Performance of top-quark and $W$-boson tagging with ATLAS in Run 2 of the LHC,” Eur. Phys. J. C, vol. 79, no. 5, p. 375, 2019\.
  * (26) S. Egan, W. Fedorko, A. Lister, J. Pearkes, and C. Gay, “Long Short-Term Memory (LSTM) networks with jet constituents for boosted top tagging at the LHC,” 2017.
  * (27) J. Pearkes, W. Fedorko, A. Lister, and C. Gay, “Jet constituents for deep neural network based top quark tagging,” 2017.
  * (28) D. Krohn, J. Thaler, and L.-T. Wang, “Jet trimming,” Journal of High Energy Physics, vol. 2010, Feb 2010.
  * (29) A. Andreassen, I. Feige, C. Frye, and M. D. Schwartz, “JUNIPR: a Framework for Unsupervised Machine Learning in Particle Physics,” Eur. Phys. J., vol. C79, no. 2, p. 102, 2019.
  * (30) S. Catani, Y. Dokshitzer, M. Seymour, and B. Webber, “"longitudinally-invariant $k_{t}$-clustering algorithms for hadron-hadron collisions",” Nuclear Physics B, vol. 406, no. 1, pp. 187 – 224, 1993\.
  * (31) A. Andreassen, I. Feige, C. Frye, and M. D. Schwartz, “Binary JUNIPR: an interpretable probabilistic model for discrimination,” Phys. Rev. Lett., vol. 123, no. 18, p. 182001, 2019.
  * (32) L. Bottou, “From machine learning to machine reasoning,” CoRR, vol. abs/1102.1808, 2011.
  * (33) G. Louppe, K. Cho, C. Becot, and K. Cranmer, “QCD-Aware Recursive Neural Networks for Jet Physics,” JHEP, vol. 01, p. 057, 2019.
  * (34) T. Cheng, “Recursive Neural Networks in Quark/Gluon Tagging,” Comput. Softw. Big Sci., vol. 2, no. 1, p. 3, 2018.
  * (35) K. Fraser and M. D. Schwartz, “Jet Charge and Machine Learning,” JHEP, vol. 10, p. 093, 2018.
  * (36) D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by jointly learning to align and translate,” in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (Y. Bengio and Y. LeCun, eds.), 2015.
  * (37) A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” CoRR, vol. abs/1706.03762, 2017.
  * (38) J. Devlin, M. Chang, K. Lee, and K. Toutanova, “BERT: pre-training of deep bidirectional transformers for language understanding,” in Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL-HLT 2019, Minneapolis, MN, USA, June 2-7, 2019, Volume 1 (Long and Short Papers) (J. Burstein, C. Doran, and T. Solorio, eds.), pp. 4171–4186, Association for Computational Linguistics, 2019.
  * (39) M. Zaheer, S. Kottur, S. Ravanbakhsh, B. Póczos, R. Salakhutdinov, and A. J. Smola, “Deep sets,” CoRR, vol. abs/1703.06114, 2017.
  * (40) P. T. Komiske, E. M. Metodiev, and J. Thaler, “Energy flow networks: deep sets for particle jets,” Journal of High Energy Physics, vol. 2019, Jan 2019\.
  * (41) “Deep Sets based Neural Networks for Impact Parameter Flavour Tagging in ATLAS,” Tech. Rep. ATL-PHYS-PUB-2020-014, CERN, Geneva, May 2020.


††institutetext: 1 Department of Physics, Yantai University, Yantai 264005, P.
R. China
2 School of Physics and Microelectronics, Zhengzhou University, ZhengZhou
450001, P. R. China
3 CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, P. R. China
4 School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, P. R. China

#
A joint explanation of W-mass and muon g-2 in 2HDM

Xiao-Fang Han1, Fei Wang2, Lei Wang1, Jin Min Yang3,4, Yang Zhang2

###### Abstract

Since both $W$-mass and muon $g-2$ can be affected by the mass splittings
among extra Higgs bosons $(H,~{}A,~{}H^{\pm})$ in a 2HDM, we take a model with
$\mu$-$\tau$ LFV interactions to examine the two anomalies reported
respectively by CDF II and FNAL. We obtain the following observations: (i)
Combined with theoretical constraints, the CDF $W$-mass measurement disfavors
$H$ or $A$ to degenerate in mass with $H^{\pm}$, but allows $H$ and $A$ to
degenerate. The mass splitting between $H^{\pm}$ and $H/A$ is required to be
larger than 10 GeV. The $m_{H^{\pm}}$ and $m_{A}$ are favored to be smaller
than 650 GeV for $m_{H}<120$ GeV, and allowed to have more large values with
increasing of $m_{H}$. (ii) After imposing other relevant experimental
constraints, there are parameter spaces that simultaneously satisfy (at
$2\sigma$ level) the CDF $W$-mass, the FNAL muon $g-2$ and the data of lepton
universality in $\tau$ decays, but the mass splittings among extra Higgs
bosons are strictly constrained.

##  1 Introduction

The CDF collaboration presented their new result for the $W$-boson mass
measurement cdfmw

$\displaystyle\Delta m_{W}=80.4335\pm 0.0094{\rm GeV}.$ (1)

The experimental central value has an approximate $7\sigma$ discrepancy from
the Standard Model (SM) prediction, $80.357\pm 0.006$ GeVsmmw . On the other
hand, there has been a long-standing discrepancy between the SM prediction and
experiment for the muon anomalous magnetic moment (muon g-2). The combined
result of the FNAL E989 experiment mug2-exp  and the BNL experiment fermig2
has an approximate $4.2\sigma$ deviation from the SM prediction smg2-1 ;
smg2-2 ; smg2-3

$\displaystyle\Delta a_{\mu}=a_{\mu}^{exp}-a_{\mu}^{SM}=(25.1\pm 5.9)\times
10^{-10}.$ (2)

Both the deviations strongly imply existence of new physics beyond SM. Some
plausible explanations have already performed for the CDF $W$-mass 2204.03693
; 2204.03796 ; 2204.03996 ; 2204.04183 ; 2204.04191 ; 2204.04202 ; 2204.04204
; 2204.04286 ; 2204.04356 ; 2204.04514 ; 2204.04559 ; 2204.04770 ; 2204.04805
; 2204.04834 ; 2204.05024 ; 2204.05031 ; 2204.05085 ; 2204.05260 ; 2204.05267
; 2204.05269 ; 2204.05283 ; 2204.05284 ; 2204.05285 ; 2204.05296 ; 2204.05302
; 2204.05303 ; 2204.04688 ; 2204.05728 ; 2204.05760 ; 2204.05942 ; 2204.05962
; 2204.05965 ; 2204.05975 ; 2204.05992 .

Among various new physics models, two-Higgs-doublet models (2HDMs) are rather
simple extensions of the SM (for a recent review, see, e.g., stu4 ). A 2HDM
introduces a second $SU(2)_{L}$ Higgs doublet and thus predicts two neutral
CP-even Higgs bosons $h$ and $H$, one neutral pseudoscalar $A$ and a pair of
charged Higgs boson $H^{\pm}$ t.lee . This kind of models can give additional
corrections to the masses of gauge bosons via the self-energy diagrams
exchanging extra Higgs fields. In addition, if the extra Higgs bosons have
appropriate couplings to the leptons, the muon $g-2$ can be simply explained.
Since both the $W$-mass and muon $g-2$ can be affected by the mass splittings
among $H$, $A$, and $H^{\pm}$, we in this note take a 2HDM with the
$\mu$-$\tau$ lepton flavor violation (LFV) interactions to study the
possibility of a simultaneous explanation of both anomalies. In our analysis
we will intensively examine the parameter space of this model by considering
various relevant theoretical and experimental constraints. For the single
explanation of muon g-2 using the Higgs doublet field with the $\mu$-$\tau$
LFV interactions, see e.g., taumug2-3 ; 0207302 ; 10010434 ; 150207824 ;
151108544 ; 151108880 ; 1610.06587 ; 160604408 ; 190410908 ; 1907.09845 ;
1908.09845 ; 2104.03242 ; lfv11 ; lfv12 ; lfv13 ; 1711.08430 ; 2010.04266 .

The paper is organized as follows. In Sec. II we will introduce the 2HDM with
the $\mu$-$\tau$ LFV interactions. In Sec. III and Sec. IV we study the
W-boson mass and muon g-2 anomaly, considering relevant theoretical and
experimental constraints. Finally, we give our conclusion in Sec. V.

##  2 A two-Higgs-doublet model with $\mu$-$\tau$ LFV interactions

The 2HDM with the $\mu$-$\tau$ LFV interactions may be derived from a general
2HDM by take specific parameters. Also it can be naturally obtained by
introducing an inert Higgs doublet $\phi_{2}$ under a discrete $Z_{4}$
symmetry, and the $Z_{4}$ charge assignment is displaced in Table I 190410908
.

Table 1: Assignment of $Z_{4}$ charge in the 2HDM with $\mu-\tau$-philic Higgs doublet. |  $\phi_{1}$ |  $\phi_{2}$ |  $Q_{L}^{i}$ |  $U_{R}^{i}$ |  $D_{R}^{i}$ |  $L_{L}^{e}$ |  $L_{L}^{\mu}$ |  $L_{L}^{\tau}$ |  $e_{R}$ |  $\mu_{R}$ |  $\tau_{R}$
---|---|---|---|---|---|---|---|---|---|---|---
Z 4 | $1$ | -1 | 1 | 1 | 1 | 1 | $i$ | $-i$ | 1 | $i$ | $-i$

The Higgs potential with the $Z_{4}$ symmetry is given by

$\displaystyle\mathrm{V}$ $\displaystyle=$ $\displaystyle
Y_{1}(\phi_{1}^{\dagger}\phi_{1})+Y_{2}(\phi_{2}^{\dagger}\phi_{2})+\frac{\lambda_{1}}{2}(\phi_{1}^{\dagger}\phi_{1})^{2}+\frac{\lambda_{2}}{2}(\phi_{2}^{\dagger}\phi_{2})^{2}$
(3)
$\displaystyle+\lambda_{3}(\phi_{1}^{\dagger}\phi_{1})(\phi_{2}^{\dagger}\phi_{2})+\lambda_{4}(\phi_{1}^{\dagger}\phi_{2})(\phi_{2}^{\dagger}\phi_{1})+\left[\frac{\lambda_{5}}{2}(\phi_{1}^{\dagger}\phi_{2})^{2}+\rm
h.c.\right].$

Here all the parameters are real. Although $\lambda_{5}$ is the only possible
complex parameter, it can be rendered real with a phase redefinition of one of
the two Higgs fields. The two Higgs doublets $\phi_{1}$ and $\phi_{2}$ are
expressed by

$\phi_{1}=\left(\begin{array}[]{c}G^{+}\\\
\frac{1}{\sqrt{2}}\,(v+h+iG^{0})\end{array}\right)\,,\ \ \
\phi_{2}=\left(\begin{array}[]{c}H^{+}\\\
\frac{1}{\sqrt{2}}\,(H+iA)\end{array}\right).$

The $\phi_{1}$ field obtains a nonzero vacuum expectation value (VEV), $v$=246
GeV, while the $\phi_{2}$ field has zero VEV. The parameter $Y_{1}$ can be
determined from the minimization condition for the Higgs potential,

$Y_{1}=-\frac{1}{2}\lambda_{1}v^{2}.$ (4)

The fields $G^{0}$ and $G^{+}$ indicate Nambu-Goldstone bosons eaten by the
gauge bosons. The fields $A$ and $H^{+}$ represent the mass eigenstates of the
CP-odd Higgs boson and charged Higgs boson, whose masses are written as

$m_{H^{\pm}}^{2}=Y_{2}+\frac{\lambda_{3}}{2}v^{2},~{}~{}~{}m_{A}^{2}=m_{H^{\pm}}^{2}+\frac{1}{2}(\lambda_{4}-\lambda_{5})v^{2}.$
(5)

There is no mixing between the two CP-even Higgs bosons $h$ and $H$, and their
masses are

$m_{h}^{2}=\lambda_{1}v^{2}\equiv(125~{}{\rm
GeV})^{2},~{}~{}~{}m_{H}^{2}=m_{A}^{2}+\lambda_{5}v^{2}.$ (6)

We obtain the masses of fermions via the Yukawa interactions with $\phi_{1}$,

$-{\cal
L}=y_{u}\overline{Q}_{L}\,\tilde{{\phi}}_{1}\,U_{R}+y_{d}\overline{Q}_{L}\,{\phi}_{1}\,D_{R}+y_{\ell}\overline{L}_{L}\,{\phi}_{1}\,E_{R}+\mbox{h.c.},$
(7)

where $E_{R}$, $U_{R}$ and $D_{R}$ stand for the three generations of right-
handed fermion fields for the charged leptons, up-type quarks and down-type
quarks, and we define $L_{L}=(v_{L_{i}},\ell_{L_{i}})^{T}$,
$Q_{L}=(u_{L_{i}},d_{L_{i}})^{T}$ and
$\widetilde{\Phi}_{1}=i\tau_{2}\Phi^{*}_{1}$ with $i$ representing generation
indices. Under the $Z_{4}$ symmetry, the lepton Yukawa matrix $y_{\ell}$ is
diagonal and therefore the lepton fields ($L_{L}$, $E_{R}$) are mass
eigenstates.

Under the $Z_{4}$ symmetry, the $\phi_{2}$ doublet is allowed to have
$\mu$-$\tau$ interactions 190410908 ,

$\displaystyle-{\cal L}_{LFV}$ $\displaystyle=$
$\displaystyle\sqrt{2}~{}\rho_{\mu\tau}\,\overline{L^{\mu}_{L}}\,{\phi}_{2}\,\tau_{R}\,+\sqrt{2}~{}\rho_{\tau\mu}\,\overline{L^{\tau}_{L}}\,{\phi}_{2}\,\mu_{R}\,+\,\mbox{h.c.}\,.$
(8)

The additional Higgs bosons $H$, $A$ and $H^{\pm}$ only have $\mu-\tau$
flavor-violating Yukawa couplings. On the other hand, the SM-like Higgs boson
$h$ has exactly same couplings to the gauge bosons and fermions as the SM
Higgs, with no $\mu-\tau$ flavor-violating couplings at the tree level.

##  3 The $S,~{}T,~{}U$ parameters and $W$-mass

The model can produce main corrections to the masses of gauge bosons via self-
energy diagrams exchanging extra Higgs fields. The oblique parameters
$(S,~{}T,~{}U)$ stu ; w-stu  represent radiative corrections to the two-point
functions of gauge bosons. Most effects on precision measurements can be
described by these parameters. Recently, Ref. 2204.03796  gave the values of
these parameters from an analysis of precision electroweak data including the
CDF new result of the $W$-mass,

$S=0.06\pm 0.10,~{}~{}T=0.11\pm 0.12,~{}~{}U=0.14\pm 0.09.$ (9)

The correlation coefficients are given by

$\rho_{ST}=0.9,~{}~{}\rho_{SU}=-0.59,~{}~{}\rho_{TU}=-0.85.$ (10)

The $W$-boson mass can be inferred from the following relation w-stu ,

$\Delta m_{W}^{2}=\frac{\alpha
c_{W}^{2}}{c_{W}^{2}-s_{W}^{2}}m_{Z}^{2}(-\frac{1}{2}S+c_{W}^{2}T+\frac{c_{W}^{2}-s_{W}^{2}}{4s_{W}^{2}}U).$
(11)

In our analysis, we adopt 2HDMC 2hc-1  to calculate the 2HDM corrections to
$S,~{}T,~{}U$ parameters, and perform a global fit to the predictions of
$S,~{}T,~{}U$ parameters. As the global fit results are presented on two-
dimension planes, a limit of $\chi^{2}<\chi^{2}_{\rm min}+6.18$ is set to
obtain profiled $2\sigma$ favoured regions, where $\chi^{2}_{\rm min}$ is the
minimum of $\chi^{2}$ corresponding the best fit point.

We scan $m_{H}$, $m_{A}$ and $m_{H^{\pm}}$ parameters in the following ranges:

$\displaystyle
80{\rm~{}GeV}<m_{H^{\pm}}<1000{\rm~{}GeV},~{}~{}65{\rm~{}GeV}<m_{A}<1000{\rm~{}GeV},$
$\displaystyle
10{\rm~{}GeV}<m_{H}<120{\rm~{}GeV},~{}~{}130{\rm~{}GeV}<m_{H}<1000{\rm~{}GeV}.$
(12)

Figure 1: For $m_{H}>$ 130 GeV, the samples explaining the CDF II results of
the $W$-mass within $2\sigma$ range while satisfying the constraints of the
oblique parameters and theoretical constraints. Figure 2: Same as Fig. 1, but
for $m_{H}<$ 120 GeV.

In Fig. 1 and Fig. 2, we show the samples explaining the CDF II $W$-boson mass
measurement within $2\sigma$ range while satisfying the constraints of the
oblique parameters and theoretical constraints from the perturbativity, vacuum
stability and unitarity. From Fig. 1 and Fig. 2, we see that $H$ or $A$ is
disfavored to exactly degenerate in mass with $H^{\pm}$, but their masses are
allowed to be degenerate. The mass splitting between $H^{\pm}$ and $H/A$ is
imposed upper and lower bounds, and required to be approximately larger than
10 GeV. When one of $m_{H}$ and $m_{A}$ is close to $m_{H^{\pm}}$, the other
is allowed to have sizable deviation from $m_{H^{\pm}}$. The $m_{H^{\pm}}$ and
$m_{A}$ are favored to be smaller than 650 GeV for $m_{H}<120$ GeV (see Fig.
2), and allowed to have more large values with increasing of $m_{H}$ (see Fig.
1).

Now we analyze the reason. In the 2HDM, the corrections to $T$ is much larger
than $S$ and $U$, and the $T$ parameter is expressed as

$\displaystyle T\;$ $\displaystyle=$ $\displaystyle\;\frac{1}{16\pi
M_{W}^{2}s_{W}^{2}}\,\Biggl{\\{}\biggl{[}\mathcal{F}(M_{H^{\pm}}^{2},M_{H}^{2})-\mathcal{F}(M_{H}^{2},M_{A}^{2})+\mathcal{F}(M_{H^{\pm}}^{2},M_{A}^{2})\biggr{]}\;\Biggr{\\}}\,,$
(13)

where the $\mathcal{F}$ function is stu1 ; stu2 ; stu3 .

$\mathcal{F}(m_{1}^{2},m_{2}^{2})=\frac{1}{2}\,(m_{1}^{2}+m_{2}^{2})-\frac{m_{1}^{2}m_{2}^{2}}{m_{1}^{2}-m_{2}^{2}}\;\log{\left(\frac{m_{1}^{2}}{m_{2}^{2}}\right)}\,.$
(14)

This expression shows that the parameters $T$ is sensitive to the mass
splittings among $H,~{}A$ and $H^{\pm}$. The $T$ parameter will be zero for
$m_{H}=m_{H^{\pm}}$ or $m_{A}=m_{H^{\pm}}$, but takes a non-zero value for
$m_{A}=m_{H}$. Therefore, the corrections of 2HDM to the oblique parameters
tend to become small as one of $m_{H}$ and $m_{A}$ approaches to
$m_{H^{\pm}}$. However, in order to accommodate the $W$-mass reported by the
CDF II collaboration, the 2HDM needs to produce an appropriate value of $T$,
which excludes $m_{H}=m_{H^{\pm}}$ or $m_{A}=m_{H^{\pm}}$.

##  4 Muon $g-2$, $\tau$ decays, and other relevant constraints

The model gives additional corrections to the muon $g-2$ anomaly ($\Delta
a_{\mu}$) via the one-loop diagrams involving the $\mu$-$\tau$ LFV couplings
of $H$ and $A$ taumug2-3 ; 0207302 ; 10010434 ,

$\displaystyle\Delta
a_{\mu}=\frac{m_{\mu}m_{\tau}\rho^{2}}{8\pi^{2}}\left[\frac{(\log\frac{m_{H}^{2}}{m_{\tau}^{2}}-\frac{3}{2})}{m_{H}^{2}}-\frac{\log(\frac{m_{A}^{2}}{m_{\tau}^{2}}-\frac{3}{2})}{m_{A}^{2}}\right].$
(15)

Here we find $\Delta a_{\mu}>0$ for $m_{A}>m_{H}$.

Because the extra Higgs bosons have the $\mu$-$\tau$ LFV interactions, the
model can affect the lepton flavor universality (LFU) in the $\tau$ lepton
decays. The HFAG collaboration tests the LFU from ratios of the partial widths
of a heavier lepton. They obtains tauexp

$\displaystyle\left(\frac{g_{\tau}}{g_{\mu}}\right)=1.0011\pm
0.0015,\left(\frac{g_{\tau}}{g_{e}}\right)=1.0029\pm
0.0015,\left(\frac{g_{\mu}}{g_{e}}\right)=1.0018\pm 0.0014,$ (16)

using pure leptonic processes, namely

$\displaystyle\left(\frac{g_{\tau}}{g_{\mu}}\right)^{2}$
$\displaystyle\equiv\frac{\overline{\Gamma}(\tau\to
e\nu\overline{\nu})}{\overline{\Gamma}(\mu\to e\nu\overline{\nu}},$ (17)
$\displaystyle\left(\frac{g_{\tau}}{g_{e}}\right)$
$\displaystyle\equiv\frac{\overline{\Gamma}(\tau\to\mu\nu\overline{\nu})}{\overline{\Gamma}(\mu\to
e\nu\overline{\nu}},$ (18) $\displaystyle\left(\frac{g_{\mu}}{g_{e}}\right)$
$\displaystyle\equiv\frac{\overline{\Gamma}(\tau\to\mu\nu\overline{\nu})}{\overline{\Gamma}(\tau\to
e\nu\overline{\nu}},$ (19)

with $\overline{\Gamma}$ representing the partial width which is normalized by
the corresponding SM value. With the two semi-hadronic processes,

$\displaystyle\left(\frac{g_{\tau}}{g_{\mu}}\right)_{h}^{2}$
$\displaystyle\equiv\frac{{\rm Br}(\tau\to h\nu)}{{\rm
Br}(h\to\mu\overline{\nu})}\frac{2m_{h}m_{\mu}^{2}\tau_{h}}{(1+\delta_{h})m_{\tau}^{2}\tau_{\tau}}\left(\frac{1-m_{\mu}^{2}/m_{h}^{2}}{1-m_{h}^{2}/m_{\tau}^{2}}\right)^{2},$
(20)

where $h$ indicates $\pi$ or $K$, they measure

$\displaystyle\left(\frac{g_{\tau}}{g_{\mu}}\right)_{\pi}=0.9963\pm
0.0027,~{}~{}~{}~{}~{}\left(\frac{g_{\tau}}{g_{\mu}}\right)_{K}=0.9858\pm
0.0071.$ (21)

The statistical correlation matrix for the five fitted coupling ratios is

$\displaystyle\left[\begin{array}[]{ccccc}1&53\%&-49\%&24\%&12\%\\\
53\%&1&48\%&26\%&10\%\\\ -49\%&48\%&1&2\%&-2\%\\\ 24\%&26\%&2\%&1&5\%\\\
12\%&10\%&-2\%&5\%&1\end{array}\right]$ (27)

In this model, we can calculate

$\displaystyle\bar{\Gamma}(\tau\to\mu\nu\bar{\nu})=(1+\delta_{\rm
loop}^{\tau})^{2}~{}(1+\delta_{\rm loop}^{\mu})^{2}+\delta_{\rm tree},$
$\displaystyle\bar{\Gamma}(\tau\to e\nu\bar{\nu})=(1+\delta_{\rm
loop}^{\tau})^{2},$ $\displaystyle\bar{\Gamma}(\mu\to
e\nu\bar{\nu})=(1+\delta_{\rm loop}^{\mu})^{2},$ (28)

where the tree-level correction $\delta_{\rm tree}$ is from the charged Higgs
contribution to $\tau\to\mu\nu\overline{\nu}$,

$\delta_{\rm tree}=4\frac{m_{W}^{4}\rho^{4}}{g^{4}m_{H^{\pm}}^{4}},$ (29)

which can give a positive correction to $\tau\to\mu\nu\bar{\nu}$. $\delta_{\rm
loop}^{\mu}$ and $\delta_{\rm loop}^{\tau}$ are the corrections to vertices
$W\bar{\nu_{\mu}}\mu$ and $W\bar{\nu_{\tau}}\tau$ from the one-loop diagrams
containing $A$, $H$, and $H^{\pm}$, respectively. As we assume
$\rho_{\mu\tau}=\rho_{\tau\mu}$ in the lepton Yukawa matrix, the two loop
corrections are identical 190410908 ; tavv-1 ; mu2h16 ,

$\delta_{\rm loop}^{\tau}=\delta_{\rm loop}^{\mu}={1\over
16\pi^{2}}{\rho^{2}}\left[1+{1\over
4}\left(H(x_{A})+H(x_{H})\right)\right]\,,$ (30)

where $H(x_{\phi})\equiv\ln(x_{\phi})(1+x_{\phi})/(1-x_{\phi})$ with
$x_{\phi}=m_{\phi}^{2}/m_{H^{\pm}}^{2}$. Meanwhile, for the semi-hadronic
processes, we have

$\left(g_{\tau}\over g_{\mu}\right)=\left(g_{\tau}\over
g_{\mu}\right)_{K}=\left(g_{\tau}\over g_{\mu}\right)_{\pi}.$ (31)

In our study, we perform a fit for predictions of these five ratios. Note that
there is a vanishing eigenvalue in the covariance matrix constructed from Eq.
(16), Eq. (21) and Eq. (27), and therefore such a degree of freedom is
removed. With remaining four degrees of freedom, the $2\sigma$ confidence
level region is obtained by adopting a limit of $\chi^{2}_{\tau}<9.72$. Thus,
the surviving samples are much more consistent with the experimental results
than the SM, which has a $\chi^{2}_{\tau}$ of $12.25$. Also the model can
affect the LFU in the $Z$-boson decays zexp . The constraints from the
$Z$-boson decays are generally weaker than those from the $\tau$ decays.

In the model, the fields $H$, $A$ and $H^{\pm}$ have no couplings to quarks,
and therefore they are produced at the LHC mainly via the electroweak
processes, $pp\to W^{\pm*}\to H^{\pm}A/H$, $pp\to Z^{*}\to HA$ and $pp\to
Z^{*}/\gamma^{*}\to H^{+}H^{-}$. The final state signal mainly includes multi-
leptons, and therefore the multi-lepton event searches at the LHC can impose
stringent constraints which require $m_{H}$ to be larger than 560 GeV
1908.09845 . Also a very light $H$ may escape the constraints of the direct
searches at the LHC 2104.03242 . In this case, the explanation of muon $g-2$
anomaly requires a very small $\rho$, which leads that the contributions of
the model to the $\tau$ decays are too small to explain the data of LFU in the
$\tau$ decays. Therefore, in this work we discuss the scenario of $m_{H}>560$
GeV. In addition, to respect the perturbativity we choose the Yukawa coupling
parameter $\rho<1$.

At the tree level, the 125 GeV Higgs has the same couplings to the SM
particles as in the SM. The $h\to\gamma\gamma$ decay will be corrected by the
one-loop diagram of the charged Higgs. We consider the bound of the diphoton
signal strength smmw ,

$\mu_{\gamma\gamma}=1.11^{+0.1}_{-0.09}.$ (32)

Figure 3: The bullets can explain the muon $g-2$ within $2\sigma$ range, while
the squares can explain the muon $g-2$ and the $\tau$-decays within $2\sigma$
ranges while satisfying the constraints of $Z$-decays. Other relevant
constraints from the theory, the oblique parameters, the CDF II $W$-mass and
the diphoton signal data of the 125 GeV Higgs are also satisfied).

In Fig. 3, we project the surviving samples explaining the muon $g-2$ anomaly
and the LFU in $\tau$-decays within $2\sigma$ ranges, with other relevant
constraints from the theory, the oblique parameters, the CDF II $W$-mass and
the diphoton signal data of the 125 GeV Higgs being satisfied. Eq. (15) shows
that $H$ and $A$ respectively give positive and negative contributions to the
muon $g-2$, and their contributions are suppressed by their masses. Therefore,
the explanation of muon $g-2$ requires $m_{A}>m_{H}$ and $\rho$ to become
large with increasing of $m_{H}$ and decreasing of $(m_{A}-m_{H})$ (see the
upper-left panel).

The ratio $\left(g_{\tau}\over g_{e}\right)$ in the $\tau$ decays has
approximately $2\sigma$ deviation from the SM. Therefore, enhancing
$\Gamma(\tau\to\mu\nu\bar{\nu})$ may give a better fit to the data of LFU in
the $\tau$ decays. The decay $\tau\to\mu\nu\nu$ obtains a positive correction
from the $\delta_{\rm tree}$ term of Eq. (4). Such a term is from the tree-
level diagram mediated by the charged Higgs, and proportional to
$\rho^{4}/m^{2}_{H^{\pm}}$. Therefore, the upper-right panel of Fig. 3 shows
that the value of $\chi_{\tau}^{2}$ tends to become large with increasing of
$m_{H^{\pm}}$ and decreasing of $\rho$. In addition, because of the
constraints of the oblique parameters and $W$-boson mass, the upper-left panel
of Fig. 1 shows that $m_{H^{\pm}}$ tends to increase with $m_{H}$, especially
for a large $m_{H}$, which implies that the value of $\chi_{\tau}^{2}$ becomes
large with increase of $m_{H}$. Therefore, the upper-left panel of Fig. 3
shows that a simultaneous explanation of the muon $g-2$ and LFU in $\tau$
decays favors a large $\rho$ which increases with $m_{H}$.

From the lower panel of Fig. 3, we see that the mass splittings among $H$, $A$
and $H^{\pm}$ are stringently constrained in the region simultaneously
explaining the $W$-mass, muon $g-2$ and LFU in $\tau$ decays, i.e., 10 GeV
$<m_{A}-m_{H}<$ 75 GeV, 65 GeV $<m_{H^{\pm}}-m_{A}<$ 100 GeV, 85 GeV
$<m_{H^{\pm}}-m_{H}<$ 125 GeV ( -150 GeV $<m_{H^{\pm}}-m_{A}<$ -85 GeV, -105
GeV $<m_{H^{\pm}}-m_{H}<$ -55 GeV).

##  5 Conclusion

We examine the CDF II $W$-boson mass and the FNAL muon $g-2$ in the 2HDM in
which the extra Higgs doublet has the $\mu$-$\tau$ LFV interactions. Imposing
the theoretical constraints, we found that the CDF II $W$-boson mass disfavors
$H$ or $A$ to degenerate in mass with $H^{\pm}$, and the mass splitting
between $H^{\pm}$ and $H/A$ is favored to be larger than 10 GeV. The
$m_{H^{\pm}}$ and $m_{A}$ are favored to be smaller than 650 GeV for
$m_{H}<120$ GeV, and allowed to have more large values with increasing of
$m_{H}$. Considering other relevant experimental constraints, we found that
the mass splittings among $H,~{}A$ and $H^{\pm}$ are stringently restricted in
the parameter space which can simultaneously explain the CDF II $W$-mass, the
FNAL muon $g-2$, and the data of LFU in the $\tau$ decays.

## Acknowledgment

We thank Lei Wu for helpful discussions. This work was supported by the
National Natural Science Foundation of China under grants 11975013, 12105248,
11821505, 12075300, 12075213, by Peng-Huan-Wu Theoretical Physics Innovation
Center (12047503), by the CAS Center for Excellence in Particle Physics
(CCEPP), and by the Key Research Program of the Chinese Academy of Sciences,
Grant NO. XDPB15.

## References

  * (1) CDF Collaboration, Science ${\bf 376}$, 6589 (2022).
  * (2) P. A. Zyla et al. [Particle Data Group], “Review of Particle Physics,” PTEP 2020, (2020) 083C01.
  * (3) Muon g-2 Collaboration, Phys. Rev. Lett. 86, (2001) 2227; Phys. Rev. D 73, (2006) 072003.
  * (4) B. Abi et al. [Fermilab Collaboration], Phys. Rev. Lett. 126, (2021) 141801.
  * (5) T. Aoyama, N. Asmussen, M. Benayoun, J. Bijnens, T. Blum, Phys. Rept. 887, (2020) 1-166.
  * (6) T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Phys. Rev. Lett. 124, 132002 (2020).
  * (7) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, P. Stoffer, JHEP 03, 101 (2020).
  * (8) Y.-Z. Fan, T.-P. Tang, Y.-L. S. Tsai, L. Wu, arXiv:2204.03693.
  * (9) C.-T. Lu, L. Wu, Y. Wu, B. Zhu, arXiv:2204.03796.
  * (10) P. Athron, A. Fowlie, C.-T. Lu, L. Wu, Y. Wu, B. Zhu, arXiv:2204.03996.
  * (11) G.-W. Yuan, L. Zu, L. Feng, Y.-F. Cai, arXiv:2204.04183.
  * (12) A. Strumia, arXiv:2204.04191.
  * (13) J. M. Yang, Y. Zhang, arXiv:2204.04202.
  * (14) J. d. Blas, M. Pierini, L. Reina, L. Silvestrini, arXiv:2204.04204.
  * (15) X. Du, Z. Li, F. Wang, Y. K. Zhang, arXiv:2204.04286.
  * (16) T.-P. Tang, M. Abdughani, L. Feng, Y.-L. S. Tsai, Y.-Z. Fan, arXiv:2204.04356.
  * (17) G. Cacciapaglia, F. Sannino, arXiv:2204.04514.
  * (18) M. Blennow, P. Coloma, E. Fernandez-Martinez, M. Gonzalez-Lopez, arXiv:2204.04559.
  * (19) K. Sakurai, F. Takahashi, W. Yin, arXiv:2204.04770.
  * (20) J. J. Fan, L. Li, T. Liu, K.-F. Lyu, arXiv:2204.04805.
  * (21) X. Liu, S.-Y. Guo, B. Zhu, Y. Li, arXiv:2204.04834.
  * (22) H. M. Lee, K. Yamashita, arXiv:2204.05024.
  * (23) Y. Cheng, X.-G. He, Z.-L. Huang, M.-W. Li, arXiv:2204.05031.
  * (24) Huayang Song, Wei Su, Mengchao Zhang, arxiV:2204.05085.
  * (25) E. Bagnaschi, J. Ellis, M. Madigan, K. Mimasu, V. Sanz, T. You, arXiv:2204.05260.
  * (26) A. Paul, M. Valli, arXiv:2204.05267.
  * (27) H. Bahl, J. Braathen, G. Weiglein, arXiv:2204.05269.
  * (28) P. Asadi, C. Cesarotti, K. Fraser, S. Homiller, A. Parikh, arXiv:2204.05283.
  * (29) L. D. Luzio, R. Grober, P. Paradisi, arXiv:2204.05284.
  * (30) P. Athron, M. Bach, D. H. J. Jacob, W. Kotlarski, D. Stockinger, A. Voigt, arXiv:2204.05285.
  * (31) J. Gu, Z. Liu, T. Ma, J. Shu, arXiv:2204.05296.
  * (32) J. J. Heckman, arXiv:2204.05302.
  * (33) K. S. Babu, S. Jana, P. K. Vishnu, arXiv:2204.05303.
  * (34) B.-Y. Zhu, S. Li, J.-G. Cheng, R.-L. Li, Y.-F. Liang, arXiv:2204.04688.
  * (35) Y. Heo, D.-W. Jung, J. S. Lee, arXiv:2204.05728.
  * (36) X. K. Du, Z. Li, F. Wang, Y. K. Zhang, arXiv:2204.05760.
  * (37) K. Cheung, W.-Y. Keung, P.-Y. Tseng, arXiv:2204.05942.
  * (38) A. Crivellin, M. Kirk, T. Kitahara, F. Mescia, arXiv:2204.05962.
  * (39) M. Endo, S. Mishima, arXiv:2204.05965.
  * (40) T. Biekotter, S. Heinemeyer, G. Weiglein, arXiv:2204.05975.
  * (41) R. Balkin, E. Madge, T. Menzo, G. Perez, Y. Soreq, J. Zupan, arXiv:2204.05992.
  * (42) L. Wang, J. M. Yang, Y. Zhang, arxiv:2203.07244.
  * (43) T. D. Lee, Phys. Rev. D 8, (1973) 1226.
  * (44) Y. Zhou, Y.-L. Wu, Eur. Phys. Jour. C 27, (2003) 577.
  * (45) K. Adikle Assamagan, A. Deandrea, P.-A. Delsart, Phys. Rev. D 67, (2003) 035001.
  * (46) S. Davidson, G. J. Grenier, Phys. Rev. D 81, (2010) 095016.
  * (47) Y. Omura, E. Senaha, K. Tobe, JHEP 1505, (2015) 028.
  * (48) R. Benbrik, C.-H. Chen, T. Nomura, Phys. Rev. D 93, (2016) 095004.
  * (49) Y. Omura, E. Senaha, K. Tobe, Phys. Rev. D 94, (2016) 055019.
  * (50) M. Lindner, M. Platscher, F. S. Queiroz, Phys. Rept. 731, (2018) 1-82.
  * (51) L. Wang, S. Yang, X.-F. Han, Nucl. Phys. B 919, (2017) 123-141.
  * (52) Y. Abe, T. Toma, K. Tsumura, JHEP 1906, (2019) 142.
  * (53) S. Iguro, Y. Omura, M. Takeuchi, JHEP 1911, (2019) 130.
  * (54) L. Wang, Y. Zhang, Phys. Rev. D 100, (2019) 095005.
  * (55) H.-X. Wang, L. Wang, Y. Zhang, Eur. Phys. Jour. C 81, (2021) 1007.
  * (56) S. Iguro, Y. Omura, M. Takeuchi, JHEP 09, (2020) 144.
  * (57) A. Crivellin, D. Müller, C. Wiegand, JHEP 06, (2019) 119.
  * (58) A. Das, T.i Nomura, H. Okada, S. Roy, Phys. Rev. D 96, (2017) 075001.
  * (59) P. S. B. Dev, R. N. Mohapatra, Y. Zhang, Phys. Rev. Lett. 120, 221804 (2018).
  * (60) T. Nomura, P. Sanyal, JHEP 05, (2021) 232.
  * (61) M. E. Peskin and T. Takeuchi, Phys. Rev. Lett. 65, 964 (1990).
  * (62) M. E. Peskin and T. Takeuchi, Phys. Rev. D 46, (1992) 381-409
  * (63) D. Eriksson, J. Rathsman, O. Stål, Comput. Phys. Commun. 181, (2010) 189.
  * (64) H.-J. He, N. Polonsky, S. Su, Phys. Rev. D 64, (2001) 053004.
  * (65) H. E. Haber, D. ONeil, Phys. Rev. D 83, (2011) 055017.
  * (66) A. Celis, V. Ilisie, A. Pich, JHEP 07, (2013) 053.
  * (67) Y. Amhis et al. [Heavy Flavor Averaging Group (HFAG) Collaboration], arXiv:1412.7515.
  * (68) T. Abe, R. Sato and K. Yagyu, JHEP 1507, (2015) 064.
  * (69) X.-F. Han, T. Li, L. Wang, Y. Zhang, Phys. Rev. D 99, (2019) 095034.
  * (70) S. Schael et al. [ALEPH and DELPHI and L3 and OPAL and SLD and LEP Electroweak Working Group and SLD Electroweak Group and SLD Heavy Flavour Group Collaborations], Phys. Rept. 427, (2006) 257.