# From Flow to Jamming:
Lattice Gas Automaton Simulations in Granular Materials

Mohamed Gaber  Minerva University, 14 Mint Plaza, Suite 300, San Francisco, CA
94103, USA<EMAIL_ADDRESS>Raquel H. Ribeiro  Janek Kozicki  Faculty of
Applied Physics and Mathematics, Gdańsk University of Technology, 80-233
Gdańsk, Poland  Advanced Materials Center, Gdańsk University of Technology,
80-233 Gdańsk, Poland

###### Abstract

We introduce the first extension of a Lattice Gas Automaton (LGA) model to
accurately replicate observed emergent phenomena in granular materials with a
special focus on previously unexplored jamming transitions by incorporating
gravitational effects, energy dissipation in particle collisions, and wall
friction. We successfully reproduce flow rate evolution, density wave
formation, and jamming transition observed in experiments. We also explore the
critical density at which jamming becomes probable. This research advances our
understanding of granular dynamics and offers insights into the jamming
behavior of granular materials.

††journal: Physical Review Letters

##  1  Introduction

Granular materials, composed of solid particles like sand or grains, challenge
conventional matter classifications [1, 2, 3, 4]. Unlike regular substances,
they mimic solid [5, 6], liquid [2, 4, 7], or gas states [8] depending on the
environment conditions [4, 9, 5], even forming a distinct matter phase (or a
transitional bridge) between the primary phases [3]. Their unique attributes
(such as insensitivity to temperature changes and nonlinear friction [10])
challenge predictive modeling [11], and yet offer the promise of innovative
applications in the energy and pharmaceutics sectors [3].

The boundary between solid- and fluid-like states in granular flow is marked
by the jamming transition, a well-studied shift from a flowing to a static
state [12, 13, 14, 15, 16], triggered by critical control parameters: density,
temperature, and shear stress [15]. This process resembles the glass
transition in amorphous materials, from dynamic to ordered states [15].
Increased density naturally halts dynamics in granular materials, leading to
the static, “jammed” state [14, 15, 16, 17, 18, 19]. Teitel et al. linked
density to viscosity, showing viscosity rises with density until a critical
point, triggering the solid jammed state [16]. Shear stress and temperature
further influence this state—when either falls below a critical threshold, the
reduced mobility leads to the system becoming “jammed.”

Density waves emerge as spatial density variations that reveal flow dynamics
and mechanical traits [5, 20], and are therefore especially challenging to
model. They are triggered by vibrations and shear forces imprinting intricate
patterns and structures [21, 22], and collectively contributing to particle
motion.

This study leverages Lattice Gas Automaton (LGA) simulations to extend prior
experimental work on density waves and jamming transition phenomena. Our novel
contribution lies in the simulation of the jamming process by specifying
appropriate model rules, addressing an unexplored gap in prior work and
triggering potential practical applications [23].

##  2  Model Description

Modeling fluid flow dynamics requires solving the Navier–Stokes equations—a
well-known computationally intensive task [24, 25, 26, 27]. In contrast, the
Lattice Gas Automaton (LGA) model introduced by Frisch et al. [28] offers a
more efficient approach. Unlike traditional methods that require the
definition of macroscopic variables and solving partial differential
equations, LGA uses rules valid at microscopic scales to predict macroscopic
behavior. This bottom-up approach simplifies continuous equations into
discrete rules, rendering LGA an accurate model for fluid simulation [29].
Notably, LGA effectively models granular materials, aligning closely with
experimental observations [22, 30, 31, 32]. In this section, we present a
modified LGA variant tailored for simulating density waves and jamming
transitions in granular materials.

###  2.1  Model Configuration

The LGA model comprises a $2D$ lattice with $L\times L$ hexagonal cells, each
having six neighboring cells. Each cell can be occupied by one or more
particles, empty, or act as a wall. The particles possess velocity vectors
indicating their movement directions. The velocity vector consists of six
binary elements, with only one element set to $1$ representing movement in the
direction of that element (Fig. 1).

Figure 1:  The cells are tagged with the particle’s velocity. A value of $0$
indicates a particle at rest, while values $1-6$ represent particles in
movement.

###  2.2  Adding Gravitational Effects

We follow the approach by Kozicki et al. [30] and introduce a parameter $g$,
ranging from $0$ to $1$, which represents the probability of particles
changing their velocity vector direction toward gravity at each iteration. For
example, if the particle’s direction is $0$ (Fig. 1), its new direction value
will be $4$ with probability $g$. If the particle’s direction were $2$, its
new direction would be $3$ with probability $g$. This modification allows for
parabolic behavior in individual particle velocities, consistent with the
findings of Kozicki et al. [30].

###  2.3  Adding Energy Dissipation

In LGA models, particle interactions are typically assumed to be elastic
collisions, which works well to model ideal gases but not granular materials
[33]. Several research experiments have found that granular particle collision
can be elastic or show energy dissipation due to particle crushability or
grain roughness that varies from one particle to another [22, 23, 31]. To
address this possibility, we adapt the collision model proposed by Herrmann et
al. [6] and introduce an energy dissipation parameter, $p$, ranging from $0$
to $1$—a value of $0$ corresponds to perfectly elastic collisions, while a
value of $1$ signals fully inelastic collisions with complete energy
dissipation. Since the roughness of particles is not known in advance, an
additional parameter $p$ accounts for probabilistic energy dissipation in
collisions, following the rules depicted in Figs. 2-4. Additionally, we
implement mirror deflection to simulate the system’s behavior when particles
collide with walls.

Figure 2:  Head-on particle collisions result in one of three deflection
outcomes, each with a 1/3 probability. Figure 3:  Collisions between particles
at an angle yield different outcomes depending on the energy dissipation
parameter, $p$. For $p=0$, collision conservation laws govern energy and
momentum. With $p>0$, one of the particles loses velocity through energy
dissipation. Note that given the hexagonal grid, the collision angles are
equal to the rebound angles. Figure 4:  A three-particle collision event can
lead to different outcomes depending on the value of $p$. Collisions conserve
energy and momentum in the absence of energy dissipation ($p=0$). However,
when p $>$ 0, there are three possible outcomes, each with a probability of
$p/3$. In each outcome, one of the particles experiences energy dissipation.
Note that the hexagonal lattice configuration limits the direction of movement
to the given array of possibilities.

###  2.4  The Onset of the Jamming Transition

Experimental studies have consistently demonstrated the critical role of the
formation of an arch structure at the narrow opening of a hopper in initiating
jamming transitions [34]. This arch structure results from the friction
between the particles and the hopper walls, creating a barrier that obstructs
the material flow.

To model friction within a lattice gas system, we have carefully designed a
friction rule that governs the interaction between particles and their
surrounding environment. Firstly, the walls exert a force in the opposite
direction to the movement of particles. Secondly, particles can transmit
frictional forces when they come into contact. However, for the friction force
to take effect, particles must be compressed from both sides, creating a state
of compression. Figure 5a shows particles arranging themselves into an arch-
like structure, with arrows indicating the upward-pointing frictional forces
due to the compression of particles toward the walls. The friction force is
transmitted upwards through the compressed particles, ultimately leading to
the structure acting as a barrier, impeding the flow of particles.

Figure 5:  Illustration of an arch-like structure formation as a result of the
upward frictional forces caused by the compression from both sides of the
walls in 5a. In contrast, 5b illustrates a scenario where the particle does
not experience friction from the wall due to the absence of compression from
the left side.

##  3  Simulation Results

###  3.1  Entropy of the LGA Model

Table I: Parameter Values for Each Simulation | $g$ | $p$ | Grid H$\times$W | Density
---|---|---|---|---
Particles in a box (Figs. 6–7) | 0 | 0 | 100$\times$100 | 10%
Density Waves (Figs. 8–9) | 0.2 | 0.1 | 200$\times$50 | 70%
Narrow Pipe (Figs. 10–11) | 0.2 | 0.1 | 200$\times$10 | 85%
Jamming (Figs. 12–13) | 0.2 | 0.1 | 200$\times$20 | 5%–100%

In the absence of gravitational and frictional effects, the model is suitable
for simulating standard fluid flows and is expected to adhere to the
principles of thermodynamics. Notably, the second law of thermodynamics should
apply, indicating that the system’s entropy, when out of equilibrium, should
increase with time. This implies that the number of possible microstates
available to the system should progressively grow. We initiated the simulation
with a low entropy state to verify this, as shown in Fig. 6a. The parameters
used for simulation in this section and for the following sections are listed
in Tab. I.

Figure 6:  The time evolution of the particles after being confined to a small
region in a box. A video of the simulation can be consulted here.

The system, depicted in Figure 6, starts with a low entropy state, where all
particles are placed in one corner of the box. As the system evolves, the
particles will collide and occupy a larger portion of the available volume
(Figs 6 b-f). Consequently, the number of possible microstates at equilibrium
will significantly exceed the restricted microstates confined to the $20\times
20$ region, increasing entropy. To quantify the entropy of this system, we
employ the Shannon entropy, defined as follows:

$H(x)=-\sum_{i=1}^{n}P(x_{i})\log(P(x_{i}))$ (1)

In the context of information entropy, $x$ represents a random variable, while
$x_{i}$ represents the possible outcomes associated with that random variable.
The term $P(x_{i})$ denotes the probability of a specific outcome $x_{i}$. In
the present study, the equation is applied by making $x$ akin to a random
indicator variable: $1$, indicating the cell is occupied by a particle, or
$0$, showing the cell is empty. To compute the probability $P(x_{i})$, we
conduct 50 distinct simulations, each spanning 500 steps. This generates an
ensemble of 50 microstate configurations for each iteration. The probability
$P(x_{i})$ for a given cell $x_{i}$ in a particular iteration is then computed
as the average occupancy frequency of that cell over the 50 distinct
simulations. This average is determined by dividing the number of times cell
$x_{i}$ is occupied by the particle by the total number of simulations (50).
After computing the probabilities $P(x_{i})$ for all cells, they are summed
using equation (1). This procedure is repeated for each iteration in the
model’s evolution.

Figure 7: The entropy over time for the system initialized with a non-
equilibrium state (Fig. 6).

In Fig. 7, we show the time evolution of entropy in the system that
accompanies the stages depicted in Fig. 6. A transient decline in entropy is
allowed—rapid entropy shifts can occur in out of equilibrium states. We
further confirm the finding by Tribel and Boon that the monotonous increase in
LGA entropy is not guaranteed [35]. This fluctuating behavior persists
equilibrium is reached, and a homogeneous distribution of particles across the
entire box and a high entropy ensue.

###  3.2  Density Waves in a hopper flow

Behringer et al. initially documented density waves in 1989 [5]. Employing a
Plexiglas hopper with adjustable width, the team captured gravity-driven sand
flow using X-ray imaging through digital fluoroscopy. They examined two sand
types: rough (resulting in inelastic collisions due to crushability and shear
strain) and smooth (with predominant elastic collisions). These experiments
established insights into density wave formation and behavior in granular
materials.

Figure 8:  Snapshots of the sand flow in a hopper taken from the model
simulation. The three consecutive snapshots, separated by a few iterations,
illustrate the temporal evolution of the particle discharge process. A video
of the simulation can be consulted here.

We successfully reproduced the phenomenon of density waves observed in the
experimental study through simulations. As depicted in Fig. 8, the flow forms
a central low-density region positioned directly above the outlet. This
central region serves as the primary pathway for the flow, while the adjacent
regions near the walls exhibit limited movement. We also observed the
emergence of low-density regions propagating upwards between the central
region and the walls, contrary to the downward sand flow. These propagating
density waves persist and are followed by additional waves until they reach
the system’s top, eventually collapsing. This simulation-based observation
aligns with the previous experimental findings on the propagation of density
waves [5].

The flow rate evolution at the hopper outlet, measured by counting the number
of particles leaving the outlet every iteration, is depicted in Fig. 9. It is
characterized by a consistent flow rate during the particle discharge,
followed by a non-linear reduction towards the end that can be attributed to
the decreasing particle density towards the end. The observed flow rate
behavior aligns with experimental findings [9, 36].

Figure 9:  Flow rate of the granular materials measured leaving the hopper
outlet. The data is averaged over ten iterations, showing a roughly constant
flow rate (indicated by the black horizontal line) until a decrease occurs
towards the end.

###  3.3  Density Waves in narrow pipes

In 1994, Peng et al. provided a distinctive perspective on density wave
propagation, relevant to industrial applications [37]. Their study, utilizing
narrow pipes instead of conventional hoppers, examined smooth and rough sand.
Density wave propagation exclusively occurred with rough sand. The underlying
reason is the inelastic collisions experienced by rough sand particles during
descent. These collisions lead to energy dissipation, causing certain
particles to lose velocity. Consequently, these decelerated particles lag
behind, triggering successive collisions and more energy loss. This cascading
effect resembles shock waves in traffic, where one car’s abrupt stop triggers
a chain reaction of deceleration among preceding vehicles.

Figure 10:  The set of figures, from left to right, shows the evolution of the
flow in a narrow pipe and under periodic boundary conditions. Note the
formation of shock waves. A link to the simulation video can be consulted
here.

Our model successfully replicated the results when we adjusted the simulation
settings to a narrow pipe and replaced the hopper outlet with periodic
boundary conditions. This adjustment allowed particles that flowed through the
bottom of the pipe to re-enter the system through the top, ensuring a
continuous flow. The simulation results are presented in Fig. 10.

Fig. 10 illustrates a notable phenomenon—the emergence of a shock wave of high
density and propagation in an upward direction. This behavior is prominently
observed in narrow pipes, as the occurrence of inelastic collisions results in
the halting of numerous particles from flowing, consequently leading to the
formation of the density wave (or shock wave). Figure 11 provides a
comprehensive visualization of the temporal propagation of density waves,
following a methodology similar to that employed by Peng et al. [37]. Multiple
snapshots were taken from the narrow pipe at different time intervals to
generate the figure. Each snapshot involved dividing the pipe into small
vertical segments and calculating the density within each segment. The density
values were then represented using a grayscale scheme, with darker regions
indicating higher densities. These snapshots were sequentially arranged from
left to right, creating the impression of multiple narrow pipes aligned
horizontally, with time progressing from left to right.

