Searching for Monopoles via Monopolium Multiphoton Decays

We explore the phenomenology of a model of monopolium based on an electromagnetic dual formulation of Zwanziger and lattice gauge theory. The monopole is assumed to have a finite-sized inner structure based on a 't Hooft-Polyakov like solution, with the magnetic charge uniformly distributed on the surface of a sphere. The monopole and anti-monopole potential becomes linear plus Coulomb outside the sphere, analogous to the Cornell potential utilised in the study of quarkonium states. Discovery of a resonance feature in the diphoton channel as well as in a higher multiplicity photon channel would be a smoking gun for the existence of monopoles within this monopolium construction, with the mass and bound state properties extractable. Utilising the current LHC results in the diphoton channel, constraints on the monopole mass are determined for a wide range of model parameters. These are compared to the most recent MoEDAL results and found to be significantly more stringent in certain parameter regions, providing strong motivation for exploring higher multiplicity photon final state searches.


Introduction
The existence of magnetic monopoles was first postulated by Dirac, with their existence found to provide an explanation for the quantisation of electric charge [1][2][3]. Monopoles have a rich phenomenology with connections to new physics scales as they are generic predictions of spontaneously broken Grand Unified Theories (GUT). Thus, discovering monopoles would provide a window through which to explore higher energy scales and new physics both terrestrially and cosmologically. Possible early universe phase transitions associated with them could also lead to gravitational waves, providing an additional avenue for discovery. There are important cosmological implications for monopoles due to their stability, ensured by magnetic charge conservation, such as electroweak Baryogenesis and Dark matter [4][5][6][7].
Monopoles have been extensively searched for experimentally, but have yet to be discovered [8,9]. One important way to search for monopoles is through production at colliders which are able to probe monopoles of low mass relative to that usually expected from GUT models, such as current efforts by CMS, ATLAS, and MoEDAL collaborations. No evidence of monopoles has been found to date, with the current collider lower bounds on the monopole mass reaching as high as a few TeV, with dependence upon the nature of the monopoles [10,11]. These direct production searches face many technical and theoretical difficulties due to monopoles strong couplings and their highly ionising nature.
One exciting prospect for observing monopoles is through monopole-anti-monopole bound states, known as monopolium, and their associated decay products [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. The monopolium decay to diphotons provides a clean signal region in which to observe possible decays, compared to the types of signals expected from monopole-anti-monopole pair production. Additionally, the potential large negative binding energy of the monopolium states opens the window to probing higher energy scales beyond that which is possible through direct production at colliders. The large binding energy means the production of monopolium is highly favoured over pair production, while allowing a probe of the higher energy physics scale of the constituent monopoles. Such searches can provide complementary avenues for discovery to those done at MoEDAL and other collider searches, and illuminate the dynamics of the monopoles and their possible bound state formation. The details of the monopolium bound state are not clearly understood because of the strong couplings and subsequent non-perturbative nature. Various analyses has been done, but many uncertainties remain with mainly qualitative results able to be gleaned [14][15][16]18]. By considering cross-section constraints in the diphoton channel, it may be possible to derive constraints on the monopolium states [29][30][31]. In future, the measurement of higher multiplicity photon final states may provide additional avenues for discovery because the Standard Model background is expected to be small, and the strong coupling of the monopoles to photons potentially leading to large signals.
In a previous work, we presented a possible method for constructing a model of a monopoleanti-monopole bound state [18]. The model treated the monopolium analogously to a quarkonium state, deriving a modified Cornell potential [32], while taking into account the possible shielding of the internal structure of the monopole at a finite scale. We utilise a strong coupling expansion in lattice gauge theory, and accordingly a linear term is included in the monopole-anti-monopole potential, which allows for large binding energies to be studied. The phenomenological aspect of our work was initially motivated by the tentative evidence for a new resonance in the diphoton channel and the possible bound state nature [33][34][35], and thus had a narrow phenomenological scope with the main focus the construction of the bound state. In this work, we make improvements on various aspects in the construction and investigation of the model, and provide a more thorough exploration of the possible collider phenomenology within this framework.
The paper is structured as follows; firstly, in Section 2 we briefly review the monopoleanti-monopole bound state model proposed in our previous work [18]. In Section 3, we will discuss the properties of the monopolium state, including the binding energy, bound state wave function, and the diphoton decay properties. Section 4 derives the expected diphoton cross section for the monopolium states at the LHC, including the possible contributions of higher energy eigenstates. The implications to the collider phenomenology upon including higher multiplicity photon decays is explored in Section 5. After which we present a discussion of the results and concluding remarks, including comments on future prospects.

