# Hybrid Quantum Benders’ Decomposition For Mixed-integer Linear Programming

Zhongqi Zhao1,2, Lei Fan2, and Zhu Han1 Dept. of Electrical and Computer
Engineering1 and Dept. of Engineering Technology2
University of Houston, TX, USA

###### Abstract

The Benders’ decomposition algorithm is a technique in mathematical
programming for complex mixed-integer linear programming (MILP) problems with
a particular block structure. The strategy of Benders’ decomposition can be
described as a strategy of divide and conquer. The Benders’ decomposition
algorithm has been employed in a variety of applications such as
communication, networking, and machine learning. However, the master problem
in Benders’ decomposition is still NP-hard, which motivates us to employ
quantum computing. In the paper, we propose a hybrid quantum Benders’
decomposition algorithm. We transfer the Benders’ decomposition’s master
problem into the quadratic unconstrained binary optimization (QUBO) model and
solve it by the state-of-the-art quantum annealer. Then, we analyze the
computational results and discuss the feasibility of the proposed algorithm.
Due to our reformulation in the master problem in Benders’ decomposition, our
hybrid algorithm, which takes advantage of both classical and quantum
computers, can guarantee the solution quality for solving MILP problems.

###### Index Terms:

Benders’ Decomposition, Mixed-integer Linear Programming, Optimization,
Quantum Computing, Communication, Networking

##  I Introduction

In recent years, mixed-integer linear programming (MILP) is widely employed in
many fields includes but is not limited to planning and scheduling [1],
transportation and telecommunications [2], energy and resource management
[3][4], network design [5], and vehicle routing [6]. However, because MILP is
an NP-hard problem [7] [8], it is desirable to have powerful solutions to
solve large-scale MILP problems.

Benders’ decomposition algorithm [9] provides an efficient way to solve the
MILP problem. In recent years, the Benders’ decomposition method has become a
popular algorithm. It exploits the structure of the MILP problem and
successfully reduces the computation workload. Benders’ decomposition divides
a MILP problem into a master problem and a subproblem. The subproblem is a
linear programming model that has strong duality. During the iterative solving
process, the optimality cuts and feasibility cuts obtained from the solution
of the subproblem will be added to the master problem. The iterative process
will stop when the gap between the upper bound provided by the master problem
and the lower bound provided by the subproblem is sufficiently small. However,
the master problem in Benders’ decomposition is still NP-hard, which motivates
us to employ quantum computing.

Due to the characteristics of parallel computing, the quantum computer can
simultaneously calculate the final result with every possible model input.
Therefore, quantum computers are often considered to have surpassed the
computational speed of traditional computers. The quantum advantage has been
proved by quantum algorithms, including Deutch-Jozsa Algorithm [10], Grover’s
Algorithm [11], Shor’s Algorithm [12]. In recent years, leading quantum
computer companies such as D-Wave, IBM, Google, IonQ, Honeywell have made
significant progress in hardware design. Normally, their quantum computers are
in one of the two groups (i.e., analog quantum model, universal quantum gate
model). D-Wave currently provides the quantum computer with the largest number
of qubits on the market.

By deploying the Ising model, the D-Wave’s quantum annealer computer can solve
the problem formulated by the quadratic unconstrained binary optimization
(QUBO) model. It depicts the energy state with coupling qubits interaction and
externally applied fields [13].

The power of quantum computers and challenges in MILP problems inspire us to
design a hybrid quantum Benders’ decomposition algorithm by jointly using
quantum computing and classical computing techniques. However, there are
several obstacles when we try to use Benders’ composition to combine quantum
and traditional computing. The first difficulty is how to convert the problem
into an integer linear programming (ILP) problem recognized by the quantum
computer. The Second difficulty is how to convert the NP-hard ILP problem into
a QUBO model as an input to a quantum computer. In addition, the third
difficulty is how to reformat the output of quantum computers from the binary
solution to the decimal numeral system.

To overcome the above challenges, this paper reformulates the master problem
as an ILP model in Benders’ decomposition under the MILP problem. We transfer
the master problem, which is an ILP model, into a QUBO form. Although the
master problem in the Benders’ decomposition is still NP-hard, we overcome
this issue by proposing a hybrid Benders’ decomposition algorithm with the
D-Wave hybrid quantum computer for the MILP problem. Finally, we study a
simple case to evaluate the performance of the D-Wave hybrid solver in solving
the MILP problem. The contributions of this paper are summarized as follows.

  * •

We propose a hybrid quantum Benders’ decomposition algorithm to find the
solution for the MILP problem. Our hybrid quantum Benders’ decomposition
algorithm converges and returns the correct final result as the classical
algorithm does.

  * •

We propose an ILP model for the master problem in Benders’ decomposition for
the MILP problem. We reformulate the ILP model as a QUBO model recognized by
the quantum annealing machine.

  * •

We employ the quantum computer provided by D-Wave to solve the MILP problem by
their Leap™ quantum cloud service. Our experiments show the possibility of
using quantum computing to solve the MILP.

The rest of this paper is organized as follows. Section II introduces the
basics about the MILP problem and the Benders’ decomposition. Section III
illustrates our hybrid quantum Benders’ Decomposition algorithm. Section IV
validates our algorithm by showing the corresponding simulation. Finally,
Section V concludes the whole paper.

##  II MILP and Benders’ Decomposition Basics

###  II-A Mixed-integer Linear Programming

MILP has been widely adopted optimization problems that include but are not
limited to communication and networks. In the field of communication, for
instance, it has been employed to the resource allocation in wavelength
division multiple access (WDMA) [14], multi-access edge computing (MEC) [15],
offloading for fog computing [16], etc. Furthermore, it plays an essential
role in the field of network research, such as minimizing the network delay
[17], scheduling virtual network re-configurations [18], and improving the
supply chain network [19]. Besides civilian applications, it is no doubt that
MILP has been adopted for military purposes such as assisting military pilot
training [20]. These demonstrate that MILP is a powerful tool for either
classical problems or future popular topics such as communication networks and
resource allocation. Consider a MILP model as follows:

$\displaystyle\max_{x,y}$ $\displaystyle c^{\intercal}x+h^{\intercal}y$ (1)
s.t. $\displaystyle Ax+Gy\leq b,$ $\displaystyle x\in X,\
x\in\left\\{0,1\right\\}^{n},\ y\in\mathbb{R}^{p}_{+}.$

Here $x$ is a vector with binary variables and $A$ is its corresponding
coefficient matrix in constraints. $y$ is a vector with non-negative
continuous variables and $G$ is its corresponding coefficient matrix in
constraints. $c^{\intercal}$ is the coefficient vector of variable $x$ in the
objective function and $h^{\intercal}$ is the coefficient vector of variable
$y$ in the objective function.

Figure 1: The structure of hybrid Benders’ decomposition algorithm for MILP
using both quantum and traditional computing

###  II-B Benders’ Decomposition for MILP

We first provide a brief introduction of the Benders’ Decomposition. In the
MILP problem (1), for each possible choice of $\bar{x}\in X$, we find the best
choice for $y$ by solving a linear program. So we regard $y$ as a function of
$x$. Then we replace the contribution of $y$ to the objective with a scalar
variable representing the value of the best choice for a given $\bar{x}$. We
start with a crude approximation to the contribution of $y$ and then generate
a sequence of dual solutions to tighten up the approximation. Hence, the
original MILP problem in (1) can be written equivalently to a master problem
as follows.

$\displaystyle\mbox{Master Problem:}\ \max_{x,t}$ $\displaystyle
c^{\intercal}x+t$ (2) s.t.
$\displaystyle\left(b-Ax\right)^{\intercal}u^{k}\geq t\quad\textrm{for}\
k\in\hat{K},$ $\displaystyle\left(b-Ax\right)^{\intercal}r^{j}\geq
0\quad\textrm{for}\ j\in\hat{J},$ $\displaystyle x\in X,\
x\in\left\\{0,1\right\\}^{n},\ t\in\mathbb{R}.$

We denote $K$ and $J$ as the extreme points $u^{k}$ and extreme rays $r^{k}$
of the dual polyhedron $Q=\left\\{u\in\mathbb{R}^{m}_{+}\mid
G^{\intercal}u\geq h\right\\}$ generated by the linear programming’s duality,
which is called the subproblem. The subproblem is as follows.

$\displaystyle\mbox{Subproblem:}\ z_{LP}(x)$ $\displaystyle=\min_{u}\
\left(b-Ax\right)^{\intercal}u$ (3) s.t. $\displaystyle G^{\intercal}u\geq h,$
$\displaystyle u\in\mathbb{R}^{m}_{+}.$

In the subproblem, if the inner product between $(b-Ax)$ and any dual ray
$r^{j^{\prime}}$ is negative then $z_{LP}(x)=-\infty$. Equivalently, in this
situation, the dual problem of Problem (3) is infeasible. Then, $x$ does not
allow a feasible solution to the original mixed-integer problem (1). Thus, we
have a new feasibility cut, i.e.,

$(b-Ax)^{\intercal}r^{j^{\prime}}\geq 0,\quad\ j^{\prime}\in J.$ (4)

If $x$ satisfies (4), then we yield an extreme point $u^{k^{\prime}}$ and the
value of $z_{LP}(x)$ is given by

$z_{LP}(x)=\min_{k^{\prime}\in K}\
\left(b-Ax\right)^{\intercal}u^{k^{\prime}}.$ (5)

Thus, problem (1) is written equivalently as problem (2) by denoting
$z_{LP}(x)$ as a continuous real number variable $t$. Equation (5) posts a new
optimality cut in master problem. In addition to the extreme points and
extreme rays, we use $\hat{K}$ and $\hat{J}$ to denote the current known
extreme points and extreme rays of $Q$, respectively, where $\hat{K}\subseteq
K$ and $\hat{J}\subseteq J$.

##  III Hybrid Quantum Benders’ Decomposition

Quantum computer is a powerful tool to solve NP-hard integer problems.
Accordingly, we introduce the quantum computer to solve the master problem in
(2) since it has a special formulation that can convert to an ILP problem. The
subproblem is a linear programming model which can be solved well on the
classical computer. As Fig. 1 shows, we design a hybrid algorithm in which the
master problem is solved by the quantum computing techniques and the
subproblem is solved by the classical computing techniques. The two different
computers are communicating with each other. The subproblem returns new
optimality and feasibility cuts to tighten the bounds of the master problem.
The integer solution returned from the master problem gives the subproblem a
direction to find new optimality and feasibility cuts.

TABLE I: Table of Common Constraint-Penalty Pairs Constraint | Equivalent Penalty
---|---
$x_{1}+x_{2}=1$ | $P(x_{1}+x_{2}-1)^{2}$
$x_{1}+x_{2}\geq 1$ | $P(1-x_{1}-x_{2}+x_{1}x_{2})^{2}$
$x_{1}+x_{2}\leq 1$ | $P(x_{1}x_{2})$
$x_{1}+x_{2}+x_{3}\leq 1$ | $P(x_{1}x_{2}+x_{1}x_{3}+x_{2}x_{3})$

###  III-A Quantum Formulation

Quantum annealers are able to solve the optimization problem in a QUBO
formulation. To leverage state-of-art quantum annealers provided by D-Wave,
the ILP problem has to be converted to the corresponding QUBO formulation. The
definition of QUBO are (6) and (7). Let
$f$:$\left\\{0,1\right\\}^{n}\rightarrow\mathbb{R}$ be a quadratic polynomial
over binary variables $x$ of length $n$ and $q_{ij}$ is corresponding cost
coefficient of $x_{i}x_{j}$. We get equation (6)

$f_{Q}\left(\mathbf{x}\right)=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}q_{ij}x_{i}x_{j}=\mathbf{x}^{\intercal}\mathbf{Q}\mathbf{x}.$
(6)

The QUBO problem consists of finding a binary vector $\mathbf{x}^{*}$ that is
minimal with respect to $f$ among all other binary vectors, namely,

$\mathbf{x}^{*}={\arg\min}_{\mathbf{x}}\ f\left(\mathbf{x}\right).$ (7)