The graph serves as a clear demonstration of the formation and movement of
density waves. The initial state, depicted on the left side of the graph,
portrays a semi-uniform density distribution. As time elapses from left to
right, density waves’ emergence and upward propagation become increasingly
evident. It is worth noting that once the density wave reaches the top of the
graph, it reappears at the bottom due to the periodic boundary representation
implemented in the simulation. After $10,000$ iterations near the right side
of the figure, the number of density waves reduces as the voids merge into
larger voids.

Figure 11:  Sequential snapshots illustrating density wave propagation in a
narrow pipe over time. Each snapshot represents a vertical segment of the
pipe, with density depicted using grayscale colors (darker shades indicating
higher densities).

###  3.4  Jamming Transition

As mentioned above, the jamming transition occurs when an arch structure forms
at the narrow opening of a hopper, obstructing the particle flow. This
structure originates from the friction between the particles and the hopper
walls, effectively blocking the granular flow. Multiple experimental studies,
such as the one by Kiwing et al. [34] have explored the underlying mechanisms
governing this process highlighting the role of the specific hopper structure.
The emergence of this configuration is stochastic in nature and dependent on
particle’s phase space.

Furthermore, these investigations shed light on the role of density in jamming
transition. Higher densities are associated with increased jamming transition,
as they entail more particles and a more extensive set of possible
microstates, which significantly enhances the probability that one of the
configurations will lead to the formation of an arch and subsequent jamming.
Consequently, the jamming transition is not solely contingent upon reaching a
critical density value but becomes more likely above this critical density
threshold. This distinction sets the jamming transition apart from traditional
phase transitions, such as the water-ice transition, typically occurring at
specific critical values of control parameters like temperature or pressure.
Our simulations identified a critical density threshold above which jamming
transition is more likely to occur, as illustrated in Fig. 13.

Figure 12:  Flow evolution in a narrow-outlet hopper. Fig. 12a is initialized
with $15\%$ packing density for which no jamming occurs, while Fig. 12b is
initialized with $30\%$ packing density for which jamming occurs. A video of
the simulation can be consulted here.

The simulation depicted in Fig. 12 involved initializing a narrow hopper with
varying packing densities of particles, followed by measuring the resulting
flow rate. To achieve this, a predetermined number of particles corresponding
to each density was randomly distributed within the hopper. The simulation was
then executed for 5000 iterations. At each density, the flow rate was
quantified as the average number of particles exiting through the hopper
outlet per unit of time. To ensure the robustness of the results, this process
was repeated $30$ times for each density. The reported flow rate in Fig. 13 is
the average flow rate from the $30$ simulations for each density value, with
error bars indicating the 95% confidence intervals.

As expected, Fig. 13 demonstrates a persistent linear relationship between the
particle density and the flow rate until reaching a density of approximately
$20$%, beyond which the system contains enough particles to increase the
probability of arch formation, reducing the flow rate. As the density
approaches $100$%, the box becomes densely packed with particles, resulting in
a near-certain arch formation and a complete halt in particle flow. It should
be noted that, statistically, the flow rate will not reach zero as a few
particles may still exit the system before the arch forms.

Figure 13:  The graph illustrates the flow rate measured at different
densities and the onset of the jamming transition. The model was run for 30
different simulations at each density, each running for 5000 iterations. The
average flow rate was reported with 95% confidence intervals.

We also explored jamming and flowing outcomes for multiple particle densities.
Jamming is relatively rare for densities below $20$%, with infrequent flowing
events for densities above $20$%. This probabilistic nature of jamming
indicates that higher densities increase the likelihood of its occurrence,
again corroborating previous experimental findings [34].

##  4  Conclusion

This study reports on the first successful simulation of granular materials
using an LGA framework. The model accurately replicates critical empirical
observations by incorporating gravitational effects, energy dissipation in
particle collisions, and frictional interactions.

The main novelty of this research lies in the reproduction of the jamming
transition, previously only observed experimentally. We recognized the
probabilistic nature of the critical density threshold for arch formation
demonstrating how it is more likely at higher densities. Our tailored LGA
model also succeeds in reproducing the flow rate evolution and density wave
formation at the hopper outlet. The simulated non-linear reduction in flow
rate towards the end of the discharge process suggests the role of density
fluctuations in arch formations.

These findings further our understanding of granular dynamics, providing
valuable insights into the complex behavior of granular materials. Based on
these results, one natural direction to follow up would be to delve into how
the jamming transition in granular materials might affect the propagation and
behavior of density waves. Investigating the interplay between the jamming
transition and the emergence of density waves could unveil intriguing
connections between these two phenomena, potentially shedding light on the
underlying mechanisms that govern both processes. This venue of research could
contribute to a more comprehensive understanding of the complex behaviors
exhibited by granular materials near the jamming transition and their
implications for various practical applications.

ORCID
Mohamed Gaber: 0009-0004-9562-1244
Raquel H. Ribeiro: 0000-0002-5434-252X
Janek Kozicki: 0000-0002-8427-7263

## References

  * [1] T. Shinbrot, H. M. Jaeger, P. B. Umbanhowar, Does the granular matter?, Proceedings of the National Academy of Sciences 97 (24) (2000) 12959\.
URL https://www.pnas.org/content/97/24/12959#ref-6

  * [2] S. B. Savage, Flow of granular materials, Theoretical and Applied Mechanics (1989) 241–266doi:10.1016/b978-0-444-87302-6.50022-3.
  * [3] F. F. Chung, S. Liaw, M. C. Ho, Energy and phase transition in a horizontally vibrating granular system, Granular Matter 12 (4) (2010) 369–374.  doi:10.1007/s10035-010-0184-9.
  * [4] F. in Granular Materials, Built on Sand, The MIT Press, 2020.  doi:10.7551/mitpress/12202.003.0013.
  * [5] G. W. Fagert, R. P. Behringer, T., G. A. Johnson, Pattern formation in flowing sand, Physical Review Letters 62 (24) (1989) 2825–2828.  doi:10.1103/physrevlett.62.2825.
  * [6] H. Herrmann, On the shape of the tail of a two dimensional sand pile, Physical Review Letters (1996).
  * [7] L. P. Kadanoff, S. R. Nagel, L. Wu, S. Zhou, Scaling and universality in avalanches, Physical Review A 39 (12) (1989) 6524–6537.  doi:10.1103/physreva.39.6524.
  * [8] G. Mesri, B. Vardhanabhuti, Compression of granular materials, Canadian Geotechnical Journal 46 (4) (2009) 369–392.  doi:10.1139/t08-123.
  * [9] R. P. Behringer, T. Fagert, G. A. Johnson, Pattern formation and time-dependence in flowing sand, in: Two Phase Flows and Waves The IMA Volumes in Mathematics and Its Applications, Springer, 1990, pp. 1–29.  doi:10.1007/978-1-4613-9022-0_1.
  * [10] Z. X. Yang, J. Yang, L. Z. Wang, On the influence of inter-particle friction and dilatancy in granular materials: A numerical analysis, Granular Matter 14 (3) (2012) 433–447.  doi:10.1007/s10035-012-0348-x.
  * [11] H. M. Jaeger, C. Liu, S. R. Nagel, Relaxation at the angle of repose, Physical Review Letters 62 (1) (1989) 40–43.  doi:10.1103/physrevlett.62.40.
  * [12] J. P. W. M. E. Cates, Jamming and static stress transmission in granular materials, in: Jamming and Rheology, Pubmed, 2001, pp. 149–160.  doi:10.1201/9781482268171-14.
  * [13] M. Sadjadpour, C. S. Campbell, Granular chute flow regimes: mass flowrates, flowrate limits and clogging, Powder Technology 186 (2) (2008) 96–107.  doi:10.1016/j.powtec.2008.01.028.
  * [14] V. Trappe, V. Prasad, L. Cipelletti, P. N. Segre, D. A. Weitz, Jamming phase diagram for attractive particles, Nature 411 (6839) (2001) 772–775.  doi:10.1038/35081021.
  * [15] C. S. O’Hern, L. E. Silbert, A. J. Liu, S. R. Nagel, Jamming at zero temperature and zero applied stress: The epitome of disorder, Physical Review E 68 (1) (2003).  doi:10.1103/physreve.68.011306.
  * [16] S. Teitel, P. Olsson, Critical scaling of shear viscosity at the jamming transition, Physical Review Letters 99 (17) (2007).  doi:10.1103/physrevlett.99.178001.
  * [17] I. Pajic-Lijakovic, M. Milivojevic, Viscoelasticity and cell jamming state transition, The European Physical Journal Plus (2021).  doi:10.1101/2021.03.19.436195.
  * [18] M. Jean, Numerical simulation of granular materials, in: Micromechanics of Granular Materials, ISTE Ltd, 2009, pp. 149–315.  doi:10.1002/9780470611616.ch4.
  * [19] G. Giraldi, Animation of gas-liquid systems through lattice gas cellular automata and smoothed particle hydrodynamics, ArXiv (2005).
  * [20] T. Riethmüller, L. Schimansky-Geier, D. Rosenkranz, T. Pöschel, Langevin equation approach to granular flow in a narrow pipe, Journal of Statistical Physics 86 (1-2) (1997) 421–430.  doi:10.1007/bf02180213.
  * [21] L. E. Silbert, J. W. Landry, G. S. Grest, Granular flow down a rough inclined plane: Transition between thin and thick piles, Physics of Fluids 15 (1) (2003) 1–11.  doi:10.1063/1.1521719.
  * [22] M. El Shamy, C. Denissen, Microscale characterization of energy dissipation mechanisms in liquefiable granular soils, Computers and Geotechnics 61 (2014) 347–363.  doi:10.1016/j.compgeo.2014.06.007.
  * [23] R. Mukwiri, Energy dissipation in granular materials in triaxial tests, Materials Today: Proceedings 4 (2) (2017) 3840–3847.  doi:10.1016/j.matpr.2017.03.048.
  * [24] S. F. Judice, Lattice Gas Cellular Automata for Fluid Simulation, Springer, 2018\.
  * [25] K. Matsushima, S. Takanashi, Naiver-stokes simulation of transonic flows about a space-plane, in: 12th Applied Aerodynamics Conference, 1994, pp. 132–164.  doi:10.2514/6.1994-1864.
  * [26] Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element, Dover Publications, 1998.
  * [27] M. Perić, Finite-Volume Methods for Navier-Stokes Equations, Springer, n.d.
  * [28] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice-gas automata for the navier-stokes equation, Physical Review Letters 56 (14) (1986) 1505–1508.  doi:10.1103/physrevlett.56.150.
  * [29] B. Chopard, Cellular automata modeling of physical systems, in: Cellular Automata, Cambridge University Press, 2009, pp. 657–689.  doi:10.1007/978-1-4939-8700-9_57.
  * [30] J. Kozicki, Application of a cellular automaton to simulations of granular flow in silos, Granular Matter 7 (1) (2005) 45–54.  doi:10.1007/s10035-004-0151-3.
  * [31] Z. Zhao, C. Liu, B. Brogliato, Energy dissipation and dispersion effects in granular media, Physical Review E 78 (3) (2008).  doi:10.1103/physreve.78.031307.
  * [32] A. D. Fitt, P. Wilmott, Cellular-automaton model for segregation of a two-species granular flow, Physical Review A 45 (4) (1992) 2383–2388.  doi:10.1103/physreva.45.2383.
  * [33] A. S. Kim, H.-J. Kim, Dissipative dynamics of granular materials, in: Granular Materials, InTech Open Science, 2017, pp. 34–48.  doi:10.5772/intechopen.69196.
  * [34] T. Kiwing, K. A. Dahmen, Jamming patterns in granular flow through a horizontal pipe, Physical Review Letters 94 (15) (2005).  doi:10.1103/physrevlett.94.158002.
  * [35] O. Tribel, J. P. Boon, Entropy and correlations in lattice-gas automata without detailed balance, International Journal of Modern Physics C 08 (04) (1997) 641–652.  doi:10.1142/s0129183197000552.
URL http://dx.doi.org/10.1142/S0129183197000552

  * [36] G. W. Baxter, Cellular automata models of granular flow, Experimental and Computational Techniques in Soft Condensed Matter Physics (1991) 209–229doi:10.1017/cbo9780511760549.008.
  * [37] G. Peng, Herrmann, Density waves of granular flow in a pipe using lattice-gas automata, Physical Review E 49 (3) (1994).


aainstitutetext: LUPM, CNRS, Université Montpellier Place Eugene Bataillon,
F-34095 Montpellier, Cedex 5, Franceccinstitutetext: Physik-Institut,
Universität Zürich, Winterthurerstrasse 190, CH–8057 Zürich,
Switzerlandddinstitutetext: Paul Scherrer Institut, CH–5232 Villigen PSI,
Switzerland

# A viable $L_{e}-L_{\mu}$ model with $\mu\to e$ violation

Marco Ardu c,d Fiona Kirk<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We extend the Standard Model gauge group by $U(1)_{L_{e}-L_{\mu}}$ and
introduce two scalars, a doublet and a singlet, that are charged under this
new group and have lepton flavour violating couplings. Since in this model
$\mu\to e$ processes can only be mediated by $\mu\to\tau\times\tau\to e$
interactions, bounds from $\mu\to e$ transitions can be avoided while allowing
for accessible new physics. We consider a $Z^{\prime}$ boson with a mass of
$M_{Z^{\prime}}\simeq 10~{}$GeV and a gauge coupling $g^{\prime}\simeq
10^{-4}$, which is in reach of Belle-II. Neutrino masses and mixing angles can
also be accounted for if sterile neutrinos are added to the spectrum.

††preprint: PSI-PR-22-12 ZU-TH-17/22

##  1 Introduction

The success of the Standard Model of particle physics (SM) in describing the
experimental data suggests that New Physics (NP) is either heavy or weakly
coupled, or that it preferentially interacts with those sectors of the SM that
are least constrained by experiment (e.g. the third generation of matter
fields, in the case of Minimal Flavour Violation DAmbrosio:2002vsn ).

