Simulating magnetic monopole-defect dynamics

We present simulations of one magnetic monopole interacting with multiple magnetic singularities. Three-dimensional plots of the energy density are constructed from explicit solutions to the Bogomolny equation obtained by Blair, Cherkis, and Durcan. Animations follow trajectories derived from collective coordinate mechanics on the multi-centered Taub-NUT monopole moduli space. We supplement our numerical results with a complete analytic treatment of the single-defect case.


Introduction and summary
A standard course in introductory physics can't help but suggest a false dichotomy, with particle mechanics -or more generally the mechanics of rigid bodies -on the one side, and field theory on the other side. One learns that charged particles create fields and fields apply forces to charged particles, but one is not presented with a complete description of a coupled particle-field system until later. Furthermore, the fields are typically singular in the vicinity of the particles sourcing them, so that the dichotomy is only truly resolved through quantum field theory. This is the situation, at least, for electromagnetism and its charged particles.
Solitons provide a fascinating alternative if one's goal is to see how particle dynamics can emerge from, and be completely embedded in, the framework of classical field theory in a nonsingular way. Static solitons are represented by self-supporting localized field configurations and occur in theories that admit a topological charge. 1 Solitons exist in JHEP04(2021)286 a variety of theories; classic examples include kinks (or domain walls) for theories in one spatial dimension, vortices for theories in two dimensions, and magnetic monopoles for theories in three dimensions. See [1] for a modern review.
The existence of solitons typically relies on nonlinearity in the field equations, and one cannot linearly superpose two one-soliton solutions to construct a two-soliton solution. Nevertheless, soliton solutions do come in smooth families with a number of parameters, or moduli, parameterizing the family. Focusing on the case of magnetic monopoles in Yang-Mills-Higgs theory [2,3], solutions representing a single monopole occur in a fourdimensional family. Three moduli represent the position of the monopole in R 3 while the fourth parameter is a circle coordinate, with momentum along this circle corresponding to electric charge [4]. This four-dimensional space of solutions is referred to as the onemonopole moduli space.
In this paper we work in the context of a special type of Yang-Mills-Higgs theory -namely, one which has (four-dimensional N = 2 extended) supersymmetry [5][6][7][8][9]. In supersymmetric Yang-Mills-Higgs theory, the net force vanishes between two stationary monopoles and hence there exist static n-monopole solutions coming in smooth 4ndimensional families for any positive integer n. 2 These monopole moduli spaces have been intensely studied by mathematicians and physicists alike since their invention in the late 1970's. They inherit a natural Riemannian metric induced from the energy functional of the parent Yang-Mills-Higgs theory that carries rather special geometric structures. Specifically, these moduli spaces are hyperkähler manifolds admitting a number of isometries; see [12] for details. In [12], Atiyah and Hitchin explotied these special structures to pin down the metric on the two-monopole moduli space, despite the fact that the full family of two-monopole solutions was not known at the time. 3 Furthermore, the collective coordinate paradigm of Manton [15] reaches its full brilliance in the context of multi-monopole moduli spaces. Collective coordinates are the ultimate example of the physicist's ball-rolling-on-a-hill. In this analogy, the ball is the multi-monopole configuration, the terrain of peaks and valleys is the infinite-dimensional space of field configurations, and monopole moduli space is a minimum-energy valley where the ball can roll without change in kinetic energy. In fact the analogy is rigorous. It can be proven that time-dependent solutions to the full field equations are well-approximated by allowing the moduli to become time-dependent -i.e. promoted to collective coordinates -so that they trace out specific trajectories in moduli space [16].

JHEP04(2021)286
The trajectories in moduli space are determined by a specific form of Newton's Laws, and hence the emergence of particle mechanics from field theory. In ordinary Yang-Mills-Higgs theory in the Bogomolny-Prasad-Sommerfield (BPS) limit, Newton's Laws imply that the trajectories are geodesics [15]. However, in the supersymmetric extension considered here, there is an additional force due to a secondary Higgs field that modifies the trajectories in a way compatible with the special structure of the moduli space [17,18]. 4 Motion on the two-monopole moduli space was studied in [12], where it was shown to predict rather beautiful and dramatic phenomena for the scattering of two monopoles in real space. For example, in a head-on collision, the two monopoles -represented by spherical blobs of energy when they are far apart -deform as they get close to each other. The monopoles' individual identities disappear momentarily as they overlap and form an axially symmetric ring of energy. The spherical blobs then re-emerge from the collision region traveling away from each other on a line rotated by 90 degress from the line of incidence. A simulation of the collision [19], constructed in the late 1980's on an IBM supercomputer, can still be found on YouTube. Very recently, an interactive applet has been constructed for the two monopole solution based on the new analytic results of [14]. See the final appendix of Reference [14] for discussion and links.
Fascinating n-monopole collisions with n > 2 have been studied for special initial conditions such that the configuration maintains some specific symmetry throughout the evolution; see e.g. [20,21]. The reason for this symmetry restriction is that the full moduli space geometry, required to simulate collisions with generic initial conditions, is not known for n > 2.
Our goal in this work is to explore and simulate multi-monopole interactions in a different limit -namely, when all but one of the monopoles are infinitely heavy and immobile, while the remaining one can move in accordance with the appropriate moduli space force law. The heavy monopoles can be placed at arbitrary fixed positions in three dimensional space and are modeled as magnetic singularities known as (supersymmetric) 't Hooft defects [22,23]. They can indeed be viewed as infinite-mass limits of ordinary monopoles in a precise sense described in [24,25]. This allows us to utilize a relatively recent set of analytic solutions obtained by Blair, Cherkis, and Durcan (BCD) [26][27][28][29], describing one ordinary mobile monopole in the presence of any number of fixed 't Hooft defects. We compute the energy density and construct three-dimensional plots, using Mathematica, for the smooth monopole in arbitrary position relative to the defects. Like the authors of reference [14], we plot several level sets of the energy density with varying opacity, so as to allow one to see inside the monopole configuration. See figure 1 below.
The moduli space metric for one monopole in the presence of k minimally charged singularities is also known and given by the k-centered Taub-NUT manifold [30]. This is a four-dimensional manifold constructed over an R 3 base. Over each point on the base is a fiber that is generically a circle. However the size of the circle varies, and the circle shrinks to a point over the k "nut" points of the base, whose positions are specified by k fixed three-JHEP04(2021)286 dimensional vectors. The moduli parameterizing the three-dimensional base correspond to the position of the smooth monopole in physical three-space, and the nut points represent the positions of the singularities. Momentum along the circle fiber represents electric charge carried by the monopole. There is also a potential energy function on the moduli space due to the secondary Higgs field in the supersymmetric model that we consider.
We numerically integrate Newton's Laws with Mathematica to determine the trajectory on moduli space for any given initial position and velocity. We then use Manton's collective coordinate ansatz and the explicit BCD solutions to create the resulting simulations of the monopole interacting with the defects in real space. There exist both bound orbits (generically non-repeating, but closed for special initial conditions) and unbound trajectories. When the smooth monopole passes near the defect we observe significant but transitory deformations of both the monopole and defect shape. Complete, momentary screening of the defect by the monopole is also observed when the defect carries the same charge as the monopole.
When only a single defect is present, the equations of motion on moduli space can be integrated analytically for generic initial conditions. We carry out this analysis as well, since the analytic results offer valuable insight into the more complex scenarios with multiple defects. This system is mathematically equivalent to several related systems that have been studied over the years, starting with work of Zwanziger [31][32][33]. In these references it was found that trajectories are conic sections. However, the plane of motion does not contain the defect when the monopole carries electric charge. We review and extend some of these results to include the period of bound orbits and a new and elementary approach to the classical scattering problem.
The structure of the paper is as follows. In section 2 we review the theoretical background for the emergence of monopole-defect solutions and dynamics in supersymmetric Yang-Mills-Higgs theory. Then in section 3 we apply these ideas to the explicit BCD solutions and illustrate them with three-dimensional energy density plots and simulations based on moduli-space dynamics. Section 4 contains our analysis of the single-defect case. We conclude in section 5 with a brief summary and description of future directions.
Several simulations are highlighted in subsections 3.4, 4.2, and 4.3. These simulations, and the Mathematica code required to create such simulations, are included as supplementary material. High resolution simulations can be created in a few hours to a couple days on current commercial laptops, depending on the number of defects included. Low resolution simulations can be made in a matter of minutes.