Here $\mathbf{x}$ is a vector of binary variables where the length of $n$, and
$\mathbf{Q}$ is either an upper-diagonal matrix or a symmetric matrix. Since
(6) is an unconstrained optimization model, we need to reformulate our
constrained ILP as unconstrained QUBO by using penalties. Next, we get the
optimal solution by finding the best penalty coefficients of the constraints.
The principles of transforming classical constraints to their equivalent
penalties are displayed in the TABLE I. Here $x_{1}$, $x_{2}$ and $x_{3}$ are
binary variables. $s_{i}$ is a binary slack variable. $a_{l}$ is the
coefficient for the corresponding slack variable. $b$ is a constant. $P$ is a
user-defined penalty coefficient.

###  III-B Variable Representation

Now consider problem (2), the initial binary decision variables make up the
vector $\mathbf{x}\in X$ with length of $n$. $\bar{t}\in\mathbb{R}$. In order
to reformulate the master problem into the QUBO formulation, we needs to
represent the continuous variable using binary bits. We use a binary vector
$\mathbf{w}$ with length of $M$ bits to replace continuous variable $\bar{t}$.
In general, $\bar{t}$ requires the binary numeral system assigning $M$ bits to
replace continuous variable $\bar{t}$. Then we can recover the $\bar{t}$ by

$\displaystyle\bar{t}$
$\displaystyle=\sum^{\bar{m}_{+}}_{i=-\underline{m}}2^{i}w_{\left(i+\underline{m}\right)}-\sum^{\bar{m}_{-}}_{j=0}2^{j}w_{j+1+\underline{m}+\bar{m}_{+}}$
(8) $\displaystyle=\bar{t}\left(\mathbf{w}\right).$

In (8), $\bar{m}_{+}$ is the number of bits that assigned to represent
positive integer part, $\underline{m}$ is the number of bits assigned to
represent the positive decimal part, and $\bar{m}_{-}$ is the number of bits
that represent the negative integer part.

###  III-C QUBO Setup

Let’s reconsider problem (2) and replace $\bar{t}$ in (8). The new problem is
only depend on binary decision variables $\mathbf{x}$ and $\mathbf{w}$. Then,
the new master problem is

$\displaystyle\max_{\mathbf{x},\mathbf{w}}$ $\displaystyle
c^{\intercal}\mathbf{x}+\sum^{\bar{m}_{+}}_{i=-\underline{m}}2^{i}w_{\left(i+\underline{m}\right)}-\sum^{\bar{m}_{-}}_{j=0}2^{j}w_{j+1+\underline{m}+\bar{m}_{+}}$
(9) s.t.
$\displaystyle\left(b-A\mathbf{x}\right)^{\intercal}u^{k}\geq\bar{t}\left(\mathbf{w}\right),\quad\textrm{for}\
k\in\hat{K},$ $\displaystyle\left(b-A\mathbf{x}\right)^{\intercal}r^{j}\geq
0,\quad\textrm{for}\ j\in\hat{J},$ $\displaystyle\mathbf{x}\in
X,\quad\mathbf{x}\in\left\\{0,1\right\\}^{n},$ $\displaystyle\mathbf{w}\in
W,\quad\mathbf{w}\in\left\\{0,1\right\\}^{M}.$

Then we create a new binary decision variable collection
$\mathbf{x}^{\prime}=\left\\{\mathbf{w},\mathbf{x}\right\\}$ to set up our
QUBO matrix, where the length of the binary vector $\mathbf{x}^{\prime}$ is
$n+M$. According to the rule of setting up QUBO, we will transfer our
objective and constraints one by one to be a symmetric matrix.

####  III-C1 Objective Function

Because the quantum computer only accepts a quadratic polynomial over binary
variables. Following the principle of the QUBO formulation, we convert the
objective function to be a QUBO matrix as follows.

$\displaystyle
c^{\intercal}\mathbf{x}+\sum^{\bar{m}_{+}}_{i=-\underline{m}}2^{i}w_{\left(i+\underline{m}\right)}-\sum^{\bar{m}_{-}}_{j=0}2^{j}w_{j+1+\underline{m}+\bar{m}_{+}}$
(10) $\displaystyle\Rightarrow$
$\displaystyle\mathbf{Q}_{obj}=\sum^{\bar{m}_{+}}_{i=-\underline{m}}w_{\left(i+\underline{m}\right)}2^{i}w_{\left(i+\underline{m}\right)}+\mathbf{x}^{\intercal}\textrm{diag}\left(c\right)\mathbf{x}$
$\displaystyle-\sum^{\bar{m}_{-}}_{j=0}w_{j+1+\underline{m}+\bar{m}_{+}}2^{j}w_{j+1+\underline{m}+\bar{m}_{+}}.$

####  III-C2 Optimality Cuts

Similarly, following the principle of constraint-penalty pairs in TABLE I, we
not only introduce the penalty but also convert the optimality cuts constraint
to be a QUBO matrix as follow.

$\displaystyle\bar{t}\left(\mathbf{w}\right)+\left(u^{k}\right)^{\intercal}A\mathbf{x}\leq
b^{\intercal}u^{k},\ \textrm{for}\ k\in\hat{K}.$ $\displaystyle\Rightarrow$
$\displaystyle
P_{k}\left(\bar{t}\left(\mathbf{w}\right)+\left(u^{k}\right)^{\intercal}A\mathbf{x}+\sum_{l=0}^{\bar{l}^{K}}2^{l}s^{K}_{kl}-b^{\intercal}u^{k}\right)^{2},$
where
$\displaystyle\bar{l}^{K}=\left\lceil\log_{2}\left(b^{\intercal}u^{k}-\min_{\mathbf{w},\mathbf{x}}\left(\bar{t}\left(\mathbf{w}\right)+\left(u^{k}\right)^{\intercal}A\mathbf{x}\right)\right)\right\rceil.$

####  III-C3 Feasibility Cuts

Similarly, following the principle of constraint-penalty pairs in TABLE I, we
convert the feasibility cuts constraint to be a QUBO matrix with penalties as
follows:

$\displaystyle\left(r^{j}\right)^{\intercal}A\mathbf{x}\leq
b^{\intercal}r^{j},\quad\textrm{for}\ j\in\hat{J}.$ $\displaystyle\Rightarrow$
$\displaystyle
P_{j}\left(\left(r^{j}\right)^{\intercal}A\mathbf{x}+\sum_{l=0}^{\bar{l}^{J}}2^{l}s^{J}_{kl}-b^{\intercal}r^{j}\right)^{2},$
where
$\displaystyle\bar{l}^{J}=\left\lceil\log_{2}\left(b^{\intercal}r^{j}-\min_{\mathbf{x}}\left(\left(r^{j}\right)^{\intercal}A\mathbf{x}\right)\right)\right\rceil.$

###  III-D Proposed Algorithm

We deploy the D-Wave solver in our proposed heuristic algorithm to solve the
model. In addition, we need to carefully tune the penalties for a decent QUBO
model. An extremely large penalty may cause quantum annealer malfunctioning
since it will explode the coefficients. Similarly, a soft penalty may make
quantum annealer ignore the corresponding constraints. As long as the penalty
is well-tuned, the quantum solver will give the right answer with a relatively
high probability. Therefore, we summarize our proposed hybrid quantum Benders’
decomposition algorithm as Algorithm 1.

Algorithm 1 Hybrid Quantum Benders’ Decomposition Algorithm

Initial sets $\hat{K}$ of extreme points and $\hat{J}$ of extreme rays of $Q$

$\bar{t}$ $\leftarrow$ $+\infty$

$\underline{t}$ $\leftarrow$ $-\infty$

while $\mid\bar{t}-\underline{t}\mid\geq\epsilon$ do

$\mathbf{P}$ $\leftarrow$ Appropriate penalties numbers or arrays

$\mathbf{Q}$ $\leftarrow$ Reformulate both objective and constraints in (2)
and construct the QUBO formulation by using corresponding rules

$\mathbf{x}^{\prime}$ $\leftarrow$ Solve problem (6) by quantum computer.

$\bar{t}$ $\leftarrow$ Extract $\mathbf{w}$ and replace the $\bar{t}$ with
$\bar{t}\left(\mathbf{w}\right)$ (8)

$z_{LP}(\mathbf{x})$ $\leftarrow$ Solve the problem (3)

$\underline{t}$ $\leftarrow$ $z_{LP}(\mathbf{x})$

if $z_{LP}(\mathbf{x})=-\infty$ then

An extreme ray $j$ of $Q$ has been found.

$\hat{J}$ = $\hat{J}$ $\cup$ $\left\\{j\right\\}$

else if $z_{LP}(\mathbf{x})<\bar{t}$ and $\bar{t}\neq+\infty$ then

An extreme point $k$ of $Q$ has been found.

$\hat{K}$ = $\hat{K}$ $\cup$ $\left\\{k\right\\}$

end if

end while

return $\bar{t}$, $\mathbf{x}$

##  IV NUMERICAL VALIDATION

We validate the proposed algorithm by running on a hybrid D-Wave quantum
processing unit (QPU). Because the input of interest to practice is too large
to fit onto current-model QPUs and be solved directly by quantum annealing. We
choose to use the hybrid model instead of the pure QPU. The reason for using
the hybrid model is that the hybrid solver overcomes a lot of input size
barriers and allows the QPU to accept a large input. More details can be found
in [21]. We accessed the D-Wave system by Leap™ quantum cloud service.

###  IV-A Example Setup

In our simulation, we consider a simple problem to test our proposed quantum
algorithm, where $\mathbf{x}\in X$, $\mathbf{x}\in\left\\{0,1\right\\}^{2}$,
$\mathbf{y}\in Y$, $\mathbf{y}\geq 0$

$A=\begin{bmatrix}0&0\\\ 0&0\\\ 0&0\\\ 0&0\\\ -1&-1\\\ -1&0\\\ -1&0\\\ 0&-1\\\
0&-1\\\ \end{bmatrix},\ G=\begin{bmatrix}1&0&1&0\\\ 1&0&0&1\\\ 0&1&1&0\\\
0&1&0&1\\\ 0&0&0&0\\\ 1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\\\
\end{bmatrix},\ b=\begin{bmatrix}1\\\ 1\\\ 1\\\ 1\\\ -1\\\ 0\\\ 0\\\ 0\\\ 0\\\
\end{bmatrix},$ $h^{\intercal}=\begin{bmatrix}8&9&5&6\\\ \end{bmatrix},\quad
c^{\intercal}=\begin{bmatrix}-15&-10\\\ \end{bmatrix}.$

###  IV-B Result

Fig. 2 shows how our proposed hybrid quantum Benders’ decomposition algorithm
adds optimality and feasibility cuts to get the solution and corresponding
$\bar{t}$. In the example, the algorithm takes 4 rounds to let $t$ converge.
Four cuts are added iteratively in the space and tighten the feasible region.
The algorithm can find the optimal solution in each round. Therefore, Fig. 2
demonstrates that our algorithm is reliable and efficient. In Fig.
LABEL:fig:tutl, the hidden dashed line denotes that the lower bound in the
corresponding round is negative infinity. As we can see in the graph, the
upper bound and lower bound converge. Our algorithm only takes two rounds to
find the non-negative infinity lower bound. This result proves that our
proposed algorithm is mathematically consistent with the classic Benders’
decomposition algorithm. In other words, as long as the classic Benders’
decomposition algorithm can solve it, our algorithm can at least achieve the
same.

Figure 2: Benders’ Decomposition Process

##  V Conclusion

In this paper, we not only reformulate the master problem of Benders’
decomposition for the MILP problem as an ILP model but also successfully
convert it to a QUBO model. our hybrid quantum Benders’ decomposition
algorithm converges and returns the correct final result as the classical
algorithm does. In addition to that, our algorithm also guarantees the
solution quality for solving the MILP problem. In our simulation, we solved an
NP-hard MILP problem by using the hybrid quantum computer provided by D-Wave.
Therefore, we can conclude that the quantum computer can potentially replace
the classical computer in solving the NP-hard master problem of Benders’
decomposition algorithm in MILP problems.

