# Floating-base manipulation on zero-perturbation manifolds

Brian A. Bittner, Jason Reid, Kevin C. Wolfe Johns Hopkins University Applied
Physics Laboratory, Laurel, MD, USA<EMAIL_ADDRESS>
###### Abstract

To achieve high-dexterity motion planning on floating-base systems, the base
dynamics induced by arm motions must be treated carefully. In general, it is a
significant challenge to establish a fixed-base frame during tasking due to
forces and torques on the base that arise directly from arm motions (e.g. arm
drag in low Reynolds environments and arm momentum in high Reynolds
environments). While thrusters can in theory be used to regulate the vehicle
pose, it is often insufficient to establish a stable pose for precise tasking,
whether the cause be due to underactuation, modeling inaccuracy, suboptimal
control parameters, or insufficient power. We propose a solution that asks the
thrusters to do less high bandwidth perturbation correction by planning arm
motions that induce zero perturbation on the base. We are able to cast our
motion planner as a nonholonomic rapidly-exploring random tree (RRT) by
representing the floating-base dynamics as pfaffian constraints on joint
velocity. These constraints guide the manipulators to move on zero-
perturbation manifolds (which inhabit a subspace of the tangent space of the
internal configuration space). To invoke this representation (termed a
perturbation map) we assume the body velocity (perturbation) of the base to be
a joint-defined linear mapping of joint velocity and describe situations where
this assumption is realistic (including underwater, aerial, and orbital
environments). The core insight of this work is that when perturbation of the
floating-base has affine structure with respect to joint velocity, it provides
the system a class of kinematic reduction that permits the use of sample-based
motion planners (specifically a nonholonomic RRT). We show that this allows
rapid, exploration-geared motion planning for high degree of freedom systems
in obstacle rich environments, even on floating-base systems with nontrivial
dynamics.

##  I Introduction

Figure 1: A thirteen-link Purcell swimmer with knowledge of its kinematic
chain via (4) (left) is unable to track (red proximal to end effectors) the
shapes (blue) sequentially attempted right to left due to perturbation induced
on the base (red proximal to base) during tracking. The thirteen-link Purcell
swimmer with knowledge of its dynamics as well as its kinematic chain written
as pfaffian constraints via (5) (right) is able to succeed given the same
task. This tracking result is explained in detail in Section IV. The
performant tracker leverages a surface in the joint velocity space which we
term a zero-perturbation manifold (ZPM). This surface is the core element of a
new planning framework covered in Section V. Accompanying video in
Supplementary Materials.

Floating-base manipulation platforms offer the potential to perform critical
tasking in underwater [1, 2, 3, 4], aerial [5, 6, 7], and orbital [8, 9]
environments. While manipulation has succeeded in high complexity environments
[10, 11, 12], floating-base systems have yet to achieve the same ability to
precisely manipulate objects in similar environments. The inability to
physically anchor any link of the kinematic chain can induce complex dynamics
on the ungrounded robot. These dynamics present differential constraints that
are often challenging for the motion planner to accommodate when pursuing
precision tasking. If the floating-base manipulator assumes a kinematic
operating environment, perturbations on the base can cause imprecision in
movement that is unpredictable and unsafe, especially in obstacle rich
environments. By considering these perturbation dynamics, one can pursue more
precise positioning of the end effector(s).

Typically in floating-base manipulation, the dynamics are represented as a
second-order mechanical system. Such a representation stipulates that given a
path in the joint configuration, the rate along that path has a variational
effect on the evolution of the robot in its workspace. These timing-based
concerns induce substantial complexity on the motion planning problem; these
planning problems are classified as kinodynamic. Generally, these kinodynamic
path planning problems are formulated as trajectory optimization problems.
These approaches typically approximate the Hamilton-Jacobi-Bellman equation,
which recursively formulates the optimal control law for a physical system.
These exploitative methods excel at finding local refinements to seeded
trajectories that leverage the agile capabilities of the system [13, 14, 15],
as long as the underlying model maintains high fidelity to the realized system
and environment. A separate trajectory optimization approach involves policy
gradient methods [16, 17], which are performant at rapidly finding collision
free paths, but again are typically used to find local refinements to
approaches (due to the exploitative nature of gradient descent).111While
stochastic formulations are designed to specifically avoid local minima, we
emphasize that the nature of a policy gradient search is far more localized
than sampling-based path planning methods. Additionally, this class of
solutions suffer from poor sampling scalability with the dimension of the
configuration – termed the curse of dimensionality. This notion of scalability
is encountered in many practical scenarios, such as adding joints to improve
the reachable configurations of a robot in key task domains (e.g. pick and
place in a cluttered environment). Exploitative trajectory optimization
methods thus tend to be insufficient for accomplishing many critical tasks due
to their locality and sample inefficiency. However, as is typically the case
in mobility architectures, a higher level planner might successfully hand off
a performant, coarse plan intended for final refinement by the trajectory
optimizer.

The main contribution of this work is the justification, design, and
implementation of a framework that decouples perturbations during floating-
base manipulation within a sampling-based motion planner (specifically an
RRT). We show that a kinematic reduction of the floating-base dynamics can be
applied to systems through mechanical design, environment selection, and
thruster-aided pose control. This reduces the motion planning problem from
kinodynamic to kinematic, making sampling-based motion planning algorithms
accessible. Empirically, such methods extend to high dimensional problems with
practical sample efficiency [18]. Sample efficiency for large dimensions is
critical, since as described earlier, higher numbers of joints will be
required for complex tasking in cluttered environments. Gaining this property
(sample efficiency for platforms with significant internal complexity –
required for many critical application domains) while respecting the dynamics
is the core capability that we pursue through this work. Using the method
detailed in Section V, we structurally expect (A) higher precision end
effector tracking and placement than planners with no dynamics knowledge as
well as (B) rapid computability of valid plans in many-joint, obstacle rich
scenarios. We show in obstacle rich scenarios that for emplacement and
extraction tasks, the proposed method can shift induced perturbation to a
trivial moment during the plan, allowing the ZPM to fascilitate a high
precision engagement with a target configuration.

The class of kinematic systems we assume stipulates that the body velocity is
a joint-defined linear mapping of joint velocity. These pfaffian constraints
can be used at any configuration to identify a subspace of joint velocities
where instantaneous joint motion yields zero perturbation (a sufficient
condition being that the dimension of the internal configuration is greater or
equal to the dimension of the position space). These pfaffian constraints are
termed a perturbation map, and this map provides a nonheuristic manifold that
can be used to guide the system to rapidly find collision-free plans. No
objective function designing or tuning is intended to be a part of
implementing the methods motivated in this work.

In Section II, we provide background on kinematic representations of floating-
base systems in application-oriented physical scenarios. We then define and
discuss key properties of the ZPM in Section III, followed by an example of an
end-effector tracking shapes without inducing motion on the base in Section
IV. In Section V we cover the proposed method for kinematic planning on
floating-base systems followed by Section 2 where we demonstrate global path
planning using the proposed method. Notably, the method outperforms (A)
kinematic algorithms that do not leverage dynamical information and (B)
trajectory optimizations that find locally optimal, low performance solutions.
Finally, Sections VII and VIII include discussion and conclusion of the
results.

##  II Background: kinematic reductions of Lagrangian systems

In this work, we seek to compute plans for floating-base systems when their
motion model can be written as a first-order (kinematic) system expressed as
pfaffian constraints connecting velocity on the shape space (joint space) $B$
to velocity on the position space $G$. These constraints are expressed as

$\dot{x}_{b}=P(\theta)\dot{\theta}$ (1)

yielding a linear mapping, $P$ (as a function of $n$ joints $\theta\in B$)
from joint velocity $\dot{\theta}\in TB$ to body velocity $\dot{x}_{b}\in TG$.
This structure will be shown to be an exact representation of the physics in
conventional orbital environments and scenarios governed by Rayleigh
dissipation. We’ll show that violations of this physical structure can be
related to intuitive phenomena and suggest that in many cases these complexity
inducing features can be designed out of the problem through mechanical
design, environment selection, and control.

###  II-A Obtaining Environmental Symmetry

Gravity directly violates environmental symmetry, inducing position and
orientation dependent forces on the floating-base. This will, for instance,
cause the robot to move when the joints are still, clearly violating our
assumed model at $\dot{\theta}=0$. Gravity is however negligible for orbital
robots as well as neutrally buoyant underwater systems222Some platforms such
as RE2s sapien class system are designed to be neutrally buoyant to avoid
dynamical complexity.. In less convenient scenarios, one might supply a model-
based control law to compensate for gravity with thrusters, through which one
could regain environmental symmetry of the floating-base under this thruster-
governed control law.

When we can assert symmetry of the dynamics over $G$ for Lagrangian systems
via a powerful reframing of the dynamics that provides three distinct
equations written for the generalized momentum, body velocity, and joint
acceleration of the system. This result from the field of geometric mechanics
is known as the reduced Lagrangian [19].

###  II-B Mitigating Interaction Between Inertia and Friction

When the generalized momentum starts at zero (as is the imagined case for many
floating-base tasks), the dynamics of the body are governed exclusively by
conservation of zero net angular and linear momentum. The resulting dynamics
take the form of (1) and is referred to as the mechanical connection [20]. A
similar phenomena occurs when the motion model becomes dominated by Rayleigh
dissipation. Here, friction-governed forces applied along the body yield a
viscous interaction with the fluid, again identical to the form in (1) and
termed the viscous connection [21].

Coupling of forces driven by inertia and friction result in second-order
behaviors. Take the free-floating astronaut as a mechanical connection. A
brief friction-driven interaction of brushing a satellite imparts an
uncancellable impulse on the astronaut. Post interaction it has obtained a net
momentum that it cannot cancel through internal motion, providing a drift that
violates the first order structure of the mechanical connection. Take a micro-
swimmer governed by a viscous connection. If we iteratively increase the
inertia of the swimmer, it will eventually noticeably drift as a function of
accumulated momentum while swimming. This drift again violates the first order
structure of the viscous connection.

Systems can be engineered to reduce interaction of these forces. Underwater
robots might move slowly to reduce excitation of momentum and leverage the
structure of a drag-dominated environment (assuming gravity compensation).
Aerial robots might be designed to avoid arm-related drag to encourage
isolation of momentum driven dynamics (assuming gravity compensation). Orbital
robots experience no environmental friction unless contacting objects in
space.

###  II-C Approximating Nonviscous Friction

In general, friction can certainly act in non-viscous ways, violating the
structure of (1). One can think of the viability of the Rayleigh dissipation
model as the viability of a first-order Taylor series approximation of the
more general friction-governed motion model
$\dot{x}_{b}|\theta=f(\dot{\theta})$. 333In many cases this can be a complex,
non-smooth function (e.g. in the case of coulomb friction) motivating the use
of stochastic linearization to obtain useful model parameters. Isotropy of the
friction can also have a significant impact on behavior [22, 23], further
motivating the use of stochastic linearizations to build linearizations that
support prediction over a desired range of velocities. However, it has been
observed in practice that the assumption of a Rayleigh dissipation model
provides a strong fit across a surprisingly broad class of robots [24, 25].

In this section we have detailed the motion model of interest and the relevant
assumptions for dynamics and control to arrive at its structure. We now
outline a method which allows us to use sampling-based motion planning while
respecting a broad class of nontrivial floating-base dynamics.

##  III Zero Perturbation Manifolds