Processes that are suppressed or forbidden in the SM, such as lepton flavour
violation (LFV) in the charged sector, are smoking gun signals of NP and
constitute valuable probes for Beyond Standard Model scenarios. Charged LFV
has not yet been observed but is to be expected, since the discovery of
neutrino masses and neutrino oscillations represents direct evidence for
lepton flavour changing transitions. The current upper limits on the rates of
rare muon processes $\mathrm{Br}(\mu\to e\gamma)<4\times 10^{-13}$ MEG:2016leq
, $\mathrm{Br}(\mu\to 3e)<10^{-12}$ SINDRUM:1987nra  and $\mathrm{Br}(\mu A\to
eA)<7\times 10^{-13}$ ($\mu\to e$ conversion in nuclei) SINDRUMII:2006dvw
provide the most stringent constraints in models that allow for $\mu\to e$
transitions. The next generation of $\mu\to e$ experiments promise an
improvement in sensitivity of up to four order of magnitude MEGII:2018kmf ;
Blondel:2013ia ; COMET:2009qeh . The bounds on LFV processes involving $\tau$
leptons are less constraining, with
$\mathrm{Br}(\tau\to\ell\gamma),\mathrm{Br}(\tau\to
3\ell)\lesssim\text{few}\times 10^{-8},\ell=e,\mu$ BaBar:2009hkt ;
Hayasaka:2010np . Here the experimental sensitivities are expected to improve
by a factor $\sim 10$ in the near future Belle-II:2018jsg .

A particularly simple way to extend the SM is to enlarge the SM gauge group by
a new spontaneously broken abelian $U(1)^{\prime}$ group. This leads to new
interactions mediated by a neutral massive vector boson, commonly referred to
as a $Z^{\prime}$ gauge boson. $U(1)^{\prime}$ groups may be remnants of
larger non-abelian groups, such as $SO(10)$ and $E_{6}$, which feature in
certain Grand Unified Theories (GUTs) Langacker:2008yv  or they can be
considered as standalone extensions, as in the case of
$U(1)_{L_{e}-L_{\mu}},U(1)_{L_{\mu}-L_{\tau}},U(1)_{L_{e}-L_{\tau}}$ He:1991qd
or $U(1)_{B-L}$ CARLSON1987378 , which are anomaly-free global symmetries of
the SM.

LFV $Z^{\prime}$ couplings are introduced in such models if the mass and gauge
eigenstates of the leptons are misaligned. The stringent bounds from $\mu\to
e$ transitions Langacker:2000ju ; Murakami:2001cs ; Chiang:2011cv  can be
avoided in $U(1)^{\prime}$ models that allow only for
$\tau\leftrightarrow\ell$, flavour changes, $\ell=e,\mu$ Heeck:2011wj ;
Heeck:2014qea ; Heeck:2016xkh , which are less constrained. However, if both
$\tau\leftrightarrow e$ and $\tau\leftrightarrow\mu$ couplings are present,
there is no symmetry forbidding $\mu\leftrightarrow e$ processes mediated by
the product of $\mu\to\tau$ and $\tau\to e$ interactions. Considering that the
future sensitivities of $\mu\leftrightarrow e$ and $\tau\leftrightarrow\ell$
searches approximately satisfy the inequality

$\mathrm{Br}(\mu\to e)\lesssim\mathrm{Br}(\tau\to\mu)\mathrm{Br}(\tau\to e),$
(1)

$\mu\to e$ observables can probe products of $\mu\leftrightarrow\tau$ and
$\tau\leftrightarrow e$ couplings which are beyond the reach of direct
$\tau\to\ell$ searches. If the NP is heavy, model independent
$\mu\to\tau\times\tau\to e$ contributions can be calculated in the Standard
Model Effective Field Theory (SMEFT) Ardu:2022pzk . Since these correspond to
the combination of two dimension 6 operators, they only arise at dimension 8.
In models, the combination of $\mu\to\tau$ and $\tau\to e$ can be larger than
dimension 8 effects, making it more accessible to $\mu\to e$ experiments.

In this article, we propose a $Z^{\prime}$ model where $\mu\to e$ transitions
are mediated by the product of $\mu\leftrightarrow\tau$ and
$\tau\leftrightarrow e$ couplings only. In this way we suppress $\mu\to e$
rates enough to respect the current upper limits, while being in reach of
future $\mu\to e$ experiments. We gauge $L_{e}-L_{\mu}$111In models with
gauged $L_{\mu}-L_{\tau}$ or $L_{e}-L_{\tau}$, where both
$\mu\leftrightarrow\tau$ and $e\leftrightarrow\tau$ couplings are present,
gauge invariance does not forbid $\mu\leftrightarrow e$ couplings., aiming for
a feebly-coupled $Z^{\prime}$ with a mass below the electroweak scale, and
extend the scalar sector by an extra doublet and a new singlet, both of which
are charged under the new $U(1)^{\prime}$. At some unknown high energy, the
vacuum expectation value (VEV) of the singlet, $v_{S}$, breaks
$U(1)_{{L_{e}}-L_{\mu}}$, gives the $Z^{\prime}$ boson a mass and generates
Majorana masses for the sterile neutrinos that are charged under
$U(1)_{{L_{e}}-L_{\mu}}$, setting the stage for a type I seesaw neutrino mass
model. LFV is introduced via the Yukawa interactions with the new scalars, and
the $Z^{\prime}$ boson receives flavour changing couplings when the
$U(1)^{\prime}$ charged doublet acquires a non-vanishing vacuum expectation
value (VEV).

A large region of the model’s parameter space is in reach of future
$Z^{\prime}$ and LFV searches, and will be tested in the upcoming years.

This article is organized as follows: In Section 2 we present the particle
content of the model and define our notation. In Section 3 the phenomenology
is studied. We briefly review the constraints on non-SM gauge interactions of
electrons and muons and give the rates of $\tau\to\ell$ and $\mu\to e$
processes.

##  2 The Model

We extend the Standard Model gauge group $SU(3)_{c}\otimes SU(2)_{L}\otimes
U(1)_{Y}$ by the abelian anomaly-free $U(1)_{L_{1}-L_{2}}\equiv U(1)^{\prime}$
and consider the particle content summarised in Table 1.

|  |  | $U(1)^{\prime}$
---|---|---|---
$L_{1}$ (1, 2, $-\frac{1}{2}$) |  $e_{1}$ (1, 1, -1) |  $N_{1}$ (1, 1, 0) | +1
$L_{2}$ (1, 2, $-\frac{1}{2}$) |  $e_{2}$ (1, 1, -1) |  $N_{2}$ (1, 1, 0) |  $-$1
$L_{3}$ (1, 2, $-\frac{1}{2}$) |  $e_{3}$(1, 1, -1) |  $N_{3}$ (1, 1, 0) | 0
$H$ (1, 2, $\frac{1}{2}$) |  |  | 0
$\phi$ (1, 2, $\frac{1}{2}$) |  |  |  $-$1
$S$ (1, 1, 0) |  |  |  $-$1
Table 1: Field content of the model. In parentheses we give the
representations under the SM gauge group ($SU(3)_{c}$, $SU(2)_{L}$,
$U(1)_{Y}$), in the last column the charges under $U(1)^{\prime}$.

Then the most generic renormalisable Lagrangian is

$\mathcal{L}=\mathcal{L}_{\rm kin}-\mathcal{L}_{\rm Yuk}-V(H,\phi,S)$ (2)

where

$\displaystyle\mathcal{L}_{\rm kin}$
$\displaystyle=-\frac{1}{4}B^{\alpha\beta}B_{\alpha\beta}-\frac{1}{4}W^{a\alpha\beta}W^{a}_{\alpha\beta}-\frac{1}{4}B^{\prime\alpha\beta}B^{\prime}_{\alpha\beta}+$
$\displaystyle+\sum_{\Phi=H,\phi,S}(D_{\alpha}\Phi)(D^{\alpha}\Phi)^{\dagger}+\sum_{\psi}i\overline{\psi}\not{D}\psi,$
(3) $\displaystyle\mathcal{L}_{\rm Yuk}$
$\displaystyle=\sum_{i\in\\{1,2,3\\}}y_{ii}\bar{L}_{i}e_{i}H+y^{\nu}_{ii}\bar{L}_{i}N_{i}\tilde{H}+\mathrm{h.c.}$
$\displaystyle+y_{31}\bar{L}_{3}e_{1}\phi+y_{23}\bar{L}_{2}e_{3}\phi+y_{13}^{\nu}\bar{L}_{1}N_{3}\tilde{\phi}+y_{32}^{\nu}\bar{L}_{3}N_{2}\tilde{\phi}+$
$\displaystyle+\frac{1}{2}y^{\nu}_{13}\bar{N}_{1}^{c}N_{3}S+\frac{1}{2}M^{N}_{33}\bar{N}_{3}^{c}N_{3}+\frac{1}{2}M^{N}_{12}\bar{N}_{1}^{c}N_{2}+\mathrm{h.c.},$
(4) $\displaystyle V(H,\phi,S)$
$\displaystyle=m_{H}^{2}H^{\dagger}H+m_{\phi}^{2}\phi^{\dagger}\phi-
m_{S}^{2}S^{\dagger}S$ $\displaystyle-m_{\phi
HS}\left(\left(\phi^{\dagger}H\right)S+S^{\dagger}\left(H^{\dagger}\phi\right)\right)$
$\displaystyle+\frac{1}{2}\lambda_{H}\left(H^{\dagger}H\right)^{2}+\frac{1}{2}\lambda_{\phi}\left(\phi^{\dagger}\phi\right)^{2}+\frac{1}{2}\lambda_{S}\left(S^{\dagger}S\right)^{2}$
$\displaystyle+\lambda_{H\phi}\left(H^{\dagger}H\right)\left(\phi^{\dagger}\phi\right)+\tilde{\lambda}_{H\phi}\left(H^{\dagger}\phi\right)\left(\phi^{\dagger}H\right)$
$\displaystyle+\lambda_{HS}\left(H^{\dagger}H\right)\left(S^{\dagger}S\right)+\lambda_{\phi
S}\left(\phi^{\dagger}\phi\right)\left(S^{\dagger}S\right)\,.$ (5)

Disregarding $SU(3)_{c}$, we use the following conventions for the covariant
derivative

$D_{\mu}=\partial+ig_{1}YB_{\mu}+ig_{2}\frac{\tau^{a}}{2}W_{\mu}^{a}+ig^{\prime}Q^{\prime}B_{\mu}^{\prime}\,,$
(6)

where $g_{1},g_{2}$ and $g^{\prime}$ are the gauge couplings of
$U(1)_{Y},SU(2)_{L}$ and $U(1)^{\prime}$, respectively, while $\tau^{a}$ are
the Pauli matrices. The label $\psi$ in Eq. (3) runs over all fermions of the
model. We identify the gauge eigenstates $1,2,3$ by their diagonal Yukawa
couplings with the $U(1)^{\prime}$ neutral doublet $H$.

$U(1)^{\prime}$ gauge invariance allows for off-diagonal Yukawa interactions
of type $2\leftrightarrow 3$ and $1\leftrightarrow 3$ with the doublet $\phi$,
but forbids interactions of type $1\leftrightarrow 2$ among the charged
leptons. The off-diagonal Yukawas $y_{31},y_{23}$ are the only parameters of
the charged lepton sector that introduce LFV. As we discuss in Section 2.1,
the mass eigenstates $e,\mu,\tau$ are nearly aligned with the interaction
eigenstates $1,2,3$, and as a result $\mu\to e$ transitions are controlled by
the product $y_{31}\times y_{23}$. Due to the approximate alignment of gauge
eigenstates and mass eigenstates, we will use the term _lepton flavour_ in
both bases and refer to $U(1)_{L_{1}-L_{2}}$ as $U(1)_{L_{e}-L_{\mu}}$.

The scalar potential parameters can be such that all scalars acquire VEVs,
$\expectationvalue{H}=v_{H}/\sqrt{2}$,
$\expectationvalue{\phi}=v_{\phi}/\sqrt{2}$,
$\expectationvalue{S}=v_{S}/\sqrt{2}$. $v_{H}$ and $v_{\phi}$ spontaneously
break the electroweak gauge symmetry and thus must satisfy

$\sqrt{v_{H}^{2}+v_{\phi}^{2}}=v\qquad\text{with}\ v=246\ {\rm GeV}\,,$ (7)

whereas either $v_{S}$ or $v_{\phi}$ can be the VEV that breaks
$U(1)^{\prime}$ and provides the dominant contribution to the $Z^{\prime}$
mass. Since we are gauging the lepton flavor difference $L_{e}-L_{\mu}$,
resulting in a $Z^{\prime}$ boson that couples to electrons and muons, several
experiments can constrain the gauge coupling $g^{\prime}$. Values larger than
$g^{\prime}\sim\text{few}\times 10^{-4}$ are excluded in a vast region of the
$Z^{\prime}$ mass vs. $g^{\prime}$ coupling plane. Since, as discussed in
Section 3, we are aiming at a $Z^{\prime}$ mass in the GeV range, we will
focus on the limit $v_{S}\gg v$.

###  2.1 Lepton Mass Basis

Upon spontaneous symmetry breaking, the Yukawa interactions contribute as
follows to the mass matrix of the charged leptons:

$\displaystyle\sum_{i,j\in\\{1,2,3\\}}\bar{L}_{i}\mathcal{M}_{ij}e_{j}+\mathrm{h.c.}=$
$\displaystyle\begin{pmatrix}\bar{L}_{1}&\bar{L}_{2}&\bar{L}_{3}\end{pmatrix}\begin{pmatrix}\frac{v_{H}}{\sqrt{2}}y_{11}&0&0\\\
0&\frac{v_{H}}{\sqrt{2}}y_{22}&\frac{v_{\phi}}{\sqrt{2}}y_{23}\\\
\frac{v_{\phi}}{\sqrt{2}}y_{31}&0&\frac{v_{H}}{\sqrt{2}}y_{33}\end{pmatrix}\begin{pmatrix}e_{1}\\\
e_{2}\\\ e_{3}\end{pmatrix}+\mathrm{h.c.}$ (8)

Without loss of generality, we can take the Yukawa couplings to be real, as it
is possible to absorb the complex phases in the field definitions. As a
result, the symmetric matrices $\mathcal{MM}^{T}$ and
$\mathcal{M}^{T}\mathcal{M}$ can be diagonalised by the orthogonal matrices
$O^{L}$ and $O^{R}$, respectively,