Model of Monopolium Bound State
In our previous work, we introduced a model of monopolium based on the Zwanziger formalism. The resultant potential for the monopolium state, which we shall investigate in this work, is analogous to the Cornell potential utilised in studies of quarkonium states [32]. We assume that there is only one kind of spin 1/2 monopole M in our world, that is electrically neutral, having magnetic charge g. Within the Standard Model we have a number of magnetically neutral spin 1/2 fields with electric charge Q = ±e (e is the unit of the electron charge), and the monopole satisfies the Schwinger quantization condition with these charged particles. This gives The difficulty of highly charged particles is the estimation of the potential between them. One photon exchange is insufficient, and so we will use a lattice gauge theory approach with a finite lattice constant a, in which a strong coupling expansion is possible [36][37][38]. We will apply this to the magnetic U (1) part of the manifestly electromagnetic dual formulation of Zwanziger [39,40]. However, U (1) gauge theories are not well defined in lattice gauge theory, due to the existence of a first-order phase transition between the weak coupling perturbative region and the strong coupling confinement region, and so the continuum limit of a → 0 cannot be consistently taken [41]. On the other hand, the continuum limit is properly taken for SU (2) [42] and other asymptotically free gauge theories. So, we consider that a non-Abelian structure is revealed when we go inside the finite sized U (1) monopole, with the gauge group expected to be enhanced from U (1) to SU (2) as a → 0, and the 't Hooft-Polyakov-like structure will appear [43,44]. In the 't Hooft-Polyakov monopole, the U (1) magnetic charge is located at the origin as a point-like singularity. Instead, we consider that the U (1) magnetic charge is distributed non-locally, with magnetic charge distributed uniformly on the surface of a sphere with radius R, for which a solution exists. Inside the sphere (r < R) there is no magnetic force, and so the potential becomes flat, while the potential between the monopole and anti-monopole becomes linear plus Coulomb outside the sphere. Below we provide a summary of the details of the monopolium construction postulated in our previous work.

Zwanziger's manifestly dual formulation of gauge theory
The manifestly dual formulation of U (1) gauge theory by Zwanziger [39,40] is given by the following action, Here, a constant unit vector η µ , denoting the direction of Dirac strings [1][2][3], is displayed in parallel along the space-like direction, and where fermions with electric charge e i and magnetic charge g i are introduced, and ε 0123 = 1. From the consistency condition of the finite Lorentz transformation, the following Dirac, or Schwinger type, quantization condition [1][2][3] appears, where N ij is an integer. Here, we consider a flat space-time g µν = η µν = (1, −1, −1, −1) and ε 0123 = −ε 0123 = 1. In the Zwanziger formulation, the degrees of freedom are doubled by the introduction of electric and magnetic vector potentials, but are halved by the projection to the η µ direction. If the axial gauge is taken, no ghost fields appear, and the Feynman rules are obtained as follows, The kinetic terms of the gauge field are complicated and depend on η µ , but the 2 × 2 matrix form of the propagators is simple and satisfies, where V 1 µ = A µ , V 2 µ = B µ , andD abµν is the differential operator for the gauge fields in the action.

