text:
# 1. Introduction

## 1. Introduction

The idea of a spontaneously broken phase as a ‘condensate’ is now widely accepted. For instance, in the physically relevant case of the Standard Model, the situation can be summarized saying  that ”What we experience as empty space is nothing but the configuration of the Higgs field that has the lowest possible energy. If we move from field jargon to particle jargon, this means that empty space is actually filled with Higgs particles. They have Bose condensed.”

Here, by ‘Bose condensation’ one means the phenomenon of spontaneous particle creation in the same quantum state ($`𝐤=0`$ in some reference frame) from the empty vacuum $`|o`$ of perturbation theory. In this way, the translation from ‘field jargon to particle jargon’, amounts to establish a well defined functional relation (see ref. and Sect.2) $`n=n(\varphi ^2)`$ between the average density of scalar quanta, the ‘phions’, and the average value of the scalar field $`\mathrm{\Phi }=\varphi `$. Thus, Bose condensation is just a consequence of the minimization condition of the effective potential $`V_{\mathrm{eff}}(\varphi )`$. This has absolute minima at some values $`\varphi =\pm v0`$ for which $`n(v^2)=\overline{n}0`$ .

Of course, in order this picture to be consistent, spontaneus symmetry breaking should occur for physical (i.e. real and non-negative) values of the phion mass $`m_\mathrm{\Phi }`$. In other words, Bose condensation requires the phase transition in (cutoff) $`\lambda \mathrm{\Phi }^4`$ theories to be first order. While in the presence of gauge bosons this can be shown perturbatively , the use of perturbation theory in a pure $`\lambda \mathrm{\Phi }^4`$ theory leads to contradictory results between even and odd orders . Therefore, one has to analyze the effective potential beyond perturbation theory. By relying on the assumed exact ‘triviality’ property of the theory in 3+1 space-time dimensions , one is driven to consider the class of ‘triviality compatible’ approximations  to the effective potential, say $`V_{\mathrm{eff}}(\varphi )=V_{\mathrm{triv}}(\varphi )`$. This includes the one-loop potential, the gaussian and post-gaussian  calculations of the Cornwall-Jackiw-Tomboulis  effective potential for composite operators, i.e. all approximations where the fluctuation field is governed by an effective quadratic Hamiltonian. In all such cases, the lowest energy state of the massless theory at $`m_\mathrm{\Phi }=0`$ corresponds to a broken-symmetry phase as first shown in ref. for the gaussian approximation. Therefore the phase transition, occurring earlier for small positive values of $`m_\mathrm{\Phi }`$, is first order. However, there is a subtlety since it corresponds to an infinitesimally weak first-order phase transition. In fact, it is first order for any finite ultraviolet cutoff $`\mathrm{\Lambda }`$ but becomes asymptotically second order  in the continuum limit $`\mathrm{\Lambda }\mathrm{}`$.

In any case, for all finite values of $`\mathrm{\Lambda }`$, a particle-gas picture of the underlying scalar system is possible. This can be formulated in terms of two basic quantities: the equilibrium phion density $`\overline{n}`$ and the phion-phion scattering length $`a`$

$$a=\frac{\lambda }{8\pi m_\mathrm{\Phi }}$$
(1)
defined, in the limit of zero-momentum scattering, from the dimensionless scalar self-coupling $`\lambda `$ and the phion mass.

As shown in ref., these two quantities combine to produce all relevant length scales of the system

$$a\frac{1}{\sqrt{\overline{n}a}}\frac{1}{\overline{n}a^2}$$
(2)
that decouple for an infinitely dilute system where $`\overline{n}a^30`$. In this situation, which corresponds to approaching the continuum limit of quantum field theory , one discovers an unexpected result: the hierarchical nature of the scalar condensate.

At the same time, a particle-gas picture of the broken-symmetry phase raises several questions. For instance, in condensed media the properties of the system over vastly different scales are described by different branches of the energy spectrum. Are there similar transitions in the scalar condensate ? Also, what about the coexistence of exact Lorentz covariance and vacuum condensation in effective quantum field theories ? Can the violations of locality, at the energy scale fixed by the ultraviolet cutoff, induce non-Lorentz-covariant modifications of the infrared energy spectrum that depends on the vacuum structure  ?

To indicate this type of infrared-ultraviolet connection, originating from vacuum condensation in effective quantum field theories, Volovik  has introduced a very appropriate name: reentrant violations of special relativity in the low-energy corner. These occur in a small shell of three-momenta, say $`|𝐤|<\delta `$, that only vanishes in the strict local limit where $`\mathrm{\Lambda }\mathrm{}`$ and an exact Lorentz-covariant energy spectrum is re-obtained in the whole range of momenta.

By denoting $`M_H`$ as the typical energy scale associated with the Lorentz-covariant part of the energy spectrum, the ‘reentrant’ nature of the violating effects means that $`𝒪(\frac{\delta }{M_H})`$ vacuum-dependent corrections are equivalent to $`𝒪(\frac{M_H}{\mathrm{\Lambda }})`$ effects (see below). The $`1/\mathrm{\Lambda }`$ terms, that have always been neglected when discussing  how Lorentz covariance emerges in effective theories when removing the ultraviolet cutoff, although infinitesimally small, can play an important role over distances larger than $`\mathrm{\Lambda }/M_H^2`$, i.e. infinitely larger than the typical elementary particle scale $`\xi _H=1/M_H`$.

As discussed in ref., the basic ingredient to understand the nature of the ‘reentrant’ effects in the scalar condensate consists in a purely quantum-field-theoretical result: the two-valued nature of the connected zero-four-momentum propagator $`G^1(k=0)`$ in the broken phase . In fact, besides the well known massive solution $`G_a^1(k=0)=M_H^2`$, one also finds $`G_b^1(k=0)=0`$.

The b-type of solution corresponds to processes where assorbing (or emitting) a very small 3-momentum $`𝐤0`$ does not cost any finite energy. This situation is well known in a condensed medium, where a small 3-momentum can be coherently distributed among a large number of elementary constituents, and corresponds to the hydrodynamical regime of density fluctuations whose wavelengths $`2\pi /|𝐤|`$ are larger than $`r_{\mathrm{mfp}}`$, the mean free path for the elementary constituents.

In this sense, the situation is very similar to superfluid 4He, where the observed energy spectrum is due to the peculiar transition from the ‘phonon branch’ to the ‘roton branch’ at a momentum scale $`|𝐤_o|`$ where $`E_{\mathrm{phonon}}(𝐤_o)E_{\mathrm{roton}}(𝐤_o)`$. The analog for the scalar condensate amounts to an energy spectrum with the following limiting behaviours :

i) $`E(𝐤)E_b(𝐤)c_s|𝐤|`$ for $`𝐤0`$

ii) $`E(𝐤)E_a(𝐤)M_H+\frac{𝐤^2}{2M_H}`$ for $`|𝐤|\delta `$

where the characteristic momentum scale $`\delta M_H`$, at which $`E_a(\delta )E_b(\delta )`$, marks the transition from collective to single-particle excitations. This occurs for

$$\delta 1/r_{\mathrm{mfp}}$$
(3)
where

$$r_{\mathrm{mfp}}\frac{1}{\overline{n}a^2}$$
(4)
is the phion mean free path, for a given value of the phion density $`n=\overline{n}`$ and a given value of the phion-phion scattering length $`a`$. In terms of the same quantities, one also finds (see  and Sect.2)

$$M_H^2\overline{n}a$$
(5)
giving the anticipated trend of the dimensionless ratios ($`\mathrm{\Lambda }1/a`$)

$$\frac{\delta }{M_H}\frac{M_H}{\mathrm{\Lambda }}\sqrt{\overline{n}a^3}0$$
(6)
in the continuum limit where $`a0`$ and the mass scale $`\overline{n}a`$ is held fixed .

Now, deducing the detailed form of the energy spectrum that interpolates between $`E_a(𝐤)`$ and $`E_b(𝐤)`$ is a formidable task. To have an idea, the same problem in superfluid 4He, after more than fifty years and despite the efforts of many theorists, notably Landau and Feynman, has not been solved in a satisfactory way. Therefore, by taking into account the above remark, one can simply approximate  the expansion of the scalar field in the broken phase by two separate branches as (phys=‘physical’)

$$\mathrm{\Phi }_{\mathrm{phys}}(x)=v_R+h(x)+H(x)$$
(7)
with

$$h(x)=\underset{|𝐤|<\delta }{}\frac{1}{\sqrt{2𝖵E_k}}\left[\stackrel{~}{h}_𝐤\mathrm{e}^{i(𝐤.𝐱E_kt)}+(\stackrel{~}{h}_𝐤)^{}\mathrm{e}^{i(𝐤.𝐱E_kt)}\right]$$
(8)
and

$$H(x)=\underset{|𝐤|>\delta }{}\frac{1}{\sqrt{2𝖵E_k}}\left[\stackrel{~}{H}_𝐤\mathrm{e}^{i(𝐤.𝐱E_kt)}+(\stackrel{~}{H}_𝐤)^{}\mathrm{e}^{i(𝐤.𝐱E_kt)}\right]$$
(9)
where $`𝖵`$ is the quantization volume and $`E_k=c_s|𝐤|`$ for $`|𝐤|<\delta `$ while $`E_k=\sqrt{𝐤^2+M_H^2}`$ for $`|𝐤|>\delta `$. Also, $`c_s\delta M_H`$.

Eqs.(7)-(9) replace the more conventional relations

$$\mathrm{\Phi }_{\mathrm{phys}}(x)=v_R+H(x)$$
(10)
where

$$H(x)=\underset{𝐤}{}\frac{1}{\sqrt{2𝖵E_k}}\left[\stackrel{~}{H}_𝐤\mathrm{e}^{i(𝐤.𝐱E_kt)}+(\stackrel{~}{H}_𝐤)^{}\mathrm{e}^{i(𝐤.𝐱E_kt)}\right]$$
(11)
with $`E_k=\sqrt{𝐤^2+M_H^2}`$. Eqs.(10) and (11) are reobtained in the limit $`\frac{\delta }{M_H}\frac{M_H}{\mathrm{\Lambda }}0`$ where the wavelengths associated to $`h(x)`$ become infinitely large in units of the physical scale set by $`\xi _H=1/M_H`$. In this limit, where for any finite value of $`𝐤`$ the broken phase has only massive excitations, one recovers an exactly Lorentz-covariant theory.

Now, as anticipated, the interpretation of the gap-less branch in terms of density fluctuations places no particular problem . In fact, ”Any quantum liquid consisting of particles with integral spin (such as the liquid isotope 4He) must certainly have a spectrum of this type…In a quantum Bose liquid, elementary excitations with small momenta $`𝐤`$ (wavelengths large compared with distances between atoms) correspond to ordinary hydrodynamic sound waves, i.e. are phonons. This means that the energy of such quasi-particles is a linear function of their momentum” . In this sense, a superfluid vacuum provides for $`𝐤0`$ a universal picture. This result does not depend on the details of the short-distance interaction and even on the nature of the elementary constituents. For instance, the same coarse-grained description is found in superfluid fermionic vacua  that, as compared to the Higgs vacuum, bear the same relation of superfluid 3He to superfluid 4He.

On the other hand, a full analogy with 4He would also require to establish the interpretation of the massive branch for $`|𝐤|\delta `$ as a ‘roton’, i.e. in terms of a suitable vortical motion in the superfluid, and this is by no means obvious. For instance, Landau’s roton spectrum

$$E_{\mathrm{roton}}(𝐤)=\mathrm{\Delta }+\frac{𝐤^2}{2\mu }$$
(12)
depends on two parameters $`\mathrm{\Delta }`$ and $`\mu `$ that, in superfluid 4He, are vastly different. In fact, (in units $`\mathrm{}=c=1`$, where $`\mathrm{\Delta }`$ and $`\mu `$ have the same physical dimensions) one finds $`\mathrm{\Delta }710^4`$ eV and $`\mu 610^8`$ eV while $`E_a(𝐤)=\sqrt{𝐤^2+M_H^2}`$ depends on a single mass parameter $`M_H`$. Under which conditions can a superfluid medium exhibit rotons with $`\mathrm{\Delta }=\mu `$ ? Moreover, even if $`\mathrm{\Delta }=\mu `$, does this value agree with the Higgs mass parameter $`M_H`$ obtained in quantum field theory ?