Supersymmetric Yang-Mills-Higgs with 't Hooft defects
We study a field theory on Minkowski space, R 1,3 with coordinates x µ = (t, r), consisting of a non-abelian gauge field A µ = (A 0 , A i ), i = 1, 2, 3, two adjoint-valued scalars, X, Y , and a pair of adjoint-valued Weyl fermions. Although the fermions are crucial for supersymmetry, they will not be utilized in the following and so we suppress their contribution to the action and Hamiltonian below. We work with the simplest nonabelian gauge group, G = SO(3).

JHEP04(2021)286
The three generators of the Lie algebra, g = so (3), are denoted T a , a = 1, 2, 3. We use anti-Hermitian generators satisfying [T a , T b ] = ab c T c and normalized such that Tr(T a T b ) = 1 2 δ ab . 5 Each field can be expanded in this basis; A i = A ia T a , X = X a T a , etc. The covariant derivative and non-abelian field strength tensor are The magnetic field is B i = 1 2 ijk F jk , and we work in mostly plus conventions such that the electric field is Gauge transformations act on the fields according to where g (x) ∈ SO(3). Taking g = exp( a T a ), the infinitesimal form of these transformations is Two field configurations related by a local gauge transformation, i.e. a gauge transformation with g → 1 (or → 0) as r ∈ R 3 is sent to infinity, are physically equivalent. In contrast, global gauge transformations -those that do not approach the identity at spatial infinity -generate symmetries that can be used to simplify asymptotic boundary conditions and generate conserved Noether charges. We probe this theory with magnetic defects known as 't Hooft defects [22], which can be thought of as magnetic duals to the Wilson lines of external electrically charged particles. 't Hooft defects are a type of disorder operator, in that they are defined not in terms of the local fields in the theory but rather in terms of singular boundary conditions on the fields. This construction was extended to the supersymmetric context in [23]. Such supersymmetric "line defects" have played a central role in many of the new theoretical developments of the past decade for supersymmetric Yang-Mills-Higgs theory, beginning in large part with the work of Gaiotto, Moore, and Neitzke [34]. Reference [35] analyzed the semiclassical description of magnetic defects in supersymmetric gauge theory and the connection to singular monopoles. We refer the reader to [35] for details and further discussion of the results reviewed here.
A supersymmetric 't Hooft defect placed at position ν J ∈ R 3 is specified by a charge P J and defined by imposing the singular boundary conditions where the ellipses represent subleading terms. By making local gauge transformations, the charge P J can be taken to be a constant element of the Lie algebra, valued in a Cartan subalgebra. For so(3) we take this Cartan subalgebra to be generated by T 3 , so that P J is specified by a single integer: Here "Tr" denotes a positive-definite Killing form on the Lie algebra. We work in the minimal twodimensional representation for so (3) where T a = − i 2 σ a with σ a the Pauli matrices. Then Tr is the negative of the ordinary matrix trace.

JHEP04(2021)286
Hence 't Hooft defects may be thought of as Dirac monopoles, where the 't Hooft charge, P , specifies an embedding of the magnetic charge into the non-abelian gauge group. Dirac quantization restricts p J to be an integer. Since gauge transformations can be used to send P J → −P J , it is only |p J | that is physical.
The action for supersymmetric Yang-Mills-Higgs theory in the presence of some number of 't Hooft defects placed at positions { ν J } is a sum of two types of terms. The first set of terms is referred to as the "vanilla" action in [35], and comprises those terms that would ordinarily be present for the theory without defects. The vanilla action depends on two parameters, the Yang-Mills coupling g ym and the theta angle θ ym . The second set of terms are boundary terms supported on infinitesimal two-spheres, S 2 J , surrounding the defects. These terms are required for consistency of the variational principle and preservation of supersymmetry. Letting ε J denote the radius of S 2 J and Ω J the solid angle, one has (2.9) The integrals over space in S van should be taken to exclude the infinitesimal balls bounded by the S 2 J 's surrounding the defect insertions. The θ ym term can be written as a total derivative and ordinarily does not contribute to the dynamics, but in the presence of defects the additional boundaries enable this term to become dynamical, as we will see below. When θ ym = 0, the 't Hooft defect is a source for the electric field and Y scalar as well, with these fields behaving as [35]

Magnetic monopoles and the BPS bound
Setting the fermions to zero, the Hamiltonian, or energy functional, associated with the action (2.6) is

JHEP04(2021)286
with local energy density Here U is R 3 with the infinitesimal balls around the defects removed. This result for the Hamiltonian holds provided that the Gauss Law constraint (or A 0 equation of motion), is imposed. This constraint is preserved by the time evolution as a consequence of gauge invariance.
The conditions for energy-minimizing field configurations are exposed by rewriting the Hamiltonian (2.12) as a sum of squares. Using integration by parts, cyclicity of the trace, the Gauss Law constraint (2.14), and the Bianchi identity µνρσ D ν F ρσ = 0, one finds that (2.12) can be written as 6 where M is a boundary term receiving contributions from the two-sphere at spatial infinity: (2.16) Here, the V (J) def terms in (2.12) cancel boundary terms generated from integration by parts on the infinitesimal two-spheres surrounding the defects, so that only the two-sphere at infinity contributes to M .
Asymptotic boundary conditions can be chosen to ensure finiteness of the energy. We require the magnetic and electric fields to fall off like O(1/r 2 ) while the Higgs fields must become covariantly constant, mutually commuting, and must also commute with the O(1/r 2 ) terms of the electric and magnetic fields. By a suitable gauge transformation we can assume where γ m , γ phys e , X ∞ , Y ∞ are all constants, valued in the same Cartan subalgebra. The boundary term (2.16) can be evaluated in terms of the asymptotic data: then the mass-squared of the W -boson is m 2 = m 2 X +m 2 Y . Meanwhile, γ phys e is the (physical) electric charge in the system as measured by the flux of the electric field through the two-sphere at infinity. The notation γ e is reserved for the charge of the Noether current associated with global gauge transformations that preserve the vacuum. The two are different when θ ym = 0: see [35] for details. The magnetic charge, γ m , may include contributions from ordinary monopoles in addition to the magnetic singularities. The presence of such monopoles requires m X = 0. The allowed values of γ m are constrained by both topology and dynamics. The condition for any simple Lie group and set of 't Hooft defects was determined in [56], building on earlier works [38][39][40][41], and here we state the result for SO(3). Letting T 3 be the generator of the Cartan subalgebra defined by X ∞ , such that we have where the p J determine the 't Hooft charges, (2.5). The non-negative integer n m is the number of ordinary monopoles present in the system. The expression (2.15) implies the lower bound on the energy functional, for a given set of asymptotic boundary conditions. The bound is saturated when all of the squares in (2.15) vanish, leading to the Bogomolny-Prasad-Sommerfield (BPS) equations A solution to these equations and the Gauss Law constraint will automatically solve the full equations of motion. A convenient gauge choice for studying solutions to (2.24) is A 0 = Y , in which case the last three equations imply A i , X, Y are time-independent. This leaves only the first equation, which we recognize as Bogomolny's equation for magnetic monopoles [10], and the Gauss Law constraint. The constraint can be rewritten, using the latter three of (2.24), as a linear equation for Y in a background (A i , X) that solves the Bogomolny equation:

