text:
# Q-ball Dynamics

## 1 Introduction

One of the most fascinating areas of inter-disciplinary research at the interface between mathematics and physics is the study of solitons. This word has as many definitions as there are people who study them, but in general terms they are stable, localized energy distributions. From a purely mathematical perspective, solitons are described as extended solutions to a set of hyperbolic or parabolic partial differential equations, which can travel without dissipation at a uniform velocity and maintain, at least asymptotically, their shape during collisions; often these properties of solitons are attributable to the existence of an infinite number of conserved quantities, connected with the notions integrability, and radiation-free soliton collisions can be constructed.

In the context of particle physics, which is our main interest here, the usage of the word soliton is less rigorous and any kind of localized energy distribution falls under this broad umbrella. Only in very specialised circumstances will soliton collisions not generate radiation and solutions which radiate substantially during, for example, a highly relativistic 2-soliton collision are included. Of course very little exact analytic progress is possible in these more general settings since radiative processes are notoriously difficult to model analytically, thus numerical simulations are necessary to probe more complicated situations. The usual development of understanding in this subject follows an intricate interplay between the two, with analytic work putting on solid ground more qualitative observations from simulations. The domain of validity of analytic approximations can then be checked and extended by further simulations. Examples of the classes of solitons which have been investigated in this way are vortices , monopoles  and skyrmions , whose existence and stability is essentially due to conserved topological currents and charges, along with energy bounds related to the charge and stable scaling laws. In these examples the topological features constrain the amount of radiation produced in low energy collisions and allows approximations to be applied which ignore the radiative effects.

The subject of this paper is a particular class of solitons known as Q-balls . These are different in many ways to the topological solitons mentioned above. Firstly, they are time-dependent with a rotating internal phase. Secondly, the conserved charge associated with their stability, Noether charge (Q), is not topological and therefore their stability is also a dynamical issue. As we will see these two features lead to a much more complicated variety of interaction properties than seen in the study of topological solitons. The main difference is that the charge is quantized in topological models, it usually being scaled to be an exact integer, whereas we shall see that the charge Q can take any value (in a specified range) allowing for the possibility of charge transfer between solitons and/or fission during the interaction process.

Although the concepts associated with Q-balls are extremely general and they are likely to exist in a wide variety of physical contexts (for example, see ref. ), the main motivation for the current study is the recent realization that they are a generic consequence of the Minimal Supersymmetric Standard Model (MSSM)  due to the existence of D-flat directions in the effective potential created by tri-linear couplings. In this context the conserved charge is that associated with the U(1) symmetries of Baryon and Lepton number conservation and the relevant U(1) fields correspond to either squark or slepton particles. Therefore, the Q-balls can be thought of as condensates of either a large number of squarks or sleptons. It has been suggested that such condensates can be involved in baryogenesis via the Affleck-Dine mechanism  after an epoch of inflation in the early universe. If this is the case then there are two interesting possibilities. If the Q-balls can avoid evaporation into lighter, stable particles such as protons, then it might be possible for them to be important cosmologically as cold dark matter . Whereas if they are unstable, they would decay in a non-trivial way into baryons and could lead to observable isocurvature baryon fluctuations .

Underpinning these interesting suggestions are assumptions as to how Q-balls actually interact and it is our intention here to make an exhaustive study of this issue. Our approach will be to identify numerically the important dynamical processes that can occur in general situations of two interacting Q-balls, which we will then explain qualitatively, leaving a more detailed analytic exposition of the dynamics  and the cosmological implications to subsequent papers. In the next section we will discuss the basic properties of static U(1) Q-balls. Then we will present a detailed and extensive study of Q-balls on the line. As we will see, many of the properties of interest such as fission and charge transfer are observed in one-dimension and given the simplicity of simulations, it seems sensible to make the most exhaustive study there. In the subsequent sections on planar Q-balls and fully three-dimensional Q-balls we will show to what extent the one-dimensional simulations can be used to understand the dynamics in higher dimensions and what effects are clearly of higher dimensional origin. In a penultimate section we will discuss the interactions of Q-balls with anti-Q-balls which have an equal and opposite rotation, before making a concluding summary in the final section.

We should note that there is a disparate literature  on Q-balls in which some (but by no means all) of the processes we will discuss have already been noted, but not necessarily completely understood. In particular we should mention recent work  which presented results for Q-balls in one and two dimensions. There it was suggested that the right-angle scattering of solitons seen in two-dimensional topological soliton models is also prevalent in these non-topological models. At the relevant points we will point out in what ways we disagree with their explanation of this phenomena, and demonstrate that it is by no means general.

## 2 U(1) Q-balls

Given our motivation for studying Q-balls it seems sensible to work with the U(1) Goldstone model, although Q-balls can exist in a variety of field theoretic models. To be precise, the model we consider is that of a single complex scalar field $`\varphi `$ in $`D=1,2,3`$ spatial dimensions with a potential $`U(|\varphi |).`$ Explicitly, the Lagrangian is

$$=\frac{1}{2}_\mu \varphi ^\mu \overline{\varphi }U(|\varphi |),$$
(2.1)
with the key feature being the fact that the potential is only a function of $`|\varphi |`$. The model has a global $`U(1)`$ symmetry and the associated conserved Noether current $`J_\mu `$ is given by

$$J_\mu =\frac{1}{2i}\left(\overline{\varphi }_\mu \varphi \varphi _\mu \overline{\varphi }\right),$$
(2.2)
whose covariant conservation $`^\mu J_\mu =0`$ leads to the existence of the conserved Noether charge Q, given by

$$Q=\frac{1}{2i}\left(\overline{\varphi }_t\varphi \varphi _t\overline{\varphi }\right)d^Dx=\text{Im}(\overline{\varphi }_t\varphi )d^Dx.$$
(2.3)
A stationary Q-ball solution has the form

$$\varphi =e^{i\omega t}f(r),$$
(2.4)
where $`f(r)`$ is a real radial profile function which satisfies the ordinary differential equation

$$\frac{d^2f}{dr^2}=\frac{(1D)}{r}\frac{df}{dr}\omega ^2f+U^{}(f),$$
(2.5)
with the boundary conditions that $`f(\mathrm{})=0`$ and $`\frac{df}{dr}(0)=0.`$

This equation can either be interpreted as describing the motion of a point particle moving in a potential with friction , or in terms of Euclidean bounce solutions ; in each case the effective potential being $`U_{\mathrm{eff}}(f)=\omega ^2f^2/2U(f)`$. This leads to constraints on the potential $`U(f)`$ and the frequency $`\omega `$ in order for a Q-ball solution to exist. Firstly, the effective mass of $`f`$ must be negative. If we consider a potential $`U(f)`$ which is non-negative and satisfies $`U(0)=U^{}(0)=0`$, $`U^{\prime \prime }(0)=\omega _+^2>0`$, then one can deduce that $`\omega <\omega _+`$. Furthermore, the minimum of $`U(f)/f^2`$ must be attained at some positive value of $`f`$, say $`0<f_0<\mathrm{},`$ and existence of the solution requires that $`\omega >\omega _{}`$ where

$$\omega _{}^2=2U(f_0)/f_0^2.$$
(2.6)
Hence, Q-ball solutions exist for all $`\omega `$ in the range $`\omega _{}<|\omega |<\omega _+`$. Note (i) that solutions exist for positive and negative values of $`\omega `$, the negative ones being termed anti-Q-balls, (ii) it is often interesting to think of the Q-balls as being akin to charged bubbles; their profiles being very similar.

The classical stability of the solutions is a more sensitive issue. For sufficiently large $`Q`$ these solutions are guaranteed to be stable, as can be seen using the ‘thin wall limit’ , where the profile function can be approximated by a smoothed-out step function. For small Q it is necessary to perform a full stability analysis using the second variation of the action. In general, the results depend on the details of the potential, but it can be shown that arbitrarily small Q-balls are stable for certain potentials . For a rigorous approach to the classical stability of Q-balls see ref.  and references therein. From the quantum mechanical point of view, the solutions are always stable for large enough $`Q`$ since the energy per unit charge approaches $`\omega _{}`$, which is always less than that for the $`\varphi `$ particle itself, which is $`\omega _+`$.

In choosing a simple potential which admits Q-balls there are three natural classes which have been considered, although there are obviously many other possibilities,

$`I`$ $`:U(f)=\alpha _1f^2+\alpha _2f^4+\alpha _3f^6,`$ (2.7)
$`II`$ $`:U(f)=\beta _1f^2+\beta _2f^3+\beta _3f^4,`$ (2.8)
$`III`$ $`:U(f)=\gamma _1f^2(1\gamma _2\mathrm{log}(\gamma _3f^2))+\gamma _4f^{2p}.`$ (2.9)
In each case it is possible to remove two of the parameters by rescaling the units of energy and time. Note therefore that the potentials of type I and II have one free parameter, while potentials of type III have two free parameters for a fixed value of $`p`$.

The type I potential is the simplest allowed potential which is a polynomial in $`f^2`$, while type II is the simplest which is polynomial in $`f`$. Finally, those of type III mimic the D-flat direction in the MSSM. Here $`p6`$ is some integer that ensures the growth of the potential for large $`f`$, but does not destroy the flatness property for intermediate values of $`f.`$ We should note that none of these types of potential are the kind which might be associated with a renormalizable quantum field theory, but are typical of effective theories incorporating radiative or finite temperature corrections to a bare potential.

In this paper we will mainly be concerned with the type I potential, although we have also studied the type II case. We should note that although the existence and stability properties of Q-balls with these potentials are somewhat different, the qualitative features of the dynamics appears to be almost independent of the potential. The reason for this is that the main interaction processes that we will describe, charge transfer and fission, are due mainly to the time-dependent nature of the solution, rather than the precise profile function.

Explicitly, we shall choose our potential so that

$$U(f)=f^2(1+(1f^2)^2),$$
(2.10)
and therefore in terms of the earlier notation we have that $`w_+=2`$ and $`w_{}=\sqrt{2}`$, so that stable Q-balls exist for $`\sqrt{2}<\omega <2`$. To illustrate the important features of Q-ball solutions we shall focus on the case of one dimension where the profile function equation (2.5) can be solved exactly  to give

$$f_\omega (r)=\sqrt{\frac{4\omega ^2}{2+\sqrt{2\omega ^24}\mathrm{cosh}(2r\sqrt{4\omega ^2})}}.$$
(2.11)
The associated energy $`E_\omega `$ and charge $`Q_\omega `$ can then be computed to be

$$Q_\omega =\sqrt{2}\omega \mathrm{tanh}^1\left(\frac{2\sqrt{2\omega ^24}}{\sqrt{2}\sqrt{4\omega ^2}}\right),E_\omega =\frac{\sqrt{4\omega ^2}}{2}+\frac{(\omega ^2+2)}{2\omega }Q_\omega .$$
(2.12)
In figure 1 we plot the charge $`Q_\omega `$ and energy $`E_\omega `$ for $`\omega `$ in the allowed range $`\sqrt{2}<\omega <2.`$ and we see that both $`Q_\omega `$ and $`E_\omega `$ are monotonically decreasing functions of $`\omega .`$ From (2.11) we can deduce that

