# SDEX: Monte Carlo Simulation of Stochastic Differential Equations on
Memristor Crossbars

Louis Primeau1, Amirali Amirsoleimani2, and Roman Genov1
1Department of Electrical and Computer Engineering, University of Toronto,
Toronto, Canada
2Department of Electrical Engineering and Computer Science, York University,
Toronto, Canada
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Here we present stochastic differential equations (SDEs) on a memristor
crossbar, where the source of gaussian noise is derived from the random
conductance due to ion drift in the memristor during programming. We examine
the effects of line resistance on the generation of normal random vectors,
showing the skew and kurtosis are within acceptable bounds. We then show the
implementation of a stochastic differential equation solver for the Black-
Scholes SDE, and compare the distribution with the analytic solution. We
determine that the random number generation works as intended, and calculate
the energy cost of the simulation.

###### Index Terms:

Memristor, Crossbar, Stocasticity, Stochastic Differential Equation, Vector-
Matrix Multiplication, Black-Scholes.

##  I Introduction

Memristor crossbar arrays are an emerging CMOS alternative for acceleration of
problems requiring compute-intensive matrix multiplication. Vector-matrix
multiplication (VMM) operations can be done in $\mathcal{O}(1)$ time
complexity using Ohm’s law to sum analog currents. Despite their promise as
VMM accelerators, current devices suffer from nonidealities stemming from
immature fabrication techniques [1], preventing perfectly accurate VMM
operations. Because of this, several previous works have cleverly taken
advantage of these nonidealities [2][3] to create analog computing devices
beyond VMM accelerators.

Stochastic differential equations (SDEs) are used in many areas of applied
mathematics, modeling systems with inherent stochastic terms or systems in the
presence of noise [4]. Their wide use is despite the fact that there are no
analytic solutions for most models of interest, meaning that numerical methods
are necessary. Monte Carlo methods are one method of tackling the problem,
where many realizations of the process are simulated, and then quantities of
interest are computed from the resulting distributions.

However, Monte Carlo simulation suffers from two bottleneck operations:
solution of nonlinear equations for certain types of solvers and generation of
gaussian random variables [5]. Therefore, memristor crossbars can be used to
solve the latter problem, where gaussian random numbers can be generated
directly in memory, reducing time and energy cost compared to conventional
algorithms implemented on CMOS such as the Ziggurat algorithm [6][7] or the
Box-Muller transform.

Previous works on hardware acceleration of Monte Carlo methods for SDEs have
been focused on FPGAs [8] or fully analog computation [9][10], but to date no
work has been written evaluating the potential in this space for memristor
crossbars. There have been tangentially related works on MCMC algorithms on
in-memory computing platforms, such as [11] (on SRAM) and [2] (on memristor
crossbars), which are not directly applicable but provide motivation for this
work.

In this paper we examine generation of random numbers on an in memory
computing device, allowing the acceleration of Monte Carlo simulation of
stochastic differential equations. The inherent noisiness of crossbar
operations is shown to be not a barrier to their use as accelerators, but an
advantage which allows even more efficient calculation of Monte Carlo
simulations.

##  II Background

In this section we describe the characteristics of memristor crossbars as well
as basic mathematical results about stochastic differential equations and
numerical integration techniques.

###  II-A Stochasticity in Memristor Crossbars

Memristor crossbars are a type of resistive RAM that use programmable
resistance grids to perform highly efficient vector-matrix multiplication.
Because fabrication techniques are still in development for such devices, the
resolution of analog memristive devices is quite low compared to conventional
computers. The low “resolution” is due to limited programming tuning accuracy,
with $99$% of memristors within $4$% of the target conductance in the most
advanced crossbars [12].

Figure 1: (a) The conductance of a memristor in its low conductance state is
normally distributed with mean equal to the target conductance and standard
deviation equal to some percentage $v$ of the target conductance. This can be
exploited to create a gaussian random variable from the difference of the
currents in the bit lines (in the figure $i_{2}-i_{1}$). The difference of the
currents is used to encode negative numbers, and makes summation with non
random variables below simpler.(b) An example mapping of a random vector
block. The random source memristors are reprogrammed to generate random
vectors, and the outputs are multiplied by a covariance matrix to correct
differences in variance. (c) A Wiener process, the solution to the simplest
SDE $dX=dW$. The samples dW are taken from a pair of memristors on a crossbar,
which are added up over time to generate the full trajectory. The trajectory
was created by taking $1000$ time steps with value distributed as N(0,dt),
which could be drawn from the conductance distribution of the memristors.

In OxRAM devices, the stochasticity in the cycle-to-cycle variability can be
explained by variations in the conductive filaments. When the device is set
and reset, the size and location of the filament changes, causing a change in
the conductance state. Because of this random property, the high conductance
state (HCS) displays a gaussian distribution [13], with standard deviation
equal to anywhere from 10% to 40% of the mean, which can be controlled via the
write current [2].

Although read operations can cause drifts in neighboring memristors in the
same rows and columns, we neglect this effect because it is tolerable for a
single read operation[14] even at larger ($>$1V) voltages that use up the
random sample drawn from the memristor conductance distribution.

The crossbar model used in this paper was written in PyTorch [15]. It supports
variable input, output, and line resistance. The programming is modeled as
linear with finite precision enforced by perturbing the write conductance by
some percentage error, in this work a pessimistic 5%. It should be noted that
this variation is distinct from the cycle to cycle variability described
above, which occurs in a region where the IV curve is not linear and accurate
programming is very difficult.

The read step is modeled using a 1-bit DAC, i.e. floating point vector inputs
are converted into n-bit fixed precision binary inputs, and the n outputs are
scaled and added together off the crossbar [16]. All simulations in this paper
used a line resistance of $5$ $\Omega$, an input and output resistance of $1$
K$\Omega$[17], unless stated otherwise.

###  II-B Stochastic Differential Equations

Stochastic differential equations are a kind of differential equation which
involves a random term, used to model stochastic processes such as ion drift,
financial markets, and neuronal activity. The Itô stochastic differential
equation is

$dX=r(t,X)dt+\sigma(t,X)^{T}dW,$ (1)

where $r$ and $\sigma$ are vector valued functions and $dW$ is the derivative
of Brownian motion. The time evolution of the probability distribution $p$ of
$X$ follows the 2nd-order Fokker-Planck partial differential equation

$\frac{\partial p}{\partial t}=-\frac{\partial}{\partial
x}(rp)+\frac{1}{2}\frac{\partial^{2}}{\partial x^{2}}(\sigma^{2}p)$ (2)

which can be used to find analytic solutions to certain stochastic
differential equations.

###  II-C SDE Integrators

Only a small class of SDEs can be solved analytically. Just as with their
deterministic counterparts, numerical methods must be used for all but the
simplest SDEs. The one step approximation to the differential equation (1)
yields the Euler-Maruyama method:

$X_{n}-X_{n+1}=r(t_{n},X_{n})\Delta t+\sigma(t_{n},X_{n})\Delta W$ (3)

which is the simplest numerical stochastic integration method. $\Delta W$ is a
normal random vector with zero mean and standard deviation $\Delta t$. The
method has an order of convergence of $O(\Delta t^{1/2})$ where $\Delta t$ is
the time step size.

##  III Solution

In this section we describe the implementation of a stochastic integrator over
the crossbar, including generation and characterization of noise, as well as
the solver itself. A memristor crossbar simulator is used, including
input/output and line resistance, tuning variability, quantized DAC inputs,
and stuck on and stuck off nonidealities. The crossbar is modeled as a semi-
passive array, with $8\times 8$ tiles of memristors.

###  III-A Random Sampling

The integrator described above needs a source of gaussian random vectors. The
crossbar can provide this through the uncertainty in the low conductance state
of a memristor, which follows a gaussian distribution centered around the
target conductance with variability proportional to the conductance. The
variability is a characteristic of the memristor and the programming current.
A line of memristors can be programmed to the same conductance in order to
generate a gaussian random vector with a certain mean and standard deviation,
which can be rescaled to a unit normal vector. Although line resistance can
cause voltage degradation in larger passive arrays, we do not observe this to
have an effect on the distributions of the outputs, shown in Figure 2. The
pair of memristors at the end of the word-line produce the same output
variance, regardless of the line resistance.

Figure 2: A 32x32 simulated crossbar was used to generate $500$ 16-length
random gaussian vectors using the first word line. (a) Shows the change in the
distribution due to line resistance in the first pair of memristors on the top
word line. (b) Shows the same for the last pair of memristor. The mean and
standard deviation are not changed by the increase in line resistance. (c)
Shows the skew and kurtosis of the data in the top two plots. The kurtosis has
3 subtracted since the kurtosis of the unit normal distribution is 3. These
higher order moments show no correlation with the position on the word line,
indicating that line resistance does not negatively affect the properties of
the gaussian sources. For reference, the standard deviation of the sample
distribution of the kurtosis$-3$ for 500 normally distributed variables is
approximately 0.2, although the distribution is not symmetric. Therefore, most
of the samples fell within one standard deviation of the of the sample mean of
their test statistic.

For generation of gaussian random vectors with arbitrary covariance, the unit
normal random vector can be multiplied by the lower triangular matrix of the
Cholesky decomposition of the covariance matrix. This matrix multiplication
has limited precision, resulting in a gaussian distribution close to but not
exactly equal to the desired gaussian spectra. This is useful for modeling
correlated processes.

###  III-B Function Evaluation

The crossbar is not capable of computing general nonlinear functions such as
$r(t,X)$ and $g(t,X)$. Therefore, they must be evaluated on a normal computer,
using a computer-in-the-loop architecture. Linear functions can be evaluated
as usual on the crossbar, and the allocation of computational tasks to the
crossbar should be decided based on the problem. Recent work on neural
networks in the framework of stochastic differential equations are promising
for memristor crossbar acceleration because they involve many linear weight-
stationary operations [18].

##  IV Applications

###  IV-A Solving the Black-Scholes SDE