Monopolium Formulation
The bound state of the monopoles is formed by the exchange of magnetic photons γ M . However, the coupling of the monopole to the magnetic photon is very strong, so that we adopt a lattice gauge theory approach with a lattice constant a as a UV cutoff [36][37][38]. The space-time is considered to be a square lattice n = (n 0 , n 1 , n 2 , n 3 ), where n µ are integers, with a lattice constant a. The link variables U nμ are introduced as usual for the electric and magnetic photons A µ (x) and B µ (x), respectively, The minimum Wilson loop is given for the boundary curve C nµν ≡ ∂P nµν of the minimum region P nµν = rectangular (n, n +μ, n +μ +ν, n +ν, n). Then, the Zwanziger action can be written as the lattice gauge theory action in the Euclidean metric, where the dual Wilson loop reads In the action, we assume that the magnetic coupling g is strong, while the electric coupling e is weak and perturbative. So, in estimating the expectation value of the large Wilson loop W (B) [C], the strong coupling expansion is used [36][37][38]. We choose C to be a rectangle of length T in time and length r in the space-like direction µ; then we have where T · r denotes the minimum area of the rectangle C, so that the Wilson's area law is realized only if µ is in the η-direction. We define the potential between a heavy monopole and its anti-monopole separated by a distance r to be V (r). Then, the potential has a linear term in r, if the monopole and anti-monopole are separated in the η-direction; This clarifies the meaning of the special direction η that appears in the Zwanziger formulation, giving the direction of the Dirac string starting from the monopole. Therefore, the monopole and anti-monopole are connected by the string, starting from the monopole and ending at the anti-monopole, which contributes to the linear potential between them.
In addition to this strong coupling contribution, we will add the usual perturbative weak coupling contribution, and so the potential at this stage is This matches with high precision QCD calculations in which the potential between a quark and anti-quark pair is well approximated by the linear plus Coulomb potential [45][46][47]; this potential is known as the Cornell potential and is well-defined for point-like quarks. The monopole, however, may not be point-like, and may have an internal structure. The U (1) lattice gauge theory is usually considered not to be well defined, since there exists a firstorder phase transition between the confinement phase and the perturbative phase [41], and it obstructs the continuum limit of a → 0. One way out from this difficulty is to lift the U (1) theory to SU (2) theory, or other asymptotically free theory, when we approach the short distance region. As is shown by 't Hooft and Polyakov [43,44], the U (1) monopole was given as a classical solution of SU (2) gauge theory, which is broken to a U (1) by a triplet Higgs field φ a (x) (a = 1, 2, 3). Therefore, if we go inside the monopole, the non-Abelian gauge theory may appear, in which case, we may take the continuum limit properly.
The 't Hooft-Polyakov monopole is a classical solution of the SU (2) gauge theory with a triplet Higgs, based on, and the following ansatz: where we define a dimensionless parameter ξ = evr. If we take the limit λ → 0 while keeping v = 0, the solution, called the BPS solution, is given analytically [48,49]. Then, the equations of motion become first order, which have the following solutions, We want to know the distribution of the magnetic charge inside the monopole, ξ < 1 or r < 1/ev. The SU (2) gauge potential and the Higgs field are properly reduced by factors of 1 − K(ξ) and H(ξ), but the monopole charge is unfortunately not smeared even inside the monopole. This can be understood from the U (1) field strength proposed by 't Hooft. This gauge invariant expression can be rewritten as follows [50], whereφ a = φ a /|φ|. From this expression the magnetic charge is found to be a topological number, where Φ m is the magnetic flux, and N is the winding (wrapping) number of the sphere of the Higgs fields (φ 2 = v 2 ) by the sphere in space (x 2 = 1). We will consider the scenario in which the magnetic charge is distributed uniformly on the surface of the sphere, that is one of the solutions that exists. Let's take the following modified gauge and Higgs fields in which the contribution from r < R is cutoff, where θ(r − R) is a step function 0 for r < R, and 1 for r > R. Accordingly, K and H are modified: This modification inside the sphere (r < R) is allowed, since K(r) = 1 and H(r) = 0 is the solution.
The distribution of the U (1) magnetic charge is easily understood. Since φ a = 0 inside r < R, the magnetic flux on the surface of the radius r sphere reads, which shows that the magnetic charge 4π/e is distributed uniformly on the surface of the monopole sphere with radius R. The solution satisfies the following equations, which are modified only on the surface, where the magnetic charge is distributed. Therefore, the magnetic force between the monopole and anti-monopole vanishes when one goes inside either respective sphere (r < R). We will consider that the linear potential is the dominant contribution and the Coulomb potential plays the role of lowering the potential. Then, the corresponding potential is, where α g = g 2 4π = 1 α = 137 . We shall utilise this potential to describe the monopolium bound state in what follows.