JHEP04(2021)286
Thus we see that the Bogomolny equation arises as an energy minimizing condition. One is interested in solutions to these equations modulo gauge transformations that preserve the condition A 0 = Y . These are simply the time-independent gauge transformations. Finding solutions to the Bogomolny equation, with singularities of the form (2.4) and asymptotics of the form (2.17), is a well-studied problem going back to [39]. While explicit solutions are rare, a great deal is known about general properties of the space of gauge-inequivalent solutions. This space is referred to as the moduli space of singular monopoles. For G = SO(3), with 't Hooft defects carrying charges P J , (2.5), and an asymptotic magnetic charge given by (2.22), the moduli space of singular monopoles will be denoted M(n m , m X ; {p J , ν J }). We will often use the shorthand M when the context is clear. 7 M(n m , m X ; {p J , ν J }) is a 4n m -dimensional space that carries a natural Riemannian metric. The geometry of M will be discussed in the next subsection. Here we note that the dimension can interpreted as follows. A point in M represents a nonlinear superposition of n m ordinary monopoles in the presence of the defects. Each monopole has four moduli associated to it: three for its position in physical three-space and a fourth whose conjugate momentum corresponds to an electric charge that each monopole can carry. We let R n , n = 1, . . . , 4n m denote local coordinates on M, and we write the family of solutions to the Bogomolny equation as Given a solution (A i , X) to the Bogomolny equation and a boundary value Y ∞ , there will be a unique solution to the secondary BPS equation, [35]. Hence, by (2.24), the electric field and therefore the electric charge, γ phys e , will be determined. As we move in M, the electric charge will change. In other words, the electric charge is a function on monopole moduli space, determined by Y ∞ . Thus, if we fix an electric charge, solutions to (2.24) carrying that charge lie in a subspace of the moduli space defined by a level set γ phys e = γ phys e (R n ). Next we turn to a discussion of moduli space geometry and monopole dynamics via the collective coordinate ansatz. As we will explain, in this context it is more natural to proceed from solutions to the Bogomolny equation only, without imposing the secondary BPS equation for Y . The effects of Y will instead be felt through a potential energy on the moduli space, and possible values of the electric charge will be realized as constants of motion for the moduli space dynamcs.

Moduli space geometry
The kinetic terms in the action (2.7) for (A i , X) specify a metric on the infinite-dimensional space of field configurations through the identification of a tangent vector at the point (A i , X). This metric is the standard flat one: (2.27) and it induces a metric on the moduli space of singular monopoles. (The factor of 1/(2π) turns out to be a convenient normalization; see [35] for details.) To determine the moduli space metric we need a set of tangent vectors to solutions (A i , X) of the Bogomolny equation that will generate motion along the moduli space. These tangent vectors should therefore correspond to By differentiating the Bogomolny equation with respect to the moduli R n appearing in (2.26), one finds that (∂ n A i , ∂ n X), where ∂ n ≡ ∂ ∂R n , solves the linearized equations. However we must additionally require that the tangent vector (δA i , δX) be orthogonal to local gauge transformations, since the moduli space is the space of gauge-inequivalent solutions. This is achieved by requiring g(δ, δ ) = 0 for all , where δ = (−D µ , [ , X]) is an infinitesimal gauge transformation. The configuration (∂ n A i , ∂ n X) can be adjusted to solve this constraint by shifting it by a local gauge transformation: This condition ensures that g(δ n , δ ) = 0. Then the components of the moduli space metric with respect to the local coordinates R n are If one has an explicit family of solutions to the Bogomolny equation, (2.26), one can in principle compute the tangent vectors (2.28) and determine the metric directly from this definition. This was carried out in [42] for the family of solutions describing a single SO(3) monopole in the presence of defects studied in this paper. Typically, however, in cases where the metric is known, it is obtained from other mathematical representations of the moduli space, and there is a great literature on the subject.
Away from singular points the metric is hyperkähler. Co-dimension four singularities can exist in the moduli space of singular monopoles and are related to the phenomenon of monopole bubbling [41], in which an 't Hooft defect emits or absorbs a smooth monopole, changing the asymptotic magnetic charge of the system. In all known examples the singularities are of a fairly benign orbifold type. Furthermore, if all 't Hooft defects are taken to be minimally charged, |p J | = 1, then monopole bubbling does not occur, and the moduli space is smooth.
Another geometric construction on moduli space we will require is the Killing vector fields induced by global gauge transformations that preserve the Bogomolny equation and asymptotic data. Such Killing fields are in one-to-one correspondence with the Cartan subalgebra, t, of the gauge group. These vector fields generate isometries -in fact they

JHEP04(2021)286
are tri-holomorphic, generating isometries that preserve the hyperkähler structure as well the metric. The map is constructed as follows. Let H be the unique solution to Since we are restricting to gauge group G = SO(3) in this paper, the Cartan subalgebra is one-dimensional and there will be a single linearly independent Killing field generated corresponding to the action of global gauge transformations. Specifically, since exp(2πT 3 ) = 1 ∈ SO(3), G(T 3 ) will be a Killing field that generates a 2π-periodic isometry.