The answer to this type of questions can only be obtained by combining a field-theoretical description of the condensation phenomenon with the basic ingredients of quantum hydrodynamics. This analysis represents the main content of this paper and will be presented in the following. In Sect.2, we shall first review the formalism of ref. with the hierarchical pattern of scales that is established in the scalar condensate. Further, in Sect.3, we shall use the formalism of quantum hydrodynamics and discuss the interpretation of the massive branch as a roton. Finally, Sect.4 will contain a summary together with other possible consequences of our approach.

## 2. The Higgs condensate and a hierarchy of scales

We shall now first resume the main results of ref. in the case of a one-component $`\lambda \mathrm{\Phi }^4`$ theory, a system where the condensing quanta are just neutral spinless particles, the ‘phions’.

One starts by quantizing the scalar field $`\mathrm{\Phi }(𝐱)`$ in terms of $`a_𝐤`$, $`a_𝐤^{}`$, the annihilation and creation operators for the elementary phions whose ‘empty’ vacuum state $`|o`$ is defined through $`a_𝐤|o=o|a_𝐤^{}`$=0.

The phion system is assumed to be contained within a finite box of volume $`𝖵`$ with periodic boundary conditions. There is then a discrete set of allowed modes $`𝐤`$. In the end one takes the infinite-volume limit and the summation over allowed modes goes over to an integration: $`_𝐤𝖵d^3k/(2\pi )^3`$. In this way, the scalar field is expressed as

$$\mathrm{\Phi }(𝐱)=\underset{𝐤}{}\frac{1}{\sqrt{2𝖵E_k}}\left[a_𝐤\mathrm{e}^{i𝐤.𝐱}+a_𝐤^{}\mathrm{e}^{i𝐤.𝐱}\right],$$
(13)
where $`E_k=\sqrt{𝐤^2+m_\mathrm{\Phi }^2}`$, $`m_\mathrm{\Phi }`$ being the physical, renormalized phion mass.

Bose condensation means that in the ground state there is an average number $`𝖭_0`$ of phions in the $`𝐤=0`$ mode, where $`𝖭_0`$ is a finite fraction of the total average number $`𝖭`$

$$𝖭=\underset{𝐤}{}a_𝐤^{}a_𝐤$$
(14)
At zero temperature, if the gas is dilute, almost all the particles are in the condensate; $`𝖭_0(T=0)𝖭`$. In fact, the fraction which is not in the condensate (‘depletion’)

$$D=1\frac{𝖭_0}{𝖭}=𝒪(\sqrt{na^3})$$
(15)
is a phase-space effect that becomes negligible for a very dilute system  where

$$ϵ=\sqrt{na^3}1$$
(16)
In Eqs.(15) and (16) we have introduced the phion density

$$n=\frac{𝖭}{𝖵}$$
(17)
and the phion-phion scattering length Eq.(1).

Therefore, for a very dilute system, where, to a first approximation, one neglects the residual operator part of $`a_{𝐤=0}`$, one gets $`a_{𝐤=0}^{}a_{𝐤=0}𝖭`$ and so, $`a_{𝐤=0}`$ can be identified with the c-number $`\sqrt{𝖭}`$ (up to a phase factor). In this way

$$\varphi =\mathrm{\Phi }=\frac{1}{\sqrt{2𝖵m_\mathrm{\Phi }}}(a_{𝐤=0}^{}+a_{𝐤=0})\sqrt{\frac{2𝖭}{𝖵m_\mathrm{\Phi }}}.$$
(18)
or

$$n(\varphi ^2){\scriptscriptstyle \frac{1}{2}}m_\mathrm{\Phi }\varphi ^2,$$
(19)
With this identification, setting $`a_{𝐤=0}^{}=a_{𝐤=0}=\sqrt{𝖭}`$ is equivalent to shifting the quantum field $`\mathrm{\Phi }`$ by a constant term $`\varphi `$. Further, using Eq.(19), the energy density $`=(n)`$ can be translated into the effective potential

$$V_{\mathrm{eff}}(\varphi )=(n)$$
(20)
Now, by exploring the limit $`m_\mathrm{\Phi }0`$ in the class of ‘triviality compatible’ approximations to the effective potential  one discovers non-trivial absolute minima $`\varphi =\pm v0`$ of $`V_{\mathrm{eff}}(\varphi )=V_{\mathrm{triv}}(\varphi )`$ and Bose condensation with an average density $`\overline{n}=n(v^2)`$.

The basic relations of the broken phase are

$$M_H^2\lambda v^2\overline{n}a$$
(21)
and

$$m_\mathrm{\Phi }^2\lambda ^2v^2ϵ^2\overline{n}a$$
(22)
The key-ingredient to understand why the phase transition occurs for a still positive value of $`m_\mathrm{\Phi }`$, consists in the observation  that the phion-phion interaction is not always repulsive. Besides the contact $`+\lambda \delta ^{(3)}(𝐫)`$ potential there is an induced attraction $`\lambda ^2\frac{e^{2m_\mathrm{\Phi }r}}{r^3}`$ from ultraviolet-finite parts of higher order graphs (see also ref.) that, differently from the usual ultraviolet divergences, cannot be re-absorbed into a standard redefinition of the tree-level coupling. For small values of $`\lambda `$ and for sufficiently small values of $`m_\mathrm{\Phi }`$ the corresponding graphs, when taken into account consistently in the effective potential, can compensate for the effects of both the short range repulsion and of the non-zero phion mass. In this situation, the perturbative empty vacuum state $`|o`$, although locally stable, is not globally stable and the lowest energy state becomes a state with a non-zero density of phions that are Bose condensed in the zero-momentum state.

We emphasize that this weakly first-order scenario of symmetry breaking is discovered in a class of approximations to the effective potential, just those that are consistent with the assumed exact ‘triviality’ property of the theory in 3+1 space-time dimensions . In any case, it can be objectively tested against the standard picture based on a second-order phase transition. To this end one can run numerical simulations near the phase transition region and compare the predictions of refs. with the conventional existing two-loop or renormalization-group-improved forms of the effective potential. When this is done, the quality of the fits to the existing lattice data  favours unambiguosly the first-order scenario of refs..

Let us now consider the range of momenta associated with the two different branches of the energy spectrum. In condensed matter, the transition between acoustic branch and single-particle spectrum corresponds to their matching at a momentum scale set by the inverse mean free path for the elementary constituents. As anticipated in the Introduction, in the scalar condensate, the matching condition corresponds to a momentum scale $`\delta M_H`$ where $`E_a(\delta )E_b(\delta )`$ or

$$c_s\delta \sqrt{\delta ^2+M_H^2}M_H+\frac{\delta ^2}{2M_H}$$
(23)
with $`\delta \frac{1}{r_{\mathrm{mfp}}}`$, $`r_{\mathrm{mfp}}`$ being the phion mean free path for $`n=\overline{n}`$ in Eq.(4).

Now, the scattering length $`a`$ can be used to define a far ultraviolet scale

$$\mathrm{\Lambda }1/a$$
(24)
up to which phions can be treated as ‘hard spheres’. Using Eqs.(3), (4), (21), and (23), this yields

$$t=\frac{\mathrm{\Lambda }}{M_H}\sqrt{\frac{1}{\overline{n}a^3}}$$
(25)
and

$$\frac{1}{c_s}\frac{\delta }{M_H}\sqrt{\overline{n}a^3}$$
(26)
Therefore, in the continuum limit where $`t\mathrm{}`$ one gets an infinitely dilute Higgs condensate where $`ϵ=\sqrt{\overline{n}a^3}0`$ and the hierarchy of scales

$$\delta M_H\mathrm{\Lambda }$$
(27)
related as in Eq.(6).

Finally, by using Eq.(1), the condition for spontaneous symmetry breaking Eq.(22) can also be expressed as

$$m_\mathrm{\Phi }\overline{n}a^2\frac{1}{r_{\mathrm{mfp}}}$$
(28)
so that the phion mean free path in the condensate is of the same order as the phion Compton wavelength $`1/m_\mathrm{\Phi }`$.

Notice that the order of magnitude of the sound velocity

$$c_s\frac{M_H}{\delta }\frac{\mathrm{\Lambda }}{M_H}\frac{1}{ϵ}$$
(29)
is much larger than unity (in units of the light velocity $`c=1`$). Actually, $`c_s`$ must diverge in the continuum limit where Lorentz-covariance becomes exact . In this limit, where the vacuum acquires an infinite rigidity, the condensate becomes incompressible and the massive branch of the energy spectrum remains valid down to $`𝐤=0`$.

The presence of superluminal sound in the scalar condensate has different motivations. First of all we observe that, on a general ground, ”..it is an open question whether $`\frac{c_s}{c}`$ remains less than unity when nonelectromagnetic forces are taken into account” . For this reason, several authors  have considered the possibility of media whose long-wavelength compressional modes for $`𝐤0`$ have phase and group velocity $`\frac{E}{|𝐤|}=\frac{dE}{d|𝐤|}=c_s>c`$.

This possibility depends on the approximate nature of locality in cutoff-dependent quantum field theories where the elementary quanta are treated as ‘hard spheres’. In this case, in fact, a hard-sphere radius is known  to imply a superluminal propagation within the sphere boundary. Now, in the perturbative empty vacuum state (with no condensed quanta) such superluminal propagation is restricted to very short wavelengths, smaller than the inverse ultraviolet cutoff. However, in the condensed vacuum, the hard spheres can ‘touch’ each other so that the actual propagation of density fluctuations in a hard-sphere system might take place at a superluminal speed. This intuitive idea is at the base of the ‘macroscopic’ violations of locality discussed in ref. (and of the ‘reentrant’ violations of special relativity in the low-energy corner  mentioned above).

In some case, superluminal sound is known to arise  when a large negative bare mass and a large positive self-energy combine to produce a very small physical mass, just the situation expected for the quanta of the scalar condensate. In this way, the physical origin of the superluminal sound is traced back to the asymmetric role of mass renormalization: it subtracts out self-interaction energy without altering the tree-level interparticle interactions that contribute to the pressure.

On the other hand, following refs., superluminal sound is also consistent with the equation of state of a perfect fluid whose energy density has a minimum at some given value of the particle density, as it happens in the scalar condensate. Just for this reason, near the minimum, long-wavelength density fluctuations represent nearly instantaneous effects that can propagate at arbitrarily large speed. In this sense, the scalar condensate rensembles an elastic medium near the inconpressibility limit where the Poisson ratio $`\nu 1/2`$ and the propagation speed of the longitudinal waves of dilatation diverges in units of the propagation speed of the transverse waves of distortion .

Summarizing the previous results, we find that in the local limit of the theory, where $`\mathrm{\Lambda }/M_H\mathrm{}`$, one also finds $`\delta /M_H0`$ so that the energy spectrum $`E(𝐤)`$ reduces to $`E_a(𝐤)=\sqrt{𝐤^2+M_H^2}`$ in the whole range of $`|𝐤|`$. In the cutoff theory, however, one should expect infinitesimal deviations in an infinitesimal region of three-momenta. For instance, assuming $`\mathrm{\Lambda }=10^{19}`$ GeV and $`M_H=250`$ GeV, a scale $`\delta \frac{M_H^2}{\mathrm{\Lambda }}10^5`$ eV, for which $`\frac{\delta }{M_H}410^{17}`$, might well represent the physical realization of a formally infinitesimal quantity. If this were the right order of magnitude, the collective density fluctuations of the Higgs vacuum described by $`E_b(𝐤)`$ have wavelengths $`>\frac{2\pi }{\delta }`$ thus ranging from about a centimeter up to infinity .

On the other hand, for $`|𝐤|\delta `$, the excitation spectrum describes single-particle states of mass $`M_H\sqrt{\overline{n}a}`$. In the following section, we shall show that these states can be interpreted as elementary excitations associated with a vortical motion.

## 3. The ‘roton’ picture of the massive branch

In his theory, Landau suggested that in a superfluid medium there must be elementary vortex excitations, the ‘rotons’, whose energy has the form in Eq.(12). In his original papers , Landau did not work out an explicit derivation of Eq.(12). This was, however, deduced subsequently by Ziman  whose formalism we shall briefly resume in the following (for the convenience of the reader we shall adopt in this section the same notations of ref.).

Ziman’s starting point is the form of the Hamiltonian density of a fluid

$$H=\frac{1}{2}\rho 𝐮^2+\rho W(\rho )$$
(30)
where $`\rho `$ is the mass density and $`W`$ the internal energy whose mimimum is obtained for $`\rho \overline{\rho }`$. The fluid velocity