We now shift focus to explaining how the representations motivated in Section
II populate a manifold that can be exploited for motion planning that induces
zero perturbation on the floating-base. The linear mapping $P(\theta)$ has a
null space computable at all shapes $\theta$

$null(P(\theta))\subset TB$ (2)

which inhabits a subspace of the shape velocity space $TB$. This null space
provides a basis in the joint velocity space for exploration from a shape that
will infinitesimally invoke zero perturbation on the base. This null space
changes as a function of $\theta$, and so during a plan, the allowable motions
available to the joints are subject to change. When the null space changes
continuously with respect to shape, the null space populates a connected
submanifold of joint velocity space. We know this since the notion of
connectedness is preserved under continuous mappings. However, the null space
is provably not continuous in this sense at singular $P$ where the rank of the
nullspace increases over that of non-singular $P$. Thus, the zero perturbation
manifold (ZPM), $null(P)$, has a basis of size $dim(B)-dim(G)$ and
continuously changes over its connected regions. It is however a larger
manifold that is not continuous, e.g. at locations of singular $P$.

Continuous or not, the null space indicates at all internal configurations
instantaneous directions to move the joints that induce zero base
perturbation. Continuity of the null space is however important for generating
$C^{1}$ smooth arm motions. If a null space basis vector disappears while
moving a joint in that direction, it will appear as a nonsmooth path in joint
space (if persisting along the zero perturbation manifold). This manifold is
not unlike leveraging the null space of components of the manipulator Jacobian
to stabilize elevation, roll and pitch when moving a weight along a table, as
is done in work that in-part inspired this approach [26].

In many classical manipulation scenarios, additional joints provide
exploitable redundancy. Likewise, additional joints directly increase the
dimension of the zero perturbation manifold. We make a final note about the
challenge of connecting two specific points in the shape space. The
probability of the null space directly aligning with two randomly selected
shapes is approximately zero:

$\displaystyle\forall\delta\in
B:\textbf{p}(null(P(\theta)))^{T}null(P(\theta+\delta)=1)\approx 0.$ (3)

where p provides our single-case use of notation for probability. While it is
easy to grow a tree along the null space from a known joint location, it can
be hard to grow the tree along the null space toward a known joint location.
One might think of an Ackermann steering vehicle (a system whose maneuvers are
governed by pfaffian constraints) that is placed arbitrarily close to a final
location, then asked to arrive at that location directly, without knowledge of
Lie brackets or ”parallel parking” maneuvers. In many scenarios, there is not
a direct route to that location. In higher-degree of freedom floating base
manipulation problems, it is almost never the case that such as direct route
exists. To our knowledge, it is an open problem whether structures exist to
procedurally connect such proximal shape configurations along the zero
perturbation manifold. Lie brackets offer an interesting avenue for further
inspection, but are not the subject of this work.

##  IV Local End Effector Tracking along the Zero Perturbation Manifold

Here we leverage the exploitable redundancy property discussed in Section III
to show that a many jointed, bimanual swimming robot can complete a welding
task with a closed form solution. This solution will implicitly leverage the
structure of the zero perturbation manifold, which grows in dimension with
each joint made available to the system.

We compare two platforms. A baseline platform has access to its kinematic
chain, so induced disturbances on the body cannot be accounted for:

$\dot{x}_{ee}=J(\theta)\dot{\theta}.$ (4)

It takes a pseudoinverse of the Jacobian to compute $\dot{\theta}$ given a
desired velocity for the end effector $\dot{x}_{ee}$. The second platform has
access to both its kinematic chain and the model of the viscous fluidic model
that it inhabits:

$\begin{bmatrix}\dot{x}_{ee}\\\
\dot{x}_{b}\end{bmatrix}=\begin{bmatrix}J(\theta)\\\
P(\theta)\end{bmatrix}\dot{\theta}.$ (5)

where $\dot{x}_{b}$ can be set to zero.

This toy system is a simple adaptation of the Purcell swimming model [27, 28],
where we magnify the ratio of across-to-along body forces by a factor of 5,
mimicking a more fin-like relationship with the water to excite larger
perturbations. In the trial, $\dot{x}_{ee}$ is set by the residual between the
location of the hand (with respect to the base) and the desired hand location
(with respect to the base). The desired hand location is a function of time
that traces each shape (from right to left)). The closest hand is driven to
track the desired hand location, the other being free to use its neighboring
links as a tail-like appendage if necessary to stabilize the base. Fig. 1
shows that the thirteen-link swimmer was capable of exactly tracing the
desired shapes with a closed form solution when leveraging the pfaffian form
of the dynamics. Meanwhile, the platform embedded with only kinematic-chain
knowledge both drifted and was unable to draw shapes resembling those desired.
These results were not intended to capture phenomena such as collision-
avoidance including self-collisions, but these topics are addressed in the
next section.

##  V RRT on Zero Perturbation Manifolds

In Sections II and III we overviewed a key class of floating-based dynamics
and its relevance for distinction of a zero-perturbation manifold that
occupies $TB$. Here, we will show how such a manifold can be used effectively
within sampling-based motion planners for rapid exploration based planning on
high degree of freedom floating-base robots in obstacle rich settings.
Equation (3) shows us that connecting two points in the shape space can be
challenging. This suggests that objects like probabilistic roadmaps (PRMs)
will be structurally challenging to generate, since the inter-connectedness of
nodes requires direct confrontation of this issue. Turning our attention to
single-query planning approaches, RRT offers the ability to grow quickly from
a known location. As we have established, knowing the final location does not
make it easy to find a plan that arrives exactly at the final location. One
might think of a car placed at a close distance but incorrect orientation to a
desired pose. While Lie brackets or primitive-based motion planning may offer
a solution for the car, it is not obvious that such approaches are
computationally tractable in the floating-base manipulation setting.

We observe three options, which include growing the RRT from an initial
location to a desired location, growing a bi-RRT from both locations towards
each other, and growing from the desired location toward the initial location.
We accept the third option as suboptimal but highly practical. An RRT will be
designed to rapidly grow from a location of interest toward the robots initial
configuration, all the while restricting itself to the zero perturbation
manifold. The optimal path selected from this tree is the core contribution of
this work. The robot will start by connecting itself from initial
configuration to the ZPM. Once it reaches the trajectory along the ZPM, we
assert that it will stop while the thrusters mitigate any induced
perturbation. Now, after regulation, the floating-base system is afforded a
zero perturbation path, through clutter, to the object of interest.

Algorithm 1 $ZPMRRT(\theta_{s},\theta_{g})$

1:T.init($\theta_{g}$)

2:while $TimeRemaining()$ do

3: $\theta_{sample}\leftarrow sample(i,\theta_{s})$

4: $\theta_{nearest}\leftarrow nearest(T,\theta_{sample})$

5: $\theta_{extend}\leftarrow
ExtendAlongZPM(T,\theta_{nearest},\theta_{sample})$

6: if $||\theta_{extend}-\theta_{s}||<\alpha$ then

7: return $Smooth(ExtractPath(T,\theta_{s}))$

8: end if

9: return $null$

10:end while

Algorithm 2 $ExtendAlongZPM(T,\theta_{near},\theta_{sample})$

1:$\theta_{extend}\leftarrow\theta_{near};\texttt{ }spacing\leftarrow 0$

2:while $Timeout()$ or  $Collision(\theta_{extend})$ do

3: $\delta\theta\leftarrow\theta_{extend}-\theta_{sample}$

4: $u,s,v\leftarrow svd(P(\theta_{extend}))$

5: $P_{null}\leftarrow v[dim(G):]$

6: if $P_{null}\perp\delta\theta$ then

7: return $\theta_{extend}$

8: else

9: $\theta_{extend}\leftarrow\theta_{extend}+unit(\delta\theta
P_{null}^{T}P_{null})dt$

10: $spacing\leftarrow spacing+dt$

11: if $spacing>L$ then

12: $T.addVertex(\theta_{near},\theta_{extend})$

13: $\theta_{near}\leftarrow\theta_{extend};\texttt{ }spacing\leftarrow 0$

14: end if

15: end if

16:end while

17:return $\theta_{extend}$

Figure 2: Using a classical RRT (top), the thirteen-link Purcell swimmer is
not able to accurately obtain a desired pose during emplacement (taken from 1
of 30 trials detailed in Fig. 3). We show snapshots of the swimmer at the
beginning (first column), midpoint (second column), and end (final column),
approaching hand locations (pairs of green x’s) while avoiding obstacles (grey
circles). Using ZPMRRT (bottom), the swimmer is able to reach the desired
endpoint. We show snapshots of the behavior, with columns carrying the same
timestamps with the exception of the middle, which marks the transition point
to the ZPM. Accompanying video in Supplementary Materials.  Figure 3: This
Fig. show the displacement of the Purcell swimmer in distance (left $\Delta d$
in bodylengths) and orientation (right $\Delta\psi$ in radians) over 30 trials
for emplacement with one obstacle (top row, obstacle numbers identified in
Fig. 2), two obstacles (middle row), and three obstacles (bottom row). The
mean (solid) and bounds (transparent band) of performance are plotted for the
RRTZPM (blue, Algorithm 1) motivated in this work as well as the classical RRT
(red) which does not incorporate any dynamics knowledge.

Here we will outline our approach for floating-base manipulation along zero
perturbation manifolds. Algorithm 1 is provided a start and goal internal
configuration for the robots ($\theta_{s},\theta_{g}$). As described above, we
will start at the goal location, where there is a higher need for precision
(such as in emplacement and extraction tasking). We initialize a tree $T$ with
this node and query a sample function at each iteration. This will provide a
random internal configuration or the start configuration (with probability
$0.5$ in this work). We then attempt to grow the tree toward this sampled
configuration, while adhering to the zero perturbation manifold. The
computation details of this process are available in Algorithm 2, which
involves computation of the null space via singular value decomposition (SVD)
and projection of the error between goal and current point along the computed
subspace of the internal configuration that carries zero eigenvalues. If we
avoid singular configurations (as is assumed in this algorithm), then the null
space is consistently of $rank(null(P))=dim(B)-dim(G)$ and can be indexed
directly as the corresponding bottom rows of V. We can project the growth
direction onto this manifold (see line 9 of Algorithm 2 where $unit$ produces
a unit vector), and stop extending in the event of a timeout, collision, or
perpendicular orientation of the null space and remaining error. During
extension, we grow the tree by adding a node after such a length has been
extended from a prior node in the tree. Once the search has timed out or
achieved sufficient proximity to the goal configuration, the algorithm returns
a smoothed path. In this work, smoothing involves attempting to find a
shortcut between the start and end point of the path by direct projection onto
the ZPM, then recursively repeating the procedure on two randomly split
segments of the path (until those segments reduce to two points). This
smoothing process is repeated three times in this work.

##  VI Results

Here we will ask the thirteen-link Purcell swimmer to plan in constrained
environments while executing an emplacement task. This task will require that
both hands reach a particular location in the workspace. Through expectation
that the base does not move, we stipulate a shape that will achieve the
desired end effector locations and provide this desired shape as a goal to the
planner. As mentioned before, the planner will grow from the goal towards the
starting location. We identify that there are a variety of viable starting
configurations from which to approach the emplacement task, and encode this
within the planner. When trying to connect to starting shapes $\theta_{s}$, we
will provide a randomly selected, generalized $U$-shape for the swimmer. 444A
$U$ includes a negative quarter rotation of a joint left of the base and a
positive quarter rotation of a joint right of the base. We can provide this
freedom of choice since we will assume that the system can start from a
further displaced position, change configuration to any other $U$-shape, and
then trivially re-engage at the position identified in the planner555Any such
shape can translate in the coordinate perpendicular to the base link to exit
or enter the starting configuration, making it a viable entry point to start
the emplacement plan from.. Thus, we claim to be solving a broader planning
problem (starting from a disengaged, distal position with any $U$-shaped
initial joint configuration) even though we only explicitly solve the non-
trivial component of the planning problem (once engaged in an obstacle rich
environment).

