Discovery prospects for long-lived multiply charged particles at the LHC

In this work, we aim to provide a comprehensive and largely model independent investigation on prospects to detect long-lived multiply charged particles at the LHC. We consider particles with spin 0 and $\frac{1}{2}$, with electric charges in range $1 \le |Q/e| \le 8$, which are singlet or triplet under $SU(3)_C$. Such particles might be produced as particle-antiparticle pairs and propagate through detectors, or form a positronium(quarkonium)-like bound state. We consider both possibilities and estimate lower mass bounds on new particles, that can be provided by ATLAS, CMS and MoEDAL experiments at the end of Run 3 and HL-LHC data taking periods. We find out that the sensitivities of ATLAS and CMS are generally stronger than those of MoEDAL at Run 3, while they may be competitive at HL-LHC for $3 \lesssim |Q/e| \lesssim 7$ for all types of long-lived particles we consider.


Introduction
One of the motivations behind the design, construction and operation of the Large Hadron Collider (LHC) [1] is the discovery of new particles beyond the Standard Model, or the (partial) exclusion of their potential existence. Priorities during the LHC Run 1 and Run 2 included scans for peaks in invariant-mass spectra, searches for excesses of missing transverse energy in events with various final states, and precision tests of Standard Model (SM) predictions. Apart from the discovery of a SM-like Higgs boson, no search has led to evidence of a new, still unobserved particle.
As conventional searches kept on showing no deviation from the SM predictions, the interest of the experimental -and theoretical -community started shifting towards less orthodox final states, such as those by long-lived particles a e-mail: r.maselek@uw.edu.pl (corresponding author) (LLPs) [2][3][4]. In theory, long lifetimes may be due to a symmetry (leading to stable particles), narrow mass splittings, small couplings, or a heavy mediator. Here, we focus on multiply charged particles that either are pair produced and leave a high ionising trace in the detectors, or they form a bound state which decays into two photons. In this article we are interested in the lifetime that is long enough such that the particles can be treated as stable within the detectors. For shorter lifetimes of the order of mm and cm, other exotic signatures, such as displaced vertices and disappearing tracks, are available. However, these are beyond the scope of the paper.
ATLAS and CMS experiments have been actively searching for multiply charged LLPs [5][6][7][8] with electric charge Q ranging 2 ≤ |Q/e| ≤ 7 with e being the electric charge of the proton. Their searches exploit the fact that multi-charged LLPs give rise to highly ionising tracks with anomalously large energy loss, d E/dx. So far, their analyses focus on colourless fermions and include only Drell-Yan production channels. As pointed out in Refs. [9,10], however, production channels involving photons in the initial state may become important for |Q| > 1e and should be included in the analysis. In this paper we recast the large d E/dx search [8] for various types of LLPs with spin 0 and 1 2 and with and without colour charges. Our analysis includes production channels with photonic initial states, and take two important effects into account. One is underestimation of the transverse momentum, p T , of charged tracks for |Q| > 1e. It is caused by the estimation of track p T 's based on the track curvatures (r ∝ Q/ p T ) assuming Q = 1e. This leads to a loss of efficiency in p T cuts. The second effect comes from the fact that produced charged particles lose kinetic energy while traveling inside the detector via strong electromagnetic interactions with detector materials. For large |Q|, this effect may be so strong that for a large fraction of events charged LLPs arrive at the muon system too late after the next bunch cross-ing. Since such events are rejected, this negatively affects the overall sensitivity.
Pair-produced charged particles can also form a bound state, B. Those bound states are generally short-lived and decay into the Standard Model particles. The best way to look for such bound states is to focus on the γ γ decay mode and search for a bump in diphoton invariant mass distributions. We compute the event rate of the pp → B → γ γ process for various types of charged LLPs and estimate the current limit as well as the projected sensitivities achievable at Run 3 and HL-LHC.
The MoEDAL experiment, which is designed to search for highly ionising avatars [11], has published results on isolated magnetic charges, namely magnetic monopoles [12][13][14]. Rather recently, it also constrained high electric charges either indirectly through dyons [15], which possess both electric and magnetic charge, or as highly electrically charged LLPs [16]. Its capability to probe heavy highly electrically charged LLPs has been demonstrated quantitatively in previous studies [17][18][19], which have focused on |Q| ≤ 4e. These studies indicated that MoEDAL's sensitivity is generally less significant than those obtained from large d E/dx analyses by ATLAS and CMS, apart from some exceptional cases [17]. This is largely because the data collected at MoEDAL is roughly 10 times smaller than those at ATLAS and CMS. Also, MoEDAL can register only slow-moving charged particles with β < 0.15 · |Q/e|, whilst those particles are often produced with relativistic velocities and the above condition is not satisfied in majority of events for |Q| 4e.
The situation may, however, be different for |Q| 4e. While the signal efficiency of large d E/dx analyses drops for larger |Q| due to the two aforementioned obstacles (underestimation of p T and late arrival time at the muon system due to the velocity loss), MoEDAL's detection rate only increases with |Q|. In this paper we estimate, for the first time, MoEDAL's sensitivity to multiply charged LLPs with |Q| > 4e including photonic initial states and compare the sensitivities obtained from ATLAS and CMS analyses at Run 3 and HL-LHC.
The paper is structured as follows. In the next section we define the relevant theoretical models for this study. In Sect. 3 the open production modes for charged LLPs are studied. Our prescription to estimate the effect of hadronisation for coloured LLPs is also outlined. In Sect. 4, we discuss the large d E/dx analysis and the methodology to recast the CMS analysis. In Sect. 5, we briefly describe the MoEDAL detector and outline the estimation of the MoEDAL's sensitivity to the charged LLPs. The calculation of the bound state production and decay are discussed in Sect. 6, where the sensitivity from the diphoton resonance searches is also studied. The results of our numerical studies are summarised in Sect. 7. Section 8 is devoted to conclusions.

