text:
# Neural Network with Semiclassical Excitation

## I Introduction

The classical integrate-and-fire neural network has been studied both in the simpler zero-width HH1  action potential and the more involved finite-width FS1 cases. In these works every neuron integrates the current coming from neighboring neurons and when the threshold for firing is exceeded, it too sends off an action potential to its neighbors. Hopfield and Herz had found that there is a simple relation between the contributions $`A`$ from the neighbors and an external current $`I`$, with the time period of the firing of the network when phase-lock is established:

$$\tau =(1A)/I$$
(1)
It has been shown FS1  that there is some modification of this formula when the action potential is not exactly a delta function but is spread over time, which is, of course, a more realistic assumption, in biological as well as physical contexts.

In view of the recent great interest in quantum computing we think it is worthwhile to investigate the changes if any that would result from converting such networks to a quantum model. As a first step, we have here tried to see the changes when the action potential acts like a quantum radiation to the neighbors, instead of a current, and the firing of the neuron is replaced by induced radiation from it to the neighbors. The quantity integrated therefore gives not a classical state with deterministic firing on reaching the threshold, but the superposed quantum amplitude, the square of whose magnitude is related to the probability of firing.

Hence, instead of the deterministic approach of the classical case, we have a quantum stochastic interpretation of the dynamics of the network. It is obvious that if the neurons are allowed to emit according to time dependent perturbation theory at random, there will be no exact phase-locking, but it may still be possible to find average periods of firing of a typical neuron which may depend on the strengths of the various parameters of the model, as in the classical case. We study this important possibility in this paper, because it may throw light on some of the most fundamental similarities and differences between classical and quantum networks.

This model is not a fully quantized network where the states of all the neurons are entangled and the operations are all unitary, as in quantum computing models. Our interest is not in the exploitation of the enlargement of the capacity of the net in the quantum form, or of the speed of information processing in such a system, but to make an investigation of the significant changes that the nonlinearity and the stochasticity may be expected to bring about. The problems of measurement make the attractive possibilities of a fully quantized networked system of gates and other devices (see e.g. CN1 ) unreachable in the near future. Hence, as a first step it is possible that quantum computers will make a beginning as semiclassical devices with all links not fully entangled by a unitary operation. How dissipation destroys the memory of the initial state and how quickly, and what periodicity, if any, emerges, depending on which parameters, are important questions that need to be studied.

## II The Model

We have simulated a square lattice of such neurons which can communicate with neighbors with quantum signals and have varied their durations, in analogy with the action potential currents of varying widths in the classical case. We have also retained a constant potential background in analogy with the constant external current.

We know

$$|t>=U(t,t_0)|t_0>$$
(2)
with

$$U(t,t_0)=exp[i𝑑tV]1+i𝑑tV$$
(3)
So that the transition rate

$$\mathrm{\Gamma }|𝑑tV|^2$$
(4)
In our case we have retained a constant $`V_o`$ background and have added the contributions $`𝑑tf(t_i)v`$ from the $`ith`$ neighbor. The integration is spread over the duration of the quantum pulse since the last triggering.

We might expect in this case

$$\tau =1/|𝑑t(V_o+\mathrm{\Sigma }f_iv)|^2$$
(5)
For small perturbations this would look like

$$\tau =k^{}(1|\mathrm{\Sigma }𝑑tf_iv|^2)/|𝑑tV_o|^2$$
(6)
which is analogous to the classical formula FS1 ; HH1 .

However, we must realize that in the classical case it is possible to have $`1=A`$ (in the standardized units used there) and get a zero period, i.e. an always saturated net. This is impossible in this quantum case because the more exact expression (1) cannot give a zero .

If we make the further assumption that the pulse width is a duration greater than the average period, then we have to take a proportionate amount of contribution from the neighbors. In this case we get the relation, for a square lattice network:

$$1/\tau =k^{}(V_0\tau +4v\tau /w)^2$$
(7)
which is a cubic equation in the period and can be solved in terms of the system parameters $`V_0`$ , width $`w`$ of the pulse and the pulse size $`v`$. The nonlinearity of the quantum version of the network mentioned earlier is responsible for this form.

Unlike the classical case we have a randomness in the contribution from the neighbors. We may try to account for this by introducing a parameter $`q`$ with the interaction $`v`$ between neighbors. This gives:

$$1/\tau =k(V_0+4qv/w)^{(2/3)}$$
(8)
## III Input Dependence

A system that does not respond to the environment is useless. So while the system period found above has its intrinsic interest, we need also study the behavior of the system to different inputs, i.e. initial states of pre-assigned nodes to serve as the interface with the environment. Obviously the system behavior will depend on the geometry and connectivity of the net and the pattern of the input. However, in this semi-quantum model we foresee a loss of memory as the nodes fire and collapse from their original states. This may be a desirable characteristic if we are interested in a system with only a short-term memory, because the automatic erasure of memory after a time may save us the software steps or the hardware needed to erase.

We shall principally study simulation experiments, as the details of memory loss are expected to depend in a complicated way on the system used and any simple closed analytic expression is unlikely. However, we shall also try to see of a simple formula such as Eq. 8 can give agreements with the simulation for the average periods of oscillation after the initial memory is effaced.

## IV Results of Simulation

Our simulations show that:

(1) for a 40X40 lattice (as taken previously  in the classical case) with periodic boundary conditions - to make it look like an infinite lattice - for about 40,000 simulations for a given set of parameters but different initial states, we get a constant average for the neurons in the lattice. So in an average sense we do have a period (Fig. 1). The system is not chaotic. For a single neuron, if we follow its history we see a fairly linear relation between the cumulative number of firings and time (Fig. 2). If we consider the pooled sum of all nodes of the lattice, the stochasticity almost disappears and we see a practically straight line (Fig. 3).

(2) Indeed, as the strength of the signal from the neighbor increases, the time period decreases (Table 1), but it does not appear to go to zero. Unlike the classical case, here we can make the signals arbitrarily strong without worrying about running into a singularity or negative periods. The parameter $`q`$ was fitted to the first value and then Eq. 8 seems capable of predicting all the other periods with remarkable accuracy.

(3) The width of the pulse is important here too. If the pulse is spread out, the average period increases (Table 2). Again we get excellent agreement of Eq. 8 with the simulation results for the periods, with the same value of $`q`$ used in Table 1.

(4) The input dependence is observed by averaging over 100 simulation runs. The following four types of inputs were used:

(a) All peripheral nodes in state $`|1>`$; all body nodes in state $`|0>`$. In this case we see (Fig. 4) a smooth transition from a state with an initial firing rate proportional to the number of initially excited nodes dying down quickly as the system forgets the input and lets the system parameters take over with a noisy pattern, despite the averaging over the runs. It is remarkable that the initial few cycles with the memory show virtually no noise.

(b) Peripheral nodes alternately in states $`|1>`$ and $`|0>`$; body nodes in state $`|0>`$ (Fig. 5). In this case we begin with a smaller number of firings on account of halving the initially excited nodes, and there are a few kinks in the initial cycles, possibly due to the conversion of the spatial lack of symmetry to temporal. As expected the system moves to the common noisy asymptotic behavior after forgetting the input.

(c) Peripheral nodes in random states between $`|0>`$ and $`|1>`$. Here too we see fairly prominent kinks (slightly less than the oscillating $`|1>|0>`$ pattern in the previous case) from the interaction of the noncoherent randomness of input in the neighboring peripheral nodes until the short-term memory disappears (Fig. 6).

(d) In this case we take the whole lattice to be initially in state $`|0>`$. However because of the external driving potential the system soon develops into a noisy final state from a smooth initial state with no excitation.

## V Discussion

We have seen that this semiclassical neural network which is designed to mimic in some ways the neural network with an integrate-and-fire model with classical action potentials indicates the presence of an average time period reminiscent of the classical model. The following differences stand out: 1) We cannot have a system with a zero time period (infinitely fast). 2) We can have arbitrarily high potentials linking the neurons.

We have seen that the average system period depends in a simple way on the system parameters as given in Eq. 8. A single parameter choice predicts a wide range of potential pulse strengths and durations.

An interesting possible extension of the ideas presented here would be to use multi-state neurons (multibits) with no classical analogue and of course also to use coherent complex quantum action to investigate the possibility of phase-locking from coherence.

We have not considered Hebbian learning by introducing dissipation. Altaisky AL1  and Zak ZK1  have investigated the effects of dissipation in quantum neural networks. In our model also dissipation is inevitable as the neurons are all allowed to fire independently, with no coherence or entanglement. However, unlike these authors, we have considered large systems where the role of adaptivity in the quantum context becomes too complex for a preliminary investigation. The system, nevertheless, can serve as short-term memory with a built-in mechanism for effacing all input dependence in the long run, which may be a desirable characteristic where learning and unlearning must go together. As the system we have considered is very basic, it is possible that more complex systems may be designed for greater functionality.

###### Acknowledgements.

I would like to thank Andrew Tan at UCSF for useful discussions.


format: text

text:
# Neutralino spectrum in top-down models of UHECR

## I Introduction

Over one hundred cosmic ray events with energy higher than $`4\times 10^{10}`$ GeV have been detected by past and present observatories. About 30 events have energy in the range $`(13)\times 10^{11}`$ GeV (see  for an experimental review). In top-down production mechanisms of Ultra High Energy Cosmic Rays (UHECR) the decay products of very massive particles $`X`$ with $`M_X>10^{12}`$ GeV account for these events. Particle X can be either a GUT scale particle produced by the decay of a cosmological network of topological defects created in a phase transition in the early universe , or a cold dark matter particle clustered in the galactic halo whose abundance was generated during or shortly after inflation . Whatever the nature of $`X`$ is, it decays into partons which hadronise and give the UHECR primaries that we observe on the Earth. Therefore, in order to calculate the predicted spectrum and composition of UHECRs one has to calculate the nonperturbative Fragmentation Functions (FFs) of partons into hadrons at the energy scale $`M_X`$. Several approaches have been taken to solve the fragmentation process and calculate the spectrum of UHECRs: Montecarlo generators , analytical estimates (MLLA)  and DGLAP evolution  (see  for a comparison of these different approaches).

As long as supersymmetry (SUSY) is a low energy symmetry of nature, the decay of $`X`$ will produce supersymmetric particles. If R-parity is conserved, the Lightest Supersymmetric Particle (LSP) is stable and should be among the primary cosmic rays that reach the Earth. Theoretical motivations favour the neutralino as the LSP, hence we shall concentrate on neutralinos***In some SUSY models the gravitino, the gluino and even the sneutrino could be the LSP. An exotic hadron state containing a stable gluino has been proposed to explain the UHECR events .. In the Minimal Supersymmetric Standard Model (MSSM) with R-parity conserved, neutralinos are a mixture of higgsinos and neutral gauginos (bino, $`\stackrel{~}{B}`$, and neutral wino, $`\stackrel{~}{W}^0`$). The precise composition of the neutralinos in terms of the interaction eigenstates depends on the particular SUSY scenario. For simplicity, we shall concentrate on the Constrained MSSM, which assumes that the soft SUSY-breaking scalar masses, the gaugino masses and the trilinear parameters are all universal at the GUT scale. In this scenario, it can be shown that in the region of the parameter space where the relic density of neutralinos is of cosmological significance, the LSP is mainly a bino .

Weakly interacting particles like neutralinos or neutrinos have cross sections with ordinary matter too small to be detected with present day UHECR observatories. Baryons, nuclei and perhaps photons may account for all the events detected so far . Forthcoming detectors will be sensitive to weakly interacting UHECRs (see Refs.  and references therein). It is therefore important to obtain a good estimate of the neutralino flux expected in top-down models.