## References

  * [1] S. P. Canto, “Application of Benders’ decomposition to power plant preventive maintenance scheduling,” _European Journal of Operational Research_ , vol. 184, no. 2, pp. 759–777, Jan. 2008.
  * [2] A. M. Costa, “A survey on benders decomposition applied to fixed-charge network design problems,” _Computers & Operations Research_, vol. 32, no. 6, pp. 1429–1450, Jun. 2005.
  * [3] X. Cai, D. C. McKinney, L. S. Lasdon, and D. W. Watkins, “Solving large nonconvex water resources management models using generalized benders decomposition,” _Operations Research_ , vol. 49, no. 2, pp. 235–245, Apr. 2001.
  * [4] J. L. Zhang and K. Ponnambalam, “Hydro energy management optimization in a deregulated electricity market,” _Optimization and Engineering_ , vol. 7, no. 1, pp. 47–61, Mar. 2006.
  * [5] B. Fortz and M. Poss, “An improved Benders decomposition applied to a multi-layer network design problem,” _Operations Research Letters_ , vol. 37, no. 5, pp. 359–364, Sep. 2009.
  * [6] A. I. Corréa, A. Langevin, and L.-M. Rousseau, “Scheduling and routing of automated guided vehicles: A hybrid approach,” _Computers & Operations Research_, vol. 34, no. 6, pp. 1688–1707, Jun. 2007.
  * [7] A. Bulut and T. K. Ralphs, “On the complexity of inverse mixed integer linear optimization,” _arXiv:2104.09002 [cs, math]_ , Apr. 2021, arXiv: 2104.09002. [Online]. Available: http://arxiv.org/abs/2104.09002
  * [8] C. H. Papadimitriou and M. Yannakakis, “The complexity of facets (and some facets of complexity),” in _Proceedings of the fourteenth annual ACM symposium on Theory of computing_ , ser. STOC ’82. New York, NY: ACM, May 1982, pp. 255–260.
  * [9] A. M. Geoffrion, “Generalized benders decomposition,” _Journal of Optimization Theory and Applications_ , vol. 10, no. 4, pp. 237–260, Oct. 1972.
  * [10] D. Deutsch and R. Jozsa, “Rapid solution of problems by quantum computation,” _Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences_ , vol. 439, no. 1907, pp. 553–558, Dec. 1992.
  * [11] L. K. Grover, “A fast quantum mechanical algorithm for database search,” in _Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing_ , ser. STOC ’96. New York, NY: ACM, Jul. 1996, p. 212–219.
  * [12] P. W. Shor, “Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer,” _SIAM Journal on Computing_ , vol. 26, no. 5, pp. 1484–1509, Oct. 1997.
  * [13] E. Ising, “Beitrag zur Theorie des Ferromagnetismus,” _Zeitschrift für Physik_ , vol. 31, no. 1, pp. 253–258, Feb. 1925.
  * [14] A. S. Elgamal, O. Z. Alsulami, A. A. Qidan, T. E. H. El-Gorashi, and J. M. H. Elmirghani, “Q-learning algorithm for resource allocation in WDMA-based optical wireless communication networks,” _arXiv:2105.10729 [cs, eess]_ , May 2021, arXiv: 2105.10729. [Online]. Available: http://arxiv.org/abs/2105.10729
  * [15] Y. Yu, X. Bu, K. Yang, H. Yang, X. Gao, and Z. Han, “UAV-Aided Low Latency Multi-Access Edge Computing,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 5, pp. 4955–4967, May 2021.
  * [16] Z. Wu, B. Li, Z. Fei, Z. Zheng, B. Li, and Z. Han, “Energy-Efficient Robust Computation Offloading for Fog-IoT Systems,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 4, pp. 4417–4425, Apr. 2020\.
  * [17] W. Xuan, Z. Zhao, L. Fan, and Z. Han, “Minimizing Delay in Network Function Visualization with Quantum Computing,” _arXiv:2106.10707 [cs, eess, math]_ , Jun. 2021, arXiv: 2106.10707. [Online]. Available: http://arxiv.org/abs/2106.10707
  * [18] X. Pan, S. Zhao, H. Yang, S. Tang, and Z. Zhu, “Scheduling Virtual Network Reconfigurations in Parallel in Hybrid Optical/Electrical Datacenter Networks,” _Journal of Lightwave Technology_ , vol. 39, no. 17, pp. 5371–5382, Sep. 2021.
  * [19] P. Santander, F. A. Cruz Sanchez, H. Boudaoud, and M. Camargo, “Closed loop supply chain network for local and distributed plastic recycling for 3D printing: a MILP-based optimization approach,” _Resources, Conservation and Recycling_ , vol. 154, p. 104531, Mar. 2020.
  * [20] V. Mak-Hau, B. Hill, D. Kirszenblat, B. Moran, V. Nguyen, and A. Novak, “A simultaneous sequencing and allocation problem for military pilot training: Integer programming approaches,” _Computers & Industrial Engineering_, vol. 154, p. 107161, Apr. 2021.
  * [21] “D-Wave Hybrid Solver Service: An Overview.” [Online]. Available: https://www.dwavesys.com/resources/white-paper/d-wave-hybrid-solver-service-an-overview/


# Two-to-two processes at an electron-muon collider.

Antonio O. Bouzas and F. Larios
Departamento de Física Aplicada, CINVESTAV-Mérida,
A.P. 73, 97310 Mérida, Yucatán, México
abouzas<EMAIL_ADDRESS>corresponding author.

###### Abstract

Based on a recent proposal to build an electron-muon collider, we study two-
to-two production processes $e^{-}\mu^{+}\rightarrow f\bar{f}$, $\gamma\gamma$
that originate from dimension 6 and 8 operators. We compare the sensitivity to
those effective couplings obtained at the collider with that of low energy
measurements of $\mu\rightarrow e\gamma$, $\mu\rightarrow e\bar{e}e$ and
$\mu\rightarrow e$ conversion that have recently been reported in the
literature. Whereas for the production of first family fermions the
sensitivity of the collider processes is much weaker, for the second and third
familiy fermions it is similar or stronger than that of low-energy processes.
In the case of $e^{-}\mu^{+}\rightarrow\gamma\gamma$, the sensitivity to a
dimension 8 contact operator turns out to be the strongest in comparison.

##  1 Introduction

With the main objective of measuring Higgs boson properties, there has been
renewed interest in muon colliders [1, 2, 3]. At the low energy level, there
have been numerous works concerning $(g-2)_{\mu}$ from the theoretical and
experimental point of view [4, 5, 6, 7, 8]. In addition, the recent
confirmation of a significant deviation from the Standard Model (SM) value of
$(g-2)_{\mu}$ at Fermilab [9] has also encouraged the proposal of new physics
models [10]-[17] as well as the calculation of higher order QED contributions
to the $e\mu\to e\mu$ scattering process [18, 19]. The idea of electron-muon
collisions has been analyzed in past decades [20]-[25], and recently it has
been brought up again in a study of an extra $Z^{\prime}$ boson with generic
couplings that could be searched for through a Lepton Flavor Violating (LFV)
process like $e^{-}\mu^{+}\to e^{+}\mu^{-}$ [26]. In a subsequent article, the
construction of a high energy $e^{-}\mu^{+}$ collider has been proposed [27].
Besides the clear capability of probing muon LFV effects like the $He\mu$
coupling, this machine could even test the Higgs-bottom quark coupling [27].
One important advantage of such a machine is the very low level of SM
background processes as they come mostly from vector boson fusion. Other than
the elastic $e^{-}\mu^{+}\to e^{-}\mu^{+}$ there is no two-to-two process to
contend with and this is what motivates our study. We are interested in
fermion pair and two-photon production that can come from contact terms
$e^{-}\mu^{+}f\bar{f}$ and $e^{-}\mu^{+}\gamma\gamma$.

LFV operators are already being strongly constrained from low energy
measurements like muon decays and $\mu-e$ transitions [28, 29, 30]. We will
show the comparison of the potential limits estimated here with these
precision measurements. For $e^{-}\mu^{+}\to f\bar{f}$ we consider the four-
fermion dimension-6 operators of the Standard Model Effective Field Theory
(SMEFT) as given by the well known Warsaw basis [31]. We will also consider
some dimension 8 operators that have been constrained in [29]. Given the
chiral structure of the SMEFT, it is straightforward to obtain the
corresponding amplitudes in the helicity basis. For the two-photon production
$e^{-}\mu^{+}\to\gamma\gamma$ we address the contribution from the LFV dipole
operator, that is very strongly constrained by $\mu\to e\gamma$, and we
consider as well a dimension 8 contact operator that is bound more strictly by
our results than by low-energy measurements.

The paper is organized as follows. In section 2 we analyse the two-body
processes induced by four-fermion $e\mu ff$ operators. We explicity provide
the corresponding helicity amplitudes and cross sections. Then, we estimate
limits on the Wilson coefficients for 13 dimension 6 and 5 dimension 8 four-
fermion operators. In section 3 we consider $\gamma\gamma$ production
originating from two different contributions: a trivalent dimension 6 dipole
and a dimension 8 contact operator. We write down the sum of squared
amplitudes and the total cross sections. We estimate limits on the
coefficients in the same way we did in the previous section, but then we carry
out a detailed Monte Carlo analysis of signal and background processes. In
this way we show that our simplified strategy to obtain limits is realistic.
Finally, in section 4 we summarize our results.

##  2 Four fermion operators.

Since we are dealing with massless chiral fermions, the amplitudes are most
conveniently written in the helicity basis (see for example Ref. [32]).
Assuming some degree of polarization in the incoming beams, the general form
of the cross section can be divided in four terms:

$\displaystyle\sigma_{Pe^{-}P\mu^{+}}$ $\displaystyle=$
$\displaystyle\frac{1+\mathcal{P}_{e^{-}}}{2}\frac{1+\mathcal{P}_{\mu^{+}}}{2}\;\sigma_{++}\;+\;\frac{1-\mathcal{P}_{e^{-}}}{2}\frac{1-\mathcal{P}_{\mu^{+}}}{2}\;\sigma_{--}$
(1) $\displaystyle+$
$\displaystyle\frac{1+\mathcal{P}_{e^{-}}}{2}\frac{1-\mathcal{P}_{\mu^{+}}}{2}\;\sigma_{+-}\;+\;\frac{1-\mathcal{P}_{e^{-}}}{2}\frac{1+\mathcal{P}_{\mu^{+}}}{2}\;\sigma_{-+}\;.$

In this section we will address four fermion contact operators and we will
obtain sensitivity limits with the following method: first, for each operator
we will assume a total polarization in the cross section. For instance, in
$e_{R}^{-}\mu_{R}^{+}\to u_{R}{\bar{u}}_{R}$ that comes from the right-
chirality $Q_{eu}$ operator the ${\cal M}(+-+-)$ amplitude gives rise to the
$\sigma_{+-}$ cross section term, and we assume $\mathcal{P}_{e^{-}}=+1$,
$\mathcal{P}_{\mu^{+}}=-1$. Second, we neglect backgrounds and obtain
sensitivities in optimal conditions. Third, we will assume an integrated
luminosity of order $1$ ${\rm ab}^{-1}$ (as in Ref. [27]) and require a
minimum cross section of $0.04{\rm fb}$ from the operator contribution. Such a
cross section yields a significant amount of $40$ events. In this way, since
we are working with ideal conditions we will avoid being too optimistic. There
are $13$ dimension 6 four-fermion operators and 9 of them only contribute to
one of the four polarized cross sections. In their case, if there is only
partial beam polarization the actual number of events could be reduced by
about half or even by one quarter in case of unpolarized beams. On the other
hand, we do not expect the presence of background to reduce the sensitivity
significantly. In the next section we will corroborate that this is true for
$\gamma\gamma$ production; where we perform a full signal vs background
analysis.

###  2.1 Helicity amplitudes.