Wave Function and Binding Energy of Monopole in Bound State
Using the potential V (r), we can obtain the wave function and energy eigenvalues of the monopolium states. Given that we have a piecewise potential, we solve the Schrödinger's equation requiring continuity at the point r = R in the radial wavefunction, where we will assume spherical symmetry. We define the potential as follows and derive the radial wavefunction for both potentials, where we denote the string tension by κ = ln(8πα g )/a 2 . The radial wavefunction is defined by ψ(r) = χ(r)/r, such that χ(r) satisfies the onedimensional Schrödinger equation. From the potentials given in Eqs. (32) and (33), the corresponding radial wavefunction solutions are described by sinusoidal and Airy functions, respectively, where the boundary conditions χ 1 (0) = 0 and lim r→∞ χ 2 (r) = 0 have been imposed, and the nth energy eigenvalue is, with corresponding eigenstate mass, for n ≥ 1 , where x n = x n (p, q, m) is a function of p, q and m, that is to be determined numerically.
In order to determine C 1 and C 2 , and the corresponding energy eigenvalues, we require that the following three conditions are satisfied, for continuity, and to normalise the wavefunction. These can be solved numerically, and the energy eigenvalues derived. By requiring 0 < M n < 2m, the list of possible spherically symmetric bound state can be compiled for given input parameters, R, a and m.
In the analysis that follows we will make the identification that R = q m and a = p m . The energy scales in this scenario are given by both the core size of the monopole, and necessarily when the QED description would break down embodied in the lattice scale a. Furthermore, motivated by the usual classical monopole radius we will assume that R and a are proportional to 1 m , with constants of proportionality defined as q and p [51]. The most motivated choice of q is 137 from the definition of the approximate monopole core radius, there may be some ambiguity in this denominator so we will allow it to vary for our analysis. In this parametrisation, the monopole mass parameter m factors out for many of the calculations, allowing determination of the monopolium properties by varying the parameters q and p, each exhibiting unique bound state signatures and phenomenology.
We consider that the monopole and anti-monopole collide when the distance r = R in our model, so that the wave function is evaluated at r = R.
The eigenvalue masses are determined by the value of x n , which increases with n. That is, for n ≥ 1, and subsequently the lower bound on the monopolium mass M * , for a given q, is given by, and is depicted as a function of q and p in Figure 1. For large p, the smallest allowed q value can be derived to be approximately 68.5 < q, above which M n > 0 is always satisfied. Although smaller q values can be accessed by choosing q > p, in what follows we explore q ≥ 70 to ensure that the monopolium masses considered always satisfy M 0 > 0. We will require 2m > M n > 0, for which the monopolium states exhibit negative binding energies, opening the possibility to probing energy scales much greater than production of the constituent particles would allow, particularly for large binding energy scenarios. In the following we will confine our analysis to q ∈ [70, 200] and p ∈ [60,2000], which takes into all major regions of interesting behaviour in this model. p values that are too small can lead to monopolium masses greater than 2m, while for q this would give M < 0. We do not consider combinations in which q is so small and p so large that M → 0, which we ensure be considering q > 70. The large negative binding energies we obtain allow us to probe monopole masses well beyond those probed by pair production. Now we can investigate the properties of the wavefunction for the corresponding energy eigenvalues. Upon numerically solving the above equations, the number of eigenvalues obtained in the region M 0 < M < 2m, for the range of q and p parameters to be considered, varies between 100-600 states. Below we will also consider the stability and decay properties of these different energy level states, and the phenomenological implications.

Diphoton decay width of Monopolium States
After obtaining the energy and wave function of the monopolium, we can estimate the decay properties of the monopolium states. The decay rate of of the monopolium state to diphotons can be estimated from the cross section σ(m + m → 2γ) as, where since the monopoles have a finite core size R, we use the average probability density of finding the monopole and anti-monopole at a relative distance R or less in the bound state. We assume that a distance of r ≤ R indicates collision of the monopole-anti-monopole pair, and hence annihilation. The factor 4 σv rel gives the reaction rate, with the factor 4 coming from the possible combination of monopole and anti-monopole spin states. In this calculation, the polarization vector of the magnetic photon m (k) µ , has to be converted to the usual electric photon's polarization γ (k) ν , by using the off-diagonal propagator The polarization sum of photons, pol γ (k) * λ γ (k) ρ = −g λρ , is performed, then we have γ s pol.
if the photon is on the mass shell, k 2 = 0, and the Ward identity is used, k µ and k ν terms vanish in the amplitude of monopoles on the mass shell, and hence γ s pol. m (k) * µ m (k) ν = −g µν . Therefore, the expression of σ(m + m → γγ) is obtained from σ(e + + e − → γγ), by replacing the mass m e by m and the coupling constant e by g.
Now we can obtain the full form of the two photon decay rate as, which is the formula valid at low energy, in the case of √ŝ = M ≈ 2m. This approximation will be used in our subsequent analysis.
Due to the large coupling associated with the monopole coupling, this decay width can be very large depending upon the relation between the average probability density and the monopolium mass. We will consider as a general requirement throughout our work that the total decay width is less than 10% of the corresponding monopolium mass, in order to not strongly violate the narrow width approximation.
In Figure 2 we depict the dependence of the two photon decay width to monopolium mass ratio for a range of p and q parameters. The decrease of this ratio with increasing q is attributable to the enlargement of the monopolium mass and decrease in the average probability density. Values of q 130 lead to large decay widths for p > 40, so we will consider q > 130 for the analysis of the lowest energy level monopolum states.
We must also consider the contribution to the decay rate coming from the emission of multi-photons due to the monopoles strong coupling to photons. The total decay rate will become, to calculate this total decay width we will follow [20], which utilises an analogy to multiphoton positronium decays to estimate the branching ratios of higher multiplicity photon final states.
In the next section we will confine ourselves to the two photon final state only, with the subsequent section introducing higher multiplicity photon decays and their collider phenomenological implications.