$$f_\omega (0)=\sqrt{(4\omega ^2)/(2+\sqrt{2\omega ^24)}},$$
(2.13)
and therefore $`f_\omega (0)`$ increases with the charge $`Q_\omega `$, since it is a decreasing function of $`\omega .`$ In figure 2 we display the energy per unit charge $`E_\omega /Q_\omega `$ as a function of the charge $`Q_\omega `$, from which it can be seen that $`E_\omega /Q_\omega `$ is a decreasing function of the charge. Recall that the asymptotic limit as $`Q_\omega \mathrm{}`$ is $`E_\omega /Q_\omega =\omega _{}=\sqrt{2}.`$ Thus, these Q-balls are stable against decay into a number of smaller Q-balls preserving the total charge.

Most of the general properties of Q-balls in any dimension and for differing choices of the potential are captured by this one-dimensional example, where explicit formulae are available. However, there are some slight differences, for example if $`D=3`$ then there is a lower bound on the charge of a Q-ball, whereas in the $`D=1`$ case considered above Q-balls can have an arbitrarily small charge. This is not a generic feature of every potential and, for example, it has been shown that for $`D=3`$ arbitrarily small Q-balls can be found with a potential of type II using a ‘thick wall limit’ in ref. . These kind of details are easily determined by solving the profile function equation (2.5) numerically and a complete treatment of these issues can be found in ref. . But, as we have already noted, we don’t believe that they are important for the dynamical processes which we focus on in the subsequent sections.

## 3 Q-ball dynamics on the line

In this section we shall investigate the dynamics of Q-balls in one-dimension; even for $`D=1`$ we shall see that multi-Q-ball dynamics is a complicated issue, there being a rich variety of phenomena associated with the non-quantization of the charge and the time-dependent phase.

To investigate the dynamics of Q-balls we numerically solve the field equations which follow from the Lagrangian (2.1) with the potential (2.10), namely

$$\ddot{\varphi }^2\varphi +2\varphi (24|\varphi |^2+3|\varphi |^4)=0,$$
(3.1)
which is valid for any value of D. The numerical methods we use are simple finite difference schemes involving either second or fourth order accurate spatial derivatives and a second order leapfrog algorithm for the time evolution with 1000 points, the spatial step size $`\mathrm{\Delta }x=0.1`$ and the time step size $`\mathrm{\Delta }t=0.02`$. We apply absorbing boundary conditions, which allows radiation to leave the grid and therefore simulates an infinite domain (see refs.  for details on how to apply these boundary conditions).

As initial conditions to describe two well-separated Q-balls we use the ansatz

$$\varphi =e^{i\omega _1t+i\alpha }f_{\omega _1}(|x+a|)+e^{i\omega _2t}f_{\omega _2}(|xa|),$$
(3.2)
in one dimension, which can be trivially generalized to higher dimensions. This ansatz describes a Q-ball with frequency $`\omega _1`$ at the position $`x=a`$ and a second Q-ball with frequency $`\omega _2`$ at the position $`x=a.`$ The $`U(1)`$ symmetry of the theory means that for a single Q-ball the phase of $`\varphi `$ can be set to zero at $`t=0`$ without loss of generality. However, for multi-Q-ball configurations only the initial overall phase can be removed and so for a 2-Q-ball configuration the relative phase, $`\alpha `$, remains as an important parameter.

The total charge of this configuration is

$$Q=Q_{\omega _1}+Q_{\omega _2}+(\omega _1+\omega _2)\mathrm{cos}\alpha _{\mathrm{}}^{\mathrm{}}f_{\omega _1}(|x+a|)f_{\omega _2}(|xa|)𝑑x.$$
(3.3)
The final term in the above expression is exponentially small in the separation parameter $`a`$, since the profile functions have an exponential fall-off. However, we see that the relative phase $`\alpha `$ does affect the value of the total charge. Thus, strictly speaking, it is not valid to substitute this ansatz into the energy functional to determine how the energy depends on the relative phase, since this would involve comparing configurations with differing values of $`Q.`$ The same remark also applies to any attempt to determine how the potential energy depends upon the relative separation $`a,`$ which would help to determine the nature of the interaction force between Q-balls. This issue and its resolution using the methods of ref.  will be discussed in ref. .

To begin with we consider configurations in which the two Q-balls have the same charge, that is $`\omega _1=\omega _2=\omega .`$ In figure 3, and in all subsequent figures illustrating the dynamics of Q-balls, we plot the charge density (the integrand in equation (2.3)) for the initial conditions and at later times111We could have made similar plots of the energy density which are not equivalent. We believe the charge density gives a better representation of the dynamics.. The parameter values used for this simulation are $`\omega =1.5`$ and $`a=4`$, with the Q-balls initially in phase so that $`\alpha =0.`$ The two Q-balls slowly attract and coalesce to form one larger Q-ball which has a charge which is slightly less than the sum of the charges of the two original Q-balls; the charge deficit being carried away by the fission of two additional Q-balls which, in figure 3, can just be seen moving away from the almost stationary large Q-ball at the origin. The attraction of the two Q-balls is simple to explain; it being a consequence of the ratio $`E/Q`$ decreasing as $`Q`$ increases. However, the process of fission is less intuitive and is a novel concept to those who might have studied the dynamics of topological solitons in an attractive potential. In the topological case the charge is an integer and so the fission of higher charge solitons can only be achieved when the solitons are moving sufficiently fast for the kinetic energy to overcome the attraction and release a soliton. But here the charge of an isolated Q-ball can have any value, arbitrarily close to zero, and so one might imagine that the energy barrier to fission at low interaction speeds is small, particularly when the charge is high.

The fission of Q-balls is a process which can occur when a Q-ball suffers a large distortion, for example, during a collision. This can be demonstrated by taking a single Q-ball and squashing it by applying the scale transformation

$$x\lambda x,\varphi \sqrt{\lambda }\varphi ,$$
(3.4)
where $`\lambda >1`$ is a scale factor. The scaling of the $`\varphi `$ field is required in order that the charge of the Q-ball is not changed by the scaling. For small enough values of $`\lambda `$ the Q-ball oscillates, but then eventually returns to the original un-squashed Q-ball corresponding to $`\lambda =1`$, which is to be expected since these Q-balls are stable. However, for a sufficiently large distortion the Q-ball will break up into smaller Q-balls. This is illustrated in figure 4, where we take a Q-ball with charge $`Q=8.4`$ and perturb it with a scale $`\lambda =1.6.`$

To consider the efficiency of the fission process as a function of the charge, we define the quantity

$$\mathrm{\Delta }(Q)=(2E(Q/2)E(Q))/E(Q),$$
(3.5)
where $`E(Q)`$ denotes the energy of a Q-ball with charge $`Q.`$ $`\mathrm{\Delta }(Q)`$ is the fractional increase in energy required to allow a charge $`Q`$ Q-ball to fission into two charge $`Q/2`$ Q-balls. This is a monotonically decreasing function of $`Q`$ with $`\mathrm{\Delta }(\mathrm{})=0`$, with the limit $`Q\mathrm{}`$ being a Bogomolny-like limit in which the energy is proportional to the charge. Thus we expect that the fission of Q-balls is more easily stimulated when the charge is large. To verify this we apply the same distortion factor $`\lambda =1.6`$ as displayed in figure 4, but this time we take a smaller Q-ball with charge $`Q=5.6`$, and the results are displayed in figure 5. The Q-ball performs a breather-like motion in which two structures initially begin to form but then recombine. This motion persists for many cycles with a slowly decreasing amplitude until it eventually settles down to a configuration which is very close to the original Q-ball without a distortion (the configuration at $`t=3000`$ is almost identical to the stationary $`Q=5.6`$ Q-ball).

The above expectations of the fission of Q-balls are confirmed by performing simulations, as in figure 3, with two stationary Q-balls and varying the charge (increasing or decreasing the value of $`\omega `$). The charge of the additional Q-balls produced decreases as the charge of the initial Q-balls is reduced and for small enough Q-balls no fission takes place; rather the sole Q-ball formed oscillates for some time, with a decreasing amplitude.

If the two Q-balls are initially Lorentz boosted toward one another, each with a velocity $`v`$ say, then if $`v`$ is large enough the two Q-balls can be made to pass through each other. In figure 6 we display the charge density for a simulation with $`\omega _1=\omega _2=1.5`$, $`a=10`$ and $`v=0.3.`$

In this case the two Q-balls pass through each other, although they do lose some charge via radiation during the interaction process and their velocities are reduced. For a slightly lower value of the velocity the two Q-balls pass through each other, but do not escape to infinite separation. Rather they subsequently recombine, forming a stationary Q-ball at the origin and producing two additional Q-balls in the same manner as described above. In figure 7 we plot the positions of the two main Q-balls (determined as the location of the maximum of the charge density) as a function of time for the velocity $`v=0.28.`$

To study how the relative phase affects the interaction we consider the same initial configuration used to produce figure 3 except that we set $`\alpha =\pi `$, so that the two Q-balls are exactly out of phase. In this case the resulting evolution is very different and the two Q-balls simply drift apart with no change in their shape or charge, even though the crude arguments based on $`E/Q`$ suggest that they should attract. The effects of changing the overall phase are similar in many ways to the overall isospin rotations possible in 2-skyrmion interactions. There it is possible to make the skyrmions attract or repel by an internal SU(2) rotation about the line joining the two soliton centres. Another way of understanding the relative phase is as a current between two charged bubbles, if the Q-balls are thought of as bubbles. The results of this repulsive interaction channel are displayed in figure 8.

As we have seen, for $`\alpha =0`$ the Q-balls attract and for $`\alpha =\pi `$ they repel. However, for intermediate values of $`\alpha `$ the dynamics is much more complicated. We display the evolution for $`\alpha =\pi /9`$ in figure 9, illustrating a novel process which we shall call charge transfer. The Q-balls initially move very slightly towards each other, but they eventually repel and separate off to infinity. The most interesting aspect is that the charge of the first Q-ball has clearly decreased and that of the second Q-ball has increased, despite the fact that the Q-balls remain two distinct objects with only a very small overlap throughout the time evolution. As they separate the smaller Q-ball moves at a greater speed than the larger one.

In figure 10 we plot the total charge $`Q`$ in the right half of the line, $`x>0,`$ as a function of time for the simulation displayed in figure 9, where $`\alpha =\pi /9`$, and also for the cases $`\alpha =\pi /4`$ and $`\alpha =\pi /2.`$ We see that as the relative phase is decreased the rate at which the charge transfer initially takes place is reduced but the total charge transfered is increased. The in-phase limit $`\alpha =0,`$ where the two Q-balls form a single larger Q-ball, is a smooth limit if we interpret it as a total charge transfer. In the out-of-phase limit $`\alpha =\pi `$, as we have already seen, there is no charge transfer. If $`\alpha <0`$ then the same amount of charge is transferred as in the case of a phase $`\alpha `$, but this time it is the Q-ball in the left half of the interval which increases in charge.

It may seem surprising that there is charge transfer, but this result can be understood, at least qualitatively, by considering a simplified mechanical analogue of the field dynamics associated with two well-separated Q-balls. Consider two equal charge Q-balls, fixed at the positions $`x=\pm a`$, then the ansatz (3.2) may be written in the form

$$\varphi =e^{i\theta _1}f(|x+a|)+e^{i\theta _2}f(|xa|),$$
(3.6)
where $`\theta _1\theta _1(t),\theta _2\theta _2(t)`$ are the time-dependent phases of the two Q-balls and $`f`$ is a profile function. To leading order in the separation $`a`$, corresponding to the limit of large separation, the contribution to the Lagrangian is given by

