text:
# Singularity formation in the Gross-Pitaevskii Equation and Collapse in BEC

## I Introduction

Collapse phenomena in Bose-Einstein condensates received a great deal of interest in the last few years Refs.pit ; x ; xx ; xxx ; xxxx ; xxxxx ; xxxxxx ; sackett . The experimental technique based on a Feshbach resonance allows to tune, through application of an external magnetic field, the strength of the interatomic interaction and even switch between repulsive and attractive interactionscornish . The realization of a collapse controlled by the magnetic field on a gas of $`{}_{}{}^{85}\text{Rb}`$ reveals a rich dynamical properties of the collapsing condensates Ref.donley2001 .

When the number of atoms becomes sufficiently large, the interatomic energy overcomes the quantum pressure and the condensate implodes. In the course of the implosion stage the density increases in the small vicinity of the trap center. When the density approaches certain critical values a fraction of the atoms gets expelled. In a time period of an order of few milliseconds the condensate stabilizes. There are two observable components in the final stage of the condensate dynamics: remnant and burst atoms. The remnant atoms are the atoms which remain in the condensate. The burst atoms have the energy which is much larger then the energy of the condensed atoms. There is also a fraction of atoms, which is not observable. This fraction is referred to as missing atoms. To account for dissipative losses in the collapse phenomena, the most part of theoretical studies reflected in existing literature try to modify the mean field theory through addition of damping terms modelling the mechanism of dissipative losses in the system. In particular, the results of experiment of Ref.donley2001  have been analyzed by numerical solution of the Gross-Pitaevskii (GP) equation with three body recombination term Refs.saito2001 ; santos2001 ; abhi2002 . A mechanism of losses through elastic collisions was discussed in Ref.stoof .

The References saito2001 ; santos2001  follow the ideas of experimental work donley2001  and obtain the burst atoms fraction fitting the Gaussian distribution to the condensate wave function. The essence of the approach of Refs. saito2001 ; santos2001  is the idea that the final state of the condensate is comprised of the remnant fraction surrounded by a very dilute fraction of burst atoms. Hence, the only fraction which is expelled from the condensate is that of the missing atoms, while the burst atoms fraction is formed in the course of the condensate dynamics.

Effective workings of a nonlinear mechanism of dissipation in description of collapse in BEC are inherently connected to the structure of singularity zakh1991  which is formed in the course of the nonlinear dynamics as described by the Gross-Pitaevskii equation pit . This equation reads

$$i\mathrm{}\mathrm{\Psi }_t+\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\mathrm{\Psi }g\mathrm{\Psi }|\mathrm{\Psi }|^2V(𝐱)\mathrm{\Psi }=0.$$
(1)
Here $`\mathrm{\Psi }(𝐱,t)`$ is the wave function of the condensate, the external potential $`V(𝐱)`$ models the wall-less confinement (the trap), $`m`$ is the mass of an individual atom, $`g=4\pi \mathrm{}^2m^1a_s`$ is the effective interaction strength (the ’coupling constant’), $`a_s`$ is the scattering length, and $`\mathrm{\Delta }=_i\frac{^2}{x_i^2}`$ is the Laplace operator. A convenient as well as practical choice for the confining trap is the paraboloidal potential $`V=\frac{m}{2}_i^3\omega _i^2x_{i}^{}{}_{}{}^{2}`$. Notice that so far very little is known about the singularity structure of the $`3D`$ GP equation Eq.(1) outside treatments of spherically symmetric case (see the monograph sulem  and references therein). In this work we will try to gain an insight into the problem through reformulation of the Gross-Pitaevskii equation Eq.(1) in a different form. Our approach resembles, at least formally, the de Broglie-Bohm formulation of quantum mechanics holland ; bohm .

We introduce quantum trajectories ($`q`$-trajectories) of the GP equation Eq.(1). These trajectories satisfy a system of ordinary differential equations, viz.

$$\frac{r_i}{t}=\frac{1}{m}\frac{\varphi }{x_i}(𝐫(𝜼,t),t);r_i(𝜼,0)=\eta _i,$$
(2)
with $`i=1,2,3`$ in the $`3D`$ case. Here $`\varphi `$ is the phase defined by the representation $`\mathrm{\Psi }=\sqrt{\rho }e^{i\varphi /\mathrm{}}`$ of the condensate field. The fields $`r_i(𝜼,t)`$ can now be interpreted as trajectories of the fictitious particles with the initial points located at $`\eta _i`$. These trajectories satisfy the Newton equations of motion for a set of particles,

$$m\frac{^2r_i}{t^2}=\frac{V}{x_i}\frac{U}{x_i},$$
(3)
where the ’pressure potential’ $`U(𝐱,t)`$ is

$$U(𝐱,t)=\frac{\mathrm{}^2}{2m}\rho ^{1/2}\mathrm{\Delta }\rho ^{1/2}+g\rho .$$
(4)
In addition, the Liouville formula for Eq.(2) can be written in the form

$$det\left(\frac{r_i}{\eta _k}\right)\rho (𝐫,t)=\rho _0(𝜼),$$
(5)
where $`\rho _0`$ is the initial density of the condensate. The system of equations Eqs.(3), (5) is equivalent to the original GP equation Eq.(1).

As was indicated above the $`q`$-trajectories are used in the de Broglie-Bohm causal interpretation of quantum mechanics Refs.holland ; bohm . In our previous work Ref.rybin2001  we employed these trajectories to define singularities of the GP equation Eq.(1). Indeed, Eq.(5) immediately suggests that a set of singular points $`𝐫_s`$ of the GP equation Eq.(1) at which $`\rho (𝐫_s,t_{})=\mathrm{}`$ for some fixed $`t_{}`$ in terms of $`q`$-trajectories, describes either the caustic or focal points of the trajectories. Notice that in the case of linear  Schrödinger equation the probability density $`\rho `$ is always finite and caustics do not exist.

The Newton equations Eqs.(3) define the flux of particles. The appearance of the burst atoms in the mean field theory can be qualitatively understood through analysis of the energy distribution near the singularity in the flux of particles as shown in Figure 1. In this picture the peak with positive energy density arises due to the acceleration of the particles by the singularity. This acceleration is proportional to $`g\rho `$ and results in increase of the kinetic energy of the particles surrounding the singularity. In Figure 1 the curves with greater peaks correspond to the moments of time closer to $`t_{}`$. The burst atoms can be formed only if the dissipation losses occur faster than the collapse domain is refilled by the flux. In this case one observes intermittent implosions as reported in Ref.saito2001 .

The burst atoms cannot appear in the framework of the mean field theory if the damping parameter is chosen to be large enough (cf. Ref.abhi2002 ). In this case the density cannot increase to sufficiently large values and the peak of the positive energy density is not formed.

The dynamics of the $`q`$-trajectories which begin at the points located far from the origin in the $`𝜼`$-space cannot be significantly affected by possible dissipative term. This is because the compression time of the breathing mode of the harmonic oscillator, which is $`1/4`$ of the period of the oscillator, is of the same order of magnitude as the time of the collapse including the time of dissipation. This means that the remote trajectories cannot reach the collapse domain.

In the presence of the burst atoms the condensate dynamics outside the dissipation domain is a quite complicated interference of the ’outgoing’ burst atoms and ’incoming’ atoms. Should we neglect this interference, the fraction of the remaining atoms could be found through the analysis of the particles flux falling into the singularity and described by the Newton equations Eqs.(3). Notice that these equations are still valid in the presence of dissipative term. However in this case the relation Eq. (5) does not hold.

The analysis of flux of particles falling into the singularity shows that the fraction of atoms remaining in the condensate is consistently constant for the wide range of parameters. This effect is supported by an experimental evidence donley2001 . In this work we deal with the general case of asymmetric trap. In our previous work Ref. rybin2001  we reported for the symmetric trap an effect similar to the described above. Our computation was based on a Gaussian trial wave function which was used to estimate a part of atoms focusing on the caustic. Here we discuss different types of collapse which can exist in asymmetric trap.

The paper is organized as follows. In the next section we discuss different types of collapse which can occur in the asymmetric trap. In the section 3 we obtain the collapsing fraction and show that this fraction in the broad range of parameters consistently remains almost constant.

## II Breathing mode of the condensate

In this section we consider types of collapse which can be observed in an asymmetric trap. We use as the initial condition the Gaussian profile. For this profile

$$\rho _0=\frac{N}{\pi ^{3/2}}\underset{i=1}{\overset{3}{}}\frac{1}{\mathrm{w}_i}e^{x_i^2/\mathrm{w}_i^2}$$

the harmonic oscillator equation has explicit solution, viz.

$$r_i(𝜼,t)=\sigma _i(t)\eta _i,$$
(6)
where

$$\sigma _i(t)=\sqrt{\mathrm{cos}^2(\omega _it)+\frac{a_i^4}{\mathrm{w}_i^4}\mathrm{sin}^2(\omega _it)},a_i=\sqrt{\frac{\mathrm{}}{m\omega _i}}.$$

It seems natural, at least for Gaussian initial conditions, to expect that the breathing mode described by Eq. (6) gives the main contribution into the action

$$S(𝐫)=𝑑t𝑑𝜼\rho _0(𝜼),$$

where

$$=\underset{i=1}{\overset{3}{}}\frac{m}{2}\dot{r}_i^2\frac{\mathrm{}^2}{2m\rho }(\sqrt{\rho })^2V\frac{g}{2}\rho .$$

This provides us with the motivation to employ the variational principle and to seek for the solution of the GP equation Eq. (1) in the form $`r_i(𝜼,t)=\tau _i(t)\eta _i`$. Substituting this into the action we obtain

$$\ddot{\tau }_i+\omega _i^2\tau _i\frac{\mathrm{}^2}{m^2}\frac{1}{\tau _i^3}\frac{gN}{(2\pi )^{3/2}m}\frac{1}{_k\tau _k}\frac{1}{\tau _i}=0$$
(7)
The density as obtained through the relation Eq. (5) reads

$$\rho (x,t)=\frac{N}{\pi ^{3/2}}\underset{i=1}{\overset{3}{}}\frac{1}{\tau _i}e^{x_i^2/\tau _i^2}$$
(8)
Here we have included $`\mathrm{w}_\mathrm{i}`$ in $`\tau _i`$ and used $`\rho _0`$ with $`\mathrm{w}_\mathrm{i}=1`$.

The result Eq. (8) is well-known Gaussian trial wave function Ref.gar . This function describes the breathing (monopole) mode of the condensate.

In what follows we measure the time in units of $`1/\omega `$, where $`\omega =(\omega _x\omega _y\omega _z)^{1/3}`$, and the distance in units of the the oscillator length $`a_{HO}=\sqrt{\mathrm{}/m\omega }`$. For the axial symmetric equation we have $`\omega _x=\omega _y=\omega _r`$ and $`\tau _x=\tau _y=\tau _r`$. Then the system Eqs.(7) is reduced to the following equations