The Black-Scholes SDE (Eq 4) for geometric brownian motion is a classic
stochastic differential equation [19]. The case we consider is for a market
with a fixed interest rate $r$ and volatility $\sigma$, since this simple case
has a closed form solution.

$dX=rXdt-\sigma XdW\textrm{ subject to }X(t=0)=X_{0}$ (4)

Because it has a closed form solution (Eqn. 5), we can use it to test the
accuracy of our solver by comparing the distribution with the distribution of
the final values of the trajectories simulated on the crossbar.

$X=X_{0}\exp(\sigma W+(r-\sigma^{2}/2)t)$ (5)

The parameters $r$ and $\sigma$ were mapped to the crossbar, and a 3rd pair of
memristors set to their high resistance state (100 k$\Omega$) was used as the
random source. The variability of the HRS was set at 25%. Because the problem
is 1-D, no covariance matrix is needed. Figure 3(a) shows the layout of the
parameters on the crossbar. All unused memristors were set to their low
resistance state (10 k$\Omega$), with 10% standard deviation.

The results of Monte Carlo simulation using the Euler-Maruyama method are
shown in Figure 3. For calculation of $rX$ and $\sigma X$, the input was
encoded in 16 bits. The simulation involving only generation of random numbers
on the crossbar is in good agreement with the analytic solution. However, the
simulation involving computation of $r$ and $\sigma$ displays a slight skew.
This is due to the crossbar model perturbing the conductances to model
imperfect writing. When this parameter is turned down the skew disappears.
Although the variability is small, over $100$ time steps the effect on the
distribution can be significant. Because these values are so important, more
care is required in programming them so that they land closer to their target
conductance.

###  IV-B Static Power Consumption

Gaussian random sampling is expensive on digital devices because several
calculations are needed to transform from a uniform distribution for which
pseudo-random number generators exist to normal random variables. Because of
this, the implicit gaussian distribution on the crossbar allows for some
savings in computation despite requiring a write operation each time a new
sample is needed.

The static power consumption of the memristor array was calculated following
the model for the eVM synaptic array in NeuroSim [21], although their software
was not used. The memristors were modeled as taking $100$ of $1$ V write
pulses of $200$ ns to get from low to high resistance state, as well as a $1$
ms $0.2$ V read operation between each write pulse. This is a pessimistic
estimate because writing to the HRS does not require so many read operations
at first because it is unlikely that the resistance will overshoot the HRS.
Unselected resistors were set at the low resistance state.

For the 2-memristor random variable shown here, this resulted in an average
writing energy of 0.8 $\mu$J, which was dominated by the energy consumed by
the read operations. The full simulation of the Black-Scholes equation
involved roughly $200,000$ write operations, resulting in an energy of
consumption 0.16J. Because all unused memristors were set at the low
resistance state, the power consumption is a worst-case scenario. The energy
consumed by read operations was 3$\mu$J. Scaling to higher dimensional
problems would increase this figure linearly.

Figure 3: Simulation of the Black-Scholes equation. The equation was simulated
from t=0 to t=1, with 100 time steps. The parameters were r=0.1,
$\sigma$=$0.2$. (a) Heatmap of resistances of the devicess in a simulated 8x8
array, showing the location of $r$, $\sigma$, and the random source for
$\Delta W$. (b) Shows the first 16 bitline currents for read operations. The
first two operations are single read operations used to model the energy
consumption during writing of $r$ and $\sigma$. The subsequent read operations
are from the random source, and the stochasticity in the output currents can
be observed. In (c), the blue line shows the distribution of the final
positions of the trajectories using the crossbar random number generator and
solver weights. The black line shows the analytic solution given that W is
brownian motion. In (d) only the random numbers were generated on the
crossbar. Note the slight skew in the top plot. The orange line shows the
final positions of the trajectories calculated using fully digital
calculations.

##  V Conclusion

In this paper, simulation of a simple SDE was shown over memristor hardware,
where the white noise was drawn from the intrinsic uncertainty in programming.
The effects of crossbar nonidealities such as line resistance were examined in
order to determine their effects on the shape of the gaussian distribution.
Despite nonidealities associated with inaccurate correction for variance, the
solver produces a good fit to the analytic solution. This paper demonstrates
that programming variability can be taken advantage of for solutions of SDEs
while still being useful for vector matrix multiplication operations.

There are several directions this work can be taken in the future. Integration
of a more complete crossbar model including 1T1R configurations and stochastic
memristor models would improve the characterization of distributions resulting
from crossbar operations. Although not discussed in this paper, the methods
used here could be applied to higher dimensional SDEs where no analytical
solutions exist, where the VMM acceleration of the crossbar would be
significant.

## References

  * [1] A. Amirsoleimani, F. Alibart, V. Yon, J. Xu, M. R. Pazhouhandeh, S. Ecoffey, Y. Beilliard, R. Genov, and D. Drouin, “In‐memory vector‐matrix multiplication in monolithic complementary metal–oxide–semiconductor‐memristor integrated circuits: Design choices, challenges, and Perspectives,” Advanced Intelligent Systems, vol. 2, no. 11, Aug. 2020.
  * [2] T. Dalgaty, N. Castellani, C. Turck, K.-E. Harabi, D. Querlioz, and E. Vianello, “In situ learning using intrinsic memristor variability via Markov chain Monte Carlo sampling,” Nature Electronics, vol. 4, no. 2, pp. 151–161, Jan. 2021.
  * [3] R. Zhang, H. Jiang, Z. R. Wang, P. Lin, Y. Zhuo, D. Holcomb, D. H. Zhang, J. J. Yang, and Q. Xia, “Nanoscale diffusive memristor crossbars as physical unclonable functions,” Nanoscale, vol. 10, no. 6, pp. 2721–2726, Jan. 2018.
  * [4] B. Oksendal, Stochastic differential equations: An introduction with applications. Berlin, Germany: Springer, 2013.
  * [5] H. Gilsing and T. Shardlow, “SDELab: A package for solving stochastic differential equations in MATLAB,” Journal of Computational and Applied Mathematics, vol. 205, no. 2, pp. 1002–1018, 2007.
  * [6] G. Marsaglia and W. W. Tsang, “The Ziggurat Method for Generating Random Variables”, J. Stat. Soft., vol. 5, no. 8, pp. 1–7, Oct. 2000.
  * [7] Jianing Su and Jun Han, ”An improved Ziggurat-based hardware Gaussian random number generator,” 2016 13th IEEE International Conference on Solid-State and Integrated Circuit Technology (ICSICT), 2016, pp. 1606-1608, doi: 10.1109/ICSICT.2016.7998817.
  * [8] F. Martini, M. Piccardi, N. B. Liberati and E. Platen, ”A hardware generator for multi-point distributed random variables,” 2005 IEEE International Symposium on Circuits and Systems, 2005, pp. 1702-1705 Vol. 2, doi: 10.1109/ISCAS.2005.1464934.
  * [9] Y. Huang, N. Guo, M. Seok, Y. Tsividis, and S. Sethumadhavan, “Analog Computing in a modern context: A linear algebra accelerator case study,” IEEE Micro, vol. 37, no. 3, pp. 30–38, Jun. 2017.
  * [10] Y. Huang, N. Guo, S. Sethumadhavan, M. Seok and Y. Tsividis, ”A Case Study in Analog Co-Processing for Solving Stochastic Differential Equations,” 2018 IEEE 23rd International Conference on Digital Signal Processing (DSP), 2018, pp. 1-5, doi: 10.1109/ICDSP.2018.8631831.
  * [11] P. Shukla, A. Shylendra, T. Tulabandhula and A. R. Trivedi, ”MC2RAM: Markov Chain Monte Carlo Sampling in SRAM for Fast Bayesian Inference,” 2020 IEEE International Symposium on Circuits and Systems (ISCAS), 2020, pp. 1-5, doi: 10.1109/ISCAS45731.2020.9180701.
  * [12] H. Kim, M. R. Mahmoodi, H. Nili, and D. B. Strukov, “4K-memristor analog-grade passive crossbar circuit,” Nature Communications, vol. 12, no. 1, Aug. 2021.
  * [13] C. Li, M. Hu, Y. Li, H. Jiang, N. Ge, E. Montgomery, J. Zhang, W. Song, N. Dávila, C. E. Graves, Z. Li, J. P. Strachan, P. Lin, Z. Wang, M. Barnell, Q. Wu, R. S. Williams, J. J. Yang, and Q. Xia, “Analogue signal and image processing with large memristor crossbars,” Nature Electronics, vol. 1, no. 1, pp. 52–59, Dec. 2017.
  * [14] Y. Yilmaz and P. Mazumder, ”A Drift-Tolerant Read/Write Scheme for Multilevel Memristor Memory,” in IEEE Transactions on Nanotechnology, vol. 16, no. 6, pp. 1016-1027, Nov. 2017, doi: 10.1109/TNANO.2017.2741504.
  * [15] A. Paszke, S. Gross, and F. Massa, “PyTorch: An Imperative Style, High-Performance Deep Learning Library,” NEURIPS2019, pp. 8024–8035, 2019.
  * [16] A. Shafiee, A. Nag, N. Muralimanohar, R. Balasubramonian, J. P. Strachan, M. Hu, R. S. Williams, and V. Srikumar, “Isaac: A convolutional neural network accelerator with in-situ analog arithmetic in Crossbars,” 2016 ACM/IEEE 43rd Annual International Symposium on Computer Architecture (ISCA), Jun. 2016.
  * [17] K. Pham, S. Tran, T. Nguyen, and K.-S. Min, “Asymmetrical Training Scheme of Binary-Memristor-Crossbar-Based Neural Networks for Energy-Efficient Edge-Computing Nanoscale Systems,” Micromachines, vol. 10, no. 2, p. 141, Feb. 2019.
  * [18] B. Tzen and M. Raginsky, “Neural Stochastic Differential Equations: Deep Latent Gaussian Models in the Diffusion Limit,” CoRR, vol. abs/1905.09883, May 2019.
  * [19] R. U. Seydel, Tools for Computational Finance, 5th ed. New York, NY: Springer, 2012.
  * [20] A. Nag et al., ”Newton: Gravitating Towards the Physical Limits of Crossbar Acceleration,” in IEEE Micro, vol. 38, no. 5, pp. 41-49, Sep./Oct. 2018, doi: 10.1109/MM.2018.053631140.
  * [21] P. Chen, X. Peng and S. Yu, ”NeuroSim: A Circuit-Level Macro Model for Benchmarking Neuro-Inspired Architectures in Online Learning,” in IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 37, no. 12, pp. 3067-3080, Dec. 2018, doi: 10.1109/TCAD.2018.2789723.