$$=\frac{1}{2}M(\dot{\theta }_1^2+\dot{\theta }_2^2)ϵ^2\mathrm{cos}(\theta _1\theta _2)4M,$$
(3.7)
where $`M=_{\mathrm{}}^{\mathrm{}}f(|x|)^2𝑑x`$ is treated as a constant moment of inertia and $`ϵ^2=4_{\mathrm{}}^{\mathrm{}}f(|x+a|)f(|xa|)𝑑x`$ is a small interaction coefficient. To derive this Lagrangian we have made the assumption that the profile function is time independent, which obviously has a very limited range of validity as we shall discuss further below; it is nonetheless instructive. The equations of motion which follow from (3.7) are

$$\ddot{\theta }_1+\ddot{\theta }_2=0,\ddot{\theta }_1\ddot{\theta }_2=\frac{2ϵ^2}{M}\mathrm{sin}(\theta _1\theta _2).$$
(3.8)
The first of these equations simply represents the fact that the sum of the rotation frequencies $`\dot{\theta }_1+\dot{\theta }_2`$ is conserved, and the second equation determines the dynamics of the relative phase. There are symmetric solutions, $`\theta _1=\theta _2`$, and $`\theta _1=\theta _2+\pi `$, where the phase difference remains constant, corresponding to the two Q-balls being exactly in-phase or exactly out-of-phase for all time, but for general values of the initial relative phase, $`\alpha =\theta _1(0)\theta _2(0)`$, there will be a non-trivial time dependence. For all $`\alpha (0,\pi )`$ there will be an initial positive acceleration in the relative phase, so that $`\dot{\theta }_1>\dot{\theta }_2`$ for small $`t>0.`$ Thus the first Q-ball will have a higher frequency than the second Q-ball and, since we know that the charge of a Q-ball decreases with increasing frequency, then this corresponds to the charge of the first Q-ball decreasing and the charge of the second Q-ball increasing. This simple analysis also predicts that the initial rate of charge transfer will be greatest for a relative phase $`\alpha =\pi /2`$ and will decrease as $`\alpha `$ decreases. This agrees with the observation we made earlier by an examination of the plots in figure 10 for small times.

However, what we are clearly not able to study with our simple restricted mechanical model is the whole charge transfer process for later times. One reason for this is that we assumed that the profile function $`f`$ was fixed when of course we know that it is highly dependent on the rotation frequency (see equation (2.11)). In particular this dependence constrains the rotation frequencies to satisfy $`\omega _{}<\dot{\theta }_1,\dot{\theta }_2<\omega _+`$ and as either of these limits are approached our simple model breaks down. One might be tempted to improve our simple model to deal with this issue by including the known frequency dependence of the profile function, but it is not obvious how to do this since the profile function depends upon the frequency $`\dot{\theta }`$ so using such an ansatz in the Lagrangian would lead to a Lagrangian for a mechanical system with second order derivatives $`\ddot{\theta }`$ and hence a fourth order equation of motion. Furthermore, we have assumed that the positions of the Q-balls are fixed when in fact the results of the full simulations show that they eventually drift apart. This effect will also serve to cut-off the relative phase dynamics since it will correspond to reducing the $`ϵ^2`$ coefficient in our simple mechanical model.

In summary, we have shown that a simple mechanical model is useful in understanding the qualitative features of the charge transfer process, but a more sophisticated analysis is required to explain the quantitative behaviour found. The analysis of relative phase dynamics in mechanical systems, such as discrete breathers, has been studied in some detail and the phase space trajectories are well understood . These methods can be extended to study the more complicated relative phase dynamics, and hence charge transfer, of Q-balls .

So far we have only considered initial conditions in which the two Q-balls have the same charge. For two Q-balls which have different charges the initial relative phase does not have the same importance as for Q-balls of the same charge, due to the fact that the Q-balls have different frequencies and so the initial relative phase will not be preserved, even with no interaction. This can easily been seen using the simple mechanical system above.

In figure 11 we plot the charge density for the initial conditions $`\omega _1=1.8,\omega _2=1.5,a=3,\alpha =0.`$ It can be seen that the two Q-balls repel and there is virtually no charge transfer since the solitons never get close enough for the charge transfer process to become important. Similarly, if a non-zero relative phase is introduced then virtually no charge transfer takes place, although the rate of separation does vary very slightly. However, if the Q-balls are Lorentz boosted towards each other with a sufficiently large velocity $`v`$ so that they collide and pass through each other, it is possible to induce charge transfer as illustrated in figure 12 for the parameters $`\omega _1=1.8,\omega _2=1.5,a=6,\alpha =0,v=0.2.`$ The amount of charge transferred depends on the value of the relative phase as the Q-balls collide, as can be verified by changing the initial phase. It can be checked that this is equivalent to varying the initial separation, since the time to collision is then altered and hence the relative phase is different by an amount equal to the change in collision time multiplied by the frequency difference. Just using the simple mechanical analogy, one might have naively expected that an initial difference in the rotation speeds would be on a similar footing to an initial phase difference, but this is clearly not the case. It is evident that there is a non-trivial interaction between the relative dynamics of the Q-balls and the charge transfer process.

## 4 Planar Q-balls

The main features of one-dimensional Q-balls which we have described in the previous section, such as charge transfer and the dependence of the interaction force on the relative phase, carry through to the two-dimensional case. We demonstrate this by again performing numerical simulations of the field equations via an equivalent finite difference scheme to the one dimensional case. We find that a grid containing $`200^2`$ points with $`\mathrm{\Delta }x=0.2`$ and $`\mathrm{\Delta }t=0.05`$ gives an accurate representation of the dynamics in this case. In contrast to the one-dimensional case an exact solution is not known for the profile function in two-dimensions, but it is a simple matter to numerically obtain the profile function using a standard shooting method.

One might assume that head-on collisions of Q-balls with a small charge (for example, $`\omega =1.6`$) are equivalent to those in one-dimension with attraction, repulsion and charge transfer taking place as before and indeed this is the case as we will discuss later. The phase-dependent force is, however, a more general concept which we will illustrate in fully two dimensional interactions with a non-zero impact parameter. In particular, we consider the collision of two equal charge Q-balls ($`\omega _1=\omega _2=1.6`$) with a non-zero impact parameter, where the initial positions of the Q-balls are $`(x_1,x_2)=\pm (6,3)`$, each is boosted along the $`x_1`$-axis with a velocity $`v=0.05`$, and the relative phase $`\alpha `$ is set to zero. In figure 13 we plot the charge density for this interaction at $`t=0,104,112,3200.`$

As in the one-dimensional case (where collisions are head-on) the two Q-balls attract and form a single larger Q-ball, although the large Q-ball has some angular momentum due to the fact that the collision was not head-on. Figure 14 displays the results of the same simulation except that the two Q-balls are exactly out-of-phase, that is $`\alpha =\pi ,`$ where the Q-balls clearly repel. Finally, in figure 15 we display the simulation with a relative phase $`\alpha =\pi /4`$, where there is an initial attraction, followed charge transfer and finished off by a repulsion which forces the two Q-balls apart. We conclude, therefore, that while the head-on collisions of small Q-balls in two-dimensions can be thought of as being effectively one-dimensional, the same dynamical processes are active in the case of a non-zero impact parameter.

Next, we turn our attention to head-on collisions of Q-balls with higher charge, where — based on the intuition of the one-dimensional interactions — one would expect things to be slightly different. We take two Q-balls with $`\omega _1=\omega _2=1.5`$ at positions $`(x_1,x_2)=\pm (10,0)`$ with each Lorentz boosted towards the other with a velocity $`v=0.4.`$ Figure 16 displays the charge density at $`t=0,20,24,28,32,36,40,44,52`$ for the in-phase case, $`\alpha =0.`$ As can be seen from the figure, there is a very complicated interaction process involving the charge being strongly deformed and the emission of some radiation. Eventually, two Q-balls emerge from the interaction region at right angles to the initial direction of approach. Naively one may think that this is a simple $`90^{}`$ scattering process as seen in a number of topological soliton models  and suggested for Q-balls in ref. 222We should note that the work presented in ref.  uses a potential of type II, not type I as used in our work, but that the qualitative nature of this process is independent of the choice of potential.. However, the scattering of Q-balls is a complicated dynamical issue rather than being topological, and the underlying mechanisms are very different. In particular, there is no associated geometry of a moduli space which forces the Q-balls to scatter at right angles. Rather, during collisions the Q-matter becomes highly deformed with huge charge densities and it is this deformation, and its associated pressures, that lies at the heart of the interaction process and the fission of Q-balls in the plane perpendicular to the incident direction. As we demonstrated in the one-dimensional case, a sufficient distortion of a Q-ball will induce fission, and it this same process which is responsible for this more complicated phenomena in two-dimensions.

This point can be illustrated immediately by considering the same scattering process, but this time we set the Q-balls to be exactly out-of-phase, that is $`\alpha =\pi `$ with the results displayed in figure 17. Although the initial conditions look identical in figures 16 and 17, the evolution is clearly very different. As the two Q-balls are now in a repulsive phase the Q-matter gets distorted in a very different way. Rather than forming a single structure, as in figure 16c, the two individual Q-balls never actually coalesce because of the repulsive interaction, getting squashed separately and this distortion induces the fission of each. Thus, each Q-ball splits into two and the two pairs repel each other, producing four Q-balls in all, which are clearly visible in figure 17g. For this particular set of initial parameters the Q-ball pairs do not have enough energy to escape each others attraction and eventually recombine leaving two Q-balls which move off to infinity. By, for example, increasing the initial velocity it is found that the four Q-balls can be produced in such a way that they all separately move off to infinity without any subsequent recombination. If the initial velocity is small enough then the distortion is not large enough to induce fission and the two Q-balls eventually repel keeping their individual structure intact.

In figure 18 we investigate the same simulation as in figure 16 but with an initial relative phase $`\alpha =\pi /2.`$ In this case we expect that the distortion will also be accompanied by a charge transfer, and indeed this is what we find, as the first Q-ball loses charge to the second one and then each Q-ball undergoes fission to produce a total of two small Q-balls and two large Q-balls. The two small Q-balls move away from the interaction region at a faster rate than the large Q-balls, and they do not recombine. The two large Q-balls move away from the interaction region at only a very slow speed and in fact they eventually recombine after a time of around $`t=200`$ (the final plot in figure 18 is only at $`t=52`$).

As we discussed in the one-dimensional case, large Q-balls are more susceptible to fission than smaller Q-balls, so we expect that the scattering processes we have described above will vary depending on the charge of the initial Q-balls. For example, our reasoning predicts that the fission process in figure 16, which produced two Q-balls moving at right-angles to the initial line of approach, will be more difficult to reproduce for smaller charges. To test this we perform the same simulation, with $`v=0.4`$ again, but decrease the charge of the initial Q-balls by taking $`\omega =1.6`$ rather than $`\omega =1.5.`$ The resulting evolution is shown in figure 19 and confirms that now the distortion is not sufficient to liberate two Q-balls. The configuration oscillates for some time before settling down to a single larger Q-ball, after a small amount of charge has been dissipated through radiation. Fission can be produced for these smaller Q-balls by increasing the impact velocity, but this also has the result that some of the charge passes straight through the interaction region producing small Q-balls which continue to travel along the direction of approach. A collision at increased velocities is also a more violent process and more charge is lost to radiation in these circumstances. In figure 20 we display the evolution for the case where the velocity is increased to $`v=0.6`$, with all other parameters kept the same as in figure 19. Just visible in figure 20d are the four very small Q-balls which are produced by this collision together with a ring of charge carried away by the radiation generated. If the collision velocity is increased further then the two Q-balls have less time to interact and their momentum carries them through the collision process with no deflection. This is demonstrated in figure 21 where $`v=0.8`$ and no additional Q-balls are produced. In summary, as the charge is reduced an increased velocity is required in order for a sufficient deformation to be generated to produce fission, but this also results in more of the charge being carried straight through the collision. Thus for small charge there is a very limited window of velocities for which collisions of the form displayed in figure 16 may occur.