Collective coordinate dynamics on monopole moduli space
The metric (2.30) and Killing fields (2.31) play a central role in the collective coordinate description of monopole dynamics. The basic idea of the collective coordinate ansatz [15] is that time-dependent solutions of the full field theory describing monopole dynamics should be well-approximated by motion on moduli space -that is, by allowing the moduli to become functions of time (collective coordinates). Intuitively, this should be true provided the collective coordinate velocities are small, since the moduli space is a minimum-energy surface in the space of field configurations. Quantifying this condition requires some care. The reason the question is subtle is that field fluctuations around the monopole include modes of arbitrarily long wavelength for components of the fields along the u(1) preserved by the vev. In other words, there is no mass gap in the spectrum of fluctuations, and so energy can freely leak into radiation. Nevertheless, radiation is sourced by accelerating charges and one might expect energy loss to be small if the time variation of the collective coordinates is small.
In the context of classical time-dependent solutions in ordinary Yang-Mills-Higgs theory, the following mathematical result has been obtained by Stuart [16]. Suppose the collective coordinates are slowly varying, such that time derivatives behave as ∂ n t R = O( n ) for n = 1, 2, 3, with a small parameter. Given an initial field configuration that is close to the moduli space, such that the distance from the moduli space with respect to the metric (2.27) is O( 2 ), then the exact time-dependent solution to the Yang-Mills-Higgs equations will stay O( ) close to a model trajectory, This result is consistent with a physical estimate of the energy lost to radiation over the same time scale. Following [43], we evaluate the fields on the collective coordinate ansatz (2.26) with the model trajectory R (0) (t) and consider the time-dependence of the asymptotic multipole expansion for the massless u(1) components. In general, since the monopole terms -both magnetic and electric -are time-independent, the leading con-JHEP04(2021)286 tribution comes from dipole radiation. 8 The electric and magnetic dipole moments of the asymptotic fields can depend on the collective coordinates and are thus time-dependent. Since dipole radiation has a total radiated power of order P rad ∼ (∂ 2 t d) 2 , where d is either the magnetic or electric dipole moment, (see e.g. [44]), the rate of energy loss is O( 4 ). One can view the effects of this energy loss as a radiation reaction force of O( 3 ) acting on the system. Over time scales T = O(1/ ), one expects such a force to cause a deviation in the trajectory of O( ), and this is consistent with the theorem in [16].
A third point of view on the limits of the collective coordinate approximation arises in the context of quantum Yang-Mills-Higgs. In the semiclassical approximation to soliton states in quantum field theory, it is natural to take the collective coordinate velocities to be the same order as the Yang-Mills coupling, g ym , which is assumed to be small. Hence g ym plays the role of in the above discussion. With this identification one ensures that quantum corrections from field fluctuations around the monopole are suppressed relative to the leading collective coordinate dynamics and can be treated perturbatively. This perspective goes back to the original work on soliton quantization (see e.g. [45,46] for the diagrammatic approach). One can then define an effective Hamiltonian for the collective coordinates in the n-monopole sector by path-integrating out the field theoretic fluctuation fields around the background configuration (2.26).
In the leading saddle-point approximation to this path integral, one is solving the classical equation of motion for the fluctuation field and inserting this solution back into the field theory action to arrive at an effective action for the collective coordinates. This process can be carried out order by order in the small velocity expansion. At zeroth order in time derivatives of the collective coordinates one finds the (classical) mass of the soliton, which is O(g −2 ym ). There are no terms at first order in time derivatives because the static soliton is an exact solution. At second order in time derivatives, corresponding to O(g 0 ym ) terms in the Hamiltonian, one recovers the standard two-derivative collective coordinate Hamiltonian, whose equations of motion reproduce the model trajectory R (0) (t). One also recovers the first quantum correction to the soliton mass, which is independent of the collective coordinates. The first effects of the coupling between collective coordinates and radiation modes enter the effective Hamiltonian for the collective coordinates at third order in time derivatives, corresponding to O(g ym ). It is these terms that provide the explicit radiation reaction force discussed above, 9 and the same conclusion applies: over time scales T = O(1/ ) these terms will lead to a deviation from the model trajectory of O( ).
All of these approaches consistently show that, in the slowly-varying regime, the trajectories resulting from Manton's collective coordinate approximation for Yang-Mills-Higgs theory remain O( ) close to the true trajectories through times of O(1/ ). The advantage of the third approach, based on the collective coordinate effective Hamiltonian, is that it has JHEP04(2021)286 been extended to supersymmetric Yang-Mills-Higgs with the secondary Higgs field Y and its relation to electric charge [17,18,[49][50][51][52][53][54][55], and with the inclusion of 't Hooft defects [35]. We recall the key insights and results of this extension now.
As noted previously, the secondary BPS equation in (2.25), and the equation E i = D i Y , imply that the electric charge depends on the vev m Y and the point in moduli space. A analysis of this constraint shows that having configurations with electric charges q O(g −1 ym ) requires a hierarchy of scales m Y /m X ∼ O(g ym ). Hence, one should treat Y, A 0 on the same footing as the collective coordinate velocities. Specifically, one makes the ansatz and solves the remaining equations of motion for A 0 and Y in this background, working perturbatively in g ym , under the assumption thaṫ Upon inserting these expressions back into the action, one can integrate over space and, using the definition of the metric (2.30), one finds that the field theory action reduces to a particle mechanics action for the collective coordinates, R n (t).
This calculation was carried out in detail in [35] allowing for the presence of 't Hooft defects, and here we simply quote the results. In fact, we will only give part of the results since we are not considering the dynamics of the fermionic degrees of freedom in this paper. Index theory can be used to show that the fermions also carry 4n m massless real degrees of freedom in the monopole background; these are the superpartners of the bosonic collective coordinates. While this structure is essential for understanding the correct quantum mechanical model for the collective coordinates, it plays no role in the classical dynamics. Hence, setting the fermionic degrees of freedom to zero, one finds the following expansion for the field theory Lagrangian (2.6) around the monopole background: The terms in the first line, (2.34), are organized by scaling in g ym . In units of the vev, Here we are using (2.33).
are the bosonic pieces of two separate supersymmetry invariants identified in [35]. The terms comprising L (0) c.c. form the bosonic part of a collective coordinate Lagrangian that was first obtained for monopoles without defects in supersymmetric Yang-Mills-Higgs theory in [17,18]. In particular, they feature a potential energy term given by the norm-squared of the Killing field G(Y ∞ ). Thus we see how the Y Higgs field gives rise to a potential energy on the moduli space.
The final term, L c.c. , and its fermionic completion, were first obtained in [35]. This term is only dynamical in the presence of 't Hooft defects. When defects are absent,

JHEP04(2021)286
the Killing field G(X ∞ ) can be shown to be covariantly constant and L (θym) c.c. becomes a total time derivative. Since this term is O(g ym ), we should either drop it or write all contributions to the effective Hamiltonian at this order, which include the first higherderivative corrections to L (0) c.c. . In [35], the focus was on certain BPS trajectories and their quantum analogs, where supersymmetry can be used to argue that the higher-derivative corrections are inessential.
In this paper our interest is in generic collective coordinate motion, so we cannot make the same argument. We will nevertheless keep the terms in L (θym) c.c. . The reason is that, on the one hand, their effects are innocuous -modifying the definition of the canonical momenta below and adding a correction to the parameter that controls the strength of the moduli space potential energy. On the other hand, keeping these terms makes it easier to compare with [35], where they are important for matching onto predictions from the Seiberg-Witten description of BPS states [9]. 10 In particular, we recognize the constant term in L (θym) c.c.
as the remaining magneticcharge contribution to the BPS mass, (2.18), once γ phys e is expressed in terms of γ m and γ e using (2.20). The γ e contribution will instead be obtained from conserved momenta in the collective coordinate dynamics. We set and we find that the Lagrangian (2.34) leads to the conjugate momenta Note that the momenta are O(g −1 ym ) since the velocities are O(g ym ) while the mass is O(g −2 ym ). We note that the Hamiltonian can also be written in the form where we have introduced the combination and used the linearity of the G-map. The momenta and Hamiltonian are subject to O (1) and O(g ym ) corrections, respectively, coming from higher-derivative terms. As discussed JHEP04(2021)286 above, these higher derivative terms originate from the coupling of the collective coordinates to radiation modes in the full field theory. Thus, under the scaling assumptions (2.33), the coupling to radiation continues to be suppressed as it is in ordinary Yang-Mills-Higgs theory. This strongly suggests there should exist a direct analog of Stuart's theorem [16] in the supersymmetric context, with or without 't Hooft defects, for the moduli space with potential approximation. For such an extension of the theorem, the collective coordinate ansatz for all of the fields in the presence of 't Hooft defects would be the one given in subsection 4.3.1 of [35]. One does not expect the presence of defects to cause additional difficulties in the analysis since the linearized fluctuation operator controlling the radiation spectrum is sufficiently regular at the defect points: no special boundary conditions are required, and the modes are locally L 2 in a neighborhood of the defect points. Indeed, this was a key point in the analysis of [56] determining the dimension of the moduli space from a Callias index theorem. 11 The same procedure of solving the equations of motion perturbatively, as described under (2.32), leads to an expression for the electric charge, γ e , as function on moduli space [35,53]: In the semiclassical quantization of the collective coordinate dynamics, q is constrained to take integer values, since G(T 3 ) m generates a 2π-periodic isometry and the corresponding momentum eigenvalues are quantized. 12 In the classical theory, however, q can be any real number. Since the momenta are O(g −1 ym ) in the scaling regime we work in, it is natural to consider charges q O(g −1 ym ). Note such charges still lead to an electric field that is O(g ym ) according to (2.17) and hence suppressed compared to the magnetic field.
In the remainder of this paper we will analyze a class of solutions to the Bogomolny equation, describing one smooth monopole in the presence of any number of 't Hooft defects. We will construct simulations of monopole-defect interactions based on the corresponding moduli space geometry and collective coordinate Hamiltonian.