Production Cross Section and LHC Constraints on Monopoles from Monoplium Diphoton Decays
Now that we have obtained the states that may contribute to resonance like collider signatures, we can derive the current constraints on the monopolium model parameters and corresponding monopole mass. The current strongest cross section limits in the diphoton channel are provided by the ATLAS collaboration with ∼139 fb −1 of data [31], and CMS collaboration with ∼39 fb −1 at 13 TeV [30]. The derived monopole mass constraints can then be compared to those derived from the MoEDAL experiment which has been designed to search for pair produced monopoles. We will not consider states for which Γ M > 0.1M , as the smeared resonance structures exhibited by these states are weakly constrained by current experimental bounds. Despite this, for cases where the energy difference between the energy levels of the monopolium and the total decay widths are of comparable size, it may be important to consider possible interference effects. To calculate these effects would require analysis of the many possible states, which is beyond the scope of the current analysis. Our approximate assumption that these interference effects may not be detrimental will be discussed below.
In calculating the predicted diphoton cross sections, we will begin by following the works of [15] [16],σ where the total decay width includes the sum of all possible contributions, including multiphoton channels which will be discussed in the next section. Photon fusion processes are made up of three key processes, inelastic, semi-elastic and elastic, which must each be included to obtain the total predicted cross section. From [52] we can compute the inelastic contribution to the production cross-section, that in which neither of the protons remains intact after exchange of a photon. This is given by the following integral, where the energy of the subprocess is given byŝ = x 1 x 2 z 1 z 2 s, √ s th = M is the threshold centre of mass energy, and the factors F p 2 (x, Q 2 ) represent the deep inelastic proton structure function. The value of Q 2 is chosen throughout to beŝ/4. Subsequently, we use the NNPDF2.3QED package to compute the partonic densities in the proton, in order to calculate the cross section.
In the semi-elastic process, only one of the protons remains intact after photon exchange. The associated production cross section is given by, whereŝ = x 1 z 1 z 2 s. For the photon spectrum f γ (z), we use where Q 2 1 =ŝ/4 − M 2 /4 and Q 2 2 = 1 GeV 2 . To calculate the elastic photon spectrum we use an approximate analytic form utilised in [52], which are based upon the Weizsäcker-Williams approximation [53][54][55]. It is of the following form, where and where in the limit s m 2 p , we can approximate this as Finally, the elastic contribution, where both protons remain intact after exchanging a photon, is given by whereŝ = z 1 z 2 s. The total monopolium production cross section will be given by the sum of these three contributions. From this we can now begin our analysis of the diphoton channel predictions for our model at the LHC. To start we will consider the scenario in which the monopolium only has one decay path, into two photon, and consider only the lowest energy eigenvalue for given p and q values. Thus we use the two photon decay width given in Eq.

Predicted Cross Section to Diphoton Final State
We can now calculate the predicted diphoton cross sections for our monopolium model at the 13 TeV LHC searches. In calculating these we have included the elastic, inelastic and semi elastic contributions to the photon fusion processes. These results can then be compared to the latest constraints determined by the ATLAS and CMS collaborations, to see which parameter regions of our model have already been probed. The diphoton channel constraints we utilise are the 95% CL constraints from ATLAS and CMS, across the S = 0 state mass ranges [500, 3000] GeV and [2250, 4600] GeV, respectively. The range of masses considered by the ATLAS Collaboration, [500, 3000] GeV, is smaller than that given by CMS, [500, 4600] GeV. Although ATLAS utilises 139 fb −1 [31], compared to 39 fb −1 at 13 TeV [30] for CMS, hence providing more stringent constraints for the mass range [500, 3000] GeV. Thus, we include both analyses in order to extend the monopolium masses probed to a maximum of 4.6 TeV. Given the prevalence of large decay widths in our formulation, it is important to note the constraints applied by the ATLAS collaboration can probe widths as large as 10% of the monopolium mass, see Table 1 Table 1: Current ATLAS diphoton cross section constraints for both the Narrow Width Approximation (NWA), and broad width resonances [31].
The dependence of the predicted diphoton cross section on p and q is depicted in Figure 3. In each figure we fix q to different values and vary p within p ∈ [60,2000]. The variation of p seen here indicates a continuous decrease in the predicted cross section with decreasing p, and equivalently increasing monopolium mass. Larger p values decrease the monopolium mass until the p dependent term in the monopolium mass becomes negligible, after which further decreases in p have minimal effects. On the other hand, increasing q values lead to larger monopolium masses for given p and m.
By comparing different fixed q values we can see the monopole mass constraints across the range of possible p. The variation observed for the full range of parameters makes it difficult to apply explicit limits on the allowed monopole masses. In Figure 4, we show the allowed parameter regions for fixing only the monopole mass. The parameters p and q are varied across a range of values for which Γ 2γ < 0.1M , namely p ∈ [60,2000] and q ∈ [130,200]. Due to the continuous behaviour observed for varying p and m for a fixed q, we obtain a polygon which subtends the region of predicted cross sections for a given monopole mass, within our model.
In the context of our monopolium construction, we can thus constrain monopole masses to be greater than 1500 GeV for the entire model parameter space. Constraints can be applied on different sets of input parameters up to 4000 GeV. The most theoretically motivated choice of the q parameter, is likely q ∼ 137 due to the relation to the classical monopole radius. For this value of q the monopole mass can be constrained up to ∼ 3500 GeV, which is larger than the current constraints from the MoEDAL experiment.
The shape of the predicted cross section region is approximately conserved with changing monopole mass, which is expected from the factoring out of m dependence in the energy The latest 95% CL constraints from the CMS (purple) [30] and ATLAS (brown) [31] collaborations are included. eigenvalue calculation. The region is then defined by the range of p and q values that are taken. Once p q, there is no increase in the size of this region because the p contribution to the monopolium binding energy gets increasingly suppressed. An increase in the region is only possible through considering smaller p and q values, but these have minimum values fixed by the requirement that the monopolium mass satisfies 0 < M < 2m. In the next section, we include model parameters with smaller q values that still maintain M > 0, for which the lowest energy eigenvalues that exhibited large decay widths. We do this by considering higher energy levels of the bound states, which can have decay widths satisfying Γ < 0.1M .