## 5 3D Q-balls

In the previous two sections we have built up a picture of the dynamics of Q-balls in one and two dimensions. In going from the extensive study in one dimension to two dimensions we have noted a number of subtle effects associated with the extra dimension. However, the basic processes involved are the same: attraction, repulsion and charge transfer. In this section we will apply the same numerical techniques to the case of three dimensions. To begin with we conducted an extensive study of the dynamics on grids containing $`100^3`$ points and have once again found that in many cases the dynamics are very similar to those in one-dimension. At the risk of labouring the point we found that for small Q-balls, if they were initially in-phase they coalesced, while if they were out-of-phase they repelled, and if they had any other phase they engaged in charge transfer. However, we did find some extremely complicated interactions which are related to those seen in the case of two dimensions. As was pointed out in the previous section when the Q-balls have a large charge, their interactions can have some interesting variants in two dimensions and it is these particular cases in three dimensions on which we will focus in this section.

In particular we have focussed on the analogues of figures 16, 17 and 18 in which complicated dynamical phenomena were identified. In each of these cases we have performed the analogous simulations on a grid of $`300^3`$ points with $`\mathrm{\Delta }x=0.3`$ and $`\mathrm{\Delta }t=0.03`$ in order to accurately simulate the complicated dynamical processes. In each of the three cases we start with 2 Q-balls each with $`\omega =1.5`$ at $`\pm (15,0,0)`$, Lorentz boosted toward each other with a velocity $`v=0.4`$. The results of the simulations are displayed in figures 22, 23 and 24.

The in-phase case (figure 22) has some marked similarities to the equivalent case in two dimensions (figure 16) if one looks at the two-dimensional slice through the centre of the Q-balls333It should be noted that the two interactions are not equivalent even when the parameters are almost identical since the relationship between the charge $`Q`$ and the frequency $`\omega `$ is not the same in two and three dimensions.. However, the extra dimension has one remarkable effect: it allows for the production of a loop in the plane perpendicular to line joining the two incident Q-balls. This phenomena is the three dimensional analogue of the right-angled fission process described in the two dimensional case. But in three dimensions the fission takes place symmetrically in all directions in the plane while respecting the cylindrical symmetry of the initial configuration. The loop expands leaving some charge in the centre which later expands to create a second, much smaller loop. Later both the loops collapse and Q-balls emerge back along the incident direction.

The formation of the loop in this three dimensional simulation adds further weight to our earlier discussion of the two dimensional case where we pointed out that right-angled fission of two Q-balls was not of topological origin, nor was it even related to the topological interactions of, for example, vortices , monopoles  and skyrmions . The reason being that in a topological interaction in three dimensions one would have expected there to have been a preferred direction (this is because in the case of topological solitons, for example skyrmions, the field configuration of a single soliton is not spherically symmetric, although the change in the field due to a spatial rotation can be undone by acting with a symmetry of the theory, which means quantities such as the energy density are spherically symmetric). But as we have already pointed out the formation of a loop is reliant on all directions being on an equal footing as far as the fission process is concerned, which of course is due to the fact that the field itself is spherically symmetric for a single Q-ball.

This explanation of the formation of a loop during the interaction of two large Q-balls in three dimensions is compatible with the results of the out-of-phase case (figure 23) and that of a general relative phase (figure 24). In both cases, the two dimensional slice through the interaction region is very similar to that of the equivalent two dimensional interactions (figures 17 and 18 respectively). In the out-of-phase interaction, two identical loops form which are then repelled back along the direction from which they came. As they move away they begin to collapse, the final outcome being a series of symmetrically placed Q-balls along the incident direction. The interaction for $`\alpha =\pi /2`$ is similar, except that, as expected, charge transfer takes place during the interaction and the two loops created have very different charge.

Just to finish this section off and by way of illustrating that this process is reasonably generic, we have performed an equivalent simulation to figure 22 with a much higher speed of incidence $`(v=0.8)`$. Due to the Lorentz contraction of the initial conditions, this requires a smaller value of $`\mathrm{\Delta }x=0.15`$ and consequently $`\mathrm{\Delta }t=0.015`$, and the results of this simulation are displayed in figure 25. We see the production of a big loop at the point of interaction plus two others which are repelled from the centre along the line of interaction. At the end of the simulation the loops are still expanding in size and are also getting close to the size of the discrete grid. It is an interesting question as to whether loops can be stable, and this question will be addressed in a separate publication .

## 6 Q-ball Anti-Q-ball Dynamics

In the preceding sections we have studied in detail 2-Q-ball interactions. The charge $`Q`$, however, can also be negative; this being achieved in the Q-ball solution by replacing $`\omega `$ with $`\omega `$ and these solutions are known as anti-Q-balls. In this section we will study the interactions of Q-ball/anti-Q-ball pairs in two and three dimensions.

Intuitively, one would expect slow soliton/anti-soliton interactions with equal and opposite charges to result in annihilation into radiation. However, only for a very small range of parameters does this annihilation take place in the case of Q-balls due to the complicated nature of the time-dependent interaction potential which we have highlighted in the case of 2-Q-ball interactions.

In general a Q-ball/anti-Q-ball interaction will result in either the two solitons bouncing back, or them passing through each other. In both cases, the charge is partially annihilated, but only for a very limited region of the interaction parameter space can it be thought of as being complete. This absence of the annihilation can be attributed to the concept of charge transfer which we have discussed in 2-Q-ball interactions444Here, as in the case of two Q-ball interactions with different charges, the initial relative phase is a less important concept. But charge transfer can still take place since the difference between the two rotation speeds of the Q-balls is maximal. The main difference between the charge transfer process for 2-Q-ball interactions and the Q-ball/anti-Q-ball interactions under consideration here is the two solitons now have opposite charge and hence when charge is transfered it results in annihilation. However, we have showed that the charge transfer never takes place fully in 2-Q-ball interactions and hence annihilation also never takes place fully in a single interaction. When there is a lower bound on the charge of a Q-ball it is more likely that sufficient charge transfer can take place for complete annihilation, whereas in models where arbitrarily small Q-balls can exist complete annihilation is likely to be much more difficult.

The process of partial annihilation, via charge transfer, is illustrated in figure 26 where we have displayed the charge density at $`t=0,15,20,50`$ for the collision of a Q-ball and an anti-Q-ball in two dimensions. The solitons were initially at $`\pm (6,0)`$, with $`\omega _1=\omega _2=1.8`$ and were Lorentz boosted together with a velocity $`v=0.3`$. It is clear that the momentum of the Q-ball carries it through the interaction region and that there is some annihilation of the charge; the maximum charge density of the outgoing Q-ball being lower than that for the incoming one. As we have already discussed this is generically what takes place during a Q-ball/anti-Q-ball collision. A variant on this kind of interaction is that the Q-balls bounce back, again partially annihilating charge, which takes place at low incident velocities for this particular charge.

This picture of partial annihilation with bounce back at low velocities and the solitons passing through each other at high, suggests that there exists some critical incident velocity $`v_\mathrm{c}`$, a function of $`\omega `$ at which annihilation takes place, and indeed this is what we find555In fact, there will exist a small range of velocities around $`v_\mathrm{c}`$ for which complete annihilation takes place.. Figure 27 illustrates this by displaying the charge density at $`t=0,25,50,100,150,250,300,350,400`$ for a Q-ball/anti-Q-ball collision, with the solitons initially positioned at $`\pm (6,0)`$, charges $`\omega _1=\omega _2=1.5`$, Lorentz boosted together with velocity $`v=0.3`$. These are the same parameters as in figure 26, except that the charge is much larger. It can be clearly seen that annihilation eventually takes place, but even in this case the mechanism is complicated, involving a number of oscillations of the system before it is achieved. The charge here is much higher than for the example above which had $`\omega =1.8`$, and in the Q-ball interactions we saw complicated fission processes for interactions involving solitons with this charge. This is also the case in Q-ball/anti-Q-ball interactions as is illustrated in figure 28, which uses the same parameters as in figure 27 except that the incident velocity is now much higher, $`v=0.6`$. The figures shown are at times $`t=0,10,15,25,30,35,40,45,50.`$ One can see, after some partial annihilation of charge, that the fission of the incident Q-balls takes place in the direction perpendicular to the line of incidence and that now there is an effective bounce back of the incident solitons with a reduced charge.

These processes are also prevalent in both one and three dimensions. By way of illustration we have also included two examples in three dimensions. In figure 29 the solitons initially at $`\pm (15,0,0)`$, are Lorentz boosted together with a velocity $`v=0.3`$ and have $`\omega _1=\omega _2=1.5`$. In this particular interaction there appears to be very little annihilation and the solitons effectively pass through each other. Figure 30 has the same initial configuration, except that the solitons are boosted together with an initial velocity of $`v=0.6`$. The subsequent interaction is complicated involving first the formation of two loops, which is the equivalent of the right-angled fission observed in two dimensions, and then what appears to be almost total annihilation once the loops collapse. Two small Q-balls are emitted along the line of incidence, along with much radiation.

## 7 Summary and conclusions

We have identified the key parameters in the interactions of Q-balls to be the relative phase, the incident velocity and the charge of Q-balls, with the resulting interactions being strongly dependent on these parameters. The generic interaction involves attraction if the relative phase $`\alpha =0`$ and repulsion if $`\alpha =\pi `$ when the two Q-balls have the same charge. If they have an initial relative phase other than $`\alpha =0`$ or $`\alpha =\pi `$, or if the charges of the Q-balls are different, the dynamics of the Q-balls results in charge transfer, a phenomena which is analogous to that observed in discrete breather systems (see, for example, ref.  and references therein), although this behaviour often has to be induced by making the solitons come together via a Lorentz boost. With no Lorentz boost such breather systems naturally repel and this is also seen after the charge transfer process is complete. If the incident velocity is extremely high Q-balls can be made to pass through each other with very little interaction taking place.