$$𝐮=\phi \frac{i}{2\rho }(\mathrm{\Psi }^{}\mathrm{\Psi }\mathrm{\Psi }\mathrm{\Psi }^{})$$
(31)
is expressed in terms of the three Clebsch potentials , $`\phi `$, $`\mathrm{\Psi }`$ and $`\mathrm{\Psi }^{}`$. Notice that, although the fluid is described in terms of the 4-component field $`(\rho ,\phi ,\mathrm{\Psi },\mathrm{\Psi }^{})`$, there are no fundamental ‘charges’ and all dynamical effects derive from the possible states of motion of the fluid.

In quantum hydrodynamics $`(\rho ,\phi )`$ and $`(\mathrm{\Psi },\mathrm{\Psi }^{})`$ are pairs of canonically conjugated variables, i.e.

$$[\rho (𝐫),\phi (𝐫^{})]=i\delta ^3(𝐫𝐫^{})$$
(32)
and

$$[\mathrm{\Psi }(𝐫),\mathrm{\Psi }^{}(𝐫^{})]=\delta ^3(𝐫𝐫^{})$$
(33)
all other commutators being zero. In this way, one obtains Landau’s relations for the commutators of the velocity components

$$[𝐮_x(𝐫),𝐮_y(𝐫^{})]=\frac{1}{i\rho }\delta ^3(𝐫𝐫^{})\zeta _z$$
(34)
where

$$\zeta =\mathrm{x}𝐮$$
(35)
is the vorticity.

In the incompressibility limit, where $`\rho =\overline{\rho }`$ and the phase $`\phi `$ are constant throughout the volume of the fluid, the fluid Hamiltonian density reduces to its ‘roton’ part

$$H_{\mathrm{roton}}=\frac{1}{8\overline{\rho }}(\mathrm{\Psi }^{}\mathrm{\Psi }\mathrm{\Psi }\mathrm{\Psi }^{})^2$$
(36)
In a cubic box of volume $`𝖵`$, one can expand in plane waves

$$\mathrm{\Psi }=\frac{1}{\sqrt{𝖵}}\underset{𝐤}{}b_𝐤e^{+i𝐤𝐫}\mathrm{\Psi }^{}=\frac{1}{\sqrt{𝖵}}\underset{𝐤}{}b_𝐤^{}e^{i𝐤𝐫}$$
(37)
with

$$[b_𝐤,b_𝐥^{}]=\delta _{\mathrm{𝐤𝐥}}$$
(38)
so that by integrating the Hamiltonian density over the whole volume one gets quadratic $`b_𝐤^{}b_𝐤`$ terms. These can be used to define a free-roton Hamiltonian

$$_{\mathrm{roton}}^{(o)}=\frac{1}{8\overline{\rho }𝖵}\underset{\mathrm{𝐤𝐥}}{}(𝐤^2+𝐥^2)b_𝐤^{}b_𝐤$$
(39)
that, after converting the discrete sum into an integral, gives finally

$$E_{\mathrm{roton}}(𝐤)=\frac{1}{8\overline{\rho }}\frac{d^3𝐪}{(2\pi )^3}(𝐤^2+𝐪^2)$$
(40)
Eq.(40) is formally divergent, so that, to extract the relevant values of $`\mathrm{\Delta }`$ and $`\mu `$ in Eq.(12), one has to find out a suitable cutoff momentum for the single-roton excitations, say $`|𝐪|=q_{\mathrm{max}}`$. If this is done, we obtain

$$\mathrm{\Delta }=\frac{q_{\mathrm{max}}^5}{80\pi ^2\overline{\rho }}$$
(41)
and

$$\mu =\frac{24\pi ^2\overline{\rho }}{q_{\mathrm{max}}^3}$$
(42)
As a consequence of the introduction of $`q_{\mathrm{max}}`$, the uncertainty associated with any cutoff procedure will only allow an order of magnitude estimate of $`\mathrm{\Delta }`$ and $`\mu `$ and of the regime of parameters associated with the Lorentz-covariant condition $`\mathrm{\Delta }=\mu `$ .

To obtain an estimate of $`q_{\mathrm{max}}`$, we observe that using the operators $`b_𝐤^{}`$ and $`b_𝐤`$ one can construct the roton Fock space labelled by $`N_𝐤`$, the eigenvalue of the roton number operator $`b_𝐤^{}b_𝐤`$, the single-roton states corresponding to $`N_𝐤=1`$. Now, by computing the value of the vorticity vector in a single-roton state, one gets the idea of the roton ”… as a steady rotational motion of the fluid, capable of moving as a ‘vortex’ through the liquid” .

This observation arises by inspection of the single-roton wave functions expressed in cylindrical co-ordinates $`(r,\theta ,z)`$. In this case, by introducing the squared transverse momentum $`\kappa ^2=k_x^2+k_y^2`$, the azimuthal quantum number $`\nu `$, and the z-component of the momentum $`k_z`$, one finds the equivalent form for the free-roton Hamiltonian

$$_{\mathrm{roton}}^{(o)}=\frac{1}{8\overline{\rho }𝖵}\underset{\kappa ,\kappa ^{},\nu ,\nu ^{},k_z,k_z^{}}{}(\kappa ^2+k_z^2+(\kappa ^{})^2+(k_z^{})^2)b_{\kappa \nu k_z}^{}b_{\kappa \nu k_z}$$
(43)
whose single-particle eigenvalues

$$E_{\kappa \nu k_z}=\mathrm{\Delta }+\frac{\kappa ^2+k_z^2}{2\mu }$$
(44)
have the same $`\mathrm{\Delta }`$ and $`\mu `$ as in Eqs.(41) and (42) with $`q_{\mathrm{max}}^2=(\kappa ^2+k_z^2)_{\mathrm{max}}`$. This can be understood by first noticing that the quantities $`\mathrm{\Delta }`$ and $`\mu `$ can be expressed in terms of the number of quantum states with $`|𝐪|q_{\mathrm{max}}`$, say $`g(q_{\mathrm{max}})`$, as

$$\mathrm{\Delta }q_{\mathrm{max}}^2\frac{g(q_{\mathrm{max}})}{\overline{\rho }𝖵}$$
(45)
and

$$\mu \frac{\overline{\rho }𝖵}{g(q_{\mathrm{max}})}$$
(46)
Further, using the results of ref., one finds the same leading behaviour

$$\frac{g(q_{\mathrm{max}})}{𝖵}\frac{q_{\mathrm{max}}^3}{6\pi }$$
(47)
by switching from cubical to cylindrical geometry. Finally, averaging over all possible orientations gives $`(k_x^2)_{\mathrm{max}}=(k_y^2)_{\mathrm{max}}=(k_z^2)_{\mathrm{max}}`$ and one obtains $`q_{\mathrm{max}}^2\frac{3}{2}(\kappa ^2)_{\mathrm{max}}`$.

In cylindrical coordinates, the single-particle wave functions can be expressed as

$$\psi =\psi _{\kappa ,\nu ,k_z}(r,\theta ,z)=𝒩J_\nu (\kappa r)e^{i\nu \theta }e^{ik_zz}$$
(48)
where $`J_\nu (\kappa r)`$ are Bessel functions and $`𝒩`$ a normalization factor. In this case, one finds $`u_r=u_z=0`$ with the only non vanishing component of the velocity being

$$u_\theta =\frac{1}{2ir}(\psi ^{}\frac{\psi }{\theta }\psi \frac{\psi ^{}}{\theta })=𝒩^2\frac{\nu J_\nu ^2(\kappa r)}{r}$$
(49)
Analogously, $`\zeta _r=\zeta _\theta =0`$ and the only non-zero component of the vorticity is

$$\zeta _z=\frac{1}{r}\frac{(ru_\theta )}{r}=𝒩^2\frac{\nu }{r}\frac{d}{dr}J_\nu ^2(\kappa r)$$
(50)
By introducing the vortex radius $`R`$ as the value at which $`u_\theta (r=R)=0`$, and using the relation

$$_0^Rr𝑑rJ_\nu ^2(\kappa r)=\frac{R^2}{2}J_{\nu +1}^2(\kappa R)$$
(51)
(with $`J_\nu (\kappa R)=0`$) we can set finally $`𝒩=\sqrt{\frac{2}{R}}\frac{1}{J_{\nu +1}(\kappa R)}`$.

Therefore, a reasonable value of the cutoff momentum for single-roton states is

$$q_{\mathrm{max}}\kappa _{\mathrm{max}}\frac{1}{R_{\mathrm{min}}}$$
(52)
where $`R_{\mathrm{min}}`$ denotes the minimum transverse size of the thinnest vortices that can be established in the superfluid. These are the so called ‘vortex filaments’, whose transverse size can be obtained from ref. in terms of the particle density $`n`$ and of the scattering length $`a`$. In this case, for a dilute (‘almost ideal’) Bose condensate, the filament ‘core’ is

$$R_{\mathrm{min}}r_{\mathrm{core}}\frac{1}{\sqrt{na}}$$
(53)
so that, for $`n=\overline{n}`$, we find

$$q_{\mathrm{max}}\frac{1}{r_{\mathrm{core}}}\sqrt{\overline{n}a}$$
(54)
In this way, assuming for the mass of the fluid constituents the same Eq.(28) established for the quanta of the scalar condensate, one finds

$$\overline{\rho }=m_\mathrm{\Phi }\overline{n}\overline{n}^2a^2$$
(55)
so that replacing relations (54) and (55) into Eqs.(41) and (42), one gets the order of magnitude estimate

$$\mathrm{\Delta }\mu \sqrt{\overline{n}a}$$
(56)
Therefore, by comparing with Eq.(21), we have identified the relativistic regime $`\mathrm{\Delta }=\mu =M_H`$ foreseen in the Introduction. In fact, $`M_H`$ coincides with the energy-gap $`\mathrm{\Delta }`$ for vortex formation in a superfluid medium possessing the same density and the same type of constituents as the scalar condensate.

More generally, it is interesting to compare the possible regimes of a Bose superfluid, made up of particles with mass $`m`$, in terms of the two dimensionless parameters

$$x=\frac{m}{na^2}$$
(57)
and

$$ϵ=\sqrt{na^3}$$
(58)
In terms of $`x`$, using Eq.(54) for arbitrary $`n`$, Eqs.(41) and (42) give

$$\mathrm{\Delta }\frac{\sqrt{na}}{x}$$
(59)
and

$$\mu x\sqrt{na}x^2\mathrm{\Delta }$$
(60)
so that, as anticipated, the Lorentz-covariant condition $`\mathrm{\Delta }\mu `$ corresponds to $`x1`$.

Notice that the momentum $`\delta 1/r_{\mathrm{mfp}}na^2`$, related to the transition from the phonon branch to the roton-like excitations of Eq.(12) can be used to obtain the sound velocity from the relation

$$c_s\delta \mathrm{\Delta }$$
(61)
so that (in units of $`c`$)

$$c_s\frac{1}{xϵ}$$
(62)
Using Eq.(60), the non-relativistic limit, where $`\mathrm{\Delta }\mu `$, is seen to correspond to very large values of $`x`$ such that also $`1xϵ`$. In this case, the sound velocity becomes

$$c_s\frac{1}{xϵ}\frac{\sqrt{na}}{m}1$$
(63)
which is the Lee-Yang-Huang result for a dilute hard-sphere Bose gas . On the other hand, a Lorentz-covariant form of the massive branch requires a value $`x1`$ in Eq.(60). This, when replaced in Eq.(62), produces the anticipated vastly superluminal sound velocity Eq.(29)

$$c_s1/ϵ$$
(64)
that diverges in the limit where $`ϵ0`$ and the scale $`\mathrm{\Delta }\mu \sqrt{na}=`$ is kept fixed.

The above relations can also be compared with superfluid 4He. Although this is not a dilute system (with typical values $`a2.710^8`$ cm and $`n10^{23}`$ cm-3, one gets $`ϵ1.4`$), we find, nevertheless, a good agreement with our picture. In fact, using the experimental values  $`(\mu )_{\mathrm{exp}}=0.16m_{\mathrm{He}}`$ and $`(\mathrm{\Delta })_{\mathrm{exp}}=7.410^4`$ eV in Eq.(60) one can obtain an experimental value of $`x`$, say

$$x_{\mathrm{exp}}\sqrt{\frac{(\mu )_{\mathrm{exp}}c^2}{(\mathrm{\Delta })_{\mathrm{exp}}}}10^6$$
(65)
that, if used in Eq.(62), produces a value $`c_s0.710^6`$ in good agreement with the experimental result for the sound velocity $`(c_s)_{\mathrm{exp}}=239`$ m/sec . Finally, the theoretical input prediction from Eq.(57)

