text:
# Large Scale Simulations of Jets in Dense and Magnetised Environments

## 1 Introduction

Several billion years ago, at redshifts in excess of two, the centers of galaxy clusters, typically hosting already a relaxed elliptical galaxy with an old population of stars, were usually equipped with a powerful radio jet. There are many different lines of evidence for this (Carilli et al. 2001). The host galaxies of radio jets are the brightest ones at their redshift. Also, there have been found dozens of line emitting objects around five high redshift radio galaxies, so far (Venemans et al. 2003). This is a significant overdensity. Furthermore, the space density of galaxy clusters at low redshift agrees with that of high redshift radio galaxies. Hence it is clear that radio galaxies pinpoint the most massive structures in the high redshift universe.

Contrary to the situation at high redshift, the brightest cluster galaxies (BCGs) in the local universe are generally associated only with weak radio jets. The only exception being Cygnus A (Fig. 2a). This classical double radio galaxy has an outstanding power, only reached by sources with redshifts in excess of roughly unity. The affected gas surrounding the radio jet can be observed in great detail in the X-ray regime by the Chandra satellite. Therefore, Cygnus A can serve as a model for the interaction of the jet with the intergalactic medium (IGM), for powerful, classical double radio jets. Much work has been done on the propagation of extragalactic jets (compare e.g. Krause & Camenzind 2001; Krause 2003, and references therein). Of particular interest for the present study are simulations that include all the external gas that is affected by action of the jet (Clarke et al. 1997; Reynolds et al. 2001, 2002; Saxton et al. 2002a, b; Krause & Camenzind 2002b, c; Krause 2002a, 2003; Zanni et al. 2003). These simulations did not yet reach the actual size of the jet in Cygnus A, but could be extrapolated to derive a hydrodynamical model of the jet (Krause 2002b; Krause & Camenzind 2002c). In this model, the jet in Cygnus A resides in a stratified galaxy cluster atmosphere, with a central density contrast (jet/IGM) of roughly $`10^4`$. We have tested this model with a large scale simulation that we report in sect. 2. The radio emission of jets has been applied to study the large scale coherence of magnetic fields in galaxy clusters (e.g. Dolag et al. 2002). We propose the idea that the magnetic field is randomly oriented, but with roughly the same plasma $`\beta `$ everywhere, before the jet passes. We have already shown, that a slab jet with similar parameters as used here, produces a coherent field structure in the shocked ambient gas (unpublished), with reversal scales of 10 to 50 $`R_\mathrm{j}`$, when simulated up to a diameter of 400 $`R_\mathrm{j}`$. In order to check the validity of this result in 3D, we tried to accomplish such a simulation. The results are described in sect. 3.

### 1.1 Numerics

For the computations in this contribution, the magneto-hydrodynamic (MHD) code Nirvana\_CP was employed. The main part of this code (NIRVANA) was written by Udo Ziegler Ziegler & Yorke (1997). In that version, it solves the MHD equations in three dimensions (3D) for density $`\rho `$, velocity $`\stackrel{}{v}`$, internal energy $`e`$, and magnetic field $`\stackrel{}{B}`$:

$`{\displaystyle \frac{\rho }{t}}+\left(\rho \stackrel{}{v}\right)`$ $`=`$ $`0`$ (1)
$`{\displaystyle \frac{\rho \stackrel{}{v}}{t}}+\left(\rho \stackrel{}{v}\stackrel{}{v}\right)`$ $`=`$ $`p\rho \mathrm{\Phi }+{\displaystyle \frac{1}{4\pi }}\left(\stackrel{}{B}\right)\stackrel{}{B}{\displaystyle \frac{1}{8\pi }}\stackrel{}{B}^2`$ (2)
$`{\displaystyle \frac{e}{t}}+\left(e\stackrel{}{v}\right)`$ $`=`$ $`p\stackrel{}{v}`$ (3)
$`{\displaystyle \frac{\stackrel{}{B}}{t}}`$ $`=`$ $`\times (\stackrel{}{v}\times \stackrel{}{B}),`$ (4)
where $`\mathrm{\Phi }`$ denotes an external gravitational potential.

NIRVANA can be characterised by the following properties:

1. explicit Eulerian time–stepping,
2. operator–splitting formalism for the advection part of the solver,
3. method of characteristics–constraint–transport algorithm to solve the induction equation and to compute the Lorentz forces;
4. artificial viscosity has been included to dissipate high–frequency noise and to allow for shock smearing in case the flow becomes supersonic.

The code was vectorised and parallelised by OpenMP like methods, and successfully tested on the SX-5 (Krause & Camenzind 2002a). All the significant loops could be vectorised. The number crunching part scales without significant performance loss. This is also true for the MHD part of the solver. We show a typical profile output below (Tables 1 and 2), for a run without data output, which indicates an optimum in vectorisation and parallelisation efficiency. The average performance in the 2D simulation was only 434 cumulative MFLOPS with eight processors, probably because in this run, about 500 GB of data had to be dumped to the hard disk, which is a serial process.

## 2 Simulation of a Very Light Jet to Large Scale

We have performed a bipolar axisymmetric simulation – in the following called run A – of a very light jet in a King type galaxy cluster atmosphere with a final jet size of $`>200`$ jet radii. The simulation was run for 6200 CPU hours on eight processors of the SX-5 at the HLRS.

### 2.1 Simulation Setup

The jet was injected in both directions (bipolar) in the center of a cylindrical grid with 8300 and 2000 points in axial (Z) and radial (R) direction, respectively. The basic length scale of the problem is the jet radius which is represented by 20 points. With that resolution global parameters like the bow shock velocity on the Z-axis or energy and momentum conservation are accurate to $``$ 10% (Krause & Camenzind 2001). As our unit of length we choose the observed jet radius in Cygnus A, i.e. the jet radius is set to $`R_\mathrm{j}=0.5`$ kpc. The total grid size is therefore $`[207.5\times 50]`$ kpc. The environmental density ($`\rho _\mathrm{e}`$) is given by an isothermal King profile:

$$\rho _\mathrm{e}(R,Z)=\rho _{\mathrm{e},0}\left(1+\frac{r^2}{a^2}\right)^{3\beta /2},$$
(5)
where $`r^2=R^2+Z^2`$. This means that the density is constant up to a core radius a that was set to $`a=10`$ kpc. Then it starts to decrease, asymptotically reaching $`\rho _\mathrm{e}r^{3\beta }`$ ($`\beta =3/4`$). The central density was set to $`\rho _\mathrm{e}=m_\mathrm{H}\mathrm{cm}^3`$, with the hydrogen mass $`m_\mathrm{H}=1.67\times 10^{24}`$ g. The temperature was set to $`T=30`$ Mio. K. This atmosphere is kept in hydrostatic equilibrium by the gravity of a dark matter halo:

$$\mathrm{\Phi }=\frac{3\beta kT}{2\mu m_\mathrm{H}}\mathrm{log}\left(1+\frac{r^2}{a^2}\right).$$
(6)
$`\mu `$ is the number of particles per proton mass. Here, $`\mu =0.5`$ for an ionised medium. In order to brake the symmetry, density perturbations were included, i.e. with 10% probability, the density in a cell was increased by a random factor between 0 and 40%. The jet density was set to $`\rho _\mathrm{j}=\eta _0\rho _\mathrm{j}`$, with the density contrast $`\eta _0=10^4`$. The sound speed in the jet was set to $`20\%c`$, $`c`$ being the the speed of light, and the jet’s Mach number to $`M=3`$.

The simulation was run for 20 Myr, in total. During that time, the jet reached an extention of 110 kpc, i.e. 220 jet radii, on the axis.

### 2.2 Results

We present logarithmic density plots of the simulation results for four different simulation times ($`5,10,15,20`$ Myr) in Fig. 3. The morphology that appears in these figures is a continuation of previous simulations that could not reach the size shown here. Fig. 3a shows the state that was reached by Krause (2003), and extensively discussed therein. In this early phase, the bow shock is spherical, its radius following an expansion law given by the force balance equation which can be integrated to yield, for arbitrary mass distribution $`(r)`$ and energy injection $`E(t)`$:

$$_0^r(r_1)r_1𝑑r_1=2_0^t𝑑t_1_0^{t_1}E(t_2)𝑑t_2.$$
(7)
For the given matter profile (5), (7) can be integrated numerically. We only discuss the asymptotic power law parts here. For a power law density distribution ($`\rho =\rho _0(r/r_0)^\kappa `$) and constant energy injection ($`E=Lt`$), the solution is:

$$r=\sqrt[\kappa +5]{\frac{(\kappa +3)(\kappa +5)r_0^\kappa Lt^3}{12\pi \rho _0}}$$
(8)
The density profile used here has the asymptotic power law approximations: $`lim_{r0}(\rho )=\rho _0`$, and $`lim_r\mathrm{}(\rho )r^{9/4}`$. Therefore, the bow shock should expand with $`rt^{0.6}`$ in the beginning, steepening towards $`rt^{1.09}`$, at least as long as it remains spherical. The radial bow shock position was determined every 0.4 Myr (compare Fig 2c).

It was fitted with a function $`r=a+bt^c`$, where a,b, and c were simultaneously varied. At $`5`$ Myr, where the bow shock is still almost spherical (aspect 1.2), $`c`$ is 0.82. The expected exponent, using a local power law approximation for the density and the analytic approximation, is: $`c_{\mathrm{theo}}=0.79`$. This confirms the analytic approximation. But the bow shock continues to follow this expansion law far beyond the spherical phase: at $`17.6`$ Myr, the fit gives $`c=1.06`$ versus 1.00 from the analytic model. Even a global power law, with exponent 0.86 fits quite good.

The propagation in axial direction is shown in Fig 2d. Besides the early times, the bow shock always accelerates. Its acceleration ($`vt^{0.2}`$) is always significantly higher than predicted by self-similar modelling (asymptotically: $`vt^{0.09}`$) (compare e.g. Carvalho & O’Dea 2002, and references therein), but does not reach the prediction for narrow beams of constant radius that would be super-exponential for $`\kappa <2`$.

This reflects the state of the beam plasma. Figure 3 shows a sound beam, reaching to the tip of the bow shock, only for $`t<10`$ Myr. At later phases , the beam is disrupted towards the head. Thereby, the beam thrust is distributed over a greater area, leading to lower bow shock velocity compared to the case of a narrow beam with constant radius.

The cocoon transforms gradually from conical to cylindrical at later times. Its width does not grow much for the second half of the simulation. Towards the center, there is a turbulent region, where the ambient gas is entrained into the cocoon by the action of gravity, and mixed with the shocked beam plasma.

### 2.3 Comparison to Observations