The neutralino spectrum in top-down models has already been estimated using Montecarlo generators by Berezinsky&Kachelrieß . In the Montecarlo approach a primary parton with energy $`M_X/2`$ produced in the decay of $`X`$ initiates a parton cascade which proceeds until a specific minimal virtuality is reached. The squarks and gluinos in the parton shower decay into the LSP and ordinary partons once they reach the scale $`M_{\mathrm{SUSY}}`$, the universal mass of squarks and gluinos. All quark and gluons evolve down to the hadronization scale where a phenomenological model is employed to bind them into hadrons.

The Montecarlo simulations did not include the electroweak radiation of neutralinos by off-shell partons in the shower. In the next section we show how the DGLAP formalism employed to calculate the neutrino, photon and baryon spectra of UHECR  can be extended to calculate the spectrum of neutralinos. We concentrate on the electroweak radiation of neutralinos by off-shell partons to complement the Montecarlo simulations. We find that at large $`x`$ electroweak radiation of neutralinos is important. The large $`x`$ region is the most relevant for future cosmic ray observatories since in this region the neutralino flux is larger than the neutrino flux and will not be swamped by the neutrino signal.

## II Evolution equations for neutralino fragmentation functions

The main channels of neutralino production are the hadronic decays of $`X`$. Hence the flux of neutralinos is given by the fragmentation functions from partons, $`D_a^\chi (x,M_X^2)`$, where $`a`$ is any quark $`q_k`$ and squark flavour $`s_k`$, a gluon $`g`$ or a gluino $`\lambda `$. The flux is proportional to the sum over $`a`$ of $`D_a^\chi (x,M_X^2)`$ each one weighted with the branching ration for the decay of $`X`$ into parton $`a`$. The variable $`0<x<1`$ is the fraction of the momentum of $`a`$ carried off by $`\chi `$. These nonperturbative functions depend on the energy scale $`\mu `$; in our case $`\mu =M_X`$. The rate of change of $`D_a^\chi (x,\mu ^2)`$ with $`\mu ^2`$ is given by the sum of two terms. The first one takes into account the ordinary SUSY QCD branching of partons and gives rise to the so-called Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) equations . The second term takes into account the emission of neutralinos in the parton shower, which stems from the tree level weak SUSY coupling among quarks, squarks and neutralinos. For any (anti)quark and (anti)squark we obtain, to lowest order in the strong coupling $`\alpha _S`$ and the weak couplings $`\alpha _W^{\chi ,k}`$,

$`\mu ^2_{\mu ^2}D_{q_k}^\chi (x,\mu ^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _S(\mu ^2)}{2\pi }}\left[{\displaystyle \underset{l}{}}P_{q_lq_k}(x)D_{q_l}^\chi (x,\mu ^2)+\mathrm{}\right]+{\displaystyle \frac{\alpha _W^{\chi ,k}(\mu ^2)}{2\pi }}P_{\chi q_k}(x)D_\chi ^\chi (x,\mu ^2),`$ (1)
$`\mu ^2_{\mu ^2}D_{s_k}^\chi (x,\mu ^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _S(\mu ^2)}{2\pi }}\left[{\displaystyle \underset{l}{}}P_{q_ls_k}(x)D_{q_l}^\chi (x,\mu ^2)+\mathrm{}\right]+{\displaystyle \frac{\alpha _W^{\chi ,k}(\mu ^2)}{2\pi }}P_{\chi s_k}(x)D_\chi ^\chi (x,\mu ^2),`$ (2)
where the dots inside the square brackets stand for the remaining leading order SUSY QCD termsThe DGLAP equations do not include soft gluon emission coherence. Coherence is important at low $`x`$. Our results will only hold for $`x`$ much larger than the position of the Gaussian peak produced by coherence, at $`x_p\sqrt{\mathrm{\Lambda }_{\mathrm{QCD}}/M_X}`$. As pointed out in Sec. I we are only interested in the large $`x`$ region. The DGLAP equations can be modified to include coherence, see for example ., see . For gluons and gluinos, to lowest order in $`\alpha _W`$, there are no additional terms to the usual SUSY QCD ones. The FF of neutralinos from neutralinos is $`D_\chi ^\chi (x,\mu ^2)=\delta (1x)+O(\alpha _W)`$. Note that the functions $`D_a^\chi (x,\mu ^2)`$ are of order $`\alpha _W^{\chi ,a}`$. The splitting functions of quarks and squarks into neutralinos are given by $`P_{\chi q_k}=1x`$ and $`P_{\chi s_k}=1`$. We make use of the convolution operator

$$A(x)B(x)_x^1\frac{\mathrm{d}z}{z}A(z)B(\frac{x}{z}).$$
(3)
Analogous equations to Eqs. (12) employed to study photon structure or photon production by parton showers can be found in .

It is useful to introduce the following evolution variable

$$\tau \frac{1}{2\pi b}\mathrm{ln}\frac{\alpha _S(\mu _0^2)}{\alpha _S(\mu ^2)},$$
(4)
$`b`$ being the coefficient in the leading order $`\beta `$-function governing the running of the strong coupling: $`\beta (\alpha _S)=b\alpha _S^2`$.

For simplicity we shall only consider gaugino production. We shall argue later that the final results are not substantially altered if higgsinos are included in the calculation. Also for clarity we shall concentrate on the evolution of quark and squark singlet functions, defined as the following sum over flavours

$`D_q^\chi `$ $``$ $`{\displaystyle \underset{k}{}}D_{q_k}^\chi +D_{\overline{q}_k}^\chi ,`$ (5)
$`D_s^\chi `$ $``$ $`{\displaystyle \underset{k}{}}D_{s_k}^\chi +D_{\overline{s}_k}^\chi .`$ (6)
The singlet functions are coupled to the gluons and gluinos by means of the following $`4\times 4`$ matrix integro-differential equation

$$_\tau \left(\begin{array}{c}D_q^\chi \hfill \\ D_g^\chi \hfill \\ D_s^\chi \hfill \\ D_\lambda ^\chi \hfill \end{array}\right)=\left(\begin{array}{cccc}P_{qq}& 2n_\mathrm{F}P_{gq}& P_{sq}& 2n_\mathrm{F}P_{\lambda q}\\ P_{qg}& P_{gg}& P_{sg}& P_{\lambda g}\\ P_{qs}& 2n_\mathrm{F}P_{gs}& P_{ss}& 2n_\mathrm{F}P_{\lambda s}\\ P_{q\lambda }& P_{g\lambda }& P_{s\lambda }& P_{\lambda \lambda }\end{array}\right)\left(\begin{array}{c}D_q^\chi \hfill \\ D_g^\chi \hfill \\ D_s^\chi \hfill \\ D_\lambda ^\chi \hfill \end{array}\right)+\frac{\alpha _W^\chi }{\alpha _S}\left(\begin{array}{c}P_{\chi q}\\ 0\\ P_{\chi s}\\ 0\end{array}\right),$$
(7)
where $`\alpha _W^\chi _k\alpha _W^{\chi ,k}`$ and $`n_\mathrm{F}`$ is the number of flavours. These will be the only relevant equations if one assumes flavour universality in the decay of $`X`$. If one wishes to study the evolution of each flavour separately then one needs the complete set of DGLAP equations as given in  to which one has to add the neutralino radiation terms on the right hand side of Eqs. (12).

We can write Eq. (7) in the simple matrix form

$$_\tau xD^\chi =xPxD^\chi +xX^\chi ,$$
(8)
where $`xD^\chi `$ and $`xX^\chi `$ are $`4`$vectors and $`xP`$ is a $`4\times 4`$ matrix whose elements are all the leading order SUSY QCD splitting functions $`P_{ab}`$ times $`x`$. These functions were calculated in . We have multiplied Eq. (7) by $`x`$ to improve numerical stability.

Equation (8) is a linear inhomogeneous equation. Its associated homogeneous equation is the usual DGLAP equation for the coupled evolution of the quark and squark singlets, gluons and gluinos. The inhomogeneous or source term arises from the neutralino radiation by partons. The general solution to Eq. (8) is the sum of the general solution to the associated homogeneous equation plus a particular solution

$$xD^\chi (x,\tau )=E(x,\tau \tau _0)xD^\chi (x,\tau _0)+_{\tau _0}^\tau d\tau ^{}E(x,\tau \tau ^{})xX^\chi (x,\tau ^{}).$$
(9)
The evolution operator $`E(x,\tau )`$ is the solution to the associated homogeneous equation with the initial condition $`E(x,0)=\delta (1x)I_4`$ (see )

$$E(x,\tau )=e^{xP(x)\tau }\delta (1x)I_4+xP(x)\tau +\frac{1}{2!}xP(x)xP(x)\tau ^2+\mathrm{}$$
(10)
We take as initial energy scale $`\mu _0=M_{\mathrm{SUSY}}`$, $`\tau _0\tau (\mu _0)=0`$, the typical scale for all sparticle masses.

Equation (9) gives the total neutralino FF as the sum of two terms. The first one, the solution to the homogeneous equation, is the contribution from partons that decay on-shell into neutralinos (assuming a universal mass for all squarks and gluinos $`M_{\mathrm{SUSY}}`$). This contribution was studied with Montecarlo generators in Refs. . The second term is the contribution stemming from the electroweak radiation of neutralinos by partons. From now on, we concentrate on this second term and take

$$xD^\chi (x,\tau )=_0^\tau d\tau ^{}E(x,\tau \tau ^{})xX^\chi (x,\tau ^{}).$$
(11)
This expression is our main equation. It will allow us to estimate the relevance of electroweak radiative emission of neutralinos.

## III Numerical calculation and results

In order to calculate numerically $`xD^\chi (x,\tau )`$ as given in Eq. (11) we expand the evolution operator in Laguerre polynomialsLaguerre polynomial expansions to study scaling violations in QCD were first introduced in . $`L_n(x)`$

$`E(x,\tau )`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}E_n(\tau )L_n(\mathrm{ln}x),`$ (12)
$`E_n(\tau )`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{4}{}}}e^{\lambda _i\tau }{\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \frac{\tau ^k}{k!}}B_{i,n}^k.`$ (13)
The scalars $`\lambda _i`$ are the eigenvalues of $`xP_0`$, the first coefficient in the Laguerre expansion of the matrix $`xP(x)`$. Energy conservation gives $`\lambda _1=0`$. All other eigenvalues are negative. The $`4\times 4`$ matrices $`B_{i,n}^k`$ are given by recursive relations that were calculated in . Likewise we expand the source term in Laguerre polynomials

$`xX^\chi (x,\tau )`$ $`=`$ $`{\displaystyle \frac{\alpha _W^\chi (\tau )}{\alpha _S(\tau )}}f(x),`$ (14)
$`f(x)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}f_nL_n(\mathrm{ln}x),`$ (15)
$`f_n^T`$ $`=`$ $`(\left({\displaystyle \frac{1}{2}}\right)^{n+1}{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{2}{3}}\right)^n,0,\left({\displaystyle \frac{1}{2}}\right)^{n+1},0).`$ (16)
The Laguerre expansions for the neutralino FFs are

$$xD^\chi (x,\tau )=\underset{N=0}{\overset{\mathrm{}}{}}\stackrel{~}{C}_N^\chi (\tau )L_N(\mathrm{ln}x),$$
(17)
where $`\stackrel{~}{C}_0^\chi C_0^\chi `$ and $`\stackrel{~}{C}_N^\chi C_N^\chi C_{N1}^\chi `$ for $`N>0`$. The $`4`$vectors $`C_N^\chi `$ are calculated substituting all Laguerre expansions in Eq. (11)

