Construction of a Model of Monopolium and its Search via Multi-Photon Channels at LHC

A model of monopolium is constructed based on an electromagnetic dual formulation of Zwanziger and lattice gauge theory. To cope with the strong coupling nature of the magnetic charge, for which the monopole is confined, $ U(1) $ lattice gauge theory is applied. The monopole is assumed to have a finite-sized inner structure based on a 't Hooft-Polyakov like solution in which the magnetic charge is uniformly distributed on the surface of a sphere. The monopole and antimonopole potential becomes linear plus Coulomb outside the sphere and is constant inside. Numerical estimation gives two kinds of solutions: One which has a small binding energy, and hence the para-($ J=0 $) and ortho-($ J=1 $) monopoliums have degenerate masses. For the parameter choices considered, they both have $ \mathcal{O}(1-10) $ TeV masses and are very short-lived. The other solution has a small monopole mass and large binding energy, with an illustrative example of parameter choices giving a 750 GeV para-monopolium and 1.4 TeV ortho-monopolium. The production rate of the former is one order of magnitude smaller than the announced enhancement, but they may be the target of future LHC searches and the 100 TeV colliders.


Introduction
Last year, an interesting enhancement was reported at 750 GeV by the ATLAS [1] and CMS [2] groups at LHC. 1 The enhancement X was observed in the γγ channel and was not observed in the other di-boson channels such as ZZ, W + and W − : σ(p + p → X → γγ) = 10 ± 3 fb (ATLAS) and 6 ± 3 fb (CMS), (1) at 13 TeV (Run2) experiment, but no such enhancement was reported at 8 TeV (Run1). Upon combining the ATLAS and CMS data and including 8 TeV data the following cross-section estimate is obtained [3]: σ(pp → X → γγ) ≈ (4.6 ± 1.2) fb. (2) The total width Γ tot has not been well determined, but appears to be large ∼ 45 GeV. Various models have been proposed for X [4], but if X couples exclusively to γγ with a large coupling constant, then it may be possible to consider the enhancement as a monopolium, or a monopole-antimonopole dumbbell as described by Nambu [5]. A recent paper by M. Yamada et al. [6] considers a similar idea, but they identify the enhancement as a magnetic Higgs particle where it is dual to the electric Higgs particle in a hidden U(1) H gauge theory. Ours is more directly a bound state system of a monopole and an antimonopole. The detection of monopolium has already been studied by L. I. Epele et al., aiming for observation at Tevatron and LHC [7,8]. The main difference between their paper and the following analysis is that we adopt a strong coupling expansion in lattice gauge theory and accordingly a linear term is included in the monopoleantimonopole potential. Therefore, the large binding energy case can be studied. This work is motivated by papers [9,10] written by one of the authors (N. D. B.) on explaining the diphoton excess via photon fusion production of X, in particular where X is identified as a leptonium of highly charged leptons having a charge Q = (5 − 7) e. Higher charge can be naturally understood in the context of monopoles due to the large magnetic charge, which is opposite to the small electric charge. Also, the stability of the monopole is less crucial than the lepton, since our monopole is confined and can't exist individually.
The difficulty of highly charged particles is the estimation of the potential between them. One photon exchange is not enough, and so we will use lattice gauge theory with a finite lattice constant a, in which a strong coupling expansion is possible [11,12,13]. We will apply this to the magnetic U(1) part of the manifestly electromagnetic dual formulation of Zwanziger [14,15]. In lattice gauge theory, however, there exists another difficulty, that is, U(1) gauge theories are not well defined; a first-order phase transition exists between the weak coupling perturbative region and the strong coupling confinement region, and so the continuum limit of a → 0 can't be taken properly [19]. On the other hand, the continuum limit is properly taken for SU(2) [20] and other asymptotically free gauge theories. So, we consider that the welcoming non-Abelian structure reveals when we go inside the finite sized U(1) monopole, that is, at some point of taking a → 0, the gauge group is expected to be enhanced from U(1) to SU(2) and the 't Hooft-Polyakov [29,30,31] like inside structure will appear. Even in the 't Hooft-Polyakov monopole, the U(1) magnetic charge is unfortunately located at the origin as a point-like singularity. To relax this situation, we look for a solution in which the U(1) magnetic charge is distributed nonlocally. Fortunately, there exists a solution in which the magnetic charge is distributed uniformly on the surface of a sphere with radius R. Inside the sphere (r < R) there is no magnetic force, and so the potential becomes flat, while the potential between the monopole and antimonopole becomes linear plus Coulomb outside the sphere.
We estimate the monopolium masses and wave functions using this potential. In this model, the monopole and antimonopole are confined like quark and antiquark, so that they are not observed individually.
If the monopolium is found and the world of monopoles is opened in reality, then various interesting phenomena are expected. One of these is that the neutrino mass as well as possible lepton and baryon number violation effects may be attributed to the monopolium or the monopole-antimonopole dumb bell system [32].