$\displaystyle(O^{L})^{T}\mathcal{M}\mathcal{M}^{T}O^{L}=\mathcal{M}_{\rm
diag}^{2}\,,$ (9)
$\displaystyle(O^{R})^{T}\mathcal{M}^{T}\mathcal{M}O^{R}=\mathcal{M}_{\rm
diag}^{2}\,,$ (10)

where $\mathcal{M}_{\rm
diag}^{2}=\text{diag}\begin{pmatrix}m_{e}^{2}&m^{2}_{\mu}&m^{2}_{\tau}\end{pmatrix}$
has non-negative diagonal entries corresponding to the squared charged lepton
masses.

Assuming the flavour off-diagonal Yukawas to be smaller than the flavour
conserving ones, the angles that rotate the gauge eigenbasis into the mass
eigenbasis can be treated perturbatively, and the orthogonal matrices
$O^{L,R}$ can be written as $O^{L,R}=R^{L,R}_{13}R^{L,R}_{23}$, where

$\displaystyle R_{13}^{L,R}=\begin{pmatrix}1&0&\theta_{13}^{L,R}\\\ 0&1&0\\\
-\theta_{13}^{L,R}&0&1\end{pmatrix},\qquad
R^{L,R}_{23}=\begin{pmatrix}1&0&0\\\ 0&1&-\theta_{23}^{L,R}\\\
0&\theta_{23}^{L,R}&1\end{pmatrix}\,.$ (11)

At leading order in the ratio of flavour off-diagonal Yukawa couplings and
flavour diagonal Yukawa couplings, the rotation angles read

$\displaystyle\theta_{13}^{L}$
$\displaystyle\simeq\frac{v_{\phi}}{v_{H}}\frac{y_{31}y_{11}}{y_{33}^{2}}$
$\displaystyle\theta_{23}^{L}$
$\displaystyle\simeq-\frac{v_{\phi}}{v_{H}}\frac{y_{23}}{y_{33}}$
$\displaystyle\theta_{13}^{R}$
$\displaystyle\simeq\frac{v_{\phi}}{v_{H}}\frac{y_{31}}{y_{33}}$
$\displaystyle\theta_{23}^{R}$
$\displaystyle\simeq-\frac{v_{\phi}}{v_{H}}\frac{y_{23}y_{22}}{y_{33}^{2}}$
(12)

and the charged lepton masses are given by

$\text{diag}\begin{pmatrix}m_{e}&m_{\mu}&m_{\tau}\end{pmatrix}\simeq\frac{v_{H}}{\sqrt{2}}\
\text{diag}\begin{pmatrix}y_{11}&y_{22}&y_{33}\end{pmatrix}\,,$ (13)

where $v_{H}$ is the VEV of the $U(1)^{\prime}$ neutral Higgs doublet.

###  2.2 Gauge Sector

The mass matrix of the neutral gauge bosons can be obtained from the kinetic
terms of the scalars in unitary gauge:

$\displaystyle\left(D_{\alpha}H^{\dagger}\right)\left(D^{\alpha}H\right)+\left(D_{\alpha}\phi^{\dagger}\right)\left(D^{\alpha}\phi\right)+\left(D_{\alpha}S^{\dagger}\right)\left(D^{\alpha}S\right)\supset$
$\displaystyle-\frac{1}{8}\begin{pmatrix}B_{\alpha}&W_{\alpha}^{3}&B_{\alpha}^{\prime}\end{pmatrix}\begin{pmatrix}g_{1}^{2}\left(v_{H}^{2}+v_{\phi}^{2}\right)&-g_{1}g_{2}\left(v_{H}^{2}+v_{\phi}^{2}\right)&-2g^{\prime}g_{1}v_{\phi}^{2}\\\
-g_{1}g_{2}\left(v_{H}^{2}+v_{\phi}^{2}\right)&g_{2}^{2}\left(v_{H}^{2}+v_{\phi}^{2}\right)&2g^{\prime}g_{2}v_{\phi}^{2}\\\
-2g^{\prime}g_{1}v_{\phi}^{2}&2g^{\prime}g_{2}v_{\phi}^{2}&4g^{\prime
2}\left(v_{S}^{2}+v_{\phi}^{2}\right)\end{pmatrix}\begin{pmatrix}B^{\alpha}\\\
W^{3\alpha}\\\ B^{\prime\alpha}\end{pmatrix}\,.$ (14)

If $M_{Z^{\prime}}\lesssim 10~{}$GeV, $g^{\prime}$ must be $\lesssim 10^{-4}$
to avoid the constraints on electron fifth-force (see Section 3.1), leading to
suppressed mixing, $\propto g^{\prime}$, between the $B^{\prime}$ boson and
the SM gauge bosons. In the limit $v_{S}\gg v$, we approximate the mass matrix
as follows:

$\displaystyle\frac{1}{8}\begin{pmatrix}B_{\alpha}&W_{\alpha}^{3}&B_{\alpha}^{\prime}\end{pmatrix}\begin{pmatrix}g_{1}^{2}v^{2}&-g_{1}g_{2}v^{2}&0\\\
-g_{1}g_{2}v^{2}&g_{2}^{2}v^{2}&0\\\ 0&0&4g^{\prime
2}v_{S}^{2}\end{pmatrix}\begin{pmatrix}B^{\alpha}\\\ W^{3\alpha}\\\
B^{\prime\alpha}\end{pmatrix}\,.$ (15)

The mass eigenstate $Z^{\prime}$, with $M_{Z^{\prime}}=g^{\prime}v_{s}$, is
aligned with $B^{\prime}$, while the photon $A$ and the $Z$ boson (with
$M_{Z}=v/2\sqrt{g_{1}^{2}+g_{2}^{2}}$) are related to the interaction basis
via a rotation by the Weinberg angle $\theta_{W}$, with $\tan\theta_{W}\equiv
g_{1}/g_{2}$, just as in the SM.

$\begin{pmatrix}A_{\alpha}\\\ Z_{\alpha}\\\
Z^{\prime}_{\alpha}\end{pmatrix}=\begin{pmatrix}\cos\theta_{W}&\sin\theta_{W}&0\\\
-\sin\theta_{W}&\cos\theta_{W}&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}B_{\alpha}\\\ W^{3}_{\alpha}\\\
B^{\prime}_{\alpha}\end{pmatrix}$ (16)

Note that kinetic mixing of the abelian fields, $\epsilon
B^{\mu\nu}B^{\prime}_{\mu\nu}$, is not forbidden by any symmetry and
introduces a new coupling $\epsilon$ that leads to $g^{\prime}$ independent
interactions of the $Z^{\prime}$ boson with the SM fermions. If the
$U(1)^{\prime}$ is the remnant of a larger spontaneously broken non-abelian
gauge group, kinetic mixing is absent at tree-level. However, kinetic mixing
can always be generated via loops involving fermions that are charged under
both abelian groups. These lead to the finite and calculable contribution
$\epsilon_{1-\rm loop}\sim g^{\prime}g_{1}/(16\pi^{2})\log(m_{\mu}/m_{e})$. In
order to obtain a more predictive model, we consider vanishing tree-level
$\epsilon$ but loop-induced kinetic mixing.

The lepton mass eigenstates are related to the gauge eigenstates via
$L_{i}=O^{L}_{i\alpha}L_{\alpha}$, $e_{i}=O^{R}_{i\alpha}\alpha$ with
$\alpha\in\\{e,\mu,\tau\\}$ and $i\in\\{1,2,3\\}$. Whereas the flavour
universal photon and $Z$ couplings are unaffected by these unitary
transformations, the lepton flavour non-universal $Z^{\prime}$ couplings are
sensitive to the misalignment of the lepton gauge and mass eigenbases and
receive flavour changing contributions. The $Z^{\prime}$ interactions with the
charged leptons take the form

$\displaystyle\mathcal{L}^{\rm
LFV}_{Z^{\prime}}=-g^{\prime}Z^{\prime}_{\rho}\bigg{[}$
$\displaystyle\theta^{L}_{13}(\overline{L}_{e}\gamma^{\rho}L_{\tau}+\tau\leftrightarrow
e)+\theta^{L}_{23}(\overline{L}_{\mu}\gamma^{\rho}L_{\tau}+\tau\leftrightarrow\mu)$
$\displaystyle+\theta^{L}_{13}\theta^{L}_{23}(\overline{L}_{e}\gamma^{\rho}L_{\mu}+\mu\leftrightarrow
e)+\theta^{R}_{13}(\overline{e}\gamma^{\rho}\tau+\tau\leftrightarrow e)$
$\displaystyle+\theta^{R}_{23}(\overline{\mu}\gamma^{\rho}\tau+\tau\leftrightarrow\mu)+\theta^{R}_{13}\theta^{R}_{23}(\overline{e}\gamma^{\rho}\mu+\mu\leftrightarrow
e)\bigg{]}$ (17)

where $\theta^{L,R}$ are the rotation angles given in Eq. (12). The $\mu\to e$
couplings are proportional to the product of the off-diagonal $y_{31}\times
y_{23}$ Yukawas that parametrise the $\tau\leftrightarrow e$ and
$\tau\leftrightarrow\mu$ flavour changes. Note that the products
$\theta_{13}^{L}\theta_{23}^{L}$ and $\theta_{13}^{R}\theta_{23}^{R}$ in the
tree-level $\mu\to e$ $Z^{\prime}$ interactions are suppressed by the SM
electron and muon Yukawas respectively. Only the combination
$\theta_{13}^{R}\theta_{23}^{L}$ is sizeable.

###  2.3 Scalar Sector

We supplement the SM Higgs sector by two new scalars, an $SU(2)_{L}$ doublet
$\phi$ and a singlet $S$, both of which are charged under $U(1)^{\prime}$.
Assuming the singlet to be heavier than the doublet, we integrate it out,
reducing the scalar sector to that of a regular two Higgs Doublet Model (2HDM)
with the doublets $H$ and $\phi$. Since the VEV of the singlet breaks
$U(1)^{\prime}$, the scalar potential now features all $SU(2)_{L}\otimes
U(1)_{Y}$ gauge invariant potential terms, including $U(1)^{\prime}$ breaking
interactions, which are generated via singlet VEV insertions:

$\displaystyle\mathcal{V}^{\mathrm{eff}}=$ $\displaystyle
M_{HH}^{2}H^{\dagger}H+M_{\phi\phi}^{2}\phi^{\dagger}\phi-
M_{H\phi}^{2}\left(H^{\dagger}\phi+\phi^{\dagger}H\right)$
$\displaystyle+\frac{1}{2}\Lambda_{H}\left(H^{\dagger}H\right)^{2}+\frac{1}{2}\Lambda_{\phi}\left(\phi^{\dagger}\phi\right)^{2}$
$\displaystyle+\Lambda_{H\phi}\left(H^{\dagger}H\right)\left(\phi^{\dagger}\phi\right)+\tilde{\Lambda}_{H\phi}\left(H^{\dagger}\phi\right)\left(\phi^{\dagger}H\right)$
$\displaystyle+\frac{1}{2}\Lambda_{a}\left(\left(H^{\dagger}\phi\right)^{2}+\left(\phi^{\dagger}H\right)^{2}\right)$
$\displaystyle+\Lambda_{b}\left(H^{\dagger}H\right)\left(\left(H^{\dagger}\phi\right)+\left(\phi^{\dagger}H\right)\right)+\Lambda_{c}\left(\phi^{\dagger}\phi\right)\left(\left(H^{\dagger}\phi\right)+\left(\phi^{\dagger}H\right)\right)\,.$
(18)

It is particularly convenient to rotate into the so-called Higgs basis, where
only one of the doublets acquires a vacuum expectation value:

$\begin{pmatrix}H_{1}\\\
H_{2}\end{pmatrix}=\begin{pmatrix}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{pmatrix}\begin{pmatrix}H\\\
\phi\end{pmatrix}\qquad\text{with}\ \tan\beta\equiv v_{\phi}/v_{H}$ (19)

$\expectationvalue{H_{1}}=v/\sqrt{2}$, with $v=\sqrt{v_{\phi}^{2}+v_{H}^{2}}$,
and $\expectationvalue{H_{2}}=0$. We relabel the potential parameters as
follows

$\displaystyle V(H_{1},H_{2})$
$\displaystyle=m^{2}_{11}H^{\dagger}_{1}H_{1}+m^{2}_{22}H^{\dagger}_{2}H_{2}-m^{2}_{12}(H^{\dagger}_{1}H_{2}+\text{h.c})+\frac{\lambda_{1}}{2}(H^{\dagger}_{1}H_{1})^{2}+\frac{\lambda_{2}}{2}(H^{\dagger}_{2}H_{2})^{2}$
$\displaystyle+\lambda_{3}(H^{\dagger}_{1}H_{1})(H^{\dagger}_{2}H_{2})+\lambda_{4}(H^{\dagger}_{1}H_{2})(H^{\dagger}_{2}H_{1})$
$\displaystyle+\left(\frac{\lambda_{5}}{2}(H^{\dagger}_{1}H_{2})^{2}+\lambda_{6}(H^{\dagger}_{1}H_{1})(H_{1}H^{\dagger}_{2})+\lambda_{7}(H^{\dagger}_{2}H_{2})(H_{1}H^{\dagger}_{2})+\text{h.c}\right),$
(20)

and expand the doublets in terms of the canonically normalized scalar fields

$H_{1}=\begin{pmatrix}G^{+}\\\
\frac{1}{\sqrt{2}}(v+H^{0}_{1}+iG^{0})\end{pmatrix}\,,\qquad
H_{2}=\begin{pmatrix}H^{+}\\\
\frac{1}{\sqrt{2}}(H^{0}_{2}+iA)\end{pmatrix}\,.$ (21)

Once the Goldstone bosons $G$ are eaten by the electroweak gauge bosons, the
spectrum contains one charged scalar $H^{+}$, two CP even scalars
$H^{0}_{1},H^{0}_{2}$ and one CP odd neutral scalar $A$. If the potential
parameters are real, the physical states are CP eigenstates and only the two
CP even neutral scalars, $H_{1}^{0}$ and $H_{2}^{0}$, can mix. The mass matrix
of $H_{1}^{0}$ and $H_{2}^{0}$ is diagonalised by the angle $\beta-\alpha$. We
identify two scalar mass-eigenstates