Models
In this study four types of electrically charged long-lived particles are considered: (spin 0 and 1 2 ) × (colour-singlet and triplet). In order to write a simple gauge invariant Lagrangian, we introduce the corresponding quantum field, φ, for scalars and, ψ, for fermions, and postulate that they are SU (2) L singlet and have the hypercharge Y = Q/e, where Q is the electric charge measured in the proton charge unit, e. With these fields, we extend the Standard Model Lagrangian, L SM , as (2.1) for scalars and μ for colour-triplet fields with B μ and G a μ being the U (1) Y and SU (3) C gauge fields, respectively, and T a being the generator of SU (3) C . After the electroweak symmetry breaking, the interaction with the B μ field should be rewritten in terms of the photon ( A μ ) and Z -boson (Z μ ) fields as where θ W is the weak mixing angle. We denote the scalar (fermion) particle φ +Q (ψ +Q ) and its antiparticle φ −Q (ψ −Q ). We also write ξ +Q (ξ −Q ) to indicate the charged particle (antiparticle) when not specifying the spin. Without the terms denoted by · · · in Eqs. (2.1) and (2.2), these charged particles are absolutely stable. Since the presence of stable charged particles is problematic in cosmology, we postulate extra terms responsible for decays in · · · in the Lagrangian. In favour of model-independence of the analysis, we do not specify these terms as well as the details of the decay. We instead treat the lifetime of ξ ±Q as a free parameter when studying the sensitivity at the MoEDAL experiment. In estimating the sensitivities of ATLAS and CMS, we assume the charged particles are meta-stable, i.e. their typical decay length is much larger than the size of the detectors.
The models described above are implemented in MadGraph 5 [20] with the help of the FeynRules [21] package. MadGraph 5 is used for later analyses in calculating tree-level production cross sections and generating Monte Carlo events.

Open production mode
In the models described in the previous section, long-lived charged particles may be produced at the LHC in a particleantiparticle pair and travel through the detectors before decaying. We call this type of production process the open production mode. In this case, the production modes are similar to those for other exotic particles such as the magnetic monopoles [22].

Colour-singlet particles
For colour-singlet particles, the following processes are possible (see Fig. 1): • Drell-Yan: qq initial state with s-channel γ /Z exchange • Photon fusion: γ γ initial state with a t-channel interaction (and a 4-point interaction for scalars) As can be seen in the diagrams of Fig. 1, the production rate is proportional to Q 2 for the Drell-Yan qq initiated process, while the γ γ initiated process has the production rate proportional to Q 4 . Although the latter process is generally suppressed by the parton distribution function (PDF) of photons, it can be dominant for large values of |Q|. 1 In Fig. 2 we show the leading order cross sections of open production modes for colour-singlet scalars (left) and fermions (right). The cross sections are calculated with MadGraph 5 [20,23,24] using LUXqed PDF set (LUXqed17_plus_PDF4LHC15_nnlo_100) [25,26], which includes improved calculation for the photon PDF. We see that around m = 100 GeV the cross section varies more than two orders of magnitude from |Q/e| = 1 to 8. The impact of Q on the cross section is stronger for larger masses and for scalars since the photon fusion channel is relatively more important than the Drell-Yan channel for those cases. For scalars with m = 2 TeV, the cross section with Q = 8e is more than three orders of magnitude larger than that with Q = 1e. For Q 1e, we expect that the electromagnetic higher order correction is sizeable. However, the next-to-leading order cross sections for multicharged particles have not been well studied and including this effect is beyond the scope of this study. Figure 3 shows the relative contributions from the Drell-Yan and photon fusion channels to the total cross section for scalars (left) and fermions (right) with m = 1 TeV. We see that the photon fusion is subdominant for Q = 1e due to the photon PDF suppression. Moving to larger Q, the relative importance of photon fusion rapidly grows sincê σ γ γ /σ qq ∝ Q 2 as mentioned above, and becomes dominant for Q/e 3 (5) for scalars (fermions). This growth is faster for scalars due to the following reason. The Drell-Yan production (with off-shell s-channel γ /Z exchange) from the qq initial state necessarily has J = 1 (p-wave). In order for scalars to have non-zero cross sections, the final state configuration must have a non-zero angular momentum, which indicates the cross section must be suppressed by the positive power of the velocity of the final state particles at the centre of mass frame. This p-wave suppression is absent for the photon fusion channel, which makes this contribution relatively more important and quickly dominates over the Drell-Yan for the scalar case. The p-wave suppression is absent in the fermion production both for Drell-Yan and photon fusion channels.

Colour-triplet particles and a hadronisation model
For colour-triplet particles, the following processes contribute to the open production: • Drell-Yan: qq initial state with off-shell g/γ /Z exchange in s-channel • gluon fusion: gg initial state with QCD interaction • gluon-photon fusion: gγ initial state with mixed QCD and QED t-channel interaction • photon fusion: γ γ initial state with pure QED t-channel interaction (and the 4-point interaction for scalars) Figure 4 shows the leading order cross section for the open production channel pp → ξ +Q ξ −Q for colour-triplet scalars (left) and fermions (right) as a function of the mass. Compared with the colour-singlet case, the Q-dependence of the cross section is much weaker particularly for small values of m. This is because the production channels with QCD interaction (Drell-Yan with gluon s-channel exchange, gg and gγ fusions) take up a large fraction of the total cross section for all values of Q/e (up to 8), as can be seen in Fig. 5, and they are independent or only weakly dependent on Q. The QCD interaction is relatively more important for lighter ξ ±Q due to the enhancement of gluon PDF at a small x (momentum fraction) region. As mentioned above, electromagnetic higher order corrections for Q > 1e particles are not available to date. For consistency, we use the leading-order cross sections also for the QCD productions throughout this paper.
If colour-triplet particles have a sufficiently long lifetime, they draw quark-antiquark pairs from the vacuum so that they dress them and hadronise into colour-singlet states. The electric charge of the colour-singlet heavy hadrons Q gets modified by Q from the original charge Q of ξ +Q due to the extra quarks; Q = Q + Q. We estimate this effect using a simple hadronisation model, assuming that a colour-triplet particle goes into one of the mesonic states with probability k and one of the baryonic states otherwise. The model assigns an equal probability to every state among the mesonic (baryonic) states. The colour-singlet Colourless Fermion   Tables 1 and  2 for a scalar (φ +Q ) and a fermion (ψ +Q ) particle, respectively, together with the charge shift Q and the probability p assigned to each state. In estimating the sensitivity of the ATLAS, CMS and MoEDAL experiments, we vary the free parameter k from 0.3 to 0.7 for a ball-park estimation of the uncertainty originated from the hadronisation effect.