We start with one obstacle and provide the RRTZPM 30 separate trials to find
connections to $U$-shapes along the zero perturbation manifolds. We then ask
the RRT to connect the same two configurations. In each trial we allow the
planner to replan 4 times after 10s intervals if a collision free path is not
found. Paths are termed collision free if they would not cause a collision
with a static base. In every trial both planners succeed during the first 10s
interval. Fig. 3 shows the displacement as a function of a fraction of path
length, spanning the start and finish of the attempt. It is clear that the
proposed method accrues a brief initial displacement to access the zero
peturbation manifold (distinguishing the effect of the properties following
from (3). However, once on the ZPM, the platform is capable of achieving the
desired configuration with zero displacement. Meanwhile, the classical RRT
accrues error, critically during the final steps of engagement of the desired
internal configuration. This classical approach persistently accrues error as
it approaches the final location, exactly when the proposed method accrues
zero error by structure.

We see this relationship persist as a second and third obstacle are added.
Notably, during trials with 3 obstacles, the RRT unsuccessfully attempts to
find a collision free path in 20 of 30 trials. Across all numbers of
obstacles, the proposed method (for each obstacle) experiences an initial
perturbation of the base followed by an approximately zero perturbation
approach to the final configuration. We observe this difference in performance
in Fig. 2, where we can see the drift over the duration of the trial for the
classical RRT on the top row, as well as the lack of drift for the example on
the bottom row.

Figure 4: Using iterative TVLQR (top), the thirteen-link Purcell swimmer is
not able to accurately obtain a desired pose during emplacement (taken from 1
of 30 trials in Fig. 5). We show snapshots of the swimmer at the beginning
(first column), midpoint (second column), and end (final column), approaching
hand locations (pairs of green x’s) while avoiding obstacles (grey circles).
Using RRTZPM (bottom), the swimmer is able to reach the desired endpoint. We
show snapshots of the behavior, with columns carrying the same timestamps with
the exception of the middle, which marks the transition point to the ZPM.
Accompanying video in Supplementary Materials.  Figure 5: This Fig. shows the
displacement of the purcell swimmer in distance (left $\Delta d$ in
bodylengths) and orientation (right $\Delta\psi$ in radians) over 30 trials
with no obstacles. The mean (solid) and bounds (transparent band) of
performance are plotted for RRTZPM (blue, Algorithm 1) motivated in this work
as well as the iterative TVLQR which is subject to the locality of
exploitation-driven solutions.

In obstacle-free environments we compare this approach to an iterative time
varying linear quadratic regulator (TVLQR). This provides the opportunity to
observe the exploration-driven solutions available via RRTZPM with respect to
this exploitation-driven approximation of the Hamilton-Jacobi-Bellman
equation. The iterative TVLQR is seeded with a direct route from initial to
final shape for the linear system:

$\begin{bmatrix}\dot{x}_{b}\\\
\dot{\theta}\end{bmatrix}=\begin{bmatrix}P(\theta)\\\
I\end{bmatrix}\dot{\theta}.$ (6)

where we take the shape velocity $\dot{\theta}$ as the control input. We
discretize each trajectory at 50 samples. For $n=dim(B)$ and $m=dim(G)$ take

$\displaystyle Q=\begin{bmatrix}I^{mxm}&&0\\\ 0&&0\end{bmatrix}$
$\displaystyle Q_{f}=\begin{bmatrix}I^{mxm}&&0\\\ 0&&I^{nxn}\end{bmatrix}$
$\displaystyle R=Ie^{-2}$ (7)

where $Q_{f}$ is the final element of $Q$. By iteratively supplying the output
of a TVLQR solution to itself, we provide the opportunity to further refine
the solution with the context of the dynamics centered at the most recent
solution. At the end of 20 iterations we have a trajectory that has been
locally refined to find the desired joint location while avoiding
perturbations. We ask the swimmer to achieve the previously identified
emplacement shape configuration from randomly selected $U$-shapes (this time
the swimmer cannot select a $U$-shape.) Fig. 4 shows that over 30 trials we
observe the RRTZPM is capable of accessing the ZPM, although it accessed later
in the maneuver than in Fig 3 when it had the choice of starting $U-shape$
locations. The iterative TVLQR critically accrues error as it arrives at the
emplacement location, repeating the precision mitigating behavior seen in the
RRT in Fig. 3. We can see impact of this error at the emplacement location in
Fig. 5, where the error accrued through iterative TVLQR renders the platform
unable to complete the emplacment task.

##  VII Discussion

The active use of the other arm to decouple the generated drag force in Fig. 1
reminded us of the use of tails in many animals. This may provide inspiration
for design of systems in similar dynamic environments. Overall, the end
effector tracking results of Fig. 1 are not entirely surprising (prior work
exists in null space projection for swimming robots [4]). The example
reinforces what is known in manipulation, that redundancy offers a greater
variety of options to achieve the desired result. Observing this property in
the dynamic case while tracking shapes informs us that this property, since it
exists for a local tracking problem, may be possible to embed in a more global
planning problem with more constraints.

The emplacement performance of the RRTZPM method in Fig. 3 and Fig. 2 shows us
that it was possible to embed this property in a high degree of freedom
platform with three obstacles. It is clear from comparison to the RRT planner
(that neglects dynamics) that the ZPM is providing significant value. The
practical relevance of this result is made clear in Fig. 2, where the RRT
platform can be seen colliding and drifting from the desired location. Such
precision is unacceptable for even coarse emplacement tasking.

This relationship persists when we relate the motivated method to iterative
TVLQR in an obstacle-free environment. The exploitation driven algorithm can
locally refine the initial trajectory, but is unable to find a more global
solution that permits acceptable precision for the emplacement task. This
supports the perceived value of our method which can respect dynamics while
finding inventive solutions rapidly.

The measured perturbation rejection on these simulated systems can be impacted
largely by the time-step of the integration. In principle, simulating at a
small enough time step would permit a numerically zero perturbation when
integrating along the ZPM. For real systems, feasibility of precisely tracking
shape evolution along the ZPM may significantly impact the quality of the
overall plan.

##  VIII Conclusion

The results suggest that the proposed methods provide a novel class of
algorithms that can address the weaknesses of classical sampling-based motion
planners and optimization-based approaches. We showed that we were able to
maintain the ”sampling-scalability at high dimension” feature of RRTs while
providing them the dynamical context to avoid precision-related failures. We
showed that the proposed method maintained the exploration-geared
inventiveness of RRTs, permitting viable solutions where exploitation-driven
optimization schemes will fail.

Our examples are restricted to a planar swimming simulated system. Further
questions will remain pertinent in extensions of these results to other
systems, simulated or real. How viable are pfaffian constraint approximations
for systems that do not exist exactly in high or low Reynolds domains? How
many degrees of freedom are typically needed to achieve similar results in
three-dimensional scenarios? How easy is it to execute trajectories along this
manifold with conventional hardware? In future work we plan to explore open
questions directly on hardware in underwater and orbital domains.

## References

  * [1] James C Kinsey, Ryan M Eustice and Louis L Whitcomb  “A survey of underwater vehicle navigation: Recent advances and new challenges”  In _IFAC conference of manoeuvering and control of marine craft_ 88, 2006, pp. 1–12  Lisbon
  * [2] Enrico Simetti et al.  “Floating underwater manipulation: Developed control methodology and experimental validation within the trident project”  In _Journal of Field Robotics_ 31.3 Wiley Online Library, 2014, pp. 364–385
  * [3] Patryk Cieslak, Pere Ridao and Mariusz Giergiel  “Autonomous underwater panel operation by GIRONA500 UVMS: A practical approach to autonomous underwater manipulation”  In _2015 IEEE International conference on robotics and automation (ICRA)_ , 2015, pp. 529–536  IEEE
  * [4] Gerald Brantner and Oussama Khatib  “Controlling Ocean One: Human–robot collaboration for deep-sea manipulation”  In _Journal of Field Robotics_ 38.1 Wiley Online Library, 2021, pp. 28–51
  * [5] Fabio Ruggiero, Vincenzo Lippiello and Anibal Ollero  “Aerial manipulation: A literature review”  In _IEEE Robotics and Automation Letters_ 3.3 IEEE, 2018, pp. 1957–1964
  * [6] Marin Kobilarov  “Nonlinear trajectory control of multi-body aerial manipulators”  In _Journal of Intelligent & Robotic Systems_ 73.1 Springer, 2014, pp. 679–692
  * [7] Moju Zhao, Kei Okada and Masayuki Inaba  “Versatile articulated aerial robot DRAGON: Aerial manipulation and grasping by vectorable thrust control”  In _The International Journal of Robotics Research_ SAGE Publications Sage UK: London, England, 2022, pp. 02783649221112446
  * [8] Peter Kazanzides et al.  “Teleoperation and Visualization Interfaces for Remote Intervention in Space”  In _Frontiers in Robotics and AI_ 8 Frontiers Media SA, 2021
  * [9] Evangelos Papadopoulos, Farhad Aghili, Ou Ma and Roberto Lampariello  “Robotic manipulation and capture in space: A survey”  In _Frontiers in Robotics and AI_ Frontiers, 2021, pp. 228
  * [10] Eric Krotkov et al.  “The DARPA robotics challenge finals: Results and perspectives”  In _The DARPA Robotics Challenge Finals: Humanoid Robots To The Rescue_ Springer, 2018, pp. 1–26
  * [11] Jeffrey Mahler et al.  “Learning ambidextrous robot grasping policies”  In _Science Robotics_ 4.26 American Association for the Advancement of Science, 2019, pp. eaau4984
  * [12] Dale McConachie, Andrew Dobson, Mengyao Ruan and Dmitry Berenson  “Manipulating deformable objects by interleaving prediction, planning, and control”  In _The International Journal of Robotics Research_ 39.8 SAGE Publications Sage UK: London, England, 2020, pp. 957–982
  * [13] Sébastien Kleff et al.  “High-frequency nonlinear model predictive control of a manipulator”  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 7330–7336  IEEE
  * [14] Tomasz Rybus, Karol Seweryn and Jurek Z Sasiadek  “Control system for free-floating space manipulator based on nonlinear model predictive control (NMPC)”  In _Journal of Intelligent & Robotic Systems_ 85.3 Springer, 2017, pp. 491–509
  * [15] Dongjae Lee, Hoseong Seo, Dabin Kim and H Jin Kim  “Aerial manipulation using model predictive control for opening a hinged door”  In _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 1237–1242  IEEE
  * [16] Nathan Ratliff, Matt Zucker, J Andrew Bagnell and Siddhartha Srinivasa  “CHOMP: Gradient optimization techniques for efficient motion planning”  In _2009 IEEE International Conference on Robotics and Automation_ , 2009, pp. 489–494  IEEE
  * [17] Mrinal Kalakrishnan et al.  “STOMP: Stochastic trajectory optimization for motion planning”  In _2011 IEEE international conference on robotics and automation_ , 2011, pp. 4569–4574  IEEE
  * [18] Mohamed Elbanhawi and Milan Simic  “Sampling-based robot motion planning: A review”  In _Ieee access_ 2 IEEE, 2014, pp. 56–77
  * [19] Jim Ostrowski and Joel Burdick  “The geometric mechanics of undulatory robotic locomotion”  In _The international journal of robotics research_ 17.7 Sage Publications Sage CA: Thousand Oaks, CA, 1998, pp. 683–701
  * [20] Jerrold E Marsden and Jim Ostrowski  “Symmetries in motion: Geometric foundations of motion control”  National Academies Press, 1998
  * [21] Scott D Kelly and Richard M Murray  “Geometric phases and robotic locomotion”  In _Journal of Robotic Systems_ 12.6 Wiley Online Library, 1995, pp. 417–431
  * [22] Daolin Ma and Alberto Rodriguez  “Friction variability in planar pushing data: Anisotropic friction and data-collection bias”  In _IEEE Robotics and Automation Letters_ 3.4 IEEE, 2018, pp. 3232–3239
  * [23] Silas Alben  “Efficient sliding locomotion with isotropic friction”  In _Physical Review E_ 99.6 APS, 2019, pp. 062402
  * [24] Henry C Astley et al.  “Surprising simplicities and syntheses in limbless self-propulsion in sand”  In _Journal of Experimental Biology_ 223.5 The Company of Biologists Ltd, 2020, pp. jeb103564
  * [25] Dan Zhao et al.  “Walking is like slithering: A unifying, data-driven view of locomotion”  In _Proceedings of the National Academy of Sciences_ 119.37 National Acad Sciences, 2022, pp. e2113222119
  * [26] Dmitry Berenson, Siddhartha S Srinivasa, Dave Ferguson and James J Kuffner  “Manipulation planning on constraint manifolds”  In _2009 IEEE international conference on robotics and automation_ , 2009, pp. 625–632  IEEE
  * [27] Leif E Becker, Stephan A Koehler and Howard A Stone  “On self-propulsion of micro-machines at low Reynolds number: Purcell’s three-link swimmer”  In _Journal of fluid mechanics_ 490 Cambridge University Press, 2003, pp. 15–35
  * [28] Ross L Hatton and Howie Choset  “Geometric swimming at low and high Reynolds numbers”  In _IEEE Transactions on Robotics_ 29.3 IEEE, 2013, pp. 615–624


# Precision Muon-Related Observables as a Tool
to Constrain New Physics Models

Gabriel Lourenço${}^{1\,}$111E-mail<EMAIL_ADDRESS>, Andre
Milagre${}^{2\,}$222E-mail<EMAIL_ADDRESS>, Rui
Santos${}^{1,3\,}$333E-mail<EMAIL_ADDRESS>, João P.
Silva${}^{2\,}$444E-mail<EMAIL_ADDRESS>,
1Centro de Física Teórica e Computacional, Faculdade de Ciências,
Universidade de Lisboa, Campo Grande, Edifício C8 1749-016 Lisboa, Portugal
2 CFTP, Departamento de Física,
Instituto Superior Técnico, Universidade de Lisboa,
Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal
3ISEL - Instituto Superior de Engenharia de Lisboa,
Instituto Politécnico de Lisboa 1959-007 Lisboa, Portugal

(August 28, 2024)

###### Abstract

We propose a set of precision muon-related observables that serve as a tool to
constrain new physics models. Using LEP’s precision measurements on the
$Z$-boson pole, we derive bounds on the new physics quantum contributions to
the decay $Z\to\mu^{+}\mu^{-}$. We show that the new precision observables
have a real impact on two specific models that solve the $g-2$ anomaly and
provide a sound dark matter candidate.

##  1 Introduction

Until 2022 there was an observed anomaly in the semileptonic $B$ meson decay
rates that suggested a violation of lepton flavor universality. However, the
recent reinterpretation performed by the LHCb collaborations washed out the
discrepancy between the SM prediction in the $b\to s\mu^{+}\mu^{-}$ transition
[1, 2] and the experimentally measured value. The anomaly was present for a
long time in the measurements of the ratios of the exclusive branching
fractions, $R(K^{(*)})={\cal B}(B\to K^{(*)}\mu^{+}\mu^{-})/{\cal B}(B\to
K^{(*)}e^{+}e^{-})$, obtained by the LHCb Collaboration [3, 4, 5] and by the
Belle Collaboration [6, 7] although with larger errors. In an attempt to solve
this problem, several new models were proposed among which some focused on
adding a number of new fields that would solve the problems via one-loop
contributions [8, 9, 10, 11, 12, 13, 14, 15, 16, 17] from particles living in
a new Dark Sector. Some of these models [11, 16, 17] also solved the muon
$g-2$ anomaly [18, 19, 20, 21, 22] while providing a Dark Matter (DM)
candidate compatible with all present constraints. It is important to note
that although the anomaly in the semileptonic $B$ meson decay rates vanished,
the latter models are still in agreement with all experimental data and are
able to explain the $g-2$ anomaly and still provide a scalar dark matter
candidate.

Triggered by exploring in more detail the constraints on this class of models
we thought it was timely to use the precision measurement of muon-related
observables to constrain any new physics contributions coming from the
proposed extensions of the Standard Model (SM). The idea is inspired by the
use of precision observables related to $Z\to b\bar{b}$ [23, 24, 25, 26] to
constrain models with an extended scalar sector of the SM. In fact, the
precision measurements performed at LEP on the $Z$-boson pole can be used for
other $Z$ decay channels such as $Z\to\mu^{+}\mu^{-}$ and
$Z\to\tau^{+}\tau^{-}$. If a new Dark Sector contributes to a new observable
like for instance the muon $g-2$, it will certainly appear in other
observables related to the muon. The need to have large couplings to muons
will then manifest itself in the $Z$ decays to muon allowing to restrict the
parameter space of the model. We will further show that the results for the
tau lepton are obtained directly from the muon and that because the
experimental measurements related to the tau have larger error bars, the
results are less restrictive in constraining the parameter space of a new
model.

The paper is organized as follows. In Sec. 2 we present and discuss the
$Z\to\mu^{+}\mu^{-}$ observables to be used and in Sec. 3 we present the
calculation of the radiative corrections. In Sec. 4 we present the results and
apply the constraints to two models of the new physics. Our conclusions are
given in Sec. 5.

##  2 Observables

In the SM, the interaction Lagrangian between the $Z$ and the muons can be
written as

$\mathcal{L}_{Z\mu^{+}\mu^{-}}=-\frac{g}{c_{W}}Z^{\lambda}\bar{\mu}\,\gamma_{\lambda}\left(g_{\mu
L}P_{L}+g_{\mu R}P_{R}\right)\mu,$ (2.1)

where $P_{L}$ and $P_{R}$ are the left and right-handed chirality projectors.
At tree-level, the $Z\mu^{+}\mu^{-}$ couplings in the SM are

$\displaystyle g_{\mu L}$ $\displaystyle=-\frac{1}{2}+s_{W}^{2},$
$\displaystyle g_{\mu R}$ $\displaystyle=s_{W}^{2},$ (2.2)

where $s_{W}$ and $c_{W}$ are the sine and cosine of the Weinberg angle,
respectively. These quantities were measured at LEP [27] and the most precise
values for the corresponding SM couplings read [27]

$\displaystyle g_{\mu L}^{\textsc{SM}}$ $\displaystyle=-0.26929,$
$\displaystyle g_{\mu R}^{\textsc{SM}}$ $\displaystyle=0.23208.$ (2.3)

As previously discussed, quantum loop corrections make the couplings deviate
from their tree-level values. These corrections are partly generated by SM
particles and partly from the new physics contributions.

The muon couplings can be derived from two $Z$-pole observables related to the
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ process: the asymmetry parameter
($\mathcal{A}_{\mu}^{0}$) and the muon hadronic decay rate ($R_{\mu}$).
Henceforth, we assume all observables to be corrected for initial-state
radiation but make the distinction between observables that include effects of
final-state radiation by placing a ”0” superscript on the ones that do not.

The asymmetry parameter can be calculated from the muon forward-backward
asymmetry $A_{\textup{FB}}^{0}$. This is a measure of the angular asymmetry in
the distribution between final state particles produced in the forward
direction (defined by the angle they make with the initial electron beam, i.e.
$\theta<\frac{\pi}{2}$) and the ones that are produced backward
($\theta>\frac{\pi}{2}$). Its formal definition reads

$A_{\textup{FB}}^{0}=\frac{\sigma_{F}^{0}-\sigma_{B}^{0}}{\sigma_{F}^{0}-\sigma_{B}^{0}}=\frac{3}{4}\mathcal{A}_{e}^{0}\mathcal{A}_{\mu}^{0},$
(2.4)

with,

$\displaystyle\sigma_{F}^{0}$
$\displaystyle=2\pi\int_{0}^{1}d\cos\theta\,\frac{d\sigma^{0}}{d\,\cos\theta},$
$\displaystyle\sigma_{B}^{0}$
$\displaystyle=2\pi\int_{-1}^{0}d\cos\theta\,\frac{d\sigma^{0}}{d\,\cos\theta}.$
(2.5)

which are the cross sections of the muons detected in the forward and backward
directions, respectively. Since the center of mass energy range is near the
$Z$ resonance ($s\sim m_{Z}^{2}$), $\mathcal{A}_{\mu}^{0}$ can be written in
terms of the tree-level left and right $Z\mu^{+}\mu^{-}$ couplings as

$\mathcal{A}_{\mu}^{0}\equiv\frac{\sqrt{1-4\mu_{\mu}}(g_{\mu L}^{2}-g_{\mu
R}^{2})}{(1-\mu_{f})(g_{\mu L}^{2}+g_{\mu R}^{2})+6\mu_{\mu}\,g_{\mu
L}\,g_{\mu R}}\ \ \ ,\ \ \ \mu_{\mu}\equiv\frac{m_{\mu}^{2}}{m_{Z}^{2}}.$
(2.6)

where $m_{\mu}$ is the mass of the muon and $m_{Z}$ that of the $Z$ boson.

The muon hadronic decay rate is defined as the ratio between the partial decay
width of the $Z$ to hadrons and to muons,

$\displaystyle R_{\mu}^{0}$ $\displaystyle\equiv\frac{\Gamma^{0}(Z\to
had)}{\Gamma^{0}(Z\to\mu^{+}\mu^{-})},$ $\displaystyle\Gamma^{0}(Z\to had)$
$\displaystyle=\sum_{q=u,d,c,s,b}\Gamma^{0}(Z\to q\bar{q}),$ (2.7)

which depends on the couplings between the quarks and the $Z$. As will be
explained in Section 3, although cross-sections are defined experimentally
without final-state corrections, partial widths are not [27]. To avoid
performing QCD-related computations for the final state radiation of quarks,
instead of working with $R_{\mu}^{0}$ directly we define the quantity 555This
observable was used to constrain the parameter space of a Leptophilic two-
Higgs doublet model [28]. $R_{e\mu}=R_{e}/R_{\mu}$ which, at tree level, and
without final state radiation factors, can be written as

$R_{e\mu}^{0}\equiv\frac{R_{e}^{0}}{R_{\mu}^{0}}=\frac{\Gamma^{0}(Z\to\mu^{+}\mu^{-})}{\Gamma^{0}(Z\to
e^{+}e^{-})}=\frac{g_{\mu L}^{2}+g_{\mu R}^{2}-\mu_{\mu}(g_{\mu L}^{2}+g_{\mu
R}^{2}-6\,g_{\mu L}\,g_{\mu R})}{g_{eL}^{2}+g_{eR}^{2}}.$ (2.8)

There is good agreement between the predicted and experimental results of both
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ since the central SM value is in the
uncertainty interval of the corresponding experimental result. Table 1
presents the corresponding SM prediction and fitted values for $R_{e}$,
$R_{\mu}$, $\mathcal{A}^{0}_{\mu}$, and $R_{e\mu}$. The uncertainty in
$R_{e\mu}$ was calculated from the quadratic sum of the uncertainties in
$R_{e}$ and $R_{\mu}$, since these observables are poorly correlated [27].

Observable | Exp | SM
---|---|---
$R_{e}$ | $20.804\pm 0.050$ | $20.736\pm 0.010$
$R_{\mu}$ | $20.784\pm 0.034$ | $20.736\pm 0.010$
$\mathcal{A}^{0}_{\mu}$ | $0.142\pm 0.015$ | $0.1468\pm 0.0003$
$R_{e\mu}$ | $1.0009\pm 0.0029$ | $1.0000\pm 0.0007$
Table 1: Experimental (Exp) value and Standard Model (SM) predictions for
$R_{e}$, $R_{\mu}$, $\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ [29]. The SM
predictions include the one-loop contributions, as well as the dominant two-
loop contributions. $R_{e\mu}$ is obtained from the hadronic decay rates of
the electron and the muon in the first two lines.

The SM predictions include one-loop contributions and the dominant two-loop
contributions as described in the review Electroweak Model and Constraints on
New Physics [29].

##  3 Radiative corrections

When calculated at the $Z$ pole, the radiative corrections to
$Z\to\mu^{+}\mu^{-}$ are linked to the radiative corrections to the process
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ at one loop level. This allows for the use of
the measurements performed at LEP of the latter to constrain the higher-order
contributions of new physics processes to the former. In all calculations of
one-loop processes, we will use the on-shell renormalization procedure
developed in [30, 31]. This scheme conveniently uses parameters with clear
physical meaning that can be measured experimentally and also allows for the
separation of QED, and Weak corrections into two distinct gauge invariant
classes. In this classification, QED corrections correspond to diagrams with a
virtual photon exchanged in the loop and diagrams with real photon emission in
the form of bremsstrahlung radiation originating from the initial or final
state particles. The one-loop QED-type diagrams contributing to the processes
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ and $Z\to\mu^{+}\mu^{-}$ are represented in
Figures 1 and 2, respectively. Summing over these diagrams eliminates IR
divergences and the result is UV finite by itself after renormalization [32].
The remaining one-loop diagrams correspond to the Weak corrections and include
the corrections to the vertices, corrections to the photon and $Z$
propagators, and box diagrams with massive boson exchange. These are
represented for each process in Figs. 3 and 4. The QED corrections depend on
the experimental setup, whereas the Weak corrections depend only on the
underlying structure of the theory. In the case of colored final states,
additional QCD corrections have to be included. These are formally identical
to the QED corrections [32].

Figure 1: QED-type diagrams contributing to the $e^{+}e^{-}\to\mu^{+}\mu^{-}$
process.

Figure 2: QED-type diagrams contributing to the $Z\to\mu^{+}\mu^{-}$ process.

Figure 3: Weak-type diagrams contributing to the $e^{+}e^{-}\to\mu^{+}\mu^{-}$
process.

Figure 4: Weak-type diagrams contributing to the $Z\to\mu^{+}\mu^{-}$ process.

As previously stated, partial widths can only be compared with experiment once
the effects of final-state radiation are included. In the case of the $Z$
decaying into leptons ($l$), the final state QED corrections are accounted for
by introducing a multiplicative factor as

$\Gamma(Z\to l^{+}l^{-})=\Gamma^{0}(Z\to l^{+}l^{-})F_{\textup{QED}}^{l},$
(3.1)

where,
$F_{\textup{QED}}^{l}=1+\frac{3}{4}Q_{l}^{2}\frac{\alpha(m_{Z}^{2})}{\pi}+\mathcal{O}\left(\alpha^{2}\right)$
[32]. At leading order in $\alpha$ this factor is independent of the lepton
species. As a consequence, these corrections effectively cancel out in the
measured $R_{e\mu}$ value which, in turn, allows one to identify
$R_{e\mu}\approx R_{e\mu}^{0}$.

The weak corrections are expected to change the chiral structure of the theory
at one loop level. However, the matrix elements that do not preserve the tree-
level structure can be neglected [31]. The propagator corrections at one loop
level come from the mixing of the photon and $Z$ propagators and from the $Z$
and photon exchange diagrams. These enter as redefinitions of the QED and
neutral current coupling strengths, $Z$ width, and sine of the Weinberg angle.
The vertex corrections are defined as the $Z\mu^{+}\mu^{-}$ three-point
functions at one loop level and can be interpreted as additive contributions
to the left and right-handed coupling constants. In the on-shell
renormalization scheme, the counterterms needed to renormalize the vertex
depend on the bare self-energy of the mixed $\gamma-Z$ propagator and on the
external legs’ self-energy [33].

###  3.1 New Physics Contributions

In order to be able to constrain the parameters stemming from the new physics
Lagrangian, it is important to divide the different pieces of the one-loop
corrections into the two and three-point function related to the
$Z\mu^{+}\mu^{-}$ vertex.

One-loop corrections to the bosonic propagator are expected in the presence of
new physics. These can be absorbed into the definition of the effective
Weinberg angle and are parameterized by the Peskin-Takeuchi $S$, $T$, and $U$
oblique parameters [34]. The deviation from the SM prediction can be written
in the form [35]

$\delta
s_{W}^{2}=s_{W}^{\text{eff}\,2}-s_{W}^{\text{SM}\,2}=\frac{\alpha}{c_{W}^{2}-s_{W}^{2}}\left(\frac{1}{4}S-s_{W}^{2}c_{W}^{2}T\right),$
(3.2)

where $s_{W}^{\text{eff}}$ is the sine of the effective Weinberg angle, which
accounts for all known and unknown corrections, and $s_{W}^{\text{SM}}$ is the
sine of the SM effective Weinberg angle. In the following discussion, we use
the experimental value $s_{W}^{2}=[s_{W}^{\text{eff 2}}]_{\textup{EXP}}\approx
0.23148$ [29]. The most recent experimental data shows the oblique parameters
agreeing with the SM value of zero [29]. In both models studied in Section 4,
we follow the procedure outlined in [17] to constrain our parameter space
accordingly and neglect new physics’ direct contributions to the propagator.
We do, however, take into account how new physics affects the $T$ parameter,
which is the most relevant to constrain our parameter space, but do so in the
scan prior to computing the corrections to the muon-related observables.
Additional contributions coming from box diagrams would have to be included,
but neither model we chose to work with has new fields coupling to the first-
generation fermions. It is important to note that whatever the model under
study is, since we can write the corrections as a function of parameters that
are very constrained, the oblique parameters, these corrections will always be
negligible.

New physics corrections to the vertex are absorbed into a redefinition of the
left and right $Z\mu^{+}\mu^{-}$ couplings. At one-loop, the renormalized
couplings can be spit as
$\hat{g}_{L,R}=g_{L,R}^{\textup{SM}}+\delta\widetilde{g}_{L,R}$, where the
first term contains the tree-level and all SM one-loop contributions, and
$\delta\widetilde{g}_{L,R}$ all one loop vertex contributions from new
physics. From now on we drop the muon lower index since there will be no
source of confusion about which lepton we are referring to. The new physics
term is defined as [33]

$\delta\widetilde{g}_{L,R}=-\frac{c_{W}s_{W}}{e}\widetilde{\Gamma}^{Z\mu\mu}_{L,R}+g_{L,R}\delta\widetilde{Z}_{L,R}-g_{L,R}\frac{c_{W}}{s_{W}}\frac{\widetilde{\Sigma}^{\gamma
Z}(0)}{m_{Z}^{2}}+s_{W}c_{W}\frac{\widetilde{\Sigma}^{\gamma
Z}(0)}{m_{Z}^{2}}.$ (3.3)

where a $+\delta s_{W}^{2}$ factor would have to be added to the effective
coupling if the oblique parameters were not neglected [24]. The first term in
(3.3) is computed from the bare three-point function containing the new
particles of the model. The decomposition into left and right components is
done as666Where we neglect terms proportional to the anomalous weak-electric
and magnetic dipole moments (i.e. $\propto\sigma_{\mu\nu}k^{\nu}$, with
$k^{\nu}$ the momentum of the $Z$ boson) because of a $m_{\mu}/m_{Z}$
suppression when compared to $\widetilde{\Gamma}_{L,R}^{Z\mu\mu}$ [36].

$i\widetilde{\Gamma}_{\lambda}^{Z\mu\mu}=\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure22.pdf}}}}}\approx\ \
-\frac{ie}{s_{W}c_{W}}\gamma_{\lambda}\left(\widetilde{\Gamma}_{L}^{Z\mu\mu}P_{L}+\widetilde{\Gamma}_{R}^{Z\mu\mu}P_{R}\right).$
(3.4)