Zwanziger's manifestly dual formulation of gauge theory
The manifestly dual formulation of U(1) gauge theory by Zwanziger [14,15] is given by the following action, Here, a constant unit vector η µ , denoting the direction of Dirac strings [33,34,35] 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. We can add scalars if necessary. From the consistency condition of the finite Lorentz transformation, the following Dirac (or correctly Schwinger type) quantization condition [33,34,35] appears where N ij is an integer. Here, we consider a flat space-time g µν = η µν = (1, −1, −1, −1) and ε 0123 = −ε 0123 = 1. If we want to obtain the dual gravity theory, then consider that the metric g µν (x) is formed as a collective excitation of a pair of gauge bosons, [36]. 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. Such a formulation in a non-Abelian gauge theory is more complicated but is possible [16,17,18]. 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 satisfieŝ where V 1 µ = A µ , V 2 µ = B µ , andD abµν is the differential operator for the gauge fields in the action. 2

Production and decay of monopolium at LHC
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. We have, of course, a number of magnetically neutral spin 1/2 fields having 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 g = 4π e N (N = integer). (13) 2 Propagators that Zwanziger derived have no η 2 kµkν (k·η) 2 term in the numerator [15]. These propagators are for a gauge-fixing Lagrangian L G = − 1 2 (∂ µ (η · A)) 2 + (∂ µ (η · B)) 2 . Two kinds of propagators with or without η 2 kµkν (k·η) 2 term are, of course, equivalent for the on-shell amplitudes, if the Ward identity is taken into account.
We denote the mass of the monopole simply as m, and assume that the monopole M and its antimonopole M form a J = 0 bound state X 0 = (MM ) ( 1 S 0 , paramonopolium) and a J = 1 bound state of X 1 = (MM ) ( 3 S 1 , ortho-monopolium) via the exchange of magnetic photons, which we will denote as γ M and are described by the B µ fields. The mass of the monopolium so obtained will be denoted generally as m X , or specifically m X 0 and m X 1 .
After obtaining the energy and wave function of the monopolium, we can estimate the various partial decay widths Γ(X 0 → f ) for X 0 or Γ(X 1 → f ) for X 1 , and the total decay widths Γ tot (X 0 → all f ′ s) or Γ tot (X 1 → all f ′ s), where f denotes a final state.
In this section we use perturbation theory even for the strong coupling monopole. Our strategy is to store the effects into the bound state problem or the estimation of the monopole and antimonopole potential, and study the multi-photon emission processes perturbatively. This is because, if we fix the number of emitted photons to n, then even if photons convert to magnetic photons and couple to the monopolium strongly, the order of coupling is O(g n ), so that perturbative calculation at this order can be applied. This is valid, however, under the restriction that radiative corrections due to monopole loops are ignored.

Production cross section
At LHC, to produce X we consider that a parton a inside one proton collides with parton b inside another proton to produce X as a resonance scattering, where the parton can be anything produced from the proton, namely, quark, gluon, or photon. To do this we introduce the parton distribution function inside the proton f a/p (x), where x is the momentum fraction. Then the production cross section reads Here, the patrons are labeled a and b, andσ is the cross section of the two colliding patrons, a and b, with the center of mass (CM) energy √ŝ , τ ≡ x 1 x 2 , andŝ = τ s, where √ s is the CM energy of the colliding protons.
If we introduce the "parton luminosity" for a and b as dL ab (τ )/dτ , and define it by [21] where c a and c b are the color degrees of freedom for a and b (c a = 3 and 8 for quark and gluon, respectively), then, the cross section becomes If we replaceσ by the resonance scattering formula, or its small width approximation 3 , From these, we have where the factor F = 1 for a = b, but is 1 2 for a = b, since a + b is originally the initial state, but it is converted to the final state in the formula. Now the production cross section is expressed as the product of two factors, the first is given by the parton luminosity function and s −1 , while the second one is a model-dependent number determined by the decay widths and mass of the resonance X. If we understand that s −1 = (13 TeV) −2 and (8 TeV) −2 correspond to 2.4 pb and 6.2 pb, respectively, in order to obtain 3 -6 fb (to be consistent with Eq. (2)), we have to find the suppression factor of 10 −3 from the product of dL ab (τ ) dτ , Γ(X → a + b)/m X , and Br(X → f ). This is one of the reasons why it is difficult to understand the 750 GeV enhancement.