Including Higher Energy Eigenstates
For the model parameters within the range 70 < q < 130 the lowest energy level consists of a total decay width that is too large to be probed by current collider results, which assume narrower resonances. In order to test this region of parameter space we consider the lowest energy state for which Γ M = Γ 2γ < 0.1M is satisfied. From the production cross section extracted from these energy eigenstates we can then attempt to constrain them.
The justification for this choice can be seen in Figure 5, which depicts all the energy eigenvalues with Γ M < 0.1M for the q = p = 100 model, and their predicted cross sections. It can be seen that the state with the largest cross section is that which has the lowest energy eigenvalue, and hence smallest monopolium mass, with all subsequent states are suppressed. This suppression is caused by the combination of larger monopolium masses, and the decrease of the probability density factor ρ for increasingly excited energy states. Although it may be possible to probe other states in this series, we will consider only the one with the largest cross section. One possible concern with this analysis is the potential for interference effects between the different energy eigenstates. Although this could be an important effect for the large decay widths considered here, for simplicity we assume the relative cross sections is sufficiently suppressed for the higher energy states that there contributions can be ignored. A full analysis of the possible implications of this interference is beyond the scope of the current work, with further investigation of interest in future work.
By opening our analysis to this higher energy states we can probe smaller values of q, namely we consider 70 < q < 130, in addition to n = 1 states for 130 < q < 200. As can be seen in Eq. (40), such model parameter selections leads to larger negative binding energies, and thus smaller monopolium masses for a given monopole mass. This gives an opportunity to probe even larger monopole masses than would otherwise be possible be monopole-anti-monopole pair production.
In Figures 6 and 7, we provide the expected cross sections for the monopolium diphoton decays for the model parameter range q ∈ [70, 200] and p ∈ [60,2000]. The monopole mass is chosen for various values between 1500 GeV and 14000 GeV for illustrative purposes. By  comparing these results to those in Figure 4, for the ground state case, we see that including excited states (n = 2, 3,...) has resulted in enlarging the predicted cross section region. This is due to the inclusion of monopolium models consisting of q < 130 parameter values, which were originally excluded as their ground state decay widths violated our requirement that Γ M < 0.1M . The excited energy states of q < 130 monopoliums exhibit larger cross sections and smaller bound state masses, than the ground states of q > 130 models. Furthermore, where in the n = 1 case, no parameter constraints could be found for monopole masses greater than 4000 GeV, in the n ≥ 1 case, it is possible to constrain monopole masses as large as 12000 GeV from Figure 7, for certain input parameter. This points to the key advantage and motivation of considering the formation of monopole-anti-monopole bound states.