This picture is almost independent of the number of dimensions if the charge of the Q-balls is small (for example, $`\omega 1.6`$, but if the charge is much larger fission can take place during the interaction process due to the compression of the charge. In one dimension this process results in the emission of Q-balls during the slow interaction of two large Q-balls. In higher dimensions complicated, but analogous, phenomena are observed in which fission takes place in the direction perpendicular to the line of incidence. This leads, under specialized circumstances to the right angled fission of Q-balls in two dimensions and the production of loops in three. These fission effects can also be coupled with those of attraction, repulsion and charge transfer to produce complicated compound phenomena.

Interestingly, the naive expectation that an Q-ball/anti-Q-ball should annhilate into radiation is modified by the complicated breather-type interactions which take place. Since the difference between the two charges is large this can be thought of a special case of the charge transfer process, leading to the phenomenon of partial charge annihilation, the case of complete annihilation being very special. Fission of the Q-ball and anti-Q-ball can also take place if the charge is high.

Our original motivation was to understand the microphysical interactions of Q-balls formed by the Affleck-Dine mechanism for baryogenesis within the framework of the MSSM. The potential that we have concentrated on in this paper is different to that expected in the MSSM, but as we have pointed out the main interaction processes that we have identified are related to the phase. Therefore, we expect our result to be qualitatively independent of potential, and hence our results have some bearing on this case. We have repeated several of the simulations described in this paper for a potential of type II and found the same qualitative results. The next step in this research is an analytic description of the dynamics in terms of a slow manifold approach , before their application to the problem of Q-ball formation in the Early Universe.


format: text

text:
# 1 Introduction

## 1 Introduction

Long time ago, Dolgov and Zakharov  have introduced an alternative approach to the axial anomaly , based on the study of the triangular diagram, with two vector and one axial vertices, by means of dispersion relations. From this approach, the interpretation of the axial anomaly as an infrared phenomenon follows. The infrared aspect of the axial anomaly, arised in this paper, can be seen as complementary to the more familiar ultraviolet one, which emerges from the renormalization procedure. It allows to shed light upon the physical meaning of the chiral symmetry breaking as related to a non conservation of helicity.
In  the absorbitive part of the triangular diagram is computed by making an unitary cut and results proportional to the product of two amplitudes. The first one is relative to the production of a fermion-antifermion pair by an axial source and the second one corresponds to the subsequent annihilation of the pair into two real photons. In both these processes there occur helicity flips, thus the chirality is not conserved in the zero mass limit. For what follows, it is important to examine in particular the second process contributing to the absorbitive part of the triangular diagram ().
By taking the $`m0`$ limit of the absorbitive part of the triangular diagram, Dolgov and Zakharov find a finite result :

$$\mathrm{Im}g_1(q^2)\pi \delta (q^2)\mathrm{as}m0,$$
(1)
with $`g_1(q^2)`$ an invariant scalar function and $`q^2`$ the momentum transfer. Eq. (1) shows that the $`m0`$ limit is not smooth .
This behaviour is connected with a non conservation of helicity, in the massless limit, a non conservation of chirality. Thus the axial anomaly emerges by studying the properties of the amplitude in the infrared region. The purpose of this work is to analyse the cancellation of mass singularities in polarized amplitudes. We will not discuss the physical implications of the zero fermion mass limit.
As stated in refs. , the finite result is a consequence of the singularity occurring in the fermion propagator as $`m0`$, that exactly cancels the suppression factor $`m^2`$ of the product of the amplitudes.
The connection between the anomaly and the non conservation of physical quantities has riceved a lot of attention in the literature. In a seminal work , Gribov has described the source of the anomalies as a collective motion of particles with arbitrarily large momenta in the vacuum. Related to this work is ref. . The Dolgov and Zakharov approach has been also considered in refs. .
We attempt here to interpret some simple processes, as the $`\pi ^+`$ and $`Z^0`$ polarized radiative decays, as a manifestation of the infrared nature of the axial anomaly . In particular, we calculate the relative decay rates, for the cases in which the outgoing leptons are in a definite helicity state and we examine in some detail the structure of the mass singularities and their cancellation. We study how the Kinoshita-Lee-Nauenberg (KLN) theorem  applies to these cases and we consider the analogies and the differences with respect to the corresponding unpolarized decay rates.

## 2 Polarized decay of the pion

In the second process contributing to the absorbitive part of the triangular diagram a charged fermion with fixed helicity undergoes an helicity flip by emitting a photon . This is a forbidden process in the massless limit, but in this limit the absorbitive part doesn’t vanish.
We consider helicity changing reactions, studied for the first time by Lee and Nauenberg , by using the approach of Dolgov and Zakharov.
We have found in the literature other works dealing with this kind of processes. For example in ref.  the authors consider the process

$$e^{}(p_{},\lambda )+A(p)e^{}(p_{}^{},\lambda ^{})+\gamma (k,\lambda _\gamma )+B(q_i),$$
(2)
where $`A`$ is a target, $`\gamma `$ is a bremsstrahlung photon, assumed almost collinear with respect to the direction of the incident electron and $`B`$ is a set of particles produced in the reaction.
For $`\lambda \lambda ^{}`$, the electron makes an helicity flip. The helicity flip cross section, computed to the leading order in the electron mass, is given by:

$$\frac{\mathrm{d}\sigma _{hf}}{\mathrm{d}x}=\sigma _0\left(s(1x)\right)\frac{\alpha }{2\pi }x,$$
(3)
where $`s=(p_{}+p_{}^{})^2`$, $`x=k_0/E`$, $`E`$ is the energy of the incoming electron and $`\sigma _0`$ is the cross section for the Born process

$$e^{}(p_{}k,\lambda )+A(p)e^{}(p_{}^{},\lambda )+B(q_i).$$
(4)
We see that the expression (3) does not vanish in the massless limit .
We now consider the leptonic decay of the pion. At the Born level, the total decay rate relative to the non radiative pion decay

$$\pi ^+l^++\nu _l,$$

where $`l^+`$ is an antilepton ($`e^+`$ or $`\mu ^+`$) and $`\nu _l`$ is the associated neutrino, is given by:

$$\mathrm{\Gamma }_0(\pi ^+l^+\nu _l)=\frac{G^2f_\pi ^2}{8\pi }V_{ud}^2\frac{m_l^2}{m_\pi ^3}(m_\pi ^2m_l^2)^2,$$
(5)
where $`G`$ is the Fermi coupling constant, $`f_\pi `$ is the pion decay constant, $`V_{ud}`$ is the CKM matrix element, $`m_l`$ and $`m_\pi `$ are the lepton and pion masses, respectively.
$`\mathrm{\Gamma }_0`$ is proportional to $`m_l^2`$; this factor is due to the fact, that, due to angular momentum conservation, the pion produces a left-handed lepton, while the structure of the weak coupling requires the $`l^+`$ to be right-handed for $`m_l=0`$.

The total decay rate for the process in which the lepton is polarized can be written as:

$`\mathrm{\Gamma }_0^{pol}(\pi ^+l^+\nu _l)={\displaystyle \frac{G^2f_\pi ^2}{16\pi }}V_{ud}^2{\displaystyle \frac{m_l^2}{m_\pi ^3}}(m_\pi ^2m_l^2)^2\left(1{\displaystyle \frac{𝐩_l𝐬_l}{|𝐩_l|}}\right),`$ (6)
here $`𝐩_l`$ and $`𝐬_𝐥`$ are the linear momentum and the spin vector of the lepton, respectively, giving $`𝐩_l𝐬_l/|𝐩_l|=\pm 1`$, for right-handed and left-handed lepton, respectively. Eq. (6) indicates that the lepton is mandatory left-handed, as requested by angular momentum conservation.
Let us consider the radiative dacay

$$\pi ^+l^++\nu _l+\gamma .$$

The relative amplitude can be written as the sum of two contributions

$$M(\pi ^+l^+\nu _l\gamma )=M_{IB}+M_{SD},$$
(7)
where $`M_{IB}`$ is the Inner Bremsstrahlung amplitude and $`M_{SD}`$ is the Structure Dependent one.
We calculate the probability that the lepton flips its helicity and becomes right-handed, by emitting a real photon. The relevant part, for the problem we are considering, is the Inner Bremsstrahlung contribution described by the diagrams of fig. 1. The $`IB_3`$ diagram gives the so called contact term, introduced to ensure gauge invariance (see for example ). For the moment, we consider only tree diagrams. At the first order in perturbation theory, we retain terms of all powers in the lepton mass. We will argue that the manifestation of the axial anomaly is strictly connected to these terms.
The Inner Bremsstrahlung amplitude is :

$`M_{IB}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}IB_i`$ (8)
$`=`$ $`ie{\displaystyle \frac{G}{\sqrt{2}}}f_\pi V_{ud}m_l\overline{u}(p_\nu )(1+\gamma _5)\left[{\displaystyle \frac{pϵ}{pk}}{\displaystyle \frac{\overline{)}k\overline{)}ϵ+2p_lϵ}{2p_lk}}\right]v(p_l,s_{lR}),`$
where $`u`$ and $`v`$ are the Dirac spinors for the neutrino and the lepton, $`p`$ is the pion momentum, $`p_l`$ and $`s_{lR}`$ are the momentum and the polarization vector of the right-handed lepton, $`k`$ and $`ϵ`$ are the momentum and the polarization vector of the photon, respectively.
We see that $`M_{IB}`$ is proportional to the lepton mass $`m_l`$, thus the decay rate is proportional to $`m_l^2`$. As we have said above, this factor is a consequence of the structure of the weak current. We remove it by normalizing the radiative decay rate with respect to the non radiative one.
The differential Inner Bremsstrahlung contribution for right-handed lepton is:

$`{\displaystyle \frac{1}{\mathrm{\Gamma }_0}}{\displaystyle \frac{\mathrm{d}\mathrm{\Gamma }_{IB}^R}{\mathrm{d}y}}={\displaystyle \frac{\alpha }{4\pi }}{\displaystyle \frac{1}{(1r)^2}}{\displaystyle \frac{1}{A(1y+r)}}\{2A[1+y^22A+r(2A+r6)]`$
$`+\left[(A+2r)(1+r)^2+Ay(y4r)+2ry(1+y)y(1+y^2+5r^2)\right]\mathrm{ln}{\displaystyle \frac{y+A}{yA}}`$
$`+(1y+r)^2(y2rA)\mathrm{ln}{\displaystyle \frac{y+A2}{yA2}}\}.`$ (9)
The dimensionless variable $`y`$ is defined as

$$y=\frac{2E_l}{m_\pi }$$

and

$$r=\frac{m_l^2}{m_\pi ^2},A=\sqrt{y^24r}.$$

The physical region for $`y`$ is:

$$2\sqrt{r}y1+r.$$
(10)
We see that, in eq. (9), there is a term independent of the lepton mass, the one in the first square brackets. Owing to this term, the differential decay rate doesn’t vanish in the limit $`m_l0`$. We have:

$$\frac{1}{\mathrm{\Gamma }_0}\frac{\mathrm{d}\mathrm{\Gamma }_{IB}^R}{\mathrm{d}y}\frac{\alpha }{2\pi }(1y)\mathrm{as}r0.$$
(11)
This indicates that there occurs an helicity flip and thus a chirality non conservation in the limit of zero lepton mass. According to the interpretation given above, this term can be interpreted as connected to the axial anomaly. It corresponds to the anomalous term present in the divergence of the axial current.
Since the radiative process with the right-handed $`l^+`$ is forbidden in the limit $`m_l0`$ by the chiral invariance of the massless QED Lagrangian, the result different from zero, in this limit, indicates the occurrence of a cancellation mechanism, analogous to the one acting in the absorbitive part of the triangular diagram. Indeed, to obtain the decay rate in eq. (9), we have integrated over the emission angle terms containing the lepton propagator squared; this produces power collinear singularities, exactly cancelling the chiral suppression factor $`m_l^2`$ and thus producing the term related to the axial anomaly. The integration over the lepton propagator produces logarithmic collinear singularities.

## 3 Mass singularities

In this section we discuss the structure of mass singularities and their cancellation in the decay rate for the radiative pion and $`Z^0`$ decays and how the KLN theorem applies to these cases.

### 3.1 The pion case

It is useful to separate the cases of unpolarized, right-handed and left-handed outgoing lepton.
Let us consider the mass singularity cancellation for the familiar case of unpolarized radiative $`\pi ^+`$ decay to the first order in $`\alpha `$. The expression for the Inner Bremsstrahlung contribution, differential with respect to the lepton energy, is given by:

$`{\displaystyle \frac{1}{\mathrm{\Gamma }_0}}{\displaystyle \frac{\mathrm{d}\mathrm{\Gamma }_{IB}}{\mathrm{d}y}}`$ $`=`$ $`{\displaystyle \frac{\alpha }{4\pi }}{\displaystyle \frac{1}{(1r)^2}}{\displaystyle \frac{1}{(1y+r)}}\{4A(r1)+[(1+r)^2+y(y4r)]\mathrm{ln}{\displaystyle \frac{y+A}{yA}}`$ (12)
$``$ $`(1y+r)^2\mathrm{ln}{\displaystyle \frac{y+A2}{yA2}}\}.`$
One can easily see that the eq. (12) is divergent both in the collinear and in the infrared limits. The coefficients of the collinear logarithms don’t go to zero in the limit $`r0`$. There are also infrared divergences, since if we let $`y`$ reach its kinematical limit $`y^{MAX}=1+r`$, corresponding to the photon energy going to zero, the expression (12) diverges.
The decay rate is made free from mass singularities in the ordinary way: the divergences cancellation occurs in the total inclusive decay rate, when we add all the first order contributions to the perturbative expansion, i.e. those relative to real and virtual photon emission.
The diagrams describing the real photon emission were already given in fig. 1; the diagram for the virtual correction is drawn in fig. 2.

The expression for the lepton energy spectrum, including the Inner Bremsstrahlung contribution and the virtual photon one, calculated to the leading order in $`m_l/m_\pi `$, is given by :

$$\frac{1}{\mathrm{\Gamma }_0}\frac{\mathrm{d}\mathrm{\Gamma }}{\mathrm{d}y}=D(y,r)\left[1+\frac{\alpha }{\pi }K_l(y)\right].$$
(13)
$`D(y,r)`$ is the lepton distribution function, given, to the first order in $`\alpha `$, by:

$$D(y,r)=\delta (1y)+\left[\frac{\alpha }{\pi }(L1)+O(\alpha ^2)\right]P^{(1)}(y),$$
(14)
where $`L`$ is the logarithm

$$L=\mathrm{ln}\frac{m_\pi }{m_l},$$
(15)
diverging in the collinear limit; if the lepton is an electron, $`L5.6`$. $`P^{(1)}(y)`$ is the Gribov-Lipatov-Altarelli-Parisi kernel , which can be expressed in the form:

$$P^{(1)}(y)=\frac{1+y^2}{1y}\delta (1y)_0^1dz\frac{1+z^2}{1z}.$$
(16)
$`K_l(y)`$ is a finite term, free from infrared and collinear singularities, which has the expression:

$$K_l(y)=1y\frac{1}{2}(1y)\mathrm{ln}(1y)+\frac{1+y^2}{1y}\mathrm{ln}y.$$
(17)
The differential decay rate to order $`\alpha `$ therefore becomes:

$$\frac{1}{\mathrm{\Gamma }_0}\frac{\mathrm{d}\mathrm{\Gamma }}{\mathrm{d}y}=\delta (1y)+\frac{\alpha }{\pi }(L1)P^{(1)}(y)+\frac{\alpha }{\pi }K_l(y).$$
(18)
To calculate the inclusive decay rate, we have to integrate the expression in eq. (18) over $`y`$; to the leading order in the lepton mass, the physical region for y is $`0y1`$. The kernel $`P^{(n)}`$ has the property that:

$$_0^1dyP^{(n)}(y)=0;$$
(19)
thus, when we calculate the inclusive decay rate, the coefficient of the collinear logarithm vanishes and the resulting expression is finite in the zero mass limit.
Carring out the integration over $`y`$, we obtain the well known inclusive decay rate to order $`\alpha `$

$$\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }_0}=1+\frac{\alpha }{\pi }[\frac{15}{8}\frac{\pi ^2}{3}].$$
(20)
As expected, the expression (20) is finite in the collinear limit and is also free from infrared divergences, because, as usual, the infrared divergences present in the soft photon contribution and in the virtual photon contribution have cancelled each other.
Let us now discuss the mass singularities in the case of the right-handed Inner Bremsstrahlung contribution, given in eq. (9). We have already seen in section 2 that $`\mathrm{d}\mathrm{\Gamma }_{IB}^R/\mathrm{d}y`$ is finite in the limit $`r0`$, i.e. it is free from collinear singularities. In this limit the coefficients of both the collinear logarithms vanish. Indeed, as we have discussed above, only the term related to the axial anomaly survives in the zero mass limit, that is the part of the decay rate indipendent of the lepton mass.
We observe that the right-handed Inner Bremsstrahlung contribution is free from infrared divergences as well. If we make the lepton energy $`y`$ reach its kinematical limit $`y^{MAX}`$, we obtain a finite result:

