Dynamical evolution of $U(1)$ gauged Q-balls in axisymmetry

We study the dynamics of $U(1)$ gauged Q-balls using fully non-linear numerical evolutions in axisymmetry. Focusing on two models with logarithmic and polynomial scalar field potentials, we numerically evolve perturbed gauged Q-ball configurations in order to assess their stability and determine the fate of unstable configurations. Our simulations suggest that there exist both stable and unstable branches of solutions with respect to axisymmetric perturbations. For solutions belonging to the stable branch, the gauged Q-balls respond to the perturbations by oscillating continuously or weakly radiating before returning to the initial configuration. For the unstable branch, the solutions are eventually destroyed and can evolve in several ways, such as dispersal of the fields to infinity or fragmentation into smaller gauged Q-balls. In some cases, we observe the formation of ring-like structures which we call"gauged Q-rings". We also investigate the stability of gauged Q-balls when the gauge coupling is small, finding that the behaviour of these configurations closely resembles that of ordinary (non-gauged) Q-balls.


INTRODUCTION
Solitons are a fundamental prediction of many physical theories. They are characterized as stable, localized solutions to non-linear field equations that behave in many ways like particles. Broadly speaking, solitons can be classified as either topological or non-topological. Topological solitons owe their existence and stability to the specific topological constraints of a given model. Nontopological solitons, in contrast, can arise simply due to the balance of attractive and repulsive self-interactions in the theory. In addition, the stability of non-topological solitons is often associated with the presence of conserved charges which emerge from the theory's underlying symmetries (though one can also construct solitonic configurations in the absence of such charges [1]).
Perhaps the simplest example of a non-topological soliton in field theory is the Q-ball: a stable, localized solution of a complex scalar field theory with a non-linear attractive potential and a global or gauge U (1) symmetry. In recent years, Q-balls have attracted significant attention due to their prevalence in supersymmetric theories [2] and their possible cosmological consequences. In particular, it has been suggested that Q-balls may be relevant for baryogenesis [3,4], cosmological phase transitions [5,6], and the dark matter problem [7,8]. The formation of Q-balls could also lead to detectable gravitational wave signatures [9]. However, regardless of their physical applications, Q-balls are also interesting from a theoretical perspective as stable, particle-like objects that can be constructed from smooth classical fields and that have vanishing topological charge.
The properties of Q-balls under a global U (1) symmetry have been studied extensively in the literature. Starting with the work of Rosen [10], Q-ball solutions have * mikin@physics.ubc.ca † choptuik@physics.ubc.ca been found in a variety of physically-motivated models using various scalar field potentials (see [11] for a recent review). For some special potentials, the equations can be solved exactly [12][13][14][15], but in the general case one must use approximations or numerical methods in order to determine the characteristic features of Q-balls. Associated with each solution in a given model, there is a conserved energy E and a conserved Noether charge Q (from which the Q-ball gets its name) corresponding to the particle number. Each solution is also characterized by an internal oscillation frequency ω which can be interpreted as the chemical potential of the configuration [16]. In addition to ordinary (ground state) Q-balls, one can construct excited Q-balls which have additional radial nodes or non-zero angular momentum [17][18][19][20][21]. The basic theory has also been extended by coupling Q-balls to gravity [22][23][24], by introducing a massless or massive gauge field [25][26][27], and by considering non-spherical configurations such as Q-tubes [28], Q-rings [29], and composite systems of Q-balls [30].
When the global U (1) symmetry of the theory is gauged, the Q-balls acquire an electric charge and are known as gauged Q-balls [26]. Gauged Q-balls have properties that can differ significantly when contrasted to their global (non-gauged) counterparts. The presence of a massless gauge field introduces a long-range repulsive force that can destabilize the solutions for large gauge couplings. This repulsive force can give rise to novel scalar field configurations such as Q-shells [31][32][33], but it can also place limits on the maximum size and charge of gauged Q-balls for some scalar field potentials [32,34]. The existence of this maximal charge corresponds with the limits of the allowed range of the frequency ω, and in general the gauged Q-ball configurations cannot be uniquely characterized by the value of ω [34]. Despite these differences, there exists a correspondence which allows for the properties of gauged Q-balls to be approximated from the properties of non-gauged (global) Qballs, which are often much simpler [35]. In addition, when the interaction between the scalar field and gauge field is weak, gauged Q-balls are expected to closely resemble their non-gauged counterparts [36].
One of the essential properties of Q-balls relates to their dynamical stability. In order to be physically viable, solitons must be robust against perturbations. However, the problem of establishing the stability of solitons is often complicated by the non-linear nature of the governing field equations. In some cases, linear perturbation analyses and stability theorems can be applied to determine the expected regions of stability and instability.
For non-gauged Q-balls, it has been shown that the simple relation serves as an effective criterion for establishing regions of stability [17,37]. However, the case of gauged Q-balls is more complicated due to the presence of the repulsive gauge field. It was recently shown in Ref. [38] that the sign of dQ/dω cannot be used to assess the stability of gauged Q-balls in the general case. In the absence of such a criterion, one can still analyze the stability of these solutions using (among other alternatives) a numerical approach: dynamically evolving perturbed configurations through direct solution of the non-linear equations of motion. This method was applied in Ref. [38] to show that gauged Q-balls in several models can be stable with respect to spherical perturbations. However, it remains an open question as to whether gauged Q-balls can be classically stable against decay from more general perturbations beyond spherical symmetry. In addition, the instability mode for non-gauged Q-balls is always spherical [39], but it is not known whether gauged Q-ball decay can be mediated by non-spherical modes.
In this paper, we make progress toward understanding some aspects of gauged Q-ball dynamics by performing fully non-linear numerical simulations of the field equations in axisymmetry. There are two main questions we shall explore: (i) what is the range of stability of gauged Q-balls in axisymmetry? And (ii), what is the final fate of those configurations which are unstable? To answer these questions, we construct spherical gauged Q-ball initial data using a numerical shooting technique. We then assess the stability of these configurations by dynamically perturbing the system and observing the subsequent behaviour.
Numerical results presented below suggest that there exist both stable and unstable branches of solutions in axisymmetry. We find that stable gauged Q-balls, when perturbed, can survive over timescales which are long compared to the dynamical time with no evidence of measurable growing modes which destroy the configuration. These solutions respond to perturbation by oscillating continuously or weakly radiating before returning to the initial configuration. Unstable gauged Q-balls, in contrast, are typically short-lived and can decay in one of several ways. Some unstable solutions break apart into many smaller gauged Q-balls or shed scalar field until they relax into a smaller stable configuration. Other unstable solutions fragment into non-spherical ring-like structures which propagate away from the axis of symmetry and can survive for some time. In addition, for the case of a logarithmic potential we observe that the maximum magnitude of the scalar field can grow without bound. We interpret this behaviour as a consequence of the potential being unbounded from below. Finally, we test the effect of the gauge coupling strength on the stability, finding that gauged Q-balls closely resemble their non-gauged counterparts when the coupling is small. This paper is organized as follows: in Sec. II, we present the equations of motion of the theory. In Sec. III, we discuss the procedure for obtaining axisymmetric initial data which is used in the numerical evolutions. In Sec. IV, we briefly discuss the types of perturbations that are applied to the system. In Sec. V, we present the results for several representative evolutions. We conclude with some final remarks in Sec. VI.
Throughout this work, we use natural units where c = = 1 and employ the metric signature (−, +, +, +). We focus on unexcited gauged Q-ball solutions (those for which the scalar field modulus attains only one maximum). For brevity, we will use the term "Q-ball" interchangeably with "gauged Q-ball" when the distinction between the gauged and non-gauged solutions is obvious by context.