$$C_N^\chi (\tau )=\underset{n=0}{\overset{N}{}}\underset{k=0}{\overset{n}{}}\underset{i=1}{\overset{4}{}}\frac{B_{i,n}^kf_{Nn}}{k!}_0^\tau d\tau ^{}e^{\lambda _i(\tau \tau ^{})}(\tau \tau ^{})^k\frac{\alpha _W^\chi (\tau ^{})}{\alpha _S(\tau ^{})}.$$
(18)
In order to present our results we multiply FFs by $`x^3`$ since the measured spectrum is usually multiplied by $`E^3`$ to highlight its structure. We plot in Fig. (1) our numerical calculation for the bino $`x^3D_a^{\stackrel{~}{B}}(x,M_X^2)`$ when $`M_X=10^{12}`$ GeV and $`M_{\mathrm{SUSY}}=400`$ GeV. The major contribution comes from squarks and quarks, which can radiate neutralinos at order $`\alpha _W`$. The contribution from gluinos and gluons is much smaller since they do not have order $`\alpha _W`$ coupling to neutralinos, therefore to lowest order they can only generate neutralinos through mixing with squarks and quarks, see Eq. (7).

The fraction of available energy $`M_X=10^{12}`$ GeV carried by binos is $`xD_a^{\stackrel{~}{B}}(x,M_X)=C_0^{\stackrel{~}{B}}(\tau (M_X))=(0.047,0.002,0.109,0.003)`$ for $`q,g,s,\lambda `$, respectively (divide by $`2n_\mathrm{F}`$ to get fraction per quark and squark flavour). From Eqs. (11) and (14) one can check that for the quark and squarks singlets the order of magnitude is given by $`xD^\chi \tau \alpha _W/\alpha _S`$.

In Fig. (2) we show how $`\stackrel{~}{B}`$ distributions change for different values of $`M_X`$. Each curve is the sum of the quark and squark singlets, gluon and gluino contributions. Note that $`_\tau D^\chi >0`$ for all $`x`$ because of the inhomogeneous term in Eq. (7). If the sparticle mass scale $`M_{\mathrm{SUSY}}`$ is around the electroweak scale, the final result depends feebly on $`M_{\mathrm{SUSY}}`$ .

The decay of $`X`$ will produce neutral winos as well. The $`\stackrel{~}{B}`$ and $`\stackrel{~}{W}^0`$ curves have similar shapes, the main difference between them being that, for a common scale $`M_X`$, $`\stackrel{~}{W}^0`$ FFs are always slightly larger than $`\stackrel{~}{B}`$ FFs, since the coupling $`\alpha _W^\chi `$ is slightly stronger for $`\stackrel{~}{W}^0`$ than for $`\stackrel{~}{B}`$. The wino will eventually decay into the bino, the LSP. Its lifetime times its Lorentz factor is much smaller than 1 kpc/$`c`$, therefore neutral winos produced in the galactic halo or further will disintegrate into the LSP before reaching the Earth. The $`\stackrel{~}{B}`$ spectrum on the Earth will be the sum of the $`\stackrel{~}{B}`$ contribution produced at the decay spot of $`X`$ plus the $`\stackrel{~}{B}`$ contribution stemming from $`\stackrel{~}{W}^0`$ decay into $`\stackrel{~}{B}`$.

The decay sequence of a neutral wino into a bino depends on the details of the scenario considered. Even within the framework of the Constrained MSSM, the freedom is still large. Nevertheless, some benchmark points have been proposed  for study at the Tevatron collider, the LHC and $`e^+e^{}`$ colliders. These points are consistent with different experimental constraints, as well as cosmology, and can be regarded as generic in the whole parameter space. We shall be only concerned with the benchmark points A-D and G-M in Ref. , where the LSP is mainly the bino and the next-to-LSP is the neutral wino, with some admixture of higgsinos. The two remaining points correspond to the so-called “focus-point” region at large scalar masses, that we shall not consider. In any case, our analysis covers a large region of the allowed CMSSM parameter space.

The neutral wino decays mostly into sleptons and leptons, and it is followed by a decay of the slepton into the LSP. Therefore, only two-body decays are relevant. Given the generic case $`AB+C`$, we want to calculate which is the distribution of the decay product $`\mathrm{d}n_C/\mathrm{d}x`$ once the distribution for the decaying particle $`\mathrm{d}n_A/\mathrm{d}x`$ is known. From the phase space of the disintegration process we obtain

$$\frac{\mathrm{d}n_C(x)}{\mathrm{d}x}=\frac{m_A^2}{\sigma _2}_{y_{min}}^{y_{max}}\frac{\mathrm{d}y}{y}\frac{\mathrm{d}n_A(y)}{\mathrm{d}y},$$
(19)
where we have defined

$`y_{max}`$ $``$ $`\text{min}[{\displaystyle \frac{\sigma _1+\sigma _2}{2m_C^2}}x,1],`$ (20)
$`y_{min}`$ $``$ $`{\displaystyle \frac{\sigma _1\sigma _2}{2m_C^2}}x,`$ (21)
$`\sigma _1`$ $``$ $`m_A^2+m_C^2m_B^2,`$ (22)
$`\sigma _2`$ $``$ $`\sqrt{\lambda (m_A^2,m_B^2,m_C^2)},`$ (23)
$`\lambda (a,b,c)`$ $``$ $`a^2+b^2+c^22ab2ac2bc.`$ (24)
In addition there is the following kinematical constraint

$$0<x<\frac{2m_C^2}{\sigma _1\sigma _2}.$$
(25)
The LSP flux reaching the atmosphere is the sum of the $`\stackrel{~}{B}`$ contribution produced directly in the decay of $`X`$ plus the $`\stackrel{~}{B}`$ contribution produced by the decay of $`\stackrel{~}{W}^0`$ on its way to the Earth. The two-body decay of winos into binos pushes the momenta to lower values of $`x`$, making the direct $`\stackrel{~}{B}`$ component dominate over the component produced by $`\stackrel{~}{W}^0`$ decay, for $`x>0.1`$.

So far we have ignored higgsino production, so some words are now in order. The magnitude of the gaugino and higgsino couplings to partons are comparable, hence their FFs will be also similar. For the benchmark points that we have analysed, higgsinos will decay on their way to the Earth into winos and subsequently into binos. Therefore, for the same kinematical reason as the wino case discussed above, the total number of binos on Earth coming from higgsino decay will be a small correction to the direct bino contribution, for $`x>0.1`$. For the same reason, we have ignored binos coming from the decay of the charginos that are also produced in the parton shower.

We compare in Fig. 3 the spectra of baryons, neutrinos and LSPs expected in the case that UHECRs are produced by the slow decay of a population of superheavy dark matter particles $`X`$ clustered in the galactic halo. We take as dark matter particle mass $`M_X=10^{12}`$ GeV and as scale at which SUSY switches on $`M_{\mathrm{SUSY}}=400`$ GeV. The baryon and neutrino curves are obtained from Ref. . The neutralino spectrum shown corresponds to benchmark point C in Ref. . The final neutralino spectrum does not depend significantly on which benchmark point we choose, since the dominant contribution in the region of interest is given by the bino component produced directly by $`X`$.

## IV Conclusions

We have used the DGLAP formalism to calculate the spectrum of neutralinos produced by the decay of particles with mass $`M_X>10^{12}`$ GeV. We have concentrated on the electroweak radiation of neutralinos by partons. This contribution has to be added to the on-shell contribution calculated with Montecarlo generators. We find that at large $`x`$ the radiative contribution is slightly larger.

For $`x>0.3`$, LSPs dominate over baryons, photons and neutrinos and their flux on Earth is non-negligible. If UHECR are produced by the decay of superheavy particles and SUSY is a low energy symmetry of nature, then forthcoming observatories sensitive to weakly interacting UHECRs should detect a flux of ultra high energy neutralinos.

## Addendum

As this manuscript was being finished, Ref. appeared. This paper also discusses neutralino production by the decay of super-heavy dark matter particles using DGLAP evolution.

## Acknowledgements

We would like to thank Venya Berezinsky, Celine Boehm, Graham Ross and specially Subir Sarkar for discussion. A.I. would like to thank PPARC for financial support. R.T. was supported by a Marie Curie Fellowship No. HPMF-CT-1999-00268 and by the Leverhulme Trust.


text:
# Comment on “Time-like flows of energy-momentum and particle trajectories for the Klein–Gordon equation”

## Abstract

Horton, Dewdney, and Nesteruk  have proposed Bohm-type particle trajectories accompanying a Klein–Gordon wave function $`\psi `$ on Minkowski space. From two vector fields on space-time, $`W^+`$ and $`W^{}`$, defined in terms of $`\psi `$, they intend to construct a timelike vector field $`W`$, the integral curves of which are the possible trajectories, by the following rule: at every space-time point, take either $`W=W^+`$ or $`W=W^{}`$ depending on which is timelike.

This procedure, however, is ill-defined as soon as both are timelike, or both spacelike. Indeed, they cannot both be timelike, but they can well both be spacelike, contrary to the central claim of . We point out the gap in their proof, provide a counterexample, and argue that, even for a rather arbitrary wave function, the points where both $`W^+`$ and $`W^{}`$ are spacelike can form a set of positive measure.

Let $`\psi =\mathrm{e}^{P+\mathrm{i}S}`$ (where $`P`$ and $`S`$ are real) solve the Klein-Gordon equation, $`\mathrm{}\psi =m^2\psi `$. Set $`P_\mu =_\mu P`$, $`S_\mu =_\mu S`$, and

$$\theta =\text{sinh}^1\frac{P^\mu P_\mu S^\mu S_\mu }{2P^\mu S_\mu }.$$

That $`P_\mu `$ and $`S_\mu `$ are orthogonal is an exceptional case that we neglect, like the authors of . For $`W_\mu `$ one is supposed to take either $`W_\mu ^+=\mathrm{e}^\theta P_\mu +S_\mu `$ or $`W_\mu ^{}=\mathrm{e}^\theta P_\mu +S_\mu `$, depending on which is timelike; they cannot both be timelike since they are orthogonal. The question is, could they both be spacelike?

The authors of  declare that $`W_\mu ^+`$ and $`W_\mu ^{}`$ cannot both be spacelike and argue like this: otherwise there exists a Lorentz frame such that $`W_0^+=0`$ and $`W_0^{}=0`$, thus $`\mathrm{e}^\theta P_0=\mathrm{e}^\theta P_0`$, from which they conclude $`\mathrm{e}^\theta =\mathrm{e}^\theta `$, which is impossible.

It is correct that any two orthogonal spacelike vectors span a spacelike 2-plane (corresponding to $`x^0=0`$ in the appropriate Lorentz frame), but no contradiction arises since in this case $`P_0`$ would be $`0`$ (in this frame). This is the mistake in the proof.

Together with $`W_0^+=0`$ (or $`W_0^{}=0`$), $`P_0=0`$ implies $`S_0=0`$. Hence, for $`W_\mu ^+`$ and $`W_\mu ^{}`$ to be spacelike it is necessary and sufficient that $`P_\mu `$ and $`S_\mu `$ span a spacelike 2-plane.

Can this case occur? Clearly: since the Klein–Gordon equation is of second order, one may choose $`\psi `$ and $`_0\psi `$ ad libitum on the $`x^0=0`$ hyperplane. Can it also occur for the first-order Klein–Gordon equation $`\mathrm{i}_0\psi =\sqrt{m^2\mathrm{\Delta }}\psi `$, or, equivalently, for functions from the positive-energy subspace? Here is an example: let $`\psi `$ be a superposition of three222Two will not suffice for an example since $`P_\mu `$ and $`S_\mu `$ are linear combinations of the $`k_\mu `$ vectors. plane waves

$$\psi (x)=\underset{i=1}{\overset{3}{}}c_i\mathrm{e}^{\mathrm{i}k_\mu ^{(i)}x^\mu }$$