$`{\displaystyle \frac{1}{\mathrm{\Gamma }_0}}{\displaystyle \frac{\mathrm{d}\mathrm{\Gamma }_{IB}^R}{\mathrm{d}y}}0\mathrm{as}yy^{MAX}.`$ (21)
The result (21) shows that the emission of soft photons does not contribute to the radiative $`\pi ^+`$ decay with a right-handed lepton. Indeed, the soft photon contribution factorizes with respect to the Born decay rate, but this vanishes in the case of right-handed $`l^+`$ (see eq. (6)). Physically, eq.(21) is a consequence of the fact that soft photons don’t carry spin, thus they cannot contribute to the angular momentum balance; therefore the process with the right-handed lepton emitting a soft photon, is forbidden by angular momentum conservation.
For the same reason of angular momentum conservation, in the right-handed case also the virtual contribution is zero. The virtual correction is factorized with respect to the Born decay rate, but, as we have already seen, if the lepton is right-handed, this is identically zero.
In the right-handed case, the mass singularity cancellation occurs through a mechanism different from the one working in the unpolarized decay rate. The infrared and the collinear limits give, separately, a finite result. In particular, the coefficient of the collinear logarithms is the lepton mass, instead of the usual correction factor coming from the soft and the virtual photon contributions, as in eq. (18).
This particular mass cancellation mechanism is the consequence of the combination of two constraints: the angular momentum conservation in the pion vertex and the helicity flip in the photon-lepton vertex.
The situation is completely different if we consider the radiative process with the outgoing left-handed lepton, i.e. the process without helicity flip. The left-handed Inner Bremsstrahlung contribution is given by:

$`{\displaystyle \frac{1}{\mathrm{\Gamma }_0}}{\displaystyle \frac{\mathrm{d}\mathrm{\Gamma }_{IB}^L}{\mathrm{d}y}}={\displaystyle \frac{\alpha }{4\pi }}{\displaystyle \frac{1}{(1r)^2}}{\displaystyle \frac{1}{A(1y+r)}}\{2A[r(2Ar+6)y^212A]`$
$`+\left[(A2r)(1+r)^2+Ay(y4r)2ry(1+y)+y(5r^2+y^2+1)\right]\mathrm{ln}{\displaystyle \frac{y+A}{yA}}`$
$`+(1y+r)^2(2rAy)\mathrm{ln}{\displaystyle \frac{y+A2}{yA2}}\}`$ (22)
The expression (22) contains collinear singularities, since the coefficients of the collinear logarithms don’t vanish in the limit $`r0`$. It is also infrared divergent, as one can verify by taking the limit $`yy^{MAX}`$. In this case we don’t have the constraint constituted by the helicity flip in the photon-lepton vertex and in the pion vertex the angular momentum is conserved for soft and virtual photon emission. Thus, the mass singularity cancellation mechanism occurs in the ordinary way, as in the unpolarized case, i.e. in the total inclusive decay rate, obtained by adding all the order $`\alpha `$ contributions.
Let us show how the cancellation takes place. As we have already seen (eq. (6)), in the Born $`\pi ^+`$ decay the outgoing lepton is left-handed, due to angular momentum conservation. Thus the Born decay rates with unpolarized and left-handed $`l^+`$ coincide:

$$\mathrm{\Gamma }_0^L=\mathrm{\Gamma }_0.$$
(23)
Because of the factorization with respect to the Born decay rate, also the unpolarized and the left-handed virtual contributions are equal:

$$\mathrm{\Gamma }_v^L=\mathrm{\Gamma }_v.$$
(24)
Expressing the left-handed Inner Bremsstrahlung contribution in terms of the unpolarized and the right-handed ones, the total contribution to order $`\alpha `$ to the left-handed process is given by:

$$\mathrm{\Gamma }_{TOT}^L=(\mathrm{\Gamma }_0+\mathrm{\Gamma }_v+\mathrm{\Gamma }_{IB})\mathrm{\Gamma }_{IB}^R.$$
(25)
The expression (25) is finite both in the infrared and in the collinear limit, because the mass singularities present in the terms between brackets cancel each other, as we have seen (eq. (20)) and $`\mathrm{\Gamma }_{IB}^R`$ is free from mass singularities.
The presence of mass singularities is a conseguence of the fact that the states of a theory containing massless particles are highly degenerate. The KLN theorem states that the mass singularities disappear from the transition probability when we average it over the ensemble of degenerate states. It results that imposing on the outgoing lepton a polarization opposed to the one prescribed by the interaction taking place before the photon emission, implies a reduction of the degeneration subspace . Thus in this case we have a particular application of the KLN theorem.

### 3.2 The $`Z^0`$ case

Let us start the discussion on the $`Z^0`$ radiative decay by saying that it differs from the pion case. For radiative pion decay, due to the angular momentum conservation in the pion vertex, there is no room for a right-handed lepton. For such a channel, soft and virtual photon contribution give a contribution equal to zero. This result is valid independently of the lepton mass. Let us now consider a more general process, by loosing the value of the angular momentum of the decaying state. As an example, we study the radiative $`Z^0`$ decay in a lepton-antilepton $`(l^{}l^+)`$ pair, in which the lepton is in a definite helicity state.
The $`Z^0`$-leptons vertex is:

$$i\frac{M_Z}{\sqrt{2}}\left(\frac{G}{\sqrt{2}}\right)^{1/2}\gamma _\mu (g_vg_a\gamma _5),$$

with

$$g_v=14\mathrm{sin}\theta _{W}^{}{}_{}{}^{2}g_a=1$$

where $`\theta _W`$ is the Weinberg angle and $`M_Z`$ is the $`Z^0`$ mass.
If we set $`g_v=g_a=1`$, we require that in the limit of zero lepton mass, the $`Z^0`$ couples to a left-handed lepton. We calculate the decay rate for the process in which the outgoing lepton is right-handed.
At the Born level this is given by:

$`\mathrm{\Gamma }_0^R`$ $`=`$ $`{\displaystyle \frac{GM_Z^3}{48\sqrt{2}\pi }}\left\{\sqrt{14r}\left[(g_v^2+g_a^2)(1r)+3r(g_v^2g_a^2)\right]2g_vg_a(14r)\right\}.`$ (26)
If in eq. (26) we set $`g_v=g_a=1`$, $`\mathrm{\Gamma }_0^R`$ vanishes in the chiral limit, since there isn’t the term related to the axial anomaly.
Let us now study the decay process with the lepton emitting a real photon (see fig. 3)