The prime target to compare the simulation data to is the X-ray data of Cygnus A (compare Fig. 2a). The simulation was run to almost the size of the real source (110 kpc and 140 kpc, for simulation and observation, respectively). It is evident from Fig. 3 that this is necessary in order to get the correct cocoon morphology. Many other details can be understood by this new simulation: Fig. 2a shows that the radio cocoon is bordered by the brightest X-ray emission. We show the line-of-sight integrated X-ray emission for the simulation together in Fig. 3d. These can be compared directly to the observations in Fig. 2. The simulation clearly shows the peaks in emission next to the radio cocoon. This was already predicted analytically, in self-similar models (Alexander 2002). Outwards of this peak, the emissivity is increased with respect to the undisturbed cluster emission, by a constant factor in the plane of symmetry. The emission falls suddenly at the bow shock. This detail cannot be seen in the observational data, because there are not enough photons (but compare Fig. 4 in Smith et al. (2002)). Fig. 3d shows that the emission of the shocked ambient gas has been shaped by the bow shock in an elliptical way. Smith et al. (2002) have shown that the elliptical isophotal fit is better than the spherical fit inside a radius of $`66/(h/0.7)`$ kpc, which should indicate the bow shock’s radial position in Cygnus A.

The simulation shows significant emission in the mixing region in the central parts of the cocoon. However, this is less than observed. The reason is numerical mixing with cocoon gas in that region. Also, the observation shows dominantly non-axisymmetric modes, which cannot be represented in our axisymmetric simulation. We note that this region has been found to be non-axisymmetric in 3D simulation (Krause & Camenzind 2002a).

It has been suggested that repeated jet episodes could heat the surrounding gas to the observed X-ray temperatures during cosmological timescales. A lower limit to the heat injection in the cluster gas was derived by Krause (2003), equ. (21). For the simulation parameters here, this lower limit amounts to $`3.19\times 10^{58}`$ erg. The cluster gas was followed by a passive tracer that was set to unity. Counting only cells where the tracer variable is above 0.1, we measure $`8.17\times 10^{59}`$ erg internal energy injected by the jet into the cluster gas by the end of the simulation. The cluster gas has also gained $`2.37\times 10^{60}`$ erg of potential energy. In total, 69% of the energy injected by the jet has been put into the ambient gas. This is enough to power the X-ray emission of the cluster in Cygnus A for $``$500 Mio. years. This makes it entirely plausible that the cluster gas is heated by that mechanism.

We have proposed to measure jet parameters based on the radial bow shock propagation. The radial bow shock velocity in Cygnus A can be constrained from the shock temperature measured by Chandra (Krause & Camenzind 2002c). This procedure gives a jet power of $`L=8\times 10^{46}`$ erg/s and an age of $`t=27`$ Myr for Cygnus A. In order to demonstrate the validity of the procedure, we determine the jet parameters for the simulation in the same way. The simplest approximation for the external density is taking it to be constant. Then, from (8) it follows with the bow shock velocity of 1220 km/s, the radial bow shock position 30.17 kpc and the external density of $`10^{25}`$ g/cm3: $`L=9.96\times 10^{45}`$ erg/s and $`t=19.7`$ Myr. This is to be compared to the true jet power, $`L_{\mathrm{true}}=\pi R_\mathrm{j}^2\rho _\mathrm{j}v_\mathrm{j}^3=7.72\times 10^{45}`$ erg/s and the true jet age of 20 Myr. This very good agreement demonstrates that the exact shape of the cluster atmosphere is not critical in determining jet parameters. The agreement can even be improved by taking better approximations to the density profile.

The simulation result shows an important difference to the observation: In the simulation, the beam is very unstable, reaching the tip of the bow shock not even once, after 10 Myr. The reason for this is the low Mach number in the beam that is a consequence of the strong interaction with the cocoon and the entrained shocked ambient gas therein. Higher Mach numbers can only be reached by a relativistic jet. However, the beam could also be stabilised by an appropriate, significant magnetic field. The magnetic field is also demanded in order to preserve the contact discontinuity near the tip of the bow shock from the action of the Kelvin-Helmholtz-instability, because the strong disruption found in the simulation can not be found in the radio data. This result is in good agreement with with magnetic field determinations in Cygnus A’s hot spots via the self-synchrotron-Compton assumption (Wilson et al. 2000) and our earlier suggestion, based on a jet power argument, that the jet’s mean Lorentz factor is $`\mathrm{\Gamma }20`$.

### 2.4 Comparison to High Redshift Radio Galaxies

So far, high redshift radio galaxies have not been observed long enough in order to study the cluster gas emission in great detail. But the region where the shocked ambient gas could be expected is typically bright in emission lines. The emission line regions often have the same cone shaped structure as the X-ray emission in Cygnus A (compare Fig 2a and 2b). The line emission is brightest in the region, corresponding to the mixing region in the simulation. This could be interpreted in the way that in these objects the line emission is caused by material that was entrained into the radio cocoon.111The alternative interpretation is that cooling is significant in the whole shocked ambient gas region. This would lead to very large cocoon width and turbulent mixing of radio plasma and emission line gas (Krause 2002a). Most recent X-ray data for 4C 41.17 (Scharf et al. 2003), showing probably inverse Compton cocoon emission, point in this second direction for that particular object. This gas is subject to the combined thermal and Rayleigh-Taylor instability, which may cool some gas to the appropriate temperatures (Basson 2002; Krause 2002b).

## 3 A Jet in a Randomly Magnetised Environment

We have accomplished a 3D simulation of a jet in a randomly magnetised environment. The simulation run for 500 CPU hours on eight processors of the SX-5. Unfortunately, the timestep became too low, before significant evolution of the jet could be seen.

### 3.1 Numerical Setup

We injected the jet in the center of a Cartesian grid ($`[X\times Y\times Z]=[170\times 74\times 74]\mathrm{kpc}=[512\times 224\times 224]`$) in the X direction. The jet radius was set to$`R_\mathrm{j}=1`$ kpc and resolved with 3 points. The magnetic field at the jet inlet was set to zero for the $`Y`$ and the $`Z`$ direction, and $`B_\mathrm{X}/x=0`$ on the jet nozzle. In the surrounding King atmosphere ($`\rho _{\mathrm{e},0}=m_\mathrm{p}/10`$, $`a=10`$ kpc , and $`\beta =2/3`$, $`T=3\times 10^7`$ K), a random magnetic field was established. The vector potential was randomly determined, folded with the King distribution in order to get an average plasma $`\beta `$ of $`\beta =8\pi p/B^2=12`$, i.e. sub-equipartition fields. The initial density contrast was set to $`\eta =10^3`$ and the Mach number to $`M=5`$.

### 3.2 Results and Discussion

A slice of the density distribution and the plasma $`\beta `$ is shown in Fig. 4. In that phase, the bow shock is already weak. Weak MHD shocks reduce the field strength and refract the field towards the shock normal. The simulation shows a decrease of the magnetic field of typically an order of magnitude, behind the bow shock. yet. Unfortunately the simulation could not be run for longer time, since the timestep became too low. In that phase, no amplification of the reversal scale could be found, as in a slab simulation (compare sect. 1). This probably happens only when the shocked ambient gas expands to several times the cocoon diameter, like e.g. in Fig. 3d. If this effect could be confirmed, it would be a good candidate for the explanation of the high reversal scale in galaxy clusters, i.e. the coherence of the field is created by the weak bow shock.

## 4 Conclusions

Employing the NEC SX-5 supercomputer at the HLRS, we could simulate an axisymmetric jet to the large extention of more than 200 jet radii. On that scale, the results can be directly compared to observational data of Cygnus A. Many details, like elliptically shaped emission, morphology of cocoon and shocked ambient gas, or low aspect, are well reproduced. This confirms the hypothesis that the jet is underdense with respect to its environment by a factor of $`10,000`$, and the consequences on the jet parameters discussed above. We also find differences to the observational data, concerning stability of the beam and the contact discontinuity, which we ascribe to the missing magnetic fields and the disregard of relativistic physics. We speculate that emission line halos of high redshift radio galaxies might be identified with the centrally concentrated mixing region in the cocoon.

We also tried to simulate the influence of the jet on a randomly magnetised environment. Unfortunately, the timestep became too low, and the simulation had to be stopped, before significant progress has been made.

## Acknowledgements

This work was also supported by the Deutsche Forschungsgemeinschaft (Sonderforschungsbereich 439).

## Index


format: text

text:
# Emergence of pulled fronts in fermionic microscopic particle models

## Abstract

We study the emergence and dynamics of pulled fronts described by the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation in the microscopic reaction-diffusion process $`A+AA`$ on the lattice when only a particle is allowed per site. To this end we identify the parameter that controls the strength of internal fluctuations in this model, namely, the number of particles per correlated volume. When internal fluctuations are suppressed, we explictly see the matching between the deterministic FKPP description and the microscopic particle model.

The study of diffusion-limited reaction processes has shown the important role of internal or microscopic fluctuations in low dimensions marro ; avraham ; hinrichsen . Mean-field approximation for those processes assume that diffusive mixing is much stronger than the influence of correlations produced by reactions. However, diffusive mixing is not strong enough in low-dimensional systems and fluctuations might modify the dynamics or induce nonequilibrium phase transitions. While this behavior is observed in different situations, there is an special interest in the problem of front propagation in reaction-diffusion systems avraham ; doering1 ; riordan ; moro ; saarloos ; brunet ; breuer ; kessler ; avraham1 ; panja ; tripathy ; bramson86 ; kerstein ; kerstein88 ; lemarchand . In this paper we concentrate in microscopic lattice reaction-diffusion models whose mean-field approximation is given by the FKPP equation fisher

$$\frac{\rho }{t}=D\frac{^2\rho }{x^2}+k_1\rho k_2\rho ^2,$$
(1)
where $`\rho (x,t)`$ is the local concentration of particles. Such an equation display traveling-wave solutions of the form $`\rho =\rho (\xi )`$ with $`\xi =xvt`$ which invade the unstable phase $`\rho =0`$ from the stable phase $`\rho =k_1/k_2`$ and travel with velocity $`vv_0=2\sqrt{Dk_1}`$. For steep enough initial conditions, the solution selected for large times is the one with minimal velocity, $`v_0`$, which is known to be a pulled front, since it is essentially “pulled along” by the growth and spreading of small perturbations in the leading edge where $`\rho 1`$ ebert ; saarloos . Microscopic fluctuations are expected to modify macroscopic properties of pulled fronts at two levels: (i) first because the deterministic description (1) breaks down at small densities $`\rho 1/N`$ where $`N`$ is the number of particles, which introduces an effective cutoff in the FKPP equation. Due to the importance of the tail development in pulled fronts, several front features are dramatically affected by this effective cutoff brunet : for example, the selected velocity converges as $`\mathrm{ln}^2N`$ to the mean-field value $`v_0`$, (ii) second, because internal fluctuations are present and could interplay with or even destroy pulled front development and dynamics doering1 ; avraham1 ; moro .

One of the most studied microscopic models is the reversible reaction model $`AA+A`$ on a lattice doering1 ; riordan ; avraham1 ; breuer ; kessler ; panja ; bramson86 ; kerstein ; kerstein88 ; lemarchand ; tripathy ; moro . In the bosonic version of this model breuer , the number of possible particles per site is unbounded and thus the balance between birth and coagulation gives an average number of particles per site $`N`$. If $`N\mathrm{}`$ the reaction is well stirred within each site and the front dynamics is described by the mean-field approximation (1). For very large $`N`$ discreteness effects remain and produce the predicted velocity correction $`vv_0\mathrm{ln}^2N`$ brunet ; kessler . In the case of the fermionic version of the $`AA+A`$ model only a particle is allowed per site. The main reason to consider exclusion is that, for some values of the parameters, the model is analytically tractable avraham ; doering1 ; avraham1 ; bramson86  and/or simulations are easier than in the bosonic version. Exact results are available for the two interesting regimes in the model: reaction-limited regime bramson86 , where coarse-grained density front profiles are described by the mean-field FKPP equation, and diffusion-limited regime avraham ; doering1 ; avraham1  in which internal fluctuations dominate front propagation and the mean-field approximation (1) is not valid. Our purpose in this paper is to put these two results in a general framework that can describe the emergence of pulled fronts in this fermionic model. This is done by identifying the control parameter that modulates the effect of internal fluctuations on the front propagation model. As we will see, this parameter also controls the development of the tail front and establishes the appearance of pulled fronts.