Furthermore, $\widetilde{\Sigma}^{\gamma Z}$ can be calculated from the sum of
diagrams containing new particles that contribute to the mixed $\gamma-Z$ bare
self-energy777The longitudinal term proportional to $k^{\alpha}k^{\beta}$ can
be neglected because it is suppressed by a $m_{f}^{2}/m_{Z}^{2}$ factor [33].

$-i\widetilde{\Sigma}^{\gamma
Z}(k^{2})g^{\alpha\beta}=\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure23.pdf}}}}}$ (3.5)

calculated at zero momentum transfer $k^{2}=0$. Finally,
$\delta\widetilde{Z}_{L,R}$ is the counter-term for the left or right-handed
component of the muon wave function,

$\delta\widetilde{Z}_{L,R}=-\widetilde{\Sigma}_{L,R}(m_{\mu}^{2})-m_{\mu}^{2}\left[\widetilde{\Sigma}^{\prime}_{L}(m_{\mu}^{2})+\widetilde{\Sigma}^{\prime}_{R}(m_{\mu}^{2})+2\widetilde{\Sigma}^{\prime}_{S}(m_{\mu}^{2})\right],$
(3.6)

where the functions $\widetilde{\Sigma}_{L,R,S}$ come from the decomposition
of the muon bare self-energy into