The BCD solutions
The solutions presented here were first obtained in [26,27] using a form of the Nahm transform [60] for singular monopoles developed in [30]. Later, the solutions were recovered in [28] from a modified Nahm transform referred to as the bow construction and developed in [29,61]. 11 No choice of self-adjoint extension is needed as it was in [57][58][59]. The difference between those references and the situation considered in [56] is that the background Higgs field on which the linearized fluctuation operator depends also has a singularity in the presence of an 't Hooft defect. This leads to a cancelation in the leading singularity of the operator analyzed in the earlier references. 12 For the supersymmetric Yang-Mills-Higgs theory discussed here, q would take on only even integer values due to the fact that all fields transform in the adjoint representation of the gauge group. In the notation of [35], q = −2ne.

JHEP04(2021)286
As above, r denotes the general position vector in R 3 , ν J are the positions of the singularities indexed by J = 1, 2, . . . , k t , and r J ≡ r − ν J . We let k t denote the total number of 't Hooft defects. Each defect is taken to be minimally charged, |p J | = 1, since a non-minimally charged singularity can be obtained by letting some of the ν J coincide. We denote the moduli corresponding to the smooth monopole's position by R. It will also be convenient to define R J ≡ R − ν J as the smooth monopole position relative to the J th defect and z ≡ r − R as the observation point relative to the smooth monopole. 13 These vectors are not all independent; in particular, r J − R J = z, for each J. The same letter without the arrow notation always denotes the magnitude of the vector: R J = | R J |, etc. The following quantities appear regularly in the following and are given special names: In terms of these quantities, the Blair-Cherkis-Durcan solutions are with the functions f, g given by The solutions are written in a hedgehog-type gauge, where spatial directions indexed by i, j = 1, 2, 3 are correlated with directions in the Lie algebra indexed by a, b = 1, 2, 3. We use Einstein summation conventions for repeated indices of type i, j and type a, b, but we always write the sum over defects explicitly. ijk is the totally antisymmetric symbol with 123 = 1. In the limit where all defects are sent to infinity, ν J → ∞, one sees that v → 0 and the terms involving f J , g J vanish in (3.2). Furthermore, the terms from the sums over J in f 0 , g 0 vanish, and (A ai , X a ) reduces to the Prasad-Sommerfield solution for the smooth monopole. As we approach the J th defect, the leading singularity of the Higgs field is evident from the 1/(2r J ) term in g 0 . 13 To compare with [27], let r → t, rJ → tJ , R → − T , RJ → − TJ . Additionally, due to a different normalization convention for the generators T a , (Aai, Xi) here = 2(Aai, Xa) there . For the same reason, m here X = 2λ there , and we set v here = 2zα there .

JHEP04(2021)286
Meanwhile, the asymptotic behavior of the Higgs field as r → ∞ can be extracted from the g 0 term: By making patchwise gauge transformations on the two-sphere at infinity, −ẑ a T a can be rotated to T 3 . Comparing with the asymptotic form, (2.17), (2.22), we see that n m = 1.
Hence this solution represents a single smooth monopole in the presence of the defects, as advertised.
The full family of solutions with n m = 1 depends on four moduli. The solution (3.2) exhibits dependence on three of these parameters, R ∈ R 3 , corresponding to the smooth monopole's position. Dependence on the fourth modulus, R 4 , can be implemented by acting on the configuration (A ai , X a ) with an asymptotically nontrivial gauge transformation, g = exp(R 4 X∞ ), where X∞ is the unique solution to D i D i X∞ + [X, [X, X∞ ]] = 0 satisfying lim r→∞ X∞ = X ∞ and regular in the interior. R 4 is a circle coordinate since the gauge group is compact. We will not need to carry this out explicitly, however. The reason is that the local energy density in the fields is a gauge-invariant quantity and therefore will be independent of R 4 . We turn to the computation of the energy density next.

Magnetic field and energy density
Since we are treating the effects of the second Higgs field and electric charge as a perturbation, we only consider the leading order contribution to the energy density, (2.13), due to the magnetic and primary Higgs field: where in the second and third steps we used the Bogomolny equation, B i = D i X. In this subssection we outline the computation of the components of the magnetic field, B ai = (D i X) a = ∂ i X a + abc A b i X c , and the magnetic energy density, E m , for the BCD solutions. The derivative ∂ i ≡ ∂ ∂r i acts on z and the r J , whereas R J is a constant. Thus, for example, The partial derivatives of the f 's and g's with respect to z and r K can be straightforwardly evaluated, but we suppress it here. Our main focus in the computation is to express ∂ i X a + abc A b i X c in a minimal set of tensor structures for the free a and i indices, and determine the scalar functions multiplying each of those tensor structures. We use the identities to eliminate all symbols. A minimal set of tensor structures can be taken as since we can use r J = z + R J to eliminate all appearances of r Ji and r Ja .

JHEP04(2021)286
The coefficient functions associated with the tensor structures (3.8) are denoted h, h 00 , h 0J , h J0 , and h JK , respectively, so that the magnetic field is (3.9) Tedious but straightforward computation yields The energy density of the BCD solution is obtained by squaring (3.9): Through the formulae of this section one thus obtains E m as a function of r ∈ R 3 , the 3(1 + k t ) parameters { R, ν J }, and the mass scale m X . As is clear from (3.3), m −1 X sets the natural length scale for the field configuration.
We have constructed a module in Mathematica that takes as input 1 + k t vectors,

JHEP04(2021)286
the surface, with the lowest value corresponding to the most transparent surface and the highest value corresponding to a completely opaque surface. This allows one to see "inside" the configuration. See figure 1 for two examples.
The code determines the values of the energy density to use based on a pre-sampling of values for the requested configuration. It attempts to ensure that the local maximum at the core of the smooth monopole lies between the values for the third and fourth surface, so that the smooth monopole remains semi-transparent. This however will not be possible if the smooth monopole is too close to a defect, such that there is not a well-isolated local maximum corresponding to its position. The energy density has a 1/r 4 J singularity as one approaches the J th defect, so the defects will always be accumulation points for the surfaces.
The code is denoted "EnergyPlot" in the Mathematica notebook included with the supplementary material of this submission. In addition to the position vector of the smooth monopole and a list (of arbitrary length) of position vectors for the defects, the code takes three further arguments: the number of initial plot points to use in the argument of Mathematica's RegionPlot3D, the size of the final image in pixels, and the position from which to view the configuration. (The output, however, can be rotated at will within Mathematica.) The examples in figure 1 used a relatively high value of 60 plot points and took 16 minutes and 5.5 hours respectively to render. In general, we expect the computation time to scale like k 4 t due to the quartic term in the last line of (3.11). Somewhat faster computations might be possible, especially for high k t values, by utilizing the identity Tr(D i XD i X) = ∂ i ∂ i Tr(X 2 ). The approach we have presented was motivated in part by the desire to have explicit expressions available that could be used to visualize the magnetic field itself.

Motion on multi-centered Taub-NUT
The motion of the smooth monopole in the presence of the defects is determined by the equations of motion following from the collective coordinate Hamiltonian, (2.40). This Hamiltonian requires the input of the metric on moduli space, g mn , and the tri-holomorphic The moduli space is the k t -centered Taub-NUT manifold [30], and the metric is known explicitly. Let H( R) be the harmonic function on R 3 \ { ν J } given by with | · | the standard Euclidean norm. 14 Let In this subsection we use the notation ∂i ≡ ∂ ∂R i . There should be no confusion as the coordinates on the physical space, r, do not appear in this subsection.