In the $`AA+A`$ model in one dimension, particles are allowed to occupy lattice sites and can undergo the following moves: (i) Diffusion to any one of its two neighbor lattice sites with a diffusion rate $`D`$ (ii) Birth: occupied sites spontaneously generate particles at neighbor lattice sites with rate $`\mu `$, (ii) Coagulation: a particle can get annihilated with death rate $`\eta `$ if one of its two neighboring filled lattice sites is occupied. The fermionic nature of the model makes diffusion and birth only possible if the neighboring lattice site is empty. The mean-field description of this model is given by the FKPP equation with $`k_1=2\mu `$, $`k_2=2(\mu +\eta )`$. Starting from an initial condition in which occupation number is only different fron zero on the right side of a site, a front develops and advances as a function on time. Operationally, the front position $`x_f(t)`$ is determined by a local average of density of particles over intervals of length $`\lambda ^1=2D/v_0`$, which is the width of the deterministic front selected for our initial condition by the FKPP equation (1) ebert . Thus, front position is the point where this coarse-grained density equals $`\rho _0/2`$ tripathy ; moro . Other definitions of the front position yield the same results bramson86 . Since results only depend on the ratios $`D/\mu `$ and $`D/\eta `$ we set $`D=1/2`$ throughout this paper. After a transient time (which could be long), the front advances linearly, i.e. $`x_f(t)=vt`$, where $`\mathrm{}`$ stands for average over different realizations. In this regime, statistical properties of the front with respect to the normalized coordinate $`\xi =xvt`$, are independent of time.

In this model, there are two special cases for which exact results are available: when $`\eta =D`$ the model is solvable using the method of inter-particle distribution functions doering1 ; avraham . In that case fronts advance with velocity $`v=\mu `$ which shows how internal fluctuations can dominate the system behavior in this diffusion-limited regime. On the other hand, when $`\eta =0`$ it was proved in bramson86  that fronts approach asymptotically the FKPP equation predictions ($`v=v_0`$) in the limit $`D/\mu \mathrm{}`$ (reaction-limited regime), while $`v=D+\mu `$ in the opposite regime $`D/\mu 0`$ (diffusion-limited regime). In Fig. 1 we show the results of our simulations for the velocity of the front $`v`$ as a function of $`\mu `$ for different values of $`\eta `$. Our results are consistent with the exact results doering1 ; avraham ; bramson86  and previous simulations of this model kerstein . For an intermediate case $`0<\eta <D`$ we observe that, for some values of $`\mu `$, the velocity seems to approach the deterministic value $`v_0`$. However, for small enough value of $`\mu `$ internal fluctuations seem to dominate and the velocity deviates strongly from $`v_0`$.

To understand this behavior, let us recall that mean field approximation (1) in the $`AA+A`$ is only valid when diffusive mixing is strong enough. Specifically, this happens when the typical distance traveled diffusively by a particle between reaction events, $`l_D`$, is much larger that the typical distance between particles, $`l_p`$ avraham . In that case, particles are well stirred within cells of size $`l_D`$ and thus, mean-field approximation is valid for the coarse-grained density of particles over cells of size $`l_D`$ as shown in bramson86 . In our model we have that $`l_D=\mathrm{min}(\sqrt{D/\mu },\sqrt{D/\eta })`$, while we approximate $`l_p`$ by the average distance between particles in the stable phase $`l_p=\rho _0^1=(\mu +\eta )/\mu `$. Since we are interested in the propagation of pulled fronts, which are only driven by the birth term, our condition to approach the mean-field approximation is then given only by $`\mu `$:

$$\sqrt{D/\mu }\rho _0^1=\frac{\mu +\eta }{\mu }.$$
(2)
Interestingly, this condition is equivalent to $`N^{}\lambda ^1\rho _01`$, where $`N^{}`$ is approximately the number of particles within an interval of length $`\lambda ^1`$. Thus condition (2) also means that the number of particles within the typical length scale of the front (its width) is large. Our results in Fig. 1 are then easily explained in terms of condition (2): when $`N^{}1`$ internal fluctuations should be unimportant within cells of size $`l_D\lambda ^1`$ and front propagation should approach asymptotically the FKPP predictions. Since

$$N^{}=\frac{\sqrt{\mu /D}}{\mu /D+\eta /D},$$
(3)
then $`N^{}1`$ only happens when $`(D/\eta )^2D/\mu 1`$ for fixed values of $`D`$ and $`\eta `$. We show the condition $`N^{}>1`$ in Fig. 2 along with the different set of parameters used in Fig. 1 and in other works doering1 ; avraham1 ; panja . Outside the region $`N^{}>1`$ internal fluctuations dominate and fronts are not described by the FKPP equation. This is the case for $`\eta =D`$ doering1 ; avraham  and $`\eta =\mu `$ panja . In the intermediate case $`0<\eta <D`$ we can have values of $`\mu `$ for which $`N^{}`$ is relatively large and fronts seem to approach the deterministic value of $`v_0`$, which explains the behavior observed in Fig. 1 for $`\eta =0.002`$. Note however, that although being in the region $`N^{}>1`$ is the minimum requirement for our model to approach the mean-field description (1) a finite value of $`N^{}`$ means that fronts are still subject to internal fluctuations and discreteness effects which produce a (strong) correction to the velocity. Only in the limit $`N^{}\mathrm{}`$ do both effects become negligible and the $`AA+A`$ system is effectively described by the FKPP equation. This is the case for $`\eta =0`$, $`\mu 0`$ note .

An interesting question is whether $`N^{}`$ plays any role like the average number of particles per site, $`N`$, in bosonic models brunet ; kessler ; breuer . In those models, it is observed that the deterministic description of a pulled front given by the FKPP equation is valid until the density drops to $`\rho N^1`$ which produces an effective cutoff in the tail of the front and modifies its velocity brunet . To check this possibility, we have measured in our simulations the average distance of the last particle from the front position, $`\xi ^{}`$, which is observed to saturate to a constant value for long enough times. It is obvious that for $`\xi >\xi ^{}`$ the continuum description of the front breaks down and we expect this to happen when there is only a particle in each coarse-grained site of length $`\lambda ^1`$, i.e. when $`\rho (\xi ^{})a\lambda `$, where $`a`$ is a constant. When internal fluctuations are irrelevant, i.e. when $`N^{}1`$, we assume that the continuum description (1) is still valid up to $`\xi ^{}`$ and taking that $`\rho \rho _0\lambda \xi e^{\lambda \xi }`$ for $`\lambda \xi 1`$ for a pulled front ebert  we obtain

$$\lambda \xi ^{}e^{\lambda \xi ^{}}=a/N^{}.$$
(4)
Solutions of this equation for $`\xi ^{}`$ with $`N^{}`$ given by (3) are compared with our simulations in Fig. 3. We see that for $`\lambda \xi ^{}1`$ Eq. (4) gives a rather accurate prediction of $`\xi ^{}`$. This corroborates our assumption that a pulled front described by the FKPP equation develops even for moderate values of $`N^{}`$ up to the point where $`\rho (N^{})^1`$.

An important consequence of Eq. (4) is that $`N^{}`$ controls not only the size of internal fluctuations but also the appearance and length of the tail in the pulled front. Thus, when $`N^{}1`$ internal fluctuations dominate and also the tail length is roughly zero, $`\lambda \xi ^{}0`$. This means that the front is basically a shock wave with height $`\rho _0`$. Actually, the exact solution when $`\eta =D`$ avraham1  shows that the particles behind the leading one remain distributed as in the stable phase ($`\rho =\rho _0`$) at all times, which confirm our picture (see Fig. 3). This shock wave shape is also observed in the case $`\eta D`$ when $`\mu \eta `$, where it is found that $`\lambda \xi 0`$. In the intermediate case in which $`\eta 0`$ we see that the front develops a tail which is described by the FKPP equation only for a given interval of values of $`\mu `$ (see Fig. 3 for $`\eta =0.002`$).

In the case $`\eta =0`$ we have studied the correction to the velocity as a function of $`N^{}`$ and observe in the inset Fig. 1 that it decays like $`(v_0v)/v_0(N^{})^{1/3}`$ which is consistent with simulations of other microscopic bosonic models lemarchand ; brunet  for moderate values of number of particles $`N`$. This results stresses the equivalence of $`N^{}`$ with the role that the number of particles plays in other microscopic models. In particular we expect the correction to be $`(vv_0)/v_0\mathrm{ln}^2N^{}`$ for very large values of $`N^{}`$ note1 .

Since the last particle is, on average, at a certain distance from the front position, their velocities coincide. This fact was used in kerstein ; panja  to estimate the velocity of the front by counting possible forward and backward hopping rates:

$$v=\mu \rho ^{}(\eta D)$$
(5)
where $`\rho ^{}`$ is the probability of having a particle behind the last one. Several approximations can be made for the value of $`\rho ^{}`$ kerstein ; panja . For example, in panja  is was taken as $`\rho ^{}\rho _0`$, i.e. $`\rho ^{}`$ is given by the probability to find a particle in the stable phase. Clearly, this approximation is only valid in the case in which fronts are like a shock-wave, i.e. when $`N^{}1`$ because then the last particle is very close to the stable phase. In the case in which a pulled front develops ($`N^{}1`$) we find that Eq. (5) still holds: since the last particle is on average at a fixed distance from the front position, we can approximate $`\rho ^{}\rho _0\lambda \xi ^{}e^{\lambda \xi ^{}}=a/\lambda `$ which is the concentration of particles at $`\xi ^{}`$. Our simulations for $`\eta =0`$ confirm the validity of this approximation (see Fig. 3) which brings out the effective matching between the continuum description given by the FKPP equation for $`\xi <\xi ^{}`$ and the microscopic character of the model for $`\xi \xi ^{}`$.