and evaluate the probability that the outgoing lepton is right-handed. The electromagnetic interaction doesn’t couple states with different chirality, hence we expect the decay rate to vanish for $`m_l0`$.
The decay rate for the process described by the diagram of fig. 3, differential with respect to the lepton energy, is given by the expression:

$`{\displaystyle \frac{\mathrm{d}\mathrm{\Gamma }_\gamma ^R}{\mathrm{d}y}}={\displaystyle \frac{G\alpha M_Z^3}{96\sqrt{2}\pi ^2}}\{{\displaystyle \frac{(y2)(1y)^2}{4(1y+r)^2}}[(g_v^2+g_a^2)A+2g_vg_a(2ry)]`$
$`+`$ $`{\displaystyle \frac{(1y)}{2(1y+r)}}\left[(g_v^2+g_a^2)A(2ry)+2g_vg_a(y^22r)\right]`$
$`+`$ $`{\displaystyle \frac{2}{(y1)}}\left[2(g_v^22g_a^2)Ar+(g_v^2+g_a^2)A+2g_vg_a(4ry^2+y1)\right]`$
$`+`$ $`[g_v^2+g_a^2+2g_vg_a{\displaystyle \frac{1}{A}}(4r^2+r(yy^2+1)y)][\mathrm{ln}{\displaystyle \frac{y+A}{yA}}\mathrm{ln}{\displaystyle \frac{y+A2}{yA2}}]\}.`$
Here

$$r=\frac{m_l^2}{M_Z^2}$$

and $`y`$ is the usual dimensionless variable:

$$y=\frac{2E_1}{M_Z}$$

where $`E_1`$ is the lepton energy and

$$A=\sqrt{y^24r}.$$

The physical region for $`y`$ is

$$2\sqrt{r}y1.$$
(28)
The result obtained, as given by the emission of the photon by a single leg, is gauge dependent. To have a gauge independent amplitude, the contribution of the diagram b) of fig. 4 must be added . For the purpose of the polarized amplitude, however, the helicity flip contribution of the diagram b) of fig. 4 gives zero in the massless limit and is therefore negligible in our discussion.
From now on we set $`g_v=g_a=1`$, to have the condition of chirality conservation in the $`Z^0`$ vertex for $`m_l0`$. Taking this limit in eq. (LABEL:Rz), we see that $`\mathrm{d}\mathrm{\Gamma }_\gamma ^R/\mathrm{d}y`$ does not vanish:

$$\frac{\mathrm{d}\mathrm{\Gamma }_\gamma ^R}{\mathrm{d}y}\frac{G\alpha M_Z^3}{24\sqrt{2}\pi ^2}(1y)\mathrm{as}r0,g_v=g_a=1.$$
(29)
The result of this limit is the contribution related to the axial anomaly, which has the same form of the one found in the pion case.
Let us now discuss the mass singularity cancellation mechanism for the $`Z^0`$ radiative decay case. Eq. (29) shows that the collinear limit gives a finite result. Thus, we conclude that, as in the case of the radiative pion decay with right-handed lepton, the collinear and infrared limits are disconnected.
If we keep the lepton mass different from zero, the chirality is not fixed by the interaction occurring before the photon emission, even if we set $`g_v=g_a=1`$. Thus, for $`m_l0`$, the soft and virtual photon contributions are different from zero and diverge in the infrared limit. If we let the lepton energy reach its kinematical limit, $`y^{MAX}=1`$, we see that $`\mathrm{d}\mathrm{\Gamma }_\gamma ^R/\mathrm{d}y`$ diverges. We expect the infrared divergences to cancel, if we add all the first order contributions, given by the diagrams of fig. 3 and 4

and calculate the totally inclusive decay rate.
If we take the limit $`y1`$ in eq. (29), we obtain:

$$\frac{\mathrm{d}\mathrm{\Gamma }_\gamma ^R}{\mathrm{d}y}0\mathrm{as}m_l0,y1\mathrm{and}g_v=g_a=1.$$
(30)
Eq. (30) indicates that, in the massless limit, the soft photon contribution is zero. Indeed, it is factorized with respect to the Born decay rate, which, for $`r0`$ and $`g_v=g_a=1`$, vanishes. As we have discussed in section 1, the presence of the term related to the axial anomaly is directly connected to the emission of the photon, hence it vanishes in the limit in which the photon is not radiated.
The virtual photon contribution vanishes in the zero mass limit, as well. Indeed, it is factorized respect to the Born decay rate, which goes to zero as $`r0`$. The virtual correction factor can produce only a logarithmic collinear singularity, not a power-like one, needed for the cancellation of the chiral suppression.
In the massless limit, also the diagram with the photon emitted by $`l^+`$ doesn’t contribute, since, as the virtual contribution, cannot produces a power type collinear singularity.
Taking the infrared limit $`x0`$ in eq. (3), gives a finite result (indeed the cross section vanishes). This is a consequence of the fact that the cross section (3) has been calculated to the leading order in the lepton mass. From the eq. (LABEL:Rz), we see that the infrared divergent term is the one given by the term:

$$\frac{4r}{(y1)}\left[(g_v^22g_a^2)A+4g_vg_a\right]$$

i.e., it is proportional to the lepton mass. Performing the calculation, neglecting the mass terms as done in ref. , means imposing the chirality conservation law in the $`Z^0`$ vertex; thus the soft photon contribution is zero and the infrared divergences disappear. To the leading order in the lepton mass, one has (as in eq. 3) only the term related to the axial anomaly, which vanishes in the infrared limit.

## 4 Conclusions

We have seen that the decay rates for processes with a fermion changing helicity by emitting a massless vector boson don’t vanish in the chiral limit. We have interpreted this behaviour as related to the axial anomaly, as first suggested by Dolgov and Zakharov.
In the processes we have considered, the cancellation of the collinear singularities occurs through a mechanism different from the one for the unpolarized case. However, we have found a difference between the pion case and the more general case of the $`Z^0`$ decay. Due to the angular momentum conservation, the virtual and soft photon contributions are zero, even if the lepton mass is kept different from zero, giving finite infrared and collinear limits.
In the case of the radiative $`Z^0`$ decay, the decay rate is not finite in the infrared limit, since, for $`m_l0`$, the real soft photon contributions are different from zero. We observe that, in order to make the collinear limit finite, it is sufficient to sum over the set of degenerate states containing the fermion making the helicity flip accompanied by a hard collinear photon .
In these processes the collinear limit results disconnected from the infrared one, in the sense that the virtual photon contribution is not needed to render this limit finite.
We conclude that the collinear singularity cancellation mechanism is controlled by the anomalous breaking of the chiral symmetry. The axial anomaly implies that taking the collinear limit gives a finite result, independent of the fermion mass.
These results can be applied to interpret the mass singularity cancellation mechanism for any polarized process. The extension of this approach to the case of Quantum Cromodynamics for polarized processes is under study.

Acknowledgments

We wish to thank E. Kuraev and V. Fadin for valuable comments and S. Forte, J. Kodaira and L. Lipatov for useful discussions.


text:
# A Very Fast and Momentum-Conserving Tree Code

## 1 Introduction

The tree code (cf. Barnes & Hut 1986, hereafter B&H) has become an invaluable tool for the approximate but fast computation of the forces in studies of collisionless gravitational dynamics. It has been applied to a large variety of astrophysical problems. The gravitational potential generated by $`N`$ bodies of masses $`\mu _n`$ and at positions $`𝑿_n`$ is

$$\mathrm{\Phi }(𝑿)=\underset{n=1}{\overset{N}{}}\mu _ng(|𝑿𝑿_n|),$$
(1)
where $`g(r)`$ denotes the greens function, i.e. for un-softened gravity $`g(r)=G/r`$. The essence of the tree code is to approximate this sum over $`N`$ terms by replacing any partial sum over all bodies within a single cell which is well-separated from $`𝑿`$ by just one term. The inner structure of the cell is partly taken into account using its multipole moments. This method reduces the overall costs for the computation of all forces from $`𝒪(N^2)`$ to $`𝒪(N\mathrm{log}N)`$.

The tree code, however, does not exploit the fact that the force due to the contents of some cell is very similar at nearby positions (even though one may use the fact that nearby bodies tend to have very similar interaction lists, cf. Barnes 1990). Exploiting this is the idea of the fast multipole method (FMM) (Greengard & Rokhlin, 1987). The FMM employs a (usually) non-adaptive structure of hierarchical grids and considers only interactions between nodes on the same grid level according to their geometrical neighbourhood. The gravitational field due to some source cell and within some sink cell is approximated by a multipole expansion in spherical harmonics, the order of which is adapted to meet predefined accuracy limits. This method has been claimed to reduce the overall amount of operations to $`𝒪(N)`$, but the tables given by Cheng, Greengard & Rokhlin (1999) do not support this claim. Capuzzo-Dolcetta & Miocchi (1998) find that the FMM needs $`𝒪(N\mathrm{log}N)`$ operations, and is significantly slower for astrophysical applications than the tree code at comparable accuracy.

Instead of using a spherical multipole expansion of adaptive order, it is actually more efficient to use a Cartesian expansion of fixed order. Moreover, by preserving the symmetry of the gravitational interaction for mutual cell-cell interactions, one can (i) reduce the computational effort and (ii) obtain a code that satisfies Newton’s third law by construction and hence results in exact conservation of momentum, a property not shared by the traditional tree code.

## 2 Description of the Code

We start as the B&H tree code with a hierarchical tree of cubic cells. Each cell has up to eight sub-nodes corresponding to its octants. A node can be either a single body or another cell. The tree-building phase (cf. B&H) also includes the computation of the cells’ masses, centers $`𝒁`$ of mass, and quadrupole moments.

### 2.1 The Opening Criterion

In order to benefit from the symmetry of the gravitational interaction, the opening criterion, which decides whether or not two nodes are well-separated so that a direct mutual interaction is acceptable, must be symmetric, too. We employ an extension of the criterion used in the tree-code: nodes A and B are well-separated if

$$|𝒁_\mathrm{A}𝒁_\mathrm{B}|>(r_{\mathrm{max}\mathrm{A}}+r_{\mathrm{max}\mathrm{B}})/\theta ,$$
(2)
where the opening angle $`\theta `$ controls the accuracy of the code. $`r_{\mathrm{max}}`$ is the radius of a sphere centered on the node’s center of mass and encircling all bodies within it. Bodies naturally have $`r_{\mathrm{max}}0`$, i.e. two bodies are always well separated, while for the interaction between a body and a cell the criterion (2) reduces to that used in the tree code. Note that if one additionally to equation (2) requires $`r_{\mathrm{max}\mathrm{A}}=0`$, the standard tree code is re-covered, but the symmetry between A and B is broken.

There exist two upper limits for the radius $`r_{\mathrm{max}}`$. One is the distance $`b_{\mathrm{max}}`$ between the cell’s center of mass, $`𝒁`$, and its most distant corner (Salmon & Warren, 1994). The other is

$$\underset{\mathrm{sub}\mathrm{nodes}i}{\mathrm{max}}\left\{r_{\mathrm{max}i}+|𝒁_i𝒁|\right\}$$
(3)
(Benz et al., 1990). After computation of both these upper limits, we take the smaller one to be $`r_{\mathrm{max}}`$. For cells with only a few bodies like cell A in LABEL:example, the latter often gives values significantly smaller than $`b_{\mathrm{max}}`$, while for cells with many bodies, like cell B in LABEL:example, $`b_{\mathrm{max}}`$ is the tighter limit.