Decay rate
The decay rate of X 0 → 2γ can be estimated from σ 0 (M + M → 2γ) and the value of the wave function of relative motion of the monopole and antimonopole at the distance R, ψ1 S 0 (r = R), that is since our model of the monopole has a finite size R, |ψ1 S 0 (R)| 2 gives the probability of finding the monopole and antimonopole at a relative distance R in the bound state, where the collision occurs. The factor 4 σ 0 v rel gives the reaction rate (the factor 4 comes from the possible combination of monopole and antimonopole spin states). Similarly, the decay rate of X 1 → 3γ is estimated from σ 1 (M + M → 3γ) and the value at r = R of the wave function of the ortho-monopolium ψ3 S 1 (r = R), that is To estimate σ 0 (M +M → 2γ) or σ 1 (M +M → 3γ), we first consider the process of emission of two or three magnetic photons, and convert the emitted magnetic photons to the electric photons. The estimation is done by using the perturbation theory, since the final number of emitted photons is restricted to two or three.
In the 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 The polarization sum of photons, pol ǫ γ (k) * λ ǫ γ (k) ρ = −g λρ , is performed, then we have 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 , by replacing the mass m e by m and the coupling constant e by g. So, (gβ) depending on the velocity β does not appear in our formulation. On this point, refer to Ref. [7,8]. Now we can obtain the decay rates of the para-and ortho-monopoliums as Γ(X 1 → 3γ) = 16(π 2 − 9) 9 where m is the monopole mass and is not the monopolium mass m X . These are the formulae valid at low energy, in the case of √ŝ = m X ≈ 2m. When √ŝ = m X ≫ 2m holds, and the binding energy is large, the following high energy limit formulae should be used: where ρ = ln(ŝ/m 2 ) = ln(m 2 X /m 2 ), and ζ(3) = ∞ n=1 n −3 =1.202 [22].
For the ortho-monopolium X 1 , it can decay to quark-antiquark, lepton, and weak boson pairs. The corresponding decay rates are where the electric charge of the quark and lepton is Q q,ℓ e. These formulae are derived for massless initial fermions. So, they have to be applied for the large binding energy case, or m ≪ m X and m W ≪ m X . It is noted that a large number eg/4π appears here, but it is not serious, since quarks, leptons and weak bosons don't carry magnetic charge.
The serious contribution to the decay rates comes from the emission of multiphotons, since the (electric or usual) photon can be converted to the magnetic photon and it couples to the monopole with a strong coupling g. Therefore, the total decay rates become It is difficult to estimate the multi-photon emission rate. However, the contribution of multi-photon emission can reduce the branching ratios Br(X 0 → 2γ) and Br(X 1 → 3γ) enormously, so we have to take them into account. One way is to consider the emitted photons, other than two or three hard photons, to be all soft, namely its energy E soft γ is less than the energy cut E cut γ . The energy resolution ∆E of photons at ATLAS is However, to analyze the data of two photon events, low energy photons are discarded by the energy cutoff E cut γ , which is much larger than the energy resolution [23,24]. We choose the energy cut as 35 GeV. The strategy is, by using the infrared technique, to sum up the emitted but undetectable photons with energy There is another cut of eliminating the photons emitted forward or backward, which is the cut on the pseudo-rapidity η = − ln tan(θ/2) (θ is the scattering angle), namely, It is rather technical, but we can convert the energy cut and the rapidity cut into nonvanishing photon masses λ, since the introduction of the photon mass prevents the infrared singularity which occurs at the zero photon energy, as well as the collinear (mass) singularity that occurs when the photon is emitted parallel to the emitting fermion (monopole in this case) 4 . The conversion is done by the following dictionary: Upon combining the two cuts, in the case of m X ≥ 750 GeV, the angle cut is stronger, so that the soft photons emitted at LHC satisfy The discussion so far has pertained to the emitted photons. The technique of summing infrared photons is as follows: The contribution of one soft photon emission to the decay width (or the cross section) has an infrared divergence, but the divergence is canceled by the other infrared divergence coming from the radiative corrections of the one soft photon exchange. As a result the decay width is multiplied by an infrared free factor Φ, called the eikonal factor. This eikonal factor can be easily summed to form an exponential factor, by including multi-photons. Virtual photons should also be soft, since the following approximation is taken in the numerator of the fermion propagator, where p is the fermion's (monopole's) momentum and k is the photon's momentum. This means Therefore, for the virtual soft photon corrections, we have 4 The lower bound of the massless photon propagator is identified with the lower bound of massive photon propagator, namely E cut , which gives Eq. (39). Similarly the denominator of the fermion propagator p · k is compared between the case of a massless photon with an angle cut and that of a massive photon without an angle cut. This gives a massless fermion, 1 − cos θ cut = 1 + (λ/|k|) 2 − 1, from which Eq. (40) is obtained.
Then, the cancelation between real emissions and virtual corrections leaves the soft photons in the following energy interval: If we assume that the monopole has four-momentum p 1 , the antimonopole has p 2 inside the monopolium, and t = (p 1 − p 2 ) 2 = 4p 2 rel = −4p 2 rel < 0, then We can rewrite Φ using the above dictionary as follows: where The last approximate formula, in Eq. (49), being valid in the high energy limit, is written as a reference, but is useful to estimate Φ roughly (see for example Ref. [25]). We will use Eq. (48), since the high energy limit |t| ≫ m 2 does not necessarily work in our case. This double logarithmic approximation in the high energy limit was used to estimate the total cross section σ(e + e − → multi-photons). Please refer to Ref. [22] and references theirein. Using this eikonal factor, the total decay widths can be written as, since we have to add an even number of soft photons to the 2γ decay or 3γ decay, the cosh Φ factor appears. Now we understand that the branching ratio is roughly given by, Derivation of this suppression factor in the branching ratios is very naive, but is indispensable for the strong coupling dynamics of the monopolium.

