Looking for magnetic monopoles at LHC
Luis N. Epele, Huner Fanchiotti, Carlos A. Garcia Canal, Vasiliki A. Mitsou, Vicente Vento
aa r X i v : . [ h e p - ph ] A p r FTUV-11-2702IFIC/11-14
Looking for magnetic monopoles at LHC
Luis N. Epele a , Huner Fanchiotti a , Carlos A. Garc´ıa Canal a , Vasiliki A. Mitsou b and Vicente Vento b,c (a) Laboratorio de F´ısica Te´orica, Departamento de F´ısica, IFLPFacultad de Ciencias Exactas, Universidad Nacional de La PlataC.C. 67, 1900 La Plata, Argentina. ( E-mail: epele@fisica.unlp.edu.ar, huner@fisica.unlp.edu.ar, garcia@fisica.unlp.edu.ar )(b)
Instituto de F´ısica CorpuscularUniversidad de Valencia and CSICApartado de Correos 22085, E-46071 Valencia, Spain. ( E-mail: vasiliki.mitsou@ific.uv.es )(c)
Departamento de F´ısica Te´oricaUniversidad de ValenciaE-46100 Burjassot (Valencia), Spain. ( E-mail: [email protected] ) Abstract
Magnetic monopoles have been a subject of interest since Dirac established therelation between the existence of monopoles and charge quantization. The intense ex-perimental search carried thus far has not met with success. We study the observabilityof monopoles at the Large Hadron Collider in the γ γ channel and show that LHC isan ideal machine to discover monopoles with masses below 1 TeV at present runningenergies and with less than 1 fb − of integrated luminosity.Pacs: 14.80.Hv, 95.30.Cq, 98.70.-f, 98.80.-kKeywords: Quantum, electrodynamics, duality, monopoles, photon. Introduction
The theoretical justification for the existence of classical magnetic poles, hereaftercalled monopoles, is that they add symmetry to Maxwell’s equations and explain chargequantization [1, 2]. Dirac showed that the mere existence of a monopole in the universecould offer an explanation of the discrete nature of the electric charge. His analysisleads to the Dirac Quantization Condition (DQC), e g = N , N = 1,2,... , (1)where e is the electron charge, g the monopole magnetic charge and we use naturalunits ¯ h = c = 1. In Dirac’s formulation, monopoles are assumed to exist as point-likeparticles and quantum mechanical consistency conditions lead to Eq.(1), establishingthe value of their magnetic charge. Their mass, m , is a parameter of the theory.Monopoles and their experimental detection have been a subject of much studysince many believe in Dirac’s statement [1], “...one would be surprised if Nature had made no use of it [the monopole].” At present, despite intense experimental search, there is no evidence of their existence[3, 4]. LHC is opening a new frontier in the search for monopoles and, as we showhere, should find them at present energies, with less than 1 fb − integrated luminosity,if their mass is below 1 TeV.Although monopoles symmetrize Maxwell’s equations in form there is a numericalasymmetry arising from the DQC, namely that the basic magnetic charge is muchlarger than the smallest electric charge. This led Dirac himself in his 1931 paper [1] tostate, “... the attractive force between two one quantum poles of opposite sign is (1372) ≈ times that between the electron and the proton. This very large force mayperhaps account for why the monopoles have never been separated.” This statement by Dirac is fundamental in our investigation. If instead of monopoleswe deal with monopole-antimonopole pairs, we expect that due to this very stronginteraction many of them annihilate into photons inside the detector. We study theproduction of a monopole-antimonopole pairs via photon fusion at LHC and theirsubsequent annihilation giving rise to two extremely energetic photons.Another proposal, influenced also by the strength of the interaction, assumes thatthe produced pair of monopole-antimonopole, before decaying, forms a monopoliumbound state [5, 6]. This bound state will also desintegrate producing two photons, aprocess which is being analyzed at present [7] in connection with LHC potentialities.1n the next section we describe the dynamics of monopoles and review the produc-tion of monopole-antimonopole pairs. Section 3 discusses the annihilation of the virtualmonopole-antimonopole pair into photons. In section 4 we describe how to incorporatethe elementary process into p − p scattering. In section 5 we present our results inthe context of the present running features of LHC and in the last section we drawconclusions of our studies. The theory of monopole interactions was initially formulated by Dirac [8] and later ondeveloped in two different approaches by Schwinger [9] and Zwanziger [10]. The formal-ism of Schwinger can be cast in functional form as a field theory for monopole-electroninteraction which is dual to Quantum Electrodynamics (QED) [11]. The monopolesare considered fermions and behave as electrons in QED with a large coupling constantas a result of the DQC. In this formulation the conventional photon field is Dirac stringdependent. Due to the large coupling constant and the string dependence, perturbativetreatments `a la Feynman are in principle not well defined. However, non perturbativehigh energy treatments, like the eikonal approximation, have rendered well definedelectron-monopole cross sections [11, 12].For the case of monopole production at energies higher than their mass the aboveprocedure is not applicable, and being the treatment non perturbative, there is nouniversally accepted prediction from field theory. However, the study of the classicalinteraction of a monopole passing by a fixed electron leads to an interaction for themonopole which is associated with the electric field felt. If one compares this interactionwith the that of a positron passing by an electron, one realizes that the differencebetween QED and dual QED is simply to change the electric charge by the magneticcharge times the velocity, i.e. e → βg. (2) gβgβ gβgβ Figure 1: Elementary processes of monopole-antimonopole production via photonfusion. 2hus a monopole behaves as a passing particle with electric charge βg [13].This idea has been used to define an effective field theory to lowest order [14] whichhas been applied to Drell-Yan like monopole, or monopole-antimonopole productionmechanisms [14, 15] and to monopole-antimonopole production by photon fusion [16,17] . In Fig. 1 we show the corresponding diagrams that have been calculated for thelatter process.We note that this approximation can be shown to be almost equivalent to the lowenergy effective theory of Ginzburg and Schiller [18, 19]. The theory of Ginzburg andSchiller was derived from the standard electroweak theory in the one loop approxi-mation leading to an effective coupling of the order of g eff ∼ εm g , where m is themonopole mass and ε a kinematical energy scale of the process which is below themonopole production threshold, thus rendering the effective theory perturbative. Ina photon fusion diagram the dynamical scale is √ E − m , where E is the center ofmass energy, thus εm ∼ √ E − m m ∼ Eβ m , (3)where β = q − m E is the monopole velocity. Therefore if E ∼ m , i.e. the kineticenergy is small, and both schemes coincide [6] .The aim here is to study possible signals of magnectic monopoles at LHC. Ac-cording to previous studies [17], the most promising mechanism is photon fusion. Theelementary diagrams contributing to pair production are those in Fig. 1, where theexplicit couplings have been shown. e - - e + m - m Ω Σ H Ω L Figure 2: Elementary photon-fusion cross section for electron-positron (solid) and thatof the monopole-antimonopole (dashed) as a function of ω , to show the importance ofthe β effect.The photon-fusion elementary cross section is obtained from the well-known QEDelectron-positron pair creation cross section [20], simply changing the coupling constant3 e → gβ ) and the eletron mass by the monopole mass m e → m , leading to σ ( γ γ → mm ) = π g (1 − β ) β m − β β log β − β ! − (2 − β ) ! , (4)where E the center of mass energy and β , a function of E , is the monopole velocity.In Fig. 2 we show the ω = E/ m dependence of Eq.(4). The solid curve correspondsto the electron-positron case, the dashed one to the monopole case which contains the β factor. One should notice the large effect associated with this factor in the vicinityof the threshold.LHC detectors, apart from the MoEDAL experiment [21], have not been designedspecifically to see monopoles and therefore even those which do not annihilate insidethe detectors will be difficult to detect. However, all detectors are well designed tomeasure low cross sections of photons, since the two photon decay is one of the favorablechannels to detect a low mass Higgs, thus we would like to benefit from this capabilityand proceed to investigate how monopole-antimonopole annihilation into photons isdescribed in our scheme. γ γ It is natural to think that the enormous strength and long range of the monopoleinteraction leads to the annihilation of the pair into photons very close to the productionpoint. Thus one should look for monopoles through their annihilation into highlyenergetic photons, a channel for which LHC detectors have been optimized.In order to calculate the annihilation into photons we assume that our effective the-ory, a technically convenient modification of Ginzburg and Schiller’s, agrees with QEDat one loop order, and therefore we apply light-by-light scattering with the appropriatemodifications as shown in Fig. 3. An interesting feature of the calculation is that the gβgβ gβgβ gβgβ gβgβ
Figure 3: Elementary processes for monopole-antimonopole production and annihila-tion into photons. 4dditional magnetic coupling will increase the cross section dramatically and thereforethis measurement should lead to a strong restriction on the monopole mass. However,as we have seen, in m − m production, the large magnetic coupling is always multipliedby the small electric one, leading to effective couplings e g , and the same will happen indetection. Thus, the effective coupling of the process is e g , and therefore has strengthssimilar to the strong interaction, not more. The extreme sensitivity of LHC detectorsto photons, as a consequence of the hope to detect the Higgs in the γ γ channel makesthem ideal for the purpose of detecting monopoles, even if the monopole mass is large,as will be shown. - - Ω X H Ω L Figure 4: The light-by-light scattering quantity X ( ω ), which determines the crosssection up to factors related to the coupling constant. The solid curve corresponds toforward scattering, the dashed one to right-angle scattering, both as a function of thedimensionless quantity ω .Light-by-light scattering was studied by Karplus and Neuman [22, 23], who repro-duced all previous low energy results by Euler [24] and high energy results by Achieser[25], and was later revisited and corrected by Csonka and Koelbig [26].In the case of monopoles, with the appropriate changes, the expression for the crosssection becomes, σ γγ ( θ, ω ) = ( gβ ) π m X ( θ, ω ) . (5)where X ( ω ) is < | M | >ω in the notation of ref. [23]. In Fig. 4 we reproduce these resultsfor forward scattering and right-angle scattering. In Fig. 5 (left) we show the ratioof the forward scattering to the right-angle cross sections and therefore show how thecross section, which is basically isotropic close to threshold, becomes anisotropic as theenergy increases. In the same figure (right) we plot X ( ω ) for several values of ω asa function of angle. We note that the forward cross section is larger and the right-5ngle one smaller, than that for any other value of the scattering angle. As the energyincreases the drop in the X function from θ = 0 to π/ Ω R a ti o @ Ω D Π Π Θ X H Ω L Figure 5: Left: ratio of the light-by-light scattering quantity X ( ω ) for forward to right-angle scattering as a function of ω . Right: angular dependence of the cross section forthree values of ω (1 (solid) , 5 (dashed), 10 (dotted)) .It is clear from these figures that close to the threshold the cross section is quiteisotropic and away from threshold the forward cross section, which is very difficult tomeasure, is much larger than the right-angle one. Since the detectors cannot detectall of the photons coming out, we take the right-angle cross section as a conservativeindication of the magnitudes to be expected. ´ Β - - Ω X H Ω L Figure 6: The effect of the β factor as a function of ω . The solid curve corresponds toelectron scattering, while the dashed curve corresponds to the monopole case.In Fig. 6 we show the effect of the β factor which diminishes greatly the crosssection close to the monopole-antimonopole threshold.6 P γγγγ P ( X ) P ( X ) Figure 7: Processes contributing to the γγ cross section. The blob contains the threecases described in the text. LHC is a proton-proton collider, therefore, in order to describe the production anddesintegration of the monopole-antimonopole pair, we have to study the following pro-ceses above the monopole threshold ( β > → E ≥ m ), p + p → p ( X ) + p ( X ) + γ + γ, (6)shown globally in Fig.7, where p represents the proton, X an unknown final state andwe will assume that the blob is due exclusively to monopoles. This diagram summarizesthe three possible processes:i) inelastic p + p → X + X + γ + γ → X + X + m + m → X + X + m + m + γ + γ ii) semi-elastic p + p → p + X + γ + γ → p + X + m + m → p + X + m + m + γ + γ iii) elastic p + p → p + p + γ + γ → p + p + m + m →→ X + X + m + m + γ + γ .In the inelastic scattering, both intermediate photons are radiated from partons(quarks or antiquarks) in the colliding protons.In the semi-elastic scattering one intermediate photon is radiated by a quark (orantiquark), as in the inelastic process, while the second photon is radiated from theother proton, coupling to the total proton charge and leaving a final state proton intact.In the elastic scattering both intermediate photons are radiated from the interactingprotons leaving both protons intact in the final state.7n the blob we incorporate the elementary subprocess shown in Fig. 3 and describedby Eq. (5).We calculate γγ fusion for monopolium-antimonopolium production following theformalism of Drees et al. [27].In the inelastic scattering, p + p → X + X + ( γγ ) → X + X + m + m + γ + γ ,to approximate the quark distribution within the proton we use the Cteq6–1L partondistribution functions [28] and choose Q = ˆ s/ p + p → p + X +( γγ ) → p + X + m + m + γ + γ , the photonspectrum associated with the interacting proton must be altered from the equivalent–photon approximation for quarks to account for the proton structure. To accommodatethe proton structure we use the modified equivalent–photon approximation of [27].The total cross section is obtaned as a sum of the three processes. The explicitexpressions for the different contributions can be found in [17].In order to obtain the differential photon-photon cross section from the above for-malism we develop a procedure which we exemplify with the elastic scattering case. Inthat case the pp cross section is given by [27], σ pp ( s ) = Z m /s dz Z m /sz dz f ( z ) f ( z ) σ γγ ( z z s ) , (7)where √ s is the center of mass energy of the pp system.We perform the following change of variables v = z z , w = z , which leads to σ pp ( s ) = Z m /s dv Z v dww f ( vw ) f ( w ) σ γγ ( vs ) . (8)Note that to fix the center of mass energy of the photons is equivalent to fix v . Forfixed v we have, dσ pp dv ( s ) = Z v dww f ( vw ) f ( w ) σ γγ ( vs ) , (9)which can be rewritten in terms of E γ , the center of mass energy of the photons, andthe elementary photon-photon cross section as, dσ pp dE ( E γ ) = 2 E γ s σ γγ ( s γγ ) Z s γγ /s dww f ( s γγ w ) f ( w ) . (10)This procedure can be generalized easily to the semielastic and inelastic cases, wherethe appropriate change of variables are 8
00 500 600 700 800 900 1000110100100010 m @ GeV D Σ @ f b D Figure 8: The total monopole-antimonopole production cross section (solid line) from γ γ fusion in fb as a function of monopole mass. The different components of the crosssection are shown: elastic (dotted line), semielastic (small dashed line) and inelastic(dashed line). v = z z x , w = z x , u = x and v = z z x x , w = z x x , u = x x , t = x , respectively. In Fig. 8 we show the total cross section for monopole-antimonopole production andthe contribution from each of the individual processes described above as a function ofmonopole mass. In order to set the mass axes we have taken into account the previouslower mass limit for the monopole of ref. [15] which was set at 360 GeV. The crosssection is of O (fb) thus we limit the high mass values to those potentially observableat present, or in the near future, by the LHC with an integrated luminosity of 1 fb − .Our result differs from that of ref. [17] around threshold.We see in the previous figure that monopoles of mass around 500 GeV should beproduced abundantly, while those of mass around 1000 GeV have a cross section whichmakes them difficult to detect in the near future, at least directly.Since, as already mentioned, the LHC detectors have not been tuned to detectmonopoles but are excellent detectors for γ ’s let us discuss next our annihilation crosssection. In Fig. 9 we show the differential cross section for forward (solid) and right-angle (dotted) scattering given in fb/GeV as a function of the invariant mass of the γ γ system. We have assumed a monopole of mass 750 GeV, chosen because the cross9ections turned out to be close to the expected magnitude of the Higgs to γ γ crosssection above background. The cross sections are wide, almost gaussian, structuresrising softly just above threshold (1500 GeV). LHC detectors are blind for forwardscattering and have black spots due to construction features in the non-forward regionswhich do not allow for a full detection of photons. Therefore in order to obtain aneducated estimation of the observable cross section we take the right-angle cross sectionand multiply it by 4 π . This differential cross section is the smallest possible. However,away from threshold, it corresponds quite well to a realistic estimate, since, as we haveseen, the elementary differential cross section drops fast with angle and moreover oneshould consider an efficiency factor for the various detectors. Γ @ GeV D d Σ B f b G e V F Figure 9: Forward (solid) and the right-angle cross sections for a monopole of mass750 GeV.In Fig. 10 we plot the right-angle differental cross section for the same monopolemass and show the different contributions to the differential cross section. The structureof the cross section is a wide structure, rising softly after threshold, 1500 GeV, andextending for almost 2000 GeV. The structure is centered about 2300 GeV. Clearlythe soft rise of the differential cross section is a signature of the two particle threshold,recall the β factor. The width of the stucture is associated to the mathematical form ofthe box diagram, as can be seen for both electron-positron annihilation and monopole-antimonopole, from the structue of the elementary cross section (Fig. 6).In Fig. 11 we compare our total γγ cross section with Higgs process obtained fromref. [31]. We have extrapolated their background to our energies using an inversepolynomial fit to their data and their exponential fit. Both procedures give a negligiblebackground for the signals obtained with monopoles masses up to 1 TeV and evenhigher. In the left figure we traslate the monopole-antimonopole threshold to the originin order to compare the two signals. In the right figure we show the actual energy scalesin a LogLog plot. The Higgs signal above background has been multiplied by 50 tomake it visible. The figures correspond to a monopole of mass 750 GeV. It is clear10
500 2000 2500 3000 3500 40000.00.51.01.52.02.53.0 E Γ @ GeV D d Σ B f b G e V F Figure 10: The right-angle scattering cross sections for a monopole mass of 750 GeV.The smallest is the inelastic cross section (dotted), next comes the elastic (dashed) andthe biggest is the semielastic (long-dashed). The total cross section, the sum of thethree, is represented by the solid curve.from the curves that monopole-antimonopole annihilation should appear as a soft riseof the cross section above the background over a large energy interval. H ® ΓΓ m - m ® ΓΓ
200 400 600 800 1000 1200 14000100200300400 E Γ @ GeV D d Σ B f b G e V F H ® Γ Γ m- m ® ΓΓ
100 200 500 1000 2000 50000.11101001000 E Γ @ GeV D d Σ B f b G e V F Figure 11: Comparison of the γγ monopole-antimonopole annihilation cross section fora monopole of mass 750 GeV with the Higgs γγ decay. The Higgs cross section abovethe background has been multiplied by 50. The monopole-antimonopole threshold hasbeen set at the origin (100 GeV) (left). The right figure represents the same crosssection drawn in a Log Log plot keeping the thresholds.Finally we study the mass dependence of the differential cross section. Fig. 12 is aLogLog plot of the cross section for three masses, 500 ,
750 and 1000 GeV. In order tohave a better comparison we have translated the thresholds to the origin. It is clearthat the magnitude of the cross section and its extent falls rapidly with mass. For alow mass monopole, the magnitude of the differential cross section is of the order of ∼
10 fb/GeV and for a heavy monopole is of the order of ∼ M = 500 GeVM = 750 GeVM= 1000 GeV1000 2000 300015000.0010.1101000 E Γ @ GeV D d Σ B f b G e V F Figure 12: The cross section for monopole-antimonopole annihilation into γγ for threemasses: 500 GeV (solid), 750 GeV dashed 1000 GeV (dotted). The existence of Dirac monopoles would modify our understanding of QED. The DQC,implying a huge magnetic coupling constant, complicates matters from the theoreticalpoint of view. Non perturbative methods are required. We have avoided the problem byusing an effective theory valid at one loop order. The physical coupling turns out to be e g , equivalent in magnitude to the pionic couplings, and the expected high monopolemasses provide convenient cut-offs which make the theory sensible. In the context ofthis theoretical scheme we are able to calculate monopole(antimonopole) productionand monopole-antimonopole annihilation. Centering our description to LHC, a protonmachine, we have described scenarios for the production of monopole-antimonopolepairs via photon fusion. We have described their annihilation into photons via theconventional box diagram, analogous to that of light-by-light scattering in conventionalQED. The only difference in our approach is the introduction of a threshold factor inthe form of a velocity. We have analyzed in detail the elementary process, light-by-lightscattering at monopole energies, and have subsequently incorporated the description ofthe photon-photon elementary scattering into the initial proton-proton collision. Onemain assumption is that no other particle contributes to the box diagram in the regionof interest. We have chosen in the present calculation a mass range which is aboveconventional particle masses and below the supersymmetric particle ranges. One couldthink of interference with elementary particle decays, however the latter would have aresonant structure which our box diagram does not have. Thus any interference wouldbe avoided by their different geometric structure.12ith all these preventions we have shown that monopoles up to 1 TeV in massshould be detected at LHC in the γ γ channel at present running energies and withluminosities in reach within the next few months. The signature is a soft rise of thecross section which should stay for a long energy range, because no resonant peakstructure is to be expected.
Acknowledgement
We thank the authors of JaxoDraw for making drawing diagrams an easy task [32].LNE, HF and CAGC were partially supported by CONICET and ANPCyT Argentina.VAM acknowledges support by the Spanish Ministry of Science and Innovation (MICINN)under the project FPA2009-13234-C04-01, by the Spanish Agency of International Co-operation for Development under the PCI projects A/023372/09 and A/030322/10and by the CERN Corresponding Associate Programme. VV has been supported byHadronPhysics2, by MICINN (Spain) grants FPA2008-5004-E, FPA2010-21750-C02-01, AIC10-D-000598 and by GVPrometeo2009/129.
References [1] P. A. M. Dirac, Proc. Roy. Soc. Lond. A (1931) 60.[2] J. D. Jackson, Classical Electrodynamics, de Gruyter, N.Y. (1982).[3] K. A. Milton, Rept. Prog. Phys. (2006) 1637 [arXiv:hep-ex/0602040].[4] W. M. Yao et al. [Particle Data Group], J. Phys. G (2006) 1.[5] L. N. Epele, H. Fanchiotti, C. A. Garcia Canal and V. Vento, Eur. Phys. J. C (2008) 87 [arXiv:hep-ph/0701133].[6] L. N. Epele, H. Fanchiotti, C. A. G. Canal and V. Vento, Eur. Phys. J. C (2009) 587 [arXiv:0809.0272 [hep-ph]].[7] L. N. Epele, H. Fanchiotti, C. A. G. Canal, V. Mitsou and V. Vento, (work inpreparation).[8] P. A. M. Dirac, Phys. Rev. (1948) 817.[9] J. S. Schwinger, Phys. Rev. (1966) 1087.[10] D. Zwanziger, Phys. Rev. D (1971) 880.[11] L. P. Gamberg and K. A. Milton, Phys. Rev. D (2000) 075013[arXiv:hep-ph/9910526]. 1312] L. F. Urrutia, Phys. Rev. D (1978) 3031.[13] M. J. Mulhearn, “A Direct Search for Dirac Magnetic Monopoles,” Ph.D. Thesis,Massachusetts Institue of Technology 2004, FERMILAB-THESIS-2004-51.[14] G. R. Kalbfleisch, K. A. Milton, M. G. Strauss, L. P. Gamberg, E. H. Smith andW. Luo, Phys. Rev. Lett. (2000) 5292 [arXiv:hep-ex/0005005].[15] A. Abulencia et al. [CDF Collaboration], Phys. Rev. Lett. (2006) 011802[arXiv:hep-ex/0508051].[16] Yu. Kurochkin, I. Satsunkevich, D. Shoukavy, N. Rusakovich and Yu. Kulchitsky,Mod. Phys. Lett. A (2006) 2873.[17] T. Dougall and S. D. Wick, Eur. Phys. J. A (2009) 213 [arXiv:0706.1042 [hep-ph]].[18] I. F. Ginzburg and A. Schiller, Phys. Rev. D (1998) 6599[arXiv:hep-ph/9802310].[19] I. F. Ginzburg and A. Schiller, Phys. Rev. D (1999) 075016[arXiv:hep-ph/9903314].[20] C. Itzykson and J. B. Zuber, New York, USA: McGraw-Hill (1980) 705P.(International Series In Pure and Applied Physics) [21] J. L. Pinfold, AIP Conf. Proc. (2010) 234; J. L. Pinfold [MOEDAL Collab-oration], CERN Cour. (2010) 19.[22] R. Karplus and M. Neuman, Phys. Rev. (1950) 380.[23] R. Karplus and M. Neuman, Phys. Rev. (1951) 776.[24] H. Euler Ann. Phys. (1936) 398.[25] A. Achieser, Physik Z. Sowjetunion (1937) 263.[26] P. L. Csonka and K. S. Koelbig, Phys. Rev. D (1974) 251.[27] M. Drees, R. M. Godbole, M. Nowakowski and S. D. Rindani, Phys. Rev. D (1934) 729.1430] C. F. von Weizsacker, Z. Phys. (1934) 612.[31] G. Aad et al. [Atlas Collaboration], Phys. Rev. Lett. (2011) 121803arXiv:1012.4272 [hep-ex].[32] D. Binosi and L. Theussl, Comput. Phys. Commun.161