Multi-photon Annihilation of Monopolium and Branching Ratio to Diphoton
In our analysis so far, we have neglected the possibility of other decay channels existing beyond the expected diphoton decays. Issues with monopoles arise when calculating the total decay width due to the non-perturbative nature of their couplings, with higher order photon decays can become enhanced leading to very large total widths. In our previous work, we attempted to control this divergent behaviour but the resultant total decay widths were well beyond the mass of the monopolium itself, making it difficult to derive any meaningful phenomenological constraints. Here we will consider the possibility of making an analogy between the monopolium multiphoton decays and those of the positronium bound states, which was recently suggested in Ref. [20]. By making this connection the divergence of the full width is constrained by phase space suppression effects. The model tends to lead to the prediction of significant multiphoton final states with multiplicity greater than two, as suggested in our previous work. This can lead to interesting phenomenology, as even if the production rate of monopolium states is large their decays may be hidden from the diphoton channel, with the majority of decays occurring into higher multiplicity final states, which are difficult to constrain at present colliders. One may be able to obtain tentative bounds on these from two photon plus missing energy searches, but a thorough analysis is needed. This is strongly motivated by the expected large signals that can be produced by monopolium decays. The branching ratio model suggested in Ref. [20] may be valued more on a qualitative than quantitative basis, but we will consider the implications for our monopolium construction here, under the assumption that this positronium analogy is approximately valid.

Branching Ratio to Higher Multiplicity Photon Decays
In [20], the ratios of the higher multiplicity photon final states were matched to those of the positronium states, to obtain a numerical approximation for higher n multiplicities. The resultant ratio of the rate of multiphoton emission is, where Γ 2nγ is the decay width of monopolium to n photon pairs, and α g = 1/α = 137 for monopoles of charge g. Note that, this approximate relation assumes that the monopole loop can be shrunk to a point due to the monopole mass being very heavy, and hence that the monopolium and multi-photons interact via a contact interaction. The total decay width is then given by the sum over all final state pairs of photons, which allows us to write, from which we can derive at what monopolium masses each higher multiplicity final state can be expected to dominate the total decay width, From Eqs. (57) and (60), we see that the decay to higher multiplicities n is favoured for larger monopolium masses, relative to the monopole mass. This has implications of varying importance for the diphoton branching ratio dependent upon the p and q values considered. Interestingly, this behaviour leads to an enhancement of the total decay width for increasingly higher order energy states, which is counter to the decrease observed in solely the two photon decay width case discussed in Section 4. This reduces the number of eigenstates that satisfy Γ M < 0.1M , while simultaneously suppressing the branching ratios of higher energy states to diphotons. The predicted diphoton cross section for the lowest energy state satisfying Γ M < 0.1M will thus dominate over the higher energy states.
In Figure 8, we can see how changing values of q and p effects the expected branching ratios to different multiplicity photon states. The key behaviour, is that for increasing values of q the decay width to two photons becomes rapidly suppressed. On the other hand, the dependence on p has a minor effect when p > q, with minimal effect for larger p values. This can be deduced by considering the importance of the monopolium to monopole mass term in the decay width ratio equation, Eq. (57), with the p and q dependence of the monopolium mass, see Eq. (40). Constraints from collider searches on higher multiplicity photon final states have not yet been explored, so monopolium mass states with smaller |∆E| are hidden from current experimental probes. To illuminate their properties, specialised analysis of higher multiplicity final states needs to take place. There may also be many regions of the model parameter space for which meaningful constraints could be provided by multiple decay channels, with a discover in two channels being a possible smoking gun for a monopolium state. This provides strong motivation for dedicated searches in the mulitphoton final states.
The inclusion of the extra possible decay channels means that the total decay widths are larger for all the monopolium states. In the previous section, we discussed the potential issues with constraining states exhibiting broad resonance structures. As such, in this scenario we will follow the methodology employed above in which we considered the higher energy eigenvalues, determining that state with the largest cross section while also maintaining Γ M < 0.1M . Unlike the case with solely diphoton decay, the higher energy eigenvalues now have rapidly increasing decay widths due to the monopolium mass dependence in the branching ratio calculation. This means that there will be less states that satisfy the requirement of a small total decay width, particularly for larger q models.
In Figure 9, we show the dependence of the predicted cross section on the energy eigenvalues for three different sets of parameters, q = 80, 100, 120, with p = 1000 and m = 2000 GeV fixed. This demonstrates the similar higher energy state cross section suppression to that depicted in Figure 5, but with fewer eigenvalues satisfying Γ tot < 0.1M due to the large decay width contribution to higher multiplicity decays. In what follows, we again assume that there is minimal interference effects between these states.