In this work we are interested in the process $e^{-}(p_{1})\mu^{+}(p_{2})\to
f(p_{3}){\bar{f}}(p_{4})$ that does not exist in the SM at tree level, but
that is generated by dimension six four-fermion operators in the SMEFT.
Hereafter, $f\bar{f}$ will stand for any of the three charged leptons $e$,
$\mu$ and $\tau$, or any of the quarks except the top quark. The list of four
fermion effective operators in Ref. [31] include arbitrary $(p,r,s,t)$ flavor
indices to take into account. However, this does not mean that there are so
many possible different helicity amplitudes as there are also Fierz identities
that relate them. For instance, the operator
$O_{ll}=\bar{l}_{p}\gamma_{\nu}l_{r}\bar{l}_{s}\gamma^{\nu}l_{t}$ gives rise
to the same amplitude $e^{-}\mu^{+}\to e^{-}e^{+}$ for any combination $2111$,
$1211$, $1121$ and $1112$. We choose to work with $2111$ that we denote as
$O^{2111}_{ll}=\bar{l}_{\mu}\gamma_{\nu}l_{e}\bar{l}_{e}\gamma^{\nu}l_{e}$.
For the processes $e^{-}\mu^{+}\to\mu^{-}\mu^{+}$ and
$e^{-}\mu^{+}\to\tau^{-}\tau^{+}$ the amplitudes from $O^{2122}_{ll}$ and
$O^{2133}_{ll}$ are clearly the same. In general, we shall take the flavor
indices as $2111$ for all the operators. This specific choice may exclude some
other non-equivalent combinations in some operators, but we have found that
this assignment, in the end, covers all the possible helicity amplitudes. The
purpose of this preliminary study, rather than being comprehensive, is to get
a first glimpse of the potential sensitivity of the $e^{-}\mu^{+}$ collider in
two-to-two processes generated by contact operators and compare with the
sensitivities of low energy measurements.

The helicity amplitudes are shown in Table 1. We see that
$(\bar{L}L)(\bar{L}L)$ operators $Q_{ll}$, $Q^{(1)}_{lq}$ and $Q^{(3)}_{lq}$
generate exactly the same helicity amplitudes as they all involve left-chiral
spinors. On the other hand, the $(\bar{L}L)(\bar{R}R)$ operator
$Q_{ll}=\bar{l}_{\mu}\gamma_{\mu}l_{e}\bar{e}_{e}\gamma^{\mu}e_{e}$, gives
rise to two possible combinations of final state chiralities:
$e^{-}_{R}e^{+}_{R}$ and $e^{-}_{R}e^{+}_{L}$. They are not equivalent as
shown in Table 1. Notice that many helicity amplitudes are actually equal up
to a phase factor related to azimuthal angles. The phase factor is important
only in the case of more than one diagram due to interference effects. We also
point out that as a preliminary analysis, we will be considering the
contribution of each operator separately. Because of the relations [32]:

$\displaystyle\langle ij\rangle=-\langle
ji\rangle\;,\;\;[ij]=-[ji]\;,\;\;\langle ij\rangle=[ji]^{*}\;,\;\;|\langle
ij\rangle|^{2}=2p_{i}\cdot p_{j}\;,$ (2)

there are really only three different amplitude structures in this study:
$[12][34]$, $[13][24]$ and $[14][23]$. Their squares are proportional to the
Mandelstam $s^{2}$, $t^{2}$ and $u^{2}$ respectively. So, they are actually
not independent. In fact, there is an identity that is easy to verify [32]:
$[12][34]=[14][32]+[13][24]$.

| ($\bar{L}L$)($\bar{L}L$) | $e^{-}(p_{1})\mu^{+}(p_{2})\to f(p_{3}){\bar{f}}(p_{4})$ | ${\cal M}$
---|---|---|---
$O_{ll}$ | $\bar{l}_{\mu}\gamma_{\nu}l_{e}\bar{l}_{e}\gamma^{\nu}l_{e}$ | $e^{-}_{L}\mu^{+}_{L}\to e^{-}_{L}e^{+}_{L}$ | $2[14]\langle 32\rangle$
$O_{lq}^{(1)}$ | $\bar{l}_{\mu}\gamma_{\nu}l_{e}\bar{q}\gamma^{\nu}q$ | $e^{-}_{L}\mu^{+}_{L}\to u_{L}{\bar{u}}_{L}$ | $2[14]\langle 32\rangle$
$O_{lq}^{(3)}$ | $\bar{l}_{\mu}\gamma_{\nu}\tau^{I}l_{e}\bar{q}\gamma^{\nu}\tau^{I}q$ |  $e^{-}_{L}\mu^{+}_{L}\to u_{L}{\bar{u}}_{L}$, $d_{L}{\bar{d}}_{L}$ | $\mp 2[14]\langle 32\rangle$
| ($\bar{R}R$)($\bar{R}R$) | $e^{-}(p_{1})\mu^{+}(p_{2})\to f(p_{3}){\bar{f}}(p_{4})$ | ${\cal M}$
$O_{ee}$ | $\bar{e}_{\mu}\gamma_{\nu}e_{e}\bar{e}_{e}\gamma^{\nu}e_{e}$ | $e^{-}_{R}\mu^{+}_{R}\to e^{-}_{R}e^{+}_{R}$ | $2\langle 14\rangle[32]$
$O_{eu}$ | $\bar{e}_{\mu}\gamma_{\nu}e_{e}\bar{u}\gamma^{\nu}u$ | $e^{-}_{R}\mu^{+}_{R}\to u_{R}{\bar{u}}_{R}$ | $2\langle 14\rangle[32]$
$O_{ed}$ | $\bar{e}_{\mu}\gamma_{\nu}e_{e}\bar{d}\gamma^{\nu}d$ | $e^{-}_{R}\mu^{+}_{R}\to d_{R}{\bar{d}}_{R}$ | $2\langle 14\rangle[32]$
| ($\bar{L}L$)($\bar{R}R$) | $e^{-}(p_{1})\mu^{+}(p_{2})\to f(p_{3}){\bar{f}}(p_{4})$ | ${\cal M}$
$O_{le}$ | $\bar{l}_{\mu}\gamma_{\nu}l_{e}\bar{e}_{e}\gamma^{\nu}e_{e}$ | $e^{-}_{L}\mu^{+}_{L}\to e^{-}_{R}e^{+}_{R}$ | $2\langle 24\rangle[31]$
$O_{le}$ | $\bar{l}_{\mu}\gamma_{\nu}l_{e}\bar{e}_{e}\gamma^{\nu}e_{e}$ | $e^{-}_{R}\mu^{+}_{L}\to e^{-}_{R}e^{+}_{L}$ | $2\langle 21\rangle[34]$
$O_{lu}$ | $\bar{l}_{\mu}\gamma_{\nu}l_{e}\bar{u}\gamma^{\nu}u$ | $e^{-}_{L}\mu^{+}_{L}\to u_{R}{\bar{u}}_{R}$ | $2\langle 24\rangle[31]$
$O_{ld}$ | $\bar{l}_{\mu}\gamma_{\nu}l_{e}\bar{d}\gamma^{\nu}d$ | $e^{-}_{L}\mu^{+}_{L}\to d_{R}{\bar{d}}_{R}$ | $2\langle 24\rangle[31]$
$O_{qe}$ | $\bar{q}\gamma_{\nu}q\bar{e}_{\mu}\gamma^{\nu}e_{e}$ | $e^{-}_{R}\mu^{+}_{R}\to u_{L}{\bar{u}}_{L},d_{L}{\bar{d}}_{L}$ | $2[24]\langle 31\rangle$
| ($\bar{L}R$)($\bar{R}L$) | $e^{-}(p_{1})\mu^{+}(p_{2})\to f(p_{3}){\bar{f}}(p_{4})$ | ${\cal M}$
$O_{ledq}$ | $\bar{l}_{\mu}e_{e}\bar{d}q$ | $e^{-}_{R}\mu^{+}_{L}\to d_{R}{\bar{d}}_{L}$ | $[12][43]$
| ($\bar{L}R$)($\bar{L}R$) | $e^{-}(p_{1})\mu^{+}(p_{2})\to f(p_{3}){\bar{f}}(p_{4})$ | ${\cal M}$
$O_{lequ}^{(1)}$ | $\bar{l}^{j}_{\mu}e_{e}\epsilon_{jk}\bar{q}^{k}u$ | $e^{-}_{R}\mu^{+}_{L}\to u_{R}{\bar{u}}_{L}$ | $\langle 12\rangle[34]$
$O_{lequ}^{(3)}$ | $\bar{l}^{j}_{\mu}\sigma_{\mu\nu}e_{e}\epsilon_{jk}\bar{q}^{k}\sigma^{\mu\nu}u$ | $e^{-}_{R}\mu^{+}_{L}\to u_{R}{\bar{u}}_{L}$ | $4([13][24]+[14][23])$
Table 1: Dimension 6 operators with common flavor indices $(2111)$ and the
helicity amplitudes generated. The $\frac{C_{O}}{\Lambda^{2}}$ factor to be
included. The $R,L$ indices refer to chiralities. Final state $e^{-}e^{+}$
also stands for $\mu^{-}\mu^{+}$ and $\tau^{-}\tau^{+}$. Similarly, $d\bar{d}$
also stands for $s\bar{s}$ and $b\bar{b}$; and $u\bar{u}$ also stands for
$c\bar{c}$. Top quarks excluded.

Now let us turn our attention to dimension 8 operators. There are multiple
structures but we will pay attention to those specific operators that have
been bounded from $\mu\to e$ processes [28, 29]. Moreover, there are dimension
8 operators that coincide with some dimension 6 operator except for an
additional $H^{\dagger}H$ term. Obviously, their amplitudes would be equal
except for some rescaling factor. The amplitudes of the dimension 8 operators
that give rise to chiral structures that do not appear at dimension 6 are
shown in Table 2.

| ($\bar{L}R$)($\bar{L}R$) | $e^{-}(p_{1})\mu^{+}(p_{2})\to f(p_{3}){\bar{f}}(p_{4})$ | ${\cal M}$
---|---|---|---
$O_{le}^{(8)}$ | $\bar{l}_{\mu}He_{e}\bar{l}_{e}He_{e}$ | $e^{-}_{R}\mu^{+}_{L}\to e^{-}_{R}e^{+}_{L}$ | $\langle 12\rangle[43]$
$O_{Tle}^{(8)}$ | $\bar{l}_{\mu}H\sigma^{\mu\nu}e_{e}\bar{l}_{e}H\sigma_{\mu\nu}e_{e}$ | $e^{-}_{R}\mu^{+}_{L}\to e^{-}_{R}e^{+}_{L}$ | $4([31][24]+[41][23])$
$O_{leqd1}^{(8)}$ | $\bar{l}_{\mu}He_{e}\bar{q}Hd$ | $e^{-}_{R}\mu^{+}_{L}\to d_{R}{\bar{d}}_{L}$ | $\langle 12\rangle[43]$
$O_{leqd3}^{(8)}$ | $\bar{l}_{\mu}H\sigma^{\mu\nu}e_{e}\bar{q}H\sigma_{\mu\nu}d$ | $e^{-}_{R}\mu^{+}_{L}\to d_{R}{\bar{d}}_{L}$ | $4([31][24]+[41][23])$
$O_{le}^{(8)}$ | $\bar{l}_{\mu}He_{e}\bar{l}_{e}He_{e}$ | $e^{-}_{R}\mu^{+}_{L}\to e^{-}_{L}e^{+}_{R}$ | $\langle 24\rangle[31]$
$O_{Tle}^{(8)}$ | $\bar{l}_{\mu}H\sigma^{\mu\nu}e_{e}\bar{l}_{e}H\sigma_{\mu\nu}e_{e}$ | $e^{-}_{R}\mu^{+}_{L}\to e^{-}_{L}e^{+}_{R}$ | $4([21][34]+[14][23])$
| ($\bar{L}R$)($\bar{R}L$) | $e^{-}(p_{1})\mu^{+}(p_{2})\to f(p_{3}){\bar{f}}(p_{4})$ | ${\cal M}$
$O_{leuq}^{(8)}$ | $\bar{l}_{\mu}He_{e}\bar{u}{\tilde{H}}^{\dagger}q$ | $e^{-}_{R}\mu^{+}_{L}\to u_{R}{\bar{u}}_{L}$ | $\langle 12\rangle\langle 43\rangle$
Table 2: Helicity amplitudes generated by dimension 8 operators with common
flavor indices $(2111)$ associated to four-fermion contact vertices that do
not appear with dimension 6 operators. A
$\frac{v^{2}}{2}\frac{C_{O}}{\Lambda^{4}}$ factor is to be included. As in
Table (1), final state $e^{-}e^{+}$ also stands for $\mu^{-}\mu^{+}$ and
$\tau^{-}\tau^{+}$; $d\bar{d}$ also stands for $s\bar{s}$ and $b\bar{b}$; and
$u\bar{u}$ also stands for $c\bar{c}$.