A model of a monopole and a monopolium based on the Zwanziger model and lattice gauge theory
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 [11,12,13]. The space-time is considered to be a square lattice n = (n 0 , n 1 , n 2 , n 3 ) (n µ is an integer) with a lattice constant a. The link variables U The minimum Wilson loop is given for the boundary curve C nµν ≡ ∂P nµν of the minimum plaquette 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 where ǫ 1234 = 1.
In the action, we assume that the magnetic coupling g is strong and , 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 [11,12,13]. 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. This point is very important. Define the potential between a heavy monopole and its antimonopole separated by a distance r to be V M M (r). Then, the potential has a linear term in r, if the monopole and antimonopole are separated in the η-direction; This clarifies the meaning of the special direction η that appears in the Zwanziger formulation. It gives the direction of the Dirac string starting from the monopole. Therefore, the monopole and antimonopole are connected by the string, starting from the monopole and ending at the antimonopole, 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 antiquark pair is well approximated by the linear plus Coulomb potential [26,27,28]; this potential is good 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 first-order phase transition between the confinement phase and the perturbative phase [19], 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 asymptotic free theory, when we approach to the short distance region. As is shown by 't Hooft and Polyakov [29,30,31], the U(1) monopole was given as a classical solution of SU(2) gauge theory, which is broken to U(1) with a triplet Higgs field φ a (x) (a = 1 − 3). Therefore, if we go inside the monopole, the non-Abelian gauge theory may appear. If this happens 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 [37,38]. Then, the equations of motion become the first order: from which the solutions read 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 [39]: 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 can treat the charged sphere as follows: If the sphere of radius R is uniformly charged, it can be viewed as a point charge Q(∞) from a long distance, but inside the sphere (r < R) the charge is reduced to Q(r) = Q(∞)(r/R) 3 . If we can apply this kind of treatment for the magnetic charge, the Coulomb potential can be relaxed at short distances inside the monopole. It is not clear, however, whether the magnetic charge can be distributed uniformly inside the sphere-like monopole.
Here, we will consider another possibility in which the magnetic charge is distributed uniformly on the surface of the sphere, fortunately such a solution 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 antimonopole vanishes when one goes inside either respective sphere (r < R). So, one possibility is to use the following monopole-antimonopole potential V which is defined by three parameters g, a and R. The idea of relaxing the Coulomb potential near the monopole is used by Epele et al. [7,8]. They follow Schiff and Goebel and relax the potential through the existence of Dirac strings [40,41]. Their choice of the potential is without the linear term.
Another possibility is to consider the linear potential to be the dominant contribution and the Coulomb potential just playing the role of lowering the potential. Then, the corresponding potential becomes = const = − g 2 4πR + ln(2g 2 ) a 2 R (for r < R).
The difference between the two potentials is = 0 (for r < R).
We can also include the spin-dependent hyperfine interaction. Then, the Hamil-tonianĤ, which we will use to study the energy and wave function of the lowest s-wave states ( 1 S 0 and 3 S 1 ) of the monopolium, iŝ where the nonrelativistic HamiltonianĤ NR and the hyperfine interactionĤ HF are defined byĤ In the case of the Coulomb potential, g 2 δ (3) (r) appears from ∇ 2 V (r) and the wave function at the origin contributes to the hyperfine interaction. If we use the potential V where the string tension is denoted by κ = ln(2g 2 )/a 2 , and (σ 1 · σ 2 ) = −1 for X 0 and +3 for X 1 . We don't include the relativistic corrections, such as the other Breit terms and the annihilation effect for the hyperfine interaction, since we intend to solve the problem semiclassically. The mass m X of the monopole-antimonopole bound state is obtained as the eigenvalue of the total HamiltonianĤ.