$\displaystyle h=\sin(\beta-\alpha)H^{0}_{1}+\cos(\beta-\alpha)H^{0}_{2}$ (22)
$\displaystyle\rho=\cos(\beta-\alpha)H^{0}_{1}-\sin(\beta-\alpha)H^{0}_{2}$
(23)

with the masses $m_{h}$ and $m_{\rho}$, respectively. The angle $\beta-\alpha$
can be written in terms of $\lambda_{6}$, which is defined in Eq. (20), as
well as the scalar masses $m_{h}$ and $m_{\rho}$ Davidson:2005cw ;
Davidson:2016utf :

$\cos(\beta-\alpha)\sin(\beta-\alpha)\equiv
c_{\beta\alpha}s_{\beta\alpha}=-\frac{\lambda_{6}v^{2}}{(m^{2}_{\rho}-m_{h}^{2})}$
(24)

The decoupling limit of the 2HDM is obtained when the mass term of the VEV-
less doublet $H_{2}$ satisfies the condition $m_{22}\equiv M\gg v$. The
doublet mass terms in the Higgs basis, $m_{11}$, $m_{12}$, $m_{22}$ (see Eq.
(20)), are related to the mass matrix in the $H,\phi$ basis (see Eq. (18)) in
the following way

$m^{2}_{11}=c_{\beta}^{2}\left(M^{2}_{HH}+M^{2}_{\phi\phi}t_{\beta}^{2}-M^{2}_{H\phi}\frac{2t_{\beta}}{\sqrt{1+t_{\beta}^{2}}}\right)$
(25)
$m^{2}_{22}=c_{\beta}^{2}\left(M^{2}_{\phi\phi}+M^{2}_{HH}t_{\beta}^{2}+M^{2}_{H\phi}\frac{2t_{\beta}}{\sqrt{1+t_{\beta}^{2}}}\right)$
(26)
$m^{2}_{12}=\frac{1}{2}\left(M^{2}_{HH}-M^{2}_{\phi\phi}\right)s_{2\beta}+M^{2}_{H\phi}c_{2\beta}.$
(27)

The potential minimum conditions imply

$m^{2}_{11}=-\frac{1}{2}\lambda_{1}v^{2}\qquad
m^{2}_{12}=\frac{1}{2}\lambda_{6}v^{2}\,.$ (28)

Consequently, if we assume $\lambda_{1}$ and $\lambda_{6}$ to be perturbative
couplings, $m_{11}$ and $m_{12}$ are sub-electroweak masses. The decoupling
condition $m_{22}\gg v$ is satisfied if $M_{HH},M_{\phi\phi},$ $M_{H\phi}\gg
v$, however, in this case the minimum potential conditions require a fine-
tuned cancellation in Eqs. (25) and (27). If, on the other hand, $\beta\ll 1$,
the mass configuration $M_{H}\sim v$ and $M_{\phi},M_{H\phi}\gg v$ leads to
natural minimum potential conditions in the decoupling limit. Once $H_{2}$ is
decoupled, we identify the light CP even scalar $h$ as the 125-GeV Higgs
boson, while the orthogonal state $\rho$ has a mass $m_{\rho}\sim m_{22}\sim
M$. In the decoupling limit, the relation Davidson:2005cw

$m_{\rho}^{2}-m^{2}_{A}=(\lambda_{5}+\lambda_{1})v^{2}-m_{h}^{2},$ (29)

implies that the pseudoscalar $A$ and the scalar $\rho$ are approximately
degenerate, with a mass-splitting of order $v^{2}$.

In the lepton mass basis, the doublet $H_{1}$, which is aligned with the
electroweak vacuum, has diagonal Yukawa interactions with the leptons, while
the $H_{2}$ couplings are in general flavour-changing. Written in terms of the
Yukawa couplings defined in Eq. (4), the lepton Yukawa sector in the Higgs
basis reads

$\mathcal{L}_{\rm
Yuk}=\sum_{i,j\in\\{e,\mu,\tau\\}}\left[\frac{\sqrt{2}\delta_{ij}m_{i}}{v}\bar{L}_{i}e_{i}H_{1}+(O^{T}_{L}Y_{2}O_{R})_{ij}\bar{L}_{i}e_{j}H_{2}+\mathrm{h.c}\right]\,,$
(30)

where

$Y_{2}=\begin{pmatrix}-y_{11}\sin\beta&0&0\\\
0&-y_{22}\sin\beta&y_{23}\cos\beta\\\
y_{31}\cos\beta&0&-y_{33}\sin\beta\end{pmatrix}.$ (31)

Defining $\mathcal{Y}\equiv O^{T}_{L}Y_{2}O_{R}$, the couplings of the neutral
scalar mass eigenstates to the leptons take the form

$\displaystyle\mathcal{L}^{\rm neut}_{\rm Yuk}=\sum_{i,j\in\\{e,\mu,\tau\\}}$
$\displaystyle\frac{h}{\sqrt{2}}\bar{e}_{i}\left[s_{\beta\alpha}\frac{\sqrt{2}\delta_{ij}m_{i}}{v}(P_{R}+P_{L})+c_{\beta\alpha}(\mathcal{Y}_{ij}P_{R}+\mathcal{Y}^{\dagger}_{ij}P_{L})\right]e_{j}$
$\displaystyle+$
$\displaystyle\frac{\rho}{\sqrt{2}}\bar{e}_{i}\left[c_{\beta\alpha}\frac{\sqrt{2}\delta_{ij}m_{i}}{v}(P_{R}+P_{L})-s_{\beta\alpha}(\mathcal{Y}_{ij}P_{R}+\mathcal{Y}^{\dagger}_{ij}P_{L})\right]e_{j}$
$\displaystyle+$
$\displaystyle\frac{iA}{\sqrt{2}}\bar{e}_{i}\left[\mathcal{Y}_{ij}P_{R}-\mathcal{Y}^{\dagger}_{ij}P_{L}\right]e_{j}\,.$
(32)

Since the rate of $h\to\tau^{+}\tau^{-}$ measured at the LHC ATLAS:2018ynr ;
CMS:2017zyp  is compatible with the Standard Model prediction, we require
$s_{\beta\alpha}\sim 1$. As a result, the flavour-changing couplings of the
Higgs boson $h$ are suppressed by

$c_{\beta\alpha}\simeq-\frac{\lambda_{6}v^{2}}{(m^{2}_{\rho}-m_{h}^{2})}\ll 1\
\to\ c_{\beta\alpha}\simeq-\frac{\lambda_{6}v^{2}}{M^{2}}.$ (33)

##  3 Phenomenology

In this section we discuss the phenomenological signatures of our model. In
Section 3.1 we briefly review the experimental constraints on $L_{e}-L_{\mu}$
gauge interactions, before focussing on the LFV phenomenology in Section 2. In
Section 3.3 we address neutrino mixing and the generation of neutrino masses.

###  3.1 $L_{e}-L_{\mu}$ Gauge Interactions

Several experiments search for $Z^{\prime}$ bosons that interact with SM
particles. In the absence of tree-level kinetic mixing, these searches
directly probe the size of the gauge coupling $g^{\prime}$. Since in our model
the $Z^{\prime}$ interactions arise from gauging the lepton flavour difference
$L_{e}-L_{\mu}$, $g^{\prime}$ is primarily constrained by bounds on
$Z^{\prime}ee$ and $Z^{\prime}\nu\nu$ Wise:2018rnb ; Bauer:2018onh .

Figure 1: Figure taken from Ref. Bauer:2018onh  showing the current bounds on
the space of $M_{Z^{\prime}}$ vs. $g^{\prime}$ in the case of gauged
$L_{e}-L_{\mu}$. In this article we consider two benchmark scenarios, one of
which involves a $Z^{\prime}$ boson with a mass $M_{Z^{\prime}}\sim 10~{}$GeV
and a coupling $g^{\prime}\sim 10^{-4}$ and is in reach of Belle-II, and
another benchmark with $M_{Z^{\prime}}\sim 15~{}$GeV, which avoids the bounds
from BaBar BaBar1 ; BaBar2 , allowing for a gauge coupling of $g^{\prime}\sim
2.5\times 10^{-3}$.

In the $Z^{\prime}$ mass range of $10\ \text{MeV}\lesssim
M_{Z^{\prime}}\lesssim 1$ GeV, the strongest constraints on the
$U(1)^{\prime}$ gauge coupling come from electron beam dump experiments
ElectronSearches1 ; electronSearches2 , neutrino oscillation experiments
Wise:2018rnb  and neutrino scattering experiments neutrinoscattering1 ;
neutrinoscattering2 , while cosmological and astrophysical limits are more
relevant for lighter $Z^{\prime}$ bosons Kamada:2015era ; Rrapaj:2015wgs ;
Knapen:2017xzo . Muonium spectroscopy can also probe for $Z^{\prime}$ bosons
in the sub-MeV mass range Aiba:2021bxe . For larger masses,
$M_{Z^{\prime}}\gtrsim 1$ GeV, colliders are the most sensitive probes of a
fifth force. Below $10$ GeV, BaBar BaBar1 ; BaBar2  set an upper limit of
$g^{\prime}\lesssim 10^{-4}$ by searching for $Z^{\prime}$ production in
combination with a single photon, $e^{+}e^{-}\to\gamma
Z^{\prime}\to\gamma\ell^{+}\ell^{-}$. Belle II is expected to push this limit
down to $g^{\prime}\sim 10^{-5}$ Belle-II:2010dht ; Inguglia:2016acz .
$Z^{\prime}$ bosons with masses beyond $M_{Z^{\prime}}\sim 10$ GeV can be
produced in the decay of an excited meson produced at the LHC, however the
process occurs via loop induced kinetic mixing and the bound on the
$U(1)^{\prime}$ gauge coupling is only $g^{\prime}\lesssim 10^{-1}-10^{-2}$
LHCb:2017trq ; Ilten:2015hya .

###  3.2 Lepton Flavour Violation

####  3.2.1 $\tau\to\mu$ and $\tau\to e$ Transitions

$\tau$$\ell$$\ell$$\ell$$Z^{\prime}$ Figure 2: $Z^{\prime}$ mediated $\tau\to
3\ell$ decay, with $\ell=e,\mu$. If
$2m_{\ell}<M_{Z^{\prime}}<m_{\tau}-m_{\ell}$, the decay can happen via the on-
shell production of the $Z^{\prime}$ and its subsequent decay into a lepton
pair.

Via the lepton flavour changing $Z^{\prime}$ couplings in Eq. (17), $\tau$
leptons can decay into final states such as $3e,3\mu,\mu 2e,e2\mu$. In Figure
2, we show the diagram for the tree-level $Z^{\prime}$ mediated decay of
$\tau$ into $3\ell$, where $\ell=e,\mu$. For a $Z^{\prime}$ lighter than the
mass difference $m_{\tau}-m_{\ell}$, $\tau$ leptons can decay into an on-shell
$Z^{\prime}$ and a charged lepton $\ell$. If in addition
$M_{Z^{\prime}}>2m_{\ell}$, the $Z^{\prime}$ boson can decay into a lepton
pair $\ell^{+}\ell^{-}$. If this happens before the $Z^{\prime}$ escapes the
detector, a three-lepton final state may be measured as a consequence of the
decay chain $\tau\to Z^{\prime}\ell\to\ell^{+}\ell^{-}\ell$. In the narrow-
width approximation, $\Gamma/M_{Z^{\prime}}\lesssim g^{\prime 2}\ll 1$, the
Breit-Wigner distribution can be approximated by a $\delta-$function,
resulting in a factorized rate

$\Gamma(\tau\to 3\ell)\simeq\Gamma(\tau\to
Z^{\prime}\ell)\times\mathrm{Br}(Z^{\prime}\to\ell^{+}\ell^{-})\,.$ (34)

If the muon and electron masses can be neglected with respect to the
$Z^{\prime}$ mass, we can take $\mathrm{Br}(Z^{\prime}\to\ell^{+}\ell^{-})\sim
1/4$ and 222Note that there is no divergence in the limit $M_{Z^{\prime}}\to
0$, because with $v_{S},v_{\phi}\to 0$, the LFV rotation angles $\theta_{i3}$
also vanish. Ibarra:2021xyk

$\Gamma(\tau\to
Z^{\prime}\ell)=\sum_{X=L,R}\frac{(g^{\prime}\theta^{X}_{i3})^{2}}{64\pi}\frac{m^{3}_{\tau}}{m_{Z^{\prime}}^{2}}\left(1-\frac{m_{Z^{\prime}}^{2}}{m^{2}_{\tau}}\right)^{2}\left(1+2\frac{m^{2}_{Z^{\prime}}}{m^{2}_{\tau}}\right)\,,$
(35)

where $i=1,2$ for $\ell=e,\mu$ respectively. (The mixing angles $\theta$ are
given in terms of the Lagrangian parameters in Eq. (12)). In the case of on-
shell $Z^{\prime}$ production, the rate of $\Gamma(\tau\to 3\ell)$ scales with
$g^{\prime 2}$ rather than with the $g^{\prime 4}$ that one would naively
expect from the tree-level $Z^{\prime}$ exchange shown in Figure 2.
Consequently, the bounds on the flavour off-diagonal Yukawas from
$\mathrm{Br}(\tau\to 3e)<2.7\times 10^{-8}$, $\mathrm{Br}(\tau\to
3\mu)<2.1\times 10^{-8}$ are stringent enough to suppress all other LFV
signals. Indeed, taking $M_{Z^{\prime}}=1~{}$GeV and $g^{\prime}=10^{-4}$, the
upper limit $\mathrm{Br}(\tau\to 3\mu)<2.1\times 10^{-8}$ implies
$y_{23}<3\times 10^{-7}(\sin\beta)^{-1}$.