JHEP04(2021)286
Here, as always, the flat Euclidean metric is used to raise/lower indices of type i, j. Then the metric on the four-dimensional multi-centered Taub-NUT space takes the form (3.14) Here R 4 is a circle coordinate with periodicity R 4 ∼ R 4 + 4π, and the manifold restricts to a circle bundle over R 3 \ { ν J }. As R ≡ | R| → ∞, the size of the circle remains finite, and the overall normalization of the metric can be fixed by comparing to the definition (2.30) in this limit [35,42]. As one approaches a nut point, R → ν J , the circle fiber shrinks to zero size such that the total space is smooth. The vector field ∂ 4 ≡ ∂ ∂R 4 that generates motion along the circle fiber is, up to rescaling, the only tri-holomorphic Killing field. It follows from the periodicity of R 4 that G(T 3 ) = 2∂ 4 . Hence by linearity of the G-map and (2.21), We then find the conjugate momenta and Hamiltonian from (2.38) and (2.40) to be and The m Y term provides a potential energy well in the vicinity of each defect, which can lead to bound motion. Meanwhile π 4 is the electric charge q, (2.42): Since m −1 X sets the natural length scale in the physical R 3 , we work with dimensionless position variables R i = m X R i and parametersν Ji = m X ν Ji so that We note that R 4 is already dimensionless. We also define a dimensionless time, dimensionless momenta, and a dimensionless parameter C according to

JHEP04(2021)286
The factors of g ym , together with (2.33), ensure that ∂ τR n ,π i , q, and C are all naturally O(1) quantities. The factors of 2 and π are for convenience. The expressions for the conjugate momenta now take the form while the dynamical equations are From the first equation we learn that the electric charge, q, is a constant of motion. The remaining equations determine the motion of the smooth monopole on R 3 . These equations are more conveniently expressed in terms of the shifted momentum variables which leads to where we used (3.13). We also note that the Hamiltonian expressed in the new variables is As we can see from Newton's equation, (3.24), there are four different types of forces at play.
• The first term on the right-hand side of (3.24) has the typical form of the magnetic force on an electrically charged particle. The magnetic field is the sum of monopole fields created by the defects and is given by the gradient of the harmonic function, H. The force is proportional to the electric charge q of the smooth monopole and vanishes if the smooth monopole carries no electric charge. Taking into account that the momentum, p i contains a factor of H, we see that this force falls off as inverse distance-squared from the defect. As the smooth monopole approaches the J th defect, it acts to push the monopole in a direction transverse to the plane containing the relative position vector R J and the smooth monopole's instantaneous velocity.
• The second term is a force due to the position-dependent effective mass of the smooth monopole, m ∼ H. It is most noticeable when q = 0 such that the smooth monopole is not prevented from reaching the defects by the electric charge potential barrier. As JHEP04(2021)286 the smooth monopole approaches a defect its inertia increases and becomes infinite at R J = 0. Conservation of H c.c. dictates that the speed of the monopole must vanish at this point, and hence it is a turning point of the motion.
• The third and fourth terms on the right-hand side of (3.24), proportional to C 2 and q 2 respectively, provide competing attractive and repulsive forces on the smooth monopole by each defect. The attractive force is mediated by the secondary Higgs field while the repulsive force is due to the electrical self-energy of the smooth monopole and originates from a coupling to the long range component of the Higgs field X [62]. We can also view this term on the same footing at the p i p i term since q is the momentum along the circle fiber of Taub-NUT.
In figure 2 we plot the value of the potential energy function, that appears in the Hamiltonian (3.25) on a two-dimensional plane containing two defects. The first plot shows the potential energy with q = 0 and the second plot shows the potential energy with q = 0. Nonzero electric charge gives rise to a potential barrier that prevents the smooth monopole from passing over a defect. The function H −2 is a nonnegative bounded function on R 3 that increases to the limiting value 1 along any ray to infinity. Hence it follows from (3.24) that bound motion can exist if and only if |C| > |q|. Furthermore, if q = 0, nontrivial bound motion -i.e. other than the static solution R = constantrequires |C| > 0. Equations (3.23) and (3.24) can be numerically integrated for the smooth monopole's trajectory once the initial position and velocity are specified. We use the trajectory together with the energy density plots described earlier to construct simulations of the smooth monopole interacting with defects. Some examples are described in the next subsection.
With only a single defect, additional symmetries enable one to integrate the equations analytically. In fact, the same set of equations was studied in a different context in [33], JHEP04(2021)286 where it was shown that the general trajectory is a conic section. We review and extend this analysis in section 4. These analytic results inform the discussion of the various forces above.

Simulations
We have written several pieces of code in Mathematica for constructing simulations of monopole motion in the presence of defects. They can be found in the Mathematica notebook included with the supplementary material in this submission. Brief descriptions of the code and an illustrated example are included in that notebook. Additionally, the supplementary material contain four high resolution simulations. Two of these simulations depict bound motion in the presence of a three-defect system, one without electric charge and one with electric charge, and are described here. The other two movies are a scattering simulation and an oscillating simulation that displays complete screening of the defect. They are described in subsections 4.3 and 4.2 respectively. We integrate the equations of motion, (3.23) and (3.24), numerically to determine the smooth monopole's position and momentum as a function of time. The inputs are the initial conditions R 0 , p 0 at τ = 0, a set of defect positions { ν J }, values for the electric charge q and the coupling constant C, and the final time τ max to integrate to. The resulting trajectory R(τ ) can be plotted or fed into the code used to produce the energy density plots in figure 1. The code producing the frames for the animation outputs a table of energy density plots. This table can then be exported as a .mov or .avi file using Mathematica's Export command. The examples included in the supplementary material are based on a sampling rate of 12 frames per unit of time.
In figure 3 we show two frames of an animation in the three-defect system with q = 0 and C = 3. The initial position of the smooth monopole is as in the top of figure 1, and the initial velocity is directed back towards the center of the lower two defects. The trajectory JHEP04(2021)286 the monopole follows is shown as well. The motion takes place in the x-z plane. Although the motion is bound, it is not periodic, as the plots in figure 4 show.
In the case of bound motion such as this, the moduli space trajectory depicted here cannot remain an accurate approximation to the true field theory dynamics for arbitrarily long times. The reason is that the motion involves continual acceleration which, in the full field theory, will lead to energy loss through radiation. This energy loss is not captured in the classical truncation to the collective coordinate mechanics that we have employed. As we discussed in subsection 2.4, the moduli space trajectory will remain O( ) close to the true trajectory for times t ≤ T ∼ O(1/ ) in units of the Higgs vev. The small parameter controls the time variation of the collective coordinates and is naturally identified with g ym in the semiclassical analysis of the quantum Yang-Mills-Higgs theory.
In terms of the dimensionless time τ = gym 2 √ π m X t introduced in (3.20), this result translates to range of times τ ≤ O(1), and one might worry whether the trajectory in figure 4 can be trusted over one approximate cycle, much less the full range indicated. This would indeed be an issue if we wished to consider the theory at a small but fixed value of g ym . However, g ym can be chosen arbitrarily small, and nothing we have done so far fixes its value. Thus, as long as we consider a fixed range of times τ ∈ [0, τ max ], where τ max does not scale with g ym as g ym → 0, then T = 2 √ π gymm X τ max is O(g −1 ym ). In this way, we may view a τ max of 20 or 100, as in figure 4, as reasonable.
Note that by sending g ym → 0 we are sending the collective coordinate velocities to zero, via (2.33). In this language, the observation of the previous paragraph can be phrased as follows. The effects of radiation on the trajectory over any fixed length of rescled time, τ ∈ [0, τ max ], can be made arbitrarily small, so long as we are willing to consider arbitrarily slowly moving monopoles with respect to the true coordinate time, t. The slowness of these monopoles does not affect the appearance of figure 4 or the animations, since they are computed with respect to the rescaled time, τ . From this point of view, one should regard the animations as highly sped-up versions of the "true" motion.
In figure 5 we show two frames of a simulation in which the monopole carries electric charge q = 1. All other parameters and initial conditions are chosen to be the same as in the first simulation. The additional magnetic force on the charged monopole causes it to JHEP04(2021)286 veer outward along the y direction. The value of C is large enough and the energy low enough, however, that the monopole is drawn back towards the defects and remains bound to them. We use a different algorithm to determine the level sets of the energy density that will be plotted for the animations versus the individual configurations like those shown in figure 1. This generally results in the smooth monopole being rendered as a single semitransparent surface in the majority of frames of the animation. The reason is the following. In order for the animation to be an accurate representation, the level set values that are used to plot the energy density surfaces should remain fixed from frame to frame. As the smooth monopole moves away from the defects the overall value of the energy density in its core decreases. The algorithm ensures that the lowest value for the energy density in the core is always below the lowest level set so that the monopole never disappears. The remaining level sets increase in regular steps such that the highest one is just below the greatest energy density value that occurs in the smooth monopole's core over the duration of the simulation. 15 Therefore, multiple surfaces in the smooth monopole tend to only be evident when it is near a defect.
Each frame in these simulations takes a little over ten minutes to render on a 2018 MacBook Pro with a 2.2 GHz processor and 16 GB of RAM. Thus, at 12 frames per time unit for 17.5 units, each of these animations took about 1.5 days to finish.