with wave vectors $`k_\mu ^{(1)}=(m,0,0,0)`$, $`k_\mu ^{(2)}=(\sqrt{27}m,\sqrt{26}m,0,0)`$, $`k_\mu ^{(3)}=(\sqrt{27}m,0,\sqrt{26}m,0)`$, and $`c_1=3`$, $`c_2=1/\sqrt{3}\mathrm{i}`$, $`c_3=\mathrm{i}`$. Then, at the coordinate origin, we find $`P_\mu =(0,\alpha ,\alpha ,0)`$ and $`S_\mu =(0,\beta ,0,0)`$ with $`\alpha =\sqrt{26}m/\gamma `$, $`\beta =\alpha /\sqrt{3}`$ and $`\gamma =31/\sqrt{3}`$. This example could also be made square-integrable by replacing the plane waves $`\mathrm{exp}(\mathrm{i}k_\mu ^{(i)}x^\mu )`$ by positive-energy $`L^2`$ Klein–Gordon functions $`\phi ^{(i)}(x)`$ with the properties $`\phi ^{(i)}(0)=1`$ and $`_\mu \phi ^{(i)}(0)=\mathrm{i}k_\mu ^{(i)}`$.

One may suspect, however, that perhaps this particular wave function $`\psi `$ is very exceptional, and perhaps even that for this special wave function the coordinate origin is a rather atypical point, so that the sort of situation just described can be ignored. After all, we would be willing to ignore the case where $`P_\mu `$ and $`S_\mu `$ are orthogonal because in the 8-dimensional space of all possible pairs of vectors $`P_\mu ,S_\mu `$ it corresponds to a subset of dimension 7, and therefore one would expect that the space-time points where this happens form a set of measure zero.

But since $`W_\mu ^+W^{+\mu }`$ and $`W_\mu ^{}W^\mu `$ are continuous functions of $`P_\mu `$ and $`S_\mu `$, the set of pairs $`P_\mu ,S_\mu `$ where both $`W^+`$ and $`W^{}`$ are spacelike is open (and nonempty) and thus has positive measure in 8 dimensions. I know of nothing precluding any pairs $`P_\mu ,S_\mu `$ from arising from a Klein–Gordon wave function, so it seems reasonable to expect that the space-time points with spacelike $`W^+`$ and $`W^{}`$ form a set of positive measure for many wave functions, perhaps for most.

Acknowledgements. I wish to thank Sheldon Goldstein for improvements and simplifications of the arguments.


text:
# Untitled Document

> TAUP 2686-2001

Remarks on the Physical Meaning of
Diffraction-Free Solutions of Maxwell Equations

E. Comay

School of Physics and Astronomy

Raymond and Beverly Sackler Faculty of Exact Sciences

Tel Aviv University

Tel Aviv 69978

Israel

PACS No: 03.50.De, 41.20.Jb Abstract:

It is proved that a source of electromagnetic radiation cannot emit a diffraction-free beam at the wave zone. A Bessel $`J_0`$ $`\phi `$-invariant beam does not hold even at the intermediate zone. These results negate claims published recently in the literature.

An idea of creating a diffraction-free beam has been published. The beam’s amplitude is cylindrically symmetric ($`\phi `$-invariant) where the $`r`$-dependence is proportional to the Bessel function of the first kind $`J_0(ar)`$ and $`a`$ is a constant having the dimension of $`L^1`$. Reference  has arisen a great interest in utilization of $`J_0`$ beams and has been cited more than 360 times. An application of  shows the central peak of the assumed $`J_0`$ beam and another one refers to its peculiar $`z`$-component wavelength. Another publication related to  claims that a superluminal propagation of light in air has been detected. Objections to  have been published. The purpose of this work is to show that one cannot construct a diffraction-free electromagnetic beam at the wave zone and that the Bessel function $`J_0(ar)`$ is unsuitable for describing diffraction free $`\phi `$-invariant wave at the intermediate zone too. This outcome proves that results of papers discussing this topic, in general, and those ascribing superluminal velocity to beams that take the form of Bessel function $`J_0`$, in particular, should be reevaluated. Units where the speed of light $`c=1`$ are used. The metric $`g_{\alpha \beta }`$ is diagonal and its entries are (1,-1,-1,-1). $`𝐮_r`$, $`𝐮_\phi `$ and $`𝐮_z`$ denote unit vectors in cylindrical coordinates and $`𝐮_x`$, $`𝐮_y`$ and $`𝐮_z`$ are unit vectors in Cartesian coordinates.

A general analysis of diffraction-free solutions of Maxwell equations has been published. Here the fields solving the problem are derived from a vector potential $`𝐀`$ that satisfies the wave equation together with the Lorentz-gauge requirement. It turns out that this work is relevant to  and some of its results are analyzed here in detail. (Another work is closely related to  and .) Let us start with the solution obtained in Example 1 (on p. 1557 of ). Using cylindrical coordinates and removing constant factors, the time dependent monochromatic electric field of this solution is obtained from the vector potential $`𝐄=𝐀/t`$

$$𝐄=\omega J_1(ar)e^{i(bz\omega t)}𝐮_\phi $$
(1)
where $`J_1`$ is the Bessel function of the first kind of order 1. The magnetic field is $`𝐁=curl𝐀`$

$$𝐁=bJ_1(ar)e^{i(bz\omega t)}𝐮_riaJ_0(ar)e^{i(bz\omega t)}𝐮_z.$$
(2)
Ignoring constant factors, one finds that the magnetic field $`(\text{2})`$ is dual to the electric field of Example 2 of . (The factor 2 in $`U_r`$ of example 2 is a misprint.) This outcome indicates that Examples 1 and 2 of  represent dual electromagnetic solutions where $`𝐄𝐁,𝐁𝐄`$ (see , p. 252).

Having the solution, let us examine the problem of a cylindrically shaped wave guide whose walls are made of a perfect conductor (see , p. 335). The length of the cylinder is much greater than both its diameter $`2R`$ and the wavelength $`\lambda =1/\omega `$ (see fig. 1). The boundary conditions along the wave guide’s walls are (see , p. 335)

$$E_{}=0,B_{}=0.$$
(3)
Thus, the solution $`(\text{1})`$ and $`(\text{2})`$ satisfies the boundary conditions provided

$$J_1(aR)=0.$$
(4)
Dynamical properties of the solution $`(\text{1})`$ and $`(\text{2})`$ are obtained from the energy-momentum tensor of the electromagnetic fields (see , p. 81 or , p. 605))

$$T_F^{\mu \nu }=\frac{1}{4\pi }(F^{\mu \alpha }F^{\beta \nu }g_{\alpha \beta }+\frac{1}{4}F^{\alpha \beta }F_{\alpha \beta }g^{\mu \nu })$$
(5)
where $`F^{\mu \nu }`$ denotes the tensor of the electromagnetic fields. Expression $`(\text{5})`$ is quadratic in the fields. Hence, one should use the real part of $`(\text{1})`$ and $`(\text{2})`$ in an evaluation of quantities belonging to it. Let us first examine the momentum density of the fields. This is the Poynting vector

$$𝐒=\frac{1}{4\pi }𝐄\times 𝐁.$$
(6)
The $`z`$-component of the momentum density and energy flux are obtained from the substitution of the appropriate real part of $`(\text{1})`$ and $`(\text{2})`$

$$S_z=\frac{b\omega }{4\pi }J_1^2(ar)cos^2(bz\omega t).$$
(7)
Expression $`(\text{7})`$ is non-negative at all points, a property which is consistent with the beam’s expected flux of energy that travels away from a localized source.

The radial component of the momentum density is obtained analogously

$$S_r=\frac{a\omega }{8\pi }J_1(ar)J_0(ar)sin[2(kz\omega t)].$$
(8)
Here one sees that, unlike the case of $`(\text{7})`$, the sign of $`(\text{8})`$ alternates periodically in the time and $`z`$-coordinates. Moreover, for any fixed value of $`t`$ and $`z`$, it changes sign along the $`r`$-axis, because zeroes of the Bessel functions $`J_0`$ and $`J_1`$ do not coincide. It follows that although the radial motion does not vanish locally, its mean value is null. This property indicates that the radial motion takes the type of a standing wave.

Now let us examine the interaction of the fields with the walls of the wave guide. Point $`P`$ at $`x=R,y=z=0`$ is used as a representation of the general case and cartesian coordinates are used. The $`x`$-component of the momentum current at $`P`$ is (see , p. 82 or , p. 605))

$$T_{xx}=\frac{1}{8\pi }(E_y^2+E_z^2E_x^2+B_y^2+B_z^2B_x^2).$$
(9)
Examining the fields $`(\text{1})`$ and $`(\text{2})`$ and the boundary value $`(\text{4})`$, one finds that only the $`z`$-component of the magnetic field makes a nonvanishing contribution. Thus, the momentum current at $`P`$ is

$$T_{xx}=\frac{a^2}{8\pi }J_0^2(aR)sin^2(bz\omega t).$$
(10)
This momentum current is absorbed by the walls, because the fields vanish in all space outside the inner part of the wave guide.

Another effect of the magnetic field $`(\text{2})`$ on the wave guide’s walls is the electric current induced in the $`\phi `$-direction. Indeed, let us evaluate the line integral along the infinitesimal rectangular closed path of fig. 1. Using vector analysis, Maxwell equations and the boundary condition $`(\text{4})`$, one finds

$$𝐁𝑑𝐥=curl𝐁𝑑𝐬=4\pi 𝐣𝑑𝐬.$$
(11)
Thus, a nonzero current $`𝐣`$ is induced on the walls, because only $`B_z`$ at the inner part of the wave guide makes a nonvanishing contribution to the line integral. This outcome proves that a time-dependent (and $`z`$-dependent) electric current flows along the $`\phi `$-direction of the wave guide’s walls and that fields of this current are part of the solution $`(\text{1})`$ and $`(\text{2})`$. This electric current sustains the $`B_z`$ related standing wave in the radial direction. The walls also counteract against local electromagnetic pressure.

The dual solution of example 2 of  behaves analogously. Using the same global factor of $`(\text{1})`$ and $`(\text{2})`$, one finds for this case

$$𝐁=\omega J_1(ar)e^{i(bz\omega t)}𝐮_\phi $$
(12)
$$𝐄=bJ_1(ar)e^{i(bz\omega t)}𝐮_r+iaJ_0(ar)e^{i(bz\omega t)}𝐮_z.$$
(13)
Hence, the boundary conditions $`(\text{3})`$ yield

$$J_0(aR)=0.$$
(14)
Since $`J_0(ar)`$ and $`J_1(ar)`$ have no common root, a nonvanishing radial electric field exists at the wave guide’s walls. It follows from Maxwell equation $`div𝐄=4\pi \rho `$ that a time dependent and $`z`$-dependent charge density is built on the inner part of the wave guide’s walls. Thus, we have also in Example 2 a current that flows on the walls and affects the fields inside the wave guide.

Let us examine an analogous experimental setup. Here the source of the radiation at $`z=L`$ is the same as that of the first experiment but the wave guide is removed. This situation is different from the wave guide case. Indeed, the fields of a closed electromagnetic system depend on charges and currents at the retarded space-time points (see , pp. 158-160 or , p. 225). Therefore, the wave guide’s solutions clearly do not hold for this case because here the current along the wave guide walls is missing.

Since in the second experiment the region at $`z=0`$ satisfies the wave zone requirements (see  p. 170 or , p. 392)

$$L\lambda ,L2R,$$
(15)
one can use the wave zone solution. Let $`𝐀`$ denote the retarded vector potential at the wave zone. Thus, one finds the fields (see  p. 171)

$$𝐁=\dot{𝐀}\times 𝐧,$$
(16)
$$𝐄=(\dot{𝐀}\times 𝐧)\times 𝐧$$
(17)
where $`𝐧`$ is a unit vector in the radial direction.