###  2.2 Cross sections and limits on coefficients.

As mentioned above, there are only three different types of helicity
amplitudes squared. Each of them gives rise to one specific expression for the
cross section. For $|[12][34]|^{2}=s^{2}$, after dividing by the energy scale
$\Lambda^{4}$ and integrating over the phase space we obtain:
$\sigma_{1234}=s/(16\pi\Lambda^{4})$. For the center-of-mass (CM) collider
energy $\sqrt{s}=1.095$ TeV ($E_{e}=100,\;E_{\mu}=3000$GeV) we obtain
$\sigma_{1234}(\sqrt{s}=1.095{\rm TeV})=36.3$fb. We will show the
corresponding cross sections for each operator in terms of this common
$\sigma_{1234}$. The other amplitudes yield
$\sigma_{1324}=\sigma_{1423}=\sigma_{1234}/3$. The amplitudes obtained here
grow with the collision energy, but our effective theory calculation is valid
for energies below the cut-off scale $\Lambda=4$TeV.

Three benchmark collision energies are proposed in [27]: (1)
$E_{e}=20,\;E_{\mu}=200$GeV; (2) $E_{e}=50,\;E_{\mu}=1000$GeV; and (3)
$E_{e}=100,\;E_{\mu}=3000$GeV that correspond to CM energies of
$\sqrt{s}=2\sqrt{E_{e}E_{\mu}}=\;0.126,\;0.447,\;1.095$TeV, respectively.
Since all cross sections are proportional to $s|C_{O}|^{2}$ computing the
bound at an energy $\sqrt{s_{2}}$ assuming we know the bound at an energy
$\sqrt{s_{1}}$ is straigthforward: we just multiply by the ratio
$\sqrt{s_{1}}/\sqrt{s_{2}}$. For instance, $0.447/1.095=0.4$ and so we see
that the limits at benchmark (3) will be more than twice stronger than
benchmark (2). We shall focus on the $E_{e}=100,\;E_{\mu}=3000$GeV benchmark
in this study. However, in section 3 we will see that the cross sections for
$e^{-}\mu^{+}\to\gamma\gamma$ are not proportional to $s$, but one is constant
in energy and the other is proportional to $s^{2}$. We will provide limits
obtained with the benchmark (3) $E_{e}=100,\;E_{\mu}=3000$GeV and for greater
energies.

The dimension 6 operators of Table 1 yield the following cross sections:

$\displaystyle\frac{\sigma_{++}}{\sigma_{1234}}$ $\displaystyle=$
$\displaystyle
4|C_{le}|^{2}+N_{c}|C_{ledq}|^{2}+N_{c}|C_{lequ}^{(1)}|^{2}+\frac{16}{3}N_{c}|C_{lequ}^{(3)}|^{2}$
$\displaystyle\frac{\sigma_{+-}}{\sigma_{1234}}$ $\displaystyle=$
$\displaystyle\frac{4}{3}|C_{ee}|^{2}+\frac{4}{3}N_{c}|C_{eu}|^{2}+\frac{4}{3}N_{c}|C_{ed}|^{2}+\frac{4}{3}N_{c}|C_{qe}|^{2}$
(3) $\displaystyle\frac{\sigma_{-+}}{\sigma_{1234}}$ $\displaystyle=$
$\displaystyle\frac{4}{3}|C_{ll}|^{2}+\frac{4}{3}|C_{le}|^{2}+\frac{4}{3}N_{c}|C_{lq}^{(1)}|^{2}+\frac{4}{3}N_{c}|C_{lq}^{(3)}|^{2}+\frac{4}{3}N_{c}|C_{lu}|^{2}+\frac{4}{3}N_{c}|C_{ld}|^{2}\;.$

Where the $\sigma_{--}$ term does not appear for the flavor assignment $2111$,
but the operators that generate $\sigma_{++}$ would also generate
$\sigma_{--}$ with the assignment $1211$.

By requiring that the value of a $C_{Q}$ coefficient be enough to yield the
minimum $0.04$fb of production cross section, we obtain the following lower
limits for $E_{e}=100,\;E_{\mu}=3000$GeV:

$\displaystyle C_{ll}\;,\;C_{ee}\;$ $\displaystyle\geq$
$\displaystyle\;2.88\times 10^{-2}$ $\displaystyle
C_{ledq}\;,C^{(1)}_{lequ}\;$ $\displaystyle\geq$ $\displaystyle\;1.92\times
10^{-2}$ $\displaystyle C_{le}\;,\;C_{lu}\;,C_{ld}\;,\;C_{qe}\;$
$\displaystyle\geq$ $\displaystyle\;1.66\times 10^{-2}$ (4) $\displaystyle
C_{eu}\;,\;C_{ed}\;,C^{(1)}_{lq}\;,\;C^{(3)}_{lq}\;$ $\displaystyle\geq$
$\displaystyle\;1.66\times 10^{-2}$ $\displaystyle C^{(3)}_{lequ}\;$
$\displaystyle\geq$ $\displaystyle\;0.83\times 10^{-2}$

We can compare with the limits from low energy processes $\mu\to e\gamma$,
$\mu\to e\bar{e}e$ and $\mu\to e$ conversion in nuclei $\mu A\to eA$ as
recently reported in Tables 6 and 7 of Ref. [28]. We would like to point out
that in almost all cases the most stringent bounds in those tables come from
$\mu A\to eA$ conversion in nuclei [33]. The exceptions are first family four-
lepton coefficients $C_{ll},C_{le},C_{ee}$ that are constrained by $\mu\to
e{\bar{e}}e$ (Table 6[28]) and $C^{(3)e\mu cc}_{lequ}$ that is constrained by
$\mu\to e\gamma$ (Table 7[28]). Moreover, these $\mu A\to eA$ bounds are in
fact around four orders of magnitude stronger than the ones from $\mu\to
e\gamma$ and $\mu\to e\bar{e}e$. This brings up one important observation:
that if we ignored $\mu A\to eA$ the conclusion would be that the $e\mu$
collider would yield much stronger constraints than any of the low energy
measurements. In Table 3 we show the limits to each operator coefficient (at
the scale $m_{W}$) and the ratio $C^{\rm Low\;energy}_{Q}/C^{\rm
Collider}_{Q}$ for each possible final state. Not surprisingly, we can observe
that for the first family fermions the limits are very stringent. However, for
the second and third family modes that get their low energy limits via loop
contributions, the bounds are weaker and about the same order of magnitude as
the potential limits from the $e^{-}\mu^{+}\to f\bar{f}$ production. Strictly
speaking, the collider sensitivity for each coefficient is taken at a scale of
order 1 TeV; about ten times higher than the electroweak scale. In a more
precise analysis one should take into account the renormalization group
dependence on the scale. However, one should bear in mind that such
corrections are usually of order a few percent. For instance, the coefficient
$C_{DL}$ associated to the operator:
$(C_{DL}m_{\mu}/m_{t}^{2}){\bar{e}}_{R}\sigma^{\mu\nu}\mu_{L}F^{\mu\nu}$. An
upper limit from its contribution to $\mu\to e\gamma$ decay is reported as
$|C_{DL}|<1.05\times 10^{-8}$ at the muon mass scale [28]. The same limit
becomes $|C_{DL}|<1.12\times 10^{-8}$ at the $m_{W}$ scale, which is only a
$7\%$ numerical variation. For another example, let us consider the dimension
6 coefficient $C_{el}$ above. In a series of articles [34, 35, 36, 37] we can
find the renormalization group evolution of the dimension 6 SMEFT operators
listed in the Warsaw basis. The running of $C_{el}$ in general depends on
several other operators, but let us take the term proportional to $C_{el}$
itself. Then, we find that $C_{le}(1TeV)\sim 0.96C_{le}(m_{W})$, that is only
a $4\%$ correction. As we are only interested in a preliminary assessment of
sensitivities and comparisons with low energy experiment limits we shall not
take corrections of scale dependence into account. Recently, limits on the
$C^{(1)}_{lq}$ and $C^{(3)}_{lq}$ coefficients have been published based on
the LHC $pp\to e\mu$ dilepton production that are approximately
$C^{(1,3)}_{lq}<0.3$ for first family quarks, $C^{(1,3)}_{lq}<2.0$ for second
family and $C^{(1,3)}_{lq}<5.0$ for $b\bar{b}$ [38]. They project that with a
hundred times more luminosity the HL-LHC could reduce these limits to one
third of the current value. From what we have found here, the $e\mu$ collider
would have at least one order of magnitude greater sensitivity than the HL-LHC
for these operators.

$C_{Q}$ |  $C^{\rm Low}_{Q}$ ($C^{\rm Low}/C^{\rm Coll}$) |  $C^{\rm Low}_{Q}$ ($C^{\rm Low}/C^{\rm Coll}$) |  $C^{\rm Low}_{Q}$ ($C^{\rm Low}/C^{\rm Coll}$)
---|---|---|---
| $e^{-}e^{+}$ | $\mu^{-}\mu^{+}$ | $\tau^{-}\tau^{+}$
$C_{ll}$ |  $4.16\times 10^{-4}$ ($1.45\times 10^{-2}$) |  $0.98\times 10^{-2}$ ($0.34$) |  $1.97\times 10^{-2}$ ($0.69$)
$C_{ee}$ |  $4.16\times 10^{-4}$ ($1.45\times 10^{-2}$) |  $0.98\times 10^{-2}$ ($0.34$) |  $1.97\times 10^{-2}$ ($0.69$)
$C_{le}$ |  $4.92\times 10^{-4}$ ($2.96\times 10^{-2}$) |  $1.99\times 10^{-2}$ ($1.20$) |  $1.98\times 10^{-2}$ ($1.19$)
| $d\bar{d}$ | $s\bar{s}$ | $b\bar{b}$
$C^{(1)}_{lq}$ |  $1.51\times 10^{-5}$ ($0.91\times 10^{-3}$) |  $1.60\times 10^{-2}$ ($0.96$) |  $2.49\times 10^{-2}$ ($1.50$)
$C^{(3)}_{lq}$ |  $2.69\times 10^{-4}$ ($1.62\times 10^{-2}$) |  $6.08\times 10^{-3}$ ($0.37$) |  $2.49\times 10^{-2}$ ($1.50$)
$C_{ld}$ |  $2.80\times 10^{-5}$ ($1.69\times 10^{-3}$) |  $1.97\times 10^{-2}$ ($1.19$) |  $2.49\times 10^{-2}$ ($1.50$)
$C_{ed}$ |  $2.86\times 10^{-5}$ ($1.72\times 10^{-3}$) |  $1.97\times 10^{-2}$ ($1.19$) |  $2.49\times 10^{-2}$ ($1.50$)
$C_{qe}$ |  $1.52\times 10^{-5}$ ($0.92\times 10^{-3}$) |  $1.58\times 10^{-2}$ ($0.95$) |  $2.50\times 10^{-2}$ ($1.51$)
$C_{ledq}$ |  $5.34\times 10^{-6}$ ($2.79\times 10^{-4}$) |  $1.11\times 10^{-4}$ ($5.80\times 10^{-3}$) |  $3.66\times 10^{-3}$ ($0.19$)
| $u\bar{u}$ | $c\bar{c}$ | $--$
$C_{lu}$ |  $3.30\times 10^{-5}$ ($1.99\times 10^{-3}$) |  $0.88\times 10^{-2}$ ($0.53$) |
$C_{eu}$ |  $3.19\times 10^{-5}$ ($1.92\times 10^{-3}$) |  $0.89\times 10^{-2}$ ($0.54$) |
$C^{(1)}_{lequ}$ |  $5.45\times 10^{-6}$ ($2.84\times 10^{-4}$) |  $0.97\times 10^{-3}$ ($5.1\times 10^{-2}$) |
$C^{(3)}_{lequ}$ |  $5.10\times 10^{-5}$ ($6.1\times 10^{-3}$) |  $5.34\times 10^{-6}$ ($6.4\times 10^{-4}$) |
Table 3: Upper limits to the dimension 6 coefficients from the low energy
experiments as reported in tables 6 and 7 of [28]. They have been multiplied
by a factor $(4/0.174)^{2}$ to account for the normalization scale $m_{t}$
used by [28] instead of the $\Lambda=4$TeV scale adopted in this study. In
addition, the ratio of limits Low/Collider is shown in parenthesis
($E_{e}=100,\;E_{\mu}=3000$GeV).