If $M_{Z^{\prime}}>m_{\tau}$, the $\tau\to 3\ell$ decay is mediated at tree-
level by an off-shell $Z^{\prime}$ boson . Even for a $Z^{\prime}$ mass of
$M_{Z^{\prime}}\sim 5$ GeV, treating the $Z^{\prime}$ exchange as a contact
interaction between four leptons is sufficient for an estimation of the decay
rate to a $\sim 10\%$ accuracy. Integrating out the $Z^{\prime}$ boson at the
scale of $M_{Z^{\prime}}$, the four-lepton operators

$\mathcal{O}_{V,XY}=(\bar{\ell}\gamma^{\alpha}P_{X}\tau)(\bar{\ell}\gamma_{\alpha}P_{Y}\ell)$
(36)

are generated with the following coefficients

$\displaystyle C^{e\tau ee}_{V,XY}=-\frac{g^{\prime
2}\theta^{X}_{13}}{M^{2}_{Z^{\prime}}}\qquad
C^{\mu\tau\mu\mu}_{V,XY}=\frac{g^{\prime
2}\theta^{X}_{23}}{M^{2}_{Z^{\prime}}}.$ (37)

The decay rate of $\tau\to 3\ell$ is then given by

$\displaystyle{\rm Br}\left[\tau\to
3\ell\right]=\frac{m_{\tau}^{5}}{1536\,\pi^{3}\,\Gamma_{\tau}}\left(2\left|C^{\ell\tau\ell\ell}_{V,LL}\right|^{2}+\left|C^{\ell\tau\ell\ell}_{V,LR}\right|^{2}+L\leftrightarrow
R\right)\,,$ (38)

where $\Gamma_{\tau}$ is the total decay width of the $\tau$. Here we are
neglecting the QED running from $M_{Z^{\prime}}$ to the $\tau$ mass.

Contributions from the scalar sector arise from diagrams similar to the one
depicted in Figure 2, but with a neutral scalar exchange. However, the
flavour-diagonal current couples to scalars via the light lepton Yukawas, and
LFV Higgs decays are more sensitive to $\tau\leftrightarrow\ell$ flavour
changing Yukawa couplings than the Yukawa suppressed $\tau\to 3\ell$ decay
AtlasHiggs . The width for the LFV Higgs decay is given by 2HDMZupan

$\Gamma(h\to\tau\ell)=\frac{\absolutevalue{\mathcal{Y}_{\ell\tau}}^{2}+\absolutevalue{\mathcal{Y}_{\tau\ell}}^{2}}{16\pi}c_{\beta\alpha}^{2}m_{h}$
(39)

with $\mathcal{Y}$ as defined in Eq. (32). The allowed region in the
$y_{31}-y_{23}$ plane is plotted in Figure 3, which shows that
$\mathrm{Br}(\tau\to 3e)<2.7\times 10^{-8}$, $\mathrm{Br}(\tau\to
3\mu)<2.1\times 10^{-8}$ and $\mathrm{Br}(h\to\tau\mu)<1.5\times
10^{-3},\mathrm{Br}(h\to\tau e)<2.2\times 10^{-3}$ AtlasHiggs  are all
compatible with perturbative Yukawa couplings. Note, however, that these
bounds are not able to constrain $y_{31}$ and $y_{23}$ to values smaller than
$y_{33}/t_{\beta}\sim 1/10$, as required in Eq. (12).

Figure 3: In the scenario with $g^{\prime}=10^{-4}$,
$M_{Z^{\prime}}=10~{}$GeV, $t_{\beta}=\frac{1}{10}$ and $M=10~{}$TeV, the
bounds from lepton flavour violation in the $\tau$ sector cannot constrain the
off-diagonal Yukawa couplings $y_{31}$ and $y_{23}$.

####  3.2.2 $\mu\to e$ Transitions

We now turn to the $\mu\to e$ phenomenology. In presence of the flavour mixing
angles of Eq. (12), the $Z^{\prime}$ boson can couple at tree-level to $\mu\to
e$ currents, however these couplings only arise at second order in the mixing
angles and are suppressed by the electron and muon Yukawas: the left-handed
and right-handed vector currents come with the products of mixing angles
$\theta^{L}_{13}\theta^{L}_{23}$ and $\theta^{R}_{13}\theta^{R}_{23}$
respectively, where

$\displaystyle\theta_{23}^{L}$
$\displaystyle=\tan\beta\frac{y_{23}}{y_{33}}\qquad{\theta^{R}_{23}=\theta_{23}^{L}\times\frac{y_{22}}{y_{33}}}$
$\displaystyle\theta_{13}^{R}$
$\displaystyle=\tan\beta\frac{y_{31}}{y_{33}}\qquad{\theta^{L}_{13}=\theta_{13}^{R}\times\frac{y_{11}}{y_{33}}}\,.$
(40)

As a result, the $\mu\to 3e$ decay rate, upon which the bound
$\mathrm{Br}(\mu\to 3e)<10^{-12}$ applies, receives contributions from the
diagram of Figure 4(a). The effective four-fermion interactions that result
from integrating out the $Z^{\prime}$ boson, are given by

$C^{e\mu ee}_{V,XY}=-\frac{g^{\prime
2}\theta^{X}_{13}\theta^{X}_{23}}{M^{2}_{Z^{\prime}}}\,.$ (41)

Contributions of the dipole, which comes with the coupling product
$\theta^{R}_{13}\theta^{L}_{23}$ and avoids the suppression by the SM electron
and muon Yukawas, can also be relevant.

$\mu$$e$$e$$e$$Z^{\prime}$ (a)

$\mu_{L}$$e_{R}$$\theta^{L}_{23}$$\theta^{R}_{13}$$\gamma$$Z^{\prime}$$\tau$
(b)

Figure 4: (a) Tree-level $Z^{\prime}$ exchange contributing to $\mu\to 3e$.
(b) $Z^{\prime}$ penguin diagram contributing to $\mu\to e\gamma$.

For instance, the penguin diagram of Figure 4(b) contributes to the $\mu\to e$
photon dipole

$\mathcal{O}^{e\mu}_{D,X}=m_{\mu}(\bar{e}\sigma_{\alpha\beta}P_{X}\mu)F^{\alpha\beta},$
(42)

with a coefficient

$(C^{e\mu}_{D,L})_{peng}=-\frac{3e}{48\pi^{2}M^{2}_{Z^{\prime}}}\left(\frac{m_{\tau}}{m_{\mu}}\right)g^{\prime
2}\theta^{R}_{13}\theta^{L}_{23}.$ (43)

In addition to the penguin, loop diagrams with scalars can give sizeable
contributions to the dipole. In the diagrams of Figure 5(a), the mass
insertion flips the chirality in the virtual $\tau$ line. These diagrams
contribute to the dipole coefficient as follows Chang:1993kw

$\displaystyle(C^{e\mu}_{D,L})_{\mathrm{1-loop}}=-\frac{e}{64\pi^{2}}\left(\frac{m_{\tau}}{m_{\mu}}\right)\bigg{[}\frac{\mathcal{Y}^{*}_{\mu\tau}\mathcal{Y}^{*}_{e\tau}}{m^{2}_{h}}c^{2}_{\beta\alpha}\left(2\log(\frac{m_{\tau}}{m_{h}})+\frac{3}{2}\right)$
$\displaystyle+\frac{\mathcal{Y}^{*}_{\mu\tau}\mathcal{Y}^{*}_{e\tau}}{m^{2}_{\rho}}s^{2}_{\beta\alpha}\left(2\log(\frac{m_{\tau}}{m_{\rho}})+\frac{3}{2}\right)$
$\displaystyle-\frac{\mathcal{Y}^{*}_{\mu\tau}\mathcal{Y}^{*}_{e\tau}}{m^{2}_{A}}\left(2\log(\frac{m_{\tau}}{m_{A}})+\frac{3}{2}\right)\bigg{]}\,,$
(44)

with $\mathcal{Y}$ as defined in Eq. (32). $C_{D,R}$ is obtained from
$C_{D,L}$ by replacing $\mathcal{Y}\leftrightarrow\mathcal{Y}^{\dagger}$. When
the loop is closed by the light Higgs $h$, each flavour changing vertex is
proportional to $c_{\beta\alpha}$, resulting in a contribution suppressed by
four powers of the heavy Higgs mass, $c^{2}_{\beta\alpha}\sim 1/M^{4}$.
Moreover, in the decoupling limit, the contributions of the heavy scalar
$\rho$ and the pseudoscalar $A$ cancel in the above equation, since their
masses are degenerate at leading order in $1/M^{2}$. The total amplitude ends
up being a sub-dominant $1/M^{4}$ effect.

$\mu$$e$$\gamma$$\Phi$$\tau$ (a)

$\mu$$e$$\gamma$$\mu$$t$$\gamma$$\Phi$ (b)

$\mu$$e$$\gamma$$\mu$$W$$\gamma$$\Phi$ (c)

Figure 5: Loop diagrams with neutral scalars contributing to the rate of
$\mu\to e\gamma$. $\Phi=h,\rho,A$

The misalignment of the lepton mass eigenbasis and the lepton gauge eigenbasis
also affects the scalar couplings. Whereas the $Z^{\prime}$ boson acquires
$\mu\leftrightarrow e$ couplings proportional to
$\propto\theta^{L}_{13}\theta^{L}_{23}$ or
$\propto\theta^{R}_{13}\theta^{R}_{23}$, the scalars acquire $\mu\to e$
couplings proportional to $\theta^{R}_{13}\theta^{L}_{23}$.
It is well known that the Barr-Zee diagrams of Figures 5(b)-5(c) give the
leading contribution to $\mu\to e\gamma$ in 2HDM models with $\mu\to e$ Yukawa
couplings PhysRevLett.38.622 . The branching ratio of $\mu\to e\gamma$ is
given by

$\mathrm{Br}(\mu\to
e\gamma)=\frac{m_{\mu}^{5}}{4\pi\Gamma_{\mu}}(\absolutevalue{C^{e\mu}_{D,L}}^{2}+\absolutevalue{C^{e\mu}_{D,R}}^{2})\,,$
(45)

where $\Gamma_{\mu}\simeq G_{F}^{2}m_{\mu}^{5}/(192\pi^{3})$ is the total
decay rate of the muon, and, to a good approximation,

$C^{e\mu}_{D,X}=(C^{e\mu}_{D,X})_{\mathrm{peng}}+(C^{e\mu}_{D,X})_{\mathrm{t-loop}}+(C^{e\mu}_{D,X})_{\mathrm{W-loop}}$
(46)

where $(C^{e\mu}_{D,X})_{\mathrm{t-loop}}$ and
$(C^{e\mu}_{D,X})_{\mathrm{W-loop}}$ correspond to the two-loop dipole
contributions of Figures 5(b) and 5(c), respectively, and are defined in
Appendix A.

For the $\mu\to 3e$ rate we consider the tree-level diagram of Figure 4(a) and
the dipole contributions arising from attaching an electron current to the
photon. We find

$\displaystyle{\rm Br}\left[\mu\to 3e\right]=$
$\displaystyle-\frac{e^{2}\,m_{\mu}^{5}}{192\,\pi^{3}\,\Gamma_{\mu}}\left(|C_{D,L}^{e\mu}|^{2}+|C_{D,R}^{e\mu}|^{2}\right)\left(4\log\left(\frac{m_{e}^{2}}{m_{\mu}^{2}}\right)+11\right)$
$\displaystyle+\frac{m_{\mu}^{5}}{1536\,\pi^{3}\,\Gamma_{\mu}}\left(2\left|C^{e\mu
ee}_{V,LL}\right|^{2}+\left|C^{e\mu ee}_{V,LR}\right|^{2}+L\leftrightarrow
R\right)$ (47)
$\displaystyle+\frac{e\,m_{\mu}^{5}}{192\,\pi^{3}\,\Gamma_{\mu}}\Big{(}{\rm
Re}\left[C_{D,R}^{e\mu*}\left(2\,C^{e\mu ee}_{V,LL}+C^{e\mu
ee}_{V,LR}\right)\right]+L\leftrightarrow R\Big{)}\,,$

with $C^{e\mu ee}_{V,X}$ as defined in Eq. (41).

LFV interactions with quarks via $Z^{\prime}$ exchange arise from loop induced
kinetic mixing or via penguin diagrams with an external quark current. The
kinetic mixing is both loop-suppressed and suppressed by the Yukawa couplings
featuring in the $Z^{\prime}\mu e$ vertex, leading to a sub-dominant
contribution. The relevant scalar contributions are the tree-level exchanges
with the $u$, $d$ and $s$ quarks in the nuclei, the Barr-Zee dipoles that
already featured in Eq. (46) and are discussed in Appendix A, as well as the
diagrams in Figure 6. The latter contribute to the gluon operator
$\mathcal{O}_{GG,X}=(\bar{e}P_{X}\mu)G^{\alpha\beta}G_{\alpha\beta}$, which
enters the $\mu\to e$ conversion rate. In order to compute these
contributions, we integrate out the heavy scalar doublet $H_{2}$ at the mass
scale $M$ and match onto the SMEFT operators $\mathcal{O}_{Ledq}^{e\mu
ii}\equiv(\bar{L}_{e}\mu)(\bar{d}_{i}q_{i})$, $\mathcal{O}_{Lequ}^{(1)e\mu
ii}\equiv(\bar{L}_{e}\mu)(\bar{q}_{i}u_{i})$, and
$\mathcal{O}_{eH}^{e\mu}\equiv(H_{1}^{\dagger}H_{1})(\bar{L}_{e}\mu H_{1})$,
which leads to the off-diagonal Yukawa couplings of the SM Higgs boson $h$,
given in Eq. (32). We include the QCD running of the scalar operators from the
scale $M$ down to the electroweak scale. At the electroweak scale, the SMEFT
operators are matched at tree-level onto the Low Energy Effective Field Theory
(LEFT) scalar contact interactions and the gluon operator
$\mathcal{O}_{GG,X}=(\bar{e}P_{X}\mu)G^{\alpha\beta}G_{\alpha\beta}$
Davidson:2020ord , which receives contributions from the diagrams with
$t$-loops. Finally, we compute the $\mu\to e$ conversion rate using Eq. (2.22)
of Davidson:2020hkf , which is written in terms of LEFT operators evaluated at
the electroweak scale and includes the running of the Wilson coefficients down
to the experiment.