### 2.2 Approximating Gravity

Consider two bodies at $`𝑿`$ and $`𝒀`$ which reside in two well-separated cells A and B with centers of mass at, respectively, $`𝒁_\mathrm{A}`$ and $`𝒁_\mathrm{B}`$ and separation $`𝑹𝒁_\mathrm{A}𝒁_\mathrm{B}`$. We may re-write $`𝑿𝒀=𝑹+(𝒙𝒚)`$ with $`𝒙𝑿𝒁_\mathrm{A}`$ and $`𝒚𝒀𝒁_\mathrm{B}`$ being small in magnitude compared to $`𝑹`$ (because the cells are well-separated, cf. LABEL:example). The Taylor expansion of $`g(|𝑿𝒀|)`$ around $`𝑹`$ reads

$$g(|𝑿𝒀|)=\underset{p}{}\frac{1}{p!}\left[\left[(𝒙𝒚)\mathbf{}\right]^pg(|𝒓|)\right]_{𝒓=𝑹}.$$
(4)
Separating powers of $`x`$ from powers of $`y`$ in equation (4) and subsequently taking the mass weighted sum over cell B yields a Cartesian multipole expansion of the potential $`\mathrm{\Phi }_{\mathrm{B}\mathrm{A}}(𝑿)`$ at any position $`𝑿`$ within cell A and due to all bodies inside cell B (Warren & Salmon, 1995). Since $`𝒁_\mathrm{B}`$ was chosen to be the center of mass of cell B, its dipole vanishes. The highest-order multipole occurring in such an expansion may actually be omitted, since it only contributes a constant to the approximation for $`g`$ and does not affect the approximation for $`\mathbf{}g`$ and hence the force. The expression of third order thus reads (without the octopole; using Einstein’s sum convention)

$$\begin{array}{cccc}\hfill \mathrm{\Phi }_{\mathrm{B}\mathrm{A}}(𝑿)& M_\mathrm{B}\{\hfill & & 𝒟^{(0)}+\frac{1}{2}\stackrel{~}{Q}_{\mathrm{B}ij}𝒟_{ij}^{(2)}\hfill \\ & & +& x_i\left[𝒟_i^{(1)}+\frac{1}{2}\stackrel{~}{Q}_{\mathrm{B}jk}𝒟_{ijk}^{(3)}\right]\hfill \\ & & +& \frac{1}{2}x_ix_j𝒟_{ij}^{(2)}+\frac{1}{6}x_ix_jx_k𝒟_{ijk}^{(3)}\},\hfill \end{array}$$
(5)
where $`M_\mathrm{B}`$ and $`\stackrel{~}{Q}_\mathrm{B}`$ are the mass and specific quadrupole moment

$$\stackrel{~}{Q}_{\mathrm{B}ij}\frac{1}{M_\mathrm{B}}\underset{𝒚_n+𝒁_\mathrm{B}\mathrm{cell}\mathrm{B}}{}\mu _ny_{ni}y_{nj}$$
(6)
of cell B, while $`𝒟^{(n)}\mathbf{}^ng(r)|_{r=|𝑹|}`$, i.e.

$$\begin{array}{ccc}\hfill 𝒟^{(0)}& =& D^0,\hfill \\ \hfill 𝒟_i^{(1)}& =& R_iD^1,\hfill \\ \hfill 𝒟_{ij}^{(2)}& =& \delta _{ij}D^1(R)+R_iR_jD^2,\hfill \\ \hfill 𝒟_{ijk}^{(3)}& =& \left(\delta _{ij}R_k+\delta _{jk}R_i+\delta _{ki}R_j\right)D^2+R_iR_jR_kD^3\hfill \end{array}$$
(7)
with

$$D^n\left(\frac{1}{r}\frac{}{r}\right)^ng(r)|_{r=|𝑹|}.$$
(8)
The symmetry between $`𝒙`$ and $`𝒚`$ at every order of the Taylor expansion in equation (4) has two important consequences. First, if this expansion is used to compute both $`\mathrm{\Phi }_{\mathrm{B}\mathrm{A}}(𝑿)`$ and $`\mathrm{\Phi }_{\mathrm{A}\mathrm{B}}(𝒀)`$, Newton’s third law is satisfied by construction. Note, that our omission of the octopole term broke the symmetry only in the zeroth order and has no effect on the forces.

Second, the expressions for $`M_\mathrm{A}\mathrm{\Phi }_{\mathrm{B}\mathrm{A}}`$ and $`M_\mathrm{B}\mathrm{\Phi }_{\mathrm{A}\mathrm{B}}`$ are very similar: the expansion coefficients of second and third order differ only by mere signs, such that computing these coefficients for both Taylor series at one time is substantially faster than computing them at different times.

The transformation, or shifting, of the expansion center to some other position is trivial compared to the analogous procedure in the FMM (see Cheng et al.).

### 2.3 The Algorithm

The standard tree code computes the forces on each body by a recursive tree walk, which visits each node exactly once as a gravity sink, and thus exhibits an inherent asymmetry between sources and sinks. The new algorithm avoids this asymmetry.

First, in the interaction phase, the Taylor series coefficients are evaluated and accumulated in data fields associated with each node. This phase is based on the concept of mutual interactions (MIs), pairs of nodes, A and B, such that bodies in node A must receive forces from all bodies in node B, and vice versa. We start by the MI describing the root-root self-interaction, and process a given MI as follows. (1) A body self-interaction is ignored; (2) a cell self-interaction is split into the MIs between the sub-nodes111 In a cubic oct-tree, these are at most 36 independent sub-MIs., and the process is continued on each of the new MIs; (3) a MI representing a well-separated pair of nodes is executed: the Taylor coefficients are computed and added to the nodes’ corresponding data fields; (4) finally, in any other case, the node with larger $`r_{\mathrm{max}}`$ is split, and up to eight new MIs are created and processed.

Secondly, in the collection phase, the Taylor coefficients are passed down the tree: the expansion center is shifted to the center of mass of the currently active cell and the coefficients are accumulated. The Taylor expansion is evaluated at the position of any body and the values for potential and acceleration are added to its data fields (which may already contain contributions accumulated during the interaction phase).

## 3 Performance Tests

We tested the new algorithm and compared it with the tree code in three typical astrophysical situations: (1) a spherical Plummer model, representing a rather homogeneous stellar system, (2) a spherical Hernquist-model (1990) galaxy, and (3) a group of five such galaxies with various masses and scale radii. We generated $`10^5`$ random initial positions from each of these cases, truncating the density at 1000 scale radii, and evaluated the exact mutual forces at all positions and the approximated forces due to the tree code (up to quadrupole order) and the new code for opening angles $`\theta `$ between 0.2 and 1. We used an optimally chosen softening with the biweight softening kernel (see Dehnen 2000), but the results are insensitive to these settings. Both approximate methods have been coded by the author222 At the same $`\theta `$, the tree code is twice as fast as a code publicly available from J. Barnes, mainly because the new opening criterion leads to fewer interactions. However, even at the same number of interactions, the author’s code was about 30% faster. and use the same opening criterion (§2.1). In order to measure the accuracy of the approximated forces, we evaluated (cf. Capuzzo-Dolcetta & Miocchi 1998)

$$\epsilon _n=|a_na_n^{PP}|/a_n^{PP},$$
(9)
where $`a_n`$ denotes the magnitude of the acceleration of the $`n`$th body due to either of the approximate methods and $`a_n^{PP}`$ that of the exact computation.

Figure 2 plots the CPU time needed for the force approximation on a Pentium III/500Mhz PC versus the mean relative error and that at the 99 percentile. For the new code, the time consumption scales almost inversely with the error, while the tree code flattens off333 This is, because of $`r_{\mathrm{max}}`$ is not proportional to the cell size, so that increasing $`\theta `$ from $`0.7`$ does not much decrease the number of interactions. at $`\theta 0.7`$. Evidently, at an acceptable level of accuracy, e.g. $`\epsilon _{99\%}=0.01`$, the new code is about four times faster than the tree code, even though it requires a smaller value of $`\theta `$ (0.5 as compared to 0.7 for the tree code).

For the test case of the Hernquist model, Figure 3 plots the CPU time consumption per body versus $`N`$ at fixed $`\theta =0.7`$ for the tree code and $`0.5`$ for the new code. The tree code shows the well-known $`\mathrm{log}N`$ scaling, while for $`N10^5`$ the new code requires only a constant amount of CPU time per body. This can be explained as follows. Arranging eight root cells to a new root box, increases $`N`$ to $`8N`$ and the number $`N_I`$ of interactions to $`8(N_I+N_+)`$, where $`8N_+`$ interactions are needed to compute the forces between the former root cells. Thus,

$$\frac{\mathrm{d}N_I}{\mathrm{d}N}\frac{N_I}{N}\frac{\mathrm{\Delta }\mathrm{log}N_I}{\mathrm{\Delta }\mathrm{log}N}=\frac{N_I+N_+/\mathrm{ln}8}{N}.$$
(10)
In the tree code, $`N_+N`$ yielding $`N_IN\mathrm{log}N`$. In the new code, $`N_+`$ may be estimated to contain two contributions, a constant term accounts for the interactions with distant nodes, and a term $`N^{2/3}`$ for those on the surface of the former root cells. Inserting this into equation (10) yields

$$N_INc_1N^{2/3}c_2$$
(11)
with constants $`c_1`$ and $`c_2`$ that depend on $`\theta `$. Thus, at large $`N`$, a linear relation is approached. Note that this argument differs from that given for the FMM by Greengard & Rokhlin (1987), who assumed that the resolution may remain fixed when increasing $`N`$.

## 4 Discussion

A new code for the approximate evaluation of gravitational forces has been presented, tested, and compared to the tree code. This new code is substantially faster than the tree code. Moreover, unlike the latter, it satisfies Newton’s third law by construction, such that any $`N`$-body code based on it will not introduce spurious net-accelerations. The new code is based on a Taylor expansion of the greens function in Cartesian coordinates and incorporates mutual cell-cell interactions. The simple algorithm is well suited for implementation on parallel computers: different mutual interactions (MIs) can be passed to different CPUs.

The scaling of the CPU time required for the mutual forces of a number $`N`$ of bodies becomes essentially linear at $`N10^5`$, so that with ever larger $`N`$ the new code is increasingly faster than the tree code, allowing for a substantial improvement in simulations employing large number of bodies. The only disadvantage is the increased requirement of memory compared to the standard tree code: 20 floating point numbers per cell are needed to hold the Taylor expansion coefficients. (By using a tree-walking algorithm instead of that given in §2.3, one can avoid this at the price of enhanced CPU time consumption.)

In spirit, the new code is similar to Greengard & Rokhlin’s (1987) fast multipole method, but is more efficient because it (i) uses a Cartesian instead of a spherical harmonic multipole expansion and (ii) fixes the order of the expansion while controlling the accuracy via the interaction condition, rather than fixing the interactions and adapting the expansion order to the accuracy.

A concern with codes based on cell-cell interactions is their performance in the presence of individual time steps. Clearly, when not all the forces are to be computed, such codes fare less favorably. However, when the forces for all bodies within some domain are desired, the new code is still a significant improvement over the tree code.

The new code has been written in C++ and will be electronically available from the author upon request.

The author thanks Tom Quinn for valuable discussions and Joshua Barnes, Lars Hernquist, Junichiro Makino, Andy Nelson, and Thorsten Naab for useful comments.