The dimension 8 operators of Table 2 yield the following cross sections:

$\displaystyle\frac{\sigma_{++}}{\sigma_{1234}}(e^{-}_{R}e^{+}_{L})$
$\displaystyle=$
$\displaystyle\frac{v^{4}}{4\Lambda^{4}}\left(|C_{le}^{(8)}|^{2}+N_{c}|C_{leqd1}^{(8)}|^{2}+N_{c}|C_{leuq}^{(8)}|^{2}+\frac{16}{3}N_{c}|C_{leqd3}^{(8)}|^{2}+\frac{16}{3}|C_{Tle}^{(8)}|^{2}\right)$
$\displaystyle\frac{\sigma_{++}}{\sigma_{1234}}(e^{-}_{L}e^{+}_{R})$
$\displaystyle=$
$\displaystyle\frac{v^{4}}{4\Lambda^{4}}\left(\frac{1}{3}|C_{le}^{(8)}|^{2}+\frac{112}{3}N_{c}|C_{Tle}^{(8)}|^{2}\right)$
(5)

A comment on the tensor operators $Q^{(3)}_{lequ}$, $Q_{leqd3}$ and $Q_{Tle}$
is in order here. $Q^{(3)}_{lequ}$, $Q_{leqd3}$ and the
$\sigma_{++}(e^{-}_{R}e^{+}_{L})$ contribution of $Q_{Tle}$ generate the same
amplitude $[13][24]+[14][23]$. Except for the $v^{4}/(4\Lambda^{4})$
additional factor, the cross sections are equal. However, for the
$\sigma_{+-}(e^{-}_{L}e^{+}_{R})$ helicity mode the $Q_{Tle}$ amplitude is
$[21][34]+[14][23]$. In the CM frame the functional dependence on the $p_{3}$
polar angle’s $\cos{\theta}$ is: $[13][24]+[14][23]\sim\cos{\theta}$. If we
compare with $[21][34]+[14][23]\sim(3+\cos{\theta})/2$ we can see that the
latter yields a much greater cross section.

$C_{Q}$ | $C^{\rm Low}/C^{\rm Coll}$
---|---
$C_{le}^{(8)}$ | $5\times 10^{-3}\;/14.0$
$C_{Tle}^{(8)}$ | $2\times 10^{-3}\;/1.3$
$C_{leqd1}^{(8)}$ | $3\times 10^{-5}\;/8.1$
$C_{leqd3}^{(8)}$ | $7\times 10^{-4}\;/3.5$
$C_{leuq}^{(8)}$ | $3\times 10^{-5}\;/8.1$
Table 4: The ratio of potential limits from low energy processes and minimum
observable values at the collider for the dimension 8 coefficients (Table 5 of
[29]; $E_{e}=100,\;E_{\mu}=3000$GeV).

As in Table 3, in Table 4 we show the ratio $C^{\rm Low}/C^{\rm Coll}$ for the
dimension 8 operators. Comparing with the dimension 6 coefficients, there is a
suppresing $v^{2}/(2\Lambda^{2})$ factor of order $2\times 10^{-3}$ and the
$C^{\rm Coll}$ minimum values have to be much bigger. In contrast, according
to Table 5 of Ref. [29] the limits from low energy experiments are still very
stringent for dimension 8 couplings.

What we have learned from Table 3 is that for first family fermions $f\bar{f}$
the sensitivity of the low energy measurement of $\mu A\to eA$ conversion in
nuclei is two or more orders of magnitude higher. Maybe all the operator
coefficients are indeed very suppressed; regardless of potential
cancellations. However, for most of the second and third family $f\bar{f}$
states the collider sensitivities are of the same order of magnitude as the
low energy limits. The electron-muon collider should be able to provide
additional and competitive limits to constrain the set of dimension 6 four
fermion LFV operators. The same may not be true for dimension 8 operators, at
least for the ones that can be constrained by the low energy experiments.

##  3 The $e^{-}\mu^{+}\to\gamma\gamma$ process.

The $e^{-}\mu^{+}\to\gamma\gamma$ amplitude can be generated by the dimension
6 flavor-changing magnetic dipole operator

$\displaystyle Q_{eA}=\bar{l}_{\mu}\sigma^{\mu\nu}e_{e}HF_{\mu\nu}\;,$ (6)

through $t$ and $u$ channel diagrams where one of the photons is emitted by
the effective coupling. Another possibility comes from the dimension 8
operator

$\displaystyle Q_{eAA}=\bar{l}_{\mu}e_{e}HF^{\mu\nu}F_{\mu\nu}\;.$ (7)

We have two chiral versions: $Q_{eAL}$ and $Q_{eAR}$ ($Q_{eAAL}$, $Q_{eAAR}$)
referring to left-handed and right-handed electron respectively. In figure 1
we show the Feynman diagrams associated to each operator.

Figure 1: The anomalous two-to-two process $e^{-}\mu^{+}\to\gamma\gamma$.
Process (a) (t-channel) as given by the trivalent vertex dipole operator
$Q_{eA}$. Process (b) as given by the contact term, dimension 8 operator
$Q_{eAA}$.

Both operators give rise to amplitudes that do not depend on angles:

$\displaystyle\sum|{\cal M}|^{2}$ $\displaystyle=$
$\displaystyle\left(|C_{eAL}|^{2}+|C_{eAR}|^{2}\right)e^{2}\frac{v^{2}}{\Lambda^{4}}48s$
(8) $\displaystyle+$
$\displaystyle\left(|C_{eAAL}|^{2}+|C_{eAAR}|^{2}\right)\frac{v^{2}}{\Lambda^{8}}4s^{3}$

They yield the total cross sections:

$\displaystyle\sigma_{--}+\sigma_{++}$ $\displaystyle=$
$\displaystyle\left(|C_{eAL}|^{2}+|C_{eAR}|^{2}\right)\frac{3e^{2}v^{2}}{2\pi\Lambda^{4}}$
(9) $\displaystyle=$
$\displaystyle\left(|C_{eAL}|^{2}+|C_{eAR}|^{2}\right)4.32{\rm fb}\;.$
$\displaystyle\sigma_{--}+\sigma_{++}$ $\displaystyle=$
$\displaystyle\left(|C_{eAAL}|^{2}+|C_{eAAR}|^{2}\right)\frac{v^{2}}{8\pi\Lambda^{8}}s^{2}$
$\displaystyle=$
$\displaystyle\left(|C_{eAL}|^{2}+|C_{eAR}|^{2}\right)0.0206{\rm fb}\;,$

where the numerical value on the second line is independent of the collision
energy, but the numerical value on the fourth line is taken at
$\sqrt{s}=1.095$ TeV.

We see, then, that coefficients $C_{eAL(R)}$ and $C_{eAAL(R)}$ of order $0.1$
and $\sqrt{2}$ respectively would give us a $\sigma=0.04$fb value that is our
minimum acceptable cross section. Let us notice that the current limit from
$\mu\to e\gamma$ is of order $5\times 10^{-6}$ for the dipole coefficient
$C_{eAL(R)}$, well below the potential collider sensitivity. On the other
hand, for the dimension 8 operator coefficient $C_{eAAL(R)}$ the potential
sensitivity from $\mu A\to eA$ transitions is six orders of magnitude less
stringent: $C_{eAAL(R)}\leq 3.2$ [29]. We point out that, in this study, this
last coefficient is the one with the highest sensitivity by the $e\mu$
collider as compared to the low energy measurements.

Seeing the relatively high sensitivity to the electromagnetic $Q_{eAA}$
operator one may wonder what about the effective vertex $e\mu GG$? We can use
the $Q_{eAA}$ calculation with the analogous gluon operator

$\displaystyle
Q_{eGG}=\frac{C_{eGG}}{\Lambda^{4}}\bar{l}_{\mu}e_{e}HG^{a\mu\nu}G^{a}_{\mu\nu}\;.$
(10)

The $e^{-}\mu^{+}\to gg$ production the cross section is now $N_{c}=8$ times
greater than $e^{-}\mu^{+}\to\gamma\gamma$ and a coefficient of order $0.5$
would yield the minimum observable production. However, according to [29] the
low energy $\mu A\to eA$ limit is three orders of magnitude more stringent in
this case: $C_{eGG}\leq 1.6\times 10^{-3}$.

As mentioned above, the amplitudes squared for $e^{-}\mu^{+}\to\gamma\gamma$
do not depend on the polar angle. That means that in terms of rapidity, for
instance $y=y^{*}_{3}$ in the CM frame:

$\displaystyle\frac{d\sigma}{dy}=\frac{d\sigma}{dc_{\theta}}\frac{dc_{\theta}}{dy}=a_{0}\frac{4\exp(2y)}{(1+\exp(2y))^{2}}\;,$
(11)

where $a_{0}$ is a constant. The shape of the rapidity distribution in the CM
frame is then centered around zero with a width of approximately $2$ units as
shown below in figure 3, right column. In the lab frame the center is shifted
towards $-1.70$ (see next subsection).

###  3.1 Monte Carlo analysis of $\gamma\gamma$ production and its SM
background.

In the SM, $\gamma\gamma$ production is given by the two-to-four process
$e^{-}\mu^{+}\to\gamma\gamma\nu_{e}\bar{\nu_{\mu}}$, involving 13 Feynman
diagrams in unitary gauge. There is also $e^{-}e^{+}$ production
$e^{-}\mu^{+}\to\nu_{e}\bar{\nu_{\mu}}e^{-}e^{+}$, involving 24 Feynman
diagrams. In figure 2 we show two representative diagrams for the
$\gamma\gamma$ and resonant $e^{-}e^{+}$ production in the SM.

Figure 2: Representative unitary-gauge Feynman diagrams for the SM processes
(a) $e^{-}\mu^{+}\to\gamma\gamma\nu_{e}\bar{\nu}_{\mu}$ and (b)
$e^{-}\mu^{+}\to e^{-}e^{+}\nu_{e}\bar{\nu}_{\mu}$.

Signal and background have very different kinematics and this makes the
separation straightforward. Our goal is to show how a basic set of cuts can
reduce the potential background dramatically. We point out here that in the
lab frame the rapidities $y$ are shifted with respect to those $y^{*}$ in the
CM frame:

$\displaystyle y=y^{*}-y_{0}\;,\;\;{\rm
with}\;\;y_{0}=\frac{1}{2}\ln\left(\frac{E_{\mu}}{E_{e}}\right)=1.70\;,$ (12)

where $y_{0}=1.7$ is the shift value for $E_{\mu}/E_{e}=30$. The muon beam
goes in the direction of $-\hat{k}$ and so the event products usually appear
on the backwards hemisphere.

Thus, to study $\gamma\gamma$ production in $e^{-}\mu^{+}$ collisions we
consider the three processes,

$\displaystyle e^{-}\mu^{+}$ $\displaystyle\rightarrow$
$\displaystyle\gamma\gamma,$ (13) $\displaystyle e^{-}\mu^{+}$
$\displaystyle\rightarrow$ $\displaystyle\gamma\gamma\bar{\nu}_{\mu}\nu_{e},$
(14) $\displaystyle e^{-}\mu^{+}$ $\displaystyle\rightarrow$ $\displaystyle
e^{-}e^{+}\bar{\nu}_{\mu}\nu_{e},$ (15)