$$\ddot{\tau }_z+\beta ^2\tau _z\frac{1}{\tau _z^3}\sqrt{\frac{2}{\pi }}\frac{\kappa }{\tau _z^2\tau _r^2}=0$$
(9)
$$\ddot{\tau }_r+\frac{1}{\beta }\tau _r\frac{1}{\tau _r^3}\sqrt{\frac{2}{\pi }}\frac{\kappa }{\tau _z\tau _r^3}=0$$
(10)
Here $`\beta =\frac{\omega _z}{\omega }`$ and $`\kappa =\frac{Na_s}{a_{HO}}.`$

For negative $`\kappa `$, such that $`|\kappa |`$ is greater than certain critical value, the system of Eqs.(9), (10) has singular solutions with $`\tau _r0`$ as $`tt_{}`$. We found three different regimes of collapse.

Spherically symmetric collapse is characterized by $`\tau _z\tau _r`$ as $`tt_{}`$. The functions $`\tau _r,\tau _z`$ near the singularity are described by the power law

$$\tau _z\tau _{}(t_{}t)^{2/5},\tau _r\tau _{}(t_{}t)^{2/5}$$

where $`\tau _{}`$ is some constant. In this regime the compression of the condensate is uniform in all directions.

Radial collapse is effectively two dimensional collapse for which the condensate contracts only in the radial direction. The behavior of $`\tau _z,\tau _r`$ near $`t_{}`$ is defined by

$$\tau _z\lambda ,\tau _r\tau _r(t_{}t)^{1/2}.$$
(11)
Here $`\lambda `$ is a constant and $`\tau _r=\sqrt{2}\left(\sqrt{\frac{2}{\pi }}\frac{|\kappa |}{\lambda }1\right)^{1/4}.`$

Axial collapse is characterized by $`\tau _z\tau _r^2`$ as $`tt_{}`$. In this case the compression of the condensate is faster in the axial direction. The function $`\tau _z`$ approaches zero through an oscillating regime. This regime cannot be described by a power law.

The actual type of singularity depends on the initial conditions for Eqs.(9),(10). To begin the analysis we specify the parameters and initial data as $`\beta =0.5325,\dot{\tau }_z(0)=0,\dot{\tau }_r(0)=0,\tau _z(0)=\tau _{z0}`$ and $`\tau _r(0)=\tau _{r0}`$. The constants $`\tau _{z0},\tau _{r0}`$ are a stationary solution of the system of Eqs.(9), (10) with $`\kappa =\kappa _{in}>0`$. This choice of parameters models experiments on Feshbach resonances when the scattering length is instantly changed from positive to negative values in the initial stage of the dynamics.

For accurate implementation of numerical simulation we use the following change of variables

$$\tau _z=e^{\nu _z(\theta )},\tau _r=e^{\nu _r(\theta )},\theta =\mathrm{ln}\frac{t_c}{t_ct},$$

where $`t_c`$ is a constant. The functions $`\nu _r,\nu _z`$ are regular in the interval $`[0,+\mathrm{})`$ if $`t_ct_{}`$ and acquire a singularity at some $`\theta _s<+\mathrm{}`$ if $`t_c>t_{}`$.

We also obtain $`\lambda `$ as a function of $`\kappa ,\kappa _{in}`$ as shown in Figure 2. In this Figure the open domain $`\mathrm{\Omega }_\lambda `$ in $`\kappa \kappa _{in}`$ plane consists of two disconnected components $`\mathrm{\Omega }_\lambda ^{(1,2)}`$ $`\left(\mathrm{\Omega }_\lambda =\mathrm{\Omega }_\lambda ^{(1)}\mathrm{\Omega }_\lambda ^{(2)}\right)`$ located on the different sides of the groove. The overall domain $`\mathrm{\Omega }_\lambda `$ corresponds to the radial collapse. The component $`\mathrm{\Omega }_\lambda ^{(1)}`$ corresponds to the simple radial collapse. In this subdomain the attraction between atoms is strong, i.e. $`|\kappa |`$ is large. This results in fast and monotonous decreasing of the Gaussian widths $`\tau _i,i=x,y,z`$. Notice that while $`\tau _{x,y}`$ vanish, $`\tau _z`$ remains finite as in Eq.(11). In the second component of $`\mathrm{\Omega }_\lambda `$, i.e. in $`\mathrm{\Omega }_\lambda ^{(2)}`$ the density $`\rho (0,t)`$ has at least one local maximum in the interval $`(0,t_{})`$. In this regime the collapse occurs after a regular compression-expansion cycle. The rates of contraction in the radial and axial directions are alternating with each other in magnitudes. At some instances an intensive contraction in one direction is compensated by expansion in the other. Notice that the final stage of this process is always the contraction in the radial direction. This is why we refer to this type of collapse as to the radial collapse. We turn next to the description of the open domain complimentary to the domain $`\mathrm{\Omega }_\lambda `$ in the $`\kappa \kappa _{in}`$ plane. This is the groove shown in Figure 2. This domain corresponds to the axial collapse. In the case of axial collapse the ultimate stage of the dynamics is always the faster contraction in the axial direction. A typical behavior in axial collapse is shown in Figure 3.

The boundary separating the open domain $`\mathrm{\Omega }_\lambda `$ and the groove, i.e. $`\mathrm{\Omega }_\lambda `$ corresponds to the spherically symmetric collapse. It is interesting to notice that even though the trap is asymmetric, the workings of nonlinear dynamics still allow the spherical collapse to occur.

Intuitively, the existence of a regular cycle in dynamics of the breathing mode for $`(\kappa ,\kappa _{in})\mathrm{\Omega }^{(2)}`$ suggests that fewer number of atoms must be removed to stabilize the condensate. We further discuss this idea in the next section where we analyze the flux of particles.

## III The flux of particles

The $`q`$-trajectories $`r_i(𝜼,t)=\tau _i\eta _i`$ derived in the previous section focus at the symmetry axis when the time approaches $`t_{}`$ for any $`𝜼`$. However, for the large $`|𝜼|`$ this is not true, since in this region the behavior of the function $`𝐫(𝜼,t)`$ is determined by the solution of the harmonic oscillator equation Eq.(6).

The scope of analysis of equations for quantum trajectories Eqs.(3) proposed in this section attempts go beyond the limitations posed by the Gaussian ansatz. In what follows we construct an effective approximation for the flux of particles. Our approximation effectively uses the notion of quantum trajectories. We gain a motivation for this approach from the analysis of the quasi-classical approximation for the Schrödinger equation. To construct the quasi-classical approximation of the Schrödinger equation we have to solve the Newton equations for the trajectories of particles, i.e. to solve $`\ddot{r}_i=\frac{V}{r_i}`$. The amplitude and the phase of the wave function are determined through Eqs.(2),(5). We can search further corrections substituting the quasi-classical amplitude into the Newton equations Eqs.(3) with $`g=0`$. Iterating this procedure and assuming that iterations converge we finally get the solution of the Schrödinger equation in the time interval $`[0,T]`$ such that the function $`𝐫(𝜼,t)`$ with $`t[0,T]`$ is invertible at each step. For instance, in the case of harmonic oscillator and Gaussian initial conditions this procedure converges in the interval defined by the inequality $`_i\mathrm{cos}(\omega _it)>0`$. In fact we construct here an iteration procedure $`F`$ of the form $`𝐫^{(n)}=F\left(𝐫^{(n1)}\right)`$, where $`n`$ is the number of iteration. These iterations start from the trajectories of the classical harmonic oscillator $`r_i^{(0)}=r_i^{cl}(𝜼,t)=\mathrm{cos}(\omega _it)\eta _i`$. Notice that the mechanism of the iteration procedure $`F`$ is disconnected from the quasi-classical approximation. This procedure is equally well applicable in both linear and nonlinear ($`g0`$) cases. Indeed, solving Eqs.(3) in the nonlinear case, we can then find from Eq.(5) a next iteration for the function $`\rho `$ which we again employ in Eqs.(3), and so forth.

Below we introduce an effective approximation for the quantum trajectories as the first step of the iteration procedure $`F`$ described above. We assume that the starting solution Eq.(8) delivered by the variational principle is reasonably close to the exact solution. Using Eqs.(3) with the density given by equation Eq.(8) and substituting Eq.(8) into Eqs.(3) we obtain

$$\frac{^2r_i}{t^2}=\frac{V_{eff}}{r_i},r_i(0)=\eta _i,\dot{r}_i(0)=0,$$
(12)
where

$$V_{eff}(𝐫,t)=\frac{1}{2}\underset{i}{}\left(\frac{\omega _i^2}{\omega ^2}\frac{1}{\tau _i^4}\right)r_i^2+4\pi \kappa \rho (𝐫,t)/N.$$

Notice that for the case of an axially symmetric trap it is sufficient to consider the behavior of the $`q`$-trajectories in $`x`$-$`z`$ plane. This can be achieved through reduction $`r_y=0`$.

At the moment of time $`t_c`$ close to $`t_{}`$ we observe two main types of trajectories: the trajectories which go away from the center and the trajectories which tend to focus on the symmetry axis. The domain in the $`𝜼`$-space corresponding to the focussing trajectories is denoted by $`\mathrm{\Omega }_s`$. We now define the collapsing fraction of atoms by an integral

$$N_e=\underset{\mathrm{\Omega }_s}{}\rho _0𝑑𝜼.$$
(13)
Here $`\rho _0`$ is the initial density. It is physically plausible to assume that the collapsing fraction is removed from the condensate due to the dissipation. Then the remaining fraction is given by the formula $`1N_e/N`$.

We analyze numerically the flux of particles falling into the singularity. To do this we solve Eqs.(9), (10), (12). At the moment $`t=t_c`$ the focussing trajectories with reasonable accuracy get inside the domain $`\mathrm{\Omega }_c`$. This domain is restricted by the ellipsoid of revolution of an ellipse with the semi $`x`$ and $`z`$-axes equal $`\tau _z(t_c)`$ and $`\tau _r(t_c)`$ respectively. For $`t_ct_{}`$ the domain $`\mathrm{\Omega }_c`$ reduces to the interval $`[\lambda ,\lambda ]`$ in the $`z`$-axis.

We now technically define the domain $`\mathrm{\Omega }_s`$ in the $`𝜼`$-space such that if $`𝐫(𝜼,t_c)\mathrm{\Omega }_c`$ then $`𝜼\mathrm{\Omega }_s`$. The results of simulations are presented in Figure 4. We find that the fraction $`N_e/N`$ is consistently around 70% for large values of $`|\kappa |`$ and depends neither on $`\kappa _{in}`$ nor on $`N`$. More exactly this is true for the subdomain $`\mathrm{\Omega }_\lambda ^{(1)}`$ corresponding to the simple radial collapse described above. The collapsing fraction decreases abruptly for $`|\kappa |<|\kappa _c|`$, where $`\kappa _c`$ is a certain threshold value specific to the given value of the parameter $`\kappa _{in}`$. These thresholds for three different values of $`\kappa _{in}`$ are shown in Figure 4. As was explained above, below the threshold $`|k_c|`$, in the domain $`\mathrm{\Omega }_\lambda ^{(2)}`$ and in the groove, we observe oscillatory dynamical patterns. In the course of these oscillations we observe losses in the collapsing fraction. For large $`\kappa _{in}`$ and in this region the fraction $`N_e/N`$ is found to be around 20% .