It turns out that the solution for the free space experiment is inherently different from the one which fits the wave guide’s inner space. In particular, in the case of free space, fields at the wave zone are perpendicular to the radius vector from the source to the field point. On the other hand, the wave guide solution contains a $`z`$-component ($`B_z`$ or $`E_z`$) which is an inherent part of the solution. As shown above, the $`B_z`$ (or $`E_z`$) field is associated with the electric current induced on the wave guide’s walls. This conclusion obviously holds for any pattern of source elements put at the same spatial region as the one used here, because the analysis does not refer to the source’s details. Thus, the results disagree with the claim of .

One can use general arguments for proving that a diffraction-free electromagnetic beam that has a nonvanishing $`z`$-component for at least one of the fields, contains transverse standing wave. Indeed, the beam carries energy and therefore $`𝐒`$ of $`(\text{6})`$ does not vanish. Hence, $`𝐄`$ is not parallel to $`𝐁`$ and, due to the $`z`$-component of the fields, $`𝐒`$ has a nonvanishing transverse component. Now, the diffraction-free property of the beam prevents energy from flowing transversally. Hence, the transverse component of $`𝐒`$ is a standing wave.

It can also be proved that all solutions of  have a nonvanishing $`z`$-component of at least one of the fields. Indeed, the vector potential $`𝐀`$ takes the form (see p. 1556 therein)

$$𝐀=\underset{n}{}(\alpha _n𝐌_n+\beta _n𝐍_n),$$
(18)
where $`\alpha _n`$ and $`\beta _n`$ are numerical coefficients of the expansion. Here

$$𝐌_n=curl[J_n(ar)e^{i(bz+n\phi \omega t)}𝐮_𝐳]$$
(19)
and

$$𝐍_n=\frac{1}{k}curl𝐌_n$$
(20)
where $`k`$ is the wave number. Now $`𝐍_n`$ contains a $`z`$-component (see p. 1557 therein). Hence, if $`\beta _n0`$ then $`𝐄=𝐀/t=i\omega 𝐀`$ has a $`z`$-component too. In other cases all $`\beta _n=0`$, which mean that for at least one $`n`$, $`\alpha _n0`$. Here the magnetic field $`𝐁=curl𝐀=\alpha _ncurl𝐌=k\alpha _n𝐍`$, which means that $`B_z0`$ and the proof is completed.

It follows that the family of solutions of  involves standing waves associated with the $`z`$-components of the solutions. This diffraction-free family of solutions may fit cylindrical wave guides but are unsuitable for the case of a free space.

Example 4 of  (see p. 1558) is the last one which is analyzed here in detail. This example contains one component which is proportional to $`J_0(ar)`$ and is $`\phi `$-invariant. Although it has a $`\phi `$-dependent $`z`$-component term which is associated with a standing wave, it looks simpler to show another problem of this solution. The vector potential of this example is given in Cartesian coordinates

$$𝐀=i\alpha [aJ_0(ar)𝐮_xi\frac{a^2}{b}J_1(ar)cos\phi 𝐮_z]e^{i(bz\omega t)}.$$
(21)
Using $`𝐄=𝐀/t`$, one finds

$$𝐄=\alpha \omega [aJ_0(ar)𝐮_xi\frac{a^2}{b}J_1(ar)cos\phi 𝐮_z]e^{i(bz\omega t)}.$$
(22)
Let us examine the $`z`$-component of the Poynting vector which represents energy current flowing along the beam’s direction, namely, the quantity which is analogous to $`(\text{7})`$ of Example 1. Examining $`(\text{22})`$, one finds that only $`B_y`$ is needed for this purpose. Thus, $`(curl𝐀)_y`$ of $`(\text{21})`$ is

$$B_y=\alpha [(ab\frac{a^3}{2b})J_0(ar)+\frac{a^3}{2b}cos2\phi J_2(ar)]e^{i(bz\omega t)}.$$
(23)
Hence, the required $`z`$-component of the Poynting vector is obtained as the product of the real parts of $`E_x`$ of $`(\text{22})`$ and $`B_y`$ of $`(\text{23})`$

$$S_z=\alpha ^2\omega [(a^2b\frac{a^4}{2b})J_0^2(ar)+\frac{a^4}{2b}cos2\phi J_0(ar)J_2(ar)]cos^2(bz\omega t).$$
(24)
Let us examine the $`z`$-component of the energy current near a point whose radial coordinate is $`\overline{R}`$ and $`J_0(a\overline{R})=0`$. In this neighbourhood $`J_2`$ is dominant and the contribution of the $`J_0^2(ar)`$ term of $`(\text{24})`$ can be ignored. The rest of $`(\text{24})`$ is proportional to $`J_0(ar)J_2(ar)cos\mathrm{\hspace{0.17em}2}\phi `$. Now, let us examine the value of $`S_z`$ on a circle whose radius is $`\overline{R}+\epsilon `$, where $`\epsilon `$ is an appropriate small quantity. Due to the factor $`cos\mathrm{\hspace{0.17em}2}\phi `$, one realizes that $`S_z`$ takes different signs on this circle. Hence, in the solution of Example 4 of , energy flows in opposite $`z`$-directions in certain regions of space. This property of Example 4 is inconsistent with the notion of a beam, where electromagnetic energy flows away from a localized source.

It is clear from the analysis carried out above that, in free space, one cannot build a diffraction free beam from the family of Bessel functions of , because these functions are unsuitable at the wave zone.

Some conclusions can be drawn for the intermediate zone too. The diffraction free $`\phi `$-invariant $`J_0(ar)`$ function proposed in  does not belong to the solutions of . Indeed, in , there are only two truly $`\phi `$-invariant solutions. They are the dual solutions of Examples 1 and 2 which are discussed above. As proved in this work, the $`z`$-component of the energy current is proportional to $`J_1^2(ar)`$. Hence, the flow of energy vanishes along the $`z`$-axis. It is also proved above that Example 4 of , where there is one $`J_0`$ term which is $`\phi `$-invariant, does not describe a beam of electromagnetic radiation and its $`z`$-component is not $`\phi `$-invariant. It follows that experiments using a $`\phi `$-invariant setup and showing a strong peak at the center (like ) should not be interpreted by means of diffraction free solutions.

References:

* Email: eli@tauphy.tau.ac.il
* J. Durnin, J. J. Miceli, Jr. and J. H. Eberly, Phys. Rev. Lett. 58, 1499 (1987).
* Due to this large number of papers, the rather short reference list presented here cannot be regarded as an adequate representation of the relevant literature. For having further reference, readers may use articles mentioned here or a citation database.
* C. A. McQueen, J. Arit and K. Dholakia, Am. J. Phys. 67, 912 (1999).
* T. Wulle and S. Herminghaus, Phys. Rev. Lett. 70, 1401 (1993).
* D. Mugnai, A. Ranfagni and R. Ruggeri, Phys. Rev. Lett. 84, 4830 (2000).
* N. P. Bigelow and C. R. Hagen, Phys. Rev. Lett. 87, 059401 (2001); H. Ringermacher and L. R. Mead, Phys. Rev. Lett. 87, 059402 (2001); Thilo Sauter and Fritz Paschke, Phys. Lett. 285, 1 (2001).
* Z. Bouchal and M. Olivik, J. Mod. Opt. 42, 1555 (1995).
* J. D. Jackson, Classical Electrodynamics (John Wiley, New York,1975). p. 220.
* S. V. Kukhlevsky, G. Nyitray and V. L. Kantsyrev, Phys. Rev. E64, 026603 (2001).
* L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields (Pergamon, Oxford, 1975). P. 81.
* M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, (U.S. Government Printing Office, Washington, 1972). p. 370.
* Due to , the roots $`(r>0)`$ of $`J_n(r)`$ and $`J_{n+1}(r)`$ are simple, do not coincide and interlace. Hence, the recurrence formula $`2J_1(r)/r=J_0(r)+J_2(r)`$ proves that positive roots of $`J_0(r)`$ and $`J_2(r)`$ do not coincide.

Figure captions:

Fig. 1:

Electromagnetic radiation is emitted from a source into a cylindrical wave guide whose radius is $`R`$. The source is at $`z=L`$ and $`L2R`$. $`O`$ denotes the origin of coordinates and the rectangle at point $`P`$ denotes a closed integration path (see text).


text:
# Temperature-dependent vibrational heterogeneities in harmonic glasses

\[

## Abstract

Numerical simulation is employed to study dynamical heterogeneities in model harmonic glasses whose atoms interact via three variants of the Lennard-Jones potential (monoatomic full Lennard-Jones, soft spheres, binary mixture). Heterogeneities are observed to exist in all three kinds of glasses, and in some cases they are observed to depend on temperature. The dimension of the heterogeneities is studied for the full Lennard-Jones case.

\]

The question as to whether the structure of glasses at the scale of tens to hundreds of interatomic distances is homogeneous or inhomogeneous has attracted recently much interest. Traditionally, following the continuous random network model suggested by Zachariasen , the homogeneous hypothesis has prevailed, also because no heterogeneity was clearly observed by small angle neutrron or X-ray scattering and electron microscopy. However, these techniques show only that there is no evident density heterogeneity, but tell nothing about the cohesion or elasticity at the nanometric scale, and thus cannot rule out vibrational dynamical heterogeneities.

The inhomogeneous cohesion of glasses at the nanometric scale was invoked to interpret the excess of vibrational density of states in the terahertz frequency range, the so-called boson peak . By this interpretation, the frequency of the boson peak is related to the inhomogeneity size.

Actually, recent molecular dynamics simulations on binary Lennard-Jones (LJ) mixtures  showed that in the supercooled state highly mobile and immobile particles are spatially correlated over a range that grows with temperature as the glass transition is approached; this finding is consistent with the observation of dynamical heterogeneities in supercooled liquids . One important issue is whether or not such heterogeneities are in some sense ”freezed down” through the glass transition, so that even in the cold, harmonic glass there exists a memory of these heterogeneities and, as a consequence, softer and harder zones.

In another recent simulation this problem was addressed by calculating the local pressure or stress in a glassy LJ system , and a strong correlation was found among atoms experiencing the same local pressure. In this reference, the elastic constants pertinent to one atom (which are related to local pressure) were determined by keeping all the other atoms fixed , and if the average local elastic constant was large (small) the atom was ranked as hard (soft). Finally, very recently the vibrational modes of 2-D LJ clusters were determined numerically ; it was found that at length scales less than approximately 40 atomic distances the elasticity is no more affine with respect the sample deformation, and that heterogeneities of force constants are extended on the same nanometric lengths .

In the present paper we investigate the possible effect of temperature on the heterogeneities which are found in harmonic glasses of the LJ family. The existence of such effects is expected in principle on the basis of the following argument: Consider a glass at low temperature; in this case only the modes of low frequency will be excited and the atoms will only perform the cooperative motions corresponding to the normal modes in question; therefore, irrespective of the magnitude of the local force constants, motion in a given direction will only occur if allowed by the excitable modes. In this sense, the heterogeneities observed in Ref.  can describe the high-temperature behaviour, and otherwise it is expected that if disomogeneities exist, their shape, dimension, quantity etc. may depend on temperature.

We considered monoatomic systems interacting through the LJ potential ($`N_0=`$2048 atoms) or the repulsive part of the LJ potential (soft spheres, $`N_0=`$2048), and the bynary LJ mixture introduced in Refs.  ($`N_0=`$2000). The potential parameters, mass and density of the LJ and soft-sphere system were those suitable for Argon ; for the binary mixture we chose the potential parameters as in Refs. . The stable configurations were produced by a fast quench (down to typically $`T_m/10`$, where $`T_m`$ is the melting temperature) of the liquid configurations starting from temperatures just above $`T_m`$. The liquid configurations were obtained by standard MD equilibration runs. The frozen configurations at $`T_m/10`$ were then further equilibrated in order to check against crystallization. Finally, a steepest descent was applied to find the (glassy) minimum configuration at $`T=0`$. For all systems we averaged the quantities of interest over 5 different realizations.