$\mu$$e$$\Phi$$q$ Figure 6: Contributions via heavy quark loops to the gluon
operator $\mathcal{O}_{GG,X}=(\bar{e}P_{X}\mu)G^{\alpha\beta}G_{\alpha\beta}$
that enters the $\mu\to e$ conversion rate. $q=c,b,t$ and $\Phi=h,\rho,A$.

As shown in Figure 7, the MEG bound, $\mathrm{Br}(\mu\to e\gamma)<4.2\times
10^{-13}$, is better at constraining the product of Yukawa couplings
$y_{23}\times y_{31}$ than the current upper limits $\mathrm{Br}(\mu\to
3e)<10^{-12}$ and $\mathrm{Br}(\mu A\to eA)<7\times 10^{-13}$. On the other
hand, future $\mu\to e$ conversion and $\mu\to 3e$ experiments expect an
impressive improvement in the branching ratio sensitivities. The COMET and
Mu2e collaborations aim at $\mathrm{Cr}(\mu\mathrm{Al}\to e\mathrm{Al})\sim
10^{-16}$ COMET:2009qeh ; Mu2e:2014fns  and Mu3e at $\mathrm{Br}(\mu\to
3e)\sim 10^{-16}$ Blondel:2013ia , both of which can surpass the expected
sensitivity of MEG-II with $\mathrm{Br}(\mu\to e\gamma)\sim 6\times 10^{-14}$
MEGII:2018kmf  (Figure 8).

Figure 7: Constraints on the Yukawa couplings $y_{31}$ and $y_{23}$ from
experimental upper bounds on $\mu\to e$ transitions. The orange region
corresponds to the allowed parameter space. In the plot we consider
$g^{\prime}=10^{-4}$, $M_{Z^{\prime}}=10$ GeV, $M=10$ TeV, $\lambda_{6}=1$,
$\tan\beta=0.1$.  Figure 8: Constraints on the Yukawa couplings $y_{31}$ and
$y_{23}$, derived from future experimental sensitivities to $\mu\to e$
transitions. The light orange region corresponds to the parameter space that
will be probed by Mu3e, whereas the purple region is the allowed parameter
space if no $\mu\to 3e$ signal is observed by Mu3e. For this plot we consider:
$g^{\prime}=10^{-4}$, $M_{Z^{\prime}}=10$ GeV, $M=10$ TeV, $\lambda_{6}=1$,
$\tan\beta=0.1$.

A $Z^{\prime}$ boson with $M_{Z^{\prime}}=10$ GeV and $g^{\prime}=10^{-4}$, as
considered in Figures 7 and 8, will be probed at Belle-II in the upcoming
years. In this scenario $\mu\to e$ rates are the result of contributions both
from the $Z^{\prime}$ and from scalar LFV interactions. The scalar
contributions dominate over the $Z^{\prime}$ contributions in this region of
parameter space.

Figure 9: Constraints on the Yukawa couplings $y_{31}$ and $y_{23}$ from
experimental upper bounds on $\mu\to e$ transitions for $g^{\prime}=2.5\times
10^{-3}$, $M_{Z^{\prime}}=15$ GeV, $M=10$ TeV, $\lambda_{6}=1$,
$\tan\beta=0.1$. Since in this case the contributions of the scalar singlet
cannot be neglected with respect to those of the scalar doublet, we only show
the vector contributions.

By considering a marginally larger mass, $M_{Z^{\prime}}=15$ GeV, the BaBar
constraint $g\lesssim 10^{-4}$ can be avoided and the most stringent upper
limit on $g^{\prime}$ becomes $g^{\prime}\lesssim 10^{-1}-10^{-2}$ from LHCb
Ilten:2015hya ; LHCb:2017trq . Taking $g^{\prime}\sim 2.5\times 10^{-3}$ and
$M_{Z^{\prime}}\sim 15$ GeV, the singlet has a mass $M_{S}\simeq v_{s}\simeq
M_{Z^{\prime}}/g^{\prime}\sim 10\ \text{TeV}$, which is of the same order as
the heavy scalar masses of the 2HDM sector, and thus cannot be integrated out.
In Figure 9 we show the $\mu\to e$ sensitivities in the $y_{31}-y_{23}$ plane
considering only vector contributions and $g^{\prime}=2.5\times 10^{-3}$,
$M_{Z^{\prime}}=15\ {\rm GeV}$, while the inclusion of the scalar diagrams
would require a more careful analysis of the singlet-doublet mixing.

###  3.3 Neutrino masses

After spontaneous symmetry breaking, the neutrino mass Lagrangian takes the
form

$\displaystyle\mathcal{M}_{\nu}$
$\displaystyle=(M_{D})_{ij}\bar{L}_{i}N_{j}+\frac{1}{2}\overline{N^{c}_{i}}(M_{N})_{ij}N_{j}+\mathrm{h.c.}\,,$

where

$M_{D}=\begin{pmatrix}m^{\nu}_{11}&&m^{\nu}_{13}\\\ &m^{\nu}_{22}&\\\
&m^{\nu}_{32}&m^{\nu}_{33}\end{pmatrix}\qquad
M_{N}=\begin{pmatrix}&M^{N}_{12}&M^{N}_{13}\\\ \cdot&&\\\
\cdot&&M^{N}_{33}\end{pmatrix}\,.$ (49)

When the Majorana masses $M_{N}$ are much larger than the Dirac masses
$M_{D}$, the mass matrix of the light active neutrinos is obtained via the
celebrated seesaw formula MINKOWSKI1977421

$M_{\nu}=-M_{D}M^{-1}_{R}M^{T}_{D}\,.$ (50)

The PMNS matrix $U$ diagonalises the neutrino mass matrix
$U^{T}M_{\nu}U=\text{diag}\begin{pmatrix}m_{1}&m_{2}&m_{3}\end{pmatrix}$ and
is canonically parametrized as

$U=\begin{pmatrix}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta}\\\
-s_{12}c_{23}-c_{12}s_{23}s_{13}e^{i\delta}&c_{12}c_{23}-s_{12}s_{23}s_{13}e^{i\delta}&s_{23}c_{13}\\\
s_{12}s_{23}-c_{12}c_{23}s_{13}e^{i\delta}&-c_{12}s_{23}-s_{12}c_{23}s_{13}e^{i\delta}&c_{23}c_{13}\end{pmatrix}\times\text{diag}\begin{pmatrix}1&e^{i\alpha_{12}}&e^{i\alpha_{31}}\end{pmatrix}\,.$
(51)

Assuming the so-called Normal Ordering (NO) $m_{1}<m_{2}<m_{3}$, a recent
global fit deSalas:2020pgw  to neutrino oscillation data gives the following
values for the mass squared differences and mixing angles:

$\displaystyle m^{2}_{2}-m^{2}_{1}=\left[6.94-8.14\right]\times 10^{-5}\
e\text{V}^{2}\qquad\absolutevalue{m^{2}_{3}-m^{2}_{1}}=\left[2.47-2.63\right]\times
10^{-3}\ e\text{V}^{2}$
$\displaystyle\sin^{2}\theta_{12}=\left[2.71-3.69\right]\times
10^{-1}\qquad\sin^{2}\theta_{23}=\left[4.34-6.1\right]\times 10^{-1}$
$\displaystyle\sin^{2}\theta_{13}=\left[2.000-2.405\right]\times
10^{-2}\qquad\delta=\left[0.71-1.99\right]\times\pi$ (52)

Here the lower and upper values define the $\pm 3\sigma$ range. The texture of
the neutrino mass matrix in our model is compatible with the observed mass
differences and mixing angles. To simplify the expressions, we consider
$m^{\nu}_{32}\to 0$, resulting in a neutrino mass matrix that reads

$M_{\nu}=\begin{pmatrix}-\frac{(m^{\nu}_{13})^{2}}{M^{N}_{33}}&-\frac{m^{\nu}_{22}m^{\nu}_{11}}{M^{N}_{12}}\left(1-\frac{m^{\nu}_{13}M^{N}_{13}}{m^{\nu}_{11}M^{N}_{33}}\right)&-\frac{m^{\nu}_{13}m^{\nu}_{33}}{M^{N}_{33}}\\\
\cdot&-\frac{(m^{\nu}_{22})^{2}(M^{N}_{13})^{2}}{M^{N}_{33}(M^{N}_{12})^{2}}&-\frac{(m^{\nu}_{22})(m^{\nu}_{33})(M^{N}_{13})}{M^{N}_{33}(M^{N}_{12})^{2}}\\\
\cdot&\cdot&-\frac{(m^{\nu}_{33})^{2}}{M^{N}_{33}}\end{pmatrix}\equiv\begin{pmatrix}\frac{a^{2}}{Z}&X&a\\\
\cdot&Y&b\\\ \cdot&\cdot&Z\end{pmatrix}\,.$

We find that with the choice of parameters

$\displaystyle a\sim 2\times 10^{-3}~{}\mathrm{eV},\;b\sim-2.75\times
10^{-2}~{}\mathrm{eV}\,,$ (53) $\displaystyle X\sim 9\times
10^{-3}~{}\mathrm{eV},\;Y\sim-2\times 10^{-2}~{}\mathrm{eV},\;Z\sim-2\times
10^{-2}~{}\mathrm{eV}\,,$ (54)

we are within the parameter ranges listed in Eq. (52), and predict the
neutrino masses

$\displaystyle m_{1}\sim 2\times 10^{-3}~{}\mathrm{eV},\quad m_{2}\sim 9\times
10^{-3}~{}\mathrm{eV},\quad m_{3}\sim 5\times 10^{-2}~{}\mathrm{eV}\,.$

Assuming that the sterile neutrinos have masses in the TeV range, which can be
probed at the LHC CMS:2018iaf ; ATLAS:2019kpx , the above values require Dirac
Yukawas that are $y^{\nu}\sim\order{10^{-7}}$.

##  4 Conclusions

We have proposed a model which extends the Standard Model gauge group by the
anomaly-free abelian group $U(1)^{\prime}\equiv U(1)_{L_{e}-L_{\mu}}$. Two
scalars that are charged under this $U(1)^{\prime}$ group are added to the SM
particle content, one $SU(2)_{L}$ doublet, $\phi$, and one singlet, $S$.

As a result of the spontaneous breaking of the $U(1)^{\prime}$ gauge group,
the $Z^{\prime}$ boson acquires a mass. Since the new scalar doublet is singly
charged under the $U(1)^{\prime}$ gauge group, $\mu\leftrightarrow\tau$ and
$\tau\leftrightarrow e$ Yukawa couplings are allowed, while
$\mu\leftrightarrow e$ couplings are forbidden. Nonetheless, $\mu\to e$ can be
mediated by the product of $\mu\to\tau\times\tau\to e$ interactions. After
electroweak symmetry breaking, the $Z^{\prime}$ couplings receive flavour off-
diagonal components due to the misalignment of the gauge eigenbasis and the
mass eigenbasis of the leptons.

In the case of $M_{Z^{\prime}}\gtrsim m_{\tau}$, $\tau\leftrightarrow\ell$
searches do not appreciably constrain the model and the $\mu\to e$ processes
lead to the most stringent limits, despite being proportional to the product
of Yukawa couplings $y_{23}\times y_{31}$. The model predicts rates that are
in reach of the upcoming $\mu\to e$ experiments and our $Z^{\prime}$ boson can
be searched for at Belle-II.

Finally, if we add sterile neutrinos that are charged under $U(1)^{\prime}$ to
the spectrum, the singlet VEV can contribute to their Majorana masses. We show
that, via a type-I seesaw mechanism, our model can accommodate for the
observed neutrino masses and mixing angles.

## Acknowledgments

We thank Sacha Davidson for her useful insights and feedback. M.A. is
supported by a doctoral fellowship from the IN2P3 and thanks the Physik-
Institut of the University of Zurich for its hospitality during the completion
of this work. F.K. acknowledges support by the grant PP00P2_176884 of the
Swiss National Science Foundation and thanks the Laboratoire Univers et
Particules de Montpellier for its hospitality. Last but not least, we would
like to thank Bella for her moral support.

##  Appendix A Barr-Zee Contributions

In this appendix, we give the most relevant two-loop Barr-Zee type
contributions to the $\mu\to e$ dipole for reference Chang:1993kw . The
corresponding diagrams are shown in Figures 5(b)-5(c).
The top loop contribution to the dipole is given by

$(C^{e\mu}_{D,L})_{\mathrm{t-loop}}=\frac{e\alpha}{12\pi^{3}}\frac{1}{m_{\mu}m_{t}}\left[\sum_{\Phi=h,\rho}F^{e\mu}_{L,\Phi}F^{tt}_{L,\Phi}f\left(\frac{m^{2}_{t}}{m^{2}_{\Phi}}\right)+F^{e\mu}_{L,A}F^{tt}_{L,A}g\left(\frac{m^{2}_{t}}{m^{2}_{\Phi}}\right)\right]$
(55)
$(C^{e\mu}_{D,L})_{\mathrm{t-loop}}=\frac{e\alpha}{12\pi^{3}}\frac{1}{m_{\mu}m_{t}}\left[\sum_{\Phi=h,\rho}F^{e\mu}_{R,\Phi}F^{tt}_{R,\Phi}f\left(\frac{m^{2}_{t}}{m^{2}_{\Phi}}\right)+F^{e\mu}_{R,A}F^{tt}_{R,A}g\left(\frac{m^{2}_{t}}{m^{2}_{\Phi}}\right)\right]\,,$
(56)

where (substituting $s_{\beta\alpha}\sim 1)$

$\displaystyle F^{e\mu}_{L,h}$ $\displaystyle=\frac{\mathcal{Y}^{*}_{\mu
e}}{\sqrt{2}}c_{\beta\alpha}\qquad
F^{e\mu}_{R,h}=\frac{\mathcal{Y}_{e\mu}}{\sqrt{2}}c_{\beta\alpha}\qquad
F^{tt}_{X,h}=\frac{m_{t}}{v}(1-c_{\beta\alpha}\tan\beta)$ (57) $\displaystyle
F^{e\mu}_{L,\rho}$ $\displaystyle=-\frac{\mathcal{Y}^{*}_{\mu
e}}{\sqrt{2}}\qquad
F^{e\mu}_{R,\rho}=-\frac{\mathcal{Y}_{e\mu}}{\sqrt{2}}\qquad
F^{tt}_{X,\rho}=\frac{m_{t}}{v}(c_{\beta\alpha}+\tan\beta)$ (58)
$\displaystyle F^{e\mu}_{L,A}$ $\displaystyle=-i\frac{\mathcal{Y}^{*}_{\mu
e}}{\sqrt{2}}\qquad F^{e\mu}_{R,A}=i\frac{\mathcal{Y}_{e\mu}}{\sqrt{2}}\qquad
F^{tt}_{X,A}=-i\frac{m_{t}}{v}\tan\beta\,.$ (59)