## IV Conclusions

In this paper we suggest a simple model Eq. (12) of the BEC collapse based on the dynamical properties of the GP equation. Our model is independent of the damping parameter and a microscopic mechanism of the atom ejection. In the framework of this model it is impossible to correctly interpret the burst atoms fraction. However, the simplicity of the model allows us to study the problem in wide range of the parameters.

Using the notion of the flux of particles we define the collapsing fraction $`N_e`$ as the part of atoms focusing at the symmetry axis. This focusing results in infinite increase of the density. Because of this we assume that total collapsing fraction is removed from the condensate. We estimate that $`N_e/N=0.7`$ for sufficiently large $`\kappa `$ and the remaining fraction in this case is insensitive to the choice of parameters. Interestingly, our results show that changing the initial aspect ratio by $`\kappa _{in}`$ it is possible to move from the domain $`\mathrm{\Omega }_\lambda ^{(1)}`$ to the domain $`\mathrm{\Omega }_\lambda ^{(2)}`$ for the same interaction strength $`\kappa `$. This means that the fraction of remaining atoms can be increased for the larger initial aspect ratio controlled by $`\kappa _{in}`$. In our opinion this effect deserves experimental investigation.

## Acknowledgements

This work has been supported by the Academy of Finland (Project No. 44875). I.P.V. is grateful to the Centre for International Mobility (CIMO) for financial support. G.G.V. is grateful to the University of Jyväskylä, Finland for financial support.


format: text

text:
# The origin of radio haloes and non-thermal emission in clusters of galaxies

## 1 Introduction

Since the late 1960s, clusters of galaxies have been known to be strong X-ray sources where the thermal emission from the intracluster medium (ICM) produces a diffuse X-ray halo. The ICM is diffuse (central electron density of $`10^3`$ cm-3) with typical temperature $`2\times 10^7`$ to $`2\times 10^8`$ K. In the 1970s, similarly diffuse radio emission was found in the Coma cluster. This radio emission was found to have a steep power law index ($`\alpha 1`$ for $`S_\nu \nu ^\alpha `$) which indicates a non-thermal origin. Despite efforts to detect similar radio haloes in other clusters, until recently few were found with diffuse radio emission (see review by Feretti & Giovannini 1996). This led to the belief that radio haloes are rare. However, the number of clusters found to have diffuse radio emission has more than doubled over the last couple of years, mainly because of the availability of moderate brightness sensitivity surveys at 1.4 GHz, e.g. the NVSS \[Giovannini et al. 1999\], as well as high brightness sensitivity observations of hot clusters to detect the Sunyaev-Zel’dovich effect (e.g. Herbig & Birkinshaw 1994; Liang et al. 2000).

Diffuse cluster radio emission is now divided into two distinct classes: haloes and relics. Haloes permeate the entire cluster and resemble the thermal X-ray emission, but relics are on the cluster periphery. Both haloes and relics exhibit steep spectral indices ($`\alpha 1`$) and are extended over $`1`$ Mpc with the haloes having lower surface brightness. Emission from radio relics is strongly polarised (up to $`20`$%) but no polarised emission has been detected in radio haloes. Radio haloes and relics are believed to emit by the synchrotron process. The non-detection of polarised emission from radio haloes is attributed to low surface brightness and to Faraday depolarisation effects due to the mixing of thermal and relativistic plasma, especially towards the cluster centres.

However, synchrotron emission from haloes requires a population of relativistic (GeV energy range) electrons and a cluster-wide magnetic field. The origins of both particles and field have been a subject of much debate since the late 1970s (e.g. Jaffe 1977; Dennison 1980; Roland 1981). Recent evidence on this question has come through the observation that radio haloes are preferentially found in high X-ray luminosity clusters (Giovannini et al. 1999), and that there may be a steep correlation between the radio power of haloes and the ICM temperature (e.g. Liang et al. 2000). No upper limits so far contradict this correlation in the temperature range 5–15 keV (Liang 2000). The relativistic particles cannot simply leak out of radio sources in the cluster, because the diffusion path length ($`10`$ kpc assuming Alfvén speed propagation) is short compared to the size of clusters ($`>1`$ Mpc). Giovannini & Feretti (2000) found that while diffuse radio sources were found in bright X-ray clusters, none of the 11 clusters selected for the presence of a tailed radio source within the central 300 kpc showed any diffuse radio emission. The presence of diffuse radio emission is therefore more closely connected with the properties of the thermal plasma in the ICM than the presence of tailed radio sources which were thought to provide relativistic particles. This gives observational evidence that particles escaping from radio sources in the cluster are unlikely to be accelerated to the relativistic energies required for the radio halo emission. Instead, we interpret the correlation between halo radio power and ICM temperature, and the similarity in appearance of the X-ray and radio haloes, as showing that the relativistic electrons may be accelerated from the thermal pool. The direct cause of this acceleration may be the merging of subclusters and turbulence in the ICM (e.g. Tribble 1993, Brunetti et al. 2001).

The idea that the same relativistic electrons generate radio emission through their synchrotron losses and X-ray or gamma-ray emission by inverse Compton (IC) scattering was first used by Felten and Morrison (1966) for the interpretation of the galactic diffuse emission. Later, Cooke et al. (1978) used it for the interpretation of the X-ray emission from the lobes of the radio galaxy Centaurus A, assuming that the electrons are scattered by the cosmic microwave background (CMB) photons. A similar model was used by Rephaeli (1979) for the interpretation of the hard X-ray emission from the Coma cluster: the relativistic electrons responsible for the radio halo emission can scatter the CMB photons up to X-ray energies through IC scattering. Efforts have been made to detect this excess of non-thermal emission above the thermal emission from the hot ICM. It is easiest to detect at hard X-ray energies (above $`20`$ keV) where the thermal spectrum is expected to cut off exponentially. Recent observations from Beppo-SAX have shown that such a hard X-ray excess does exist in some clusters, e.g. Coma and A2256 (Fusco-Femiano et al. 1999, 2000).

However, alternative theories of the origin of a hard excess have been proposed by various authors. Enßlin et al. (1999) first suggested that the hard X-ray excess observed in the Coma cluster may be the result of bremsstrahlung of supra-thermal electrons in the ICM accelerated by turbulence. Both Enßlin et al. (1999) and Sarazin & Kempner (2000) have assumed the acceleration process works on particles with energies of a few times $`kT_x`$ (i.e. a few tens of keV), and that the electron energy distribution has a sharp transition from a Maxwellian to a power law at a few times $`kT_x`$. Dogiel (2000) pointed out that the acceleration of particles is only efficient above $`900`$ keV (i.e. the power law distribution starts at $`900`$ keV), but the distortion of the Maxwellian is already significant at much lower energies (as low as $`30`$ keV) due to the flux of thermal particles running away into the region of acceleration (for details see Gurevich 1960). It is the ‘quasi-thermal’ particles just above $`30`$ keV that produce the observed hard X-ray excess through bremsstrahlung emission. However, our analysis of the in situ acceleration model given below (Sec. 6) shows that the spectrum of MHD turbulence needed to produce the radio emitting electrons is very steep which may make this model problematic.

On the other hand, a model where the X-ray and radio emitting particles are generated by electrons from an ‘external’ injected population cannot be ruled out. A possible model could be that particles in the intracluster medium are accelerated only in regions with high level of turbulence, for example near collisionless bow shocks generated by motions of hypothetical dark matter haloes (e.g., Bykov et al. 2000b). Radio observations have not shown significant spatial fluctuations of intensity in radio haloes (e.g. Giovannini et al. 1993, Liang et al. 2000) which means that we must have either in situ acceleration of these electrons everywhere in the intracluster medium or injection by numerous sources whose separation is smaller than the path length of the electrons ($`10`$ kpc) or the length corresponding to the resolution of the radio maps. In the case of the Coma cluster which is relatively bright and close by, the highest resolution VLA map of the radio halo corresponds to a length scale similar to the diffusion path length of the elecrons. Therefore, the injection model requires a high density of localised ‘external’ sources.

Below we investigate both the local and global acceleration models and determine formally the model parameters required to fit the whole observed spectrum from radio to X-rays for the Coma cluster, since this is the best-studied cluster in all wave-bands.

Our analysis is performed in the framework of a model where spatial variations of the model parameters are neglected. There are several reasons for this: a) observations show more or less uniform distribution of non-thermal emission in the halo; and b) we do not have enough information about spatial variations of the hard X-ray flux from the Coma halo, therefore we derive averaged parameters of particle acceleration based on the total flux of hard X-rays.

The cosmological parameters adopted are $`H_0=50`$ km s-1 Mpc-1, $`q_0=0.5`$, so that for the Coma cluster 1 arcmin is 40 kpc.

## 2 The Coma Cluster

Extensive data from radio to X-ray bands exist for the Coma cluster ($`z0.0232`$). The best estimates of the ICM temperature are $`T_x=8.21\pm 0.09`$ keV from the Ginga satellite (Hughes et al. 1993) and $`T_x=8.25\pm 0.10`$ keV from the XMM-Newton satellite (Arnaud et al. 2001). The new high resolution XMM-Newton data have shown that at least the central 1 Mpc of the cluster is isothermal. The gas distribution is well parameterised by a $`\beta `$-model (Cavaliere et al. 1979) with electron density distribution given by

$$n_e(r)=n_0[1+(r/r_0)^2]^{3\beta /2}$$
(1)
where $`\beta =0.75\pm 0.03`$, the core radius $`r_0=10\stackrel{}{.}5\pm 0\stackrel{}{.}6`$ (i.e. $`408h_{50}^1`$ kpc), the emission extends to $`10r_0`$, and the central electron density $`n_0=(2.89\pm 0.04)\times 10^3h_{50}^{1/2}`$ cm-3, from the ROSAT all-sky survey data (Briel et al. 1992).

An excess of hard X-rays above the thermal component has been been detected in the 20-80 keV range by Beppo-SAX satellite \[Fusco-Femiano et al. 1999\]. From Eq. 1, we estimate the average electron density and squared electron density within the Beppo-SAX field of view of $`1\stackrel{}{.}3`$ (i.e. a volume of $`V1.7\times 10^{75}`$ cm3) to be $`\overline{n}_e1.23\times 10^4`$ cm-3 and $`\overline{n_e^2}5.22\times 10^8`$ cm-6. We adopt these values in the discussions that follow, including the calculation of Fig. 1 & 2.