Following the procedure of Ref. , soft (hard) atoms, relative to a given eigenvector, are identified as those with small (large) variation of potential energy under the displacement pattern produced by the eigenvector itself. Therefore, in the harmonic approximation, the dynamical matrix was diagonalized, thus obtaining eigenvalues and eigenvectors. Subsequently, for each atom $`i`$ and for each normal mode $`p`$, we calculated the potential-energy change due to the collective motion determined by the normal mode in question:

$$\delta V_i^p(T)=\underset{ji}{}V(r_{ij}^p(T))V(0)$$

Here the sum runs on all atoms other than $`i`$, $`r_{ij}^p(T)`$ is the set of atomic distances resulting from the displacements produced by normal mode $`p`$, whose amplitude depends on temperature; $`V(0)`$ is the potential energy of the equilibrium configuration. At a given temperature, all normal modes are excited, with amplitude depending on their energies $`\mathrm{}\omega _p`$. Since the motions due to the different modes are not mutually coherent, we assume that the total potential energy variation is the sum of the weighted variations produced by the single modes:

$$\delta V_i(T)=\underset{p}{}\delta V_i^p(T)$$
(1)
The distances $`r_{ij}^p(T)`$ can be written as $`[r_i^p(T)r_j^p(T)]`$, where the atomic positions depend on the eigenvector displacements:

$$r_i^p(T)=r_i^0+Q^p(T)e_i^p,$$

and the amplitudes $`Q^p`$ are the temperature-dependent quantities. For a normal mode of frequency $`\omega _p`$, and at temperature $`T`$, the average vibrational level is

$$<n^p(\omega _p,T)>=\frac{1}{[exp(\mathrm{}\omega _p/kT)1]},$$

and the expectation value of $`Q^2`$ is calculated in the state $`|n>`$ as:

$$Q^2=<n|Q^2|n>=\mathrm{}/(m\omega _p)(2n+1),$$

and consequently $`Q^p(T)=\sqrt{Q^2}`$ is determined as a function of $`T`$ and $`\omega _p`$, and the temperature-dependent potential-energy variations are computed through Eq. (1). In what follows, temperatures are alway expressed in units of the maximum frequency $`\omega _M`$ of the respective systems, $`T=\frac{\mathrm{}\omega _M}{k}`$.

The presence of vibrational heterogeneities is detected by comparing the normalized pair correlation function of the whole system $`g_W(R)`$ to that of the $`N`$ softest and/or hardest atoms, $`g_N(R)`$ ; in fact, for example, if the nearest-neighbor or next-nearest-neighbor peaks in $`g_N`$ are higher than in $`g_W`$ it means that the $`N`$ softer or harder atoms tend to group.

In Fig. 1 we report the pair correlation functions $`g_N(R)`$ of the monoatomic full LJ system for the $`N=20,50,100`$ and $`200`$ softest atoms at $`T=0.3`$, and $`g_W(R)`$ (which is $`T`$-independent). Clear differences appear in the first peak: that of all selected groups of soft atoms is narrower, more intense than for the whole, and shifted to shorter distances, indicating that the soft atoms belong to the low-$`R`$ tail. Excess correlation is observed on the next-neighbour peaks, but only for the smallest value $`N=20`$, whose $`g_N(R)`$ is very noisy: however, the difference is so marked that such peaks might be significant.

In Fig. 2 is shown, for the same system, the temperature variation of the nearest-neighbour peak of $`g_N`$ for $`N=100`$, for soft ($`a`$) and hard ($`b`$) atoms. A marked temperature dependence is observed for the soft ones, together with a jump in correlation between $`T=0.1`$ and $`T=0.3`$. For both kinds of atoms, the peak is shifted to low $`R`$ values, but for the soft ones it becomes narrower than in $`g_W`$ at high $`T`$; this can be understood by considering that the modes of low frequency, i.e. of long wavelength, are weakly dependent on the local variation of the bonding between neighbouring atoms, the dependence increases with the frequency of the modes so that the softest atoms become better differentiated from the whole. The hard atoms show little or no temperature variation on the highest point, but contrary to soft ones the overall area under the nearest-neighbour peak is increased with respect to $`g_W`$, especially at the two highest $`T`$-values.

The nearest-neighbour peaks of soft spheres are reported in Fig. 3. In this case again there is very little $`T`$-variation, but the excess correlation is quite marked across the whole peak whose area is much larger than in $`g_W`$; the low-$`R`$ side of hard atoms seems again to be favoured.

In Fig. 4 are reported the nearest-neighbour peaks for the LJ binary mixture. The major effects concern the low-$`R`$ tail of the $`AA`$ peak and the whole $`AB`$ peak. For soft atoms, there is a strong, $`T`$-dependent increase of the former and a corresponding, $`T`$-dependent decrease of the latter. The same, but less $`T`$-independent, changes occur for hard atoms.

Summarizing the results presented so far, it can be said that for the three harmonic systems studied dynamical heterogeneities are clearly observed; the most pronounced overall correlation increases are for the binary mixture and for soft spheres, while soft inhomogeneities are observed to depend on temperature for the full LJ and binary mixture. In all cases, the most pronounced effects concern the low-$`R`$ tail of the pair correlation function.

One important question is whether with the present system size ($``$ 2000 atoms) it is possible to estimate the expected size, $`\lambda `$, of the vibrational dynamical inhomogeneities. To check this point, in the case of full LJ we evaluated $`\lambda `$ for samples with different number of atoms, $`N_0=500,1000,1500,2048`$. In order to evaluate $`\lambda `$, the height of the first point of $`g_N(R)`$ (marked by an arrow in Fig. 2) was evaluated as a function of $`N`$ and compared to the corresponding height in $`g_W`$; the value of $`N`$ at which the two become equal is taken as an estimate of the number of atoms which constitute the inhomogeneity and hence of its size. The results are reported in Fig. 5 and show that certainly for hard atoms, and probably for soft ones as well at high temperature, the heterogeneity size is still growing with $`N_0`$ at the maximum sample dimension available to us, which seems to indicate that the ”true” size is larger than we can determine. One possible exception to this is the low-temparature ($`T=`$0.01 and 0.1) soft atoms, for which there appears to be a flattening of the curve in Fig. 5 starting at $`N_0`$=1000. However, before drawing definite conclusions we think it is wiser to await for more statistics (only 1 sample was examined for $`N_0=`$500 and 1000, and 3 samples for $`N_0`$=1500). However, if the low-$`T`$ result of Fig. 5 should be confirmed, temperature would play a role also in determining the soft heterogeneity size, and the same should be expected for the binary mixture (see Fig. 4(a) ).

In summary, it is the first time that the existence of $`T`$-dependent vibrational dynamical heterogeneities in glasses is clearly shown to exist by simulation with LJ or soft spheres systems. However, it is not possible to draw definite and general conclusions about the heterogeneity size and its relation with the boson peak . For that, the number of atoms in the samples studied by simulation should be much larger than 2000; work in this direction is in progress. We also plan to use other potentials taking into account long-range interactions, like for example the BKS one  which is suitable for SiO2. We expect a more contrasted inhomogeneous elasticity or vibrational inhomogeneity, like in vitreous silica, possibly in relation with a more intense boson peak .


text:
# Constraints on primordial black holes and primeval density perturbations from the epoch of reionization

## 1. Introduction

Primordial black holes (PBHs) form from the gravitational collapse when an initially dense region on super-horizon scales enters the horizon during the early radiation-dominated era (Carr & Hawking 1974, Carr 1975). In the inflationary universe, the abundance of PBHs is extremely sensitive to the spectral index $`n`$, which specifies the dependence of the power spectrum of primordial density fluctuations on the comoving wavenumber $`k`$, $`|\delta _k|^2k^n`$. Thus, a constraint on the abundance of PBHs will also yield a constraint on the index $`n`$.

The simplest constraint on the abundance of PBHs is given by the condition $`\mathrm{\Omega }_{bh}1`$, i.e., the mass density of the PBHs should be less than the critical density of the universe. This condition results in $`n<1.5`$ regardless of whether a black hole can evaporate completely, or ends with a Planck relic (Kim & Lee 1996). More effective constraints can be derived by comparing the particle emission from PBH evaporation with observed background. For instance, the cosmic background of $`\gamma `$-ray spectrum on $`100`$ Mev places a tightest constraint on the PBH density parameter $`\mathrm{\Omega }_{bh}10^8`$ (MacGibbon & Carr 1991), which yields the index $`n<1.25`$ (Kim, Lee & MacGibbon 1999). However, these constraints depend on an additional parameter, the reheating temperature $`T_{rh}`$, or equivalently, the minimum mass of PBHs $`M_{min}`$, which is not well determined.

Recently, the epoch of the onset of the cosmic reionization is explored with absorption spectrum of high redshift ($`5.7<z<6.3`$) SDSS quasars (Fan et al. 2001, and references therein). The existence of the complete Ly$`\alpha `$ Gunn-Peterson trough in the absorption spectra of high redshift quasars indicates that the dark age probably ended at the redshift $`z68`$. The ionizing photon background in the IGM at $`z6`$ is found to be more than 20 times lower than that at $`z3`$. These results motivate us to explore the constraints on the abundance of PBHs, and then on the primordial density perturbations with the low energy photon background at the epoch of reionization. We will show that this constraint is basically independent of the uncertain parameters $`T_{rh}`$ or $`M_{min}`$. PBHs have been considered as a possible source of the cosmic reionization (Gibilisco 1996, 1997), but in these works, the mass function of PBHs is treated as free parameters. It cannot be used to constrain $`n`$.

This Letter is organized as follows. In §2, we briefly present the methods of calculating the mass function of the PBHs and the ionizing photon background from PBH’s evaporation. The numerical results of the constraints on $`n`$ from the observed high redshift ionizing background are given in §3. §4 are discussions and conclusions.

## 2. Methods

### 2.1. Mass function of PBHs

For Gaussian primeval density perturbations, the mass function $`n_{bh}(M_{bh})`$ of PBHs at redshift $`z`$ can be calculated by the same way as calculating the mass function of collapsed massive halos with the Press-Schechter formalism (Press & Schechter 1974). Accordingly, the number density of regions with density contrast from $`\delta `$ to $`\delta +d\delta `$ and mass $`M`$ to $`M+dM`$ at the time $`t_i`$ of the onset of radiation-dominated epoch is

$`n(M,\delta )dMd\delta =`$
$`\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\rho _i}{M^2}}{\displaystyle \frac{1}{\sigma _R}}\left|{\displaystyle \frac{\mathrm{ln}\sigma _R}{\mathrm{ln}M}}\right|\left|{\displaystyle \frac{\delta ^2}{\sigma _R^2}}1\right|e^{\frac{\delta ^2}{2\sigma _R^2}}dMd\delta ,`$
where $`\rho _i`$ is the cosmic mean mass density at $`t_i`$, $`M`$ the mass within $`R`$, and $`\sigma _R`$ the variance of the probability distribution function of initial Gaussian mass field smoothed by scale $`R`$. Since the mass function is normalized at $`\rho _i`$, i.e. $`Mn(M,\delta )𝑑M𝑑\delta =\rho _i`$, eq.(1) gives the physical number density of region $`(M,\delta )`$ at $`t_i`$. The comoving number density is $`n(M,\delta )_{com}=n(M,\delta )(a_i/a_0)^3`$, where $`a(t)`$ is the cosmic factor.

It has been shown (Carr 1975) that when a region $`(M,\delta )`$ enters the horizon at $`t_h`$, i.e. when $`M`$ is about the same as horizon mass at $`t_h`$, this region will collapse to a black hole if the initial density contrast $`\delta `$ satisfies