$i\widetilde{\Sigma}_{\mu}(\not{p})=\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure24.pdf}}}}}=i\not{p}\widetilde{\Sigma}_{L}(p^{2})P_{L}+i\not{p}\widetilde{\Sigma}_{R}(p^{2})P_{R}+im_{\mu}\widetilde{\Sigma}_{S}(p^{2}).$
(3.7)

with $\Sigma^{\prime}=\frac{\partial\Sigma}{\partial s}$.

###  3.2 Connection with experiment

As shown, the most significant new physics contributions come from the
$Z\mu^{+}\mu^{-}$ vertex corrections. To a good approximation,
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ can be defined at one-loop order by
replacing the tree-level couplings in (2.6) and (2.8) by
$\hat{g}_{L,R}=g_{L,R}^{\textup{SM}}+\delta\widetilde{g}_{L,R}$. In the
presence of new physics, we expect the theoretical prediction of the
observables $\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ to deviate from the SM.
This deviation can be written as

$\mathcal{A}_{\mu}^{0}=\mathcal{A}_{\mu}^{\textsc{SM}}+\delta\widetilde{\mathcal{A}}_{\mu},$
(3.8a) $R_{e\mu}=R_{e\mu}^{\textsc{SM}}+\delta\widetilde{R}_{e\mu}.$ (3.8b)