Coma-C was the first radio halo detected, and its emission has been measured between 30 MHz and 4850 MHz (Deiss et al. 1997 and references therein). The spectrum is consistent with a power law of index $`\alpha 1.25`$ between 30 and 1400 MHz. The measurement at 4850 MHz is an upper limit and the apparent steepening at 2700 MHz has been contested (Deiss et al. 1997). Radio spectral index measurements for such diffuse sources are difficult and a mismatch of the beam size or difference in the brightness sensitivity of the observations can often lead to an apparent steepening of the spectral index. Accordingly, we will disregard the measurements above 1400 MHz. The largest extent of the radio halo was found to be $`80\mathrm{}\times 45\mathrm{}`$ on an Effelsberg map at 1.4 GHz (Deiss et al. 1997); by contrast, the 1.4 GHz VLA image which is not as sensitive to the large scale emission showed a maximum extent of only $`22\mathrm{}`$ (Kim et al. 1990).

## 3 Acceleration processes

In most circumstances, the main source of charged particle acceleration in a cosmic plasma is collisions of electrons and ions with magnetic field fluctuations. This leads to a slow, stochastic, energy gain. Here we consider magneto-hydrodynamic turbulence in the strong and weak limits. We can determine which regime is most appropriate in clusters from the spectral indices of radio haloes, since the spectral signatures of the two cases differ.

### 3.1 Weak Turbulence

If the magnetic turbulence is ‘weak’, i.e. the amplitude of magnetic field fluctuations, $`\delta B`$, is much less than the strength of the large scale magnetic field $`B_0`$, $`\delta BB_0`$, then the interaction between the particles and the fluctuations has a resonance character, which has been discussed in the context of the interstellar medium (Berezinskii et al. 1990). The resonance condition is

$$\omega (k)|k_{}|v_{}\omega _H=0.$$
(2)
where the different signs correspond to the different wave polarisations, $`k`$ and $`\omega `$ are the wave-number and the wave frequency respectively, $`v_{}`$ is the component of the particle velocity along the magnetic field, and $`\omega _H`$ is the gyro-frequency of a particle with the total energy $`E_{\mathrm{tot}}`$ and the charge $`Ze`$

$$\omega _H=\frac{ZeB_0c}{E_{\mathrm{tot}}}.$$
(3)
The resonance condition (Eq. 2) can be written in the form

$$k_{}=\frac{ZeB_0}{pc(\mu \omega (k)/kv)},$$
(4)
where $`v`$ is the particle velocity, and $`\mu `$ is the cosine of the particle pitch angle (i.e. the angle between the particle momentum, $`p`$, and the large scale magnetic field, $`B_0`$). For Alfv́en waves with dispersion relation

$$\omega (k)=\pm |k_{}|v_A,$$
(5)
where $`v_A`$ is the Alfvén velocity, and for fast particles ($`vv_A`$) the resonance condition reduces to

$$k_{}=|\frac{ZeB_0}{pc\mu }|=\frac{1}{r_L|\mu |}.$$
(6)
Here $`r_L`$ is the particle gyro-radius. We see that in general electrons and protons of the same energy (especially for non-relativistic energies) are scattered by waves with different wavelength and polarisation. Therefore, in some conditions electrons may be preferentially accelerated. This may be the case in the Coma halo; and in this respect models of hard X-ray emission from Coma which take into account of only the electron acceleration are not implausible (e.g. Schlickeiser et al. 1987).

The process of acceleration may be described as momentum diffusion with a diffusion coefficient determined by the spectrum of magneto-hydrodynamic turbulence. If the spectrum of waves is a power-law

$$W(k)k^q,$$
(7)
the momentum diffusion coefficient depends on the particle momentum as (Berezinskii et al. 1990)

$$a(p)=a_0p^q.$$
(8)
It can be shown that in the simple thick target model the momentum distribution of accelerated particles is given by (see Gurevich 1960)

$$f(p)p^{(q+1)},$$
(9)
where $`f(p)`$ is the electron number density distribution per momentum volume interval, so that the total number of particles above a momentum $`p`$ is given by

$$N(>p)=_p^{\mathrm{}}4\pi f(p)p^2𝑑p,$$
(10)
and the value of the index $`q`$ can be deduced from the spectrum of non-thermal emission produced by the accelerated electrons.

Below (Sec. 6) we perform the necessary analysis in an attempt to estimate $`a_0`$ and $`q`$ and hence to define the mechanism of particle acceleration in the intracluster medium of the Coma cluster.

### 3.2 Strong Turbulence

In the case of strong magneto-hydrodynamic turbulence, where the field fluctuations are of the same order as the mean field, both plasma components, positive and negative, are accelerated. This is seen, for example, in the acceleration of particles in molecular clouds (Dogiel et al. 1987) or acceleration by supersonic turbulence, which effectively acts as an ensemble of shock waves (Bykov and Toptygin 1993). Stochastic acceleration by strong turbulence is described as momentum diffusion with diffusion coefficient $`a(p)`$ of the form (e.g. Toptygin 1985)

$$a(p)=a_0p^2,$$
(11)
which creates a relatively flat spectrum of accelerated particles

$$f(p)p^3.$$
(12)
The interpretation of the hard X-ray emission in Coma in terms of this kind of acceleration was discussed by Dogiel (2000).

A similar acceleration process takes place near shock fronts, where the magnetic field strength jumps by an amount of the order of the mean field. However, the spectrum of accelerated particles near a shock front has

$$f(p)p^4,$$
(13)
steeper than that of an ensemble of shocks (Eq. 12). This case has been analysed for the Coma cluster by Bykov et al. (2000b).

## 4 Particle Kinetic Equation

To describe the distribution function in both non-relativistic and relativistic energy ranges it is more convenient to use a canonical momentum variable, $`p`$, than the particle kinetic energy $`E`$. The electron number density distribution per unit energy interval, $`F(E)`$, is related to the momentum distribution per unit momentum volume, $`f(p)`$, by

$$F(E)=4\pi f(p)p^2\frac{dp}{dE},$$
(14)
where

$$\frac{dp}{dE}=\frac{1}{p\sqrt{m_ec^2kT_x}}\sqrt{p^2+\frac{m_ec^2}{kT_x}},$$
(15)
where $`m_e`$ is the electron rest mass, $`T_x`$ is the temperature of the thermal gas, $`k`$ is the Boltzman constant, and $`p`$ is the dimensionless particle momentum, $`p=\widehat{p}/\sqrt{m_ekT_x}`$ where $`\widehat{p}`$ is the dimensional momentum.

The generalised Fokker-Planck approximation to the kinetic equation, describing the particle distribution function $`f`$ in the intracluster medium (see Berezinskii et al. 1990), written in dimensionless variables is

$`{\displaystyle \frac{f}{t}}+u(𝐫){\displaystyle \frac{f}{𝐫}}{\displaystyle \frac{p}{3}}{\displaystyle \frac{u(𝐫)}{𝐫}}{\displaystyle \frac{f}{p}}Q(p,𝐫)=`$ (16)
$`{\displaystyle \frac{1}{p^2}}{\displaystyle \frac{}{p}}\left[𝒜(p){\displaystyle \frac{f(p)}{p}}+(p)f\right]+(\kappa (𝐫,p)f),`$
where $`u(𝐫)`$ is the particle convection, $`\kappa `$ is the spatial diffusion coefficient, and $`Q`$ describes the distribution of sources emitting fast particles into the intracluster medium. We define non-dimensional variables $`t=\widehat{t}\nu _0`$, $`r=\widehat{r}/L`$, where $`\widehat{t}`$ and $`\widehat{r}`$ are the dimensional variables, $`L`$ is the size of the cluster, and we define the characteristic frequency $`\nu _0`$ as

$$\nu _0=\frac{2\pi \overline{n}_ec^2r_e^2m_e}{\sqrt{m_ekT_x}},$$
(17)
$`\overline{n}_e`$ is the average electron number density, and $`r_e`$ is the classical electron radius $`r_e=e^2/m_ec^2`$. For the Coma halo of $`\overline{n}_e=1.23\times 10^4`$ cm-3 within a volume of $`1.7\times 10^{75}`$ cm-3 and $`T_x=8.2`$ keV (Sec. 2), $`\nu _06.4\times 10^{17}`$ s-1.

The frequency $`\nu _0`$ is proportional to the frequency $`\nu _1`$ of thermal particle Coulomb collisions at kinetic energy $`kT_x`$

$$\nu _1=2\nu _0\left(\frac{m_ec^2}{kT_x}\right)\mathrm{ln}\mathrm{\Lambda },$$
(18)
where $`\mathrm{ln}\mathrm{\Lambda }`$ is the Coulomb logarithm. The frequency $`\nu _15.3\times 10^{14}`$ s-1 for the Coma halo.

In Eq. 16, $``$ describes the particle energy gain and loss and $`𝒜`$ the momentum diffusion due to Coulomb collisions with background particles and interactions with the magnetic fluctuations,

$$(p)=p^2\left[\left(\frac{dp}{dt}\right)_{ion}\left(\frac{dp}{dt}\right)_{synIC}\left(\frac{dp}{dt}\right)_{brem}\right],$$
(19)
$$𝒜(p)=p^2\left[\left(\frac{dp}{dt}\right)_{ion}\frac{\gamma }{\sqrt{\gamma ^21}}\sqrt{\frac{kT_x}{m_ec^2}}+a_0p^q\right],$$
(20)
where $`\gamma =1+E(p)/m_ec^2`$ is the Lorentz factor, and the dimensionless parameter $`a_0`$ is related to the dimensional parameter $`\widehat{a}_0`$ as

$$a_0=\frac{\widehat{a}_0(\sqrt{m_ekT_x})^{q2}}{\nu _0}.$$
(21)
As one can see from Eq. 20 the process of momentum diffusion is determined by Coulomb collisions with background particles (the first rhs term of Eq. 20) and collisions with electromagnetic fluctuations (the second rhs term of Eq 20).

The $`dp/dt`$ terms are rates of momentum loss due to ionisation, synchrotron plus inverse Compton and bremsstrahlung emission respectively (e.g. Hayakawa 1969). The ionisation loss term is

$`\left({\displaystyle \frac{dp}{dt}}\right)_{ion}={\displaystyle \frac{1}{p}}\sqrt{p^2+{\displaystyle \frac{m_ec^2}{kT_x}}}{\displaystyle \frac{\gamma }{\sqrt{\gamma ^21}}}`$ (22)
$`\times \left[\mathrm{log}\left({\displaystyle \frac{E(p)m_ec^2(\gamma ^21)}{h^2\omega _p^2\gamma ^2}}\right)+0.43\right],`$
where the plasma frequency $`\omega _p=\sqrt{4\pi e^2\overline{n}_e/m_e}`$. The combined synchrotron and IC loss term is