Mass and wave function of monopolium
In this section we use the potential V M M (r), and obtain the mass and wave function of the monopolium semiclassically. Schrödinger's equation with the linear potential can be solved in terms of the Airy function, but the solution is essentially equal to the semiclassical treatment [42], so we use the semiclassical method.
For a given nonrelativistic energyẼ, that is an eigenvalue E NR of the nonrelativistic HamiltonianĤ NR , classical momentum p r and |p r | are given by p r = m(Ẽ − V (r)), and |p r | = m(V (r) −Ẽ).
The turning point r = r * is the point where p r * = 0. If we define the wave function by ψ(r) = χ(r)/r, then χ(r) satisfies the onedimensional Schrödinger equation, and the corresponding semiclassical (WKB) wave function reads and the Bohr-Sommerfeld quantization condition is given by where n (= 1, 2, · · ·) denotes the radial excitation number, or the principal quantum number.
If we apply this for the potential V (b) (r), then the wave function is analytically given by, where we use the notation of the string tension as κ = ln(2g 2 )/a 2 . The Bohr-Sommerfeld quantization condition becomes a third order algebraic equation; for This is the case of one real root and two mutually complex conjugate roots, and so the real root is given by the Cardano formula as where p = 3m 2 κR, and q = − 3m 2 κπ(n − 1/2).

Comparison with experiments, LHC and others
In this section, on the basis of the formulae obtained so far, we will compare the numerical estimation of the mass and the decay width of monopoliums ( 1 S 0 and 3 S 1 ) , and the production rate with experiments such as LHC and others. There are typically two kinds of Solutions, 1 and 2: Solution 1 has a large monopole mass about half of the monopolium mass, namely 2m ≈ m X . In this case the nonrelativistic energyẼ is one order of magnitude smaller than the monopole mass, and the hyperfine interaction is also negligibly small. This corresponds to (Case 2) in the previous section. See Table 1 in which a few examples of Solution 1 are listed, where J = 0 and J = 1 monopoliums are degenerate with a mass of a several TeV. Estimation of the total width Γ tot (X) ≈ Γ(X → 2γ or 3γ) cosh Φ is performed using formulae (48) and (53) with Γ(X 0 → 2γ) (28) and Γ(X 1 → 3γ) (29). However, the very large total widths and the very small lifetimes suggest that our approximation is too naive, and a more reasonable estimation should be developed. We can say here that the estimation of the total widths and branching ratios are important for the strong coupling monopole dynamics. At this stage, the states considered are difficult to observe as a resonance peak at the LHC due to their high masses, but multi-photon emission decays may be observable at future 100 TeV colliders. Other than collider experiments, they could possibly be observed in astrophysical observations, through an explosion or burst of γ rays from the monopoliums. The number distribution function f γ (n) for photons can be very roughly estimated using the eikonal factor as f γ (n) = Φ n /n! cosh Φ (n = even for X 0 , odd for X 1 ).
However, in order to give a definite statement for the number distribution function of photons, a more reasonable estimation of Φ is necessary.
The second case, Solution 2, is obtained when considering very light monopole masses, for example 870 MeV. It is found that the large energyẼ can lead to a In the estimation of the production cross section for X 0 we have used the formula of Csaki et al. [45], To attempt to replicate the 750 GeV diphoton excess we make an appropriate choice of parameters to obtain an X 0 mass of m X 0 = 750 GeV, and an X 1 mass of m X 1 = 1410 GeV. The production cross section in the 2γ channel at 13 TeV LHC is, however, at most 0.68 fb and is one order of magnitude smaller than the recent announcements. In this study the high energy limit of the decay rates Γ(X 0 → 2γ) (30), Γ(X 1 → 3γ) (31) are used, and the total widths are found to be 5 MeV, 6 MeV respectively by using Eq. (56). It is also found that the hyperfine interaction is large because of the small monopole mass, such that the X 0 and X 1 state masses are strongly separated. For the estimation of the cross section for X 1 we have used the decay width Γ(X 1 → qq), a branching ratio of X 1 to 3γ, Br(X 1 → 3γ) ∼ 1, and the parton luminosity which is derived from graphs by W.J. Stirling (MSTW2008NLO PDF) [44], As for the production cross section of X 1 , the usual fusion mechanism of quark and antiquark inside a proton is used. The value of the cross section 0.0025 ab is rather small. The parton luminosity 1 s dL qq (τ ) dτ τ = m 2 X s , however, increases when s becomes larger, so that the production cross section can be several times larger for 33 TeV. The integrated luminosity of the order of hundred (ab) −1 is necessary at 33 TeV LHC or the future 100 TeV colliders in order to detect 1.4 TeV X 1 .