$$x_{\mathrm{th}}\frac{m_{\mathrm{He}}}{na^2}310^6$$
(66)
is also in fairly good agreement with the experimental result Eq.(65), thus confirming the overall consistency of our picture.

## 4. Summary and outlook

Taking into account the two-valued nature  of the zero-4-momentum connected propagator in the broken phase, one gets the idea of a true physical medium that contains, for $`𝐤0`$, two types of excitations: a massive one, whose energy $`E_a(𝐤)M_H`$ and that corresponds to the usual Higgs boson field, and a gap-less one whose energy $`E_b(𝐤)0`$. The overall picture is similar to the coexistence of phonons and rotons in superfluid 4He that, in fact, is usually considered the condensed-matter analogue of the Higgs condensate.

Now, the gap-less branch is naturally interpreted in terms of the collective density fluctuations of the system . These dominate the physical spectrum for $`𝐤0`$ and their wavelengths are larger than $`r_{\mathrm{mfp}}`$, the mean free path for the condensed quanta.

On the other hand, the continuum limit of quantum field theory corresponds to the ideal case of an incompressible fluid so that the massive energy spectrum $`E_a(𝐤)`$ extends down to $`𝐤=0`$.

In this paper, following the original Ziman’s  approach and using the formalism of ref., we have proposed the interpretation of the massive branch as a roton, i.e. as an elementary excitation that, differently from the phonons associated with the irrotational motions of the superfluid, arises in connection with a non-zero (‘bulk’) vorticity. This interpretation requires the energy spectrum Eq.(12) to exhibit values of $`\mathrm{\Delta }`$ and $`\mu `$ such that

$$\mathrm{\Delta }\mu M_H\sqrt{\overline{n}a}$$
(67)
In turn, this relation depends on the peculiar condition Eq.(28)

$$m_\mathrm{\Phi }\overline{n}a^2$$
(68)
between the mass $`m_\mathrm{\Phi }`$ of the elementary condensing phions, their equilibrium number density $`\overline{n}`$ and their scattering length $`a`$. Relation (68) implies that the phion mean free path $`r_{\mathrm{mfp}}1/(\overline{n}a^2)`$ is of the same order as the phion Compton wavelength $`1/m_\mathrm{\Phi }`$ and is naturally found in the weakly first-order scenario of ref. where the broken phase is represented as a dilute Bose condensate for which $`ϵ=\sqrt{\overline{n}a^3}1`$, the continuum limit corresponding to $`ϵ0`$ with $`\overline{n}a=`$ fixed.

Of course, being used to consider Lorentz covariance an exact built-in requirement, an energy spectrum as in Eq.(12) with $`\mathrm{\Delta }=\mu `$ may seem more or less trivial. For instance, starting from the original annihilation and creation operators for the elementary quanta of the unphysical empty vacuum state $`|o`$, it is easily discovered within the standard covariant generalization of the Bogolubov method . However, our analysis shows that there is an alternative viewpoint. In fact, within quantum hydrodynamics the same result is far from being trivial and is only recovered in the special case (68). In this modified perspective, a Lorentz-covariant massive spectrum corresponds (with some approximations and in a certain range of wavelengths) to vortex formation in a superfluid medium possessing a well defined pattern of scales.

We conclude by observing that the proposed identification of the massive Higgs boson as a roton and of $`M_H`$ with the energy-gap for vortex formation in the superfluid vacuum is not a mere exercise. In fact, in a quantum-hydrodynamical description of the scalar condensate based on the hierarchical structure of scales Eq.(2), the single-roton states have a natural cutoff momentum at

$$q_{\mathrm{max}}\sqrt{\overline{n}a}M_H\mathrm{\Lambda }$$
(69)
(with $`\mathrm{\Lambda }1/a`$).

Therefore, accepting our interpretation, the Higgs boson momentum would be physically cut off at values that are infinitely smaller than $`\mathrm{\Lambda }`$ (by a factor $`\sqrt{\overline{n}a^3}0`$) when approaching the continuum limit. This implies, for instance, that any self-energy part that grows quadratically with the phase space, say

$$\mathrm{\Pi }_H=\underset{n}{}c_n\lambda ^nq_{\mathrm{max}}^2,$$
(70)
is much smaller than $`M_H^2`$, for a weak self-coupling $`\lambda `$. Thus $`M_H`$ emerges as a naturally intermediate, fixed energy scale associated with the condensed vacuum and, in this sense, by treating the Higgs condensate as a real physical medium, one can find a solution of the so called naturalness problem (i.e. why $`M_H`$ is so much smaller than $`\mathrm{\Lambda }`$) without any artificial fine-tuning of the basic parameters. Such fine-tuning problems, instead, appear in the standard approach where the massive Higgs boson is regarded as an ordinary elementary particle propagating in the vacuum and its maximum momentum is identified with the ultimate ultraviolet cutoff $`\mathrm{\Lambda }`$ of the theory.


format: text

text:
# Simulation of Ultra-High Energy Photon Propagation in the Geomagnetic Field

## 1 Program Summary

CPC Program Library Index: 1.1 Cosmic Rays
Title of program:    PRESHOWER 1.0     Catalog number:      Program obtainable:    from P. Homola, e-mail: Piotr.Homola@ifj.edu.pl     Computer on which the program has been thoroughly tested:    Intel-Pentium based PC     Operating system:    Linux, DEC-Unix     Programming languages used:    C, FORTRAN 77     Memory required to execute:    $`<100`$ kB     No. of bits in a word:    32     Has the code been vectorized?    no     Number of lines in distributed program:    1900     Other procedures used in PRESHOWER:    IGRF , bessik, ran2      Nature of physical problem:    simulation of a cascade of particles initiated by UHE photon passing through the geomagnetic field above the Earth’s atmosphere.         Method of solution:    The primary photon is tracked until its conversion into $`e^+e^{}`$ pair or until it reaches the upper atmosphere. If conversion occured each individual particle in the resultant preshower is checked for either bremsstrahlung radiation (electrons) or secondary gamma conversion (photons). The procedure ends at the top of atmosphere and the shower particle data are saved.     Restrictions on the complexity of the problem:    gamma conversion into particles other than electron pair hasn’t been taken into account     Typical running time:    100 preshower events with primary energy $`10^{20}`$ eV require a 800 MHz CPU time of about 50 min., with $`10^{21}`$ eV the simulation time for 100 events grows up to 500 min.

## 2 Introduction

The existence of ultra-high energy cosmic rays (UHECR) beyond the Greisen-Zatsepin-Kuzmin (GZK) cutoff energy  is one of the main puzzles in today’s astrophysics. The composition of UHECR, if known, could provide a stringent test of models of particle acceleration and propagation towards the Earth. Most theoretical efforts are made towards explaining the acceleration of protons to ultra-high energies. This is because protons seem to be the most promising candidates for the UHECR particles. There are, however, also so-called exotic models which postulate decays of supermassive “X-particles”. A common prediction of these models is that there should be a significant photon flux in the ultra-high energy region. Accurate identification of photon primaries, measurement of the ultra-high energy photon flux, or specifying the upper limit for it, will provide an excellent test for the models of cosmic-ray origin.

Ultra-high energy photons in the presence of the geomagnetic field can convert into an electron-positron pair before they enter the atmosphere. Synchrotron radiation of this pair gives rise to a cascade of photons. Some of them can have energies sufficient for the next conversion. The cascade at the top of the atmosphere thus consists mainly of photons and two or more electrons. We will call this cascade a “preshower” since it originates and develops above the upper atmosphere, i.e. before the “ordinary” shower development in the air.

Particles of a preshower induce subshowers after entering the atmosphere. A superposition of these subshowers will be seen by detectors as one extensive air shower, whose profile is expected to differ from the profile of a shower induced by a single unconverted photon.

Forthcoming experiments with large collecting apertures (Pierre Auger Observatory , EUSO ), that are designed to detect cosmic rays of ultra-high energies, give a good opportunity to test these models. A photon-induced shower might have observable features that allow to distinguish it from showers initiated by hadrons.

In this paper we describe the PRESHOWER program that allows extensive studies of properties of photon-induced showers. The program performs a simulation of the propagation of photons before they enter the Earth’s atmosphere and can be used both as an independent analysis tool as well as a part of an air shower simulation package (e.g. CORSIKA ).

With the PRESHOWER program working in stand-alone mode we performed, as an illustration, an analysis of basic properties of preshowers for the magnetic conditions of Malargüe in Argentina – the site of the Southern Pierre Auger Observatory . The results of these studies, which are presented in Sec. 5, are in general agreement with other analyses .

By linking our program to an air shower simulation code one is able to study showers initiated by various primaries. We checked this by combining our PRESHOWER code with the CORSIKA air shower simulation package which simulates air showers initiated by the preshower particles. We simulated photon-induced showers of different energies and different directions at the Southern Pierre Auger Observatory. Although the interaction in the geomagnetic field changes significantly the shower properties so that primary photon events look more similar to hadronic showers, a discrimination between primary photons and hadrons still seems possible.

The plan of the paper is the following: The relevant effects and the formulas applied in the simulation are summarized in Section 3. The structure of the program is discussed in Section 4. As an illustration, in Section 5 results obtained with the code for the case of the Auger Observatory are discussed. The paper is concluded in Section 6.

## 3 Preshower physics

At energies above 10 EeV, in the presence of the geomagnetic field, a photon can convert into an $`e^+e^{}`$ pair before entering the atmosphere. The resultant electrons will subsequently lose their energy by magnetic bremsstrahlung. The emitted photons can convert again if their energy is high enough. In this way, instead of one high energy photon at the top of atmosphere, a number of less energetic particles, mainly photons and a few electrons, will emerge constituting a preshower. A superposition of subshowers induced by these particles should be seen by detectors as one extensive air shower which reaches its maximum earlier than a shower induced by a single photon of the same primary energy. This is also due to the Landau-Pomeranchuk-Migdal (LPM) effect  which is stronger for photons of higher energies. Both photon conversion and magnetic bremsstrahlung are strongly dependent on the transverse component of the geomagnetic field. In the following, we concentrate on these effects, give the main formulas and discuss their application.

### 3.1 Geomagnetic field model

Since the geomagnetic field is responsible both for gamma conversion and for bremsstrahlung of electrons, it is very important to model it with adequate accuracy. We use the International Geomagnetic Reference Field model (IGRF)  for $`n=10`$, $`n`$ being the highest order of spherical harmonics in the scalar potential. To calculate the field according to the IGRF model we use the numerical procedures written by Tsyganenko . When choosing $`n=1`$, the output of these procedures is equivalent to a dipole model with $`|𝐌|=810^{25}`$ G$``$cm3 and with poles at 79.3N, 71.5W and 79.3S, 108.5.

### 3.2 Magnetic pair production: $`\gamma e^+e^{}`$

The number of pairs created by a high-energy photon in the presence of a magnetic field per path length $`dr`$ can be expressed in terms of the attenuation coefficient $`\alpha (\chi )`$ :

$$n_{pairs}=n_{photons}\{1\mathrm{exp}[\alpha (\chi )dr]\}$$
(1)
where

$$\alpha (\chi )=0.5(\alpha _{em}m_ec/\mathrm{})(B_{}/B_{cr})T(\chi )$$
(2)
with $`\alpha _{em}`$ being the fine structure constant, $`\chi 0.5(h\nu /m_ec^2)(B_{}/B_{cr})`$, $`B_{}`$ is the magnetic field component transverse to the direction of the photon’s motion, $`B_{cr}m_e^2c^3/e\mathrm{}=4.414\times 10^{13}`$ G and $`T(\chi )`$ is the magnetic pair production function. $`T(\chi )`$ can be well approximated by:

$$T(\chi )0.16\chi ^1K_{}^{2}{}_{1/3}{}^{}(\frac{2}{3\chi }),$$
(3)
where $`K_{1/3}`$ is the modified Bessel function of order $`1/3`$. For small or large arguments $`T(\chi )`$ can be approximated by