$`\left({\displaystyle \frac{dp}{dt}}\right)_{synIC}={\displaystyle \frac{1}{p}}\sqrt{p^2+{\displaystyle \frac{m_ec^2}{kT_x}}}{\displaystyle \frac{32\pi r_e^2}{9\sqrt{m_ekT_x}}}{\displaystyle \frac{1}{\nu _0}}`$ (23)
$`\times \left[U_{cmb}+U_{mag}\right]\left({\displaystyle \frac{E(p)}{m_ec^2}}\right)^2,`$
where $`U_{cmb}`$ is the energy density of the CMB, and $`U_{mag}`$ is the magnetic energy density. The bremsstrahlung loss term is given by

$`\left({\displaystyle \frac{dp}{dt}}\right)_{brem}={\displaystyle \frac{1}{p}}\sqrt{p^2+{\displaystyle \frac{m_ec^2}{kT_x}}}{\displaystyle \frac{\overline{n}_em_ec^2e^2r_e^2}{hc\sqrt{m_ekT_x}}}{\displaystyle \frac{1}{\nu _0}}`$ (24)
$`\times \left[\mathrm{log}{\displaystyle \frac{2(E(p)+m_ec^2)}{m_ec^2}}{\displaystyle \frac{1}{3}}\right]{\displaystyle \frac{E(p)}{m_ec^2}}.`$
In comparison to the loss terms, the average rate of acceleration due to momentum diffusion (Eq. 8) is

$`\left({\displaystyle \frac{dp}{dt}}\right)_{acc}a_0p^{q1}.`$ (25)
These loss terms and the acceleration term are plotted as a function of particle energy in Fig. 1 for the Coma cluster. At low energies ionisation loss dominates, and at high energies either the acceleration term or the combined synchrotron and IC term dominates depending on the value of $`q`$.

Equation 16 describes the particle spectrum in the energy range from thermal energies, where the spectrum is Maxwellian, up to non-thermal energies, as well as the particle spatial distribution. In general it is difficult to obtain a solution, though in some special cases it is possible to solve it analytically or numerically (see Bulanov and Dogiel 1979).

Fortunately, in many cases the spatial distribution is unimportant. Whether or not a model without spatial terms can be applied to Coma is rather questionable. At the moment we do not have enough information to prove or disprove this.

If we ignore the spatial variation, then we can integrate Eq. 16 over the volume of the object and the terms describing the spatial distribution $`f`$ converge to a single term, which can be represented as $`f/t_{\mathrm{eff}}`$ where $`t_{\mathrm{eff}}`$ is the effective escape time of the particles due to their spatial propagation. If the characteristic escape time is much greater than the energy loss timescale, we can neglect the term $`f/t_{\mathrm{eff}}`$. Then Eq. 16, rewritten for the spatially-averaged distribution function $`f(p)`$, can be presented as

$$\frac{f}{t}\frac{1}{p^2}\frac{}{p}\left[𝒜(p)\frac{f(p)}{p}+(p)f\right]=Q(p).$$
(26)
Two possible fast particle production models are analysed for the Coma halo in Sec. 6 and 7:

* in situ acceleration everywhere in the halo, where the spectrum of non-thermal particles is established by acceleration processes represented in Eq. 26 by the momentum diffusion term $`𝒜(p)`$ and balanced by the loss term, $`(p)`$; and
* injection into the intracluster medium from many relatively small regions where acceleration takes place (perhaps strong shocks, or regions with high levels of magnetic turbulence). These regions are supposed to be more or less uniformly distributed in space with no acceleration between them. The spectrum in the intracluster medium is then fully determined by energy loss processes and the injection spectrum, i.e. $`(p)`$ and $`Q(p)`$.

## 5 Processes of Non-Thermal Emission Generated by Fast Electrons

To determine the parameters for both models we need to calculate the X-ray and radio spectra expected from the non-thermal electrons and compare our results with observational data. In this section, we collect the equations describing the radiation processes.

### 5.1 Non-thermal Bremsstrahlung Emission

The bremsstrahlung flux density $`F_{brem}`$ in units of ph cm-2 s-1 keV-1 is given by

$$F_{brem}=\frac{V}{4\pi d_L^2}_{E_\gamma }^{\mathrm{}}\overline{n_e^2}\frac{F(E)}{\overline{n}_e}\sqrt{\frac{2E}{m_e}}\frac{d\sigma _x}{dE_\gamma }𝑑E$$
(27)
where $`V`$ is the emission volume, $`d_L`$ is the luminosity distance, $`E_\gamma `$ is the photon energy (not to be confused with the particle energy $`E`$), and the emission cross section is given by

$$\frac{d\sigma _x}{dE_\gamma }=\frac{16}{3}\frac{e^2}{\mathrm{}c}\frac{r_e^2m_ec^2}{EE_\gamma }\mathrm{log}\left(\frac{\sqrt{E}+\sqrt{EE_\gamma }}{\sqrt{E_\gamma }}\right)$$
(28)
\[Hayakawa 1969\]. Non-thermal bremsstrahlung radiation in the hard X-ray range comes from electrons in the energy range $`30100`$ keV.

### 5.2 Inverse Compton Emission

Here we calculate the IC flux density $`F_{IC}`$ in units of ph cm-2 s-1 keV-1 due to the scattering of CMB photons by relativistic particles. Since we are only interested in the IC scattering above 100 eV (i.e. EUV and above), we are interested only in electrons with $`E>10^5`$ keV. We will show (Fig. 2) that the electron energy distribution follows a power law with slope $`1q`$ for $`E>10^4`$ keV. We can use the analytic solution for a power law distribution in the Thomson limit ($`\sqrt{E_\gamma kT_r}m_ec^2`$) given in Blumenthal & Gould (1970; Eq. 2.65) for this power law to find

$$F_{IC}=\frac{V}{4\pi d_L^2}\frac{r_e^2}{\pi h^3c^2}K_e(m_ec^2)^{2q}(kT_r)^{\frac{q+4}{2}}\zeta (q1)E_\gamma ^{q/2},$$
(29)
where numerical values of $`\zeta (q1)`$ are given in Table 1 of Blumenthal & Gould (1970), $`T_r`$ is the temperature of the CMB (not to be confused with the electron temperature $`T_x`$), $`K_e=F(E_{min})E_{min}/(q2)`$ is the normalisation for the power law distribution and $`E_{min}`$ is typically $`10^4`$ keV. Electrons of energies $`5`$ GeV inverse Compton scatter photons of the CMB to hard X-ray energies $`60`$ keV.

### 5.3 Synchrotron Radio Emission

We calculate the synchrotron flux density in Jy from the electron energy distribution using Ginzburg & Syrovatskii (1964):

$$F_{radio}=\frac{V}{4\pi d_L^2}10^{23}_{E_{min}}^{E_{max}}\sqrt{3}r_eeBF(E)f_x𝑑E,$$
(30)
where

$$f_x=\frac{1}{2}_0^\pi \mathrm{sin}^2\theta \frac{\nu }{\nu _c}_{\frac{\nu }{\nu _c}}^{\mathrm{}}K_{5/3}(\eta )𝑑\eta 𝑑\theta ,$$
(31)
$`\nu _c=4.21B\mathrm{sin}\theta (E/m_ec^2)^2`$ MHz, and $`B`$ is in units of Gauss. Since in our case the electron energy distribution is a power law for $`E>10^4`$ keV and the synchrotron emission in the radio frequency range is given by the GeV electrons, we can safely use the analytic equation for a power law electron energy distribution given in Eq. 6.10 of Ginzburg & Syrovatskii (1964). For radio emission above $`10`$MHz, we need electrons above energies of $`1`$ GeV for $`B1\mu `$G.

## 6 Model of in situ Acceleration

If the injection term, $`Q(p)=0`$, then apart from the ionisation and radiation losses, the spectrum is completely determined by two processes: Coulomb collisions, which form the thermal and quasi-thermal part of the spectrum, and stochastic acceleration, which forms the spectrum at high energies (see Fig. 1).

A quasi-stationary solution to Eq. 26, in the sense that the presence of a ‘weak’ acceleration mechanism leads to the continuous acceleration of particles above a certain injection energy, is given by Gurevich (1960):

$$f(p)=\sqrt{\frac{2}{\pi }}\overline{n}_e\mathrm{exp}\left(_0^p\frac{(v)}{𝒜(v)}𝑑v\right)𝒢,$$
(32)
where

$$𝒢=\frac{_p^{\mathrm{}}\frac{dv}{𝒜(v)}\mathrm{exp}(_0^v\frac{(t)}{𝒜(t)}𝑑t)}{_0^{\mathrm{}}\frac{dv}{A(v)}\mathrm{exp}(_0^v\frac{(t)}{𝒜(t)}𝑑t)},$$
(33)
The energy distribution of the electron number density thus calculated is shown in Fig. 2 for 3 sets of ($`a_0`$,$`q`$). It can be seen that $`F(E)`$ starts to deviate from a pure Maxwellian at a few tens of keV and reaches a power law at $`110`$ MeV.

Radio spectra implied by these ($`a_0`$,q) parameter sets are shown in Fig. 3, and the hard X-ray spectra are shown in Fig. 4. From Eq. 32, it follows that the momentum spectral index of the accelerated particles is $`q1`$, hence $`q`$ is related to the the spectral index $`\alpha `$ of the synchrotron emission as $`q=2\alpha +2`$. The observed radio spectrum (Fig. 3) constrains the electron momentum spectral index to $`q4.5_{0.25}^{+0.2}`$. Figure 5 shows that the best fit acceleration parameters from the radio and hard X-ray data are $`a_00.43_{0.18}^{+0.27}`$, $`q4.50.2`$, where the errors correspond to 95% confidence. The corresponding magnetic field derived from adjusting the model flux density to fit the data at 326 MHz, is $`B1.51_{0.95}^{+2.31}\mu `$G. The magnetic field is well constrained in this model (with only 2 free parameters) by the X-ray and radio data, and is consistent with the limit deduced from Faraday rotation measurements, $`B5.9h_{50}^{1/2}\mu `$G on field tangling scales of $`<1`$ kpc (Feretti et al. 1995). In comparison, the equipartition magnetic field is $`B1`$$`\mu `$G assuming that there are equal number of relativistic protons and electrons and that the emission filling factor is 1. Any independent estimates of the cluster global magnetic field would provide a test for the validity of the model.

In the framework of this model the origin of the hard non-thermal X-rays is bremsstrahlung since the flux of IC X-ray photons is negligible. Figure 7 shows the expected ratio of the bremsstrahlung to IC flux density for the Coma cluster at 70 keV as a function of the momentum diffusion coefficient spectral index, $`q`$. If the GeV electrons responsible for the radio halo emission are the result of in-situ acceleration, then the observed radio halo spectrum ($`\alpha 1.25`$ or $`q4.5`$) implies that the contribution from the bremsstrahlung of 30-100 keV electrons is always greater than that of the IC scattering of CMB photons by the GeV electrons. This result is independent of the magnetic field, since neither the bremsstrahlung nor the IC emissivity depends on the magnetic field.