Another interesting property of front propagation is the wandering of the position of the front around its mean value $`\mathrm{\Delta }^2(t)=(x_f(t)x_f(t))^2`$. Several studies of different microscopic models shows that $`\mathrm{\Delta }^2(t)=2D_ft`$ for long times breuer ; brunet ; lemarchand  and that the diffusion coefficient $`D_f`$ depends on the number of particles $`N`$. Specifically, it was found that $`D_fN^{1/3}`$ for moderate values of $`N`$ breuer ; lemarchand , while $`D_f\mathrm{ln}^3N`$ for very large values of $`N`$ brunet . Our simulations for the $`A+AA`$ system shows that fronts move diffusively for all values of the parameters. In the case $`\eta =0`$, in which the model approaches asymptotically the FKPP equation, we get $`D_f(N^{})^{1/3}`$ (see Fig. 4) like in bosonic models lemarchand  which stress once again the fact that $`N^{}`$ plays the role of the number of particles in this fermionic model. Moreover, we found for small times that as $`N^{}`$ increases, the correlation between the time development of the front and internal fluctuations produces superdiffusive motion of the front position $`\mathrm{\Delta }^2(t)t^{2\nu }`$ with $`\nu 0.8`$. Once the front tail is developed (which happens at $`tk_1^1`$), front position starts to wander diffusively. Finally, our results for the diffusion of the front indicate that as the front approaches the deterministic FKPP equation, internal fluctuations make the front move diffusively at times $`t>k_1^1`$, independently of $`N^{}`$. We do not observed any signs of subdiffusive behavior conjectured by some authors tripathy ; saarloos  for pulled fronts subject to noise. This supports the idea that pulled fronts subject to internal noise belong to a different universality class than those subject to external noise moro .

In summary, we have identified the parameter that controls the strength of microscopic fluctuations for the front propagation problem in the fermionic model $`AA+A`$, namely the number of particles $`N^{}`$ per coarse-grained site of length $`\lambda ^1`$. When $`N^{}1`$, internal fluctuations are suppressed and the front becomes a pulled front like those of the FKPP equation (1). Moreover, our studies about the length of the tail, the velocity of the front and its diffusion show that $`N^{}`$ plays the same role as the number of particles in other microscopic bosonic models. Finally, it is interesting to note that in the $`AA+A`$ model, the velocity of a macroscopic object such as the front is related to the microscopic motion of the last particle, something also observed in other works brunet ; panja ; kerstein88 . We hope our results will help to understand the dynamics of fronts in microscopic fermionic reaction-diffusion models and its relevance when discussing properties of the FKPP equation subject to internal noise tripathy ; saarloos ; moro .

We are grateful to D. ben-Avraham, E. Brunet, R. Cuerno, J. Casademunt, C. Doering and W. van Saarloos for comments and discussions. Financial support is acknowledged from the European Commission through its Marie Curie program and from the Ministerio de Ciencia y Tecnología (Spain).


text:
# Multiphoton Exchange Processes in Ultraperipheral Relativistic Heavy Ion Collisions

## 1 Introduction

In ultraperipheral relativistic heavy ion collisions (UPC) very strong electromagnetic fields are present for a very short time . The parameter which characterizes the strength of the interaction is the Coulomb parameter $`\eta `$:

$$\eta =\frac{Z_1Z_2e^2}{\mathrm{}v}Z_1Z_2\alpha $$
(1)
where $`Z_1`$ and $`Z_2`$ are the charges of the nuclei and $`vc`$ is the relative velocity of the ions. For heavy nuclei like Au-Au (RHIC) or Pb-Pb (LHC) $`\eta 1`$. One can treat these systems semiclassically, with the nuclei following classical Rutherford trajectories. At high energies they are well approximated by straight lines with impact parameter $`b`$. In addition to elastic scattering, which is due to the exchange of many photons, there are various kinds of inelastic processes, like nuclear excitations, photon-nucleus interactions leading to baryonic or mesonic excitations and photon-photon processes. One can calculate the amplitudes for these processes in the same way as the corresponding ones for real photons .

Probabilities for the excitation of the giant dipole resonance (GDR) and $`e^+e^{}`$-pair production are very large, of the order of one for small impact parameter . This means that these processes will also occur simultaneously with other interesting processes, where the probabilities are generally much smaller than one. The GDR is a strongly collective mode in nuclei. The theory of relativistic heavy ion excitation of one- and multiphonon GDR states is given in . There is also mutual excitation, where each ion is excited to the (one phonon) GDR state. This process has been recently studied experimentally at RHIC ; it is important for the luminosity determination and for triggering on UPCs. The fact that the double phonon giant dipole resonance is excited strongly in a two-photon excitation mechanism is of special interest for nuclear structure physics. The properties of these new collective modes could be explored with this method .

The probability of electron-positron pair production calculated in lowest order perturbation theory exceeds the unitarity limit of one . This means that multiple $`e^+e^{}`$ pair production is expected to be appreciable for heavy ions at RHIC and LHC conditions.

The probability $`P(b)`$ of other inelastic processes is in general much smaller than one. An important process is diffractive vector meson production (especially $`\rho ^0`$-production). Probabilities for $`\rho ^0`$-production in close (grazing) collisions are about 0.01 or 0.03 for RHIC and the LHC, respectively. The probability of the multiple (double) electromagnetic production of $`\rho ^0`$ mesons is of the order of $`10^3`$ to $`10^4`$.

One example of simultaneous processes is giant dipole resonance excitation (followed by neutron decay) accompanied by $`\rho ^0`$-production. This process was recently studied experimentally at STAR.

This paper will present a general, unified treatment of multiphoton processes in heavy ion collisions, focusing on the correlations which occur because all of the processes share a common impact parameter, $`\stackrel{}{b}`$. Section 2 reviews the semiclassical treatment and then introduces multiphoton interactions. Section 3 presents some examples of multiphoton processes and gives examples of the correlations, with RHIC and the LHC used for illustration. We give our conclusions in Section 4. A preliminary account of part of the present work is given in Ref. .

## 2 Theory

We work in the rest frame of one of the nuclei (“target system”). The other nucleus moves on a straight line with impact parameter $`b>R_1+R_2`$, (where the nuclear radii are $`R_1`$ and $`R_2`$, respectively). The corresponding Lorentz-factor is $`\gamma =1/\sqrt{1v^2/c^2}`$. It is related to the Lorentz-factor $`\gamma _{coll}`$ in the collider frame by $`\gamma =2\gamma _{coll}^21`$.

Since the strength of the electromagnetic interaction decreases with increasing impact parameter $`b`$, the excitation probability also decreases, in many important cases proportional to $`1/b^2`$. Higher order effects are proportional to the product of such probabilities, and are therefore especially important for collisions with impact parameter $`b`$ close to $`b_{min}=R_1+R_2`$.

The collision time is given by $`\tau _{coll}=b/\gamma `$, so the electromagnetic spectrum will include frequencies up to

$$\omega _{max}\frac{\gamma }{b_{min}}.$$
(2)
### 2.1 Elastic Coulomb Collisions

A simple and rather accurate result is well known for classical Coulomb scattering at small angles $`\theta 1`$ . The momentum transfer is

$$\mathrm{\Delta }p=\frac{2Z_1Z_2e^2}{bv}=\frac{2\eta \mathrm{}}{b}$$
(3)
in the direction perpendicular to the motion. The momentum transfer is related to the scattering angle $`\theta `$ by $`\mathrm{\Delta }pp\theta `$. From this classical relation between the impact parameter and the scattering angle one can calculate the classical relativistic Rutherford scattering cross section as

$$\frac{d\sigma }{d\mathrm{\Omega }}=\left[\frac{2Z_1Z_2e^2}{pv}\right]^2\frac{1}{\theta ^4},$$
(4)
where $`p=\gamma M_1v`$ and $`M_1`$ is the mass of nucleus $`Z_1`$, see, e.g., p. 644 of . This result can also be obtained in a field theoretical approach, which we only sketch here. Coulomb scattering is due to the exchange of photons. The sum over all the ladders diagrams, crossed and uncrossed, can be done in the high energy limit and for forward angles . A typical diagram is shown in Fig. 1(b).

One ingredient in this procedure is the linearization of the propagator denominator

$$\frac{1}{(p+q)^2m^2}\frac{1}{2pq}$$
(5)
In order to sum up all ladder and crossed ladder graphs one uses the identity

$$\underset{\sigma }{}\frac{1}{a_{\sigma (1)}(a_{\sigma (1)}+a_{\sigma (2)})\mathrm{}(a_{\sigma (1)}+a_{\sigma (2)}+\mathrm{}+a_{\sigma (n)})}=\frac{1}{a_1a_2\mathrm{}a_n}$$
(6)
where the sum is over all permutations $`\sigma `$ of the indices $`1`$ to $`n`$.

Summing over all orders the scattering amplitude is found similar to the eikonal approximation in (non relativistic) potential scattering:

$$f_{el}(q)=\frac{p}{2\pi i}d^2b\mathrm{exp}(i\stackrel{}{q}\stackrel{}{b})\left[\mathrm{exp}(i\chi (b))1\right].$$
(7)
By integrating over the angle $`\varphi `$

$$f_{el}(q)=ip_0^{\mathrm{}}𝑑bbJ_0(qb)\left[1\mathrm{exp}\left(i\chi (b)\right)\right].$$
(8)
In this approach there is a well defined momentum transfer $`q=\mathrm{\Delta }p`$ (we set $`\mathrm{}=1`$), and the impact parameter $`b`$ is not an observable quantity. For the elastic Glauber phase $`\chi (b)`$ in the case of photon exchange, the analytic expression is well known, see e.g. :

$$\chi (b)=2\eta K_0(\lambda b),$$
(9)
where $`\lambda `$ is a small photon mass needed to regularize the integral at large $`b`$.

For $`\eta 1`$ the integrand oscillates very rapidly and large contributions to the integral are only obtained when the phase of this integrand is stationary. One can introduce the saddle point (“stationary phase”) approximation

$$e^{i\mathrm{\Phi }(t)}𝑑te^{i\mathrm{\Phi }(t_0)}\sqrt{\frac{2\pi i}{\mathrm{\Phi }^{\prime \prime }(t_0)}}$$
(10)
where $`t_0`$ is determined by $`\mathrm{\Phi }^{}(t_0)=0`$. (One may have to sum over various stationary points $`t_0`$). In our case we have a two-dimensional integration over $`\stackrel{}{b}`$; for a derivation of the semiclassical case using only a one-dimensional saddle point approximation, see . Assuming that the momentum transfer $`\stackrel{}{q}`$ points in the $`x`$-direction, the phase in the integrand of Eq. (7) is (for small $`\lambda `$)

$$\mathrm{\Phi }(\stackrel{}{b})=q_xb_x+2\eta ln(\lambda b)$$
(11)
The conditions

$$\frac{\mathrm{\Phi }}{b_y}=2\eta \frac{b_y}{b^2}=0,\frac{\mathrm{\Phi }}{b_x}=q_x+2\eta \frac{b_x}{b^2}=0$$
(12)
lead on the one hand to $`b_y=0`$, that is the dominant contribution comes from the direction, where $`\stackrel{}{b}`$ is parallel to $`\stackrel{}{q}`$. It also leads to an extremum for $`b_x=b_0`$ with

$$q=\mathrm{\Delta }p=\frac{2\eta }{b_0},b_0=\frac{2\eta }{q}.$$
(13)
The second derivatives at the saddle point are

$$\frac{^2\mathrm{\Phi }}{b_x^2}=\frac{2\eta }{b_0^2},\frac{^2\mathrm{\Phi }}{b_xb_y}=0,\frac{^2\mathrm{\Phi }}{b_y^2}=\frac{2\eta }{b_0^2}.$$
(14)
Therefore an approximation of the elastic scattering amplitude up to a phase is

$$|f_{el}(q)|=\frac{pb_0^2}{2\eta }=\frac{2p\eta }{q^2},$$
(15)
giving the same cross section as the classical calculation.