$$\gamma (M/M_i)^{2/3}\delta (M/M_i)^{2/3},$$
(2)
where parameter $`\gamma 1/3`$ and $`M_i`$ is the horizon mass at $`t_i`$, i.e. the horizon mass at the beginning of the radiation-dominated. The mass of a black hole formed from this region is equal to $`M_{bh}=\gamma ^{1/2}M_i^{1/3}M^{2/3}`$. Thus, the mass function of PBHs is

$`n_{bh}(M_{bh})={\displaystyle \frac{dM}{dM_{bh}}}{\displaystyle _{\gamma (M/M_i)^{2/3}}^{(M/M_i)^{2/3}}}n(M,\delta )𝑑\delta =`$
$`{\displaystyle \frac{3}{2}}\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\gamma ^{3/4}\rho _iM_i^{1/2}}{m_{bh}^{5/2}}}{\displaystyle _\gamma ^1}\left|{\displaystyle \frac{\mathrm{ln}\sigma _R}{\mathrm{ln}M}}\right|\left|{\displaystyle \frac{\delta _h^2}{\sigma _h^2}}1\right|e^{\frac{\delta _h^2}{2\sigma _h^2}}{\displaystyle \frac{1}{\sigma _h}}𝑑\delta _h,`$
where $`\delta _h=\delta (M/M_i)^{2/3}`$ is the density contrast when the region $`(M,\delta )`$ crosses the horizon at $`t_h`$, and variance $`\sigma _h=\sigma _R(M/M_i)^{2/3}`$. Since $`MM_i`$, the minimum mass of the PBH is $`\gamma ^{1/2}M_i`$.

For a power law initial density perturbation with spectral index $`n`$, i.e. $`\sigma _hM^{(1n)/6}`$, we have $`\mathrm{ln}\sigma _R/\mathrm{ln}M=(n+3)/6`$, and eq.(3) becomes (Kim & Lee 1996)

$`n_{bh}(M_{bh})=`$
$`{\displaystyle \frac{n+3}{4}}\left({\displaystyle \frac{2}{\pi }}\right)^{1/2}\gamma ^{7/4}\rho _iM_i^{1/2}M_{bh}^{5/2}\sigma _h^1e^{\gamma ^2/2\sigma _h^2}.`$
As usual, the variance $`\sigma _h`$ will be normalized with the COBE observation (e.g., Fang & Xu 1999, Bugaev & Konishchev 2001).

In eq.(4), only the factors $`\rho _i=3/32\pi Gt_i^2`$ and $`M_i=(4\pi /3)(ct_i)^3\rho _i`$ are dependent on $`t_i`$. Thus, $`n_{bh}(M_{bh})t_i^{3/2}`$. On the other hand, $`(a_i/a_0)^3(1+z_{eq})^3(t_i/t_{eq})^{3/2}`$, where $`t_{eq}`$ and $`z_{eq}`$ are, respectively, the time and the redshift at the end of the radiation-dominated era. Thus, the comoving mass function $`(a_i/a_0)^3n_{bh}(M_{bh})`$ is independent of $`t_i`$.

### 2.2. Ionizing photon background from PBH’s evaporation

For a Schwarzschild black hole $`m`$, the Hawking emission rate for particles with spin $`s`$, energy from $`E`$ to $`E+dE`$ and per degree of particle freedom is

$$\frac{dN}{dt}=\frac{1}{2\pi \mathrm{}}\frac{\mathrm{\Gamma }_sdE}{e^{E/k_BT(m)}(1)^{2s}},$$
(5)
where $`T(m)=\mathrm{}c^3/8\pi Gk_Bm`$ is the temperature of black hole $`m`$, and the dimensionless parameter $`\mathrm{\Gamma }_s`$ is the absorption probability of spin $`s`$ particle by a black hole $`m`$ (MacGibbon & Webber 1990; Page 1975). For low energy photon emission, i.e., $`h\nu k_bT(m)`$, $`\mathrm{\Gamma }=128G^4m^4\nu ^4/3c^{12}`$, and eq.(5) gives

$$\frac{d^2N}{d\nu dt}=f(\nu ,m)=\frac{128G^4m^4\nu ^4}{3c^{12}[e^{h\nu /k_BT(m)}1]}.$$
(6)
The evaporation leads to the decreasing of black hole mass with the rate

$$\frac{dm}{dt}=\frac{\alpha (m)}{m^2}$$
(7)
where $`\alpha (m)`$ takes into account all degrees of freedom of evaporated particles (MacGibbon 1991). Since $`\alpha (m)`$ is not a strong function of $`m`$, the time dependence of $`m`$ for a black hole with initial mass $`M_{bh}`$ at $`t_h`$ is approximated by

$$m(M_{bh},t)[M_{bh}^33\alpha (m)t]^{1/3},$$
(8)
if $`tt_h`$.

Thus, the background photon energy flux $`J(\nu ,z_{obs})`$ per unit frequency interval and per solid angle at redshift $`z_{obs}`$ is as follows:

$`{\displaystyle \frac{J(\nu ,z_{obs})}{h\nu }}={\displaystyle \frac{c}{4\pi }}{\displaystyle _{t_{eq}}^{t(z_{obs})}}𝑑t\left({\displaystyle \frac{a_{obs}}{a}}\right)\left({\displaystyle \frac{a_i}{a_{obs}}}\right)^3`$
$`{\displaystyle _{M_{min}}^{M_{max}}}𝑑M_{bh}n_{bh}(M_{bh})\mathrm{exp}[\tau _{eff}(\nu ,z_{obs},z)]`$
$`f[\nu (1+z)/(1+z_{obs}),m(M_{bh},t)],`$
$`M_{min}`$ is given by $`\mathrm{max}[\gamma ^{1/2}M_i,m(t)]`$, where $`m(t)`$ is the solution of $`m(t)=\{3\alpha [m(t)]t\}^{1/3}`$. The integral upper limit $`M_{max}`$ in eq.(9) in not important, as the mass function $`n_{bh}(M_{bh})`$ is rapidly decreasing with the increasing of $`M_{bh}`$. $`\tau _{eff}(\nu ,z_{obs},z)`$ is the opacity for photons emitted at $`z`$ and observed at $`z_{obs}`$ at the frequency $`\nu `$.

## 3. Results

Using eq.(9), we calculate the ionizing photon background at redshift $`z_{obs}6`$ \- 8. The opacity $`\tau _{eff}`$ for photons observed at the hydrogen Lyman edge ($`h\nu _0`$ =13.6 ev) is dominated by HI photoelectric absorption. We will use the $`\tau _{eff}`$ estimated by the so-called high absorption model (Meiksin & Madau, 1993). It is

$`\tau (\nu _0,z_{obs},z)`$
$`0.0097(1+z_{obs})^{1.56}[(1+z)^{1.84}(1+z_{obs})^{1.84}]`$
$`0.0068(1+z_{obs})^{1.56}[(1+z)^{0.4}(1+z_{obs})^{0.4}]`$
$`8.06\times 10^5[(1+z)^{3.4}(1+z_{obs})^{3.4}].`$
The optical depth given by eq.(10) is very short. For $`z_{obs}=6`$, the mean transmission $`e^{\tau _{eff}}`$ is less than $`10^4`$ at $`z6.5`$. That is, the effective flux $`J(\nu _0)`$ observed at $`z=6`$ is only contributed by PBHs localized within the spatial range from $`z=6`$ to 6.5. For $`z=8`$, this range is about $`z=8`$ to 8.25. Actually, eq.(10) is found from a fitting with observations at $`z3`$ \- 4. However, there is no direct observation of the hydrogen Lyman edge absorption at redshifts $`z>6`$. Nevertheless, eq.(10) is consistent with the Ly$`\alpha `$ transmission $`e^\tau =0.004\pm 0.003`$ in the redshift range $`[5.956.15]`$ (Becker et al. 2001). Hence, it would be reasonable to estimate the optical depth at $`z=6`$ \- 6.5 by eq.(10).

The latest observation of the complete Gunn-Peterson troughs in the absorption spectrum of a $`z`$=6.28 quasar indicates that the photoionization rate at $`z`$6 is lower than 8.0$`\times 10^{14}`$s-1 for Ly$`\alpha `$, or 2.0$`\times 10^{14}`$s-1 for Ly$`\gamma `$ (Fan et al. 2001). Assuming the ionizing photons are from massive stars $`J(\nu )\nu ^5`$ (Barkana & Loeb 2001), one can then find that the upper limits to the ionizing photon background $`J_{21}J(\nu )/10^{21}`$ ($`J(\nu )`$ in unit of erg cm-2 sec-1Hz-1sr-1) at photon energy 13.6ev are, respectively, $`5.3\times 10^2`$ and $`1.3\times 10^2`$ for Ly$`\alpha `$, and Ly$`\gamma `$, and at 24ev, are $`3.1\times 10^3`$, $`7.8\times 10^4`$ for Ly$`\alpha `$, and Ly$`\gamma `$, respectively.

Figure 1 presents the PBH contributed ionizing photon background $`J_{21}`$ at $`z=6`$ as a function of the index $`n`$ of the primordial fluctuation power spectrum. The upper limits to the ionizing photon background yield $`n<1.273`$ for $`E`$=13.6ev, and $`n<1.270`$ for $`E`$=24 ev.

It should be pointed out that these constraint basically are independent of $`t_i`$. Eq.(9) shows that $`J(\nu _0)`$ depends on $`t_i`$ only via $`M_{min}`$, the truncation of PBH mass function at small mass side. On the other hand, the optical depth at $`z=6`$ is very small. Only PBHs, which can survive to $`z6`$, has contribution to $`J(\nu _0)`$. Thus, $`J(\nu _0)`$ is independent of $`M_{min}`$, if $`M_{min}10^{12}`$ g, as those PBHs have already evaporated at $`z6`$ \[eqs.(7) and (8)\]. In inflationary universe, $`t_i=0.301g_{}^{1/2}M_{pl}/T_{rh}^2`$ (Kolb & Turner 1990). Thus, the $`M_{min}`$-independence leads to also $`T_{rh}`$-independence.

The $`M_{min}`$ or $`T_{rh}`$-independence of $`J(\nu _0)`$ can be also seen from Fig. 2, which plots $`dJ(\nu _0,z_{obs})/dm`$ vs. $`m`$. It shows that $`J(\nu _0)`$ at $`z=6`$ is mainly contributed by the evaporation of PBHs with mass $`m2\times 10^{14}`$ g. Therefore, $`J(\nu _0)`$ is independent of $`M_{min}`$ if it is less than $`2\times 10^{14}`$ g. The peak $`2\times 10^{14}`$ g. is also consistent with the assumption of $`h\nu k_bT(m)`$ used in eq.(6).

We calculate the flux $`J(\nu ,z)`$ at $`z=6`$ in the frequency range from $`\nu `$ = 3$`\times 10^{15}`$ to $`6\times 10^{15}`$, which corresponds to photon energy range from 13.6 to 24.6 ev (HeI ionization). In this frequency interval one can still approximately use the optical depth $`\tau _{eff}`$ of eq.(10). The result is shown in Fig.3. The flux $`J(\nu )`$ can be well approximated by a power law spectrum $`J(\nu )\nu ^3`$. This spectrum is very different from the models of massive stars ($`J(\nu )\nu ^5`$) and quasars ($`J(\nu )\nu ^{1.8}`$, see Madau, Haardt, & Rees 1999), for both of which, the flux decreases with the increasing of $`\nu `$. This is because $`J(\nu )`$ are mainly from the low energy tail of the photon emission by 1014g black holes, whose temperature $`T(m)h\nu _0/k_b`$. In the low energy tail, the intensity is increasing with $`\nu `$.