The loop functions $f$ and $g$ are defined as Chang:1993kw

$\displaystyle f(z)=$
$\displaystyle\frac{1}{2}z\int_{0}^{1}dx\frac{1-2x\left(1-x\right)}{x\left(1-x\right)z}\log\frac{x(1-x)}{z}$
(60) $\displaystyle g(z)=$
$\displaystyle\frac{1}{2}z\int_{0}^{1}dx\frac{1}{x(1-x)-z}\log\frac{x(1-x)}{z}\,.$
(61)

The $W$ loop of Figure 5(c) leads to the dipole contribution

$\displaystyle(C^{e\mu}_{D,L})_{\mathrm{W-loop}}$
$\displaystyle=\frac{e\alpha}{32\pi^{3}}c_{\beta\alpha}\frac{\mathcal{Y}^{*}_{\mu
e}}{\sqrt{2}m_{\mu}v}\left(3f\left(z_{h}\right)+5g\left(z_{h}\right)+\frac{3}{4}(g\left(z_{h}\right)+h\left(z_{h}\right))+\frac{f\left(z_{h}\right)-g\left(z_{h}\right)}{2z_{h}}\right)$
$\displaystyle-\frac{e\alpha}{32\pi^{3}}c_{\beta\alpha}\frac{\mathcal{Y}^{*}_{\mu
e}}{\sqrt{2}m_{\mu}v}\left(3f\left(z_{\rho}\right)+5g\left(z_{\rho}\right)+\frac{3}{4}(g\left(z_{\rho}\right)+h\left(z_{\rho}\right))+\frac{f\left(z_{\rho}\right)-g\left(z_{\rho}\right)}{2z_{\rho}}\right)$
$\displaystyle(C^{e\mu}_{D,R})_{\mathrm{W-loop}}$
$\displaystyle=\frac{e\alpha}{32\pi^{3}}c_{\beta\alpha}\frac{\mathcal{Y}_{e\mu}}{\sqrt{2}m_{\mu}v}\left(3f\left(z_{h}\right)+5g\left(z_{h}\right)+\frac{3}{4}(g\left(z_{h}\right)+h\left(z_{h}\right))+\frac{f\left(z_{h}\right)-g\left(z_{h}\right)}{2z_{h}}\right)$
$\displaystyle-\frac{e\alpha}{32\pi^{3}}c_{\beta\alpha}\frac{\mathcal{Y}_{e\mu}}{\sqrt{2}m_{\mu}v}\left(3f\left(z_{\rho}\right)+5g\left(z_{\rho}\right)+\frac{3}{4}(g\left(z_{\rho}\right)+h\left(z_{\rho}\right))+\frac{f\left(z_{\rho}\right)-g\left(z_{\rho}\right)}{2z_{\rho}}\right)\,,$

where we have defined $z_{\Phi}=m_{W}^{2}/m^{2}_{\Phi}$. The loop function $h$
is defined as Chang:1993kw

$\displaystyle h(z)=z^{2}\frac{\partial}{\partial
z}\left(\frac{g(z)}{z}\right)=\frac{z}{2}\int_{0}^{1}\frac{dx}{z-x(1-x)}\left(1+\frac{z}{z-x(1-x)}\log\frac{x(1-x)}{z}\right)\,.$
(62)

## References

  * (1) G. D’Ambrosio, G. F. Giudice, G. Isidori and A. Strumia, _Minimal flavor violation: An Effective field theory approach_ , _Nucl. Phys. B_ 645 (2002) 155 [hep-ph/0207036].
  * (2) MEG collaboration, _Search for the lepton flavour violating decay $\mu^{+}\rightarrow\mathrm{e}^{+}\gamma$ with the full dataset of the MEG experiment_, _Eur. Phys. J. C_ 76 (2016) 434 [1605.05081].
  * (3) SINDRUM collaboration, _Search for the Decay mu+ — $>$ e+ e+ e-_, _Nucl. Phys. B_ 299 (1988) 1.
  * (4) SINDRUM II collaboration, _A Search for muon to electron conversion in muonic gold_ , _Eur. Phys. J. C_ 47 (2006) 337.
  * (5) MEG II collaboration, _The design of the MEG II experiment_ , _Eur. Phys. J. C_ 78 (2018) 380 [1801.04688].
  * (6) A. Blondel et al., _Research Proposal for an Experiment to Search for the Decay $\mu\to eee$_, 1301.6113.
  * (7) COMET collaboration, _Conceptual design report for experimental search for lepton flavor violating mu- - e- conversion at sensitivity of 10**(-16) with a slow-extracted bunched proton beam (COMET)_ , .
  * (8) BaBar collaboration, _Searches for Lepton Flavor Violation in the Decays tau+- — $>$ e+- gamma and tau+- —$>$ mu+- gamma_, _Phys. Rev. Lett._ 104 (2010) 021802 [0908.2381].
  * (9) K. Hayasaka et al., _Search for Lepton Flavor Violating Tau Decays into Three Leptons with 719 Million Produced Tau+Tau- Pairs_ , _Phys. Lett. B_ 687 (2010) 139 [1001.3221].
  * (10) Belle-II collaboration, _The Belle II Physics Book_ , _PTEP_ 2019 (2019) 123C01 [1808.10567].
  * (11) P. Langacker, _The Physics of Heavy $Z^{\prime}$ Gauge Bosons_, _Rev. Mod. Phys._ 81 (2009) 1199 [0801.1345].
  * (12) X.-G. He, G. C. Joshi, H. Lew and R. R. Volkas, _Simplest Z-prime model_ , _Phys. Rev. D_ 44 (1991) 2118.
  * (13) E. D. Carlson, _Limits on a new u(1) coupling_ , _Nuclear Physics B_ 286 (1987) 378.
  * (14) P. Langacker and M. Plumacher, _Flavor changing effects in theories with a heavy $Z^{\prime}$ boson with family nonuniversal couplings_, _Phys. Rev. D_ 62 (2000) 013006 [hep-ph/0001204].
  * (15) B. Murakami, _The Impact of lepton flavor violating Z-prime bosons on muon g-2 and other muon observables_ , _Phys. Rev. D_ 65 (2002) 055003 [hep-ph/0110095].
  * (16) C.-W. Chiang, Y.-F. Lin and J. Tandean, _Probing Leptonic Interactions of a Family-Nonuniversal Z’ Boson_ , _JHEP_ 11 (2011) 083 [1108.3969].
  * (17) J. Heeck and W. Rodejohann, _Gauged $L_{\mu}-L_{\tau}$ Symmetry at the Electroweak Scale_, _Phys. Rev. D_ 84 (2011) 075007 [1107.5238].
  * (18) J. Heeck, M. Holthausen, W. Rodejohann and Y. Shimizu, _Higgs → $\mu$$\tau$ in Abelian and non-Abelian flavor symmetry models_, _Nucl. Phys. B_ 896 (2015) 281 [1412.3671].
  * (19) J. Heeck, _Lepton flavor violation with light vector bosons_ , _Phys. Lett. B_ 758 (2016) 101 [1602.03810].
  * (20) M. Ardu, S. Davidson and M. Gorbahn, _The sensitivity of $\mu\to e$ processes to $\tau$ flavour change_, 2202.09246.
  * (21) S. Davidson and H. E. Haber, _Basis-independent methods for the two-Higgs-doublet model_ , _Phys. Rev. D_ 72 (2005) 035004 [hep-ph/0504050].
  * (22) S. Davidson, _$\mu\rightarrow e\gamma$ in the 2HDM: an exercise in EFT_, _Eur. Phys. J. C_ 76 (2016) 258 [1601.01949].
  * (23) ATLAS collaboration, _Cross-section measurements of the Higgs boson decaying into a pair of $\tau$-leptons in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Rev. D_ 99 (2019) 072001 [1811.08856].
  * (24) CMS collaboration, _Observation of the Higgs boson decay to a pair of $\tau$ leptons with the CMS detector_, _Phys. Lett. B_ 779 (2018) 283 [1708.00373].
  * (25) M. B. Wise and Y. Zhang, _Lepton Flavorful Fifth Force and Depth-dependent Neutrino Matter Interactions_ , _JHEP_ 06 (2018) 053 [1803.00591].
  * (26) M. Bauer, P. Foldenauer and J. Jaeckel, _Hunting All the Hidden Photons_ , _JHEP_ 07 (2018) 094 [1803.05466].
  * (27) BaBar collaboration, _Search for Dimuon Decays of a Light Scalar Boson in Radiative Transitions Upsilon — $>$ gamma A0_, _Phys. Rev. Lett._ 103 (2009) 081803 [0905.4539].
  * (28) BaBar collaboration, _Search for a Dark Photon in $e^{+}e^{-}$ Collisions at BaBar_, _Phys. Rev. Lett._ 113 (2014) 201801 [1406.2980].
  * (29) S. Andreas, C. Niebuhr and A. Ringwald, _New Limits on Hidden Photons from Past Electron Beam Dumps_ , _Phys. Rev. D_ 86 (2012) 095019 [1209.6083].
  * (30) J. D. Bjorken, R. Essig, P. Schuster and N. Toro, _New Fixed-Target Experiments to Search for Dark Gauge Forces_ , _Phys. Rev. D_ 80 (2009) 075018 [0906.0580].
  * (31) S. Bilmis, I. Turan, T. M. Aliev, M. Deniz, L. Singh and H. T. Wong, _Constraints on Dark Photon from Neutrino-Electron Scattering Experiments_ , _Phys. Rev. D_ 92 (2015) 033009 [1502.07763].
  * (32) M. Lindner, F. S. Queiroz, W. Rodejohann and X.-J. Xu, _Neutrino-electron scattering: general constraints on $Z^{\prime}$ and dark photon models_, _JHEP_ 05 (2018) 098 [1803.00060].
  * (33) A. Kamada and H.-B. Yu, _Coherent Propagation of PeV Neutrinos and the Dip in the Neutrino Spectrum at IceCube_ , _Phys. Rev. D_ 92 (2015) 113004 [1504.00711].
  * (34) E. Rrapaj and S. Reddy, _Nucleon-nucleon bremsstrahlung of dark gauge bosons and revised supernova constraints_ , _Phys. Rev. C_ 94 (2016) 045805 [1511.09136].
  * (35) S. Knapen, T. Lin and K. M. Zurek, _Light Dark Matter: Models and Constraints_ , _Phys. Rev. D_ 96 (2017) 115021 [1709.07882].
  * (36) M. Aiba et al., _Science Case for the new High-Intensity Muon Beams HIMB at PSI_ , 2111.05788.
  * (37) Belle-II collaboration, _Belle II Technical Design Report_ , 1011.0352.
  * (38) G. Inguglia, _Belle II studies of missing energy decays and searches for dark photon production_ , _PoS_ DIS2016 (2016) 263 [1607.02089].
  * (39) LHCb collaboration, _Search for Dark Photons Produced in 13 TeV $pp$ Collisions_, _Phys. Rev. Lett._ 120 (2018) 061801 [1710.02867].
  * (40) P. Ilten, J. Thaler, M. Williams and W. Xue, _Dark photons from charm mesons at LHCb_ , _Phys. Rev. D_ 92 (2015) 115017 [1509.06765].
  * (41) A. Ibarra, M. Mar$́\mathrm{i}$n and P. Roig, _Flavor violating muon decay into an electron and a light gauge boson_ , _Phys. Lett. B_ 827 (2022) 136933 [2110.03737].
  * (42) CMS Collaboration collaboration, _Search for lepton-flavor violating decays of the higgs boson in the $\mu\tau$ and $e\tau$ final states in proton-proton collisions at $\sqrt{s}=13\text{ }\text{ }\mathrm{TeV}$_, _Phys. Rev. D_ 104 (2021) 032013.
  * (43) R. Harnik, J. Kopp and J. Zupan, _Flavor Violating Higgs Decays_ , _JHEP_ 03 (2013) 026 [1209.1397].
  * (44) D. Chang, W. S. Hou and W.-Y. Keung, _Two loop contributions of flavor changing neutral Higgs bosons to mu — $>$ e gamma_, _Phys. Rev. D_ 48 (1993) 217 [hep-ph/9302267].
  * (45) J. D. Bjorken and S. Weinberg, _Mechanism for nonconservation of muon number_ , _Phys. Rev. Lett._ 38 (1977) 622.
  * (46) 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 (2020) 115043 [2007.09612].
  * (47) S. Davidson, _Completeness and complementarity for $\mu\to e\gamma\mu\to e\bar{e}e$ and $\mu A\to eA$_, _JHEP_ 02 (2021) 172 [2010.00317].
  * (48) Mu2e collaboration, _Mu2e Technical Design Report_ , 1501.05241.
  * (49) P. Minkowski, _$\mu\to e\gamma$ at a rate of one out of 109 muon decays?_, _Physics Letters B_ 67 (1977) 421.
  * (50) P. F. de Salas, D. V. Forero, S. Gariazzo, P. Mart$́\mathrm{i}$nez-Miravé, O. Mena, C. A. Ternes et al., _2020 global reassessment of the neutrino oscillation picture_ , _JHEP_ 02 (2021) 071 [2006.11237].
  * (51) CMS collaboration, _Search for heavy neutral leptons in events with three charged leptons in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Rev. Lett._ 120 (2018) 221801 [1802.02965].
  * (52) ATLAS collaboration, _Search for heavy neutral leptons in decays of $W$ bosons produced in 13 TeV $pp$ collisions using prompt and displaced signatures with the ATLAS detector_, _JHEP_ 10 (2019) 265 [1905.09787].