Eq. (13) gives the classical connection of the momentum transfer $`q`$ to the impact parameter $`b_0`$, see above Eq. (3). The dominant contribution to the (pseudo-) Gaussian integral comes from values of $`\stackrel{}{b}`$ around the classical impact parameter $`b_0`$ in (opposite) direction of $`\stackrel{}{q}`$, with a width (“uncertainty of the $`b`$-value”, see Eq. (14)) of $`\frac{b_0}{\sqrt{\eta }}`$.

The classical momentum transfer $`\mathrm{\Delta }p`$ is built up from the perpendicular momenta of many exchanged photons. The perpendicular momentum of the individual photon is restricted to $`q_{}\begin{array}{c}<\hfill \\ \hfill \end{array}\frac{1}{b}`$, so $`2\eta `$ photons must be exchanged to reach the classical $`\mathrm{\Delta }p`$. This contrasts with electron or proton scattering, where $`\eta =\frac{1}{137}1`$ for $`vc`$. So, in the following, we use the semiclassical approximation to describe inelastic processes in the case $`\eta 1`$.

For $`b<R_1+R_2`$ one has to consider strong interactions. In many cases, especially for heavy ions, the “black disk approximation” is quite good. In this approximation $`\chi `$ is very large and imaginary for $`b<R_1+R_2`$. In the case of strong Coulomb elastic scattering ($`\eta 1`$), including the strong absorption of the ions for $`b<R_1+R_2`$, the diffraction pattern is given by Fresnel diffraction (rather than a Fraunhofer one for $`\eta <1`$). This is explained in detail in , see also Sec. 5.3.5 of . These grazing angles are extremely small, about $`5\times 10^5`$ (radians) for Au-Au collisions RHIC and $`2\times 10^6`$(radians) for Pb-Pb collisions at LHC (in the collider frame); they are well below the beam emittance, so the scattered ions will remain in the beam. Although the previous discussion is not entirely new, we think that it is useful to present the general ideas, which relate classical and quantal (or field-theoretical) descriptions of scattering in strong Coulomb fields, which may be a rather unfamiliar subject.

### 2.2 Perturbation Theory for Inelastic Processes

In the following we use the semiclassical approximation and treat the Coulomb field of the ion(s) as classical external fields. One ion moves on a straight line trajectory with a definite impact parameter, while the other ion is at rest. The semiclassical approximation can be justified in the same way as in the elastic case. The scattering amplitude in the inelastic case is given in the eikonal approximation by

$$f_{fi}(\stackrel{}{q})=\frac{p}{2\pi i}_0^{\mathrm{}}d^2b\mathrm{exp}(i\stackrel{}{q}\stackrel{}{b})\mathrm{exp}[i\chi (b)]a_{fi}(\stackrel{}{b}).$$
(16)
The amplitude $`a_{fi}`$ can be expanded in the form $`a_{fi}(\stackrel{}{b})=_\mu a_{fi}^\mu (b)\mathrm{exp}(i\mu \varphi )`$ where $`\mu `$ denotes the angular momentum transfer along the beam axis. After integration over the azimuthal angle $`\varphi `$,

$$f_{fi}(\stackrel{}{q})=ip_0^{\mathrm{}}𝑑b\underset{\mu }{}J_\mu (qb)i^\mu \mathrm{exp}[i\chi (b)]a_{fi}^\mu (b).$$
(17)
Using the saddle point approximation, the connection to the semiclassical case can clearly be seen and $`a_{fi}(\stackrel{}{b})`$ corresponds to the semiclassical excitation amplitude. In the semiclassical approximation, the scattering amplitude is the product of the elastic Coulomb scattering amplitude and the semiclassical excitation amplitude $`a_{fi}(b)`$. We have

$$f_{fi}(\stackrel{}{q})=f_{el}(q)a_{fi}(b_0\widehat{q})=f_{el}(q)\underset{\mu }{}(1)^\mu a_{fi}^\mu (b_0)$$
(18)
where $`\widehat{q}`$ is the unit vector pointing in the $`\stackrel{}{q}`$-direction (see Eqs. (11)–(14)). While we have not made numerical comparisons of the field theoretic solution to the semiclassical approximation in the relativistic case, such comparisons exist for the non-relativistic electromagnetic excitation. In this case the electromagnetic cross sections can be calculated exactly. They depend essentially on the Coulomb parameter $`\eta `$ which approaches infinity in the semiclassical approximation. It is found (see, e.g., Fig. (II.9) of ) that the asymptotic (semiclassical) limit is reached already for rather low $`\eta `$ values. The semiclassical approximation works better than one may expect. In Au-Au or Pb-Pb collisions $`\eta 50`$, so the semiclassical approximation should be excellent.

The electric charge of the relativistic ion gives rise to an electromagnetic potential, the Lienard-Wiechert potential $`A_\mu (\stackrel{}{r},t)`$

$`A_0(\stackrel{}{r},t)`$ $`=`$ $`\varphi (\stackrel{}{r},t)={\displaystyle \frac{Z_pe\gamma }{\sqrt{(bx)^2+y^2+\gamma ^2(zvt)^2}}}`$
$`\stackrel{}{A}(\stackrel{}{r},t)`$ $`=`$ $`{\displaystyle \frac{\stackrel{}{v}}{c}}\varphi (\stackrel{}{r},t).`$ (19)
The interaction with the target is described in terms of an electromagnetic current operator $`\widehat{j}_\mu `$. This leads to a time-dependent electromagnetic interaction of the form

$$V(t)=d^3rA_\mu (\stackrel{}{r},t)\widehat{j}^\mu (\stackrel{}{r}).$$
(20)
The dependence of $`V(t)`$ on the impact parameter is not shown explicitly. The target current $`\widehat{j}_\mu `$ describes a wealth of physics and is a quite complicated object in general. This current contains the current of the nucleus and the interaction leads to the excitation of nuclear states, with the giant dipole resonance (GDR) as an important example. It also describes photonuclear interactions like nucleon excitations, meson production and photon-photon interactions, for example lepton pair production. The lepton ($`e^+e^{}`$) states are those in the external Coulomb field of the target, i.e., the Furry picture is used here. Now one can expand the (time dependent) target state $`\mathrm{\Phi }(t)`$ in terms of all possible states $`\mathrm{\Phi }_n`$ which can be reached by the interaction $`V(t)`$. We write

$$\mathrm{\Phi }(t)=\underset{n}{}a_n(t)\mathrm{exp}(iE_nt)\mathrm{\Phi }_n,$$
(21)
where we have introduced the time-dependent amplitudes $`a_n(t)`$ for the state $`\mathrm{\Phi }_n`$. The (time-independent) states $`\mathrm{\Phi }_n`$ consist of all states which can be connected to the target nucleus ground state by the interaction $`V`$ (i.e., due to the interaction with the (virtual) photons). This includes, for example, excited nuclear states, the nucleus in its ground state and lepton pairs, or mesons, or nuclear excited states along with lepton pairs, or mesons, etc.

One can set up coupled equations for the amplitudes $`a_n(t)`$ . They are

$$i\dot{a}_n=\underset{m}{}n\left|V(t)\right|m\mathrm{exp}(i(E_nE_m)t)a_m(t)$$
(22)
The initial condition is $`a_m(t\mathrm{})=\delta _{m0}`$. The expression on the r.h.s. can be viewed as the matrix element of the operator

$$\stackrel{~}{V}(t)=\mathrm{exp}(iH_0t)V(t)\mathrm{exp}(iH_0t)$$
(23)
(“interaction representation”). A formal solution of this equation can be written down using the time-ordering operator $`𝕋`$:

$$a_n(t+\mathrm{})=n\left|𝕋\mathrm{exp}\left(i_{\mathrm{}}^+\mathrm{}𝑑t\stackrel{~}{V}(t)\right)\right|0.$$
(24)
We use perturbation theory to solve this equation. Since the electromagnetic interaction $`V`$ is in general weak, this is a good procedure. The first order amplitude is

$$a_n^{(1)}=i_{\mathrm{}}^{\mathrm{}}𝑑tV_{n0}(t)\mathrm{exp}(i\omega _{n0}t)$$
(25)
where $`V_{nm}(t)=<n|V(t)|m>`$ and $`\omega _{nm}=E_nE_m`$. This is the one-photon approximation. In second order

$$a_n^{(2)}=\frac{1}{i^2}\underset{m}{}_{\mathrm{}}^{\mathrm{}}𝑑tV_{nm}(t)\mathrm{exp}\left(i\omega _{nm}t\right)_{\mathrm{}}^t𝑑t^{}V_{m0}(t^{})\mathrm{exp}\left(i\omega _{m0}t^{}\right),$$
(26)
where the sum is over all intermediate states $`m`$. Only if the integrand is symmetric in $`t`$ and $`t^{}`$ one can extend the integral over $`t^{}`$ to infinity and factorize the amplitude. This is not possible here since the operator $`\stackrel{~}{V}`$ does not commute in general for different times $`t`$ and $`t^{}`$. (In the case of the excitation of a harmonic oscillator, the commutator of these operators is a $`c`$-number and a full analytical solution can be found; this case will be treated below in Sec. 2.4.)

In the case of the excitation of two independent modes, that is of two modes which do not mix with each other since their interaction is zero or negligible, there is such a factorization: let us assume that one excites a state $`n=(\alpha ,\beta )`$ where $`\alpha `$ denotes, e.g., a nuclear excited state and $`\beta `$ a vector meson in a state with a given momentum and helicity produced coherently. Two types of intermediate states contribute to the sum over $`m`$, see Fig. 2: $`m=(0,\beta )`$ and $`m=(\alpha ,0)`$. By summing over them the integrand becomes symmetric in $`t`$ and $`t^{}`$. Formally this can then be written in the following way: The first path goes from the state $`0=(0,0)`$ to $`(\alpha ,0)`$ and then to $`(\alpha ,\beta )`$. The contribution of this path to the second order matrix element is

$`a_{(\alpha ,\beta )}^{(2)}(1)=`$
$`{\displaystyle \frac{1}{i^2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑tV_{(\alpha ,\beta )(\alpha ,0)}(t)\mathrm{exp}\left(i\omega _{\beta 0}t\right){\displaystyle _{\mathrm{}}^t}𝑑t^{}V_{(\alpha ,0)(0,0)}(t^{})\mathrm{exp}\left(i\omega _{\alpha 0}t^{}\right).`$ (27)
and for the second path

$`a_{(\alpha ,\beta )}^{(2)}(2)=`$
$`{\displaystyle \frac{1}{i^2}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑tV_{(\alpha ,\beta )(0,\beta )}(t)\mathrm{exp}\left(i\omega _{\alpha 0}t\right){\displaystyle _{\mathrm{}}^t}𝑑t^{}V_{(0,\beta )(0,0)}(t^{})\mathrm{exp}\left(i\omega _{\beta 0}t^{}\right).`$ (28)
The independence of the two processes corresponds to the assumption that the presence of the state $`\alpha `$ (or $`\beta `$) does not influence the interaction matrix element for the production of $`\beta `$ (or $`\alpha `$) (and that a single photon cannot produce a single step transition from $`(0,0)`$ to $`(\alpha ,\beta )`$, see below):