Probing anomalous Zγγγ couplings at a future muon collider

H. AmarkhailElectronic address<EMAIL_ADDRESS>
Department of Physics, Sivas Cumhuriyet University, 58140,
Sivas, Turkey

and Department of Physics, Kandahar University, Kandahar 3801, Afghanistan

S.C. İnanElectronic address<EMAIL_ADDRESS>
Department of Physics, Sivas Cumhuriyet University, 58140,
Sivas, Turkey

A.V. KisselevElectronic address:

A.A. Logunov Institute for High Energy Physics, NRC
“Kurchatov Institute”,

142281, Protvino, Russian Federation

The sensitivity to anomalous quartic gauge couplings  (AQGCs) of the
$\gamma\gamma\gamma Z$ interaction is studied in the $\mu^+\mu^-
\rightarrow \mu^+\gamma\gamma \mu^-$ scattering at a future muon
collider with unpolarized beams. The anomalous $\gamma\gamma\gamma
Z$ vertex is described by two couplings, $\zeta_1$ and $\zeta_2$.
The differential and total cross sections are calculated for the
center-of-mass energies of 3 TeV, 14 TeV, and 100 TeV. For these
values of the collision energy the $95\%$ C.L. exclusion regions for
AQGCs are obtained depending on the systematic error. In particular,
for the 14 TeV muon collider with the integrated luminosity $L = 20$
ab$^{-1}$ the best sensitivities are derived to be $\zeta_1 = 3.1
\times 10^{-5}$ TeV$^{-4}$ and $\zeta_2 = 6.5 \times 10^{-5}$
TeV$^{-4}$. These constraints are three orders of magnitude stronger
than the bounds obtained for the 27 TeV HE-LHC with $L = 15$
ab$^{-1}$. At the 100 TeV muon collider with $L = 1000$ ab$^{-1}$
AQGCs can be probed up to $(1.64 \div 3.4) \times 10^{-8}$
TeV$^{-4}$. The partial-wave unitarity constraints on couplings
$\zeta_1$, $\zeta_2$ are evaluated. It is shown that the unitarity
is not violated in the region of the AQGCs examined in the present

§ INTRODUCTION

In the Standard Model (SM) there are only three quartic gauge
couplings which involve at least two charged $W$ bosons. In
particular, there is no contribution from the $Z\gamma\gamma\gamma$
coupling to the decay $Z \rightarrow \gamma\gamma\gamma$ which
proceeds exclusively via fermion and $W$-boson loops. First search
for anomalous couplings in $Z \rightarrow\gamma\gamma\gamma$ events
was made by the L3 collaboration at the LEP. The branching ratio
upper limit for the $Z$ boson decaying directly into three photons
was found to be $1.0 \times 10^{-5}$ [1]. The LEP 2 data,
as was shown in [2], induce limits of order
$10^{-5}$ GeV$^{-4}$ for the anomalous $Z\gamma\gamma\gamma$
couplings and of order $10^{-3}$ GeV$^{-4}$ for the $ZZ\gamma\gamma$
couplings. The best upper bound for the $Z
\rightarrow\gamma\gamma\gamma$ decay has been obtained by the ATLAS
collaboration using $pp$ collisions at $\sqrt{s} = 8$ TeV and equals
to $\mathcal{B}(Z \rightarrow \gamma\gamma\gamma) < 2.2 \times
10^{-6}$ [3]. Note that the SM expectation is
$\mathcal{B}(Z \rightarrow \gamma\gamma\gamma) = 5.4 \times

Analogously, there is no contribution from the $Z\gamma\gamma\gamma$
coupling to the $Z\gamma \rightarrow \gamma\gamma$ collision. Thus,
a diphoton production via VBF at $e^+e^-$ or muon colliders can be
regarded as a good process to probe effects induced by anomalous
quartic vertex $Z\gamma\gamma\gamma$.

The anomalous couplings for the $Z\gamma \rightarrow \gamma\gamma$
interaction is a particular case of the anomalous quartic gauge
couplings (AQGCs). They were searched for in a number of experiments
at the LHC, and upper bounds have been obtained
[4]-[9]. As for limits on AQGCs
for neutral bosons, they were derived both for the CLIC
<cit.>-[15] and LHC
<cit.>-[23]. The upper bounds on AQGCs were
also obtained for future $\gamma e$ and $\gamma\gamma$ colliders
[24]-[28], future hadron
[29]-[32], and $eh$ colliders
[33]-[37], see also
[38]-[41]. Recently expected limits on
AQGCs at a muon collider have been examined

Previously, we studied the anomalous $\gamma\gamma\gamma Z$
couplings through $\gamma Z$ production in $\gamma\gamma$ collisions
at the CLIC [12]. The goal of the present paper is to
probe these couplings in the diphoton production at a future muon
collider. The novelty of our work is that i) we study the
muon collider instead of $e^+e^-$ collider; ii) we consider a
lepton collision, while in [12] a $\gamma\gamma$ mode
of the $e^+e^-$ collider (namely, a collision of Compton
backscattered photons) was examined; iii) we address
center-of-mass energies which are significantly higher than the CLIC

On the assumption that a new physics mass scale $\Lambda$ is larger
than the collision energy $\sqrt{s}$, all new physics manifestations
can be described using an effective Lagrangian
$\mathcal{L}_{\mathrm{EFT}}$ valid for $\sqrt{s} \gg \Lambda$ (see,
for instance, [46]). For our goals, it is enough to use
dimension-8 operators
$W_{\mu\nu}W^{\nu\nu}W_{\rho\sigma}W^{\rho\varrho}$, and
$B_{\mu\nu}B^{\nu\nu}W_{\rho\sigma}W^{\rho\varrho}$, where $B_\mu$
and $W_\mu$ are the $U(1)_Y$ and $SU(2)$ gauge fields, respectively.
Note that the dual strength tensors are not considered. In a broken
phase (in which the Lagrangian is expressed in terms of physical
fields) a part of the effective Lagrangian describing anomalous
$Z\gamma\gamma\gamma$ interaction looks like [22]

\begin{equation}\label{Lagrangian}
\mathcal{L}_{Z\gamma\gamma\gamma} = \zeta_1 O_1 + \zeta_2 O_2 \;,
\end{equation}

with the operators

\begin{equation}\label{operators_B}
O_1 =  F^{\mu\nu}F_{\mu\nu} F^{\rho\sigma}Z_{\rho\sigma} \;, \quad
O_2 = F^{\mu\nu} F_{\nu\rho}F^{\rho\sigma}Z_{\sigma\mu} \;,
\end{equation}

where $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$, and
$Z_{\mu\nu} = \partial_\mu Z_\nu - \partial_\nu Z_\mu$. The
anomalous couplings $\zeta_1$ and $\zeta_2$ have dimension
TeV$^{-4}$. Instead of the  operator $O_2$ the
$F_{\mu\nu}\tilde{F}^{\mu\nu}F_{\rho\sigma} \tilde{Z}^{\rho\sigma}$
operator can be encountered, where $\tilde{F}^{\mu\nu} =
(1/2)\epsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}$ and
$\tilde{Z}^{\mu\nu} =
(1/2)\epsilon^{\mu\nu\alpha\beta}Z_{\alpha\beta}$. The corresponding
Lagrangians can be transformed into each other (up to a full
derivative) using equations of motion.

As mentioned above, two operators in (<ref>) can
naturally arise from the effective Lagrangian
$\mathcal{L}_{\mathrm{EFT}}$. However, we may not appeal to
$\mathcal{L}_{\mathrm{EFT}}$, but just to demand that our operators
must be dimension-8 Lorentz invariant CP-even operators constructed
from one $Z$ boson and three photon fields. Then we come to the same
set of two operators (<ref>). In such a case, anomalous
interactions (<ref>) can be considered independently of

In a number of papers (see, for instance,
[13]-[15]) bounds on couplings
$f_i/\Lambda^4$ of an unbroken Lagrangian (defined by the
fields $B_\mu$ and $W_\mu$) were searched for. In such an approach,
while obtaining bounds on any coupling, all other couplings are
assumed to be zero.

The anomalous couplings of the neutral gauge bosons in the
broken phase (like our couplings $\zeta_1, \zeta_2$) are
known to be linear combinations of the couplings $f_i/\Lambda^4$. In
our approach, to examine the neutral AQGCs, we suggest that only one
of the quartic couplings, $\gamma\gamma\gamma\gamma$,
$Z\gamma\gamma\gamma$, or $ZZ\gamma\gamma\gamma$, is nonzero. It is
equivalent to an assumption that only one particular linear
combinations of $f_i/\Lambda^4$ is nonzero. We think that both
approaches have the right to exist.

In our previous work [42] we addressed anomalous
couplings of the $\gamma\gamma\gamma\gamma$ vertex. The goal of the
present paper is to study the couplings $\zeta_1, \zeta_2$ of the
anomalous $Z\gamma\gamma\gamma$ vertex (<ref>).

§  PRODUCTION OF PHOTON PAIR AT MUON COLLIDER

As is already mentioned in Introduction, our goal is to study AQGCs
in the diphoton production at the muon collider,

\begin{equation}\label{process}
\mu^+\mu^- \rightarrow \mu^+ V_1 V_2 \mu^- \rightarrow
\mu^+\gamma\gamma \mu^- ,
\end{equation}