II. EQUATIONS OF MOTION
The Lagrangian density of the theory takes the form where φ is the complex scalar field, F µν = ∂ µ A ν −∂ ν A µ is the electromagnetic field tensor for the U (1) gauge field A µ , and D µ = ∇ µ − ieA µ denotes the gauge covariant derivative with coupling constant e. Here, V (|φ|) is a U (1)-invariant scalar field potential that permits Q-ball solutions in the limit e → 0. In this work, we consider the following scalar field potentials: where µ, β, m, k, and h are assumed to be positive, realvalued parameters of the potentials. The potential (3) has previously been studied in various forms in Refs. [12,32,36,[38][39][40][41] while potential (4) has been studied in Refs. [16,18,19,21,26,35,42,43]. Further details about the scalar potentials (3) and (4) will be discussed in the sections that follow. The evolution equations for the theory can be found by varying the Lagrangian density (2) with respect to the scalar and gauge fields to obtain From (6) we identify the conserved Noether current which corresponds with invariance of the theory under the U (1) transformations The conserved current (7) can be integrated to obtain a conserved Noether charge Q = j 0 d 3 x. Also associated with the theory is the energy-momentum tensor and the corresponding conserved energy E = T 00 d 3 x.
To solve the field equations of motion, we adopt the usual cylindrical coordinates (t, ρ, ϕ, z) and write the spacetime line element as In three spatial dimensions, computational constraints would limit the range of possible field configurations that could be explored. We therefore reduce the computational complexity of the problem by imposing axisymmetry on the system: no dependence of any of the fields on the azimuthal angle ϕ is assumed. This results in a system of six coupled non-linear partial differential equations which are described in App. A. Evolution of the fields is subject to the constraints where E and B are the electric and magnetic fields, respectively, and ρ c is the electric charge density Equation (13) will be trivially satisfied in axisymmetry while equation (12) can be re-expressed in terms of the gauge field components using the relation We also impose the Lorenz gauge condition to simultaneously fix the gauge and simplify the equations. It is expected that a numerical solution to the equations of motion will also satisfy the constraint equations at a given time.

III. INITIAL DATA
To generate suitable initial data, we make a spherically-symmetric ansatz for the scalar and gauge fields Inserting this ansatz into the equations of motion yields the following system of coupled equations: Here, V (f ) is the scalar potential and g(r) = ω − eA 0 (r). This system constitutes an eigenvalue problem for the parameter ω subject to the boundary conditions which are required to ensure finiteness of energy and regularity at the origin. Gauged Q-ball solutions can be found by solving the system of equations (20)-(21) using an iterative shooting technique [44] to simultaneously determine f (r) and A 0 (r). In this method, an initial choice is made for the value of g(0) and a corresponding guess is made for f (0). The equations are then integrated on a uniform grid using a fourth-order classical Runge-Kutta method out to a finite radius r 0 . Depending on the asymptotic behaviour of g(r) and f (r) at r 0 , the value of f (0) is adjusted through iterative bisection until the boundary conditions (22) and (23) are approximately satisfied at large r. Once a solution is found, the eigenvalue ω can be determined from the asymptotic value of g(r) using the boundary condition (23) and A 0 (r) can be determined as A 0 (r) = (ω − g(r))/e.
One of the main computational challenges associated with this procedure is the high numerical accuracy required in order to find satisfactory solutions. Typically, the number of digits required for convergence will exceed the capacity of double-precision floating-point numbers.
To overcome this limitation, we employ the arbitraryprecision arithmetic capabilities of Maple [45]. The software precision is adjusted and the integration is carried out until the asymptotic behaviour of f (r) is observed to decay exponentially at large r. At this point, the value of f (r) is typically very small (one part in 10 8 or smaller) and so the fields g(r) and f (r) approximately decouple in equations (20) and (21). In this asymptotic region, we fit a 1/r tail to g(r) and an exponentially-decaying tail to f (r) [34] so that the solution is determined to an arbitrarily large radius. In Fig. 1, we present the results of our shooting procedure for the logarithmic potential (3). For numerical purposes, we find it convenient to set µ = β = 1 in the model. The central scalar field value f (0) is plotted against g(0) = ω − eA 0 (0) for various values of e. When the value of e is small (representing weak gauge coupling), the curve of solutions is monotonically decreasing and single-valued, closely resembling the behaviour of non-gauged Q-balls. However, when e is increased, the curves are no longer single-valued and they begin to bifurcate with some curves ending abruptly in the solution space. These distinct endpoints generally correspond to the appearance of additional radial nodes in the solution, representing excited gauged Q-balls [17][18][19][20][21]. Also notable is the appearance of distinct curves close to the horizontal axis where f (0) is very small. These solutions correspond to Q-shells [31][32][33] which attain their maximal field values away from r = 0 and resemble shell-like concentrations of the fields. Fig. 2 are the results of our numerical shooting procedure for the polynomial potential (4) with m = k = 1 and h = 0.2. In order to clearly distinguish the curves, and following Ref. [42], we plot g(0) = ω − eA 0 (0) versus ω for various values of e. We restrict our shooting to solutions where ω ≤ 1, which is required in order to ensure that the solutions have finite energy [26,34,42]. The case of e = 0.0 (corresponding to non-gauged Q-balls) is represented by a linear line in the solution space. As e is increased, a minimal value ω min appears which separates each curve into an upper and lower branch. The value of ω min increases with e until ω min > 1, at which point no gauged Q-ball solutions can be found in the model. We note that while Q-shell solutions are known to exist for the polynomial potential [33], no such solutions are found for our choice of the potential parameters. As a basic check of our shooting procedure, we have compared our numerical solutions to those reported in previous publications on U (1) gauged Q-balls in logarithmic and polynomial models [38,42]. We find good agreement with the previously reported results.

Plotted in
In order to generate initial data that is suitable for evolution in axisymmetry, we transform the spherical solutions in (17)- (19) to cylindrical coordinates by performing a fourth-order polynomial interpolation of the spherical solution in the ρ-z plane. This provides axisymmetric initial data that will subsequently be used in our numerical simulations.

IV. DIAGNOSTICS
Here we discuss several diagnostics which are useful in characterizing the stability of each gauged Q-ball configuration. For the purposes of this work, a configuration is defined to be stable if small perturbations to the initial state remain bounded during the course of the evolution. Unstable configurations are those for which small perturbations grow exponentially on top of the solution, eventually leading to the destruction of the Q-ball (such as through fragmentation or dispersal of the fields). Note that with this definition, we classify as stable those con-figurations for which the fields may be weakly oscillating or radiating but are not destroyed by the initial perturbation.
There are several physical quantities associated with the scalar and gauge fields which are relevant when monitoring the evolution of each configuration. Principal among these are the conserved Noether charge Q and the total energy E. The Noether charge is given by and the total energy is given by Both Q and E are time-independent quantities which are expected to be conserved as long as the fields remain localized within the simulation domain.
We will now discuss how we add small perturbations to the stationary initial data. The solutions are perturbed in two ways: (i) perturbation through the inherent numerical truncation error of the finite-difference scheme, and (ii) perturbation by an auxiliary scalar field designed to explicitly excite all underlying modes of the configuration.

A. Perturbation by Numerical Truncation Error
As a first test of the stability of our gauged Q-ball configurations, we numerically evolve the fields forward in time using the axisymmetric initial data described in Sec. III. Upon evolution, the fields will be subject to small numerical perturbations due to the finite-difference discretization which is used to solve the evolution equations. This can be understood from the observation that the discrete solution of a uniform centered finite-difference scheme admits a truncation error expansion around the continuum solution in powers of the grid spacing [46]. While the exact form of this expansion is generally not known (making the perturbations effectively random), the magnitude of the associated truncation error is tied closely to the numerical resolution of the simulation and can therefore be indirectly controlled. In the sections that follow, we will refer to perturbations of this form as "Type 0".
One consequence of this type of perturbation is that any potential instabilities will take longer to manifest for higher-resolution simulations than for lower resolution ones. This is because the magnitude of the truncation error becomes smaller at higher resolutions. It is therefore necessary to evolve the configuration over sufficiently long times in order to assess stability. The notion of a "sufficiently long time" is difficult to make precise, but this timescale can generally be estimated by observing the oscillations of the scalar field modulus |φ| when subject to a perturbation. Even for small perturbations, the maximum value of |φ| will tend to oscillate at frequencies which correspond to the underlying modes of the configuration. The dynamical time can then be identified as the inverse frequency of the longest mode. For the problem at hand, we evolve each configuration with this timescale in mind so that any slowly-growing unstable modes have the opportunity to manifest. We find that this procedure provides an adequate preliminary test of stability which can be further verified using additional perturbation methods (to be discussed immediately below).

B. Perturbation by an Auxiliary Scalar Field
As a second test of stability, we dynamically perturb the gauged Q-balls by simulating the implosion of an asymmetric shell of matter onto the stationary configurations. We do this by introducing a massless real scalar field χ(t, ρ, z) that couples to the complex Q-ball field φ(t, ρ, z) in the modified theory: where U (|φ|, χ) describes the interaction potential of the fields φ and χ. We compute the modified equations of motion from (26) to obtain One can see from (27)-(29) that by choosing an interaction potential U (|φ|, χ) such that U (|φ|, χ) → 0 in the limit of χ → 0, then the modified equations (27) and (28) reduce to equations (5) and (6) (with (29) just representing an independent wave equation for χ). This means that χ and φ will elicit some mutual influence when the fields overlap, but the influence will disappear if the fields become well-separated. In this sense, χ can act as an external perturbing agent. We initialize χ as an ingoing asymmetric shell of the form where A, a, b, δ, r 0 , ρ 0 and z 0 are parameters specifying the characteristics of the initial pulse. If r 0 is made large, the field approximately vanishes in the vicinity of the gauged Q-ball at the initial time and so χ can be considered an external perturbation with a size controlled by A. Strictly speaking, the notion of an "external" perturbation cannot be made precise because gauged Q-balls do not have a finite radius. However, since the scalar field decays exponentially away from the center of the configuration, initializing the auxiliary field sufficiently far away from the center will serve as a good approximation to an external perturbation. Note also that the auxiliary field couples only to the scalar field so that the long-range behaviour of the gauge field is not a significant factor. During the evolution, the massless scalar field implodes toward the origin and collides with the gauged Q-ball. The two fields temporarily interact before the massless field passes through the origin and explodes outward to r → ∞, leaving the gauged Q-ball perturbed at the origin. Due to the asymmetry of the imploding pulse, the interaction of the two scalars is expected to excite the underlying modes of the system and induce the disruption of any unstable configurations. For our purposes, we choose where c is a coupling constant that controls the coupling strength between χ and φ. In the sections that follow, we will refer to perturbations of this form as "Type I". This technique resembles the methods of Ref. [47] to investigate the stability of boson stars. We note that configurations which are subject to perturbations of this type will inevitably also be perturbed by the inherent truncation error of the numerical simulation (Type 0). However, since Type 0 perturbations are typically very small and effectively random, Type I perturbations provide an additional level of control in determining the stability of a given configuration. For the results presented here, the simulations are repeated for various values of the Type I perturbation parameters A and c. This is done to verify that the response of the Q-ball field to the perturbation (as measured, for example, by the magnitude of the induced oscillations of the scalar field modulus |φ|) remains in the linear regime: an increase of A or c leads to a corresponding increase in the magnitude of oscillations of the perturbed |φ|.

V. NUMERICAL RESULTS
Here we present results from the dynamical evolution of gauged Q-balls in the potentials (3) and (4). For each simulation, we numerically solve the evolution equations using a second-order Crank-Nicolson finitedifference scheme implemented using the Rapid Numerical Prototyping Language (RNPL) framework [48]. Fourth-order Kreiss-Oliger dissipation is applied as a mild low-pass filter to damp poorly-resolved and potentially problematic (from a numerical stability viewpoint) high-frequency solution components. We implement a modified Berger-Oliger adaptive mesh refinement (AMR) algorithm via the PAMR/AMRD libraries [49] to increase the numerical resolution of our simulations. In all examples presented below, the base grid is taken to be 129 × 257 grid points in {ρ, z} and up to five additional levels of mesh refinement are used with a refinement ratio of 2:1. The domain is taken to be finite with outgoing boundary conditions imposed at the outer boundaries. Reflective (or anti-reflective) boundary conditions are imposed at the axis of symmetry in order to ensure regularity of the fields. We choose a Courant factor of λ = dt/ min{dρ, dz} = 0.25 and evolve the configurations until at least t ≈ 1000 to assess stability, though in many cases we evolve for longer in order to observe the late-time dynamics. Further details about the numerical implementation and code validation are given in App. B.
To illustrate the general behaviour of stable and unstable configurations, we focus on several specific solutions for the potentials (3) and (4). These solutions are listed in Table I. Configurations L1-L4 correspond to the logarithmic potential while configurations P1-P2 correspond to the polynomial potential. We note that besides L1-L4 and P1-P2, we have also performed hundreds of additional evolutions along the solution curves of Fig. 1 and Fig. 2 in order to determine the general regions of stability. This stability investigation will be discussed below.

A. V log Model
Here we consider the dynamical stability of gauged Qballs in the logarithmic model (3). Most relevant for this work are the results of Ref. [38] which conducted numerical evolutions of gauged Q-ball configurations for e = 1.1 in spherical symmetry. There it was found that both stable and unstable gauged Q-balls can exist in the model, though the classical stability criterion (1) provides little information about the stability of a given configuration in the general case. Once again, we set µ = β = 1 for numerical purposes. For brevity, and to facilitate comparison with previous work, we focus on the case of e = 1.1 where the system is fully coupled.
The relevant properties of each of the configurations L1-L4 are described in Table I along with the final result of numerically evolving the configuration forwards in time. In Fig. 3, the location of each of these configurations in the solution space is labelled with a dot. L1 corresponds to a solution on the stable branch. L2 corresponds to an unstable configuration which decays through dissipation of the fields. L3 corresponds to an unstable Q-shell solution which breaks apart into several smaller solitonic components. Finally, L4 illustrates the case of an unstable solution which responds to perturbation by growing without bound. Here, only L1 is subject to the Type I perturbation (to illustrate the dynamical stability of the configuration) while L2-L4 are all subject to Type 0 perturbations only.  (3). P1-P2 represent configurations found using the polynomial potential (4). The second column indicates the outcome of the numerical evolution. From left to right, the remaining columns give the initial central scalar field amplitude φ(0, 0), the initial central gauge field value A0(0, 0), the electromagnetic coupling constant e, the eigenfrequency ω, the total integrated energy E, the total Noether charge Q, the size of the simulation domain spanning {ρ : 0 ≤ ρ ≤ dmax} and {z : −dmax ≤ z ≤ dmax}, the type of perturbation used, and the perturbation parameters c and A (if applicable). First let us discuss the evolution of configuration L1. This evolution is run for 65000 base-grid time steps up to a final time of t = 6400. To assess the stability, we apply a Type I perturbation with parameters c = 0.1 and A = 0.1. Results for this evolution are shown in Fig. 4. The contour lines in the figure represent the scalar field modulus |φ| while the colormap represents the perturbing field modulus |χ|. The imploding pulse, which is centered around the point {ρ = 0.0, z = 0.5}, interacts with the Q-ball starting at t ≈ 20 (the second panel of the figure) and explodes through the origin, leaving the simulation domain at t ≈ 70. This induces small asymmetric distortions in the Q-ball field which can be observed as changes of the contour lines in the subsequent panels. This distortion also creates large oscillations in the maximal value of |φ| which are plotted in Fig. 5. Prior to the imploding pulse interacting with the Q-ball, the oscillations of |φ| are very small and are sourced by Type 0 perturbations. After the pulse interacts with the Q-ball at t ≈ 20, the amplitude of the oscillations grows significantly as the imploding pulse transfers energy to the configuration. It oscillates continuously around the stationary (unperturbed) solution before slowly returning toward the original configuration.
If configuration L1 was unstable, one would expect that the interaction between φ and χ would excite any unstable modes underlying the solution. Once excited, these modes should quickly grow and bring about the destruction of the configuration. However, no such behaviour is observed in our numerical experiments using different values of c and A. In addition, we have also analyzed the behaviour of the configuration when subject to Type 0 perturbations only, finding no evidence of instability. We therefore conclude that L1 is stable.
Next we consider L2, which lies on the upper branch of the solution curve in Fig. 3. We subject this configuration only to a Type 0 perturbation. The time evolution of L2 is depicted in Fig. 6. The scalar field modulus retains its initial shape for only a short time before quickly decaying and spreading radially. As the evolution proceeds, the fields continue to propagate toward the boundaries until no significant remnant of the initial configuration remains in the domain. As mentioned previously, the timescale over which the Q-ball survives before dissipating can be extended by increasing the numerical resolution of the simulation (thereby decreasing the size of the perturbation). However, even with five additional levels of refinement, the solution begins to decay within the first few oscillations of the scalar field. Since the outcome of this evolution is the total dispersal of the fields, we classify L2 as unstable.
Consider next L3, which lies on the lower branch of Fig. 3, near the bottom axis. Notably, f (0) ≈ 0 near this axis which results in |φ| attaining its maximal value away from the origin. This configuration resembles a shell-like distribution of matter and is therefore labelled a "gauged Q-shell". The evolution of this configuration subject to a Type 0 perturbation is shown in Fig. 7. At the initial time (top panel), the shell-like nature of the solution is readily apparent. As time evolves, the spher- cal for the entirety of the evolution (aside from oscillations and distortions induced by the fragmentation process) and represent smaller "child" gauged Q-balls of the initial configuration. In addition, we observe that the field also fragments into several distinct components which coalesce away from the axis of symmetry. In threedimensions, these resemble ring-like structures which we call "gauged Q-rings". Those Q-rings which are closest to the axis quickly collapse back into spherical structures (child gauged Q-balls) which remain on the axis of symmetry for the rest of the evolution. However, those rings which are initially farthest away from the axis of symmetry are observed to propagate outward and can survive for some time. The bottom panel of Fig. 7 illustrates the behaviour of the gauged Q-rings associated with the decay of L3. The largest Q-ring reaches a maximal distance from the origin of ρ ≈ 40 before turning around and collapsing back onto the axis of symmetry by t ≈ 500. We classify L3 as an unstable configuration. Moreover, we find that all Q-shell solutions on the lower branch of Fig. 3 are unstable. It is notable that this particular branch of solutions was reported to be stable in Ref. [38] under spherical symmetry assumptions. However, the formation of rings is obviously forbidden under spherical symmetry, so our current results are not inconsistent with previous findings.
The formation of gauged Q-rings does not appear to be a unique feature of the decay of L3. We observe a similar phenomenon for other Q-shells on the lower branch of Fig. 3 as well as for the decay of some unstable gauged Q-balls on the upper branch (though the resulting rings may differ in size and lifetime). We have not been able to find any gauged Q-rings which can survive indefinitely. In each case, the rings propagate outward some finite distance from the axis of symmetry before collapsing back inward and forming a gauged Q-ball. This behaviour is similar to what has been observed for non- gauged Q-balls. In Ref. [50], rings are formed through the high-energy collisions of non-gauged Q-balls which also collapse back inward at late times. Q-ring solitons with semitopological origin have also been considered in Ref. [29]. While the rings observed here do not persist indefinitely, they appear to retain their shape despite the relatively violent dynamics that occur after the decay of L3 (until eventually collapsing onto the axis of symmetry). Since these rings could potentially survive long enough to interact with other structures and produce dynamical effects, we conjecture that they may represent a new type of non-spherical solution in the model. Initially, the energy density is negative in a region surrounding the origin and positive elsewhere. As the evolution proceeds, the energy density near the origin grows; the positive region also grows to compensate. Note that the total energy integrates to a positive quantity and is conserved to within 1% over the timescales shown here.
Finally, let us discuss the evolution of L4. This configuration lies on the upper branch of Fig. 3 above L2. This configuration is subject only to the Type 0 perturbation. When evolved forward in time, we observe that the modulus of the scalar field quickly grows without bound until large field gradients are produced. These excessive field gradients cannot be numerically resolved by our code even with increasing adaptive mesh refinement, leading to termination of the evolution due to computational constraints. In Fig. 8, we plot a radial slice of the energy density of the configuration at various points during the evolution. At the initial time, the energy density of the configuration is already negative near the origin. This is likely a consequence of the logarithmic scalar field potential (3) being unbounded from below: when the value of |φ| is large enough, the scalar potential term V (|φ|) in (10) can become negative. If V (|φ|) dominates locally over the other energies in the system, then the energy density at a point in space can also become negative (even while the total integrated energy remains overall positive). When this configuration is evolved forward in time, it may become energetically favourable for the field in the negative region to grow. At the same time, the field in regions of positive energy density would have to grow to compensate and keep the total integrated energy constant. This runaway process results in the large field gradients and unbounded growth (blowup) observed in Fig. 8.
While the decay process of L4 may be unphysical, it is not entirely unexpected. Similar phenomena have been observed in other Q-ball models where local energy densities can attain negative values [51,52]. It is also possible that the decay of L4 could manifest in a different manner (such as dissipation of the fields, similar to L2) if the sign of the initial perturbation to the system could be precisely controlled. However, this level of control is not possible with the Type 0 perturbation, and the fields are found to grow too quickly for Type I perturbations to be effective. In any case, the evolution of L4 results in the destruction of the configuration, and we therefore classify L4 as unstable.
Having discussed the specific configurations L1-L4, we now turn to the general regions of stability and instability depicted in Fig. 3. The black solid line on the bottom branch indicates the regions of the solution curve which are found to be stable under both Type 0 and Type I perturbations. L1 lies in this region. At the leftmost edge of the bottom branch, we observe a turning point where the gauged Q-ball configurations suddenly become unstable. As this turning point is approached from above, the stable gauged Q-ball configurations become less robust: it becomes possible for a sufficiently large Type I perturbation to "kick" the configuration to the unstable branch, though the same configuration remains stable for smaller-sized perturbations. Due to this effect, it is difficult to exactly determine the location of the onset of instability. However, our numerical experiments suggest that it corresponds with the leftmost edge of the lower branch as depicted in the figure. The region of the curve below this point, marked by a red solid line, is found to be unstable. This region contains L3 along with other Q-shell solutions. Lastly, all portions of the curve along the upper branch (containing L2 and L4) are found to be unstable. The solutions which are found to exhibit the blowup behaviour when subject to Type 0 perturbations (including L4) are indicated by a red dashed line along this curve.
To conclude this section, let us summarize the main dynamical behaviours observed in the logarithmic model. For the stable configurations, we find that small perturbations remain bounded and the fields remain relatively close to their initial values without any sign of significant growth or decay. For the unstable configurations, we find the most common outcome to be fragmentation into a small number of "child" Q-balls which quickly propagate away along the axis of symmetry. In some cases (such as L3), this process is accompanied by the formation of Q-rings, while in other cases (such as L2), no significant Q-ball or Q-ring remnants are formed. At present, we have been unable to identify a simple criterion which can predict these changes in behaviour. In general, the fragmentation of gauged Q-balls appears to be a complicated non-linear process, with the only guarantee being the conservation of four-momentum and charge.

B. V6 Model
Here we consider the dynamical stability of gauged Qballs in the polynomial model (4). For illustrative purposes, we select two configurations P1 and P2 whose properties are listed in Table I  First we consider the evolution of P1. This configuration lies on the shortest curve of Fig. 2 with e = 0.17, which is near the maximum allowable value of e ≈ 0.182 [42]. The evolution of P1 is subject to a Type I perturbation with parameters A = 0.1 and c = 0.1 and runs for 65000 base-grid time steps up a final time of t = 6400. The maximal value of the scalar field modulus |φ| and the gauge field component A 0 for this evolution is shown in Fig. 9. The perturbative scalar field hits the Q-ball at t ≈ 20 before exploding outward and exiting the simulation domain. After the collision, the Q-ball is left oscillating at the origin around the stationary (unperturbed) configuration. However, the magnitude of this oscillation rapidly decays as the Q-ball quickly returns close to the original configuration. Similar behaviour to P1 is observed for all other solutions tested on the e = 0.17 branch depicted in Fig. 2. We therefore conclude that P1 (as well as all other solutions tested on this branch) is dynamically stable.
Finally, we consider the evolution of P2. This configuration is distinctive in that it occupies a much larger volume than any of the configurations previously considered. In addition, the scalar field profile of P2 is relatively uniform in the center of the Q-ball before dropping off rapidly to zero at a radial distance r ≈ 65. In this sense it somewhat resembles a Q-ball of the thin-wall type [16]. In Fig. 10, we show the evolution of P2 subject to a Type 0 perturbation. The distinctive flat-top shape of the configuration is apparent in the first panel of the figure. By t ≈ 525 (second panel), the original spherical shape of the configuration is lost as the field content begins to concentrate away from the axis of symmetry. At late times, these off-axis concentrations separate into two distinct Q-rings while a relic region of Q-matter remains near the origin.
In Fig. 11, we plot the growth of the scalar field modulus |φ| for configuration P2. Here, the difference ∆|φ| = |φ(t = 0, ρ, z)| − |φ(t = 225, ρ, z)| illustrates the growth of the solution between the initial time and at a point midway through the evolution (but before the dynamics have entered the non-linear regime). It is clear from the figure that the growth occurs predominantly near the edge of the Q-ball and resembles the pattern of the Y 4,0 spherical harmonic. This pattern becomes apparent in the evolution by t ≈ 100 and grows exponentially in amplitude until the Q-ball begins to break apart starting at t ≈ 500. As mentioned previously, it is well-known that the decay of unstable non-gauged Qballs is always mediated by a spherically-symmetric mode [39]. However, it remains an open question in the literature as to whether gauged Q-balls can be destroyed by the growth of non-spherical modes. Here we have found an example of an unstable gauged Q-ball where the growth of the dominant unstable mode appears to be non-spherical. Remarkably, this occurs even for a small gauge coupling value of e = 0.02. This result is suggestive (though not conclusive) that the destruction of gauged Q-balls can be mediated by non-spherical modes, in contrast to their non-gauged counterparts. However, we emphasize that we have not made any attempt to perform a full stability analysis in this work.
In Fig. 12, we plot the location of P2 on the e = 0.02 curve in the solution space. The curve can be broken down into several distinct branches: an upper unstable branch (I), a stable branch (II), and a lower unstable branch (III) which contains P2. The branch (III) is characterized by solutions which are dominated by a large, nearly-homogeneous central region and have thin surface boundaries, similar to P2. The radial extent of these solutions increases with ω along branch (III). In most cases, the gauged Q-balls on this branch are found to decay slowly into smaller gauged Q-balls or Q-rings, in contrast to the unstable branch (I) where the instability quickly manifests via complete dispersal of the fields to infinity (similar to L2 in the logarithmic model). However, for some configurations along branch (III) which are close to the transition point with branch (II), we also observe the development of large oscillations in the Qball interior which significantly disrupt the shape of the configuration but do not cause the Q-ball to immediately break apart. These oscillations are accompanied by the radiation of charge toward infinity. Since these solutions lose their resemblance to the initial configuration, we also classify them as unstable.
Before concluding, let us comment on the validity of the classical stability criterion (1) for e = 0.02. Since the gauge coupling is very small for this case, one might expect that the regions of stability should be well-predicted by the sign of dQ/dω. Indeed, we find this to be the case. The unstable branches (I) and (III) in Fig. 12 approximately correspond with (ω/Q) dQ/dω > 0 while solutions on the stable branch (II) approximately correspond with (ω/Q) dQ/dω < 0. In this sense the gauged Q-balls with e = 0.02 closely resemble their non-gauged counterparts. In contrast, we have found the entire space of solutions for e = 0.17 to be stable despite the fact that one can find both (ω/Q) dQ/dω > 0 and (ω/Q) dQ/dω < 0 for these solutions. This supports the finding that the sign of the criterion (1) provides little information on the  11. Plot of the difference in the scalar field modulus ∆|φ| = |φ(t = 0, ρ, z)| − |φ(t = 225, ρ, z)| for configuration P2 subject to a Type 0 perturbation. Here, the growth of |φ| occurs predominantly near the edge of the Q-ball and resembles the pattern of the Y4,0 spherical harmonic. The corresponding plot for the full evolution of |φ| is given in Fig. 10. classical stability of gauged Q-balls when the magnitude of the gauge coupling is appreciable [38].

VI. CONCLUSION
We have performed fully non-linear numerical evolutions of U (1) gauged Q-balls in axisymmetry to investigate their stability and dynamics. We assessed this stability in two ways: by perturbing the configurations using the inherent truncation error of the numerical grid, and by introducing an auxiliary massless real scalar field which acts as a perturbing agent designed to explicitly excite any underlying unstable modes. Our simulations suggest that both stable and unstable gauged Q-ball configurations can exist in both the logarithmic and polynomial models. For those solutions which are classified as stable, we observe no evidence of growing modes on the timescales of our evolutions. These solutions respond to perturbations by oscillating continuously or weakly radiating before returning to the original configuration. On the other hand, the decay of unstable configurations can manifest in several different ways, such as total dispersal of the solution, fragmentation into smaller gauged Qballs, or via the formation of non-spherical ring-like structures which we call "gauged Q-rings". These structures appear to be similar in appearance and behaviour to the rings observed in previous studies of non-gauged Q-ball dynamics [50]. Additionally, for some solutions governed by the logarithmic potential and which attain large field values, we observe that the configurations respond to perturbations by growing without bound. This is similar to behaviour observed in other Q-ball models that permit a negative energy density and is interpreted as a consequence of the scalar potential being unbounded from below. For the polynomial potential, we have also investigated the dynamical behaviour of gauged Q-balls when the gauge coupling is small. In this case, we find that the regions of stability and instability are well-described by the stability criterion (1).
One expected result from our study is that those configurations which were known to be unstable with respect to spherically-symmetric perturbations [38] are also unstable under axisymmetric ones. However, our results indicate that axisymmetric perturbations can lead to new regions of instability in the solution space. Furthermore, we have found that some unstable gauged Q-ball configurations can be destroyed by the growth of modes which appear to be non-spherical. These results suggest that the decay of some gauged Q-ball configurations may be mediated by non-spherical modes, in contrast to nongauged Q-balls.
While we have presented numerical evidence that gauged Q-balls can be classically stable with respect to axisymmetric perturbations, it is possible that more general perturbations may eventually destroy any gauged Qball. Addressing this issue might be accomplished with a fully three-dimensional code. Future work may also focus on trying to explicitly solve for the non-spherical gauged Q-ring configurations and numerically evolving them in order to assess their stability. Another interesting question relates to the interaction of two stable gauged Q-balls, which will be the subject of a future paper.
Lastly, we would like to emphasize that our results have addressed only the classical stability of gauged Qballs (i.e., stability of the solutions with respect to small perturbations of the fields). For a complete picture of Qball behaviour, one may also wish to consider quantum effects. For example, a Q-ball may decay through collective tunnelling or by surface evaporation when coupled to additional fields [53][54][55][56]. It is interesting to ask whether the stability of gauged Q-balls is maintained once these effects are considered, though such a question is beyond the scope of the present work. our code. Let us define the convergence factor Q c (t) as where h represents the spacing between points on the numerical grid, u n represents the solution computed with grid spacing n, and · denotes the L 2 -norm. For a second-order-accurate finite-difference scheme, it can be shown that Q c (t) → 4 as h → 0. In Fig. 13, we plot the convergence factor resulting from our test for the constraint equations and several representative fields. The rate of convergence is found to be approximately second-order (corresponding to Q c (t) = 4) which is to be expected for our second-order Crank-Nicolson finitedifference implementation. As a secondary measure, we have performed an independent residual test [57] to verify that our discrete numerical solution is converging to the true continuum solution of the underlying system (A1)-(A6). For this test, we substitute the numerical solution found via the second-order Crank-Nicolson discretization into a separate first-order forward discretization of the evolution equations. Results of this test are shown in Fig. 14. The residuals are found to approximately overlap when rescaled by factors of 2 n , indicating the expected firstorder convergence. For brevity, only the (A1) residual is presented here -other residuals are found to be similar.
In solving the equations of motion, we have used a free evolution scheme wherein a solution to the evolution equations is expected to solve the constraint equations at the initial time [58]. However, it is possible for constraint violation to grow during the course of the evolution, indicating lack of convergence. The degree to which the constraints are violated is therefore a relative measure of the error in the numerical solution. In all of our simulations, we monitor the L 2 -norm of the constraint residuals to ensure that they do not grow significantly over the timescales explored. For the results reported here, the L 2 -norm of the constraint residuals remains within O(10 −4 ). In addition, we also monitor the integrated total energy E and the charge Q during the course of the evolution to confirm that these quantities remain approximately conserved to within O(1%).