$$V_{(\alpha ,\beta )(0,\beta )}=V_{(\alpha ,0)(0,0)}=V_{\alpha 0}$$

and

$$V_{(\alpha ,\beta )(\alpha ,0)}=V_{(0,\beta )(0,0)}=V_{\beta 0}.$$

In many important practical cases this independence should be very well fulfilled. For mutual excitation the factorization is well understood, the elastic form factor of the excited ion is almost unchanged from the one of the ion in its ground state. Even the slight change in form factor takes time, as will be discussed below in connection with the validity of the sudden limit.

Similarly for photon-photon processes like lepton-pair production together with nuclear excitation, the photon-photon process depends only on the charge of the nucleus, with little sensitivity to the details of the charge distribution. The nucleus is not influenced by the additional photon-photon process and vice versa. For hadronic final states, where there may be strong interactions with the nuclei, the photon-photon process in a double equivalent photon approach will predominantly occur outside the two nuclei.

Adding the amplitudes of the two paths, Eqs. (27) and (28) one can see that the integrand is symmetric in t and $`t^{}`$ and one obtains factorization. In an obvious way this is generalized to the excitation of $`N`$ independent modes. For any function $`f(t_1,t_2,\mathrm{},t_n)`$ which is symmetric in all variables $`t_i`$

$`{\displaystyle _{\mathrm{}}^t}𝑑t_1\mathrm{}{\displaystyle _{\mathrm{}}^t}𝑑t_nf(t_1,\mathrm{}t_n)=`$
$`n!{\displaystyle _{\mathrm{}}^t}𝑑t_1{\displaystyle _{\mathrm{}}^{t_1}}𝑑t_2\mathrm{}{\displaystyle _{\mathrm{}}^{t_{n1}}}𝑑t_nf(t_1,t_2,\mathrm{}t_n).`$ (29)
One can use this formula to factorize the amplitude into a product of first order amplitudes. The factor of $`\frac{1}{n!}`$ is compensated by the $`n!`$ different ways to order the excitation of the independent modes $`\alpha ,\beta ,\gamma ,\mathrm{}`$.

Let us give here also an example of two modes, which are coupled (i.e. no longer independent): the excitation of vibrational and rotational modes of a nucleus . These modes mix due to rotation-vibration coupling. This gives rise to nuclear eigenstates of a mixed character. In the case of $`\rho ^0`$ production and GDR excitation, for example, such a coupling is expected to be so small that it can safely be neglected.

The factorization is also not fulfilled in single photon exchange processes, which are depicted in Fig. 3. One such process is photon exchange with inelastic vertices on both sides (Fig. 3(a)). This process was studied in  within an inelastic EPA approach and was found to be only of the order of 1% as compared to the elastic process. Another example is the GDR excitation together with the meson production, see Fig. 3(b). This is a special contribution to the incoherent meson production process. Again this process is expected to be small: the vector meson is mainly produced through a diffractive (absorptive) process, which is very ineffective for a GDR excitation. The magnitude of such a process could be estimated by looking, e.g., at the GDR excitation in inelastic meson scattering at high energies. Unfortunately we are unaware of any study that has been done at these high energies.

Asymmetric systems, such as $`d`$-A are also of interest, as shown by the recent $`d`$-Au run at RHIC, where UPCs were studied . In these collisions, non-factorizable processes could come into play. The Coulomb parameter is $`\eta =0.6`$, less than one and the semiclassical approximation (used, i.e., in equations like Eqs. (3), (10) or (18)) does not apply. However, the total cross-sections (i.e. integrated over $`b`$ and $`q`$, respectively) are equal in the semiclassical and Glauber (or plane wave) case . In this run, the $`\rho ^0`$\+ nuclear breakup channel was studied. Due to the $`Z^2`$ dependence the photon(s) will come predominantly from the gold nucleus. The $`\rho ^0`$ \+ breakup channel gets the main contribution from incoherent rho production, i.e. it is a one photon process of the form of Fig. 3(b). There is also a two-photon contribution from coherent $`\rho ^0`$-production accompanied by an electromagnetic deuteron breakup . This mechanism is suppressed by a factor of $`P_{breakup}(b)`$ relative to the incoherent one-photon mechanism, which is of the order of 1% in the relevant impact parameter range. Due to the selection mechanism in the experiment the incoherent rho meson production mechanism is disfavored. Therefore overall the balance between the two reaction mechanisms should be studied in more detail.

### 2.3  Multiphoton Processes in the Sudden or Glauber Approximation

The solution of the coupled equations Eq. (22) is greatly facilitated if the sudden approximation can be applied, see, e.g., . In this case the collision time $`t_{coll}=b/\gamma `$ is much smaller than the excitation time $`1/\omega `$. This condition is fulfilled in many interesting cases and we assume now that the sudden approximation can be used. This “frozen nucleus”-approximation is also used in Glauber theory. The relation between the semiclassical approach and the Glauber (or eikonal) approximation is explained in , for the non-relativistic, as well as, for the relativistic case. The excitation amplitude is

$$a_n(t\mathrm{})=<n|\mathrm{exp}(iR)|0>$$
(30)
where

$$R=_{\mathrm{}}^+\mathrm{}𝑑tV(t)=_{\mathrm{}}^+\mathrm{}𝑑tA_\mu (t)\widehat{j}^\mu .$$
(31)
The operator $`R`$ is a direct sum of operators in the space of nuclear states, the space of the nucleus-vector meson system, the nucleus-$`e^+e^{}`$ system, etc.. So,

$$R=\underset{i}{}R_i$$
(32)
where $`i`$ denotes the different final states like $`e^+e^{}`$ pairs, $`\rho `$-mesons or a GDR excited ion. One has

$$\mathrm{exp}(iR)=\underset{i}{}\mathrm{exp}(iR_i).$$
(33)
as the different $`R_i`$ commute with each other. One may expand the exponential in Eq. (30). Terms linear in $`R`$ describe, e.g., the excitation of nuclear states, like the collective giant dipole resonance (GDR), vector meson production or $`e^+e^{}`$pair production. Terms quadratic in $`R`$ give, e.g., contributions to double-phonon GDR excitation, double vector meson production, double $`e^+e^{}`$ pair production. The quadratic terms also describe, e.g., vector meson production and GDR excitation in a single collision. A contribution to the second order amplitude $`a^{(2)}`$ is, e.g.,

$$a^{(2)}=<GDR|R|0><\rho ^0|R|0>$$
(34)
where $`|0>`$ denotes the ground state of the nucleus. In the general case the $`R_i`$ commute and using Eq. (33) one gets

$$a=<GDR|\mathrm{exp}(iR_{GDR})|0><\rho ^0|\mathrm{exp}(iR_{\rho 0})|0>,$$
(35)
which has the second order term as a special limit.

For three independent processes, say GDR excitation, vector meson- and $`e^+e^{}`$ pair production, the six different time orderings compensate for the $`1/3!`$ factor in the expansion of Eq. (30). In this formalism, processes are independent and the elementary amplitudes factorize, as one would have intuitively expected. This property is used also in the experimental analysis of vector meson production with simultaneous GDR excitation, where the neutrons from the GDR decay serve as a trigger on UPC . The ion motion is not disturbed by the excitation process. The reason is that the kinetic energy of the ion is much larger than the excitation energy. Due to coherence, the quantity $`x=\frac{\omega }{E}=\frac{1}{RM_NA}`$ is very small. One has $`x<4\times 10^3,3\times 10^4`$ and $`1.4\times 10^4`$ for oxygen, tin and lead ions respectively.

The GDR excitation immediately affects the nucleon momenta, but the nucleons will take a finite time to adjust their positions. With an excitation energy of about 20 MeV, the oscillation time of the GDR can be estimated classically to be about 10 fm$`/c`$. Due to the time dilatation, this corresponds to a time in the collider frame of about 1000 fm$`/c`$ at RHIC and even longer at the LHC. This time is long compared to the time required to, for example, produce a $`\rho ^0`$.

From the general discussion above of higher order theory we have seen that this is not a necessary assumption to achieve factorization.

### 2.4  Excitation of a Harmonic Oscillator (Bosonic Modes)

An especially simple and important case is the excitation of a harmonic oscillator. It also has the virtue that it can be studied analytically. As some excitation processes in heavy ion collisions can be regarded as quasibosonic processes, such a model can serve as a good approximation. Some examples are discussed below. Although well known (see, e.g., ), we briefly present the main results of the excitation of a harmonic oscillator by an external force, this serves also to establish the notation.

In terms of the creation and destruction operators $`a^{}`$ and $`a`$ the Hamiltonian of the system is

$$H=\omega \left(a^{}a+\frac{1}{2}\right)$$
(36)
where $`\omega `$ is the energy of the oscillator. The boson commutation rules are $`[a,a^{}]=1`$ and $`[a,a]=[a^{},a^{}]=0`$. Only one mode is shown explicitly, in general one has to sum (integrate) over all the possible modes.