which goes via VBF $V_1 V_2 \rightarrow \gamma\gamma$ ($V_{1,2} =
\gamma, Z$), see Fig. <ref>. The colliding muon
beams are assumed to be unpolarized, and we sum over the photon
polarization states. This process can be regarded as an
exclusive process by requiring the outgoing muons to be
observable in the detector coverage

\begin{equation}\label{angle_cut}
10^\circ < \theta < 170^\circ \;,
\end{equation}

where $\theta$ is a scattering angle of the outgoing muons

The Feynman diagrams describing diphoton production in the
$\mu^+\mu^-$ collision via vector boson fusion.

The cross section of the process (<ref>) is

\begin{equation}\label{cs}
d\sigma = \int\limits_{\tau_{\min}}^{\tau_{\max}} \!\!d\tau
\!\!\!\int\limits_{x_{\min}}^{x_{\max}} \!\!\frac{dx}{x}
\,\sum_{V_1, V_2} \!\!f_{V_1/\mu^+}(x, Q^2) f_{V_2/\mu^-}(\tau/x,
Q^2) \,d\hat{\sigma} (V_1V_2\rightarrow \gamma\gamma) \;,
\end{equation}

\begin{equation}\label{y_z_limits}
x_{\max} = 1 - \frac{m_\mu}{E_\mu} \;, \ \tau_{\max} = \left( 1 -
\frac{m_\mu}{E_\mu} \right)^{\!2} , \ x_{\min} = \tau/x_{\max} \;, \
\tau_{\min} = \frac{p_\bot^2}{E_\mu^2} \;,
\end{equation}

$E_\mu$ is the energy of the muon beams ($4E_\mu^2 = s$), $p_\bot$
is the transverse momenta of the outgoing photons, and $m_\mu$ is
the muon mass. In (<ref>) $V_1$ and $V_2$ run over $\gamma, Z_+,
Z_-, Z_0$, where $Z_\pm$ ($Z_0$) denotes $Z$ boson with transverse
(longitudinal) helicity. The quantities $f_{V_{1,2}/\mu^\pm}(x,
Q^2)$ denote distributions of the gauge bosons in incoming muons.

For photon distributions we use the equivalent photon approximation
(EPA) [48]-[53]. In the EPA the
polarized distributions of photon inside unpolarized fermion beam
are given by [52]

\begin{align}\label{photon_spectrum}
f_{\gamma_\pm/f}(x, Q^2) = f_{\gamma/f}(x, Q^2) =
\frac{\alpha}{4\pi} \frac{1 + (1 - x)^2}{x}
\ln\frac{Q^2}{Q^2_{\min}} \;,
\end{align}

where $x = E_\gamma/E_\mu$ is fraction of the beam energy carried by
the photon, and $Q^2 = sx/2$. In our case $f = \mu^-$, $\bar{f} =
\mu^+$. Because of C and P invariance, $f_{\gamma/\bar{f}} =
f_{\gamma/f}$. For the $Z$ boson we use so-called effective $W$
approximation (EWA) [54]-[64]. It allows
to treat of massive vector bosons as partons inside the colliding
beams and reads [64]

\begin{align}
f_{Z_\pm/f}(x, Q^2) &= \frac{\alpha_Z}{4\pi} \frac{(g_V^f \mp
g_A^f)^2 + (g_V^f \pm g_A^f)^2(1 - x)^2}{x}
\ln\frac{Q^2}{Q^2_{\min}} \;,
\label{Z_trans} \\
f_{Z_0/f}(x, Q^2) &= \frac{\alpha_Z}{\pi} \frac{[(g_V^f)^2 +
(g_A^f)^2](1 - x)}{x} \;, \label{Z_long}
\end{align}

\begin{equation}\label{alpha_Z}
\alpha_Z = \frac{\alpha}{(\cos\theta_W \sin\theta_W)^2} \;, \ \
g_V^f = \frac{1}{2}(T_3^f)_L - Q^f \sin^2 \theta_W \;, \ \  g_A^f =
- \frac{1}{2}(T_3^f)_L  \;.
\end{equation}

There are CP relations $f_{Z_\pm/\bar{f}} = f_{Z_\mp/f}$,
$f_{Z_0/\bar{f}} = f_{Z_0/f}$.

Without constraints on the scattering angles $\theta$ of the
outgoing muons, $Q^2_{\min}$ has the following values