Analytic results for a single defect
Observe that this energy can be written in either of the forms which implies the lower bound The bound can only be saturated if the smooth monopole is stationary.
When only a single defect present, such that the harmonic function takes the form there are additional conserved quantities. Here we have used translation invariance to place the defect at ν = 0, and we also allow for the possibility of an arbitrary defect charge, |p| = k t . This system is equivalent to a model for dyon interactions first considered by Zwanziger [31] and has been encountered in the context of ordinary monopole moduli space dynamics [32]. Reference [33] also provides a recent and detailed treatment. The additional constants of motion are the angular momentum vector 5) and the Runge-Lenz vector Here × denotes the usual cross product for Euclidean three-vectors, andR = R/R is the unit vector in the direction of R. The angular momentum vector receives contributions from the motion of the smooth monopole and from the angular momentum in the electromagnetic field. The strength of the latter is equal to the Dirac-Schwinger-Zwanziger pairing of the electric and magnetic charges of the monopole and defect.
From the field theory perspective, the origin of the conserved angular momentum is the fact that rotations of the field configuration about the defect map a solution of the Bogomolny equation satisfying all boundary conditions to a new solution, and thus generate a corresponding set of rotational isometries on the moduli space. The origin of the Runge-Lenz symmetry is due to the extended supersymmetry inherited by the collective coordinate dynamics when the fermions are included in the analysis. However these additional symmetries can also be understood purely from the point of view of the Hamiltonian particle mechanics, where they are realized as symmetries of the associated six-dimensional phase space. This point of view is explained nicely in [33], where a relationship to the Kepler problem and its Runge-Lenz vector is also discussed. We refer the reader there for further details.

Trajectories
References [32,33] showed how the conserved charges lead to a determination of the trajectories as conic sections. We review their analysis in this subsection for completeness and since our conventions are slightly different.
First, observe from (4.5) that is a constant. This implies that motion takes place on a cone. The axis of the cone is sgn(q) J and its opening angle, θ, is given by (4.8) If q > 0 then J is along the axis of the cone, if q < 0 then − J is along the axis of the cone, and if q = 0 the "cone" is the plane orthogonal to J. Since q and J are conserved, so is the magnitude | R × p|. However the direction of R × p will in general change along the trajectory.
Next, consider the consequences of the conserved quantity K. One finds that Since R · J = Rk t q it follows that Hence, defining the conserved vector, we have that This is the equation for a plane with outward normal vector sgn(q) N and distance to the origin k t |q|(J 2 −(k t q) 2 )/N . We will compute the magnitude N as well as the angle between J and N , but first we discuss some degenerate cases.

Motion along a ray
If R × p = 0 then the motion takes place along a fixed ray. Conservation of K ∝R implies that the smooth monopole cannot pass through the defect or elseR would flip sign. As we discussed under (3.26), bound orbits can only exist if |C| > |q|. Since the asymptotic value of the potential energy is U → 4(q 2 + C 2 ), we see that the one-dimensional motion will be bounded, and hence oscillatory, when additionally E < 4(q 2 + C 2 ). If E ≥ 4(q 2 + C 2 ) the motion will have a single turning point. When q = 0 the (inner) turning point will be at the origin, where the mass of the smooth monopole becomes infinitely heavy. When |q| > 0 the turning point will be at some distance away from the origin. We will find the turning point(s) of the motion below when we analyze the generic case. These formulae will include the case of motion along a ray as a special limit.

Motion in a plane containing the defect
If q = 0 but R × p = 0, then θ = π/2 and the motion takes place in the plane orthogonal to J, which contains the defect at R = 0. This coincides with the plane defined by N since N ∝ J when q = 0. The Runge-Lenz vector K is however interesting in this case. One finds that K · J = 0 when q = 0, and hence K lies in the plane of motion. One finds that and hence, defining φ as the angle measured counterclockwise from K, we have This is the equation for a conic section, with semi-latus rectum and eccentricity given by The magnitude of K can be computed by making use of the energy equation to eliminate p 2 . This is straightforward when q = 0 and yields This results in the eccintricity Hence for E < 2C 2 the trajectory is an ellipse, for E = 2C 2 a parabola, and for E > 2C 2 a hyperbola. The turning points and time dependence follow from setting q = 0 in the general case below.

The general case
Suppose now that q = 0 = R × p. In this case the plane defined by (4.12) does not pass through the origin, and the motion takes place on the intersection of this plane with the cone defined by (4.7). Thus the trajectory is still a conic section, but it occurs in a plane that does not contain the defect. To analyze the trajectory in greater detail we need expressions for the length N and the angle δ between the outward normal to the plane, sgn(q) N , and the axis of the cone, sgn(q) J. A more tedious computation of K 2 when q = 0 gives which leads to Hence the angle is Having obtained δ, we can now determine the conditions for the different types of trajectories. See figure 6 for reference. (The same figure appears in [33]; we reproduce it here for convenience.) First, if δ ≥ π/2 + θ there will be no intersection since the plane becomes parallel to the cone at δ = π/2 + θ. Then the three types of trajectory are Note in particular that the parabolic and elliptic cases require cos δ > 0 and therefore, from (4.22), E > 4q 2 . The critical values δ = π/2 ± θ correspond to cos δ = ∓ sin θ, and therefore by (4.22),

JHEP04(2021)286
To determine whether E should be greater or less than E c in each case, consider the functions which satisfy f ± (E c ) = 0. Then f + (E) > 0 is the condition for solutions to exist while f − (E) > 0 is the condition for bounded motion. Considering f ± (E), we find that f + is a strictly increasing function while f − is a strictly decreasing function. The discussion is divided into two cases. First suppose that |C| > |q|. Then the bound (4.3) implies that E > 4q 2 and hence that f + (E) > 0. Thus, (4.3) is sufficient to guarantee solutions exist. It follows from f − (E) > 0 that closed trajectories can only exist for E < E c . In order that this condition be compatible with E > 4q 2 we require 2(q 2 + C 2 ) > 4q 2 , which is guaranteed for |C| > |q|. Hence we recover the condition discussed under (3.26) for bounded motion to exist. When E > E c , we have that f + (E) > 0 while f − (E) < 0, and therefore the trajectory is hyperbolic. Now suppose that |C| < |q|. In this case E < E c implies E < 4q 2 . Thus f − is already negative and decreases as E increases, so it must be that f + is the function passing through zero at E = E c . Therefore, in this case, solutions only exist for E ≥ E c and are open trajectories. Hence the results may be summarized as follows: If |q| > |C| the first condition is empty, and elliptic trajectories cannot occur. These results are consistent with [33].