CMB photons can also be up-scattered to the extreme-ultraviolet (EUV) energy range by electrons of energy $`2.5\times 10^5`$ keV through inverse Compton scattering. The detection of excess EUV emission above that of thermal bremsstrahlung has been reported in the Coma cluster using data gathered by the Extreme Ultraviolet Explorer (EUVE; Lieu et al. 1999 and references therein). However, this result is controversial (see e.g. Bowyer et al. 1999), hence we have chosen to predict the excess EUV flux density in the framework of our model instead of fitting it to the data. We estimate the expected excess EUV flux density at 100 eV to be $`1.9\times 10^3`$ ph cm-2 s-1 keV-1 within a $`15^{^{}}`$ radius which is nearly 2 orders of magnitude below the EUVE sensitivity. However, given the controversy over the observational data, it can not be used to rule out the model without further independent measurements.

It is interesting to note that, for $`q<4`$, IC starts to dominate over bremsstrahlung in the hard X-ray regime, but for $`q2`$ bremsstrahlung dominates again because synchrotron and IC losses exceed the gains from acceleration and causes the electron spectrum to cut off at $`1`$ GeV. To produce X-ray emission above 20 keV through IC scattering of the CMB, we need electrons with energies above $`3.5`$ GeV. For example, if the radiating particles are accelerated by strong turbulence, that is $`q=2`$ (hence $`a_0=5`$ to fit the hard X-ray data assuming a bremsstrahlung origin), then ionisation losses dominate over acceleration below 1 MeV (see Fig. 1). The gain from acceleration overwhelms particle energy losses in the energy range of 1 MeV - 1 GeV. Hence, the maximum energy of accelerated particles is $`1`$ GeV. At higher energies, the combined synchrotron and IC losses dominates over acceleration (see Fig. 1), hence strong turbulence is unable to produce radio emitting GeV electrons. Acceleration by strong turbulence can produce only relatively low energy electrons emitting in the energy range of hard X-rays to soft gamma-rays through bremsstrahlung. Another mechanism for the production of radio-emitting electrons is then necessary to explain the radio halo in Coma. This conclusion is somewhat in contradiction with the results of Dogiel (2000) who showed that in this case strong turbulence can accelerate particles up to energies of $`10`$ GeV, albeit with a much flatter spectrum than needed to produce the observed radio spectrum. The reason for this disagreement is that Dogiel (2000) used the density for central regions of the cluster, where the frequency of thermal particle collisions is high. This led to a higher acceleration efficiency than in our present treatment which uses average cluster parameters.

It is shown in Fig. 1 that for $`q<2`$, the synchrotron and IC losses win over the acceleration term at energies even lower than 1 GeV. This explains why Blasi (2000), in his study of stochastic acceleration in the ICM assuming a Kolmogorov spectrum ($`q=5/3`$) concluded that in situ acceleration could not explain the radio emission because of the synchrotron and IC loss cutoff. However, we have shown that the value $`q=4.5`$, corresponding to weak turbulence (Sec. 3.2), is consistent with the radio data (Fig. 3).

In the case of $`q>3`$, however, the acceleration term always wins over the loss terms at high energies (Fig. 1), which means a rapid decrease of the characteristic acceleration time, $`\tau _{acc}`$, with particle momentum

$$\tau _{acc}\frac{\widehat{p}^2}{\widehat{a}(\widehat{p})}=\frac{p^{2q}}{a_0\nu _0}$$
(34)
and a larger $`q`$ results in a smaller value of the injection energy at which processes of acceleration dominate. The best fit parameters, $`a_00.43`$ and $`q4.5`$, imply a cluster average $`\tau _{acc}7\times 10^6`$ yrs at the injection energy ($`200`$ keV), which means that for $`E>200`$ keV, $`\tau _{acc}`$ is much less than the cluster lifetime. Hence, the model is self-consistent in the sense that a ‘quasi-stationary’ state can easily be reached within a cluster lifetime.

We note, however, the steep spectrum of turbulence ($`q4.5`$) required for the radio data may pose a problem for the model, since it is known that turbulence with $`q>2`$ cannot effectively scatter particles at pitch angles near $`\pi /2`$ in the quasi-linear regime. Nonlinear processes have to be involved to cause effective acceleration(e.g., Berezinskii et al. 1990). Further consideration of this point is beyond the scope of our model.

Cluster thermal X-ray and radio structures are similar in morphology (Deiss et al. 1997, Feretti 1999, Liang et al. 2000, Govoni et al. 2001). This would be expected from the in situ model. We can calculate the expected relationship between the radio halo and X-ray profile, if the synchrotron-emitting relativistic particles are drawn from the thermal pool through in-situ acceleration. If the magnetic field is produced by a dynamo process (e.g. turbulence), then magnetic energy is in equipartition with the energy of turbulent motions (Kraichnan & Nagarajan 1967, De Young 1980). Assuming that the energy of turbulence is proportional to the thermal energy, then we have the magnetic energy density in a cluster medium proportional to the thermal energy density $`P_{mag}P_{th}`$, so that

$$Bn_e^{1/2}T_x^{1/2}.$$
(35)
As pointed out by Enßlin et al. (1998), this implies that the radio halo surface brightness profile should be slightly steeper than the X-ray profile. The X-ray surface brightness is given by the projection of $`n_e^2`$, which for a $`\beta `$model with $`n_e`$ given by Eq. 1 is

$$I_x(r)\left[1+(\frac{r}{r_0})^2\right]^{3\beta +1/2}$$
(36)
for fixed $`T_x`$.

Since relativistic particles are accelerated from the thermal pool, we can assume that their density is proportional to the thermal particle density $`n_e`$. Then the radio halo brightness profile is given by

$$I_{1.4\mathrm{GHz}}n_eB^{1+\alpha }𝑑l,$$
(37)
when using Eq. 35, and assuming an isothermal ICM of the form of Eq. 1

$$I_{1.4\mathrm{GHz}}\left[1+(\frac{r}{r_0})^2\right]^{3\beta (3+\alpha )/4+1/2}.$$
(38)
If the radio spectral index is exactly $`\alpha =1`$, then the radio and X-ray profiles follow each other. In the case of Coma and 1E0657-56, the spectral index is $`\alpha 1.25`$ and the radio profile is slightly steeper than the X-ray profile. For the clusters studied by Govoni et al. (2001), the radio and X-ray surface brightness profiles are very similar in shape, implying $`\alpha 1`$.

## 7 Energy Output of Non-Thermal and Quasi-Thermal Particles

In this section, for clarity we need to distinguish particles from the 3 different energy ranges: thermal, quasi-thermal and non-thermal. The Maxwellian thermal particle spectrum is shaped by Coulomb collisions, where each thermal particle collides with other background particles, exchanging energies without changing the shape of the energy spectrum. The lifetime of the thermal particle distribution is infinite if we neglect thermal bremsstrahlung emission. On the other hand, a non-thermal particle spectrum is formed by acceleration processes. If these particles are injected into a region with no further acceleration, the particles would lose their energy to the background particles through Coulomb collisions, thus heating the background plasma and become thermal particles themselves. The lifetime of non-thermal particles is equal to that of the individual particles. Unlike the thermal and non-thermal particles, the quasi-thermal particles exist only in regions of active particle acceleration. Like the thermal particles, the quasi-thermal (or suprathermal) particle distribution is formed by Coulomb collisions, but unlike the thermal particle distribution it is not in strict equilibrium (but in quasi-equilibrium) and the quasi-thermal particle distribution would dissipate like the non-thermal particles once the acceleration process is switched off. Coulomb collisions cause a net run-away of thermal particles into regions of acceleration in the spectrum, resulting in a spectral distortion with an excess above the thermal equilibrium distribution in the quasi-thermal energy range. Hence the life-time of quasi-thermal particles is shorter than the thermal particles, but longer than the non-thermal particles.

The energy output of emitting particles is determined by their energy losses: ionisation losses (Coulomb collisions with background particles) for quasi-thermal particles, and inverse Compton scattering and synchrotron radiation losses for non-thermal electrons. The necessary energy output, $`L_{cr}`$, can roughly be estimated from the observed value of X-ray luminosity $`L_x`$ as

$$L_{cr}L_x\frac{\tau _x}{\tau _{cr}},$$
(39)
where $`\tau _x`$ is the characteristic time of X-ray emission by parent particles and $`\tau _{cr}`$ is the life-time of these particles.

If Coma cluster emits X-rays by inverse Compton process, then from Eq. 39 it follows that

$$L_{cr}^CL_x^C410^{43}\text{erg s}\text{-1},$$
(40)
where $`L_x^C410^{43}`$ erg s-1 is the luminosity of hard ($`2080`$ keV) X-rays from the Coma halo. The uncertainty in the estimate $`L_{cr}^C`$ is determined by the unknown rate of synchrotron losses, since the magnetic field strength is only known within fairly wide bounds (a range that spans nearly 2 orders of magnitude).

In the case of a bremsstrahlung origin of the Coma luminosity produced by non-thermal particles, as in models of Enßlin et al. (1999) and Sarazin & Kempner (2000), the necessary energy output $`L_{cr}^C`$ is much larger than $`L_x^C`$. The reason is that these non-thermal particles lose their energy by collisions with background electrons and protons (ionisation losses) and the probability of their generating a bremsstrahlung photon is very low ($`<2.5\times 10^4`$ for photons with $`E_x<100`$ keV; Valinia and Marshall 1998). Therefore the necessary energy output in this case is huge

$$L_{cr}^C10^{48}\text{erg s}\text{-1},$$
(41)
and as pointed out by Petrosian (2001), this makes the model problematic.

Similar problems arise in the interpretation of the hard X-ray flux from the galactic ridge. The spectrum of hard X-ray emission from the galactic ridge also shows an excess above the thermal flux. However, in this case its origin cannot be interpreted as due to inverse Compton of relativistic electrons (Skibo et al. 1996), thermal emission of hot plasma or the integrated contribution of discrete unresolved source (Tanaka 2002). The only possibility to explain the origin of hard X-ray flux from the galactic disk is just bremsstrahlung radiation of sub-relativistic galactic electrons. If these electrons are non-thermal the necessary energy output should be of the order of $`10^{43}`$ erg s-1, i.e. more than can be supplied by supernovae, the most probable sources of energy in the Galaxy (see e.g. Valinia et al. 2000, and Dogiel et al. 2002a). It was shown, however by Dogiel et al. (2002b) that this energy problem can be resolved if the emission is produced by quasi-thermal electrons. The life-time of quasi-thermal electrons is longer than that of non-thermal electrons. Therefore, as follows from Eq. 39, the necessary energy output of quasi-thermal electrons is smaller than the energy output of non-thermal electrons. We showed above that in the case of in-situ acceleration, unlike the models of Enßlin et al. (1999) and Sarazin & Kempner (2000) where non-thermal electrons with power-law spectrum were responsible for the excess hard X-ray in the Coma cluster, it is the quasi-thermal electrons with exponential distribution that give rise to the excess hard X-ray emission. If our model is correct, the energy output in quasi-thermal electrons needed to explain the Coma X-ray flux is of the order of $`L_{cr}^C10^{45}`$ erg s-1 as follows from numerical calculations for the parameters $`a_0=0.43`$ and $`q=4.5`$, i.e. much less than the estimate in Eq. 41.