Detecting open production at ATLAS and CMS
Experiments with general purpose detectors, ATLAS and CMS, have been actively looking for long-lived charged particles. Their concrete analyses differ in the details but the general strategy is the same. The most stringent and general constraint on the open production mode comes from so-called "large d E/dx" searches, which look for tracks with anoma- Contributions of various initial states, qq (red), γ γ (green), gg (yellow), γ g (blue), to the production rate of coloured scalars (left) and fermions (right) with m = 1 TeV lously high energy deposition in the inner tracker. ATLAS analysed their full Run 1 data (8 TeV, 20.3 fb −1 ) and interpreted their result for spin-1 2 particles with 2 ≤ Q/e ≤ 6 [5]. This search was updated with the 13 TeV 36.1 fb −1 data in Ref. [6], where the result was interpreted for fermions with Q/e ranging from 2 to 7. CMS performed their analysis on a combined dataset of 7 and 8 TeV pp collisions with integrated luminosity of 5.0 and 18.8 fb −1 , respectively [7]. The most recent result using the 13 TeV, 2.5 fb −1 data has also been published [8]. The results of these analyses are interpreted for charged fermions with |Q/e| ≤ 2 and long-lived supersymmetric particles.
For multiply charged particles with Q/e ≥ 2, the above ATLAS and CMS analyses have performed only for colourless fermions. Moreover, their analysis considered only Drell-Yan production channels. Reference [10] has recast the CMS analysis [7] for coloured scalars and fermions and improved the original analysis by including the photonfusion 2 as well as gluon-gluon and gluon-photon initiated processes. In this article we follow the work of Ref. [10] and extend it by adding colourless scalars and estimation of projected sensitivities for Run 3 and HL-LHC.
As outlined in Ref. [10], our recasting is based on the CMS analyses [7,8], in particular the "Tracker + Time of Flight (TOF)" selection. The signal region is comprised of several kinematical cuts, such as p T and η of charged tracks, as well as more subtle requirements, e.g. conditions on the number of pixel hits and the uncertainty of 1/β measurements. We observe that the class of former cuts is sensitive to the spin and production mechanisms and its efficiency may be estimated by Monte Carlo (MC) simulations. On the other hand, the latter depends more crucially on the masses and charges of the particles. In order to obtain the final signal efficiency, , we first estimate the efficiency of the online (muon trigger) selection, on , by MC simulations. We then take the cut-flow tables for the offline selection provided in Ref. [27] for different masses and charges. The offline efficiency is given by the product of relative efficiencies of the cuts with respect to the previous ones in the cut-flow table; off (m, Q) = i i (m, Q), where i ≡ (# of events passed all cuts up to i)/(# of events passed all cuts up to i − 1). 3 We estimate the relative efficiencies of kinematical cuts by MC and replace the values in the cut-flow table. Since the relative efficiencies of non-kinematical cuts are not very sensitive to the spin, production mechanism and the collider energy, we use the same values provided in Ref. [27]. In the end the final signal efficiency is given by = on · off .
There are several caveats in estimating the efficiencies. In the ATLAS and CMS analyses the momentum of charged particles is estimated from the curvature of the track assuming |Q| = 1e. Since the curvature radius is proportional to p T /Q of the particle, this leads to underestimation of the momentum of the charged particle by a factor 1/|Q/e|. For example, the online selection requires at least one muon candidate with the measured transverse momentum larger than 50 GeV [8]. This translates to the condition p T > 50·|Q/e| GeV with p T being the true transverse momentum of the charged particle. Due to this effect, the efficiency of the p T cut tends to be smaller for particles with larger |Q|.
There is an implicit assumption in the muon trigger requirement; the track measured in the inner detector (ID) must match the one measured in the muon system (MS) in the same bunch crossing. If highly charged particles are produced and travel through the detector, they lose their velocity due to interactions with detector materials, and may arrive at the MS too late after the next bunch crossing. Following [10], we denote the minimum distance to which a particle must travel so that the momentum is measured at the MS, as x trigger , as a function of η. We then calculate, for each produced charged particle in the simulation, the amount of time needed to travel to the distance x trigger and denote it as t TOF . This calculation is performed using the procedure outlined in the appendix of [10], which is based on the Bethe-Bloch formula [28] and the information of the detector geometry and materials. Since the bunch crossing at the LHC occurs every 25 ns, we demand t TOF − x trigger /c < 25 ns to pass the muon trigger, where the left-hand-side represents the time separation between the detection of the massive charged particle and that of massless particles in the same bunch crossing. For charged particles with |η| < 1.6, the resistive plate chamber (RPC) muon trigger is applied, in which the threshold of time separation is relaxed to be 50 ns. Notice that the energy loss in the Bethe-Bloch formula is proportional to Q 2 . We therefore expect that the impact of this condition is more prominent for particles with larger |Q|.
In summary, for the online selection (muon trigger) we require at least one candidate with We then estimate and replace relative efficiencies of the following kinematical cuts in the offline selection: where β MS is the velocity measured at the MS. As mentioned above, for the relative efficiencies of non-kinematical cuts we use the values reported in the cut-flow tables in [27]. These tables provide the information only up to m = 1 TeV. Since all i (m, Q) of non-kinematical cuts, i, reach plateaus before m = 1 TeV, we use i (1 TeV, Q) for particles heavier than 1 TeV.
In Fig. 6 we show the signal efficiency, = on · off , as a function of Q. In the figure, the left (right) panels are for scalar (fermionic) particles, and the top (bottom) panels correspond to the colourless (coloured) particles. The curves with different colours represent different masses of the charged LLPs. For a fixed mass, we see a general tendency that the efficiency increases as moving from |Q/e| = 1 to 2. This is because charged tracks with |Q/e| = 2 deposit more energies in the detector than |Q/e| = 1 and lead to a distinctive signature with large d E/dx. However, the efficiency peaks around |Q/e| = 2 − 3 and decreases for larger |Q|. This can be understood as for very large |Q| (|Q/e| 4) the p T cut ( p T > 65|Q/e|) requires the particle to have a significantly large velocity, which in turn makes it difficult to pass the 1/β MS > 1.25 cut. Since the velocity is lower for larger masses at a given momentum, this also explains why the efficiency is higher for larger masses for the large Q region. On the other hand, we see that at |Q/e| = 1 the efficiency is higher for m = 500 GeV than heavier masses apart from the colourless scalar. We checked that this is because the efficiency of the charged track isolation criterion is generally smaller for heavier particles. As we will see later, the production velocity of colourless scalars is generally higher compared to the other types of particles. The efficiency is therefore still smaller for m = 500 GeV than heavier masses due to the requirement 1/β MS > 1. 25. In the upper panels of Fig. 7 we show the effective cross sections, σ eff ≡ σ · , of various colourless particles (left; scalars, right; fermions). As can be seen, the effective cross section increases significantly for larger Q. This is because the photon-fusion process is very important for colourless particles, and its production rate grows with Q 4 .
The lower panels of Fig. 7 show the effective cross sections of colour-triplet particles (left; scalars, right; fermions). We see that the σ eff is less sensitive to Q. This is because the production is dominated by the QCD Drell-Yan, gluongluon and gluon-photon initiated processes, of which the contribution is either independent of Q or grows with Q 2 for the gluon-photon fusion. We note that for m 800 GeV the effective cross section drops as a function of Q. This behaviour is due to the fact that the signal efficiency decreases for larger Q (Q 3e), and this effect is stronger for lower masses, as shown in Fig. 6.
The plots in Fig. 7 also show the (projected) modelindependent upper limits, σ UL eff , on the effective cross section. The dashed green line corresponds to the limit from the latest 13 TeV L = 2.5 fb −1 analysis [8], recasted to include all four types of the studied particles and production modes involving initial state photons. The gray and black dashed lines correspond to limits scaled to L = 300fb −1 (Run-3) and 3000 fb −1 (HL-LHC), respectively. The mass region with σ eff ≥ σ UL eff is sensitive (excluded) in the future (current) analyses.

Detecting open production at MoEDAL
In this section we outline the details of the MoEDAL detector and describe the procedure to estimate the expected number of the signal events of the open production mode observed at each stage of the LHC run by MoEDAL.

MoEDAL detector
The Monopole and Exotics Detector at the LHC (MoEDAL) [29,30] is a unique and mostly passive experiment located at the Intersection Point 8 (IP8), in the vicinity of the VErtex LOcator (VELO), in the LHCb detector cavern of the LHC. MoEDAL comprises three main sub-detector systems: nuclear track detectors, magnetic monopole trappers, and TimePix, each following a different principle of operation.
The main sub-detector is a large array of nuclear track detectors (NTDs) surrounding IP8. Each panel comprises of three layers of CR39 , three layers of Makrofol , and three layers of Lexan polymers, with 1.5 mm, 0.50 mm and 0.25 mm thickness respectively. All layers are placed inside 125 m thick aluminium bags. Currently, the Lexan foil is used only for protection, while six other polymer layers are used for detection of highly-ionising particles, such as magnetic monopoles [16]. When such particles traverse an NTD panel, they deposit energy in plastic material leading to invisible damage along particles' trajectories. After one year of exposure, NTD panels are transported to INFN Bologna, where they are dismantled and placed in an etching solution. The solvent dissolves polymer material anisotropically, leading to formation of cone-shaped etch-pits revealing a trajectory of a highly-ionising particle. After etching, NTD panels are scanned with an optical microscope in order to find etch-pits aligned in all sheets of polymer material. Signature of a signal is a string of etch-pits present in all 6 analysed layers of polymer from a given NTD panel, where all etch-pits have a shape of a cone, and are present on both sides of plastic layers. Careful analysis of signal pattern allows to assess the charge of the highly ionising particle, and can be used to discover multiply charged long-lived BSM particles.
One drawback of MoEDAL's LLP searches with respect to those at large experiments, i.e. ATLAS and CMS, is that the data available at MoEDAL is roughly 10 times smaller than that collected at the large experiments. This is because MoEDAL shares the same interaction point with the LHCb Another property of the MoEDAL experiment is its sensitivity to slowly moving particles. While large experiments are sensitive only to relativistic particles moving typically with β > 0.6, MoEDAL can detect particles with β < 0.15 · |Q/e|. Therefore, if a particle has charge Q = −3e, then it can be detected by MoEDAL only if it moves slower than 45% of speed of light. Also, the NTD panels are placed at the distance of ∼ 2 m away from the interaction point. This means that the MoEDAL is sensitive only to the charged long-lived particles with the lifetime typically τ 2[m]/c with c being the speed of light.
This leads to two interesting consequences: (a) MoEDAL provides a complementary detection phase-space when compared to large LHC experiments; (b) the acceptance and sensitivity of MoEDAL to BSM particles grows with the magnitude of electric charge. It is expected that due to large mass and significant energy loss caused by ionisation, new multiply charged BSM particles will have lower velocities than their SM counterparts. Therefore, it is possible that despite lower luminosity, MoEDAL can exhibit a comparable, or in some cases even better, sensitivity with respect to ATLAS and CMS, when searching for particles with high electric charges. Hence, we include MoEDAL in our analysis and compare its predicted sensitivity with that of large LHC experiments.
Another important property that differentiates MoEDAL from other LHC experiments is the fact that it is basically background free. Firstly, MoEDAL can directly detect only highly ionising particles. Secondly, heavy charged particles in the Standard Model are unstable and decay before reaching MoEDAL's NTD panels placed approximately two meters away from IP8. Thirdly, in the Standard Model metastable particles that can reach the NTD panels have much larger velocities and cannot satisfy the detection criteria, β < 0.15 · |Q/e|. Moreover, background radiation from the cavern surrounding the detector is actively monitored by TimePix devices. Finally, possible spallation or neutron recoil products will not produce a detectable signature, which corresponds to a particle going through whole NTD panel without much stopping. It is therefore evident, that MoEDAL can be treated as a background-free experiment. In this study we follow the approach developed in [17][18][19], where it has been proven that MoEDAL has a large potential for discovering long-lived BSM particles. The number of expected events observable at MoEDAL with the data corresponding to the integrated luminosity, L, is given by where σ is the production cross section for a pair of new particles. The efficiency, , in the above formula represents the average number of particles per event that hit one of MoEDAL's NTD panels with β < 0.15 · |Q/e|. This can be estimated with the MC simulation as where β i is the velocity of the i-th BSM particle produced in an event, (x) is the Heavyside step function ( (x) = 1 for x > 0 and 0 otherwise) and · · · MC represents the Monte Carlo average. β thr Q ≡ 0.15 · |Q/e| is the threshold value below which the particle can be registered by the MoEDAL detector. Finally, P( β i ) is the probability that the particle reaches an NTD panel, and given by where δ( β i ) is 1 if there is an NTD panel on the trajectory of the particle i and 0 otherwise, l( β i ) is the distance between the interaction point and the NTD panel along the particle's trajectory, γ ≡ 1/ 1 − β 2 i is the Lorentzian factor and τ is the lifetime of the particle. In our numerical analysis, we calculate the average in Eq. (5.2) over the Monte Carlo events generated with MadGraph 5 and evaluate the functions and P using the fast simulation code developed for [17][18][19], which contains the information of the location of NTD panels.
A couple of comments are in order regarding the efficiency calculation. First, precisely speaking, the detection efficiency depends also on the incident angle between the particle and the NTD panel. However, during Run-3 data taking period, NTD panels are placed directly facing the interaction point. In this situation incident angles cannot be large and this effect is negligible. Secondly, the condition β i < 0.15 · |Q/e| is trivially satisfied for |Q/e| ≥ 1 0.15 = 20 3 since β ≤ 1 by definition. For |Q/e| ≥ 20 3 MoEDAL can detect charged particles at any velocities as long as they hit and penetrate an NTD panel. Lastly, we assume that the charged particles traverse the LHCb VELO detector without much change of four-momentum, which is a good approximation for the range of charges we consider in this analysis. For extremely larger charges |Q/e| 8, ionisation loss becomes significant and more precise treatment is necessary to accurately estimate the efficiency.
To see how the condition β < 0.15 · |Q/e| impacts the signal acceptance of different LLPs, we show β distributions of the produced LLPs in Fig. 8. In the left (right) plot we take m = 500 (2000) GeV and the solid (dashed) histograms correspond to Q = 1e (8e). Four types of particles are shown: colourless scalar (red), colourless fermion (green) colourtriplet scalar (blue) colour-triplet fermion (yellow). Comparing the left and right plots, we see the production velocities are generally larger for smaller m. We also see that the production velocities of the colourless scalar with Q = 1 are, on average, significantly larger than those for the other types. This is because for this type of particles the dominant production channel is the qq initiated Drell-Yan, whose production rate vanishes in the β → 0 limit due to angular conservation. For Q = 8, as we have seen in Fig. 3, the dominant channel is replaced by the photon-photon fusion, in which such a property is absent. Apart from the Q = 1e colourless scalar, the velocity distributions are similar among all types of particles. The vertical dashed lines in Fig. 8 represent MoEDAL's detection threshold, β thr Q = 0.15 · |Q/e|. We see that only a small fraction of LLPs satisfy β < β thr Q for Q = 1e, while majority of produced LLPs can be registered by MoEDAL for Q = (5-6) e, depending on the masses and types of the particles.
As explained in Sect. 5.1, MoEDAL is a background free experiment, hence we consider the thresholds corresponding to N sig = 1, 2 and 3 when presenting the sensitivity. In Appendix B we show contours of those thresholds in the (m vs cτ ) planes for each type of particle with L = 30 and 300 fb −1 for Run 3 and HL-LHC, respectively.

Bound state and diphoton decay mode
For sufficiently large electric charges, the Coulomb attraction force between the ξ +Q -ξ −Q pair becomes important. At the vicinity of the threshold, the ξ +Q -ξ −Q pair can form a bound state B with the mass M B 2m, as depicted in Fig. 9. We call this production channel the closed production mode.
Depending on the initial partons, the bound state B can have different spins. For example, if the initial state is qq (or gg for coloured ξ ±Q ), the s-channel production via the off-shell γ /Z (or g) exchange allows B to have spin 1. On the other hand, the γ γ (or gg for coloured ξ ±Q ) initial state allows B to have spin 0 or 2, through the diagrams with the t-channel or the 4-point interactions. The spin of the bound state also dictates the decay. The spin 0 and 2 bound states can decay into the γ γ , Z γ and Z Z final states, as well as the gg final state if ξ ±Q is colourtriplet. On the other hand, the spin-1 bound state can decay into W + W − or any SM fermion-antifermion pair via the s-channel γ /Z (g) exchange. Among these final states, the diphoton channel, pp → B → γ γ , has by far the strongest experimental sensitivity and we therefore concentrate on this channel in this study to constrain the open production mode [10].
If the elementary particle ξ ±Q is colour-triplet, the ξ +Qξ −Q bound state may be colour-singlet or octet, since 3⊗3 = 1 ⊕ 8. However, at least QCD is concerned, the octet configuration does not lead to an attractive force, which renders the production rate subdominant compared to the production of colour-singlet states. In any case, the octet state cannot decay into the diphoton final state and we therefore do not consider it in this study.
The ξ +Q -ξ −Q system near the threshold can be treated in a good approximation by the non-relativistic quantum mechanics. The calculation of wave functions of energy eigenstates goes along the same line with the calculation for the hydrogen atom. Letting r (r = | r|) be the displacement from the position of ξ +Q to that of ξ −Q , the static potential between the two particles is given by where α s (r −1 ) is the running strong coupling evaluated at the scale r −1 , α is the structure constant of electromagnetism, C = (C 1 + C 2 − C B )/2 and C 1/2 and C B are the quadratic Casimir of SU (3) C for the two constituent particles and the bound state, respectively (C = 4/3 (0) for the colour-triplet (singlet) case). Using the reduced mass, , the Schrödinger equation can be brought into the same form as that of the hydrogen atom. The wave function nlm ( r ) of various energy eigenstates can also be found analogously to the calculation of the hydrogen atom. The decay rate of the bound state is proportional to the probability of finding the two constituent particles at the same point, | nlm (0)| 2 . This is non vanishing only for the s-wave.
The wave function ( n ≡ n00 ) at r = 0 is given by where r b is the "Bohr radius" of this bound state system and α s ≡ α s (r −1 rms ), for ground state r rms = √ 3 r b . With these inputs the probability of finding the two particles at the same point can be obtained as We see that the contributions from the radial excitations are suppressed by 1/n 3 . We therefore include only the ground state (n = 1) contribution in our analysis, (0) = 1 (0). The partial decay rates for B → γ γ and B → gg are given by where n c is the number of colours (n c = 1 for singlet and 3 for triplet) and n f is the dimension of Lorentz representation (n f = 1 for scalar and 2 for fermion). Since the gg final state is not available when ξ ±Q is colour-singlet, we set B→gg = 0 in colour-singlet case. The partonic cross section for the bound state production is related to the corresponding decay width aŝ whereŝ is the partonic centre of mass energy, c ab is a symmetric factor (c ab = 1 for a = b and 2 for a = b), D p is the dimension of the colour representation of particle p and J B is the spin of the bound state. The hadronic production cross section σ pp→B for the spin 0 bound state can be obtained by convoluting the partonic cross sectionσ ab→B (sτ ) with the luminosity function dL ab /dτ as where s is the pp collision energy and τ ≡ŝ/s. We carry out this integral with the Mathematica package ManeParse [31] using LUXqed PDF set (LUXqed17_plus_PDF4LHC15 _nnlo_100) [25,26]. Finally, the signal rate (cross section times branching ratio) for the pp → B → γ γ process can be obtained using the narrow width approximation as and (6.9)

Constraints from diphoton resonance searches
The signal rate of the diphoton channel from the bound state decay obtained in the previous subsection should be confronted with the experimental results of diphoton resonance searches conducted by ATLAS and CMS. The latest CMS result [32] is based on the 13 TeV collision data corresponding to the integrated luminosity of 36 fb −1 . On the other hand, the recent ATLAS analysis [33] uses 139 fb −1 of 13 TeV data. Since the latter analysis is based on a larger data set, we take the ATLAS analysis [33] and derive the current limit on the bound state production and project it to the Run 3 (300 fb −1 ) and HL-LHC (3000 fb −1 ). In Ref. [33], ATLAS studied the diphoton invariant mass distribution and did not find significant deviation from the Standard Model prediction. The result was interpreted for the spin 0 and 2 resonances for several values of decay widths and upper limits are placed on the production cross section times branching ratio to two photons as a function of the resonance mass. In our simple models, we have B /M B 2 % in all region of interest in the Q vs M B plane, which justifies to use the ATLAS upper limit derived using the narrow width approximation.
In Fig. 10 we show the signal cross section, σ pp→B→γ γ , for the diphoton channel as a function of the bound state mass M B = 2m for various types of particles. The upper (lower) panels correspond to the colour-singlet (triplet) ξ ±Q , and the left (right) panels represent the scalar φ ±Q (fermion ψ ±Q ) scenarios. The red dashed line corresponds to the 95 % CL upper limit obtained by ATLAS using 139 fb −1 data of 13 TeV pp collision [33] using the narrow width approximation. In this ATLAS analysis the limit is shown only up to M B = 2.5 TeV. In order to constrain the bound states with masses larger than 2.5 TeV, we adopt the same signal cross section upper limit obtained at M B = 2.5 TeV. This is justified since ATLAS observed no event with the diphoton invariant mass larger than 2.5 TeV and in this region the Standard Model background is below one event.
To estimate the projected limits for Run 3 (L Run3 = 300 fb −1 ) and HL-LHC (L HL = 3000 fb −1 ), we rescale the current upper limit σ UL 0 at L 0 = 139 fb −1 as σ UL Run3 = σ UL 0 √ L 0 /L Run3 for Run-3 and σ UL HL = σ UL 0 √ L 0 /L HL for HL-LHC, assuming the signal/background efficiencies and the systematic uncertainty do not drastically change. Since improving signal-background separation or reducing systematic uncertainties would lead to a better sensitivity, this approach is conservative. The projected signal cross section limits (95% CL) for Run 3 and HL-LHC are shown in the black and blue dashed lines in Fig. 10.

Results
In this section we collect the results of calculations outlined in the previous sections and compare the expected mass bounds for the meta-stable multi-charged particles obtained from the large d E/dx and diphoton resonance searches at ATLAS and CMS as well as from the etch-pits analysis at MoEDAL assuming the data that will be delivered at Run-3 and HL-LHC. For large d E/dx and diphoton resonance searches, Run 3 and HL-LHC correspond to the integrated luminosity of 300 and 3000 fb −1 , respectively, and we present the 95% CL limits, while for MoEDAL Run 3 and HL-LHC correspond to L = 30 and 300 fb −1 , respectively, and we show the expected signal detection N sig = 1, 2 and 3, throughout this section. We start by showing the current mass limits (at 95% CL) as functions of the electric charge in Fig. 11, where the top (bottom) panels show colourless (coloured) particles and the left (right) panels represent the scalars (fermions). The light red area enclosed by the two red curves is the region excluded by the latest CMS result of the large d E/dx analysis with the 13 TeV 2.5 fb −1 data [8]. The mass limits are re-calculated for various multi-charged particles, including the contribution of initial states involving photons. The lower boundary of the exclusion comes from the fact that for larger |Q| the p T cut ( p T > 65 · |Q/e| GeV) makes it very difficult to satisfy the slow moving requirement, 1/β MS > 1.25, for lighter particles. We see, however, that most of the region below the lower exclusion boundary is already excluded by the ATLAS diphoton resonance analysis (blue curve). We also believe that the remaining mass region below the lower boundary is excluded by the previous large d E/dx analyses with the 7 and 8 TeV data [7], where the p T cut is milder. We do not find any excluded region for colourless scalars with |Q/e| ≥ 6. This is because the production velocity of scalar particles is generally large and almost all events are rejected by the 1/β MS > 1.25 requirement. Only for colourless fermions (top right panel) we show by the purple curve the mass limits reported in the ATLAS large d E/dx analysis with the 13 TeV 36.1 fb −1 data [6]. ATLAS has interpreted their result only for colourless fermions with |Q|/e from 2 to 7. We note that those limits are computed based on leading-order cross sections excluding the photonic initial states. Although the amount of data used in the ATLAS analysis is larger than the latest CMS analysis [8], we cannot reinterpret this analysis since the signal efficiency of each cut is not publicly available.
The blue curves correspond to the mass limits obtained by recasting the ATLAS diphoton resonance analysis (13 TeV, 139 fb −1 ) [33] for the bound state decay into diphoton. We can see that for smaller charges (Q (2 − 4)e), the large d E/dx search gives tighter mass limits. Increasing |Q|, the mass bounds from the diphoton resonance rapidly grows and surpass the large d E/dx limits around Q ∼ (2 − 4)e, depending on the type of particles. In fact, as can be seen in Eqs. (6.3)-(6.5), the production cross section of the bound state grows as |Q| 10 , which is responsible for the steep rise of the mass limits as moving towards large values of |Q|. On the other hand, the limits from the large d E/dx search only mildly increase. Although the cross section of open production grows with the charge, the signal acceptance decreases as we have seen in Fig. 6, due to the underestimation of the charged track p T s and the velocity loss due to the strong elec-tromagnetic interaction with detector materials as discussed in Sect. 4. Figure 12 shows the projected mass limits for Run 3. The red and blue curves correspond to the 95% CL limits from the large d E/dx and diphoton resonance searches, respectively, while the green curves represent the expected signal detection at the MoEDAL experiment, corresponding to N sig = 1, 2 and 3. We do not show the lower boundary of the large d E/dx sensitivity since the region below it is already covered (or excluded) by the other analyses, e.g. the MoEDAL and ATLAS diphoton analyses in the plot as well as the analyses elsewhere with the Run-1 and Run-2 data. For coloured scalars and fermions, we vary the parameter, k, in our hadronisation model from 0.3 to 0.7. This in principle impacts the sensitivities of open production searches, i.e. the large d E/dx and MoEDAL analyses, since the electric charge of the final state particles are probabilistically shifted. This effect is visible in the MoEDAL's sensitivity for smaller |Q| (|Q| 3e) as a width of the green curves. However, for ATLAS γγ resonance (13 TeV, 139fb −1 ) Fig. 11 The mass limits on the multi-charged LLPs: colourless scalars (top left), colourless fermions (top right), coloured scalars (bottom left) and coloured fermions (bottom right); recasted to take into consideration photon-induced production processes. The region between the two red curves is excluded by the large d E/dx analysis with the 13 TeV 2.5 fb −1 data [8]. The purple curve in the top right plot is the mass limit reported in the ATLAS large d E/dx analysis with the 13 TeV 36.1 fb −1 data [6]. The region below the blue curve is excluded by the ATLAS diphoton resonance analysis (13 TeV, 139 fb −1 ) [33] larger charges and for the large d E/dx analysis, the effect is very small and almost invisible in the plots. This assures that the uncertainty coming from the hadronisation model is small and the mass limits are robust against this uncertainty. When extrapolating the mass bounds of large d E/dx analysis from L = 2.5 to 300 fb −1 , we follow the procedure described in [34,35] and obtained the 95% CL upper limit on the detected signal. In this procedure we assume the systematic uncertainty stays the same level (∼ 66%) relative to the total background. Our projected limits are conservative in this regard. The grey shaded regions show the currently excluded regions already presented in Fig. 11. We see that the Run-3 diphoton limit improves the current limit only mildly. This is because the current limit is calculated already with the 139 fb −1 data. This may be contrasted to the large d E/dx limits where the improvement is significant, since the current limit from the large d E/dx search in Fig. 11 is based on the 2.5 fb −1 data except for the colourless fermions with 2 ≤ |Q|/e ≤ 7, for which the limit is based on the 36.1 fb −1 data. As discussed above, the large d E/dx analysis constrains smaller charges well, while the diphoton resonance search is very constraining for large values of |Q|. At Run 3 the sensitivities from these searches become comparable for |Q| ∼ (5−6)e. The sensitivities of the MoEDAL experiment typically lie between the aforementioned two constraints. As can be seen, MoEDAL's Run-3 sensitivity with N sig = 1 may be as strong as those of the large d E/dx and diphoton resonance searches for |Q| ∼ (4-6) e.
We finally show the expected sensitivities at the HL-LHC in Fig. 13. We see that the improvement of the large d E/dx limits are very mild. We employ the same event selection used in [8], where CMS estimated the background contribution to the relevant signal region to be 0.06 with L = 2.5 fb −1 . For the HL-LHC, however, the background contami-  nation to the same signal region is expected to be 72 events, which requires a large number of detected signal events for exclusion. We expect that the limit from large d E/dx analyses can be further improved by optimising the event selection, which is beyond the scope of this paper. Compared to the large d E/dx and diphoton resonance searches, the improvement of the MoEDAL's mass limits from Run-3 to HL-LHC is significant. This is largely due to the background free nature of the MoEDAL experiment. While for ATLAS and CMS experiments the increase of luminosity results in larger number of expected signal and background events, only the expected signal rate grows for MoEDAL. We see that for 3 |Q|/e 7 MoEDAL may provide superior sensitivities than the other searches by ATLAS and CMS.

Conclusion
In this paper we have studied the prospect of observing long-lived particles with high electric charges, ranging 1 ≤ |Q/e| ≤ 8, at the LHC. In particular, we compared Run-3 and HL-LHC sensitives of three different types of searches: (a) anomalous track searches (i.e. large d E/dx analyses) at ATLAS and CMS, (b) charged LLP searches at MoEDAL and (c) diphoton resonance searches at ATLAS and CMS. The former two target the open production, pp → ξ +Q ξ −Q , with detection of charged LLPs within the detectors, whereas the latter concerns the production of the bound state followed by the decay into two photons, For the open production mode, we demonstrated the importance of the photonic initial states for the cross section, which have been omitted in the current ATLAS and CMS large d E/dx analyses. Inclusion of photonic initial states  is also crucial to correctly estimate the MoEDAL's detection sensitivity, since it has a large impact on the production velocities of the charged particles especially for scalars. In order to estimate the current sensitivities of the anomalous track searches and extrapolate them to the future luminosities, we recast the CMS large d E/dx analysis [8] for various charged LLPs including the photonic initial states. We also studied the effects of underestimation of the charged track p T as well as the velocity loss due to the (strong) electromagnetic interaction with the detector materials. Since these two effects grow with |Q| and negatively impact the signal efficiency, we found that the large d E/dx is most powerful for smaller |Q|.
On the other hand, the sensitivity of the diphoton resonance search grows rapidly with |Q|, since the event rate of pp → B → γ γ scales as |Q| 10 .
We found the MoEDAL's sensitivity is somewhere between those of the anomalous track searches and the dipho-ton resonance searches. For Run-3 MoEDAL may be competitive at best around |Q| ∼ (5 − 6)e. Due to the background free nature of MoEDAL, its sensitivity grows faster relative to that of the anomalous track and diphoton resonance searches as the luminosity increases. At the HL-LHC, therefore, MoEDAL may give superior sensitivities than ATLAS and CMS for 3 |Q| 6.
Although our conclusion is qualitatively robust, the mass bounds are derived based on the leading-order cross sections, since the electromagnetic higher order corrections for Q > 1e particles have not been calculated. Also, our estimate of Run-3 and HL-LHC sensitivities from the large d E/dx and diphoton resonance searches are based on the event selection employed in Run-2 and the conservative assumption on the systematic uncertainty. We hope this study encourages more thorough studies on the highly charged LLPs searches at ATLAS and CMS as well as the MoEDAL experiments. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.

Appendix A: Effect of Z-boson induced processes
In Sect. 3, we showed that the production of multi-charged particles from the γ γ initial state has cross sections proportional to Q 4 , hence becomes important for large Q. Since the production rate of the Z -boson initiated processes is also proportional to Q 4 , we perform a rough estimate to understand the significance of the Z -boson initiated processes. Since the Z -boson is much heavier than the proton, the parton distribution functions (PDFs) used for LHC analyses do not include the Z -boson as a parton. 4 At LHC energies, the Z -boson initiated processes should be better understood as the Z /γ fusion processes, qq → qqφ +Q φ −Q or qq → qqψ +Q ψ −Q . One of such diagrams is depicted in Fig. 14.
To evaluate the significance of the Z -boson contribution to the Z /γ fusion process with respect to the γ γ fusion, we compare the two cross sections computed with MadGraph 5 by the following commands: σ inc : generate u u > u u phiQ 4 Electroweak bosons (Z , W ± and h) should be included in the PDFs for future high energy hadron colliders with √ s ∼ 100 TeV. See e.g. [36].
The QCD=0 option prohibits QCD vertices from appearing in the diagrams. With the first command, σ inc represents the cross section of the uu → uuφ +Q φ −Q process including both Z and γ in the intermediate state. On the other hand, the second command asks for the same process but excludes Zbosons anywhere in the diagrams. σ exc therefore corresponds to the cross section to which only the γ γ fusion contributes. From these cross sections, we construct σ inc − σ exc σ exc and use it to measure the significance of the Z -boson contribution with respect to the γ γ fusion. Table 3 shows (σ inc − σ exc )/σ exc calculated for Q = 8e scalars and fermions at several mass points. For some cases we find σ inc < σ exc due to destructive interference. As can be seen, in the mass range we are interested, (σ inc −σ exc )/σ exc 15 %, which suggests that the Z -boson contribution to the Z /γ fusion process is subdominant compared to the γ γ fusion. This justifies our approximation to neglect the Zboson initiated productions in the analysis.