Explicit parameterization, turning points, and time dependence
Now let us be more explicit about the parameterization of these trajectories. By rotating our coordinate system, we can always assume that the z-axis is in the direction of sgn(q) J and that N lies in the x-z plane. Then θ, introduced in (4.7), is the polar angle in spherical coordinates, and the motion takes place at fixed θ. Writing R in spherical coordinates, the equation for the plane, (4.12), takes the form Since θ is fixed this is an equation relating R and φ which can be brought to a standard form (4.15) for a conic section. We still have the freedom of rotating our coordinate system about the z-axis by 180 degrees to make N x positive or negative. We correlate this choice with the sign of N z so that N x /N z is always positive. Having done so, we then obtain the equation with semi-latus rectum and eccentricity given by N z cos θ , e = tan θ| tan δ| . (4.29) JHEP04(2021)286 of the motion. In the case of parabolic motion, α is positive, and R = α/2 is the distance of closest approach. In the case of hyperbolic motion, e > 1 and α can have either sign. One root is positive and the other is negative, but which is which depends on the sign of α. The positive root is the distance of closest approach. Setting the results can be summarized as follows: The time dependence is then determined by elliptic: where the sign choice corresponds to the outward or inward part of the trajectory respectively, and the required indefinite integrals are As in the Kepler problem, the equation for R as a function of τ is transcendental. We can, however, obtain an analytical expression for the period of an elliptical orbit: The equations for the turning points (4.36), the time dependence (4.38), and the period (4.40) are valid in all cases, including the case of one-dimensional motion along a ray.
In particular, T gives the oscillation period in this case.
The supplementary material includes an animation of an oscillating solution with k t = 2, q = 0, and C = 2. When q = 0 the inner turning point is at R − = 0, atop the defect. Furthermore, when the defect charge is k t = 2, the smooth monopole can completely screen the defect: at R = 0 the solution becomes the trivial one with vanishing energy density. The vanishing of the asymptotic magnetic charge (2.22) when k t = 2 is consistent with the possibility of this configuration being a point in the moduli space. However, the evolution is smooth through this singularity and the monopole and defect reemerge. In figure 7 we JHEP04(2021)286 show three frames from the animation and a plot of the monopole-defect distance as a function of time.

Scattering off the defect
In this final subsection we analyze the scattering problem for the smooth monopole off of the defect at the origin. Reference [32] previously gave the differential scattering cross section for a mathematically equivalent problem by generalizing the results in [63] to the case a nonvanishing attractive potential (i.e. C = 0 in our language). We take a different and elementary approach based on the classical trajectory and a rotation between reference frames.
Scattering is easy to analyze in the adapted coordinate system of figure 6. Motion takes place at a constant value of the polar angle θ, while the projection of the trajectory into the x-y plane is a hyperbola that starts in the second quadrant and ends in the third quadrant or vice versa. The asymptotic initial and final angles, φ ± , are the two solutions to cos φ ± = −1/e , (4.41) in the range (π/2, 3π/2). See figure 8. Therefore our approach is the following. We set up some initial data -an incoming velocity and an impact parameter -for the scattering problem in the "lab" frame, compute the values of the conserved quantities, E, J, N , and rotate the incoming direction to the adapted coordinate system, where it must match up with one of n ± =    cos φ ± sin θ sin φ ± sin θ cos θ

JHEP04(2021)286
The other of these two is then the outgoing direction in the adapted frame, which we finally rotate back to the lab frame to obtain the outgoing direction in the lab frame. We will refer to the lab frame as the primed frame in the following. Without loss of generality, we take our initial data to be with v > 0 and b ≥ 0, so that the smooth monopole is coming in parallel to the z-axis, a distance b from it in the direction of the positive x-axis. Using these, we obtain the conserved quantities where we note that sgn(N x ) = sgn(q) , sgn(N y ) = sgn With the transformation from lab to adapted coordinate system in hand, we can determine the incoming direction,n in =k in the adapted coordinate system:

JHEP04(2021)286
As a check, one can use the above results to verify that (n in ) x = − 1 e sin θ . (4.50) We find thatn in =n ± , (4.42), for sgn(N x ) sgn(N y ) = ∓1 respectively. Regardless, the outgoing direction in the adapted coordinate system is thereforê (4.51) Finally we rotate this vector back to the lab frame to determine the outgoing direction of the smooth monopole after the interaction: We find the following components: . (4.53) Some comments on the result are: • If q = 0 the y-component vanishes, and this is consistent with the fact that the scattering should take place in the x-z plane, which is the plane containing the defect. The sign of q determines whether the monopole scatters above or below the initial plane of motion.
• If E − 4q 2 = 1 2 v 2 + 2C 2 − 2q 2 is positive, then the x-component is negative, meaning that the trajectory bends around the defect, while if it is negative then the trajectory bends away from the defect. Recalling that E = 4q 2 is the condition that δ = π/2, this behavior is consistent with figure 6 and the comments under (4.30).
• The sign of the the z-component can always be made negative by choosing b large enough. This makes sense: if b is large, then the trajectory shouldn't be much affected by the defect, and hence the direction of the monopole's final velocity should be close to that of the initial velocity.
• If b = 0 the outward direction simplifies ton out =k. The smooth monopole approaches the defect, comes to a stop, and reverses. This corresponds to the case of motion along a ray, discussed earlier.
The supplementary material includes an animation of a scattering trajectory. In figure 9 we show three frames from that animation and illustrate how the outward trajectory becomes parallel to the direction specified by (4.53).

Conclusions
In this work we have analyzed, both numerically and analytically, the interactions of a BPS monopole with an arbitrary number of 't Hooft defects. Our motivations were to 1. broaden our understanding of classical soliton dynamics in the presence of defect singularities, and 2. illustrate with a new class of examples the emergence of particle dynamics from field theory through the collective coordinate paradigm for solitons.
Our main numerical results consist of simulations built on two key inputs. First, the monopole and defect positions are represented in three-dimensional plots based on JHEP04(2021)286 the energy density of the fields, determined analytically from the Blair-Cherkis-Durcan solutions. We plot several level sets of the energy density with varying opacity, illustrating finiteness of the density in the core of the smooth monopole and divergences in the cores of the defects. Second, motion of the smooth monopole is generated by numerical integration of the equations of motion determined from the collective coordinate reduction to monopole moduli space. The Mathematica code developed to produce the animations, as well as several example movies have been included with this submission as supplementary material.
In section 4 we explored the case of a single defect analytically, building on the work of [32,33]. We determined the period of an elliptical orbit, and a provided a new and elementary analysis of the scattering problem.
It would be interesting to extend the numerical techniques of this paper to the case of multi-monopole interactions in models based on a higher rank gauge group. In such theories, monopoles come in different types because there are distinct types of magnetic charges they can carry -as many as the rank of the gauge group [64]. Furthermore, the classical solutions and moduli space geometry for multi-monopole configurations with constituents of distinct type are much more tractable than for multi-monopole configurations carrying only one type of magnetic charge. See, e.g. [52,65]. Simulating the dynamics for generic initial conditions should be possible, and this is an area in the field of magnetic monopoles that has not yet been explored.