The bremsstrahlung X-ray flux can also be generated by sub-relativistic protons (see Dogiel 2001). Though we do not know which specific process is responsible for particle acceleration in Coma, we note that in the case of Fermi acceleration the characteristic acceleration frequency is proportional to $`\alpha _0v/L`$ (see Toptigyn (1985)) where $`v`$ is the particle velocity and $`L`$ is the particle mean free path (an average distance between scattering centres). Since velocities of thermal electrons are $`\sqrt{M/m}`$ times larger than velocities of thermal protons, we expect that the rate of acceleration for electrons $`\alpha _0^e\alpha _0^p`$, hence background electrons may be accelerated to sub-relativistic energies more effectively than protons.

## 8 Model of ‘External’ Sources of Acceleration

As follows from Fig. 1 the particle spectrum is controlled by ionisation losses at low energies and synchrotron$`+`$IC losses at high energies, which balance the source injection spectrum and acceleration. If we assume no in situ acceleration, i.e. $`𝒜(p)=0`$ in Eq. 26, then the kinetic equation for the distribution function $`f`$ can be written as:

$$\frac{f(p)}{t}\frac{1}{p^2}\frac{}{p}(p)f(p)=Q(p).$$
(42)
The stationary solution of this equation is given by

$$f(p)=\frac{1}{(p)}p^2Q(p)𝑑p$$
(43)
i.e.

$$f(p)=\frac{K}{(p)(\eta 1)}p^{3\eta }$$
(44)
where we assume a power law injection spectrum $`Q(p)=Kp^\eta `$, or $`Q(E)E^{\gamma _0}`$ (where $`\gamma _0=\eta 2`$ as follows from Eq. 14).

The best fit to the radio spectrum gives $`\eta 4.5`$, or $`\gamma _02.5`$ (Figure 6), just as is expected for the model of particle acceleration by strong shocks (Biermann et al. 1993). In this model, the origin of the hard X-ray emission is definitely IC with no significant bremsstrahlung (Fig. 7). The model reproduces the emission spectrum from hard X-rays to radio if the strength of the magnetic field is less than $`0.10.2\mu `$G. Thus the finding from Faraday rotation measurements (e.g. Feretti et al. 1995; Clarke et al. 2001) of cluster fields of a few $`\mu `$G argues against the IC model for the hard X-ray (e.g. Enßlin et al. 1999 and Dogiel 2000). However, this is not necessarily a problem. First, since there are only a handful of radio sources behind the cluster with sufficient polarised flux to be detected, the regions of a cluster probed by these sources may not be representative of the average field of the entire cluster. Second, since as many as half of the polarised radio sources used to deduce the cluster magnetic field are cluster members rather than background sources, we expect that such measurements would be biased, as a higher $`B`$-field is likely near a radio source because of its perturbing effect on the ICM. Third, Feretti et al. (1995) in their study of the tangling scale of the magnetic field in the Coma cluster, using an extended polarised radio source in the centre of the cluster, suggested that the cluster magnetic field in Coma has a small scale ($`<1`$ kpc) component of strength $`5.9h_{50}^{1/2}\mu `$G as well as a uniform large scale ($`200`$ kpc) component of strength $`0.10.3\mu `$G which is consistent with the magnetic field of $`0.15\mu `$G deduced from the IC model.

Figure 7 shows that for this model IC emission dominates bremsstrahlung at 70 keV for the required steep spectral index, $`\eta 4.5`$. The presence of strong shocks in galaxy clusters shown by spatial and temperature images from ROSAT, ASCA and now XMM-Newton and Chandra (see Donnelly et al. 1999, Markevitch et al. 1999, Arnaud et al. 2001, Sun et al. 2001 and Boute 2001) lends observational support to this model.

## 9 Conclusions

We have shown that it is possible to explain both the excess hard X-ray data and the radio halo spectrum with a model where the relativistic electrons are weakly accelerated from the thermal pool by turbulence in the ICM. By using a quasi-stationary solution of the kinetic equation, we arrive at an electron energy distribution that deviates from a Maxwellian at $`30`$ keV and reaches a power law at $`110`$ MeV. In the process of reaching an equilibrium, the collisions between thermal particles produce the supra-thermal electrons at $`30100`$ keV required to give the excess hard X-ray emission through bremsstrahlung emission. The accelerated GeV electrons are then sufficient to produce the observed radio halo spectrum with a magnetic field close to the equipartition value. The model produces a self-consistent electron energy distribution and has only two free parameters, the acceleration parameter $`a_0`$ and the momentum spectral index $`q`$, which are well constrained by the observed excess hard X-ray spectrum and the radio halo spectral index (Fig. 5). In the hard X-ray energy range, the contribution from the inverse Compton scattering of the GeV electrons is negligible compared with the non-thermal bremsstrahlung emission for steep radio spectral indices $`\alpha >1`$.

It is also possible to explain the hard X-ray and radio spectrum by a model where accelerated particles are externally injected in numerous localised regions distributed throughout the cluster. In this case, the contribution from IC scattering dominates that of non-thermal bremsstrahlung in the hard X-ray energy range for steep radio spectral indices.

For both the above models a strong assumption was made, that the non-thermal X-ray and radio emission from Coma is generated from a single electron acceleration mechanism. An advantage of this assumption is that it limits significantly the number of free parameters in the model. Of course it is possible to have both in situ acceleration and external injection of accelerated particles (or global and local acceleration) at work in the Coma cluster. We have, however, examined the two extremes rather than such a hybrid model. This study has shown that if it is purely in situ acceleration that produces the particles necessary for the non-thermal emission in the hard X-ray and radio spectral bands, then the only explanation for the origin of the hard X-ray is bremsstrahlung of supra-thermal electrons; and consequently, to fit both X-ray and radio emission a cluster-wide magnetic field of order $`\mu `$G (i.e. close to the equipartition value) is implied. On the other hand, if the accelerated particles are purely provided by numerous, localised external sources of injection, then the only possible origin for the hard X-rays is inverse Compton scattering of CMB photons by relativistic electrons, which implies a cluster-wide magnetic field of 0.1–0.2 $`\mu `$G. Any independent means of measuring the global cluster magnetic fields would distinguish between the two models.

At present it is not easy to choose between the models of particle acceleration for Coma. Both models: in situ acceleration (global) and shock wave acceleration (local) have their advantages and problems. In favour of the in situ model is the similar morphology between the thermal X-ray and the radio emission. The radio data have not shown significant structure in the radio intensity in the Coma halo which would indicate regions of electron acceleration. Analysis of the thermal X-ray emission from the Coma cluster also did not show the strong temperature variations expected from shocks (Arnaud et al. 2001). On the other hand, several features indicating recent and ongoing mergers are found in the Coma cluster (e.g. Donnelly et al. 1999, Arnaud et al. 2001). These mergers often produce multiple large scale shocks and turbulence capable of both magnetic field amplification and in situ electron acceleration (see a recent paper by Ohno et al. 2002 on the subject). Whether or not these regions extend far enough to allow electrons to fill the whole volume of the halo, remains to be discovered. If the origins of the radio and X-ray emitting electrons are different, then a much wider range of possible models can be accommodated. For example, there has been revived interest in the class of secondary electron models where relativistic and thermal protons in the ICM collide to produce relativistic (secondary) electrons (Dennison 1980; Blasi & Colafrancesco 1999; Dolag & Enßlin 2000; and Miniati et al. 2001).

We would like to thank R. Fusco-Femiano for providing us with the data points for the PDS spectrum for Coma. We are grateful to the referee T. Enßlin for helpful comments and suggestions. VAD thanks the University of Bristol and the Institute of Space and Astronautical Science (Sagamihara, Japan) in particular Prof. H. Inoue for hospitality during his visits.


text:
# Untitled Document

Restriction Conditions of Motion Equations of Non-Abelian
Gauge Fields and Their Wave Solutions
Mei Xiaochun

(Department of Physics, Fuzhou University, Fuzhou, 350025, China, E-mail: fzbgk@pub3.fz.fj.cn)

Abstract Similar to introduce the Lorentz condition in the motion equation of electromagnetic field, the restriction condition is introduced in the motion equations of non-Ablian gauge fields so that the equations are simplified greatly and their wave solutions are obtained. Only in this way, the Feynman rule’s use become rightful in the calculation of perturbation theory.

PACS numbers: 1110, 1130

Let’s first review the motion equation of electromagnetic field and its solution. The Lagrangian of free electromagnetic field is

$$L=\frac{1}{4}F_{\mu \nu }(x)F_{\mu \nu }(x)$$
(1)
Here

$$F_{\mu \nu }(x)=_\mu A_\nu (x)_\nu A_\mu (x)$$
(2)
From Eq.(1), the motion equation of free electromagnetic field can be obtained

$$_\mu F_{\mu \nu }=^2A_\nu _\nu _\mu A_\mu =0$$
(3)
On the other hand, Eq.(3) is not a wave equation, so it is improper to be used to describe free electromagnetic field or photon field with wave nature. Free electromagnetic fields or photon fields should satisfy wave equation

$$^2A_\mu =0$$
(4)
In order to get it from Eq.(1), the Lorentz condition

$$_\mu A_\mu =0$$
(5)
is introduced. The spatial solution of Eq.(4) can be written as

$$A_\mu (x)=_{\mathrm{}}^+\mathrm{}a_\mu (\stackrel{}{k})e^{ikx}d^3\stackrel{}{k}$$
(6)
Put it into the Lorentz condition, we get

$$k_\mu a_\mu (\stackrel{}{k})=0$$
(7)
This is a restriction condition about parameters $`a_\mu (\stackrel{}{k})`$. Only three $`a_\mu (\stackrel{}{k})`$ are independent among four.

It is proved below that by introducing similar restriction condition, the motion equations of non-Ablian gauge fields and their solutions can also be simplified. The Lagrangian of non-Ablian gauge fields with zero masses is

$$L=\frac{1}{4}F_{\mu \nu }^\alpha F_{\mu \nu }^\alpha $$
(8)
in which

$$F_{\mu \nu }^\alpha (x)=_\mu A_\nu ^\alpha (x)_\nu A_\mu ^\alpha (x)+gf^{\alpha \beta \gamma }A_\mu ^\beta (x)A_\nu ^\gamma (x)$$
(9)
We can get the motion equations of free non-Abelian gauge fields from Eq.(8)