Constraints upon Inclusion of Multiphoton Decay Channels
In Figures 10 and 11, we provide the expected cross sections for the monopolium diphoton decays for the model parameter range q ∈ [70, 200] and p ∈ [60,2000]. The monopole mass is chosen for various values between 1500 GeV and 14000 GeV for illustrative purposes. The 95% CL diphoton channel constraints from the LHC collaborations, ATLAS and CMS, are included across the S = 0 state mass ranges [500, 3000] GeV and [2250, 4600] GeV, respectively. As in the previous scenarios, we can see that the monopole mass dependence has a small effect on the shape of the predicted cross section region, which instead is defined by the range of p and q values.
Upon comparing these results to those depicted in Figures 6 and 7, it can be seen that the region describing each monopole masses predicted cross section has been shifted, with the peak cross section remaining approximately the same, and the minimum being largely suppressed. Importantly, this leads to a narrowing of the predicted region, allowing easier separation of possible degeneracy between model parameters within uncertainty. Additionally, in allowing for higher multiplicity photon decays we can no longer fully rule out monopoles of mass 1500 GeV, as could be achieved when considering only diphoton decay. This is because those higher mass monopolium states now predominately decay into higher multiplicity states and as such are hidden from the diphoton channel. States that decay predominantly into higher multiplicity photon states require dedicated searches to be probed. This means their production by photon fusion is also suppressed. One may also consider multiphoton fusion processes to search for these particles, for example in heavy ion collisions.
Similar to the two photon only case, constraints can be applied on monopole mass models as high as 13 TeV, which is far beyond the current capabilities of dedicated monopole searches at the LHC. This is the key motivation of searching for monopoles through their bound state formation, and gives exciting possibilities for the monopole masses that could be probed by future colliders with larger centre of mass energies, such as 100 TeV.  GeV respectively, and the latest 95% CL LHC constraints from the CMS (red) [30] and ATLAS (green) [31] collaborations are included. The blue regions subtend the range of predicted cross sections for parameters p ∈ [60,2000] and q ∈ [70, 200], with (a) labelled with corresponding q values at the extrema, which are consistent across each figure.

Conclusions
The prospect of discovering the existence of monopoles could lead to many exciting phenomenological consequences, and also significant implications for cosmological observables, including Baryogenesis and gravitational waves. The strong coupling of monopoles leads to the possibility of monopole-anti-monopole pairs readily be produced in the form of bound states, monopolium. The large negative binding energy of monopolium states means that they can be copiously produced at energy scales well below that of the monopole mass itself. This can open the window to probing higher energy scales than could be reached by pair production alone. The monopolium states are expected to decay predominantly into multiphoton final states, which can provide a clean signal in comparison to that produced by monopole pair production. This possibility is a strong motivation for exploring further the physics of monopolium states.
In this work, we have made enhancements to and expanded upon the phenomenological details of a monopolium construction we had previously proposed. This monopolium model utilised an analogy to the strong coupling regime description of quarkonium bound states, in which the bound state is formed non-perturbatively using a strong coupling expansion in lattice gauge theory. Improvements have been made in the calculation of the properties of the bound state, and the description of higher energy eigenvalues of the monopolium state. A new method has been used to control the divergence of total decay width upon including multiphoton decays, by utilising an analogy to multiphoton positronium decays that was suggested by a recent work [20]. These improvements have allowed a thorough analysis of the phenomenology of our monopolium construction at current colliders. For certain parameter regions of our model, we have derived constraints on the spin 1/2 monopole mass that exceed the strongest constraints available from the MoEDAL collaboration, namely m > 2420 GeV [10].
Here we have constrained our analysis to the para-monopolium (S = 0) bound states with l = 0, but the existence of ortho-monopolium states (S = 1) is also expected. The signatures produced by decays to three photon final states would be an ideal complementary test to the para-monopolium decays, but at present there are minimal constraints in this decay channel [57]. We have also assumed that the possible interference effects between energy eigenstates are minimal, which will require a more detailed analysis beyond the scope of the current work. Large areas of parameter space have been neglected due to the very large decay widths that they exhibit, Γ tot > 0.1M . Future theoretical and experimental techniques may make it possible to more easily probe these states. For example, see the recent progress in searches for tt analyses which consider total decay widths as high as Γ tot ∼ 0.6M [58]. Prior to the discovery of the Higgs boson, there was discussion of the idea of the stealthy Higgs [59]. A large invisible Higgs width could lead to a very broad resonance allowing it to evade conventional methods. Such phenomenological implications considered therein could help guide how to consider such large multiphoton decay widths for the monopolium state.
The building of new colliders such as the proposed 100 TeV collider [60], will be able to probe deeply into the parameter space of our model, possibly leading to the discovery of monopoles with masses of order O(10) TeV. This in combination with new cosmological and gravitational wave experiments may provide new possibilities to search for the existence of monopoles.