where the first one is our signal process as depicted in figure 1 and (14),
(15) the SM backgrounds shown in figure 2. We implemented the effective
interactions (6), (7) in our Monte Carlo simulations by means of feynrules 2.0
[39]. We simulated the signal and background processes with madgraph 2.6 [40],
with beam energies $(E_{e},E_{\mu})=$ (100, 3000), (150, 4500) and (200, 6000)
GeV. We then have $\sqrt{s}=1.095$, 1.643 and 2.191 TeV, respectively. In the
two-body signal process (13), the CM energies of the final photons are fixed
at $E^{*}_{3}=E^{*}_{4}=\sqrt{s}/2$. This is unlike what happens with the
backgrounds (14), (15), where there is a continuous spectrum for $E^{*}_{3}$
and $E^{*}_{4}$. Similarly, for the signal process
$|\vec{p}_{3T}+\vec{p}_{4T}|=0$, but for the backgrounds
$|\vec{p}_{3T}+\vec{p}_{4T}|=\not\\!\\!E_{T}$ which has a continuous range of
values. Furthermore, we observe the final-state photons in the signal process
to be very central in the CM frame, with
$|y^{*}_{3,4}|=|y_{3,4}+y_{0}|\lesssim 2.5$ for the vast majority of events,
as expected in view of the analytical distribution (11) and as shown below in
figure 3. For the process (14), without restrictions on $E^{*}_{3,4}$, the
photon rapidity distribution in the lab frame is symmetric about
$-y_{0}=-1.70$, but very broad. If we require $E^{*}_{3,4}$ in (14) to be
large, however, the final-state photons must be very forward or backward. For
example, if $E^{*}_{3,4}\simeq 250$GeV the rapidities will have maxima at
$-y_{0}\pm\Delta y$ with $\Delta y\simeq 2$ as shown in figure 3; and $\Delta
y$ gets larger for greater values of $E^{*}_{3,4}$. For the process (15) the
electron rapidity distribution presents essentially the same features, but is
less forward-backward symmetric, as seen in the figure.

We are, thus, led to consider the following set of phase-space cuts,

$\begin{gathered}C_{0}:\;p_{3T},p_{4T}>1.0\,\mathrm{GeV},\qquad
C_{1}:\;E_{3}^{*},E_{4}^{*}>500.0\,\mathrm{GeV},\\\
C_{2}:\;p_{T}^{\mathrm{tot}}=|\vec{p}_{3T}+\vec{p}_{4T}|<20.0\,\mathrm{GeV},\qquad
C_{3}:\;|y_{3}+y_{0}|,|y_{4}+y_{0}|<1.75.\end{gathered}$ (16)

The cut $C_{0}$ is necessary to control infrared divergences in (14), (15). As
discussed above, the cuts $C_{1,2,3}$ in (16) have only small effects on the
signal cross section, but they do substantially decrease the cross section for
the backgrounds. The effect of the cuts (16) on the cross sections for the
processes (13), (14), (15) is illustrated at $(E_{e},E_{\mu})=$ (100, 3000)
GeV in table 5.

$\mathcal{P}_{\mu}$ | $\mathcal{P}_{e}$ | $\mathrm{cuts}$ | $\sigma_{\gamma\gamma}^{(3)}[\mathrm{fb}]$ | $\sigma_{\gamma\gamma}^{(4)}[\mathrm{fb}]$ | $\sigma_{\nu\nu\gamma\gamma}[\mathrm{fb}]$ | $\sigma_{\nu\nu ee}[\mathrm{fb}]$
---|---|---|---|---|---|---
$0.0$ | $0.0$ | $C_{0}$ | $2.16$ | $0.0103$ | $770.2$ | $426.6$
$0.0$ | $0.0$ | $C_{0,1}$ | $2.16$ | $0.0103$ | $0.0049$ | $0.073$
$0.0$ | $0.0$ | $C_{0,1,2}$ | $2.16$ | $0.0103$ | $0.0015$ | $0.0081$
$0.0$ | $0.0$ | $C_{0\mbox{--}3}$ | $2.04$ | $0.0097$ | $3.23\times 10^{-5}$ | $0.00108$
$+0.4$ | $+0.8$ | $C_{0,1}$ | $2.86$ | $0.0136$ | $0.0014$ | $0.081$
$+0.4$ | $+0.8$ | $C_{0,1,2}$ | $2.69$ | $0.0128$ | $0.00043$ | $0.0082$
$+0.4$ | $+0.8$ | $C_{0\mbox{--}3}$ | $2.69$ | $0.0128$ | $9.85\times 10^{-6}$ | $0.0038$
Table 5: Cumulative effects of the cuts (16) on the cross sections for (13),
(14), (15) at $E_{e}=$ 100 GeV, $E_{\mu}=$ 3 TeV.

The cross section $\sigma_{\gamma\gamma}^{(3)}$ refers to the anomalous
process (13) with only the trivalent $e\mu\gamma$ vertex and the Wilson
coefficients $C_{eAL}=C_{eAR}=1$. Similarly, $\sigma_{\gamma\gamma}^{(4)}$
refers to (13) through the $e\mu\gamma\gamma$ vertex and the coefficients
$C_{eAAL}=C_{eAAR}=1$. The numerical results agree with (9).

We expect the results for cross sections with cuts in table 5 to be quite
realistic, although detector efficiencies and acceptances have not been
allowed for in those results. However, we expect the rapidity acceptance
effects to be taken into account by the cut $C_{3}$ in (16), and we also
expect the efficiency for photon identification to be no less than 90%, so
that detector effects should be modest. The important exception to this,
however, is the background process (15), which in table 5 seems to represent
one-third of $\sigma_{\gamma\gamma}^{(4)}$, but which must actually be
adjusted for the electron-photon misidentification probability. In order to
settle this issue, we carried out a detector simulation using Delphes 3.4
[41].

An $e^{-}\mu^{+}$ collider is highly asymmetric, so we assume its detector to
have a correspondingly asymmetric design. We obtain a simple but effective
asymmetric detector simulation in the lab frame by assuming that in the CM
frame, in which the $e^{-}\mu^{+}$ collisions are forward-backward symmetric
on average, the detector possesses the same capabilities as the symmetric Muon
Collider detector implemented in Delphes 3.4111Which corresponds to the
configuration file cards/delphes_card_MuonColliderDet.tcl in the Delphes
distribution.. We simulate the collisions with madgraph 5 with a loosened
version of the cuts (16) in order to adequately populate the phase space,

$\begin{gathered}C^{\prime}_{0}=C_{0},\qquad
C^{\prime}_{1}:\;E_{3}^{*},E_{4}^{*}>250.0\,\mathrm{GeV},\\\
C^{\prime}_{2}:\;p_{T}^{\mathrm{tot}}=|\vec{p}_{3T}+\vec{p}_{4T}|<40.0\,\mathrm{GeV},\qquad
C^{\prime}_{3}=\emptyset~{}.\end{gathered}$ (17)

We run Pythia 6 [42] on these events for QED showering, followed by Delphes
3.4 with Muon Collider configuration. In figure 3 we display normalized
differential cross sections with respect to the CM rapidity $y^{*}$, the
transverse momentum $p_{T}$ and the CM energy $E^{*}$ for the final-state
photons in processes (13), (14), and the final-state $e^{\pm}$ for (15), at
three different collision energies. Those differential cross sections
correspond to detector-level events generated with the cuts (17) at the parton
level. We notice here that the differential cross sections shown in the figure
for the signal process (13) correspond to the contact interaction (7)
represented by the symbol $\sigma_{\gamma\gamma}^{(4)}$ in tables 5 and 6. For
the kinematic variables considered in figure 3, however, identical results
would have been obtained with the interaction (6). We also point out here the
good agreement of the analytical rapidity distribution (11) with the Monte
Carlo data in figure 3, right column.

Figure 3: Normalized differential cross sections for final-state photons for
processes (13), (14) and $e^{\pm}$ for process (15), with respect to: center-
of-mass energy (left column), transverse momentum (center column) and center-
of-mass rapidity (right column). The beam energies are $(E_{e},E_{\mu})=$ (100
GeV, 3 TeV) (upper row), (150 GeV, 4.5 TeV) (middle row), (200 GeV, 6 TeV)
(lower row). Blue lines: signal process (13), black line: analytical
expression (11), orange lines: background process (14), green lines:
background process (15). All differential cross section correspond to
detector-level events with cuts (17).

We then apply a preselection cut

$N_{\gamma}\geq 2,\qquad N_{e^{\pm}}=0,$ (18)

to the Delphes events. Notice that (18) implicitly includes a cut in absolute
rapidity $|y^{*}|<2.5$, corresponding to the detector rapidity acceptance
range. Finally, we apply the cuts $C_{0\mbox{--}3}$ from (16) to the
preselected events. The cross sections obtained for the detector-level events
are summarized in table 6. We notice here that the cross section for the
process (13) induced by the dim 6 operator (6), as displayed by the diagrams
in figure 1 (a), has an energy-independent cross section at the parton level,
but shows a slight decrease with increasing $\sqrt{s}$ in table 6. This is due
to the fact that the upper limit on $p_{T}^{\mathrm{tot}}$ we are using in
(16) is fixed. This effect, however, is more than compensated for by the
increase in the partonic cross section in the case of the signal process
induced by the dim 8 operator (7), as displayed in figure 1 (b), leading to a
increasing cross section also at the detector level. The cross section for the
process (14) shows a very modest growth with $\sqrt{s}$ in table 6, and
process (15) actually decreases at the highest energy.

$\mathcal{P}_{\mu}$ | $\mathcal{P}_{e}$ | $\sigma_{\gamma\gamma}^{(3)}[\mathrm{fb}]$ | $\sigma_{\gamma\gamma}^{(4)}[\mathrm{fb}]$ | $\sigma_{\nu\nu\gamma\gamma}[\mathrm{fb}]$ | $\sigma_{\nu\nu ee}[\mathrm{fb}]$
---|---|---|---|---|---
$E_{e}=$ 100 GeV, $E_{\mu}=$ 3 TeV
$0.0$ | $0.0$ | $1.72$ | $0.0082$ | $\sim 10^{-5}$ | $\sim 10^{-5}$
$+0.4$ | $+0.8$ | $2.26$ | $0.011$ | $\lesssim 10^{-5}$ | $\lesssim 10^{-5}$
$E_{e}=$ 150 GeV, $E_{\mu}=$ 4.5 TeV
$0.0$ | $0.0$ | $1.61$ | $0.039$ | $\sim 10^{-4}$ | $\sim 10^{-5}$
$+0.4$ | $+0.8$ | $2.12$ | $0.051$ | $\lesssim 10^{-4}$ | $\lesssim 10^{-5}$
$E_{e}=$ 200 GeV, $E_{\mu}=$ 6 TeV
$0.0$ | $0.0$ | $1.46$ | $0.11$ | $\sim 10^{-4}$ | $\lesssim 10^{-5}$
$+0.4$ | $+0.8$ | $1.93$ | $0.15$ | $\lesssim 10^{-4}$ | $\lesssim 10^{-5}$
Table 6: Cross sections in fb for the signal processes (13) obtained from the
operator $Q_{eA}$, see eq. (6) and fig. 1 (a), (13) obtained from the operator
$Q_{eAA}$, see eq. (7) and fig. 1 (b), and the background processes (14) and
(15). All cross sections, at the three energies shown, obtained at the
detector-simulation level with cuts $C_{0\mbox{--}3}$ from (16).

Furthemore, comparing the results in table 6 at $(E_{e},E_{\mu})=$ (100, 3000)
GeV with those in table 5, we see that detector effects result in an effective
efficiency of 84% for the $\gamma\gamma$ production processes given in the
tables by $\sigma_{\gamma\gamma}^{(3)}$, $\sigma_{\gamma\gamma}^{(4)}$. We see
also that detector efficiencies reduce the cross section for the $e^{-}e^{+}$
background process (15) to the same $O(10^{-5})\,\mathrm{fb}$ level as the
$\gamma\gamma$ background (14). Completely analogous results are obtained at
the two higher energies considered in table 6.

##  4 Conclusions