The relation between the new physics contributions to the couplings and the
deviations in the observables can be derived by taking
$|\delta\widetilde{g}_{L,R}|\ll|g_{L,R}^{\textup{SM}}|$, and expanding
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ to first order in
$\delta\widetilde{g}_{L,R}$. In the small mass limit ($\mu_{\mu}\to 0$) we
obtain

$\delta\widetilde{\mathcal{A}}_{\mu}\approx-3.6343\delta\widetilde{g}_{L}-4.2154\delta\widetilde{g}_{R},$
(3.9a)
$\delta\widetilde{R}_{e\mu}\approx-4.2619\delta\widetilde{g}_{L}+3.6744\delta\widetilde{g}_{R},$
(3.9b)

Using this formalism, the parameter space of models with new couplings to the
muon can be systematically constrained, by requiring the deviation from the SM
predictions to fall within the experimental uncertainties presented in Table
1.

##  4 Results

In order to see how these observables can be used to test and constrain
extensions of the SM, we consider the models presented in [16, 17]. We perform
a multi-parameter scan to identify the allowed parameter space obeying the
already implemented flavor constraints and the newly defined conditions
arising from the muon-related observables. The flavor-related constraints are
the new $R\left(K^{\left(*\right)}\right)$ now in agreement with the SM
predictions (as discussed in the introduction),
$\mathcal{B}\left(B_{S}\to\mu^{+}\mu^{-}\right)$, $B_{S}-\bar{B}_{S}$ mixing
and $b\to s\gamma$. We force both models to be able to explain the $g-2$
anomaly which gives a special role to the muon couplings. Besides that, we
also force the DM candidate to saturate the relic density measurement, take
into account the bounds from direct and indirect detection, and consider the
constraints imposed by collider searches (see [16, 17] for a detailed
discussion).

The models add two scalar fields, $\Phi_{2}$ and $\Phi_{3}$ (the second with a
color charge), and a vector-like fermion $\chi$. These fields are collectively
called the Dark Sector. A new $Z_{2}$ symmetry is introduced, under which the
SM fields are even and the Dark Sector fields are odd. In the one we call
Model 3, the new fermion is a singlet of $SU(2)_{L}$ and the scalars are
doublets. In Model 5, the fermion is an $SU(2)_{L}$ doublet and the scalars
are singlets. In both models, the neutral component of the $\Phi_{2}$ scalar
is split into its real and imaginary degrees of freedom
($\frac{1}{\sqrt{2}}(S+iA)$). The flavor phenomenology is the same for both
$S$ and $A$, hence, without loss of generality, we assume $m_{S}<m_{A}$,
effectively choosing $S$ to be the DM candidate. The $SU(3)_{c}$ and
$SU(2)_{L}$ representations and weak hypercharges of the Dark Sector particles
of both models are summed up in Table 2.

Model 5
---
| $SU(3)_{c}$ | $SU(2)_{L}$ | $U(1)_{Y}$
$\chi_{R}$ | 1 | 2 | $-1/2$
$\Phi_{2}$ | 1 | 1 | $0$
$\Phi_{3}$ | 3 | 1 | $+2/3$

Model 3
---
| $SU(3)_{c}$ | $SU(2)_{L}$ | $U(1)_{Y}$
$\chi_{R}$ | 1 | 1 | $-1$
$\Phi_{2}$ | 1 | 2 | $+1/2$
$\Phi_{3}$ | 3 | 2 | $+7/6$

Table 2: Representations and weak hypercharges of Dark Sector of Model 5
(left) and Model 3 (right). Adapted from [16] and [17].

LEP searches for unstable heavy vector-like charged leptons set a lower bound
on the mass of $\chi^{\pm}$ of 101.2 GeV [37]. Other more recent constraints
from the LHC do not apply, since these assume either that the vector-like
lepton couples to the tau [38] or a very small amount of missing transverse
energy in the final states [39]. Other limits on the masses come from the
flavor constraints and are detailed in [16, 17].

In the following two sections, we present the new physics contributions to
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$, and the results of the scans for the
two models. We performed the random scan and checked which conditions were
satisfied for each point. The relic density was computed with micrOMEGAs [40],
which solves the Boltzmann equations numerically. We show projections from the
parameter space with a similar color scheme to that used in [16, 17], where
all cyan points satisfy the mass constraints and explain the $B$ meson data
(including the latest reinterpretation of the LHCb collaboration) in the
2$\sigma$ range. In addition to the previous constraints, the blue points
agree with DM relic abundance observations in the 2$\sigma$ range. Similarly,
the green points verify all previous constraints, together with the
restrictions coming from DM direct detection and collider searches. Finally,
the red points satisfy all previous constraints and agree with the muon $g-2$
experiments within $3\sigma$. We updated experimental values used in these
conditions in light of recent results. In particular, the more accurate
analysis done by the LHCb Collaboration of the exclusive branching fractions,
$R(K)$ and $R(K^{*})$, [2] has implications in the $C_{9}^{NP}$ Wilson
coefficient used in [16]. We constrained all parameter space points to be
within the $2\sigma$ range around the new best-fitted value
$C_{9}^{NP}=-0.19\pm 0.06$ [41]. Furthermore, bounds on Dark Matter relic
density are from the 2018 Planck Collaboration [42], and we have included the
most recent upper bounds on the spin-independent DM-nucleon scattering cross-
section from the LZ experiment [43].

For every point that verified all previous conditions, we computed
$\delta\widetilde{\mathcal{A}}_{\mu}$ and $\delta\widetilde{R}_{e\mu}$ and
checked whether
$\mathcal{A}_{\mu}^{\textup{SM}}+\delta\widetilde{\mathcal{A}}_{\mu}$ and
$R_{e\mu}^{\textup{SM}}+\delta\widetilde{R}_{e\mu}$ were kept within their
$1\sigma$ experimental confidence levels presented in Table 1, that is

$-0.0203<\delta\widetilde{\mathcal{A}}_{\mu}<0.0097,$ (4.1a)
$-0.0020<\delta\widetilde{R}_{e\mu}<0.0038.$ (4.1b)

In the discussion that follows, parameter space points that verified this
condition were assigned the color yellow.

As noted in [17], the scalar fields in Model 3 introduce a nonzero
contribution to the oblique parameter $T$, which is given by [44]

$T=\frac{g^{2}}{64\pi^{2}m_{W}^{2}\alpha}\left[F(m_{\phi_{l}}^{2},m_{S}^{2})+F(m_{\phi_{l}}^{2},m_{A}^{2})-F(m_{S}^{2},m_{A}^{2})\right],$
(4.2)

where $g$ is the $SU(2)_{L}$ coupling constant, $m_{W}$ is the mass of the $W$
boson and $\alpha$ the fine-structure constant. The function $F(A,B)$ is
defined as

$F(A,B)=\begin{cases}\frac{A+B}{2}-\frac{AB}{A-B}\log\frac{A}{B},\text{ if
$A\neq B$.}\\\ 0\hskip 74.2617pt,\text{ if $A=B$.}\end{cases}$ (4.3)

The vector-like fermion and the colored scalars have a vanishing contribution
to the $T$ parameter. The first contributes through a diagram similar to the
photon self-energy in QED, which at zero momentum transfer vanishes, and the
colored scalars produce a term proportional to
$F(m_{\phi_{q}^{5/3}}^{2},m_{\phi_{q}^{2/3}}^{2})$, which is zero by (4.3)
since we consider $m_{\phi_{q}^{5/3}}=m_{\phi_{q}^{2/3}}$. To account for the
effect of new physics in the oblique parameter $T$, we impose a supplementary
condition that every point in the Model 3 scan must obey, which is to be
within the 2$\sigma$ experimental uncertainty interval. At present, the best
experimental value for this quantity is $T=0.03\pm 0.12$ [29]. In the case of
Model 5, the contribution to the $T$ parameter is zero [16] and it provides no
constraint on the model’s parameters.

The models were implemented both in LanHEP [45] and FeynRules [46]. We used
the packages FeynArts [47] to generate the relevant diagrams and FeynCalc [48,
49] to compute the corresponding amplitudes and perform the integration. The
results are written in terms of the Passarino-Veltman integral functions [50].
We used the package LoopTools [51] to numerically evaluate the integrals.

###  4.1 Model 5

In Model 5, the single diagram contributing to the self-energy of the mixed
$\gamma-Z$ propagator, $\Sigma^{\gamma Z}$, has the fermion $\chi$ running in
the loop. When evaluated at zero momentum transfer, its amplitude was computed
and confirmed to equal zero. As a consequence, the task of computing the NP
corrections to the $Z\mu^{+}\mu^{-}$ vertex was reduced down to the
computation of $\widetilde{\Gamma}_{\lambda}^{Z\mu\mu}$ and
$\widetilde{\Sigma}_{\mu}(\not{p})$. The relevant new diagrams contributing at
the one-loop level are

$\begin{split}i\widetilde{\Gamma}_{\lambda}^{Z\mu\mu}&=\vbox{\hbox{
\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure25.pdf}}}}}+\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure26.pdf}}}}},\\\
i\widetilde{\Sigma}_{\mu}(\not{p})&=\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure27.pdf}}}}}+\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure28.pdf}}}}}.\end{split}$ (4.4)

Besides the SM parameters, the functions
$\widetilde{\Gamma}_{\lambda}^{Z\mu\mu}$ and
$\widetilde{\Sigma}_{\mu}(\not{p})$ depend only on the new masses $m_{S}$,
$m_{A}$, $m_{\chi}$ and on the absolute value of the new Yukawa muon coupling
$|y_{\mu}|$. We confirmed that the divergences were indeed canceled in (3.3).
In addition, we found both functions to be symmetric under
$m_{S}\leftrightarrow m_{A}$ and to scale as $\sim|y_{\mu}|^{2}$, as expected.
In Model 5 the new particles couple to left-handed muons but not to their
right-handed counterparts. As a consequence, in the limit $\mu_{\mu}\to 0$,
Model 5 was found to only impact the left $Z\mu^{+}\mu^{-}$ coupling, at the
one-loop level. This behavior is evidenced by the blue line in Fig. 5, where
the result of a random scan to the parameter space forms a

Figure 5: Representation of the $\mathcal{A}_{\mu}^{0}-R_{e\mu}$ space. The
$1\sigma$ and $2\sigma$ experimental regions are shown in different shades of
orange, together with a red star marking the experimental central value. The
red diamond marks the SM prediction, and the blue points show how Model 5
impacts this prediction.

Figure 6: For $m_{\chi}=150$ GeV and for different values of $y_{\mu}$, the
left panel shows the $m_{S}-m_{A}$ projection plane with the points that keep
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ within the $1\sigma$ uncertainty
ranges. Darker shades of orange correspond to higher values for $y_{\mu}$. The
middle and right panels correspond to the parameter space slices with
$m_{\chi}=400$ GeV and $m_{\chi}=600$ GeV, respectively. The grey-shaded
region is forbidden because $m_{S}$ was chosen as the DM candidate and for
scanning purposes we have chosen $m_{S}<m_{\chi}+10$ GeV and
$m_{S}<m_{A}+10\textup{ GeV}$ with no loss of generality.