One assumes that the operator $`V`$ (see Eq. (20) is linear in the creation and destruction operators:

$$V(t)=f(t)a+f^{}(t)a^{}.$$
(37)
In this case one can calculate $`\stackrel{~}{V}`$ explicitly using the boson commutation rules given above and the expansion

$$\stackrel{~}{V}(t)=V(t)+it[H_0,V(t)]+\frac{(it)^2}{2!}[H_0,[H_0,V(t)]]+\mathrm{}$$
(38)
One finds

$$\stackrel{~}{V}(t)=f(t)e^{i\omega t}a+f^{}(t)e^{i\omega t}a^{}$$
(39)
Now one can convince oneself that the commutator of $`\stackrel{~}{V}`$ at different times $`t`$ and $`t^{}`$ is a pure $`c`$-number. In this case one can disregard the time ordering operator in Eq. (24) and obtain an exact analytical answer, up to an unimportant overall phase factor (see, e.g., ).

The excitation operator is the integral of the interaction over time

$$_{\mathrm{}}^+\mathrm{}𝑑t\stackrel{~}{V}=(ua+u^{}a^{})$$
(40)
where $`u`$ is the $`c`$-number

$$u=_{\mathrm{}}^{\mathrm{}}𝑑tf(t)\mathrm{exp}(i\omega t).$$
(41)
This leads to the excitation of a so-called coherent state, see . For the excitation of multiphonon states this is explicitly shown in . One has

$$a_n=<n|e^{i(ua^{}+ua)}|0>=\frac{(iu^{})^n}{\sqrt{n!}}e^{\frac{1}{2}uu^{}},$$
(42)
where the operator identity $`e^{A+B}=e^Ae^Be^{\frac{1}{2}[A,B]}`$ was used, which is valid for two operators $`A`$ ($`iu^{}a^{}`$) and $`B`$ ($`iua`$) for which the commutator is a $`c`$-number. See also the discussion of the forced linear harmonic oscillator in .

Let us look at a few cases where this model can be applied: Electromagnetic excitation of nuclear states, especially the collective giant multipole resonances was discussed recently by Bertulani . The possibility to excite multiphonon GDR states is studied in . Multiphonon GDR also play a role in the electromagnetic excitation of the ions in relativistic heavy ion collisions, for a recent reference see . The parameter which describes the probability $`P_{GDR}(b)`$ of GDR excitation is (see )

$$P(b)=\frac{S}{b^2}$$
(43)
with a simple parameterization for $`S`$ for the relevant high beam energies

$$S=\frac{2\alpha ^2Z_1^2N_2Z_2}{A_2m_N\omega }=5.45\times 10^5Z_1^2N_2Z_2A_2^{2/3}\text{fm}^2.$$
(44)
where $`m_N`$ denotes the nucleon mass, and the neutron-, proton-, and mass-number of the excited nucleus are $`N_2,Z_2,`$ and $`A_2`$ respectively. The excitation probability is inversely proportional to the energy $`\omega `$ $`(80\text{MeV}A^{1/3})`$ of the GDR state. As expected, the soft modes are more easily excited. Eq. (43) is based on the assumption that the classical dipole sum rule (Thomas-Reiche-Kuhn sum rule) is exhausted to 100%. For the excitation of an $`N`$-phonon state, a Poisson distribution is obtained. For heavy systems $`P(b)\frac{1}{2}`$ for close collisions ($`b\begin{array}{c}>\hfill \\ \hfill \end{array}R_1+R_2`$).

Double $`\rho ^0`$ production was studied in . In addition to the label $`m`$ for the magnetic substates of the GDR, there is a continuous label (the momenta). The probability to produce a $`\rho ^0`$-meson in a close collision is of the order 1-3% for the heavy systems. It will be interesting to study these events experimentally to see how the $`\rho \rho `$ interaction affects the simple harmonic oscillator description, especially for the production of $`\rho ^0`$’s with close enough momenta.

Multiple $`e^+e^{}`$-pairs can be produced in relativistic heavy ion collisions . This work used the sudden (or Glauber) approximation and a quasiboson approximation for $`e^+e^{}`$ pairs. Thus the Hamiltonian has the form of Eq. (36), where a sum over the quantum numbers of the lepton pair has to be included. Using a QED calculation (including Coulomb corrections in the Bethe-Maximon approach) for one pair production as an input, a Poisson distribution is obtained for multiple pair production. This is quite natural, since this problem reduces to the excitation of a harmonic oscillator (the modes are labelled by the spins and momenta of the $`e^+e^{}`$ pairs). More detailed calculations have verified that the Poisson distribution holds quite well for multiple pair production, although slight deviations are expected . A characteristic dimensionless parameter for this problem is $`\mathrm{\Xi }=[(Z_1Z_2\alpha ^2)^2]/(m_eb)^2`$, where $`m_e`$ is the electron mass and $`b>1/m_e`$. Here, $`\mathrm{log}\gamma `$ should be large.

With heavy systems like Au-Au or Pb-Pb, when $`b\begin{array}{c}>\hfill \\ \hfill \end{array}1/m_e`$, $`\mathrm{\Xi }1`$. The impact parameter dependence of $`e^+e^{}`$ pair production has been studied numerically in lowest order QED. Only recently an approximate analytical formula for the total pair production probability in lowest order $`P^{(1)}`$ was found. When $`1/m_e<b<\gamma /m_e`$

$$P^{(1)}=\frac{28}{9\pi ^2}\mathrm{\Xi }(2\mathrm{ln}\gamma ^23\mathrm{ln}(m_eb))\mathrm{ln}(m_eb)$$
(45)
The $`N`$-pair production probability decreases strongly with increasing impact parameter $`b`$ (approximately like $`b^{2N}`$, for $`b>1/m_e`$). Therefore the probabilities $`P^{(1)}(b)`$ should be known accurately for an impact parameter range of $`0<b<\text{several}1/m_e`$. Since the multiple pair production happens in a single collision, with a given impact parameter, there are correlations in the momenta of the outgoing pairs. This will be briefly discussed in Sec. 3 below.

$`e^+e^{}`$ pair production is of great practical interest. It could be useful as a trigger for UPC collisions at the LHC , but also constitutes a background for the most central parts of the trigger system of ALICE. . Electron or muon pair production might also be a way to measure the $`\gamma \gamma `$-luminosity .

Muon pair production is also of interest: in close collisions the pair production probability is appreciable, although significantly smaller than one. This can be seen from Eq. (45). The muon Compton wavelength $`1/m_\mu `$ is 1.86 fm, much smaller than the nuclear radius. For $`b2R_A`$ form factor effects also need to be considered. Total cross sections are smaller by about a factor of $`(m_e/m_\mu )^2`$ $`2\times 10^5`$ (which is still large). The results of a lowest order external field QED calculation are shown in Fig. 4. Calculations with and without a monopole form factor are shown together with the approximate expression of Eq. (54) below. Details of the calculation are given in .

## 3 Some Examples of Multi-Photon Processes

The factorization of the amplitudes has some significant, experimentally accessible implications. Although the amplitudes are independent, the processes all share a common impact parameter $`\stackrel{}{b}`$. This leads to some significant correlations, due to both the magnitude and direction of $`\stackrel{}{b}`$. In this section, we consider the effect of multiple interactions on the photon spectrum and impact parameter distribution, and on the angular distributions of the final states.

### 3.1 Impact parameter dependence of the equivalent photon spectrum

The photon spectrum at impact parameter $`b>R_A`$ is

$$N(\omega ,b)=\frac{Z^2\alpha }{\pi ^2}(\frac{\omega }{\gamma })^2\left(K_1^2(x)+\frac{K_0^2(x)}{\gamma ^2}\right)$$
(46)
where $`x=\omega b/\gamma `$. As long as $`x1`$, the photon spectrum scales as $`dN/d\omega 1/b^2`$. The mean impact parameter $`\overline{b}`$ is

$$\overline{b}=\frac{d^2bbP(b)}{d^2bP(b)}$$
(47)
where the probability $`P(b)`$ of a specific reaction occurring at impact parameter $`b`$ is

$$P(b)=N(\omega ,b)\sigma _{\gamma A}.$$
(48)
For a reaction involving a single photon with energy $`\omega `$, we can integrate over the region $`b_{min}=2R_A<|b|<b_{max}=\gamma /\omega `$. Many processes can occur over a range of photon energies; in this case, an integral over $`\omega `$ is needed. However, the details of the interaction do not matter, and

$$\overline{b}=\frac{b_{max}b_{min}}{\mathrm{ln}(b_{max}/b_{min})}.$$
(49)
For single photon processes the average impact parameter $`\overline{b}`$ is large when $`b_{max}`$ is large and insensitive to the exact value of $`b_{min}`$.

However, for two or more photon interactions, such as mutual Coulomb excitation of two nuclei, $`P(b)(d^2N_\gamma /d^2b)^n1/b^{2n}`$, with $`n`$ the number of photon interactions and the latter approximation requiring that the $`x`$ for the two photons both be less than one. Then, as long as $`b_{max}b_{min}`$,

$$\overline{b}_nb_{min}\frac{2n2}{2n3}.$$
(50)
$`b_{max}`$ drops out, and the median impact parameter only depends on $`b_{min}`$, independent of the collision energy and other details of the interaction. Of course, $`b_{max}`$ depends on the photon energy and we assume that for all photon processes $`b_{max}b_{min}`$ should be fulfilled in the relevant $`\omega `$ range. For heavy nuclei like Au $`\overline{b}_22b_{min}=4R_A=26`$ fm.

For three-photon interactions, such as vector meson production accompanied by mutual Coulomb excitation of both nuclei, $`n=3`$ and $`\overline{b}_34b_{min}/3=\frac{8R_A}{3}`$. For heavy nuclei, $`\overline{b}_318`$ fm. Of course, as $`\overline{b}_n`$ drops with increasing $`n`$, the cross sections for UPCs accompanied by hadronic interactions become more sensitive to the tails of the nuclear density profile — including these tails increases $`b_{min}`$ above the canonical $`2R_A`$.

The $`\overline{b}`$ are comparable to the median impact parameters calculated in ; these detailed calculations found that the median impact parameter is almost independent of the specifics of the reaction. This impact parameter distribution is important for studying interference between vector meson production on the two nuclei; the smaller $`\overline{b}`$, the easier it is to observe the interference .

The reduction in $`\overline{b}`$ affects the photon spectrum, since the maximum photon energy scales as $`1/b`$. One way to quantify this is to consider a two-photon reaction: excitation of one nucleus to a GDR, while the second photon produces a $`\rho `$ meson. The GDR tagging has a strong effect on the $`\rho ^0`$-production: it selects small impact parameters, where the equivalent photon spectra are harder. One can calculate the equivalent photon spectrum in the case of GDR tagging. For close collisions, which matter most, the probability of GDR excitation is given in Eqs. (43) and (44). The equivalent photon spectrum is then

$$n_{GDR}(\omega )=S\frac{2Z_1^2\alpha }{\pi }(\frac{\omega }{\gamma v})^2_{x_{min}}^{\mathrm{}}\frac{dx}{x}\left(K_1(x)^2+\frac{K_0(x)^2}{\gamma ^2}\right).$$
(51)
This expression can be evaluated numerically. Due to the high values of the Lorentz boost $`\gamma `$ at the colliders, $`K_0(x)^2/\gamma ^2`$ is very small and can safely be neglected. A simple approximation  is found by setting $`x^2K_1(x)^2=1`$ for $`x<1`$ and zero otherwise; the integral over $`x`$ is then $`\frac{1}{2x_{min}^2}\frac{1}{2}`$ where $`x_{min}=\frac{\omega R_{min}}{\gamma v}<1`$. Then

$$n_{GDR}^{BJ}(\omega )S\frac{Z_1^2\alpha }{\pi }(\frac{\omega }{\gamma v})^2\left(\frac{1}{x_{min}^2}1\right).$$
(52)
The photon spectrum scales very roughly as $`n(\omega )\omega ^0`$, in contrast to the $`n(\omega )\omega ^1`$ for untagged photons.

Using Eq. (11.3.30) of Ref.  one finds a better analytic expression

$$n_{GDR}(\omega )\frac{S}{R_{min}^2}\frac{Z_1^2\alpha }{\pi }x_{min}^2(K_1(x_{min})^2K_0(x_{min})^2).$$
(53)
Figure 5 shows the EPA spectrum $`n(\omega )`$ (without tag) and the spectrum $`n_{GDR}(\omega )`$ in the approximations of Eqs. (52) and (53) respectively.

Another example of multiphoton processes is GDR tagging of $`\gamma \gamma `$-processes: if the mass of the system produced in the $`\gamma \gamma `$ interaction is low, there is a considerable influence of the tagging, if the mass is high, the influence is less. The $`b`$-dependence of (untagged) photon-photon-processes can be found from the folding of $`b`$-dependent equivalent photon spectra, see e.g. Eq. (49) of Ref. . A simple approximate formula was developed recently for the special case of $`e^+e^{}`$ pair production. The result is given in Eqs. (0.23), (0.16) and (0.17) of .

We consider here the more general case where the invariant mass of the $`\gamma \gamma `$-system is given by $`W`$ (rather than $`2m_e`$, as in Ref. ). The lower cutoff in the impact parameter $`b_{min}`$ is $`2R_A`$, (rather than the Compton wavelength $`\frac{1}{m}`$ of the electron). We introduce the parameters $`u=\frac{b}{R}`$ and $`x=\frac{2\gamma }{WR}`$ . We reconsider Eq. (0.24) of , with different integration limits: $`b_1`$ from $`R`$ to $`b`$ and $`\omega _1`$ from $`W^2R/4\gamma `$ to $`\gamma /R`$. Generalizing Eqs. (0.16) and (0.17) of , the $`b`$-dependent $`\gamma \gamma `$-luminosity is

$$\frac{d^3L_{\gamma \gamma }}{d^2bdW}=\frac{2\pi }{WR^2}\left(\frac{Z^2\alpha }{\pi ^2u}\right)^2X.$$
(54)
The probability $`P_{AA}(b)`$ of a $`\gamma \gamma `$-process in the AA collision is given in terms of this $`\gamma \gamma `$-luminosity by $`P_{AA}(b)=2\pi b\frac{dW}{W}\frac{d^3L_{\gamma \gamma }}{d^2bdW}\sigma _{\gamma \gamma }`$. In this expression $`\sigma _{\gamma \gamma }`$ is the elementary $`\gamma \gamma `$ cross-section and

$$X=\left[\mathrm{ln}(x^2/u)\right]^2\text{for }x<u<x^2$$
(55)
(for “large ”impact parameters) . For “small” impact parameters

$$X=2\mathrm{ln}(u)(2\mathrm{ln}(x/u)+1/2\mathrm{ln}(u))\text{for }2<u<x.$$
(56)
This result is obtained by a generalization of Eq. (0.24) of .

The $`b`$-dependence is strongly influenced by the value of $`W`$: for $`WR\gamma `$ (e.g., for low mass $`e^+e^{}`$ production) there is an approximate $`1/b^2`$\- dependence over a considerable range of $`b`$-values, for $`WR\gamma `$ the value of $`x`$ is quite small and the dependence is much steeper. In this case low $`b`$-values are emphasized anyhow and the effect of tagging is less dramatic than for the case of $`\rho ^0`$-production with GDR tagging, as discussed above. This was found by Klein in  for photon-photon processes with nuclear breakup. In accordance with the present considerations, the effect of tagging is more important for low $`W`$, see his plot as a function of the invariant mass of a produced muon pair $`M_{\mu \mu }`$. A comparison of the luminosity of Eq. (54) against a full calculation is shown in Fig. 6 for different invariant masses $`W`$. In view of the crude approximations made to obtain Eq. (54) the overall agreement is quite good for small invariant masses, whereas for the highest invariant masses shown in the plot some part of the cross section come from the areas beyond the range $`u<x^2`$. This figure also shows the change in the shape of the spectrum for low and high invariant masses $`W`$.

### 3.2 Angular correlations

There can also be angular correlations since the photon polarization follows the electric field vector of the nuclear fields which in turn follows $`\stackrel{}{b}`$ . For example, the decay $`\rho _0\pi ^+\pi ^{}`$ is sensitive to the photon polarization, and hence to the direction of $`\stackrel{}{b}`$. Mutual GDR excitation is another example; the neutron transverse momentum tends to follow the electric field. The amplitude $`a_{1,2}(\stackrel{}{b})`$ for mutual excitation is

$$a_{1,2}(\stackrel{}{b})=a_1(\stackrel{}{b})a_2(\stackrel{}{b})$$

where the amplitude for single GDR excitation $`a_i`$ ($`i=1,2`$) is given, e.g., by Eq. (3.1.22) of .

In intermediate energy heavy ion scattering (in contrast to the ultrarelativistic case) it is experimentally possible to measure directly the momentum transferred to the projectile. This determines the impact parameter vector $`\stackrel{}{b}`$. The decay products show a dependence on the azimuthal angle $`\varphi `$. This dependence has been observed, e.g., in the breakup reaction $`{}_{}{}^{11}\text{Be}+\text{Pb}^{10}\text{Be}+n+\text{Pb}`$ in , see their Fig. 3. The angular distribution follows a dipole distribution, see Eq. (3.1.22) in . For large values of $`\gamma `$ the $`m=\pm 1`$ components dominate and

$$a_i\mathrm{sin}\theta \mathrm{cos}\varphi ,$$
(57)
where $`\theta `$ and $`\varphi `$ are the polar and azimuthal angles of the emitted neutrons in the system of the nucleus.

For mutual GDR excitation, the direction of $`\stackrel{}{b}`$ is the same for both excitations, producing a correlation between the emission angles of $`\theta _1,\theta _2,\varphi _1,`$ and $`\varphi _2`$ (in the system of the respective emitting nucleus). The correlation in the relative azimuthal angle $`\mathrm{\Delta }\varphi =\varphi _1\varphi _2`$ of the emitted neutrons is

$$C(\mathrm{\Delta }\varphi )=1+\frac{1}{2}\mathrm{cos}(2\mathrm{\Delta }\varphi ).$$
(58)
This correlation is shown in Fig. 7. The neutrons from the GDR decays can be detected with zero degree calorimeters which are located downstream of the collision region at relativistic ion colliders. At RHIC, the ZDC’s are located $`\pm 18`$ meters downstream; they are 10 cm wide by 18.75 cm high . The dimensions are limited by the available space. Although the sensitivity is reduced near the edges due to transverse shower leakage, coverage is good for angles $`\theta <2`$ mrad . The excitation energy for a heavy nucleus GDR is approximately 10 MeV, corresponding to a neutron momentum of about 140 MeV/c.

For the RHIC neutron longitudinal momentum of 100 GeV, $`p_T=140`$ MeV/c corresponds to a deflection of 2.5 cm, well within the ZDC acceptance. With position sensitive shower-maximum detectors to the ZDC’s, neutron position resolution of order 1 mm should be achievable , giving adequate angular resolution to study neutron-neutron correlations between the ZDCs. At the LHC, the ZDCs are located $`\pm 116`$ meters from the collision regions, giving a maximum neutron deflection of 5 mm, so determining the neutron angular direction may be more difficult.

Similar angular correlations will occur for $`\rho `$ production and decay to $`\pi ^+\pi ^{}`$ accompanied by single or mutual GDR excitation. In its rest frame the $`\rho `$ decay angular distribution is given by $`|Y_{1M}(\theta ,\mathrm{\Phi })|^2`$, where $`M`$ denotes the magnetic quantum number. The population of the different $`M`$-substates depends on the dynamics of the production process.

In the heavy ion case the population is relatively simple, since the photons are almost real, with both helicities equally populated. It is the simplest assumption that in $`\rho `$ photoproduction the $`\rho `$ has the same polarization as the initial state photon - this is $`s`$-channel helicity conservation (SCHC). This has been experimentally verified for photoproduction  and for electroproduction . In our case the photons are linearly polarized, transverse to the beam. Assuming SCHC, the final $`\rho ^0`$ state is also linearly polarized, in the same direction as the equivalent photon. It is a superposition of $`M=\pm 1`$ states. The angular distribution in the rest system is

$$Y_{11}\pm Y_{11}=\sqrt{\frac{3}{\pi }}\{\begin{array}{c}\mathrm{sin}\theta \mathrm{cos}\mathrm{\Phi }\hfill \\ i\mathrm{sin}\theta \mathrm{sin}\mathrm{\Phi }\hfill \end{array},$$
(59)
quite analogous to Eq. (57) for the GDR decay. For exclusive $`\rho `$ production, the direction of the impact parameter is unknown, and one has to integrate over all impact parameters. However, if the coherent $`\rho `$-production is accompanied by a GDR excitation, the azimuthal angles of the neutron and the $`\pi ^+\pi ^{}`$-pair would be correlated, as with mutual GDR excitation. The GDR decay provides a way to measure the photon polarization. The same correlation would be observed between two $`\rho ^0`$ that are produced together .

In $`e^+e^{}`$ pair production one also expects a dependence on the azimuthal angle $`\varphi `$ which is defined by the impact parameter vector $`\stackrel{}{b}`$. For polarized (real) photons this can be seen from the analysis of, e.g., . For linearly polarized photons the electron is most likely to be emitted in the plane of the polarization . It will be of interest to do such theoretical calculations in more detail and work on this is currently in progress. E.g. in two-pair production one might expect to see such correlations, which leave their trace on the impact parameter of the collision. The formula for one-pair production in lowest order semiclassical approximation is given e.g. in . This formula contains the correlations of the impact parameter and the momenta of the outgoing $`e^+e^{}`$ pair.

We conclude this section with a general remark on multiphoton processes. The observation of two electromagnetic processes with photon energies $`\omega _1`$ and $`\omega _2`$ gives information about the range of impact parameters of the collision. One has

$$P_{12}(b)=P_1(b)P_2(b)=N(\omega _1,b)N(\omega _2,b)\sigma _\gamma (\omega _1)\sigma _\gamma (\omega _2)$$
(60)
We use a simple approximation for the $`b`$-dependent equivalent photon spectrum $`N(\omega ,b)`$:

$$N(\omega ,b)=\frac{Z_1^2\alpha }{\pi ^2b^2}\text{for }2R<b<\frac{\gamma }{\omega }$$
(61)
Thus such an interaction takes place in the b-range of $`2R<b<\frac{\gamma }{\omega _{max}}`$, where $`\omega _{max}`$ is the larger of $`\omega _1`$ and $`\omega _2`$.

## 4 Conclusions

Due to the strong fields of highly charged ions multiphoton processes are abundant in UPCs. The important modes were identified. Quite generally, the excitation of the different modes occurs independently and the corresponding amplitudes can be factorized. They do not influence each other. This is essentially due to the very high energy of the heavy ions, which is much higher than the energy $`\omega `$ of the equivalent photon. Only a small fraction of the total kinetic energy of the nuclei is lost in the UPC. Due to the coherence condition one has $`x_{max}=\frac{\omega }{E}=\frac{1}{RM_NA}=\frac{\lambda _C(A)}{R}1`$. $`M_N`$ is the nucleon mass, and $`\lambda _C(A)`$ is the Compton wavelength of the nucleus.

The main result of this paper may seem intuitively obvious. However, the experimental implications are significant. The factorization assumption was used in detailed theoretical calculations of $`\rho `$-production and GDR excitation in  (following Table 2 it is said there that “the validity of the assumption of factorization is hard to prove rigorously …”), and  and in the corresponding experimental analysis . It is important to show also how this factorization arises in the theoretical formulation of the process. Now there is a formal basis to discuss these questions theoretically.

Ultraperipheral heavy ion collisions provide a strong source of equivalent photons up to very high energies. This offers the unique possibility to study photon-hadron (nucleus) and photon-photon processes in hitherto inaccessible regions, see . Many UPC processes like GDR excitation are now recognized as useful for practical matters like luminosity measurement and impact-parameter dependent triggers. In this case, as was already shown explicitly in  one emphasizes the hard part of the spectrum, due to the selection of low $`b`$ values.

## Acknowledgments

This work was partially supported by the US Dept. of Energy, under contract No DE-AC-03076SF00098. We would like to thank Sebastian White for interesting discussions on this subject. We thank also the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support and the organizers of the program “The first three years of heavy ion physics at RHIC”. The “UPC week” of this program, which took place April 21st–27th 2003, was vital for this paper, both due to the interaction of some of the present authors and also to the stimulating atmosphere created by all the participants.