$$_\mu F_{\mu \nu }^\alpha +gf^{\alpha \beta \gamma }A_\mu ^\beta F_{\mu \nu }^\gamma =0$$
(10)
Put Eq.(9) into (10), we obtain

$$^2A_\nu ^\alpha _\nu _\mu A_\mu ^\alpha +gf^{\alpha \beta \gamma }[(_\mu A_\mu ^\beta )A_\nu ^\gamma +A_\mu ^\beta (_\mu A_\nu ^\gamma )]$$
(11)
$$+gf^{\alpha \beta \gamma }A_\mu ^\beta (_\mu A_\nu ^\gamma _\nu A_\mu ^\gamma +gf^{\gamma \rho \sigma }A_\mu ^\rho A_\nu ^\sigma )=0$$

This is a complex non-linear equation, no solution is found up to now. Because it is not a wave equation, the non-Alebrian gauge fields seam to have no simple wave solution like Eq.(6).

On the other hand, the middle Bosons described by non-Alebrian gauge fields are also micro-particles with wave nature. From this angle, the free non-Alebrian gauge fields should also have wave solutions. In fact, in the current perturbation calculations of practical problems, we have actually supposed that the wave functions of free non-Alebrian gauge particles have the simple forms like Eq.(6). For example, in the united weak-electric interaction theory, the interaction Hamilitonian between leptons and $`W^\pm `$ particles can be written as

$$\frac{ig}{2\sqrt{2}}[W_\mu ^{}\overline{\psi }_\nu \gamma _\mu (1+\gamma _5)\psi _l+W_\mu ^+\overline{\psi }_l\gamma _\mu (1+\gamma _5)\psi _\nu ]$$
(12)
By using the Feynman rule, the vertex function of interaction between leptons and $`W^\pm `$ particles in momentum space can be written as

$$\frac{ig}{2\sqrt{2}}\gamma _\mu (1+\gamma _5)\delta ^4(p_l+p_\nu p_w)$$
(13)
In the formula, we actually acquiesce in that the wave functions of free $`W^\pm `$ particles have the simple wave solutions Eq.(6) similar to leptons. After the integral of transition matrix element is carried out in coordinate space, the formula (13) is obtained from Eq.(12). If the wave functions of free $`W^\pm `$ particles can not be written as the form of Eq.(6), we can not get vertex $`\delta `$ function, so that energy-momentum conservation in vertex processes can not be ensured and no any perturbation calculation can done.

So the difficult position we face is that in order to let perturbation calculation becoming possible, we should suppose that the wave functions of free non-Alebrian gauge particles have simple form of the wave. On the other hand, according to the motion equation of free non-Alebrian gauge fields, the wave functions of non-Alebrian gauge particles have no simple form of the wave.

This problem can be solved well by introducing a restriction condition in the motion equations of free non-Alebrian gauge fields, similar to introduce the Lorentz condition in electromagnetic theory. For this purpose, in Eq.(11) we let

$$_\nu _\mu A_\mu ^\alpha gf^{\alpha \beta \gamma }[(_\mu A_\mu ^\beta )A_\nu ^\gamma +A_\mu ^\beta (_\mu A_\nu ^\gamma )]$$
(14)
$$gf^{\alpha \beta \gamma }A_\mu ^\beta (_\mu A_\nu ^\gamma _\nu A_\mu ^\gamma +gf^{\gamma \rho \sigma }A_\mu ^\rho A_\nu ^\sigma )=0$$

The motion equation of free non-Alebrian gauge fields is simplified as

$$^2A_\mu ^\alpha =0$$
(15)
In this way, the wave functions of free non-Alebrian gauge particles can be written as simple wave’s forms

$$A_\mu ^\alpha (x)=_{\mathrm{}}^+\mathrm{}a_\mu ^\alpha (\stackrel{}{k})e^{ikx}d^3\stackrel{}{k}$$
(16)
Eq.(14) is just the restriction conditions of the motion equations of non-Alberian gauge fields. Putting Eq.(16) into (14), we get

$$k_\mu k_\nu a_\mu ^\alpha (\stackrel{}{k})e^{ikx}d^3\stackrel{}{k}+igf^{\alpha \beta \gamma }(k_\mu +k_\mu ^{})a_\mu ^\beta (\stackrel{}{k})a_\nu ^\gamma (\stackrel{}{k}^{})e^{i(k+k^{})x}d^3kd^3\stackrel{}{k}^{}$$
(17)
$$+igf^{\alpha \beta \gamma }a_\mu ^\beta (\stackrel{}{k})[k_\mu ^{}a_\nu ^\gamma (\stackrel{}{k}^{})k_\nu ^{}a_\mu ^\gamma (\stackrel{}{k}^{})]e^{i(k+k^{})x}d^3kd^3\stackrel{}{k}^{}$$

$$+g^2f^{\alpha \beta \gamma }f^{\gamma \rho \sigma }a_\mu ^\beta (veck)a_\mu ^\rho (\stackrel{}{k}^{})a_\nu ^\sigma (\stackrel{}{k}^{\prime \prime })e^{i(k+k^{}+k^{\prime \prime })x}d^3\stackrel{}{k}d^3\stackrel{}{k}^{}d^3\stackrel{}{k}^{\prime \prime }=0$$

Because $`a_\mu ^\alpha (k)`$ has nothing to do with space-time coordinate $`x`$, we integrate the formula over $`d^4x`$ and we get

$$k_\mu k_\nu a_\mu ^\alpha (\stackrel{}{k})\delta ^4(k)d^3\stackrel{}{k}+igf^{\alpha \beta \gamma }\delta ^4(k+k^{})(k_\mu +k_\mu ^{})a_\mu ^\beta (\stackrel{}{k})a_\nu ^\gamma (\stackrel{}{k}^{})d^3\stackrel{}{k}d^3\stackrel{}{k}^{}$$
(18)
$$+igf^{\alpha \beta \gamma }\delta ^4(k+k^{})a_\mu ^\beta (\stackrel{}{k})[k_\mu ^{}a_\nu ^\gamma (\stackrel{}{k}^{})k_\nu ^{}a_\mu ^\gamma (\stackrel{}{k}^{})]d^3\stackrel{}{k}d^3\stackrel{}{k}^{}$$

$$+g^2f^{\alpha \beta \gamma }f^{\gamma \rho \sigma }\delta ^4(k+k^{}+k^{\prime \prime })a_\mu ^\beta (\stackrel{}{k})a_\mu ^\rho (\stackrel{}{k}^{})a_\nu ^\sigma (\stackrel{}{k}^{\prime \prime })d^3\stackrel{}{k}d^3\stackrel{}{k}^{}d^3\stackrel{}{k}^{\prime \prime }=0$$

Take the integral of the first item about $`d^3\stackrel{}{k}^{}`$, the second and third items about $`d^3\stackrel{}{k}^{}`$ and the fourth item about $`d^3\stackrel{}{k}^{\prime \prime }`$. When $`\stackrel{}{k}\stackrel{}{k}`$, $`k_0`$ is unchanged. Let $`\overline{k}=(\stackrel{}{k},ik_0)`$, we get

$$\delta _{\nu 0}\delta (k_0)k_0^2a_0^\alpha (\stackrel{}{k})gf^{\alpha \beta \gamma }\delta (2k_0)2k_0a_0^\beta (\stackrel{}{k})a_\nu ^\gamma (\stackrel{}{k})d^3\stackrel{}{k}$$
(19)
$$+igf^{\alpha \beta \gamma }\delta (2k_0)a_\mu ^\beta (\stackrel{}{k})[\overline{k}_\mu a_\nu ^\gamma (\stackrel{}{k})\overline{k}_\nu a_\mu ^\gamma (\stackrel{}{k})]d^3\stackrel{}{k}$$

$$+g^2f^{\alpha \beta \gamma }f^{\gamma \rho \sigma }\delta (k_0+k_0^{}+k_0^{\prime \prime })a_\mu ^\beta (\stackrel{}{k})a_\mu ^\rho (\stackrel{}{k}^{})a_\nu ^\sigma (\stackrel{}{k}\stackrel{}{k}^{})d^3\stackrel{}{k}d^3\stackrel{}{k}^{}=0$$

For the gauge particles with zero masses, $`k_0=\stackrel{}{k}`$, $`k_0^{}=\stackrel{}{k}^{}`$, $`k_0^{\prime \prime }=\sqrt{\stackrel{}{k}^2+\stackrel{}{k^{}}^2}`$. Then taking differential about $`\stackrel{}{k}`$ in the formula above, we get

$$\delta _{\nu 0}\frac{}{\stackrel{}{k}}[\delta (k_0)k_0^2a_0^\alpha (\stackrel{}{k})]+gf^{\alpha \beta \gamma }\delta (2k_0)\{2k_0a_0^\beta (\stackrel{}{k})a_\nu ^\gamma (\stackrel{}{k})ia_\mu ^\beta (\stackrel{}{k})[\overline{k}_\mu a_\nu ^\gamma (\stackrel{}{k})\overline{k}_\nu a_\mu ^\gamma (\stackrel{}{k})]\}$$
(20)
$$g^2f^{\alpha \beta \gamma }f^{\gamma \rho \sigma }\delta (k_0+k_0^{}+k_0^{\prime \prime })a_\mu ^\beta (\stackrel{}{k})a_\mu ^\rho (\stackrel{}{k}^{})a_\nu ^\sigma (\stackrel{}{k}\stackrel{}{k}^{})d^3\stackrel{}{k}^{}=0$$

The formula above is non-linear differential integral equation about parameters $`a_\mu ^\alpha (\stackrel{}{k})`$ and is difficult to get its solutions. But in this way, the complexity of wave functions of non-Ablerian gauge particles is transformed into the complexity of parameters $`a_\mu ^\alpha (\stackrel{}{k})`$ and the solutions of non-Abelrian gauge fields can be written as simple forms. There are $`\alpha `$ restriction equations, so only $`3\alpha `$ are independent among $`4\alpha `$ $`a_\mu ^\alpha (\stackrel{}{k})`$. For the free non-Ablerian gauge fields with masses, the motion equation become

$$^2A_\mu ^\alpha m_\alpha ^2A_\mu ^\alpha =0$$
(21)
The form of restriction equations are unchanged, we only need to let $`\stackrel{}{k}\sqrt{\stackrel{}{k}^2+m_\alpha ^2}`$ in Eq.(20). When there exist interactions, the relative interaction item can be added into the formula, but the form of restriction equation is still unchanged.

Though the motion equations are simplified after the restriction conditions are introduced, the form of the Lagrangians of free non-Alberian gauge fields and interactions between gauge fields and other fields are unchanged. So after the solutions of motion equations of free non-Alberian gauge fields are written in the form of Eq.(16), the forms of interactions and practical calculation results are not effected. In fact, only in this way, we can use the Feynman rules in the calculation of perturbation theory, otherwise the calculations are completely impossible.