Conclusion
Being inspired by the announcement of the 750 GeV excess in the two photon channel at LHC, we have considered seriously a manageable model of monopolium, without missing the essence of it. As a result, we have a number of small but interesting findings which are summarized as follows: (1) Zwanziger's electric and magnetic dual formulation works properly. Especially, the meaning of a special direction (n µ in Zwanziger and η µ in this paper) becomes a little more manifest, since in the lattice formulation, the linear potential appears only in this direction. This supports that η µ is in the direction of the string going from the monopole and terminating at the antimonopole. 6 (2) For the strong coupling, the lattice version of Zwanziger gives a linear potential.
(3) As a finite sized monopole, the BPS equation for the 't Hooft-Polyakov monopole has another solution in which the U(1) magnetic charge is distributed uniformly on the surface of a sphere. So, the potential becomes flat inside the monopole.
The hyperfine interaction is modified depending on the potential which includes the non-perturbative effects.
(5) The mass and wave function of the monopolium (n 1 S 0 and n 3 S 1 ) can be given in a tractable manner, by applying the semiclassical method to the potential. (6) The branching ratio can be estimated very naively, by using the technique of infrared divergences, leading to Br(n 1 S 0 → 2γ) ≈ Br(n 3 S 1 → 3γ) ≈ 1/ cosh Φ, where Φ is the calculable eikonal factor that depends on the detection cut for the photon energy and angle. However, the total width calculated in this way is very big in Table 1 (not big in Table 2), so that we have to find a more reasonable way to estimate it before a definite statement can be given.
Numerical estimation based on our model of monopolium shows that there are typically two kinds of solutions: (7) One kind of solution (Solution 1) is the large monopole mass case in which the binding energy and the hyperfine splitting are very small. So, X 0 and X 1 are degenerate. A few examples considered give masses of monopoliums to be 5, 6, and 7 TeV. The signals of the multi-photon explosion or burst may be detected by 100 TeV colliders or by astronomical observations. (8) The other kind of solution (Solution 2) gives an 870 MeV monopole, and the monopole and its antimonopole form a 750 GeV J = 0 state and a 1.4 TeV J = 1 state. The solution can be understood similarly as in QCD, where the so-called small current mass of the monopole is heavily dressed with non-perturbative coats and makes the monopoliums heavy. Unfortunately this solution can't explain the 750 GeV excess at LHC, because the predicted diphoton production rate is one order of magnitude smaller than that of the observed excess. The 1.4 TeV J = 1 state may provide a target at LHC or future 100 TeV colliders. This paper has studied how to take into account strong coupling effects in monopole dynamics, especially in the dynamics of monopolium. We provide a method, in which the bound state is formed non-perturbatively using a strong coupling expansion in lattice gauge theory, while the decay rate is estimated perturbatively, since the final number of emitted photons is fixed. However, the strong magnetic coupling may give large branching ratios for the multi-photon emission processes, even if they are not intended to be detected. Overcoming this point is important and is an issue requiring future work.