We have obtained individual limits on LFV four fermion operators by looking at
the two-to-two production processes they induce at the $e^{-}\mu^{+}$
collider. For operators $e\mu f\bar{f}$ where $f$ is a second or third family
fermion the sensitivity of the collider is of the same order of magnitude as,
and for some operators even somewhat stronger than, that of the $\mu A\to eA$
conversion in nuclei. On the other hand, the $e\mu$ collider would have higher
sensitity than the other low energy measurements $\mu\to e{\bar{e}}e$ and
$\mu\to e\gamma$ even for first family fermions. In the particular case of the
Wilson coefficients $C_{lq}^{(1,3)}$, for example, the expected sensitivity at
an $e\mu$ collider would be at least as strong as that of all low-energy
measurements, and an order of magnitude larger than that projected for Drell-
Yan processes at the HL-LHC. This leads us to expect that, given the large
number of independent effective four-fermion operators, the additional
information obtained from the collider will certainly be invaluable. In the
case of $e\mu\to\gamma\gamma$ production the limits from the collider are
significantly more stringent than those from the low energy processes for the
case of the dimension 8 $e\mu\gamma\gamma$ contact operator (7).

In the SM the $f\bar{f}$ and $\gamma\gamma$ production involves an additional
pair of neutrinos and this makes the separation of signal and background
straightforward. We have made a detailed analysis of background and signal for
the case of $\gamma\gamma$ production, including detector simulation. We
observe that with appropiate cuts on the photon energies, the photon pair
transverse momentum and the photon rapidities we can dramatically lower the SM
background with very little reduction of the signal.

We point out, finally, that our conclusions are based on the assumptions of an
integrated luminosity of 1 ab-1 and beam energies $(E_{e},E_{\mu})=$ (100 GeV,
3 TeV), (150 GeV, 4.5 TeV), (200 GeV, 6 TeV), respectively. Clearly, higher
luminosities and/or beam energies would lead to stronger sensitivities to the
contact-interaction effective couplings considered here.

#### Acknowledgments

We are grateful to Georgina Espinoza Gurriz for her assistance with our
computer hardware. We acknowlegde support from Sistema Nacional de
Investigadores de Conacyt, México. We also acknowledge that a preprint has
previously been published at arXiv [27].

## References

  * [1] R. Franceschini and M. Greco, “Higgs and BSM Physics at the Future Muon Collider,” Symmetry 13 (2021) no.5, 851 [arXiv:2104.05770 [hep-ph]].
  * [2] K. Cheung and Z. S. Wang, “Physics potential of a muon-proton collider,” Phys. Rev. D 103, 116009 (2021) [arXiv:2101.10476 [hep-ph]].
  * [3] D. Acosta and W. Li, “A Muon-Ion Collider at BNL: the future QCD frontier and path to a new energy frontier of $\mu^{+}\mu^{-}$ colliders,” [arXiv:2107.02073 [physics.acc-ph]].
  * [4] T. Aoyama, N. Asmussen, M. Benayoun, J. Bijnens, T. Blum, M. Bruno, I. Caprini, C. M. Carloni Calame, M. Cè and G. Colangelo, et al. “The anomalous magnetic moment of the muon in the Standard Model,” Phys. Rept. 887, 1-166 (2020) [arXiv:2006.04822 [hep-ph]].
  * [5] P. Banerjee, C. M. Carloni Calame, M. Chiesa, S. Di Vita, T. Engel, M. Fael, S. Laporta, P. Mastrolia, G. Montagna and O. Nicrosini, et al. “Theory for muon-electron scattering @ 10 ppm: A report of the MUonE theory initiative,” Eur. Phys. J. C 80, no.6, 591 (2020) [arXiv:2004.13663 [hep-ph]].
  * [6] A. Masiero, P. Paradisi and M. Passera, “New physics at the MUonE experiment at CERN,” Phys. Rev. D 102 (2020) no.7, 075013 [arXiv:2002.05418 [hep-ph]].
  * [7] G. Abbiendi, C. M. Carloni Calame, U. Marconi, C. Matteuzzi, G. Montagna, O. Nicrosini, M. Passera, F. Piccinini, R. Tenchini and L. Trentadue, et al. “Measuring the leading hadronic contribution to the muon g-2 via $\mu e$ scattering,” Eur. Phys. J. C 77, no.3, 139 (2017) [arXiv:1609.08987 [hep-ex]].
  * [8] C. M. Carloni Calame, M. Passera, L. Trentadue and G. Venanzoni, “A new approach to evaluate the leading hadronic corrections to the muon $g$-2,” Phys. Lett. B 746, 325-329 (2015) [arXiv:1504.02228 [hep-ph]].
  * [9] B. Abi et al. [Muon g-2], Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm, Phys. Rev. Lett. 126, no.14, 141801 (2021), [arXiv:2104.03281 [hep-ex]].
  * [10] J. Ellis, J. L. Evans, N. Nagata, D. V. Nanopoulos and K. A. Olive, “Flipped $g_{\mu}-2$,” [arXiv:2107.03025 [hep-ph]].
  * [11] N. Arkani-Hamed and K. Harigaya, “Naturalness and the muon magnetic moment,” JHEP 09, 025 (2021) [arXiv:2106.01373 [hep-ph]].
  * [12] A. Greljo, Y. Soreq, P. Stangl, A. E. Thomsen and J. Zupan, “Muonic Force Behind Flavor Anomalies,” [arXiv:2107.07518 [hep-ph]].
  * [13] A. Aboubrahim, P. Nath and R. M. Syed, “Yukawa coupling unification in an SO(10) model consistent with Fermilab (g $-$ 2) result,” JHEP 06 (2021), 002 [arXiv:2104.10114 [hep-ph]].
  * [14] P. Athron, C. Balázs, D. H. Jacob, W. Kotlarski, D. Stöckinger and H. Stöckinger-Kim, “New physics explanations of $a_{\mu}$ in light of the FNAL muon $g-2$ measurement,” JHEP 09, 080 (2021) [arXiv:2104.03691 [hep-ph]].
  * [15] J. S. Alvarado, S. F. Mantilla, R. Martinez and F. Ochoa, “A non-universal $U(1)_{X}$ extension to the Standard Model to study the $B$ meson anomaly and muon $g-2$,” [arXiv:2105.04715 [hep-ph]].
  * [16] A. E. Cárcamo Hernández, C. Espinoza, J. Carlos Gómez-Izquierdo and M. Mondragón, “Fermion masses and mixings, dark matter, leptogenesis and $g-2$ muon anomaly in an extended 2HDM with inverse seesaw,” [arXiv:2104.02730 [hep-ph]].
  * [17] A. Crivellin, C. Greub, D. Müller and F. Saturnino, “Scalar Leptoquarks in Leptonic Processes,” JHEP 02 (2021), 182 [arXiv:2010.06593 [hep-ph]].
  * [18] R. Bonciani, A. Broggio, S. Di Vita, A. Ferroglia, M. K. Mandal, P. Mastrolia, L. Mattiazzi, A. Primo, J. Ronca and U. Schubert, et al. “The two-loop four-fermion scattering amplitude in QED,” [arXiv:2106.13179 [hep-ph]].
  * [19] M. Heller, “Planar two-loop integrals for $\mathbf{\mu e}$ scattering in QED with finite lepton masses,” [arXiv:2105.08046 [hep-ph]].
  * [20] S. Y. Choi, C. S. Kim, Y. J. Kwon and S. H. Lee, “High-energy FCNC search through e mu colliders,” Phys. Rev. D 57 (1998), 7023-7026 [arXiv:hep-ph/9707483 [hep-ph]].
  * [21] V. D. Barger, S. Pakvasa and X. Tata, “Are e mu colliders interesting?,” Phys. Lett. B 415 (1997), 200-204 [arXiv:hep-ph/9709265 [hep-ph]].
  * [22] J. C. Montero, V. Pleitez and M. C. Rodriguez, “Left-right asymmetries in polarized e - mu scattering,” Phys. Rev. D 58 (1998), 097505 [arXiv:hep-ph/9803450 [hep-ph]].
  * [23] G. Cvetic and C. S. Kim, “Heavy Majorana neutrino production at electron - muon colliders,” Phys. Lett. B 461 (1999), 248-255 [erratum: Phys. Lett. B 471 (2000), 471-472] [arXiv:hep-ph/9906253 [hep-ph]].
  * [24] F. M. L. Almeida, Jr., Y. do Amaral Coutinho, J. A. Martins Simoes and M. A. B. Vale, do., “Single neutral heavy lepton production at electron muon colliders,” Phys. Lett. B 494 (2000), 273-279 [arXiv:hep-ph/0008231 [hep-ph]].
  * [25] J. K. Singhal, S. Singh and A. K. Nagawat, “Possible exotic neutrino signature in electron muon collisions,” [arXiv:hep-ph/0703136 [hep-ph]].
  * [26] F. Bossi and P. Ciafaloni, Lepton Flavor Violation at muon-electron colliders, JHEP 10, 033 (2020), [arXiv:2003.03997 [hep-ph]].
  * [27] M. Lu, A. M. Levin, C. Li, A. Agapitos, Q. Li, F. Meng, S. Qian, J. Xiao and T. Yang, The physics case for an electron-muon collider, Adv. High Energy Phys. 2021, 6693618 (2021), [arXiv:2010.15144 [hep-ph]].
  * [28] S. Davidson, Completeness and complementarity for $\mu\to e\gamma\mu\to e\bar{e}e$ and $\mu A\to eA$, JHEP 02, 172 (2021), [arXiv:2010.00317 [hep-ph]].
  * [29] M. Ardu and S. Davidson, “What is Leading Order for LFV in SMEFT?,” JHEP 08, 002 (2021) [arXiv:2103.07212 [hep-ph]].
  * [30] S. Davidson, Y. Kuno, Y. Uesaka and M. Yamanaka, Probing $\mu e\gamma\gamma$ contact interactions with $\mu\to e$ conversion, Phys. Rev. D 102, no.11, 115043 (2020), [arXiv:2007.09612 [hep-ph]].
  * [31] B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek, “Dimension-Six Terms in the Standard Model Lagrangian,” JHEP 10, 085 (2010) [arXiv:1008.4884 [hep-ph]].
  * [32] M. L. Mangano and S. J. Parke, Phys. Rept. 200, 301-367 (1991), [arXiv:hep-th/0509223 [hep-th]].
  * [33] W. H. Bertl et al. [SINDRUM II], “A Search for muon to electron conversion in muonic gold,” Eur. Phys. J. C 47, 337-346 (2006)
  * [34] R. Alonso, E. E. Jenkins, A. V. Manohar and M. Trott, “Renormalization Group Evolution of the Standard Model Dimension Six Operators III: Gauge Coupling Dependence and Phenomenology,” JHEP 04 (2014), 159 [arXiv:1312.2014 [hep-ph]];
  * [35] E. E. Jenkins, A. V. Manohar and M. Trott, “Renormalization Group Evolution of the Standard Model Dimension Six Operators II: Yukawa Dependence,” JHEP 01 (2014), 035 [arXiv:1310.4838 [hep-ph]];
  * [36] E. E. Jenkins, A. V. Manohar and M. Trott, “Naive Dimensional Analysis Counting of Gauge Theory Amplitudes and Anomalous Dimensions,” Phys. Lett. B 726 (2013), 697-702 [arXiv:1309.0819 [hep-ph]];
  * [37] E. E. Jenkins, A. V. Manohar and M. Trott, “Renormalization Group Evolution of the Standard Model Dimension Six Operators I: Formalism and lambda Dependence,” JHEP 10 (2013), 087 [arXiv:1308.2627 [hep-ph]].
  * [38] A. Angelescu, D. A. Faroughy and O. Sumensari, “Lepton Flavor Violation and Dilepton Tails at the LHC,” Eur. Phys. J. C 80 (2020) no.7, 641 [arXiv:2002.05684 [hep-ph]].
  * [39] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, B. Fuks, “FeynRules 2.0 – A complete toolbox for tree-level phenomenology,” Comput. Phys. Commun. 185 (2014) 2250 [arXiv:1310.1921 [hep-ph].
  * [40] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, T. Stelzer, “MadGraph 5 : Going Beyond,” J. H. E. Phys. 06 (2011) 128 [arXiv:1106.0522 [hep-ph]].
  * [41] J. De Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, M. Selvaggi, “delphes 3, A modular framework for the fast simulation of a generic collider experiment,” J. H. E. Phys. 02 (2014) 057 [arXiv:1307.6346 [hep-ex]].
  * [42] T. Sjöstrand, S. Mrenna, P. Skands, “pythia 6.4 physics and manual,” J. H. E. Phys. 05 (2006) 026 [arXiv:hep-ph/0603175].