line in the $R_{e\mu}-\mathcal{A}_{\mu}^{0}$ plane. In this plot the $1\sigma$
and $2\sigma$ experimental uncertainty regions are represented in different
shades of orange, the red star marks the experimental central value, and the
red diamond marks the present SM prediction. The blue points show the
departing from the SM prediction when $m_{\chi}$ was varied from 101.2 to 1000
GeV, $m_{S}$ from 5 GeV to 1000 GeV, $m_{A}$ from 15 GeV to 1000 GeV, and
$y_{\mu}$ from 0 to $4\pi$ [16]. It is possible to check that the new
particles always give positive contributions to both $\mathcal{A}_{\mu}^{0}$
and $R_{e\mu}$ meaning that, in the domain of interest,
$\delta\widetilde{g}_{L}\leq 0$ and $\delta\widetilde{g}_{R}=0$. As a
consequence of the slope of the line formed by the blue points, every point
satisfying $R_{e\mu}$ at some level of confidence was found to automatically
satisfy $\mathcal{A}^{0}_{\mu}$ at that same level of confidence. The fact
that there are points outside the experimental uncertainty indicates that the
model can be constrained using these observables.

A better understanding of the observables’ impact on Model 5 is obtained by
studying the dependence of $\delta\widetilde{g}_{L}$ on the NP parameters. In
Fig. 6 we plot the parameter space regions that keep both
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ within their $1\sigma$ experimental
uncertainty range, as a function of $m_{S}$ and $m_{A}$. In the left (middle,
right) panel, $m_{\chi}$ was fixed and set equal to $150$ ($400$, $600$) GeV.
Each overlapping region has a fixed $y_{\mu}$, with darker shades of orange
corresponding to higher coupling values. The forbidden region satisfying888The
(arbitrary) choice of $+10$ GeV was originally implemented in the numerical
scan to guarantee the DM candidate ($S$) to be the lightest of the Dark
Sector’s particles. $m_{S}>m_{\chi}+10$ GeV and $m_{S}>m_{A}+10\textup{ GeV}$
is gray-shaded. We find that higher muon coupling strengths make for more
stringent constraints on the parameter space, specifically in the region of
low $m_{S}$ and low $m_{A}$. Fixing $y_{\mu}$, we conclude that these
constraints get less rigid the higher $m_{\chi}$ is. Besides the
$m_{S}\leftrightarrow m_{A}$ symmetry being manifest, we also note that in the
decoupling limit, i.e. big masses and/or small coupling, the NP corrections
behave as expected.

(a)

(b)

(c)

(d)

Figure 7: Allowed parameter space of Model 5 projected in the planes
$m_{S}-m_{A}$ (Top Left), $m_{S}-m_{\chi}$ (Top Right), $m_{S}-y_{\mu}$
(Bottom Left), and $m_{\chi}-y_{\mu}$ (Bottom Right). The color scheme
replicates the one used in [16]. Only the yellow points agree with
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ at the $1\sigma$ level while satisfying
all previous constraints.

In Fig. 7, we show the parameter space projections relevant to the
$Z\mu^{+}\mu^{-}$ observables. The color yellow was assigned to the red data
points (the ones satisfying all previous conditions) which agree with both
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ at the $1\sigma$ level. We have found
all parameter space points to deviate from $\mathcal{A}_{\mu}^{\textup{SM}}$
by less than $1\sigma$, meaning all constraints arise solely due to
$R_{e\mu}$. In these type of models, the $g-2$ discrepancy can be solved by a
large value of the coupling $y_{\mu}$. The interplay between the high
$y_{\mu}$ values required by $g-2$ and the low $y_{\mu}$ values required by
the new observables is apparent in these plots. The $m_{\chi}-y_{\mu}$
projection of Fig. 7(d) is the most interesting because of the positive
correlation between these parameters. In addition, it reveals a small strip
excluded by $R_{e\mu}$ that was accepted by previous constraints. However,
this region vanishes when we extend the threshold of the $R_{e\mu}$
requirement from $1\sigma$ to the $2\sigma$ level.

Our results (yellow points) are consistent with the previous lower bound on
$m_{A}$ of $m_{A}\gtrsim 50$ GeV, but slightly disagree with the bounds on
$m_{S}$ ($m_{S}\in[30,350]$ GeV), $m_{\chi}$ ($m_{\chi}\lesssim 600$ GeV), and
$y_{\mu}$ ($y_{\mu}\in[1.4,\sqrt{4\pi}]$) found in [16]. The initial analysis
lacked an observable that could impose an upper bound on $y_{\mu}$, and only
allowed it to range up to $\sqrt{4\pi}$. The introduction of
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ let us push $y_{\mu}$ up to its
perturbative limit of $4\pi$. The new region with $y_{\mu}>\sqrt{4\pi}$ now
allows for $m_{\chi}$ to range from $101.2$ to $1000$ GeV unboundedly, as
shown in Fig. 7(d). Due to the positive correlation between $m_{\chi}$ and
$m_{S}$ (Fig. 7(b)), the upper bound on the latter is extended to $\sim 410$
GeV, while keeping its lower bound. In sum, we argue that the bounds on the
parameters are more accurately described by

$\displaystyle m_{A}$ $\displaystyle\gtrsim 50\textup{ GeV},$ $\displaystyle
y_{\mu}$ $\displaystyle\in[1,4\pi],$ $\displaystyle m_{\chi}$
$\displaystyle>101.2\textup{ GeV},$ $\displaystyle m_{S}$
$\displaystyle\in[30,410]\textup{ GeV}.$

###  4.2 Model 3

We computed the weak corrections to the left and right-handed coupling
constants for Model 3. The $\gamma-Z$ mixing diagrams with the fermion $\chi$
and the scalars in the loop cancel in the zero momentum transfer limit. As
such, only the vertex and muon self-energy diagrams of (4.5) contribute to the
couplings. The fact that in this model the scalars are doublets of $SU(2)$
means that they couple to the $Z$ boson, which leads to two diagrams that were
forbidden in Model 5

$\begin{split}i\widetilde{\Gamma}_{\lambda}^{Z\mu\mu}&=\vbox{\hbox{
\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure29.pdf}}}}}+\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure30.pdf}}}}}+\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure31.pdf}}}}}+\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure32.pdf}}}}},\\\
i\widetilde{\Sigma}_{\mu}(\not{p})&=\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure33.pdf}}}}}+\vbox{\hbox{\leavevmode\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\immediate\hbox{{\leavevmode\includegraphics{}{}{external/zmumu-
figure34.pdf}}}}}.\end{split}$ (4.5)

Figure 8: Representation of the $\mathcal{A}_{\mu}^{0}-R_{e\mu}$ plane. The
$1\sigma$ and $2\sigma$ experimental regions are shown in different shades of
orange, together with a red star marking the experimental central value. The
red diamond marks the SM prediction. The blue points show how Model 3 impacts
the SM prediction.

In Fig. 8 we show in blue how the Model 3 points impact the SM prediction
marked by the red diamond in the $\mathcal{A}_{\mu}^{0}-R_{e\mu}$ plane. The
$1\sigma$ and $2\sigma$ experimental regions are shown in different shades of
orange, together with a red star marking the experimental central value. The
relevant parameters to calculate the corrections to the couplings are the
masses of the fermion, $m_{\chi}$, and of the scalars, $m_{S}$, $m_{A}$, and
the coupling $y_{\mu}$. For Model 3, $m_{S}$ was varied from 5 to 100 GeV and
$m_{A}$ from 100 to 1000 GeV. The limits for $m_{\chi}$ are the same as in
Model 5. $\delta\widetilde{g}_{R}\propto m_{\mu}$ and is always several orders
of magnitude smaller than $\delta\widetilde{g}_{L}$, which is the dominant
term. By varying the parameters in the allowed ranges defined in [17], we can
see that, unlike in Model 5, $\delta\widetilde{g}_{L}$ and
$\delta\widetilde{g}_{R}$ can be both positive or negative. Once again it is
possible to use $\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ to constrain the
parameter space.

Figure 9: $m_{S}-m_{A}$ projection plane with the values that keep
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ within the $1\sigma$ uncertainty ranges
for $m_{\chi}=200,350$ and $700$ GeV, respectively. Darker shades of orange
correspond to higher values for $y_{\mu}$ as indicated in the legend. Higher
values of $y_{\mu}$ bring $m_{S}$ and $m_{A}$ closer together. Smaller values
of $m_{\chi}$ decrease the allowed parameter space. The grey-shaded region is
forbidden.

In Fig. 9 we plot the 1$\sigma$ acceptance regions for $\mathcal{A}_{\mu}^{0}$
and $R_{e\mu}$ in the $m_{S}-m_{A}$ plane for different values of $y_{\mu}$
and $m_{\chi}$. The forbidden region satisfying $m_{S}>m_{\chi}+10$ GeV,
$m_{S}>m_{A}+10$ GeV and $m_{S}+m_{A}<m_{Z}$ is gray-shaded. From these plots,
we can conclude that in general the stronger the coupling, the closer the
masses of the scalars must be. Also the smaller the values of $m_{\chi}$, the
smaller the accepted region is. These are important results that show that on
one hand the smaller gap between the mass of two scalars is preferred,
especially for smaller values of $m_{\chi}$ and large values of the muon
coupling. If further restrictions are imposed on even just one of the scalar
masses limiting it to a small range, then it is expected that this would
translate to a severe restriction of the allowed parameter space.

(a)

(b)

(c)

(d)

Figure 10: Allowed parameter space of Model 3 projected in the planes
$m_{S}-m_{A}$ (top left), $m_{A}-m_{\chi}$ (top right), $m_{A}-y_{\mu}$
(bottom left) and $m_{\chi}-y_{\mu}$ (bottom right). The color scheme
replicates the one used in [16]. Only the yellow points agree with
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ at the $1\sigma$ level while satisfying
all previous constraints. Every point agrees with $T$ to 2$\sigma$. Unlike in
Model 5, the new muon-related restrictions can effectively limit the allowed
parameter space. Figure 11: Allowed parameter space of Model 3 projected in
the plane $m_{A}-m_{\phi_{l}}$. The color scheme replicates the one used in
[16]. Only the yellow points agree with $\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$
at the $1\sigma$ level while satisfying all previous constraints. Every point
agrees with $T$ to 2$\sigma$.

Finally, in Figs.10 and 11 we review our results for Model 3 in view of the
new constraints. We use the scan with the new results as already discussed for
Model 5 and show four projections of the parameter space from the random scan.
Once again the yellow points correspond to the points that verify all dark
matter and flavor-related constraints and the new muon-related constraints to
1$\sigma$. Every point shown in the plots is in accordance with the
experimental result for the $T$ parameter to 2$\sigma$. This restriction has a
visible consequence in the $m_{S}-m_{A}$ plane, where points where the
difference between $m_{A}$ and $m_{S}$ is too large to obey this condition are
excluded. In stark contrast to what we got for Model 5, in Model 3
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ can significantly limit the allowed
regions. As before, $R_{e\mu}$ provides the strongest constraint, with all
points accepted by the condition defined with $\mathcal{A}_{\mu}$ to 1
$\sigma$. As noted in [17], in Model 3 the dark matter constraints limit
$m_{S}$ to be bigger than 42 GeV and smaller than 76 GeV. This combined with
the muon-related constraints means the allowed values for $m_{A}$ are much
more restricted than indicated by the $g-2$ constraints. The dark matter
constraints end up indirectly limiting the allowed values of $m_{A}$. This can
be seen in the $m_{S}-m_{A}$, $m_{A}-m_{\chi}$ and $m_{A}-y_{\mu}$
projections, where a large portion of the red points are excluded by
$R_{e\mu}$. The yellow points are concentrated as expected in the smaller mass
regions. In the plane $m_{A}-y_{\mu}$ we can see that as $y_{\mu}$ increases,
the allowed range of $m_{A}$ gets closer to the allowed interval of values for
$m_{S}$. This is consistent with the behavior observed in Fig. 9. As in Model
5, the muon-related constraints now reject a small strip of previously
accepted space in the $m_{\chi}-y_{\mu}$ plane. These observables also provide
a much stronger constraint to the mass of $\phi_{l}$ than the $g-2$ results,
with $m_{\phi_{l}}$ now allowed to be at most 265 GeV. We can summarize the
new bounds on the parameters of Model 3 as