\begin{equation}\label{Q2_min_no_cuts}
Q^2_{\min} =
\left\{
\begin{array}{ll}
m_\mu^2 x^2/(1 - x), & \hbox{in eq.~\eqref{photon_spectrum};}
\\
m_Z^2, & \hbox{in eq.~\eqref{Z_trans}.}
\end{array}
\right.
\end{equation}

For the exclusive diphoton production (<ref>) with cuts
(<ref>) one has to put in eqs. (<ref>),
(<ref>) the quantity

\begin{equation}\label{Q2_min_with_cuts}
Q^2_{\min} = \frac{s(1-x)}{2} \!\left( \frac{1}{\cos\theta_{\min}} -
1 \right) ,
\end{equation}

where $\theta_{\min} = 10^\circ$.

In (<ref>) the VBF cross section is a sum of squared helicity

\begin{equation}\label{subprocess_cs}
\frac{d\hat{\sigma}}{d\Omega}(V_1V_2\rightarrow \gamma\gamma) =
\frac{1}{64\pi^2 \hat{s}} \sum_{\lambda_3,\lambda_4}
\!\!|M_{\lambda_1\lambda_2\lambda_3\lambda_4}|^2 ,
\end{equation}

where $\lambda_1, \lambda_2$ are helicity of the incoming bosons,
and $\lambda_3, \lambda_4$ are helicities of the outgoing photons.
its turn, each of the helicity amplitudes
$M_{\lambda_1\lambda_2\lambda_3\lambda_4}$ is a sum of an anomaly
and SM terms,

\begin{equation}\label{anom+SM}
M_{\lambda_1\lambda_2\lambda_3\lambda_4} =
M_{\lambda_1\lambda_2\lambda_3\lambda_4}^{\mathrm{anom}} +
M_{\lambda_1\lambda_2\lambda_3\lambda_4}^{\mathrm{SM}} \;.
\end{equation}

To obtain constraints directly on the couplings $\zeta_1$ and
$\zeta_2$ in the effective Lagrangian (<ref>), we assume
that there is no new quartic interactions of neutral gauge bosons
except for the $Z\gamma\gamma\gamma$ vertex. For practical
use it means that the anomalous contribution to $d\hat{\sigma}
(V_1V_2\rightarrow \gamma\gamma)$ in (<ref>) comes from the VBF
process $Z\gamma\rightarrow\gamma\gamma$ only ($V_1= \gamma, V_2 =
Z$ or $V_1 = Z, V_2 = \gamma$ in Fig. <ref>). The
explicit expressions for
$M_{\lambda_1\lambda_2\lambda_3\lambda_4}^{\mathrm{anom}}$ are given
in Appendix A.

As for a SM contribution to the cross section, all VBF processes,
like $\gamma\gamma\rightarrow\gamma\gamma$, $\gamma
Z\rightarrow\gamma\gamma$ ($Z\gamma\rightarrow\gamma\gamma$), and
$ZZ\rightarrow\gamma\gamma$, should be taken into account (see
Fig. <ref>). Each of the SM amplitudes in
(<ref>) is a sum of the fermion and $W$ boson one-loop

\begin{equation}\label{f+W_ew}
M^{\mathrm{SM}} = M^{\mathrm{SM}}_f + M^{\mathrm{SM}}_W  \;.
\end{equation}

The analytical expressions for the SM helicity amplitudes
$M^{\mathrm{SM}}_f$ and $M^{\mathrm{SM}}_W$ are well-known, see
[65]-[69]. Note that
$M^{\mathrm{SM}}_W$ dominates over $M^{\mathrm{SM}}_f$ at $\sqrt{s}
> 200$ GeV.

§ NUMERICAL ANALYSIS

The differential cross sections for the collision energy of
$\sqrt{s} = 3$ TeV, 14 TeV, and 100 TeV are presented in
Fig. <ref> versus diphoton invariant mass
$m_{\gamma\gamma}$. The cuts on the rapidity of the outgoing
photons, $|\eta| < 2.5$, and on the photon transverse momenta,
$p_\bot > 80$ GeV, were applied. For each value of $s$ we present
curves for two sets of the anomalous couplings, ($\zeta_1 \neq 0,
\zeta_2 = 0$) and ($\zeta_1 = 0, \zeta_2 \neq 0$). The SM
contributions to the cross sections are also shown. Because of the
cut (<ref>), the dominant background is the SM process
$\mu^+\mu^- \rightarrow \mu^+ \gamma\gamma \mu^{+}$ [47].

The total cross sections as functions of the minimal value of the
diphoton invariant mass $m_{\gamma\gamma}$ are given in
Fig. <ref>. We see that for $\sqrt{s} = 3$ TeV the total
cross section starts to dominate the SM one at approximately
$m_{\gamma\gamma} \gtrsim \sqrt{s}/2$. For higher collision energies
of 14 TeV and 100 TeV it takes place already if $m_{\gamma\gamma}
\gtrsim \sqrt{s}/7$.

The differential cross sections for the  $\mu^+\mu^-
\rightarrow \mu^+ \gamma\gamma \mu^-$ scattering at the future muon
collider versus diphoton invariant mass $m_{\gamma\gamma}$. The
center-of-mass energy is equal to 3 TeV (left panel), 14 TeV (middle
panel), and 100 TeV (right panel).

The total cross sections $\sigma(m_{\gamma\gamma}
> m_{\gamma\gamma,\min})$ for the $\mu^+\mu^- \rightarrow \mu^+ \gamma\gamma \mu^-$ scattering at
the future muon collider versus minimal value of the invariant mass
of the photon pair $m_{\gamma\gamma}$.

Let $s(b)$ be the total number of signal (background) events, and
$\delta$ the percentage systematic error. Then the exclusion
significance is given by [70]

\begin{equation}\label{S_excl}
S_{\mathrm{excl}} = \sqrt{ 2\left[ s - b \ln \left( \frac{b + s +
x}{2b} \right) - \frac{1}{\delta^2}\ln \left( \frac{b - s + x}{2b}
\right) - (b + s -x) \left( 1 + \frac{1}{\delta^2 b} \right) \right]
} ,
\end{equation}

\begin{equation}\label{x}
x = \sqrt{ (s+b)^2 - \frac{4s\delta^2 b^2}{(1 + \delta^2 b)} } \;.
\end{equation}

We define the regions $S_{\mathrm{excl}} \leqslant 1.645$ as the
regions that can be excluded at the 95% C.L. To reduce the SM
background, we additionally impose the lower cut on the invariant
mass of the outgoing photons. Namely, we assume that the minimal
value of $m_{\gamma\gamma}$ is equal to 2 TeV, 8 TeV, and 50 TeV for
$\sqrt{s} = 3$ TeV, 14 TeV, and 100 TeV, respectively. Following
[71, 72], we assumed that for $p_\bot > 80$ GeV
(our cut) the muon reconstruction efficiency is larger than $85\%$ .

The 95% C.L. exclusion regions for the couplings $\zeta_1,
\zeta_2$ in the unpolarized $\mu^+\mu^- \rightarrow
\mu^+\gamma\gamma \mu^-$ scattering at the future muon collider with
the systematic errors $\delta = 0\%$ (blue ellipse), $\delta = 5\%$
(red ellipse), and $\delta = 10\%$ (green ellipse). The inner
regions of the ellipses are inaccessible. The collision energy is
$\sqrt{s} = 3$ TeV, the integrated luminosity is $L = 1$ ab$^{-1}$.
The cut on the diphoton invariant mass $m_{\gamma\gamma} > 2$ TeV is

The same as in Fig. <ref>, but for $\sqrt{s}
= 14$ TeV and $L = 20$ ab$^{-1}$. The cut $m_{\gamma\gamma} > 8$ TeV
is used.

The same as in Fig. <ref>, but for $\sqrt{s}
= 100$ TeV and $L = 1000$ ab$^{-1}$. The cut $m_{\gamma\gamma} > 50$
TeV is imposed.

Our 95% C.L. exclusion regions for AQGCs for the unpolarized
process $\mu^+\mu^- \rightarrow \mu^+ \gamma\gamma \mu^-$ at the
muon collider are shown in
Figs. <ref>-<ref>. Note that in the
energy region $s \gg m_Z^2$ the squared amplitude is proportional to
the coupling combination $48\zeta_1^2 + 11\zeta_1\zeta_2 +
40\zeta_2^2$. In such a case, the exclusion regions are ellipses in
the plane $\zeta_1 - \zeta_2$, as we obtained. The best bounds on
AQGCs when one of the couplings $\zeta_1, \zeta_2$ is taken to be
zero are given in Tab. <ref>.

The 95% C.L. exclusion bounds on the anomalous
quartic couplings $\zeta_1$ and $\zeta_2$, with the integrated luminosity of 1 ab$^{-1}$,
20 ab$^{-1}$, and 1000 ab$^{-1}$ for the 3 TeV, 14 TeV, and 100 TeV muon collider.

3 TeV     14 TeV     100 TeV

$|\zeta_1|$, TeV$^{-4}$
$(\zeta_2=0)$     $\delta=0\%$
$\ \, \delta=10\%$
$\makecell{ 9.90\times10^{-3} \\ 1.18\times10^{-2} \\ 1.46\times10^{-2}}$
$\makecell{ 3.13\times10^{-5} \\ 5.75\times10^{-5} \\ 8.01\times10^{-5}}$
$\makecell{ 1.64\times10^{-8} \\ 4.79\times10^{-8} \\ 6.61\times10^{-8}}$

$|\zeta_2|$, TeV$^{-4}$
$(\zeta_1=0)$     $\delta=0\%$
$\ \, \delta=10\%$
$\makecell{ 2.07\times10^{-2} \\ 2.46\times10^{-2} \\ 3.06\times10^{-2}}$
$\makecell{ 6.54\times10^{-5} \\ 1.53\times10^{-4} \\ 1.67\times10^{-4}}$
$\makecell{ 3.43\times10^{-8} \\ 1.02\times10^{-7} \\ 1.38\times10^{-7}}$

§ UNITARITY BOUNDS ON ANOMALOUS COUPLINGS

Now let us derive constraints imposed by partial-wave unitarity. To
do it, we can use the partial-wave expansion of the helicity
amplitude in the center-of-mass system [73]

\begin{align}\label{helicity_ampl_expansion}
M_{\lambda_1\lambda_2\lambda_3\lambda_4}(s, \theta, \varphi) &=
16\pi \sum_J (2J + 1) \sqrt{(1 + \delta_{\lambda_1\lambda_2})(1 +
\delta_{\lambda_3\lambda_4})}
\nonumber \\
&\times \,e^{i(\lambda - \mu)\phi} \,d^J_{\lambda\mu}(\theta)
\,T^J_{\lambda_1\lambda_2\lambda_3\lambda_4}(s) \;,
\end{align}

where $\lambda = \lambda_1 - \lambda_2$, $\mu = \lambda_3 -
\lambda_4$, $\theta(\phi)$ is the polar (azimuth) scattering angle.
If we choose the plane $(x - z)$ as a scattering plane, then $\phi =
0$ in (<ref>). Here
$d^J_{\lambda\mu}(\theta)$ is the Wigner (small) $d$-function
[74]. In particular, we have [12]

\begin{equation}\label{d-function_00}
d^J_{00}(z) = P_J(z) \;,
\end{equation}

\begin{equation}\label{d-function_22}
d^J_{2 \pm2}(z) = (\pm 1)^J \!\left( \frac{1 \pm z}{2}\right)^{\!2}
\!{}_2F_1 \!\left( 2-J, J+3; 1; \frac{1 \mp z}{2} \right) ,
\end{equation}

where $z = \cos\theta$, $P_J(z)$ is the Legendre polynomial, and
$_2F_1(a,b;c;x)$ is the hypergeometric function. It is enough to
calculate unitarity bounds on our anomalous couplings.

Partial-wave unitarity requires that

\begin{equation}\label{parity_wave_unitarity}
\left| T^J_{\lambda_1\lambda_2\lambda_3\lambda_4}(s) \right| \leq 1
\;.
\end{equation}

Note that $T^J_{\lambda_1\lambda_2\lambda_3\lambda_4}(s)$ is related
to the partial-wave amplitude with the same helicities by formula

\begin{align}\label{parity wave_func}
T^J_{\lambda_1\lambda_2\lambda_3\lambda_4}(s) &= \frac{1}{32\pi}
\frac{1}{\sqrt{(1 + \delta_{\lambda_1\lambda_2})(1 +
\delta_{\lambda_3\lambda_4})}} \int\limits_{-1}^1 \!\!
M_{\lambda_1\lambda_2\lambda_3\lambda_4}(s, z)
\,d^{J}_{\lambda\mu}(z) \,dz \;.
\end{align}

The anomalous helicity amplitudes
$M_{\lambda_1\lambda_2\lambda_3\lambda_4}$ for process
(<ref>) are calculated in Appendix A.

In the limit $s \gg m_Z^2$ the most stringent bounds come from the
helicity amplitudes $M_{++++}$, $M_{++--}$, $M_{+-+-}$, and
$M_{+--+}$ (as well as from $M_{-+++}$, $M_{-+--}$, $M_{--+-}$, and
$M_{---+}$). On the other hand, these amplitudes coincide with
$\gamma\gamma\gamma\gamma$ anomalous amplitudes up to small
corrections $\mathrm{O(m_Z^2/s)}$, see, for instance, Appendix in
[42]. It means that unitarity bounds derived in
[42] are also valid for our anomalous couplings. They
look like

\begin{equation}\label{unitariry_bounds}
|\zeta_1| < 2\pi s^{-2} \;, \quad |\zeta_2| \leq (16\pi/3) s^{-2}
\;.
\end{equation}

As we see, the unitarity demands that the coupling $|\zeta_1|$ must
be less than $7.8 \times 10^{-2}$ TeV$^{-4}$, $1.64 \times 10^{-4}$
TeV$^{-4}$, and $6.28 \times 10^{-8}$ TeV$^{-4}$ for $\sqrt{s} = 3$
TeV, 14 TeV, and 100 TeV, respectively. Correspondingly, the
unitarity bounds for the coupling $|\zeta_2|$ are $0.21$ TeV$^{-4}$,
$4.36 \times 10^{-4}$ TeV$^{-4}$, and $1.68 \times 10^{-7}$
TeV$^{-4}$. These constraints should be compared with the bounds on
AQGCs presented in Tab. <ref>. Let us note that
in eq. (<ref>) one actually has to use the
invariant energy square $\hat{s}$ of the subprocess $V_1 V_2
\rightarrow \gamma\gamma$ instead of the collision energy $s$. To
conclude, the unitarity is not violated in the region of AQGCs
studied in the present paper.

§ CONCLUSIONS

In the present paper we have probed the quartic gauge couplings
(QGCs) of the anomalous $\gamma\gamma\gamma Z$ interaction in the
$\mu^+\mu^- \rightarrow \mu^+ \gamma\gamma \mu^-$ scattering at the
high energy muon collider with unpolarized beams. The sum over the
photon polarization states are also assumed. Our process is the
exclusive one by requiring the outgoing muons to be observable in
the detector range $10^\circ < \theta < 170^\circ$. The collision
energy of $\sqrt{s} = 3$ TeV, 14 TeV, and 100 TeV with the
integrated luminosity of $L = 1$ ab$^{-1}$, 20 ab$^{-1}$, and 1000
ab$^{-1}$, respectively, are addressed. Both differential and total
cross sections are calculated depending on the invariant mass of the
diphoton system $m_{\gamma\gamma}$. We have used the cut on the
photon rapidity, $|\eta| < 2.5$. The cut on the photon transverse
momenta, $p_\bot > 80$ GeV, has been also imposed. The 95% C.L.
exclusion regions for the AQGCs $\zeta_1$ and $\zeta_2$ are
calculated. They are ellipses in the $\zeta_1 - \zeta_2$ plane, see
Figs. <ref>-<ref>. The limits on
$\zeta_1$ (for $\zeta_2 = 0$) and $\zeta_2$ (for $\zeta_2 = 0$) are
collected in Tab. 1 for three values of the systematic error (0%,
5%, and 10%). We have investigated the constraints imposed by the
partial-wave unitarity. It has been shown that it is not violated
for the AQGCs values examine in our paper.

As was already mentioned in Introduction, bounds on the quartic
anomalous couplings were obtained for future hadron colliders
[29]-[32]. In particular, the 95% C.L.
limits on AQGCs for the LHC an FCC-hh colliders are collected in
Tab. IV in [32]. One can see that for the 14 TeV HL-LHC
with the integrated luminosity of $L = 3$ ab$^{-1}$ the best 95%
C.L. limit on AQGCs is equal to $2.1 \times 10^{-1}$ TeV$^{-4}$. The
most stringent bound for the HE-LHC with the center-of-mass energy
$\sqrt{s} = 27$ TeV and $L = 15$ ab$^{-1}$ is obtained to be $2.8
\times 10^{-2}$ TeV$^{-4}$. Finally, limits on AQGCs for the future
FCC-hh collider are derived in [32], the best one is
equal to $5.4 \times 10^{-4}$ TeV$^{-4}$. To compare, our best bound
on the coupling $\zeta_1$ ($\zeta_2$) for the 14 TeV muon collider
with $L = 20$ ab$^{-1}$ is equal to 3.13 (6.54) $\times 10^{-5}$
TeV$^{-4}$. All this demonstrates a great potential of the future
muon collider in probing anomalous interactions of the neutral gauge

§ APPENDIX A. ANOMALOUS HELICITY AMPLITUDES

First consider helicity amplitudes
$\bar{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}$ of the process with
an outgoing $Z$ boson,

\begin{equation}\label{process_Z_incoming}
\gamma(p_1,\lambda_1) \gamma(p_2,\lambda_2) \rightarrow
\gamma(p_3,\lambda_3) Z(p_4,\lambda_4) \;.
\end{equation}

They can be obtained from explicit expressions given in our paper
[12]. There the helicity amplitudes
$\tilde{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}$ were derived from
the Lagrangian

\begin{equation}\label{Lagrangian_IK}
\mathcal{L}_{\gamma\gamma\gamma Z} = g_1 F^{\rho\mu} F^{\alpha\nu}
\partial_\rho F_{\mu\nu} Z_\alpha + g_2 F^{\rho\mu} F^\nu_\mu \partial_\rho F_{\alpha\nu}
Z^\alpha  \;.
\end{equation}

They look like

\begin{equation}\label{helicity_ampl_IK}
\tilde{M}_{\lambda_1\lambda_2\lambda_3\lambda_4} = g_1
\tilde{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}^{(1)} + g_2
\tilde{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}^{(2)} \;,
\end{equation}

with the explisit expressions for
$\tilde{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}^{(i)}$ ($i=1,2$)
given in Appendix B of [12]. Then, according to
eq. (2.7) in [12], the helicity amplitude
$\bar{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}$ corresponding to
the effective Lagrangian (<ref>) is defined as

\begin{equation}\label{helicity_amlitudes_relation}
\bar{M}_{\lambda_1\lambda_2\lambda_3\lambda_4} = -2(4\zeta_1 +
\zeta_2) \tilde{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}^{(1)} +
2\zeta_2 \tilde{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}^{(2)} \;.
\end{equation}

Bose-Einstein statistics and parity invariance demand that there
exist six independent helicity amplitudes
$\bar{M}_{\lambda_1\lambda_2\lambda_3\lambda_4}$ with $\lambda_1 =
+1$ and $\lambda_4 = \pm 1$

\begin{align}\label{independent_helicity_ampl}
\bar{M}_{++++}(s,t,u) &= -\frac{(4\zeta_1 + 3\zeta_2)}{2} \,s(u+t)
\;,
\nonumber \\
\bar{M}_{+++-}(s,t,u) &= 0 \;,
\nonumber \\
\bar{M}_{++-+ }(s,t,u) &= -(4\zeta_1 + \zeta_2) \frac{tu}{t+u} m_Z^2
\;,
\nonumber \\
\bar{M}_{+-++}(s,t,u) &= -\frac{(4\zeta_1 + 3\zeta_2)}{2}
\frac{tu}{t+u} m_Z^2 \;,
\nonumber \\
\bar{M}_{++--}(s,t,u) &= -(4\zeta_1 + \zeta_2) \frac{s(t^2 + tu +
u^2)}{t+u} \;,
\nonumber \\
\bar{M}_{+-+-}(s,t,u) &= -\frac{(4\zeta_1 + 3\zeta_2)}{2}
\,\frac{su^2}{t+u} \;.
\end{align}

where $s$, $t$, $u$ are Mandelstam variables, and $s+t+u = m_Z^2$.
We also have the crossing relations

\begin{align}\label{dependent_helicity_ampl}
\bar{M}_{+--+}(s,t,u) &= \bar{M}_{+-+-}(s,u,t) = -\frac{(4\zeta_1 +
3\zeta_2)}{2} \,\frac{st^2}{t+u} \;,
\nonumber \\
\bar{M}_{+---}(s,t,u) &= \bar{M}_{+-++}(s,u,t) = -\frac{(4\zeta_1 +
3\zeta_2)}{2} \frac{tu}{t+u} m_Z^2 \;.
\end{align}

Equations (<ref>),
(<ref>) represent eight helicity amplitudes
with $\lambda_1 = +1$ and $\lambda_4 = \pm 1$. For $\lambda_1 = +1$,
$\lambda_4 = 0$ we get

\begin{align}\label{ingoing_helicity_ampl_0}
\bar{M}_{+++0}(s,t,u) &= 0 \;,
\nonumber \\
\bar{M}_{++-0}(s,t,u) &= - \frac{i(4\zeta_1 + \zeta_2)}{\sqrt{2}}
\sqrt{stu} \,\frac{t-u}{t+u} m_Z \;,
\nonumber \\
\bar{M}_{+-+0}(s,t,u) &= \frac{i(4\zeta_1 + 3\zeta_2)}{\sqrt{2}}
\frac{u\sqrt{stu}}{t+u} m_Z \;,
\nonumber \\
\bar{M}_{+--0}(s,t,u) &= \bar{M}_{+-+0}(s,u,t) = \frac{i(4\zeta_1 +
3\zeta_2)}{\sqrt{2}} \frac{t\sqrt{stu}}{t+u} m_Z \;.
\end{align}

The helicity amplitudes with $\lambda_1 = -1$ can be obtained from
them by use of parity relation,

\begin{equation}\label{parity_relations_1}
\bar{M}_{-\lambda_1\lambda_2\lambda_3\lambda_4}(s,t,u) =
\bar{M}_{\lambda_1-\lambda_2-\lambda_3-\lambda_4}(s,t,u) \;.
\end{equation}

Now consider the process with an incoming $Z$ boson,

\begin{equation}\label{process_Z_outgoing}
\gamma(p_1,\lambda_1) Z(p_2,\lambda_2) \rightarrow
\gamma(p_3,\lambda_3) \gamma(p_4,\lambda_4) \;.
\end{equation}

If time-inversion invariance holds, we have the relation

\begin{equation}\label{time_inversion_relations}
M_{\lambda_1\lambda_2\lambda_3\lambda_4}(s,t,u) =
\bar{M}_{\lambda_3\lambda_4\lambda_1\lambda_2}(s,t,u) \;.
\end{equation}

Then using eqs. (<ref>) and
(<ref>) we find that the helicity amplitudes
of the process (<ref>) are given by

\begin{align}\label{helicity_ampl_Z_outgoing}
M_{++++}(s,t,u) &= \bar{M}_{++++}(s,t,u) = -\frac{(4\zeta_1 +
3\zeta_2)}{2} \,s(t+u) \;,
\nonumber \\
M_{+++-}(s,t,u) &= \bar{M}_{+-++}(s,u,t) = -\frac{(4\zeta_1 +
3\zeta_2)}{2} \frac{tu}{t+u} m_Z^2 \;,
\nonumber \\
M_{++-+ }(s,t,u) &= \bar{M}_{+---}(s,t,u) = -\frac{(4\zeta_1 +
3\zeta_2)}{2} \frac{tu}{t+u} m_Z^2 \;,
\nonumber \\
M_{+-++}(s,t,u) &= \bar{M}_{+++-}(s,t,u) = 0 \;,
\nonumber \\
M_{++--}(s,t,u) &= \bar{M}_{++--}(s,t,u) = -(4\zeta_1 + \zeta_2)
\,\frac{s(t^2 + tu + u^2)}{t+u} \;,
\nonumber \\
M_{+-+-}(s,t,u) &=  \bar{M}_{+-+-}(s,t,u) = -\frac{(4\zeta_1 +
3\zeta_2)}{2} \,\frac{su^2}{t+u} \;,
\nonumber \\
M_{+--+}(s,t,u) &= \bar{M}_{+--+}(s,t,u) = -\frac{(4\zeta_1 +
3\zeta_2)}{2} \,\frac{st^2}{t+u} \;,
\nonumber \\
M_{+---}(s,t,u) &= \bar{M}_{++-+ }(s,t,u) = -(4\zeta_1 + \zeta_2)
\frac{tu}{t+u} m_Z^2 \;,
\end{align}

\begin{align}\label{helicity_ampl_Z_outgoing_0}
M_{+0++}(s,t,u) &= \bar{M}_{+++0}(s,t,u) = 0 \;,
\nonumber \\
M_{+0+-}(s,t,u) &= \bar{M}_{+-+0}(s,t,u) = \frac{i(4\zeta_1 +
3\zeta_2)}{\sqrt{2}} \frac{u\sqrt{stu}}{t+u} m_Z \;,
\nonumber \\
M_{+0-+}(s,t,u) &= -\bar{M}_{+--0}(s,t,u) = - \frac{i(4\zeta_1 +
3\zeta_2)}{\sqrt{2}} \frac{t\sqrt{stu}}{t+u} m_Z \;,
\nonumber \\
M_{+0--}(s,t,u) &= -\bar{M}_{++-0}(s,t,u)  = \frac{i(4\zeta_1 +
\zeta_2)}{\sqrt{2}} \sqrt{stu} \,\frac{t-u}{t+u} m_Z \;.
\end{align}

The helicity amplitudes $M_{\lambda_1\lambda_2\lambda_3\lambda_4}$
with $\lambda_1 = -1$ can be obtained from
(<ref>), (<ref>)
by use of parity relation

\begin{equation}\label{parity_relations}
M_{-\lambda_1\lambda_2\lambda_3\lambda_4}(s,t,u) =
M_{\lambda_1-\lambda_2-\lambda_3-\lambda_4}(s,t,u) \;.
\end{equation}

The amplitudes
(<ref>)-(<ref>) will be
used for our numerical calculations.

[1]
M. Acciarri et al. (L3 Collaboration), Search for anomalous
$Z\rightarrow\gamma\gamma\gamma$ events at LEP, Phys. Lett. B
345, 609 (1995).

[2]
A. Gutiérrez-Rodríguez, C.G. Honorato, J. Montaño and
M.A. Pérez, Limits on the quartic couplings
$Z\gamma\gamma\gamma$ and $ZZ\gamma\gamma$ from $e^+e^-$ colliders,
Phys. Rev. D 89, 034003 (2014) [arXiv:1304.7410].

[3]
ATLAS Collaboration, Search for new phenomena in events with
at least three photons collected in pp collisions at $\sqrt{s} = 8$
TeV with the ATLAS detector, Eur. Phys. J. C 76, 210
(2016) [arXiv:1509.05051].

[4]
G. Aad et al. (ATLAS Collaboration), Measurements of
$Z\gamma$ and $Z\gamma\gamma$ production in $pp$ collisions at
$\sqrt{s} = 8$ TeV with the ATLAS detector, Phys. Rev. D
93, 112002 (2006) [arXiv:1604.05232].

[5]
G. Aad et al. (ATLAS Collaboration), Search for new phenomena
in events with at least three photons collected in $pp$ collisions
at $\sqrt{s} = 8$ TeV with the ATLAS detector, Eur. Phys. J. C
76, 210 (2016) [arXiv:1509.05051].

[6]
ATLAS Collaboration, Measurement of electroweak
$Z\nu\bar{\nu}jj$ production and limits on anomalous quartic gauge
couplings in pp collisions at $\sqrt{s} =13$ TeV with the ATLAS
detector, JHEP 06, 082 (2023) [arXiv:2208.12741].

[7]
A.M. Sirunyuan et al. (CMS Collaboration), Measurement
of the cross section for electroweak production of a $Z$ boson, a
photon and two jets in proton-proton collisions at $\sqrt{s} = 13$
TeV and constraints on anomalous quartic couplins, JHEP
06, 076 (2020) [arXiv:2002.09902].

[8]
CMS Collaboration, Measurements of the $pp \rightarrow
W^\pm\gamma\gamma$ and $pp \rightarrow Z\gamma\gamma$ cross sections
at $\sqrt{s} = 13$ TeV and limits on anomalous quartic gauge
couplings, JHEP 10, 174 (2021) [arXiv:2105.12780].

[9]
The CMS and TOTEM Collaborations, Search for high-mass
exclusive $\gamma\gamma\rightarrow WW$ and $\gamma\gamma\rightarrow
ZZ$ production in proton-proton collisions at $\sqrt{s} = 13$ TeV,
JHEP 07, 229 (2023) [arXiv:2211.16320].

[10]
M. Köksal, Anomalous quartic $ZZ\gamma\gamma$ couplings at
the CLIC, Eur. Phys. J. Plus 130, 75 (2015)

[11]
S.C. İnan and A.V. Kisselev, Probing anomalous quartic
$\gamma\gamma\gamma\gamma$ couplings in light-by-light collisions at
the CLIC, Eur. Phys. J. C 81, 664 (2021)

[12]
S.C. İnan and A.V. Kisselev, Probing anomalous
$\gamma\gamma\gamma Z$ couplings through $\gamma Z$ production in
$\gamma\gamma$ collisions at the CLIC, JHEP 10, 121 (2021)

[13]
M. Köksal, Search for the anomalous quartic gauge
couplings through $Z\gamma$ production at $e^+e^-$ colliders,

[14]
V. Ari, E. Gurkanli, M. Köksal, A. Gutiérrez-Rodríguez,
and M.A. Hernández-Ruíz, Study of the projected
sensitivity on the anomalous quartic gauge couplings via
$Z\gamma\gamma$ production at the CLIC, Nucl. Phys. B 989,
116133 (2023) [arXiv:2112.03948].

[15]
E. Gurkanli, Sensitivities on the anomalous quartic
$\gamma\gamma\gamma\gamma$ and $\gamma\gamma\gamma Z$ couplings at
the CLIC, arXiv:2307.01326.

[16]
O.J.P. Eboli, M.C. Gonzalez-Garcia, and S.M. Lietti, Bosonic
quartic couplings at CERN LHC, Phys. Rev. D 69, 095005
(2004) [arXiv:hep-ph/0310141].

[17]
I. Sahin and B. Sahin, Anomalous quartic $ZZ\gamma\gamma$
couplings in gamma-proton collision at the LHC, Phys. Rev. D
86, 115001 (2012) [arXiv:1211.3100].

[18]
S. Fichet, G. von Gersdorff, B. Lenzi, Ch. Royond, and M. Saimpertd,
Light-by-light scattering with intact protons at the LHC: from
Standard model to new physics, JHEP 02, 165 (2015)

[19]
E. Chapon, O. Kepka and C. Royon, Probing $WW\gamma\gamma$ and
$ZZ\gamma\gamma$ quartic anomalous couplings with 10 pb$^{-1}$ at
the LHC, arXiv:0908.1061.

[20]
E. Chapon, C. Royon and O. Kepka, Anomalous quartic
$WW\gamma\gamma$, $ZZ\gamma\gamma$, and trilinear $WW\gamma$
couplings in two-photon processes at high luminosity at the LHC,
Phys. Rev. D 81, 074003 (2010).

[17]
İ. Şahin and B. Şahin, Anomalous quartic
$ZZ\gamma\gamma$ couplings in gamma-proton collision at the LHC,
Phys. Rev. D 86, 115001 (2012) [arXiv:1211.3100].

[21]
A. Senol, Anomalous quartic $WW\gamma\gamma$ and
$ZZ\gamma\gamma$ couplings in $\gamma p$ collision at the LHC, Int.
J. Mod. Phys. A 29, 1450148 (2014)

[22]
C. Baldenegro, S. Fichet, G. von Gersdorff and C. Royon,
Probing the anomalous $\gamma\gamma\gamma Z$ coupling at the
LHC with proton tagging, JHEP 06, 142 (2017)

[23]
J.-C. Yang, Y.-C. Guo, C.-X. Yue, and Q. Fu, Constraints on
anomalous quartic gauge couplings via $Z\gamma jj$ production at the
LHC, Phys. Rev. D 104, 035015 (2021) [arXiv:2107.01123].

[24]
O.J.P. Éboli, M.B. Magro, P.G. Mercadante, and S.F. Novaes,
Quartic anomalous couplings in $\gamma\gamma$ colliders,
Phys. Rev. D 52, 15 (1995) [arXiv:hep-ph/9503432].

[25]
S. Atağ and İ. Şahin, Anomalous quartic
$WW\gamma\gamma$ and $ZZ\gamma\gamma$ couplings in $e\gamma$
collision with initial beams and final state polarizations, Phys.
Rev. D 75, 073003 (2007) [arXiv:hep-ph/0703201].

[26]
O.J.P. Éboli and J. K. Mizukoshi, Probing anomalous
quartic couplings in $e\gamma$ and $\gamma\gamma$ colliders, Phys.
Rev. D 64, 075011 (2001) [arXiv:hep-ph/0105238].

[27]
İ. Şahin, Anomalous quartic $WW\gamma\gamma$ and
$WWZ\gamma$ couplings through $W^+W^-Z$ production in $\gamma\gamma$
collisions, J. Phys. G: Nucl. Part. Phys. 36, 075007
(2009) [arXiv:0807.4777].

[25]
S. Atağ and İ. Şahin, Anomalous quartic
$WW\gamma\gamma$ and $ZZ\gamma\gamma$ couplings in $e\gamma$
collision with initial beams and final state polarizations, Phys.
Rev. D 75, 073003 (2007) [arXiv:hep-ph/0703201].

[28]
M. Köksal, V. Ari and A. Senol, Search for anomalous
quartic $ZZ\gamma\gamma$ couplings in photon-photon collisioms,
Adv. High Energy Phys. 2016, 8672391 (2016)

[29]
A. Senol, O. Karadeniz, K.Y. Oyulmaz, C. Helveci, and H.  Denizli,
Sensitivity of anomalous quartic gauge couplings via
$Z\gamma\gamma$ production at future hadron-hadron colliders, Nucl.
Phys. B 980, 115851 (2022) [arXiv:2109.12572].

[30]
A. Senol, H. Denizli, and C. Helveci, Sensitivity of anomalous
quartic gauge couplings via tri-photon production at FCC-hh,

[31]
A. Gutiéerrez-Rodríguez, V. Ari, E. Gurkanli, M. Köksal,
and M.A. Hernández-Ruíz, Future projections on the
anomalous $WW\gamma\gamma$ couplings in hadron-hadron interactions
at the FCC-hh, J. Phys. G 49, 105004 (2022)

[32]
A. Senol, C.O. Karadeniz, K.Y. Oyulmaz, C. Helveci, and H. Denizli,
Sensitivity of anomalous quartic gauge couplings via
$Z\gamma\gamma$ production at future hadron-hadron colliders,
Nucl.Phys. B 980, 115851 (2022) [arXiv:2109.12572].

[33]
O.J.P. Éboli, M.C. Gonzalez-Garcia, and
S.F. Novaes,Quartic anomalous couplings in $e\gamma$
colliders, Nucl. Phys. B 411, 381 (1994)

[34]
E. Gurkanli, V. Ari, A. Gutiérrez-Rodríguez,
M.A. Hernández-Ruíz, and M. Köksal, Sensitivity
physics expected to the measurement of the quartic $WW\gamma\gamma$
couplings at the LHeC and the FCC-he, J. Phys. G: Nucl. Part. Phys.
47, 095006 (2020) [arXiv:2003.06669].

[35]
V. Ari, E. Gurkanli, A.A. Billur, and M. Köksal, Model
independent study for the anomalous quartic $WW\gamma\gamma$
couplings at future electron-proton colliders, Nucl. Phys. B
957 115102 (2020) [arXiv:1812.07187].

[36]
A. Gutiéerrez-Rodríguez, M.A. Hernández-Ruíz,
E. Gurkanli, V. Ari, and M. Köksal, Study on the anomalous
quartic $W^+W^-\gamma\gamma$ couplings of electroweak bosons in
$e^-p$ collisions at the LHeC and the FCC-he, Eur. Phys. J. C
81, 210 (2021) [arXiv:2005.11509].

[37]
E. Gurkanli, Bounds on anomalous quartic $WWZ\gamma$ couplings
in e-p collisions at the FCC-he, J. Phys. G: Nucl. Part. Phys.
50, 015002 (2023) [arXiv:2203.16247].

[38]
W.J. Stirling and A. Werthenbach,Anomalous quartic couplings
in $W^+ W^- \gamma$, $Z^0 Z^0 \gamma$ and $Z^0 \gamma\gamma$
production at present and future $e^+ e^-$ colliders, Eur. Phys. J.
C 14, 103 (2000) [arXiv:hep-ph/9903315].

[39]
G. Perez, M. Sekulla, D. Zeppenfeld, Anomalous quartic gauge
couplings and unitarization for the vector boson scattering process
$pp\rightarrow W^+W^+jjX \rightarrow l^+\nu_l l^+ \nu_l jjX$, Eur.
Phys. J. C 78, 759 (2018) [arXiv:1807.02707].

[40]
Y.-C. Guo, Y.-Y. Wang, J.-C. Yang, and C.-X. Yue, Constraints
on anomalous quartic gauge couplings via $W\gamma jj$ production at
the LHC, Chin. Phys. C 44, 123105 (2020)

[41]
Y.-C. Guo, Y.-Y. Wang, and J.-C. Yang, Constraints on
anomalous quartic gauge couplings by $\gamma\gamma\rightarrow
W^+W^-$ scattering, Nucl. Phys. B 961 115222 (2020)

[42]
H. Amarkhail, S.C. İnan, and A.V. Kisselev, Probing
anomalous $\gamma\gamma\gamma\gamma$ couplings at a future muon
collider, arXiv:2306.03653.

[43]
J.-C. Yang, X.-Y. Han, Z.-B. Qin, T. Li, and Y.-C. Guo,
Measuring the anomalous quartic gauge couplings in the $W^+W^-
\rightarrow W^+W^-$ process at muon collider using artificial neural
networks, JHEP 09, 074 (2022) [arXiv:2204.10034].

[43]
J.-C. Yang, Z.-B. Qing, X.-Y. Han, Y.-C. Guo, and T. Li,
Tri-photon at muon collider: a new process to probe the
anomalous quartic gauge coupling, JHEP 07, 053 (2022)

[44]
Y.-F. Dong, Y.-C. Mao, and J.-C. Yang, Searching for anomalous
quartic gauge couplings at muon colliders using principle component
analysis, Eur. Phys. J. C 83, 555 (2023)

[45]
S. Zhang, J.-C. Yang, and Y.-C. Guo, Using k-means assistant
event selection strategy to study anomalous quartic gauge couplings
at muon colliders, Eur. Phys. J. C 84, 142 (2024)

[46]
R.S. Gupta, Probing quartic gauge boson couplings using
diffractive photon fusion at the LHC, Phys. Rev. D 85,
014006 (2012) [arXiv:1111.3354].

[47]
T. Han, T. Li, and X. Wang, Axion-like particles at high
energy muon colliders – A white paper for Snowmass 2021,

[48]
C.F. Weizsäcker, Ausstrahlung bei Stößen sehr
schneller Elektronen, Z. Phys. 88, 612 (1934).

[49]
E.J. Williams, Nature of the high energy particles of
penetrating radiation and status of ionization and radiation
formulae, Phys. Rev. 45, 729 (1934).

[50]
S.J. Brodsky, Two-photon mechanism of particle production by
high-energy colliding beams, Phys. Rev. D 4, 1532 (1971).

[51]
H. Terazawa, Two-photon processes for particle production at
high energies, Rev. Mod. Phys. 45, 615 (1973).

[52]
V.M. Budnev, I.F. Ginzburg, G.V. Meledin, and V.G. Serbo, The
two photon particle production mechanism. Physical problems.
Applications. Equivalent photon approximation, Phys. Rep.
15, 181 (1975).

[53]
C. Carimalo, P. Kessler, and J. Parisi, Validity of the
equivalent-photon approximation for virtual photon-photon
collisions, Phys. Rev. D 20, 1057 (1979).

[54]
S. Dawson, The effective $W$ approximation, Nucl. Phys. B
249, 42 (1985).

[55]
G.L. Kane, W.W. Repko, and W.B. Rolnick, The effective $W^\pm,
Z^0$ approximation for high energy collisions, Phys. Lett. B
148, 367 (1984).

[56]
R.N. Cahn and S. Dawson, Production of very massive Higgs
bosons, Phys. Lett. B 136, 196 (1984) [Erratum-ibid. B
138, 464 (1984)].

[57]
R.N. Cahn, Production of heavy Higgs bosons: comparisons of
exact and approximate results, Nucl. Phys. B 255, 341
(1985) [Erratum-ibid. B 262, 744 (1985)].

[58]
J. Lindfors, Distribution functions for heavy vector bosons
inside colliding particle beams, Z. Phys. C 28, 427

[59]
J.F. Gunion, J. Kalinowski, and A. Tofighi-Niaki, Exact
calculation of $ff\rightarrow ffWW$ for the charged-current sector
and comparison with the effective-$W$ approximation, Phys. Rev.
Lett. 57, 2351 (1986).

[60]
G. Altarelli, B. Mele, and F. Pitolli, Heavy Higgs production
at future colliders, Nucl. Phys. B 287, 205 (1987).

[61]
J. Lindfors, Luminosity functions for $W^\pm$ and $Z^0$
initiated processes, Z. Phys. C 35, 355 (1987).

[62]
P.W. Johnson, F.I. Olness, and W.-K. Tung,
Effective-vector-boson method for high-energy collisions
Phys. Rev. D 36, 291 (1987).

[63]
I. Kuss and H. Spiesberger, Luminosities for
vector-boson-vector-boson scattering at high energy colliders,
Phys. Rev. D 53, 6078 (1996) [arXiv:hep-ph/9507204].

[64]
R. Ruiz, A. Costantini, F. Maltoni, and O. Mattelaer, The
effective vector boson approximation in high-energy muon
collisions, JHEP 06, 114 (2022) [arXiv:2111.02442].

[65]
G. Jikia and A. Tkabladze, Photon-photon scattering at the
photon linear collider, Phys. Lett. B 323, 453 (1994)

[66]
G.J. Gounaris, P.I. Porfyriadis, and F.M. Renard,
Light-by-light scattering at high energy: a tool to reveal bew
particles, Phys. Lett. B 323, 76 (1994)

[67]
G.J. Gounaris, P.I. Porfyriadis, and F.M. Renard, The
$\gamma\gamma \rightarrow \gamma\gamma$ process in the standard and
SUSY models at high energies, Eur. Phys. J. C 9, 673
(1999) [arXiv:hep-ph/9902230].

[68]
G.J. Gounaris, J. Layssac, P.I. Porfyriadis, and F.M. Renard,
The $\gamma\gamma\rightarrow\gamma Z$ process at high energies
and the search for virtual SUSY effects, Eur. Phys. J. C
10, 499 (1999) [arXiv:hep-ph/9904450].

[69]
G.J. Gounaris, J. Layssac, P.I. Porfyriadis, and F.M. Renard,
The $\gamma\gamma\rightarrow ZZ$ process and the search for
virtual SUSY effects at a $\gamma\gamma$ collider, Eur. Phys. J. C
13, 79 (2000) [arXiv:hep-ph/9909243].

[70]
Y.-J. Zhang and J.-F. Shen, Probing anomalous $tqh$ couplings
via single top production in associated with the Higgs boson at the
HE-LHC and FCC-hh, Eur. Phys. J. C 80, 811 (2020).

[71]
N. Bartosik et al., Simulated Detector Performance at
the Muon Collider, arXiv:2203.07964.

[72]
C. Accettura et al., Towards a Muon Collider,

[73]
M. Jacob and G.C. Wick, On the general theory of collisions
for particlas with spin, Ann. Phys. (N.Y.) 7, 404 (1959);
Ann. Phys. 281, 774 (2000).

[74]
E.P. Wigner, Group theory and its application to the quantum
mechanics of atomic spectra, Academic Press, New York, 1959.