$$\begin{array}{c}T(\chi )\{\begin{array}{cc}0.46\mathrm{exp}(\frac{4}{3\chi }),\hfill & \chi 1;\hfill \\ 0.60\chi ^{1/3},\hfill & \chi 1.\hfill \end{array}\end{array}$$
(4)
Representative values of $`T(\chi )`$ are listed in Table 1.

We use Eq. (1) to calculate the probability of $`\gamma `$ conversion over a small path length $`dr`$:

$$p_{conv}(r)=1\mathrm{exp}[\alpha (\chi (r))dr]\alpha (\chi (r))dr$$
(5)
and the probability of $`\gamma `$ conversion over a longer distance $`R`$:

$$P_{conv}(R)=1\mathrm{exp}[_0^R\alpha (\chi (r))𝑑r]$$
(6)
### 3.3 Magnetic bremsstrahlung

After photon conversion, the electron-positron pair propagates. The energy distribution in an $`e^+e^{}`$ pair is computed according to Ref. :

$$\frac{d\alpha (\epsilon ,\chi )}{d\epsilon }\frac{\alpha _{em}m_ecB_{}}{\mathrm{}B_{cr}}\frac{3^{1/2}}{9\pi \chi }\frac{[2+\epsilon (1\epsilon )]}{\epsilon (1\epsilon )}K_{2/3}\left[\frac{1}{3\chi \epsilon (1\epsilon )}\right],$$
(7)
where $`\epsilon `$ denotes the fractional energy of an electron and the other symbols were explained in the previous chapter. The probability of asymmetric energy partition grows with the primary photon energy and with the magnetic field. Beginning from $`\chi >10`$, the asymmetric energy partition is even more favored than the symmetric one.

Electrons traveling at relativistic speeds in the presence of a magnetic field emit bremsstrahlung radiation (synchrotron radiation). For electron energies $`Em_ec^2`$ and for $`B_{}B_{cr}`$, the spectral distribution of radiated energy is given in Ref. :

$$f(y)=\frac{9\sqrt{3}}{8\pi }\frac{y}{(1+\xi y)^3}\left\{_y^{\mathrm{}}K_{5/3}(z)𝑑z+\frac{(\xi y)^2}{1+\xi y}K_{2/3}(y)\right\}$$
(8)
where $`\xi =(3/2)(B_{}/B_{cr})(E/m_ec^2)`$, $`E`$ and $`m_e`$ are electron initial energy and rest mass respectively, $`K_{5/3}`$ and $`K_{2/3}`$ are modified Bessel functions, and $`y`$ is related to the emitted photon energy $`h\nu `$ by

$$y(h\nu )=\frac{h\nu }{\xi (Eh\nu )};\frac{dy}{d(h\nu )}=\frac{E}{\xi (Eh\nu )^2}$$
(9)
The total energy emitted per unit distance is (in CGS units)

$$W=\frac{2}{3}r_0^2B_{}^2\left(\frac{E}{m_ec^2}\right)^2_0^{\mathrm{}}f(y)𝑑y$$
(10)
with $`r_0`$ being the classical electron radius. For our purposes we use the spectral distribution of radiated energy defined as

$$I(B_{},E,h\nu )\frac{h\nu dN}{d(h\nu )dx},$$
(11)
where $`dN`$ is the number of photons with energy between $`h\nu `$ and $`h\nu +d(h\nu )`$ emitted over a distance $`dx`$. From Eqs. (8), (9), (10), and (11) we get 111 Expression (12), valid for all values of $`h\nu `$, is equivalent to Eq. (2.5a) in Ref. . A simplified form of distribution (12) is given by Eq. (2.10) in Ref. , however it can be used only for $`h\nu E`$.

$$I(B_{},E,h\nu )=\frac{2}{3}r_0^2B_{}^2\left(\frac{E}{m_ec^2}\right)^2f(h\nu )\frac{E}{\xi (Eh\nu )^2}.$$
(12)
Provided $`dx`$ is small enough, $`dN`$ can be interpreted as a probability of emitting a photon of energy between $`h\nu `$ and $`h\nu +d(h\nu )`$ by an electron of energy $`E`$ over a distance $`dx`$. In our simulations we use a small step size of $`dx=1`$ km. The total probability of emitting a photon in step $`dx`$ can then be written as

$$P_{brem}(B_{},E,h\nu ,dx)=𝑑N=dx_0^EI(B_{},E,h\nu )\frac{d(h\nu )}{h\nu }.$$
(13)
The shape of the bremsstrahlung spectrum for two different energies is shown in Fig. 1. The three curves in each plot correspond to different values of the transverse magnetic field: 0.3 G, 0.1 G and 0.03 G from the uppermost to the lowest curve respectively. The energy of the emitted photon is simulated according to the probability density distribution $`dN/d(h\nu )`$ obtained from Eq. 13.

### 3.4 Related physics topics

In this section we shortly discuss the importance of other physics effects connected to our main subject. These are:

* The geomagnetic field deflects the electrons from the initial preshower direction. To evaluate this deflection we use the equation

$$\frac{m_{rel}v^2}{R}=evB_{};m_{rel}=\frac{E}{c^2}$$
(14)
where $`v`$ is the speed of electron and $`e`$ is the elementary charge, to find the radius of curvature $`R`$. If $`R`$ is much larger than the electron path length $`L`$ and assuming constant values of $`B_{}`$, $`E`$ and $`vc`$, the linear displacement from the preshower core in the plane perpendicular to the initial direction is well approximated by:

$$\mathrm{\Delta }x\frac{L^2}{2R}=\frac{ecB_{}L^2}{E}$$
(15)
For typical values like $`B_{}=0.1`$ G, $`L=1000`$ km (see Section 5) and relatively small energy of $`E=10^{18}`$ eV, we get $`R10^{11}`$ km which gives $`\mathrm{\Delta }x0.5`$ cm. Assuming a larger value for $`L`$, one should keep in mind that at higher altitudes $`𝐁`$ is lower and $`E`$ is higher. Allowing for other possible combinations of $`B_{}`$, $`E`$ and $`L`$ one comes to a very safe conclusion that in any case $`\mathrm{\Delta }x<`$ 1 m.
* The initial angular spread of particles in both pair production and bremsstrahlung can be approximated by $`\mathrm{\Delta }\theta m_e/E`$ . Even for energies as low as $`10^{18}`$ eV one has $`\mathrm{\Delta }\theta 10^{12}`$ which even for extremely long paths like for instance 10000 km results in negligible linear deviations at the top of the atmosphere of $`\mathrm{\Delta }x10^3cm`$.
* A $`\gamma `$-ray traveling through the vicinity of the Sun can cascade in the magnetosphere of the Sun where the field is much stronger than near the Earth. The particles originated in conversion near the Sun can be significantly deflected from their initial direction by the solar magnetic field. The cascade of such particles after reaching the Earth’s atmosphere will induce an air shower of large lateral extent which might be well distinguishable from proton or iron showers. However, the predicted rate of such events for primary energies of order $`10^{19}`$ eV is very low: about 1 event per 10 years for the Pierre Auger Observatory . In our simulations we do not consider cascading in the Sun’s magnetosphere.
* The magnetic field embedded in the solar wind is of order $`10^5`$ G  and it influences the Earth’s magnetosphere at altitudes of order 8-10$`R_{}`$. Even for a primary energy as high as $`10^{21}`$ eV gamma conversion takes place within 2$`R_{}`$ above the sea level (see Section 5) where the Earth magnetic field is not less than $`10^2`$ G. Thus, the solar wind can be neglected in our studies.
* The preshower electrons traveling at speeds $`v<c`$ are delayed with respect to the photons. This delay however is very small and can be neglected. For an extreme case of a $`10^{18}`$ eV electron traveling a path of 20000 km, the delay is of order $`10^{26}`$ s. The delay of the particles not parallel to the main direction is also negligible (see the above discussion on the lateral spread of a preshower).

## 4 Structure of the program

The following subsections refer to the stand-alone mode of the PRESHOWER program. The PRESHOWER-CORSIKA mode requires only the main procedure generating events, preshw, which is described in subsection 4.3. All the input data in the PRESHOWER-CORSIKA mode are provided via the main CORSIKA input file, the standard output is printed within the CORSIKA output, and the multirun.dat file is not returned, since there is only one run per one primary trajectory allowed in CORSIKA. There is also a different random number generator used in the PRESHOWER-CORSIKA mode than in the stand-alone one. For more information on the PRESHOWER-CORSIKA mode the user can refer to the latest CORSIKA manual  or contact the corresponding author.

### 4.1 The source code and other files in the PRESHOWER packet

The PRESHOWER code consists of the following files:

contains the main procedure generating preshowers, preshw, and the other procedures in C that are called by preshw.

is used only in the stand-alone mode. It reads the input file INPUT and calls the preshower procedure preshw.

is an external procedure for calculation of the geomagnetic field according to the IGRF model (author: N.A.Tsyganenko ), it is called by preshw.

is a dummy Fortran routine used only in the stand-alone version for the proper linking of the program. In the PRESHOWER-CORSIKA mode the actual rmmard is responsible for random number generation.

For the user’s convenience there are some more files included in the PRESHOWER packet:

compilation and linking commands under linux.

short information on the content of the PRESHOWER packet and on the installation process.

typical input data read by the program working in the stand-alone mode. In the PRESHOWER-CORSIKA mode all the input data required by preshw are read directly from the main CORSIKA input file.

the example of a detailed output file containing information on the preshower particle data being the result of running the program with the input INPUT.

the example of an output file summarizing the program run with the input INPUT; this file is produced only in the stand-alone mode.

the standard output produced by the program in the stand-alone mode with the use of the input data from INPUT.

### 4.2 Input and output data

The input data are read by the main program from the file INPUT and then passed to the preshw procedure. These data are

1. the primary photon energy, \[GeV\]
2. the zenith angle of the shower trajectory in the local reference frame (see Sec. 5), \[deg\]
3. the azimuth angle of the shower trajectory in the local frame, \[deg\]
4. the top of the atmosphere – e.g. the value obtained from CORSIKA, \[km a.s.l.\]
5. geographical position of the observatory: longitude \[deg\]
6. geographical position of the observatory: latitude \[deg\]
7. the year of observation (needed for calculation of the field)
8. number of runs

An exemplary set of the input data can be found in the file INPUT.

There are two output files produced by the PRESHOWER program: part\_out.dat and multirun.dat. Both files begin with a header containing all the crucial run information, then comes the description of columns and then the data. The data in part\_out.dat are presented in two columns: particle id (1 - photons, 2 - positrons, 3 - electrons) and particle energy at the top of the atmosphere (in eV). The columns in multirun.dat are explained as follows:

1. the run number
2. altitude of the first photon conversion in km a.s.l.
3. total number of particles at the top of the atmosphere
4. number of photons
5. number of electrons
6. maximum electron energy
7. maximum photon energy
8. fraction of energy carried by electrons
9. total energy carried by preshower particles = primary photon energy
10. total energy carried by electrons
11. total energy carried by photons

The information printed as standard output begins with the same header as the output files and is followed by the specific information on each run (all energies given in eV): the run number, the altitude of each conversion together with the energy partition information, preshower summary at the top of the atmosphere: total energy carried by photons (Etot\_g), total energy carried by electrons (Etot\_e), maximum photon energy (Eg\_max), maximum electron energy (Ee\_max), number of particles (n\_part), number of photons (n\_phot), number of positrons (n\_e+) and number of electrons (n\_e-). The example of a standard output is the file out.txt included in the PRESHOWER packet.

### 4.3 Generation of events

After reading the input parameters (see Sec.4.2) the event generation procedure preshw is called, within which all the calculations are performed. In the initial stage all the necessary parameter transformations are applied and the preshower trajectory vector is computed. The propagation path length is set to be $`5R_{}`$ long, which e.g. for vertical showers means the simulation starting altitude of $`5R_{}`$ above sea level. All the data concerning the particles in the preshower are stored in the array part\_out which at the beginning has only one entry: the primary photon. In steps of $`10`$ km the transverse magnetic field is computed by the procedure B\_transverse, then the probability of conversion is found using Eq. (5). The step size has been checked to be sufficiently small. If the conversion takes place, the particle array is updated by replacing the photon with an $`e^+`$ and adding an $`e^{}`$ at the array end. The energy partition fraction is found by the function ppfrac. To assure accurate bremsstrahlung simulation, beginning from the first gamma conversion the step of 1 km is used for all particle types. In every step, for each photon in the array, we check again for the conversion effect. In case of electrons the bremsstrahlung radiation is simulated using function brem. The probability of emitting a photon is calculated from Eq. (13) with the integration starting from $`h\nu =10^{12}`$ eV and logarithmic steps of 0.01. Bremsstrahlung photons of energies lower than $`10^{12}`$ eV have a very small influence on the air shower evolution (see also Section 5) and hence they are neglected.

Once a bremsstrahlung photon is emitted, its energy is determined by the probability distribution $`dN/d(h\nu )`$ obtained again from Eq. (13), and the energy of the radiating electron is diminished correspondingly. Simulations are finished when the (adjustable) altitude of the top of the atmosphere is reached and all preshower particle parameters are written out (in the stand-alone mode) or passed to CORSIKA.

### 4.4 Compilation, linking and running the program

In order to compile, link and run PRESHOWER under Linux one should first create a file nr\_fun.c in the PRESHOWER home directory and place the header:

#include \<stdio.h\>
#include \<math.h\>

at the beginning of this file. Then, below the header, the source code of the Numerical Recipes  procedures ran2 and bessik (together with the procedures called by bessik: beschb, chebev and nrerror) should be pasted. To compile and link the program one just types “make” in the PRESHOWER home directory. For Unix compilation and linking one needs to replace the proper commands in Makefile. Once the data in the INPUT file are appropriate, the program can be started by typing “./preshower”. The compilation and linking process can be verified by running the program with the exemplary INPUT file included in the packet. If the program runs correctly, it should return the same output files as those included in the packet (part\_out.dat, multirun.dat) and the standard output as in out.txt.

## 5 Illustration of preshower simulations

The results presented in the following are obtained for the magnetic conditions of the Southern Auger Observatory in Malargüe, Argentina (35.2S, 69.2W). Other geographical positions could easily be adopted. The shower trajectories are given in the local frame where the azimuth increases in the counter-clockwise direction and $`\varphi =0^{}`$ means a shower coming from the geographical North.

To visualize how the transverse magnetic field ($`B_{}`$) can change along the shower trajectory, we have chosen some representative directions along which a particle encounters $`B_{}`$ of significantly different strengths: strong, medium and weak. The “strong field direction” is the one along which $`B_{}`$ is maximal compared to other directions with zenith angles $`\theta 60^{}`$ and is defined by the angles $`\theta =60^{}`$ and $`\varphi =177^{}`$. Similarly, the weak $`B_{}`$ trajectory (“weak field direction”) was chosen to resemble the local field direction at Malargüe of $`\theta =50^{}`$, $`\varphi =357^{}`$ . The “medium field direction” is a vertical trajectory ($`\theta =0^{}`$). In the following we will refer to these directions.

### 5.1 Geomagnetic cascading

In Fig. 2 we show how $`B_{}`$ changes with altitude for the strong, medium and weak field directions. The trajectories end at the chosen top of the atmosphere. In addition to the altitude dependence of $`B_{}`$, large directional differences are visible. Of course we expect the strong field direction to coincide with the direction of a stronger preshower effect.

In Fig. 3, the altitudes of first photon conversion are shown for different energies in the strong field direction. The percentage values shown for each curve denote the total conversion probability up to the top of atmosphere, which is adjusted here to the starting altitude of CORSIKA at 112 km. It can be seen that for the magnetic field above Malargüe, the pair production effect becomes important for energies above $`5\times 10^{19}`$ eV. For all energies of interest, even in case of photons at $`10^{21}`$ eV, the conversion altitude is not higher than 13000 km, i.e. about $`2R_{}`$ a.s.l. Thus, the range of interest of the magnetic field strength for the case of Malargüe is $`0.3`$ G $`÷`$ $`0.01`$ G (see Fig. 2).

The distribution plotted in Fig. 3 is consistent with the results of Ref. . The integrated conversion probability computed with Eq. (6) for different $`\gamma `$ energies and for the two trajectories of strong and weak field directions is plotted in Fig. 4. For photons of energies between $`5\times 10^{19}`$ eV and $`5\times 10^{20}`$ eV a change in the arrival direction, i.e. a change in the magnetic field component perpendicular to the photon trajectory, may imply a dramatic increase or decrease of the conversion probability. Thus, for these energies we expect a strong directional dependence of the preshower effect which might serve as characteristic fingerprint of photons as primary cosmic rays.

The directional dependence of the gamma conversion probability is shown in more detail in Fig. 5 for different arrival directions and primary energies. The conversion effect starts to be important for directions towards the magnetic South pole at photon energies around 50 EeV. For energies of about 150 EeV the conversion is very probable for almost every direction. The photons coming from azimuthal angles around 150 at large zenith angles convert more likely than the others. To explain this, one notes that for a given altitude the magnetic field is strongest around the geomagnetic axis and of course the field becomes stronger with decreasing altitude. Since the photons coming from the geomagnetic South, which corresponds to the azimuth about 150 in case of Malargüe, cross the very vicinity of the geomagnetic axis, they must encounter the stronger magnetic field than the others and hence their conversion probability is higher. The larger zenith angle for these photons means that the strong field region is crossed at the lower altitudes where the field is stronger and the conversion process is more effective. The plots in Fig. 5 are consistent with the values given in Ref. .

All the preshowers discussed here are presented in their final stage of development at the “top of the atmosphere” which in our simulations is assumed to be 112 km a.s.l. – where in the CORSIKA code shower simulations are started.

In Fig. 6 mean energy spectra are plotted for 1000 preshowers initiated by photons of energy $`10^{20}`$ eV arriving at the top of atmosphere from the strong field direction for Malargüe. In the top row we show the distribution of photons (left) and the distribution of photons weighted by their energies (right). We start the simulations of the bremsstrahlung radiation for the energy threshold of $`10^{12}`$ eV so the photons of lower energies are absent in the distribution. In this way we take into account the majority of the photons even if those of energies lower than $`10^{16}`$ eV may have little impact on the longitudinal profile of the resultant air shower. The total energy fraction carried by such photons is very low (see the top right plot in Fig. 6) and the induced subshowers fade very quickly. In the simulation run presented in Fig. 6, only 78 out of 1000 primary photons did not convert into $`e^+e^{}`$, these are seen in the top right plot as the narrow peak at $`E=10^{20}`$ eV.

The distributions of the electrons in the same sample of 1000 preshowers are shown in the bottom row of Fig. 6. Again, the spectrum of particles is plotted to the left and the spectrum weighted by energy to the right. Only in three cases out of 1000 a bremsstrahlung photon converted into $`e^+e^{}`$ so the total number of electrons is 1850 and not 1844. When the conversion happens at a lower altitude, the electrons radiate less bremsstrahlung before entering the atmosphere. These cases are seen in the spectrum as a tail towards high energies.

Some characteristics of the preshowers shown in Figure 6 are summarized in Table 2 together with some other exemplary preshower cases. The conversion probabilities are consistent with those in Figure 5. For higher initial energies gamma conversion takes place at higher altitudes (see Fig. 3) and the resultant electrons radiate more bremsstrahlung photons so the number of particles is larger. Additionally for energies as high as $`10^{21}`$ eV there is more than one gamma conversion per each preshower which also enlarges $`N_{part}`$. Comparing strong $`B_{}`$ and weak $`B_{}`$ directions for a given energy, the large directional dependence of preshower evolution can be noted. In case of $`E_0=10^{20}`$ eV from the weak $`B_{}`$ direction we expect unconverted photons initiating cascades with a large depth of shower maximum $`X_{max}`$ due to the LPM effect. The strong $`B_{}`$ forces the initial energy to be split into about 500 particles of which the highest energy particle typically is a photon with the energy about 5 times lower than the initial one. This means that the longitudinal profile of the resultant air shower will have the maximum higher in the atmosphere, which will make it more similar to a proton-like profile. The energy splitting in this case will compete with the LPM effect. This splitting effect for $`E_0=10^{21}`$ eV for the strong $`B_{}`$ direction is even more dramatic: the highest energy particle in the preshower is more than 16 times less energetic than the primary photon. The maximum $`\gamma `$ energy in the preshower for the weak field direction is higher than for the strong field since the conversion probability is lower in a weaker field (see Fig. 4). Detailed distributions of \<$`N_{part}`$\>, \<$`N_{e^+e^{}}`$\>, \<max $`E_\gamma `$\> and \<max $`E_e`$\> for $`E_0=10^{21}`$ eV and the strong $`B_{}`$ are shown in the top left, top right, bottom left and bottom right plots of Figure 7, respectively.

As mentioned above, the number of particles in the preshower is connected to the altitude of the first gamma conversion. This correlation is presented in Figure 8. The higher altitude of conversion gives more particles at the top of the atmosphere because the electron-positron pair radiates bremsstrahlung over a longer path. The increase in particle number is slightly larger at low altitudes because of the stronger magnetic field.

For the same reasons the total energy fraction carried by the electrons should be smaller for higher conversion altitudes. This dependence is shown in Figure 9. One notes that for high altitudes of the first conversion only about 1% of the initial energy of $`10^{20}`$ eV remains in the electrons at the top of the atmosphere, while all the rest is radiated into photons. With decreasing altitudes the energy radiated into photons decreases. The three points in excess of the general trend in Figure 9 are the cases where the primary photon converted at high altitude and one of the resulting bremsstrahlung photons converted close to the top of the atmosphere - at 228, 442 and 133 km a.s.l respectively for the points from right to left. The threshold for $`\gamma `$ conversion in this case is about $`510^{19}`$ eV, so at least one member of $`e^+e^{}`$ pair had an initial energy greater than $`10^{19}`$ eV, and since its birth altitude was relatively low, there was little time for radiation.

### 5.2 Atmospheric cascading: PRESHOWER and CORSIKA

The results presented above give a general idea about the characteristics of preshowers initiated by UHE photons before entering the atmosphere. The next step is to simulate air showers induced by different preshowers or unconverted photons and compare the results with the proton-induced showers in order to investigate the signatures of UHE photons that will be measurable by the Pierre Auger Observatory. For this purpose we connected the PRESHOWER program to CORSIKA, a commonly used air shower simulation code.

The connection between the two codes is organized as follows. First, the distribution of the preshower particles at the top of the atmosphere is simulated with the original preshower code, then the output is passed to CORSIKA and finally the resulting air shower is simulated as a superposition of the subshowers initiated by the preshower particles.

Electromagnetic interactions are simulated in CORSIKA using an adapted version  of the EGS4 code , which includes the Landau-Pomeranchuk-Migdal effect . To reduce the computational effort in CPU time, the technique of particle thinning , including weight limitation , is applied. More specifically, a thinning level of $`10^5`$ has been chosen, and particle weights are limited to $`10^5E_0`$/GeV, with $`E_0`$ being the primary energy. This keeps artificial fluctuations introduced by the contribution of individual particles sufficiently small for the analysis of shower maximum depths .

In Table 3 we compare $`X_{max}`$ and RMS of $`X_{max}`$ for gamma induced showers of different primary energies and arrival directions. The $`X_{max}`$ of proton-induced showers using the high-energy hadronic interaction model QGSJET 01  is 820 (870) g/cm2 for $`10^{20}`$ ($`10^{21}`$) eV, and using the SIBYLL 2.1 model  855 (915) g/cm2. The RMS in all cases is about 50$``$60 g/cm2. Both for $`E_0=10^{19.5}`$ eV in strong $`B_{}`$ direction and for $`E_0=10^{20}`$ eV in weak $`B_{}`$ direction, almost every photon remains unconverted when entering the atmosphere, resulting in deep $`X_{max}`$ and large fluctuations due to the LPM effect. A comparison to the longitudinal profile of hadronic primaries shows that unconverted primary photons might be well distinguishable from primary protons or nuclei even on event-by-event basis.

Almost all the photons at $`10^{20}`$ eV arriving from the strong field direction convert into $`e^+e^{}`$ pairs. As the most energetic bremsstrahlung photons have energies several times lower than the primary photon energy, the subshowers develop faster. On average, the showers reach their maxima earlier than the one induced by an unconverted photon, but still not as early as hadronic showers. For larger statistics of 1020 eV events and if a substantial fraction of UHE cosmic rays are photons, one expects to see, as described above, the directional anisotropy in $`X_{max}`$ and RMS($`X_{max}`$) and on this basis one may draw conclusions on the primary composition. The absence of this dependence would allow us to set upper limits on the photon flux.

At the primary energy of $`10^{21}`$ eV all photons convert, whatever the arrival direction. We still see the directional dependence of $`X_{max}`$ but it is not as pronounced as previously. The fluctuations of $`X_{max}`$ in this case are significantly lower (by about a factor 4) than for $`10^{20}`$ eV primaries. Table 3 shows that at $`10^{21}`$ eV it is more difficult to distinguish a photon primary on the basis of the $`X_{max}`$ value from a proton one on an event-by-event basis. It is interesting, however, to note the directional dependence of the elongation rate. From Table 3 we see that at Malargüe, for the strong $`B_{}`$ direction, the elongation rate of photon-induced showers between $`10^{20}`$ eV and $`10^{21}`$ eV is much less than for proton or iron showers. For the weak $`B_{}`$ we even have a negative elongation rate. This is because the preshowering effect for photons at $`10^{21}`$ eV splits the initial energy into energies less than $`10^{20}`$ eV and at this energy level, for the weak $`B_{}`$ direction, almost all the primary photons remain unconverted and they induce air showers with deeper $`X_{max}`$. Low or negative elongation rates could be an additional good signature of photon showers.

More detailed simulations of the signatures mentioned above will allow us to quantify the expected sensitivity of the Auger experiment to a photon fraction in UHE cosmic ray flux. Also other observables will be studied, especially the muon content and the features of the lateral distribution measured by the surface detector, since for the surface measurements the event sample will be about 10 times more numerous than for the fluorescence detector.

## 6 Summary

We presented a method for studying the signatures of UHECR gamma showers and the practical implementation of this method in the PRESHOWER program. This program was tested extensively both in stand-alone mode and as a part of CORSIKA and was checked for consistency with results quoted in other publications. The preshower option will be accessible in the future release of CORSIKA.

With the code working in stand-alone mode, we studied the preshower characteristics for different primary energies and arrival directions. The strong directional and energy dependence of the preshower effect requires a concise treatment in the simulations. Anisotropies in the preshower cascade translate into anisotropies of air shower observables such as $`X_{max}`$ and RMS($`X_{max}`$). First results obtained with PRESHOWER combined with the CORSIKA air shower simulation code for the conditions of the Auger Observatory confirm the expectation that photon showers can be distinguished from hadron ones. Further analyses including other observables such as the muon shower content are in progress.

The new experiments with large collecting apertures give promising prospects for acquisition of a substantial amount of extensive air shower data in the UHE range of cosmic-ray spectrum. An important step towards explaining the origin of UHE cosmic rays will be the observation of primary UHE photons or specifying stringent upper limits for the photon flux. This requires a detailed understanding of the UHE photon shower characteristics. The presented PRESHOWER code can serve as tool for such simulations.

## 7 Acknowledgements

We thank R. Engel, M. Kachelriess, S.S. Ostapchenko for valuable remarks and suggestions.

This work was partially supported by the Polish State Committee for Scientific Research under grant No. PBZ KBN 054/P03/2001 and 2P03B 11024 and by the International Bureau of the BMBF (Germany) under grant No. POL 99/013.

## A Compilation of procedures/functions

preshw(\*id, \*E\_gamma, \*the\_loc, \*phi\_loc, \*toa, \*glong, \*glat,
\*igrf\_year, \*iiprint, \*corsika, \*nruns, part\_out,
\*r\_zero, \*r\_fpp, \*N\_part\_out):
This is the main procedure of PRESHOWER which simulates the propagation of an UHE photon before its entering the Earth’s atmosphere. It takes into account the effect of $`\gamma `$ conversion into $`e^{+/}`$ pair and subsequent bremsstrahlung radiation of the electrons in the geomagnetic field. As a result a bunch (a preshower) of particles is obtained, mainly photons and a few electrons, instead of the primary photon. The information about all the particles in the preshower is saved (in the stand-alone mode) or returned to CORSIKA (PRESHOWER-CORSIKA mode). The output parameters r\_zero, r\_fpp and N\_part\_out are required by CORSIKA only – they are not used by the stand-alone version.
INPUT:

1. \*id: primary particle type: always id=1 (a photon)
2. E\_gamma: primary energy \[GeV\]
3. \*the\_loc, \*phi\_loc: zenith and azimuth angles (in radians) of the shower trajectory given in the local coordinate system
4. \*toa: the top of the atmosphere \[cm\]
5. \*glong: longitude position of experiment \[deg\] (Greenwich glong = $`0^{}`$, eastward is positive)
6. \*glat: latitude position of experiment \[deg\] (North Pole: glat=$`90^{}`$, South Pole: glat=$`90^{}`$)
7. \*igrf\_year: the year in which the magnetic field is to be computed – this parameter is an input for the igrf procedure
8. \*iiprint: print flag (adjustable only in PRESHOWER-CORSIKA, in the stand-alone version iiprint=1)
= 0: suppress printing from preshw;
= 1: enable printing from preshw
9. \*corsika: mode flag
= 1: PRESHOWER-CORSIKA: – random generator from CORSIKA is used (rndm);
= 0: stand-alone version: random generator from Numerical Recipes (ran2) is used
10. \*nruns: number of runs

OUTPUT:

1. part\_out: the array to store the preshower particle data, maximum number of particles: 10000
part\_out\[n\]: particle type: 1 – photon, 2 – $`e^+`$, 3 – $`e^{}`$
part\_out\[n\]: particle energy \[GeV\]
part\_out\[n\]\[2-6\]: not used
2. \*r\_zero: height above sea level of the simulation starting point given in \[km\] above sea level
3. \*r\_fpp: height of the first pair production (0 for surviving photons) given in \[km\] above sea level
4. \*N\_part\_out: number of preshower particles at the top of the atmosphere (number of entries in part\_out array)

bessik(x, xnu, \*ri, \*rk, \*rip, \*rkp):
for argument x computes modified Bessel function rk of fractional order xnu;

getrand(seed, mode):
for a given seed seed returns a random number from interval (0,1) choosing an appropriate generator according to the mode flag: ran2 for the stand-alone mode (mode=0) and rmmard for the CORSIKA-PRESHOWER (mode=1)

float kappa(x):
an auxiliary function used in the bremsstrahlung procedure brem

brem(B\_tr, E\_e, E\_sg):
computes bremsstrahlung radiation for the transverse magnetic field B\_tr, electron energy E\_e and bremsstrahlung photon energy E\_sg

B\_transverse(B, tn):
getting the geomagnetic field B component which is perpendicular to the preshower trajectory tn

ran2(seed):
the random number generator used only in the stand-alone mode;

cross(a, b, c):
returns cross product c of vectors a and b

dot(a, b):
returns dot product of vectors a and b

norm(a, n):
the input vector a is normalized and returned as vector n

normv(a):
returns the length of vector a

glob2loc(glob, loc, theta, phi):
converts global cartesian coordinates glob to the local ones (loc) for the geographical position theta (colatitude) and phi (longitude)

sph2car(sph, car):
converts spherical coordinates sph to the cartesian ones car.

car2sph(car, sph):
inverse to sph2car

locsphC2glob(locsph, glob, theta, phi, sing, cosgm, sitec):
converts local spherical coordinates locsph (CORSIKA frame of reference) to the global cartesian coordinates glob, using geographical coordinates of the observatory (theta \- colatitude, phi \- longitude), global cartesian coordinates of the observatory (sitec) and sine and cosine of the declination angle (sing, cosg)

ppp(efrac, B\_tr, E\_g0):
this auxiliary function is used within ppfrac

ppfrac(s, mode, B\_tr, E\_g0):
determines the energy partition after gamma conversion, the computations are performed according to Ref. ; s \- seed, mode \- determines the type of the random number generator, B\_tr \- transverse magnetic field, E\_g0 \- primary photon energy

igrf\_(year, n, r, theta, phi, Br, Btheta, Bphi):
external Fortran procedure written by Tsyganenko ; computes magnetic field according to the IGRF model; used within b\_igrf; year \- the year of the experiment, n \- the highest order of the spherical harmonics in the scalar potential, r, theta, phi \- spherical coordinates of the point at which the field is to be found, Br, Btheta, Bphi \- spherical coordinates of the magnetic field at given position

b\_igrf(year, r, theta, phi, bcar):
for a given year and position (r, theta, phi) computes the magnetic field bcar in global cartesian coordinates according to the IGRF model

bsph2car(colat, longit, bsph, bcar):
converts the magnetic field from spherical coordinates (bsph) into cartesian ones bcar; used only in b\_igrf

rndm():
calls CORSIKA random number generator, used only in the PRESHOWER-CORSIKA mode


text:
# About “On certain incomplete statistics” by Lima et al

## Acknowledgments

We would like to thank J.A.S. Lima for useful discussions and for sending us important references.


text:
# Probing the variability of the fine-structure constant with the VLT/UVESBased on observations made with ESO Telescopes at the La Silla or Paranal Observatories under programme ID 066.A-0212.

## 1 Introduction

Modern 10 m class telescopes equipped with instruments like the High-Resolution Echelle Spectrograph (HIRES) at the Keck Observatory or the UV-Visual Echelle Spectrograph (UVES) at the ESO Very Large Telescope (VLT) facilitate the accurate observation of QSO absorption (or emission) lines in order to study the hypothetical variability of fundamental physical constants like the fine-structure constant $`\alpha =e^2/(4\pi ϵ_0\mathrm{}c)`$ or the proton-electron mass ratio $`\mu =m_\mathrm{p}/m_\mathrm{e}`$. The interest in these studies is motivated by the unification theories incorporating varying fundamental constants (for a review see, e.g., Uzan, 2003).

From the astronomical point of view the cosmological variability of the fine-structure constant is assessed as

$$\mathrm{\Delta }\alpha /\alpha =(\alpha _z\alpha _0)/\alpha _0,$$
(1)
where $`\alpha _0`$ and $`\alpha _z`$ denote the values of the fine-structure constant in the laboratory and the specific absorption (or emission) line system at redshift $`z`$, respectively. While observational studies based on the fine-structure splitting of intergalactic alkali-doublet (AD) absorption lines (e.g., Levshakov, 1994; Murphy et al., 2001b) and intrinsic QSO emission line doublets (Bahcall et al., 2003) have provided robust upper bounds on $`|\mathrm{\Delta }\alpha /\alpha |`$, Murphy et al. (2003a, hereafter MFW) have recently detected a non-zero expectation value $`\mathrm{\Delta }\alpha /\alpha `$ for a sample of 143 complex metal absorption systems identified in QSO spectra recorded with the Keck/HIRES: $`\mathrm{\Delta }\alpha /\alpha =(5.7\pm 1.1)10^6`$ in the redshift range $`0.2<z<4.2`$. This remarkable statistical evidence for a cosmological variation of the fine-structure constant is achieved by means of the many-multiplet (MM) technique, which is a generalization of the AD method incorporating the multi-component profile decomposition of many transitions from different multiplets of different ionic species (Dzuba et al., 1999; Murphy et al., 2003b, hereafter MWF, and references therein). While the MM technique considerably improves the formal accuracy, it also shows some immanent deficiencies. Levshakov (2003) illustrates in detail that the prerequisite assumption that the spatial distribution is the same for all ionic species is not valid in typical intergalactic absorption systems. Consequently, it is more reliable to apply the MM technique to samples of absorption lines arising from only one ionic species. Furthermore, in comparison to the AD calculations the incorporation of more transitions over a wider wavelength range results in a stronger susceptibility to systematic effects. Nevertheless, sources of error like wavelength miscalibration, spectrograph temperature variations, atmospheric dispersion, and isotopic or hyperfine-structure effects do demonstrably not explain the detected non-zero expectation value (Murphy et al., 2001a, MWF, MFW). The observational discrepancy grows since the radioactive decay rates of certain long-lived nuclei deduced from geophysical and meteoritic data provide a stringent bound, $`\mathrm{\Delta }\alpha /\alpha =(8\pm 8)10^7`$, back to the epoch of Solar system formation, $`z0.45`$ (Olive et al., 2003).

In the optical spectroscopy, the most striking deficiency inherent to all decomposition techniques are non-linear inter-parameter correlations preventing the accurate optimization of the model parameters and possibly causing ambigous results. In fact, the spectral resolution attained in QSO observations at the 10 m class telescopes is still not sufficient to resolve the metal lines with an expected minimum thermal width of about 1 km s-1. Clearly, in order to solve the profile decomposition problem spectral observations with the highest possible resolution and the highest possible signal-to-noise ratio are desirable.

In this study, we present exceptional high-resolution and high signal-to-noise spectra of the notably bright intermediate redshift QSO HE 0515–4414 (Reimers et al., 1998, $`z=1.73`$, $`B=15.0`$) recorded with the VLT/UVES. The spectra reveal a multi-component complex of metal absorption lines associated with a sub-damped Lyman-$`\alpha `$ (sub-DLA) system at the redshift $`z=1.15`$ (de la Varga et al., 2000; Quast et al., 2002; Reimers et al., 2003). We analyze a homogenous subsample of Fe ii $`\lambda 1608`$, $`\lambda 2344`$, $`\lambda 2374`$, $`\lambda 2383`$, $`\lambda 2587`$, and $`\lambda 2600`$ absorption lines by means of the standard and a revised MM technique in order to assess the cosmological variability of the fine-structure constant.

## 2 Observations

HE 0515–4414 was observed with UVES during ten nights between October 7, 2000 and January 3, 2001. Thirteen exposures were made in the dichroic mode using standard settings for the central wavelenghts of 3460/4370 Å in the blue, and 5800/8600 Å in the red. The CCDs were read out in fast mode without binning. Individual exposure times were 3600 and 4500 s, under photometric to clear sky and seeing conditions ranging from 0.47 to 0.70 arcsec. The slit width was 0.8 arcsec providing a spectral resolution of about 55 000 in the blue and slightly less in the red. The ThAr lamp exposures taken immediately after each science exposure provide an accurate calibration in wavelength. The standard deviation of the wavelength versus pixel dispersion solution is about 2.0 mÅ (2.5 mÅ) in the blue (red), resulting in an absolute accuracy of about 0.15 km s-1 in radial velocity space. The raw data frames were reduced at the ESO Quality Control Garching using the UVES pipeline Data Reduction Software. The calibrated spectra were converted to vacuum wavelengths according to Edlén (1966) while the barycentric velocity correction was manually cross-checked using the ESO-Munich Image Data Analysis Software and the Image Reduction and Analysis Facility. The individual vacuum-barycentric corrected spectra were manually cleaned from cosmic rays or pixel defects, rescaled to a common median flux level, and resampled to an equidistant wavelength grid using natural cubic spline interpolation. The combined spectra (Fig. 1) show an effective signal-to-noise ratio per pixel typically better than 100.

## 3 Analysis

### 3.1 The standard MM technique

The observed spectral flux $`F(\lambda )`$ is modelled as the product of the background continuum $`C(\lambda )`$ and the absorption term convoluted with the instrumental profile, i.e., $`F(\lambda )=C(\lambda )P(\xi )\mathrm{e}^{\tau (\lambda \xi )}d\xi .`$ While the background continuum is locally approximated by an optimized linear combination of Legendre polynomials of up to second order, the instrumental profile is modelled by a normalized Gaussian given by the the spectral resolution of the instrument. Assuming pure Doppler broadening, the optical depth $`\tau (\lambda )`$ is a superposition of Gaussian functions

$$g_i(\lambda )=\frac{e^2}{4ϵ_0mc}\frac{N_if_i\lambda _i}{\sqrt{\pi }b_i}\mathrm{exp}\left[\left(c\frac{\lambda Z_i\lambda _i}{Z_i\lambda _ib_i}\right)^2\right],$$
(2)
where $`Z_i=(1+z_i)(1+u_i/c)`$ while $`\lambda _i`$, $`f_i`$, $`z_i`$, $`u_i`$, $`b_i`$, and $`N_i`$ denote, respectively, the systemic rest wavelength, the oscillator strength, the cosmological redshift, the radial velocity, the line broadening velocity, and the column density corresponding to the line $`i`$. The systemic rest wavenumber $`\omega _i=1/\lambda _i`$ is parametrized as

$$\omega _i=\omega _{0,i}+q_i(\mathrm{\Delta }\alpha /\alpha )(2+\mathrm{\Delta }\alpha /\alpha ),$$
(3)
where $`\omega _{0,i}=1/\lambda _{0,i}`$ is the wavenumber in the laboratory and $`q_i`$ is the relativistic correction coefficient (Dzuba et al., 2002).

Even though the sub-DLA system exhibits many additional metal absorption lines (Al ii, Al iii, Mg i, Mg ii, Si ii, Cr ii, Ni ii, and Zn ii) typically incorporated in the standard MM analysis, we decide to apply this technique to the Fe ii transitions only. The restriction to one ionic species avoids systematic effects if the spatial distribution is not the same for all species, and reduces systematic effects arising from isotopic line shifts. In addition, the set of Fe ii transitions between $`\lambda 1608`$ and $`\lambda 2600`$ already provides a very sensitive combination for probing the variability of the fine-structure constant.

Each absorption component is modelled by a superposition of Doppler profiles with identical radial velocities, widths, and column densities. In addition, the value of $`\mathrm{\Delta }\alpha /\alpha `$ is confined to be the same for all components. In order to find the optimal set of parameter values for both $`\tau (\lambda )`$ as well as $`C(\lambda )`$, the weighted sum of residual squares is minimized by means of an evolution strategy (ES) based on the concept of covariance matrix adaption (Hansen et al., 2003). We ignore the ensemble of Fe ii $`\lambda 2383`$ and $`\lambda 2600`$ lines associated with the central region of the sub-DLA system since the profiles are saturated and are otherwise overemphasized in the optimization.

### 3.2 The regression MM (RMM) technique

The standard MM technique can conveniently be revised to avoid the deficiencies pointed out in Sect. 1. In fact, this technique is essentially similar to the method developed by Varshalovich & Levshakov (1993) in order to infer the cosmological variability of the proton-electron mass ratio $`\mu `$ from the analysis of molecular hydrogen absorption lines. Argueing by analogy (Levshakov, 2003), in the regime $`|\mathrm{\Delta }\alpha /\alpha |1`$ we obtain the linear approximation

$$z_i=z_\alpha +\kappa _\alpha Q_i,$$
(4)
where $`z_i=\lambda _{\mathrm{obs}}/\lambda _{0,i}`$ and $`Q_i=q_i/\omega _{0,i}`$ denote, respectively, the observed redshift and the sensitivity coefficient corresponding to the line $`i`$, and the slope parameter $`\kappa _\alpha `$ is given by

$$\kappa _\alpha =2(1+z_\alpha )(\mathrm{\Delta }\alpha /\alpha ).$$
(5)
If $`\mathrm{\Delta }\alpha /\alpha `$ is non-zero, $`z_i`$ and $`Q_i`$ will be correlated and we will be able to estimate the slope $`\kappa _\alpha `$ and the intercept $`z_\alpha `$ from the linear regression analysis of the position of the line centroids in an absorption component. The accuracy of the regression analysis will be improved, if several absorption line samples are combined. In this case, the regression procedure can be generalized appropriately:

$$\zeta _i=(\mathrm{\Delta }\alpha /\alpha )(Q_i\overline{Q}),$$
(6)
where $`\overline{Q}`$ denotes the mean sensitivity coefficient of the sample, and $`\zeta _i=(\overline{z}z_i)/[2(1+\overline{z})]`$ is the normalized redshift while $`\overline{z}`$ denotes the mean redshift of an absorption component (in radial velocity space $`\zeta _i=(\overline{u}u_i)/[2(c+\overline{u})]`$).

In order to determine the central position of several selected Fe ii lines (see Sect. 4.2) we basically follow the same strategy as described in Sect. 3.1. The only differences are that we do not incorporate Eq. (3) and do not confine the radial velocities to be the same for the lines in the selected components. We point out explicitly that even though we apply a parametric profile decomposition technique to determine the position of the line centroids, the RMM analysis is not tied to any specific modelling technique. In principle, the position of the line centroids can even be determined without doing any modelling at all (cf. Levshakov, 2003). Throughout the analysis we use the atomic data listed in Table 1.

## 4 Results and discussion

### 4.1 Standard MM analysis

The optimized values of the model parameters and the standard deviations provided by the covariance matrix of the ES are listed in Table 2. The decomposition of the Fe ii absorption complex is quite evident (see Figs. 1 and 2). This contrasts with the QSO absorption systems considered in former MM studies where many components are typically unresolved and many line profiles are saturated (see MWF, Figs. 3 and 4). Clearly, the components with the most accurately defined line centroids (No. 1, 12, 20–24, 26) are the most important in the MM analysis. Expectedly, the adequate profile decomposition is reflected in the formal accuracy of the result: $`\mathrm{\Delta }\alpha /\alpha _{\mathrm{MM}}=(0.1\pm 1.7)10^6`$ is the most stringent formal bound hitherto infered from an individual QSO absorption system. The best formal accuracy achieved in former MM analyses (see MWF, Table 3) is exceeded by a factor of about three.

### 4.2 Regression MM analysis

The optimized central positions of the lines considered in the RMM analysis are listed in Table 3. We do not consider component No. 22 since the central positions of the lines in this component are strongly correlated with those in No. 21 and we ignore component No. 24, because the central position of the Fe ii $`\lambda 1608`$ line is not accurately defined. The regression line (Fig. 3) indicates no correlation between the relative displacement of the lines, $`\zeta `$, and the sensitivity coefficient $`Q`$: $`\mathrm{\Delta }\alpha /\alpha _{\mathrm{RMM}}=(0.4\pm 1.9\pm 2.7_{\mathrm{sys}})10^6`$, where the statistical errors of $`\zeta `$ and the estimated errors of $`Q`$ as well as the systematic uncertainties inherent to the wavelength calibration are propagated by means of Monte Carlo simulation. The coefficient of determination is $`r^2=0.005`$, i.e. only half a percent of the variation among the relative displacement of lines is accounted for by the difference in sensitivity. Consulting the $`t`$ statistic, the observed data support the null hypothesis $`\mathrm{\Delta }\alpha /\alpha =0`$ at a significance level of 91 percent, whereas the support for the MFW result $`\mathrm{\Delta }\alpha /\alpha =5.710^6`$ is 12 percent.

### 4.3 Systematic effects

Contributing to the discussion of potential systematic effects presented in MWF and MFW, we recall attention to another important source of error. As formerly illustrated by Levshakov & D’Odorico (1995) and Levshakov (1994), the presence of unresolved narrow lines with different optical depths can strongly affect the position of the line centroids in an ensemble of lines and result in a biased expectation value $`\mathrm{\Delta }\alpha /\alpha `$. This effect will be most noticeable in the case of optically thick lines. In fact, the expectation value $`\mathrm{\Delta }\alpha /\alpha _{\mathrm{MM}}`$ increases by $`3.610^6`$ if we do not ignore the ensemble of saturated Fe ii $`\lambda 2383`$ and $`\lambda 2600`$ profiles in the optimization, and the systematic increase of scatter in the normalized residuals with increasing optical depth (see Figs. 4 and 5 provided in the online material) may be explained by the presence of unresolved narrow lines.

## 5 Conclusions

Our results strongly support the null hypothesis of a non-varying fine-structure constant, but do not contradict the MFW result $`\mathrm{\Delta }\alpha /\alpha =5.710^6`$ at a significance level higher than $`88`$ percent. Nevertheless, we conclude: (i) The MM technique has the capability to provide stringent bounds on $`\mathrm{\Delta }\alpha /\alpha `$ even if the analysis is restricted to the Fe ii lines only. (ii) The RMM technique is illustrative and methodically more transparent than the standard MM technique. In addition, the regression analysis facilitates the consideration of systematic errors inherent to the wavelength calibration of QSO spectra. (iii) The accuracy of the individual $`\mathrm{\Delta }\alpha /\alpha `$ assessment is principally limited by systematic errors inherent to the wavelength calibration. The accuracy attainable with high-quality QSO spectra recorded with the VLT/UVES is limited to $`10^5`$. (iv) Optically thick profiles are susceptible to systematic effects biasing the expectation value $`\mathrm{\Delta }\alpha /\alpha `$. (v) The analysis of an extensive homogenous sample of Fe ii absorption lines is inevitable and will provide an independent and crucial test of the MFW result.

In particular, HE 0515–4414 is the brightest known intermediate redshift QSO in the sky and is therefore predestinated for spectroscopy with the new High Accuracy Radial velocity Planet Searcher (HARPS) operated at the ESO La Silla 3.6 m telescope. This very high-resolution spectrograph is specified to provide an efficient wavelength calibration facilitating the performance of radial velocity measurements with an accuracy of better than 1 m s-1 (Pepe et al., 2002). Presumably, HARPS will improve the accuracy of the individual $`\mathrm{\Delta }\alpha /\alpha `$ assessment by more than an order of magnitude.

###### Acknowledgements.

We kindly thank Michael Murphy for valuable comments on the manuscript. This research has been supported by the Verbundforschung of the BMBF/DLR under Grant No. 50 OR 9911 1. S. A. Levshakov has been supported by the RFBR Grant No. 03-02-17522.