$\displaystyle m_{A}$ $\displaystyle\in[100,194]\textup{ GeV},$ $\displaystyle
y_{\mu}$ $\displaystyle\in[1,4\pi],$ $\displaystyle m_{\chi}$
$\displaystyle>101.2\textup{ GeV},$ $\displaystyle m_{S}$
$\displaystyle\in[50,75]\textup{ GeV}.$ $\displaystyle m_{\phi_{l}}$
$\displaystyle\in[70,265]\textup{ GeV},$

Figure 12: $m_{S}-m_{A}$ projection plane with the values that keep
$\mathcal{A}_{\mu}^{0}$ and $R_{e\mu}$ (left panel), and
$\mathcal{A}_{\tau}^{0}$ and $R_{e\tau}$ (right panel) within the $2\sigma$
uncertainty ranges for $m_{\chi}=200$ GeV. Darker shades of orange correspond
to higher values for $y_{\mu}$ and $y_{\tau}$ as indicated in the legend. The
grey-shaded region is forbidden.

We end this section with a comment on the use of the process
$Z\to\tau^{+}\tau^{-}$. In Fig. 12 we present the $m_{S}-m_{A}$ projection
plane with the regions that keep $\mathcal{A}_{\mu}^{0}$, and $R_{e\mu}$ (left
panel), and $\mathcal{A}_{\tau}^{0}$ and $R_{e\tau}$ (right panel) within the
$2\sigma$ uncertainty ranges, for $m_{\chi}=200$ GeV. Again the darker shades
of orange correspond to higher values for $y_{\mu}$ and $y_{\tau}$ as shown in
the legend. Both calculations are the same from the theoretical side and there
is no noticeable difference coming from the different masses. However, a
direct comparison between the two plots should be done with care, since,
unlike the muon case, the SM prediction for $R_{e\tau}$ is outside the
$1\sigma$ experimental uncertainty range. Therefore, a correlation between
constraints on the tau and experimental precision cannot be made.

##  5 Conclusions

We have proposed a set of precision muon-related observables that serve as a
tool to constrain new physics models. Taking advantage of the precision
measurements performed by the LEP collaborations on the $Z$-boson pole, it is
possible to perceive the contribution of new physics to measured SM processes
via quantum corrections. The asymmetry parameter ($\mathcal{A}_{\mu}^{0}$) and
the ratio of lepton hadronic decay rates ($R_{e\mu}$) are extracted from the
$e^{+}e^{-}\to\mu^{+}\mu^{-}$ process. With these observables at hand, we can
find the limits on the variations of the $Z\mu^{+}\mu^{-}$ couplings when one-
loop contributions from new physics are included. We provide a simple recipe
that can be used for any extension of the SM.

We have discussed in detail the impact of the new bounds on two particular
models in which a new Dark Sector is added to the SM. The new sector includes
a vector like fermion and two scalars. In one scenario, dubbed Model 5 the
scalars are $SU(2)$ singlets and the fermion is an $SU(2)$ doublet. In the
other scenario, Model 3, the scalars are $SU(2)$ doublets and the fermion is
an $SU(2)$ singlet. While in the first scenario, the new constraints are very
mild, for Model 3 they constitute a real bound on the low mass region of the
scalars. This is especially important when combined with other constraints

These bounds are particularly relevant for models with an enhanced coupling to
the muons as is the case of the ones that attempt to solve the $g-2$ anomaly
via new physics loop contributions.

#### Acknowledgments

GL and RS are partially supported by the Portuguese Foundation for Science and
Technology (FCT) under Contracts no. UIDB/00618/2020, UIDP/00618/2020,
CERN/FIS-PAR/0025/2021 and CERN/FIS-PAR/0021/2021. The work of A.M. and J.P.S.
is supported in part by FCT under Contracts CERN/FIS-PAR/0008/2019, PTDC/FIS-
PAR/29436/2017, UIDB/00777/2020, and UIDP/00777/2020; these projects are
partially funded through POCTI (FEDER), COMPETE, QREN, and the EU.

## References

  * [1] LHCb, R. Aaij et al.,  Phys. Rev. D 108, 032002 (2023), 2212.09153.
  * [2] LHCb, R. Aaij et al.,  Phys. Rev. Lett. 131, 051803 (2023), 2212.09152.
  * [3] LHCb, R. Aaij et al.,  Nature Phys. 18, 277 (2022), 2103.11769.
  * [4] LHCb, R. Aaij et al.,  Phys. Rev. Lett. 122, 191801 (2019), 1903.09252.
  * [5] LHCb, R. Aaij et al.,  JHEP 08, 055 (2017), 1705.05802.
  * [6] Belle, A. Abdesselam et al.,  Phys. Rev. Lett. 126, 161801 (2021), 1904.02440.
  * [7] BELLE, S. Choudhury et al.,  JHEP 03, 105 (2021), 1908.01848.
  * [8] B. Bhattacharya, D. London, J. M. Cline, A. Datta, and G. Dupuis,  Phys. Rev. D 92, 115012 (2015), 1509.04271.
  * [9] J. Kawamura, S. Okawa, and Y. Omura,  Phys. Rev. D 96, 075041 (2017), 1706.04344.
  * [10] J. M. Cline and J. M. Cornell,  Phys. Lett. B 782, 232 (2018), 1711.10770.
  * [11] D. G. Cerdeño, A. Cheek, P. Martín-Ramiro, and J. M. Moreno,  Eur. Phys. J. C 79, 517 (2019), 1902.01789.
  * [12] B. Barman, D. Borah, L. Mukherjee, and S. Nandi,  Phys. Rev. D 100, 115010 (2019), 1808.06639.
  * [13] L. Darmé, M. Fedele, K. Kowalska, and E. M. Sessolo,  JHEP 08, 148 (2020), 2002.11150.
  * [14] G. Arcadi, L. Calibbi, M. Fedele, and F. Mescia,  Phys. Rev. Lett. 127, 061802 (2021), 2104.03228.
  * [15] M. Becker, D. Döring, S. Karmakar, and H. Päs,  Eur. Phys. J. C 81, 1053 (2021), 2103.12043.
  * [16] D. Huang, A. P. Morais, and R. Santos,  Phys. Rev. D 102, 075009 (2020).
  * [17] R. Capucha, D. Huang, T. Lopes, and R. Santos,  Phys. Rev. D 106, 095032 (2022), 2207.11556.
  * [18] Muon g-2, D. P. Aguillard et al.,  (2023), 2308.06230.
  * [19] Particle Data Group, M. Tanabashi et al.,  Phys. Rev. D 98, 030001 (2018).
  * [20] T. P. Gorringe and D. W. Hertzog,  Prog. Part. Nucl. Phys. 84, 73 (2015), 1506.01465.
  * [21] T. Aoyama et al.,  Phys. Rept. 887, 1 (2020), 2006.04822.
  * [22] Muon g-2, G. W. Bennett et al.,  Phys. Rev. D 73, 072003 (2006), hep-ex/0602035.
  * [23] J. H. Field,  Mod. Phys. Lett. A 13, 1937 (1998), hep-ph/9801355.
  * [24] H. E. Haber and H. E. Logan,  Phys. Rev. D 62, 015011 (2000).
  * [25] D. Fontes, L. Lavoura, J. C. Romão, and J. P. Silva,  Nucl. Phys. B 958, 115131 (2020), 1910.11886.
  * [26] D. Jurčiukonis and L. Lavoura,  JHEP 07, 195 (2021), 2103.16635.
  * [27] ALEPH, DELPHI, L3, OPAL, SLD, LEP Electroweak Working Group, SLD Electroweak Group, SLD Heavy Flavour Group, S. Schael et al.,  Phys. Rept. 427, 257 (2006), hep-ex/0509008.
  * [28] E. J. Chun and J. Kim,  JHEP 07, 110 (2016), 1605.06298.
  * [29] Particle Data Group, R. L. Workman et al.,  PTEP 2022, 083C01 (2022).
  * [30] W. F. L. Hollik,  Fortsch. Phys. 38, 165 (1990).
  * [31] M. Bohm, H. Spiesberger, and W. Hollik,  Fortsch. Phys. 34, 687 (1986).
  * [32] D. Albert, W. J. Marciano, D. Wyler, and Z. Parsa,  Nucl. Phys. B 166, 460 (1980).
  * [33] W. Hollik,  Adv. Ser. Direct. High Energy Phys. 14, 37 (1995).
  * [34] M. E. Peskin and T. Takeuchi,  Phys. Rev. D 46, 381 (1992).
  * [35] A. K. Grant and T. Takeuchi,  (1998), hep-ph/9807413.
  * [36] J. Bernabeu, G. A. Gonzalez-Sprinberg, and J. Vidal,  Weak dipole moments at $e^{+}e^{-}$ colliders,  in Ringberg Workshop on Perspectives for Electroweak Interactions in e+ e- Collisions, pp. 0329–342, 1995, hep-ph/9505223.
  * [37] L3, P. Achard et al.,  Phys. Lett. B 517, 75 (2001), hep-ex/0107015.
  * [38] CMS, A. M. Sirunyan et al.,  Phys. Rev. D 100, 112007 (2019), 1909.04114.
  * [39] S. Bißmann, G. Hiller, C. Hormigos-Feliu, and D. F. Litim,  Eur. Phys. J. C 81, 101 (2021), 2011.12964.
  * [40] G. Belanger, A. Mjallal, and A. Pukhov,  Eur. Phys. J. C 81, 239 (2021), 2003.08621.
  * [41] M. Algueró et al.,  (2023), 2304.07330.
  * [42] Planck Collaboration et al.,  A&A 641, A6 (2020).
  * [43] LZ collaboration,  First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment, 2022.
  * [44] W. Grimus, L. Lavoura, O. M. Ogreid, and P. Osland,  J. Phys. G 35, 075001 (2008), 0711.4022.
  * [45] A. Semenov,  Lanhep - a package for automatic generation of feynman rules from the lagrangian. updated version 3.2, 2014, 1412.5016.
  * [46] N. D. Christensen and C. Duhr,  Comput. Phys. Commun. 180, 1614 (2009), 0806.4194.
  * [47] J. Kublbeck, M. Bohm, and A. Denner,  Comput. Phys. Commun. 60, 165 (1990).
  * [48] R. Mertig, M. Bohm, and A. Denner,  Comput. Phys. Commun. 64, 345 (1991).
  * [49] V. Shtabovenko, R. Mertig, and F. Orellana,  Comput. Phys. Commun. 207, 432 (2016), 1601.01167.
  * [50] G. Passarino and M. J. G. Veltman,  Nucl. Phys. B 160, 151 (1979).
  * [51] T. Hahn and M. Perez-Victoria,  Comput. Phys. Commun. 118, 153 (1999), hep-ph/9807565.