In all the above calculations, the photons produced from the fragmentation of the PBH-emitted quarks and gluons are not considered (MacGibbon 1991). Since $`J(\nu _0)`$ is given by local PBHs, we need only to consider the fragmentation of quarks and gluons into low energy ($`10`$ \- 20 ev) photons. This problem is still quite unclear and uncertain. Since the fragmentation effect will lead to an increase of $`J(\nu _0)`$. Therefore, all above-mentioned upper limits on $`n`$ will be held.

## 4. Discussions and Conclusions

In the standard inflationary model with a power-law power spectrum of the primeval density perturbations, we calculated the ionizing background flux of the PBHs. With the latest observed ionizing photon background flux at the redshift of reionization, we find that the spectral index of the primeval density perturbations should be $`n<1.27`$, assuming that the background photon flux are emitted by massive stars. The difference is negligible if the background flux coming from quasars. If the parameter $`\gamma `$ is taken to be 0.6 (Niemeyer & Jedamzik 1998), we have $`n<`$1.29. We also calculate $`J(\nu )`$ at $`z=8`$. The result is almost the same as Fig. 3. However, there is no available data on the $`J(\nu _0)`$ at this redshift. If the reionization is really onset at $`z6`$ \- 8, $`J(\nu _0)`$ at $`z=8`$ would be significantly less than that at $`z=6`$, one may have more tight constraint. All these constraints are independent of the reheating temperature $`T_{rh}`$, if $`T_{rh}>10^{10}`$ Gev. The abundance of black hole given by this constraint is $`\mathrm{\Omega }_{bh}<10^4`$, which is much stronger than the critical density constraint $`\mathrm{\Omega }_{bh}<1`$.

The current observations of the angular power spectrum of the temperature fluctuations of cosmic microwave background radiation (CMB) are found to be able to fit with a $`n=1`$ power spectrum on angular scales from tens of degrees to sub-degrees. The PBH’s constraint $`n<1.27`$ shows that primeval density perturbations might still be scale invariant on scales much less than the direct observation of the CMB. Moreover, in the CMB data fitting, the parameters of the spectral index $`n`$ and the optical depth of the CMB photons to the last scattering surface $`\tau _c`$ are degenerated (Stompor et al. 2001). That is, even when the fitted result is $`n_{fitting}1`$, the original index may be $`n>1`$ if $`\tau _c>0`$. On the other hand, the constraint on $`n`$ given by the PBH is independent of the optical depth $`\tau _c`$.

Thanks to Hongguang Bi for stimulating discussions. PH thanks M. Gibilisco for sending him her papers and helpful discussions. PH is supported by a Fellowship of the World Laboratory.


text:
# 1 Introduction

## 1 Introduction

QCD undergoes a transition from a confined hadronic phase to a deconfined partonic phase, the quark-gluon plasma (QGP), at a temperature of $`T_c150170`$ MeV . A central quantity of matter in thermal equilibrium is the Helmholtz free energy from which the pressure $`p`$, energy density $`ϵ`$ and entropy density $`s`$ – which are important ingredients for the description of ultra-relativistic heavy-ion collisions – can be derived. A first-principles understanding of the equation of state (EOS) of hot QCD is therefore of great interest, also in order to reliably identify and calculate experimental signatures of that elusive state. Perturbative results on the EOS are available up to order $`𝒪(g_s^5)`$ . However, for temperatures of interest the strong coupling constant is presumably large: $`g_s12`$. The perturbative expansion in powers of $`g_s`$ shows bad convergence already for much smaller values of the coupling. Furthermore, in the vicinity of a phase transition, perturbative methods are in general not expected to be applicable. Non-perturbative methods such as lattice QCD calculations become mandatory. From these numerical simulations the EOS of a pure gluon plasma is known to high accuracy , and there are first estimates for the continuum EOS of systems including quarks, albeit still with unphysically large masses .

Various interpretations of the lattice data have been attempted, most prominently as the EOS of a gas of quark and gluon quasiparticles. In a phenomenological framework, quarks and gluons are simply treated as non-interacting, massive quasiparticles . More recently, a quasiparticle description of QCD thermodynamics has been derived in a more rigorous treatment using resummed hard thermal loop (HTL) perturbation theory . Employing the full HTL spectral representions, the resulting EOS can be matched to lattice data down to temperatures $`T3T_c`$; below that temperature, non-perturbative physics not amenable in an expansion in $`g_s`$ becomes important. Unfortunately, as evident from figure 1 (left panel), current experiments only probe that very temperature regime where the underlying physics, the confinement and chiral symmetry breaking mechanism, is not sufficiently understood. Phenomenological models incorporating as much physics as is known are therefore necessary. Here, we propose a new quasiparticle model of the QGP that incorporates a parametrization of confinement close to $`T_c`$, supplemented by thermal masses compatible with lattice results. Details can be found in .

## 2 Quasiparticles and confinement

The use of a quasiparticle model in QCD is based on the observation that in a strongly interacting system, the complex dynamics often rearranges itself in such a way that gross features of the physics can be described in terms of appropriate effective degrees of freedom. Consider a SU(3) gluon plasma. From asymptotic freedom, we expect that at very high temperatures the plasma consists of quasifree gluons. As long as the spectral function of the thermal excitations at lower temperatures resembles qualitatively this asymptotic form, a gluonic quasiparticle description is expected to be applicable. Their dispersion relation reads

$$\omega _k^2k^2+m_{}^2(T),$$
(1)
where $`m_{}(T)`$ acts as an effective mass generated dynamically by the interaction of the gluons with the heat bath background.

However, the picture of a simple massive gas is presumably not appropriate close to $`T_c`$ because the driving force of the transition, the confinement process, is not taken into account. This physics has to be inserted by hand. Below $`T_c`$, the relevant degrees of freedom in SU(3) gauge theory are heavy, colour singlet glueballs. Approaching $`T_c`$, deconfinement sets in and the gluons are liberated, followed by a sudden increase in entropy and energy density. Conversely, when approaching the phase transition from above, the number of thermally active degrees of freedom is reduced due to the onset of confinement. As $`T`$ comes closer to $`T_c`$, an increasing number of gluons gets trapped in glueballs which disappear from the thermal spectrum: since $`m_{GB}1.5`$ GeV and $`T_c270`$ MeV (for pure gauge theory), glueballs are simply too heavy to become thermally excited in the temperature range under consideration (up to about $`5T_c`$). So, while the confinement mechanism as such is still not understood, it is not necessary to know it in detail since we consider a statistical system. All confinement does on a large scale is to cut down the number of thermally active gluons as the temperature is lowered. This effect can be included in the quasiparticle picture by modifying the distribution function of the gluons by a temperature-dependent confinement factor $`C(T)`$: $`f_B(\omega _k)C(T)f_B(\omega _k)`$.

To become quantitative, we have to specify the thermal masses $`m_{}(T)`$ entering (1). Based on the observation that the Debye screening mass $`m_D`$ evaluated on the lattice shows approximate critical behaviour , in accordance with a weakly first order phase transition and in contrast to perturbative results, we parametrize

$$m_{}(T)G_0T\left(1\frac{T_c}{T}\right)^\beta .$$
(2)
Consider now the entropy of a gas of massive gluons with such a dropping $`m_{}(T)`$. The result for $`s(T)`$ will clearly overshoot the lattice entropy because light masses near $`T_c`$ lead to an increase in $`s(T)`$. However, since the entropy is a measure for the number of active degrees of freedom, the difference may be accounted for by the aforementioned confinement process as it develops when the temperature is lowered towards $`T_c`$. The explicit temperature dependence of the confinement factor $`C(T)`$ can be obtained simply as the ratio of the lattice entropy and the entropy calculated with a dropping input gluon mass, and $`C(T)`$ again shows near-critical behaviour: $`C(T)(1T_c/T)^\gamma `$.

Thermodynamical quantities like the energy density can now be calculated as

$$ϵ(T)=\frac{N_c^21}{2\pi ^2}\underset{0}{\overset{\mathrm{}}{}}𝑑kk^2\left[C(T)f_B(\omega _k)\right]\omega _k+B(T).$$
(3)
The function $`B(T)`$ is not an independent function, but uniquely determined by $`m_{}(T)`$, $`C(T)`$ and their $`T`$-derivatives. It is necessary to maintain thermodynamical self-consistency and can be interpreted as the thermal energy density of the vacuum. The explicit expressions for $`B(T)`$, the pressure and the entropy density can be found in . In figure 1 (right panel), we compare $`ϵ`$, $`s`$ and $`p`$ with lattice data, as a function of $`T`$. We achieve a good and economic parametrization. The so-called interaction measure, related to the trace of the energy-momentum tensor, $`T_\mu ^\mu =ϵ3p`$, is also nicely reproduced. Figure 2 (left panel) shows the function $`B(T)`$ and the spacelike plaquette expectation value $`\mathrm{\Delta }_\sigma `$, as measured on the lattice, that can be related to the thermal chromomagnetic condensate $`𝐁^2_T`$.

We find the simple relation

$$B(T)=\frac{1}{2}\mathrm{\Delta }_\sigma (T)T^4=\frac{11\alpha _s}{8\pi }𝐁^2_T+\frac{1}{4}G^2_{T=0}$$

with the zero-temperature condensate $`G^2_{T=0}`$. This relation between $`B(T)`$ and $`𝐁^2_T`$ may hint at a deeper connection between $`B(T)`$ as a carrier of non-perturbative effects, and the magnetic condensate. After all, $`B(T)`$ represents the thermal energy of the (non-trivial) Yang-Mills vacuum. Note that there is no such relation in other quasiparticle models, as evident from the dashed line in figure 2 (left panel).

## 3 Dynamical quarks

The extension of the mechanism presented so far to systems with dynamical quarks is not straightforward. Simulations of fermions on the lattice are still plagued by problems. However, when plotting the lattice pressure, normalized to the ideal gas value, for the pure gauge system and for systems with 2, 2+1 and 3 quark flavours, it is found that the QCD EoS shows a remarkable flavour independence when plotted against $`T/T_c`$. The flavour dependence is then well approximated by a term reminiscent of an ideal gas $`p(T,N_f)\left(16+10.5N_f\right)\stackrel{~}{p}(T/T_c)`$ with a universal function $`\stackrel{~}{p}(T/T_c)`$ . This hints at a confinement mechanism being only weakly flavour-dependent, and hence we assume that the function $`C(T)`$ acts in a universal way on quarks and gluons. Figure 2 (right panel) shows pressure, energy and entropy density for two massless quark flavours. Reassuringly, the pressure of the confinement model lies well within a narrow lattice estimate for the continuum EOS in the chiral limit . The picture for two light and a heavier strange quark looks similar although the approach to the ideal gas limit is slower.

The neglect of the confined particles now violates detailed balance more severely since the lightest Goldstone modes have masses comparable to $`T_c`$, hence they do contribute to the thermodynamics. Within a statistical hadronization model it may be possible to estimate their contribution. Furthermore, the calculation of susceptibilities will enable us to gain further insight in the nature of $`C(T)`$. Work along these lines is in progress.

## 4 Dilepton rates

We apply the proposed model now to dilepton production in URHIC where it enters in two ways: the EOS serves as input for the construction of the fireball that is discussed in . Second, the dilepton emissivity of a static hot spot,

$$\frac{d𝒩_{ee}}{d^4xd\omega d^3k}\frac{\rho _V(\omega ,k;T)}{\mathrm{exp}(\omega /T)1},$$
(4)
is proportional to the vector spectral function $`\rho _V`$ that is calculated from the photon self energy. By construction, the quasiparticle $`q\overline{q}`$-loop is the only contribution. In addition, the confinement mechanism reduces the probability of thermally exciting a quark by $`C(T)`$ and thus the total rate by $`C(T)^2`$. The corresponding $`\rho _V`$ agrees with qualitative lattice results .

To compare with experiment, the rate (4) has to be convoluted with the fireball expansion. A discussion of the fireball dynamics and the final dilepton rates can be found in T. Renk’s proceedings contribution .