text:
# A Fast Reconstruction Algorithm for Gene Networks

## 1 Introduction

In this paper we will consider the most simple – not trivial dynamical model of a gene network, define a reverse engineering problem111In the system and control community, this problem is commonly defined as an identification problem. We prefer here to adopt the terminology used in the biological literature. and propose a fast algorithm to tackle this problem. We will deal then with continuous–time, linear, time invariant, finite dimensional system (LTI systems) at steady state, in view of the relevant biological literature (see, for example Yeung99:2002 ). Even though this oversimplified model may not be very realistic ”far” from steady state, nevertheless it is a fundamental tool for studying and gaining insight into the basic mechanism that makes this problem an hard one, thus providing a valuable in numero testbed for gene networks reconstruction algorithms. In fact, starting from the simplest and then moving toward the more and more complex is a typical scenario in the applied sciences. Consequently, in this paper we will consider only experiments regarding steady state measurements due to ”small” input constant perturbations, so to retain linearity of the model.

More precisely, we deal with the problem of reconstructing a gene network in the typical practical situation in which the number of available data is largely insufficient to uniquely determine the network. In order to try to remove this ambiguity, we will exploit some additional biologically relevant a priori assumptions such as sparseness, as it will be expalined in the following.

The proposed algorithm has the major advantage to be very fast. In fact, its complexity is $`O(2NM^2)`$ and consequently it is particularly useful for large networks (large $`N`$) and few data (small $`M`$).

Finally we will also generate artificial data for the reconstruction algorithm. By doing so, one exactly knows the ”true” gene regulation network and then algorithms performance can be evaluated. Moreover, the algorithm will be successfully applied to real data collected from experiments regarding a nine-gene subnetwork of the SOS pathway in Escherichia coli and very recently presented in Gardner2003 .

## 2 Dynamic model of expression data and problem formulation

The idea of modelling gene expression data with differential equations has been explored by a considerable numbers of authors (see for example references Chen4:1999 , Daeseleer1999 , Yeung99:2002 , Holter98:2001 , Someren2000 , Weaver1999 , Dehoon:2003  and, for a literature review, the interested reader may refer to DeJong9:2002 ). Differential equations are used to model gene interactions under the assumption that the rate of change over time of each gene expression level is a function of the expression level of some (usually a few) other genes. Such modeling assumption is based on the reaction kinetics at the biochemical level. A fully realistic model should consider a number of relevant biological issues such as, for example, internal and external noise, time delays, specific classes of nonlinearities, and should also consider the relationships between mRNA and protein concentrations since only the first one is actually measured by microarrays. Clearly, a lot of work remains to be done in the field of gene interaction modelling. Nevertheless, a very simple linear time invariant model has proved to be useful in a number of cases (as in Daeseleer1999  for rat cervical spinal cord data) even if it is clear that nonlinearity is an unavoidable issue since it reflects also the nature of biochemical interactions. However, as in Gardner2003 , we will assume the system to behave linearly around the steady states points reached in experiments. According to the work of Yeung et al. Yeung99:2002  we consider the LTI system described by the following differential equations:

$$\dot{x}_i\left(t\right)=\lambda _ix_i\left(t\right)+\underset{j=1}{\overset{N}{}}W_{ij}x_j\left(t\right)+b_i\left(t\right)+\xi _i\left(t\right)$$
(1)
for $`i=1,2,\mathrm{},N`$, where the state variables $`x_i`$’s are the concentration of mRNA measured as a difference from the equilibrium state preceding the stimulus222Clearly, these numbers can be positive, negative or zero provided that the measured concentration levels are greater, less or equal to the concentration present prior to the stimulus., the $`\lambda _i`$’s are the self–degradation rates, the $`b_i`$’s are the external stimuli (depending on the specific experiment performed), and the $`\xi _i`$’s represent (internal) noise. The elements of the matrix $`W`$ describe the type and strength of the ”influence” of the $`j`$-th gene on the $`i`$-th gene with a positive, zero or negative sign indicating activation, no interaction and repression respectively. As described in Gardner2003 , an experiment consists in applying a prescribed (i.e. known) stimulus $`b_i\left(t\right)`$ which is a persistent perturbation (ideally a step function $`\delta _1(t)`$ of ”amplitude” $`b_i`$, i.e. $`b_i\left(t\right)=b_i\delta _1(t)`$) to $`M`$ input separately and then use a microarray to measure the response at steady state. Finally, the above situation, in view of equations (1) and taking into account the measurement noise $`w`$ introduced by the microarray device, can be formally described in the usual compact form as follows:

$`\dot{x}`$ $`=Ax+B+\xi `$ (2)
$`y`$ $`=x+w`$
where the matrix $`A`$ is a $`N\times N`$ matrix which incorporates both self–degradation rates (on its main diagonal entries) and the strength of the gene–to–gene interaction (on its off diagonal entries) and the columns of the $`N\times M`$ matrix $`B`$ are the $`b_i`$’s. We will assume standard normality properties on zero mean noises.

Steady state (equilibria) solutions are given by $`Ax_e+Bu_e+\xi =0=Ay_e+Bu_e+\xi Aw`$, and if we repeat the $`M`$ measures a sufficiently large number of times, we can calculate average values and then solve $`A\overline{y}+B\overline{u}=0`$, where $`E[y_e]=\overline{y}`$, $`E[u_e]=\overline{u}`$ and $`E[\xi Aw]=0`$ so that $`E[w]=E[\xi ]=0`$.

We can now state the problem considered in Gardner2003 :

An Identification (Reverse Engineering)

Problem for Gene Networks

Given a network described by (2), initially at rest composed of $`N`$ genes, where $`M`$ independent constant inputs333The matrix $`B`$ is known and full rank, possibly diagonal. are applied and $`M`$ noisy measurements $`y^j`$ are taken at steady state. This process is repeated a ”sufficiently” large number of times so that we can average the data and obtain the data matrix  $`\overline{Y}_{N\times M}`$, whose generic element is $`\overline{y}_i^j`$ where superscripts denote experiments (which are then repeated $`M`$ times), the overbar denotes averaging and the subscripts denote individual genes.

Then, a reverse engineering problem consists in finding the ”best” matrix $`A`$ such that

$$A_{N\times N}\overline{Y}_{N\times M}+B_{N\times M}=0_{N\times M}$$
(3)
It is important to note that the matrix $`A`$ in equation (3) is unique if and only if $`MN`$, i.e. provided that the number of experiments is equal or greater than the number of genes in the network. In what follows we will assume the typical situation in which $`MN`$ so that equation (3) is underdetermined, that is, it has many solutions yielding $`A\overline{Y}_{N\times M}+B=0`$. Thus, in order to find the best choice for the matrix $`A`$, some a priori information has to be exploited incorporating some biological knowledge into the model at hand. One possibility, as discussed in a previous section, is to try to impose the additional biological constraint that usually gene networks are sparse, i.e. that generally each gene interacts with only a small percentage of all the genes in the entire genome (see for example Jeong407:2000 ). We will therefore exploit this feature and present our algorithm in the following section.

In the next paragraph, we will first address the problem of generating artificial data consistent with the above biological assumption so to be able, after presenting a new reverse engineering algorithm, to evaluate its performances on such generated artificial data and, afterwards, also on real data taken from Gardner2003  regarding the SOS pathway of Escherichia coli.

## 3 The artificial data generation

In this section we will briefly describe how one can set up an artificial problem and the data generated by this procedure. First of all we generated data assuming a linear model (as in (1)) and the algorithms were tested on such data. However, in this preliminary work we just wanted to evaluate the performances of the algorithms presented in the next paragraphs. We would like also to stress that we applied our reconstruction algorithms to artificial data, before using real data. The reason simply being that, by doing so, one exactly knows the ”true” regulation network, and the algorithm performances can be evaluated.

Consequently, in our first artificial experiments, we generate the data with the linear model described in equation (1). Each repetition of the $`i`$–th experiment consists, in this model, in applying a constant unitary input (normalized) and then in picking a sample of the trajectory generated by equations (1) at steady state. We generate then a *sparse*, *asymptotically stable* matrix $`A`$, and, as a second step, also the entries of input matrix $`B`$ for each experiment. Finally, we obtain the artificial data by considering steady state values, and by repeating this process $`M`$ times, we get the simulated data matrix $`\overline{Y}_{N\times M}`$444As previously stated, the expression level of the genes are calculated as the difference between the current expression level (the cellular concentration of the gene product)..

## 4 A reconstruction algorithm for sparse gene networks

As discussed in the previous sections, the reverse engineering problem for gene networks when the number of measurement is (much) less than the number of genes involved in the network, leads to many alternative networks among which we have to find a way to choosing the ”best” one according to a prespecified cost function.

As stated by the reverse engineering problem, once we have collected the mRNA abundance steady state measurements $`\overline{Y}_{N\times M}`$ from microarrays with $`M<N`$, we may first find one of the many solutions $`A`$ of $`A\overline{Y}+B=0`$. A standard procedure for the case of interest with $`M<N`$, is to perform a singular value decomposition on the data matrix $`\overline{Y}_{N\times M}`$ and obtain the matrix $`A`$ with the smallest $`L_2`$ norm. The procedure goes as follows. We decomposed the data matrix $`\overline{Y}_{N\times M}`$ as $`\overline{Y}_{N\times M}=U_{N\times N}S_{N\times M}W_{M\times M}^T`$, where $`U_{N\times N}`$ and $`W_{M\times M}^T`$ are unitary matrices and the entries $`\sigma _{ij}`$ of the matrix $`S_{N\times M}`$ are such that $`\sigma _{ij}=0`$ for all $`ij`$ and $`\sigma _{11}\sigma _{22}\mathrm{}\sigma _{MM}>0`$. The numbers $`\sigma _{ii}:=\sigma _i`$ are the nonnegative square roots of the eigenvalues of $`Y_{N\times M}Y_{N\times M}^T`$ known as the singular values of $`Y_{N\times M}`$. As previously stated, the SVD provide a simple way to find a solution $`A`$ to problem (3) with minimal $`L_2`$ norm, i.e. with minimal $`A_2`$, as555Note that, in this case, we have $`A_{svd}\overline{Y}_{N\times M}+B=0`$. $`A_{svd}=BWdiag_{i=1,\mathrm{},M}(\frac{1}{\sigma _i})U^T`$. All the solutions to (3) are given by $`A=A_{svd}+C`$,where $`C`$ is any $`N\times N`$ matrix satisfying $`C_{N\times M}\overline{Y}_{N\times M}=0`$, as one can easily verify by direct substitution. Equivalently, the matrix $`C`$ is such that the columns of $`C^T`$ belong to the null space of $`\overline{Y}_{N\times M}^T`$.

A possible way of ”reasonably” choosing the ”best” matrix $`A`$ without evaluating all the feasible solutions, is given next. Our aim is to take into account the sparseness assumption, that is try to incorporate information on the structure (zero pattern) of the matrix $`A`$. In particular, we will concentrate on the case in which the number $`k`$ of nonzero entries in each row of $`A`$ is not greater than the number of available measurements, i.e. $`kM`$. Consequently, we have that $`kM<N`$ with $`kN`$. The proposed procedure, hereafter formally described, consists of three steps:

The Fast Reverse Engineering Algorithm

1. Assuming $`kM<N`$ and using the full rank data matrix $`\overline{Y}_{N\times M}`$, find the minimal $`L_2`$ norm solution $`A_{svd}=BWdiag_{i=1,\mathrm{},M}(\frac{1}{\sigma _i})U^T`$ to equation $`A\overline{Y}_{N\times M}+B=0`$ using SVD of the data matrix $`Y_{N\times M}`$.
2. For each of the $`N`$ rows of $`A_{svd}`$, determine the smallest in magnitude $`NM+\eta `$ with $`\eta 1`$ entries and consider them as ”zero”.
3. Taking into account the zero pattern detected at the previous step, find row by row the unique (sparse) matrix $`A`$ solution of $`\mathrm{min}\mathrm{arg}A\overline{Y}_{N\times M}+B`$, using the pseudo-inverse operator on the selected submatrix, and we are done.

Here, the basic idea is to exploit the fact that, at the first step, the computed matrix $`A_{svd}`$ is that of minimal $`L_2`$ norm, so that it is reasonable to assume that the smallest magnitude entries correspond to zeros in the ”true” matrix $`A`$. The main advantage of this approach is its simplicity and low computational complexity as it will discussed later in the paper. It’s worth noting that Step 3 requires only that, for each of the $`N`$ rows of $`A`$, an $`(M\eta )\times M`$ matrix to be pseudo-inverted and this can be done in a very efficient computational way. The complexity of the algorithm is then $`O(2NM^2)`$ so it is particularly useful for large networks (large $`N`$) and few data (small $`M`$).

We compare here this algorithm with another one which makes use of an heuristics to find the zero pattern for each row of the matrix. Instead of applying Step 3 (of the *Fast* algorithm), for each row $`i`$, we put to zero one entry, say $`a_{ij}`$, of the matrix $`A_{svd}`$ and denote this modified matrix by $`A_{svd,ij}`$. We evaluate then the error index $`\delta _{ij}=A_{svd,ij}Y+B_2`$. The position of the first zero in the $`i`$-th row is then that corresponding to the minimal value of $`\delta _{ij}`$. We repeat this process until we have placed all the $`NM+\eta `$ zeros. More precisely, using this criterion, we find for each row of the matrix $`A_{svd}`$ the $`NM+\eta `$ zero locations and consequently we can solve uniquely the system by selecting the $`M\eta `$ nonzero elements for each row $`i`$.

Simulation results are provided and discussed in the following section.

## 5 Results on artificial data

In this section, we will present the simulation results and the performances of the proposed reconstruction algorithm. First of all, we generate artificial data as described in the previous sections, then we estimate the matrix $`A`$, applying the reconstruction algorithms with $`\eta =1`$. We measure the error by counting the percentage of sign discrepancies666The information on the sign of the entries is the most relevant information from a biological point of view since reflect the nature of the gene-to-gene interaction. $`E=100\frac{_{i=1}^N_{j=1}^Ne_{ij}}{N^2}`$ where

$$e_{ij}=\{\begin{array}{c}1\text{if sign}\left(a_{ij}\right)=\text{sign}\left(\overline{a}_{ij}\right)\hfill \\ 0\text{otherwise}\hfill \end{array}$$

In Figure 1 we show the simulation results of a network composed by $`50`$ genes. We have generated the data using the linear model. In the figure it is clear how both the algorithms we propose produce good results.

In Figure 2 we show the computing time777The algorithms described here have been fully implemented with Matlab. They ran on a Pentium IV with a clock speed of 1.8 GHz of the two algorithms. It is clear that the proposed *Fast* algorithm grows much less than *Greedy* search. If the number of the genes of the analyzed network increases (say $`200`$, $`500`$, $`1000`$ or more genes) any exhaustive or greedy search, takes to much time to be computed. This is the strength of the algorithm we propose here.

## 6 Results on real data

As an example of application of the algorithm proposed in this paper, we consider the data obtained in Gardner2003  for a nine-gene subnetwork of the SOS pathway in Escherichia coli. The collected data, i.e. the diagonal input matrix $`B_{9\times 9}`$ and the data matrix $`Y_{9\times 9}`$, is reported in the web supplement of the paper888http://www.sciencemag.org/cgi/content/full/301/5629/102/DC1.

It is apparent from the the table in Figure 3, that the estimated matrix $`A`$ performs quite well since the numbers of correct gene interaction signs (as reported by the biological literature on the SOS pathway in E. coli) is not very different from that of Gardner et al. It is very important to note that the approach used by Gardner et al. considers, for each row, all the possible positions of the $`4`$ zeros in the dynamic matrix $`A`$, thus leading to a combinatorial approach since it evaluates every feasible solution. This is viable only for small subnetworks, while our method is suitable also for very large number of genes (e.g. the Saccaromyces Cerevisiae has about 6000 genes!) being polynomial in the number of genes.

## 7 Conclusions

In this paper we have addressed the so called ”reverse engineering” problem, that is the problem of reconstructing the gene-to-gene interactions from expression data (noisy samples of the state space trajectories) in the usual case of ”few” samples and ”many” genes. The typical situation is $`10`$ samples for gene networks with $`5000`$ genes, as in Gash2000 .

We have proposed a *fast* algorithm for finding a ”reasonable” solutions among the many feasible by exploiting the a priori biological information about sparseness of the gene networks. However, this hypothesis in too generic for selecting a ”small” set of solutions, but it is of great help. The main problem is that we only know that there are ”many” zeros, but no assumptions are made about the locations of such zeros. This makes even the problem of imposing sparseness a difficult problem since it has an intrinsic combinatorial complexity. To avoid this, in the literature many approaches have been developed mainly using a greedy search strategy, but still, the computation burden is very high. In this paper we have proposed and evaluated on artificial and real data, an algorithm with a very low computational complexity (linear in the number of genes) which can be used also for genome wide measurements. From this point of view, the algorithm presented seems very promising and future work will be devoted to explore its potential on large scale real data.


format: text

text:
# Shakura-Sunyaev Disk Can Smoothly Match Advection-Dominated Accretion Flow

## 1 Introduction

The most famous model of accretion disks is the Shakura-Sunyaev disk (SSD; Shakura & Sunyaev 1973). Since SSD was constructed exactly 30 years ago, the most important breakthrough in the field of accretion disk theory has been the proposal of advection-dominated accretion flow (ADAF; Narayan & Yi 1994; Abramowicz et al. 1995). SSD and ADAF appear to be adequate to describe the outer and the inner region of black hole accretion flows, respectively, and a phenomenological SSD+ADAF model has been quite successfully applied to black hole X-ray binaries and galactic nuclei (see Narayan, Mahadevan, & Quataert 1998 for a review). In this model, however, a smooth transition from an outer SSD to an inner ADAF was only assumed, but not proved. From the physical point of view, the question remains whether (and how) an SSD can really match an ADAF.

There have been basically three classes of answers to this question. Dullemond & Turolla (1998) and Molteni, Gerardi, & Valenza (2001) gave negative answers, arguing that a smooth transition from an SSD to an ADAF was not possible. Whereas their conclusions were under some certain conditions: Dullemond & Turolla (1998) considered only the low-viscosity case (with the viscosity parameter $`\alpha 0.1`$); and Molteni et al. (2001) referred only to the plain ADAF, i.e. that with zero cooling. The second class of answers, on the other hand, is positive. A number of authors showed that the SSD-ADAF transition was realizable if an extra heat flux caused by thermal conduction was invoked either in the radial direction (Honma 1996; Manmoto & Kato 2000; Gracia et al. 2003), or in the vertical direction (Meyer & Meyer-Hofmeister 1994; Meyer, Liu, & Meyer-Hofmeister 2000). The cost of this class of answers is, in our opinion, the involvement of an additional mechanism of energy transport, and in particular, the introduction of a new unknown parameter $`\alpha _T`$ to measure thermal conduction (Manmoto & Kato 2000). The third answer was due to Takeuchi & Mineshige (1998) and Gu & Lu (2000), who suggested that the thermal instability of a radiation pressure-supported SSD could trigger the flow to jump from the SSD state to the ADAF state. This answer is also a positive one, but without involving any extra mechanism of energy transport.

In this Letter, we present our answer to the question of SSD-ADAF transition. We demonstrate that such a transition in a smooth way is possible for flows with large values of $`\alpha `$ (different from the case of Dullemond & Turolla 1998) and non-zero radiative cooling (different from the case of Molteni et al. 2001). We do not involve any extra energy transport mechanism such as conduction, and this is different from the above mentioned second class of answers, and is similar to the third answer. We discuss in some detail the relation between our results here and those of Takeuchi & Mineshige (1998) and Gu & Lu (2000).

## 2 Equations

The dynamical equations for steady state axisymmetric accretion flows we consider here are usual in the literature (e.g. Narayan, Kato & Honma 1997). That is, the continuity, radial momentum, vertical equilibrium, and angular momentum equations read

$$\dot{M}=4\pi RH\rho \upsilon ,$$
(1)
$$\upsilon \frac{d\upsilon }{dR}=\mathrm{\Omega }^2R\mathrm{\Omega }_K^2R\frac{1}{\rho }\frac{dp}{dR},$$
(2)
$$H=\frac{c_s}{\mathrm{\Omega }_K},$$
(3)
$$\frac{d\mathrm{\Omega }}{dR}=\frac{\upsilon \mathrm{\Omega }_K(\mathrm{\Omega }R^2j)}{\alpha R^2c_s^2},$$
(4)
where $`\dot{M}`$ is the constant mass accretion rate; $`R`$ is the radius; $`H`$ is the half-thickness of the flow; $`\rho `$ is the density of the accreted gas; $`\upsilon `$ is the radial velocity; $`\mathrm{\Omega }`$ is the angular velocity; $`\mathrm{\Omega }_K`$ is the Keplerian angular velocity, and $`\mathrm{\Omega }_K^2=GM/(RR_g)^2R`$ in the well known Paczyński & Wiita (1980) potential, with $`M`$ being the mass of the central black hole, and $`R_g`$ being the gravitational radius; $`p`$ is the pressure; $`c_s`$ is the isothermal sound speed of the gas, defined as $`c_s^2=p/\rho `$ ; and $`j`$ is an integration constant that represents the specific angular momentum accreted by the black hole.

As for the energy equation, we employ the form given by Narayan & Yi (1995), which is for flows composed of two-temperature plasma with bremsstrahlung and synchrotron radiation and Comptonization, i.e. a relatively complete and complex case of black hole accretion flows:

$$Q_{vis}=Q_{adv}+Q_{Cou},$$
(5)
$$Q_{Cou}=Q_{rad}.$$
(6)
These two equations are for the energy balance of the ions and of the electrons, respectively. $`Q_{vis}`$ and $`Q_{adv}`$ are the rate of viscous heating that is primarily given to the ions and the rate of advective cooling by the ions, and are expressed by, e.g. equations (5) and (6) of Gu & Lu (2000), respectively. $`Q_{Cou}`$ is the rate of energy transfer from the ions to the electrons through Coulomb collisions, and is expressed by equation (3.3) of Narayan & Yi (1995). $`Q_{rad}`$ is the rate of radiative cooling of electrons, and is calculated using a bridging formula which is valid in both optically thick and optically thin regimes,

$$Q_{rad}=8\sigma T_e^4\left(\frac{3\tau }{2}+\sqrt{3}+\frac{8\sigma T_e^4}{Q_{br}+Q_{sy}+Q_{br,c}+Q_{sy,c}}\right)^1,$$
(7)
where $`T_e`$ is the electron temperature; $`\tau `$ is the total (election scattering plus absorption) optical depth, $`\tau =\tau _{es}+\tau _{abs}=(0.34cm^2g^1)\rho H+(Q_{br}+Q_{sy}+Q_{br,c}+Q_{sy,c})/8\sigma T_e^4`$ ; and $`Q_{br}`$, $`Q_{sy}`$, $`Q_{br,c}`$, and $`Q_{sy,c}`$ are the cooling rates of bremsstrahlung radiation, synchrotron radiation, Comptonization of bremsstrahlung radiation, and Comptonization of synchrotron radiation, and are explicitly expressed by equations (3.4), (3.18), (3.23), and (3.24) of Narayan & Yi (1995), respectively.

Finally, the equation of state is needed to close the system of equations,

$$p=p_g+p_r+p_m,$$
(8)
where $`p_g=\mathrm{}\rho (T_i+T_e)`$ is the gas pressure, $`T_i`$ is the ion temperature; $`p_r=Q_{rad}(\tau +2/\sqrt{3})/4c`$ is the radiation pressure, and $`p_m=B^2/8\pi `$ is the magnetic pressure, $`B`$ is the magnetic field, and for simplicity it is usually assumed that $`\beta _mp_m/(p_g+p_m)=\mathrm{const}`$.

Note that the dynamical equations (1)-(4) are valid for both geometrically thin and thick flows (Narayan et al. 1997), and equation (7) is a convenient interpolation between the optically thin and thick limits. When the flow is extremely optically thick, equation (7) gives $`Q_{rad}=16\sigma T_e^4/3\tau `$, which is the appropriate black body limit; whereas in the optically thin limit it gives $`Q_{rad}=Q_{br}+Q_{sy}+Q_{br,c}+Q_{sy,c}`$. Thus we expect that the above set of equations can be used to verify the possible transition from an (optically thick, geometrically thin) SSD to an (optically thin, geometrically thick) ADAF.

## 3 Numerical Solutions

There are nine equations (equations -,  plus the definition $`p=\rho c_s^2`$ and the expression of $`\tau `$) for nine unknown variables $`\upsilon `$, $`\mathrm{\Omega }`$, $`c_s`$, $`H`$, $`\rho `$, $`p`$, $`T_i`$, $`T_e`$ and $`\tau `$ as functions of $`R`$, with $`M`$, $`\dot{M}`$, $`\alpha `$, $`j`$, $`\beta _m`$ and the adiabatic index $`\gamma `$ being constant flow parameters. We use the standard Runge-Kutta method to solve the set of three differential equations (2), (4), and (5) for three unknowns $`\upsilon `$, $`\mathrm{\Omega }`$, and $`c_s`$, and then obtain other variables from the remaining algebraic equations. We integrate the three differential equations from the sonic point $`R_s`$ (where the radial velocity is equal to the sound speed) both inward and outward. The derivatives $`d\upsilon /dR`$, $`d\mathrm{\Omega }/dR`$, and $`dc_s/dR`$ at $`R_s`$, which are needed to start the integration, are evaluated by applying the l’Hôpital rule. The inward, supersonic part of the solution extends to the inner boundary of the flow, i.e. to a radius $`R_{\mathrm{in}}`$ where the no-torque condition $`d\mathrm{\Omega }/dR=0`$ (i.e. $`\mathrm{\Omega }R^2=j`$) is satisfied. More important for our purpose here is the outward, subsonic part of the solution. It should be stressed that we do not specify any ad hoc outer boundary conditions. We just observe how the outward solution evolves with increasing $`R`$. On the other hand, we obtain a standard SSD solution that is calculated from a set of purely algebraic equations (e.g. Frank, King, & Raine 2002). The given flow parameters $`M`$, $`\dot{M}`$, $`\alpha `$ and $`j`$ in the SSD solution are exactly the same as those in the above solution obtained with the Runge-Kutta method. We watch if and where the two solutions can smoothly match with each other.

Figure 1 provides an example of global solution of accretion flow, i.e. the flow quantities as functions of $`R`$. The solid line represents the solution of the nine equations in $`\mathrm{\S }`$ 2, with given parameters $`\alpha =0.7`$, $`\dot{m}=0.01`$ ($`\dot{m}\dot{M}/\dot{M}_{\mathrm{Edd}}`$, with $`\dot{M}_{\mathrm{Edd}}`$ being the Eddington accretion rate), $`j=0.742(cR_g)`$, $`\gamma =1.5`$, $`\beta _m=0.5`$, and $`R_s=2.95R_g`$; and the dashed line represents the SSD solution with the same parameters $`\alpha `$, $`\dot{m}`$, and $`j`$. Note that $`R_s`$ is not another free parameter, it is the eigenvalue of the problem, and is self-consistently determined when the constant flow parameters are given. Figure 1(a) is for the radial velocity $`\upsilon `$ and the sound speed $`c_s`$. It is seen that the solid line solution is transonic, with the sonic point being marked by a filled square; while the SSD solution is subsonic everywhere, it alone cannot describe the transonic nature of black hole accretion. Figure 1(b) shows the angular momentum $`l`$ ($`=\mathrm{\Omega }R^2`$). The SSD solution follows the Keplerian distribution $`l_K`$ ($`=\mathrm{\Omega }_KR^2`$), and the solid line solution is sub-Keplerian. Figure 1(c) draws the flow’s relative thickness $`H/R`$, which is $`0.4`$ (geometrically thick) for small $`R`$, decreases as $`R`$ increases, and reaches to $`0.005`$ (geometrically thin) of the SSD solution. Figure 1(d) is for the optical depth $`\tau `$, again with increasing $`R`$ , the flow becomes from being optically thin ($`\tau 1`$) to being optically thick ($`\tau 1`$, the SSD solution). Figure 1(e) is for the ion temperature $`T_i`$ and the electron temperature $`T_e`$. The solid line solution has $`T_iT_e`$; and as $`R`$ increases, the two temperatures drop down and become identical (the SSD solution). In Figure 1(f) one sees that the advection factor $`Q_{adv}/Q_{vis}=(Q_{vis}Q_{rad})/Q_{vis}`$ is $`1`$ (advection-dominated) for small $`R`$, and decreases dramatically with increasing $`R`$, and reaches nearly zero finally (radiative cooling-dominated, the SSD solution). From these figures we conclude that the solid line solution is an ADAF solution as it has properties of transonic radial motion, sub-Keplerian rotating, and being geometrically thick, optically thin, very hot, and of course, advection-dominated; and that this ADAF solution does match with the (dashed line) SSD solution, forming together a global solution. If the transition radius $`R_{tr}`$ is defined so that $`\tau =1`$ there, then $`R_{tr}12R_g`$ in the solution of Figure 1.

## 4 Discussion

We have shown that a smooth SSD-ADAF transition is realizable for black hole accretion flows with high viscosity ($`\alpha =0.7`$ in Figure 1) and non-zero radiative cooling. Our argument is simple and naive. The equations we solve and the numerical method we use are usual. We do not introduce any extra energy transport mechanism such as thermal conduction. Perhaps the only tool somewhat special here is the bridging formula (7) expressing the radiative cooling $`Q_{rad}`$, which we need to join the optically thick regime to the optically thin regime of the flow.

Our Figure 1 looks very similar to Figure 1 of Manmoto & Kato (2000), a representative paper of the second class of answers to the question of SSD-ADAF transition as mentioned in Introduction. However, the similar results are obtained in different ways: (1) As mentioned already, Manmoto & Kato (2000) invoked radial thermal conduction and, in particular, introduced a new unknown parameter $`\alpha _T`$ to measure this extra heat transport mechanism; while we do not. (2) They used the relaxation method to solve the differential equations, and we adopt the Runge-Kutta method. In principle, the solution obtained should not be related to the numerical method, but different methods suit solving different problems. In order to have a solution for the subsonic flow between the sonic point and the outer boundary, the relaxation method requires both the sonic point condition and the outer boundary condition, and the authors using this method imposed the SSD properties (Keplerian rotating, radiation-dominated, etc.) as the outer boundary condition of ADAF, i.e. the outer boundary of ADAF had been a priori fixed to be in a state corresponding to an SSD. The Runge-Kutta method, on the other hand, requires only one boundary condition, and is adequate for the problem we address here. We do not know a priori what the outer boundary condition of ADAF ought to be, so we do not specify any; but we know for sure that black hole accretion must be transonic, so we use the Runge-Kutta method to integrate the equations starting from the sonic point, and observe how the solution behaves as $`R`$ increases. For wrong choices of the sonic point condition and the given constant flow parameters, the outward ADAF solution does not match an SSD solution. Then we try again until a correct choice is made and the ADAF solution has an asymptotic approach to an SSD solution that corresponds to the same flow parameters as those for the ADAF solution, thus we believe that a global solution containing an SSD-ADAF transition is found. The payment for using the Runge-Kutta method is that not only the variables at the sonic point, but also their derivatives there must be supplied in order to start the integration, the calculations applying the l’Hôpital rule are troublesome, though still straightforward. (3) In Manmoto & Kato (2000) the transition radius $`R_{tr}`$ was an inputted free parameter, while in our work it is determined by the constant flow parameters and is naturally calculated.

Let us now comment on the relation between our results here and the third answer to the question of SSD-ADAF transition mentioned in Introduction. Takeuchi & Mineshige (1998) suggested firstly that for large $`\alpha 1`$ the thermal instability of radiation pressure-supported SSD could trigger the SSD-ADAF transition. They made time evolutionary calculations and obtained very interesting results: because of the dominance of radiation pressure, the SSD becomes unstable at a radius ($`3.5R_g`$ in their example solution, see their Figures 2 and 3), and the outer stable parts of the flow are disturbed and evolve towards the ADAF state; this outward propagating disturbance damps and stops at a large radius ($`160R_g`$), then a transition backward to the SSD state starts from that radius, and propagates inward until $`5R_g`$ (it is larger than the original instability radius $`3.5R_g`$); finally, a two-phased flow structure really becomes persistent, i.e. the flow stays in the ADAF state inside the transition radius $`R_{tr}5R_g`$, and in the SSD state outside it. Later, Gu & Lu (2000) made a somewhat more extensive study on the mechanism of thermal instability-triggered SSD-ADAF transition, giving an $`\alpha \dot{m}`$ parameter diagram in which the region allowing the SSD-ADAF transition is clearly seen. Both Gu & Lu (2000) and the present paper are for stationary flows, so if the work of Gu & Lu (2000) corresponds to the first step of time evolutionary calculations of Takeuchi & Mineshige (1998), i.e. it proves the cause of SSD-ADAF transition and shows how to determine the original instability radius where the transition starts to occur, then our work here corresponds to the final stage of Takeuchi & Mineshige’s evolutionary sequence, i.e. it demonstrates that a stable two-phased flow can form and exist.

In order to see more clearly the relation between the original instability radius $`R_b`$ (where the SSD solution starts to break off due to the thermal instability) and the transition radius $`R_{tr}`$ (where the SSD solution matches the ADAF solution in the final two-phased structure), we show in Figure 2 how these two radii vary with $`\dot{m}`$. In this figure the solid line for $`R_{tr}`$ is obtained by numerically solving the set of equations listed in $`\mathrm{\S }`$ 2, with $`\gamma =1.5`$, $`\beta _m=0.5`$, $`\alpha =0.7`$, and $`j=0.742(cR_g)`$ (then for each value of $`\dot{m}`$ a correctly chosen value of $`R_s`$ is required in order to obtain a solution that contains an SSD-ADAF transition); and the dashed line for $`R_b`$ is drawn by applying the instability condition $`\beta p_g/(p_g+p_r)=0.4`$ in the standard SSD theory, which gives $`R_b\dot{m}^{16/21}`$ (e.g. Kato, Fukue & Mineshige 1998). The solution of Figure 1 corresponding to $`\dot{m}=0.01`$ is marked by filled squares, which has $`R_{tr}12R_g`$ and $`R_b4.5R_g`$. It is seen that $`R_{tr}`$ is insensitive to $`\dot{m}`$, and is always larger than $`R_b`$; and the larger $`\dot{m}`$ is, the closer the two radii are. It is also clear from Figure 2 that $`R_{tr}`$ is close to the central black hole, and this is because, according to Takeuchi & Mineshige (1998) and Gu & Lu (2000), the SSD-ADAF transition is caused by the thermal instability in the radiation pressure-supported region, i.e. in the very inner part of SSD. These results about $`R_{tr}`$ are distinctive from those in other transition mechanisms. For example, in the SSD-ADAF transition model involving radial thermal conduction, $`R_{tr}`$ appears as an inputted free parameter and has a very wide range, i.e. from a few to $`10^4R_g`$ (Manmoto & Kato 2000). It is worth studying further whether the SSD-ADAF transition radius ought to be close to the central black hole or it could be far away from the hole.

This work was supported by the National Science Foundation of China under Grant No. 10233030.


text:
# Relaxation effect and radiative corrections in many-electron atoms

## Abstract

We illuminate the importance of a self-consistent many-body treatment in calculations of vacuum polarization corrections to the energies of atomic orbitals in many-electron atoms. Including vacuum polarization in the atomic Hamiltonian causes a substantial re-adjustment (relaxation) of the electrostatic self-consistent field. The induced change in the electrostatic energies is substantial for states with the orbital angular momentum $`l>0`$. For such orbitals, the relaxation mechanism determines the sign and even the order of magnitude of the total vacuum polarization correction. This relaxation mechanism is illustrated with numerical results for the Cs atom.

preprint: UNR Jan 2004-

Compared to hydrogenic one-electron systems, calculation of radiative corrections for many-electron atoms brings in an additional layer of complexity: a strong Coulomb repulsion between the electrons. The problem is especially challenging for neutral many-electron atoms, where the interaction of an outer-shell electron with other electrons is comparable to its interaction with the nucleus. At the same time, a reliable calculation of radiative corrections for a heavy neutral system is required in evaluation of the parity non-conserving (PNC) amplitude in the 55-electron 133Cs atom. Here it has been only recently realized that the sizes of radiative corrections Sushkov (2001); Johnson et al. (2001); Kuchiev and Flambaum (2002); Milstein et al. (2002); Sapirstein et al. (2003) are comparable to the experimental error bar Wood et al. (1997) of 0.35% and, together with the Breit correction Derevianko (2000), dramatically affect agreement (or disagreement Bennett and Wieman (1999)) with the Standard Model of elementary particles.

A systematic approach to the problem of radiative corrections in strongly correlated systems is to start from a Furry representation based on a self-consistent electronic potential Sapirstein (1998). This potential takes into account the fact that an electron moves in an average field created by both the nucleus and other electrons. Based on this idea, a program of calculating radiative corrections to PNC amplitudes have been put forth by Sapirstein et al. (2003). Kuchiev and Flambaum (2002) and Milstein et al. (2002) pursue a more qualitative approach using an independent-electron approximation. We believe that the question of an interplay between correlations and radiative corrections is yet to be addressed. While here we do not compute the PNC corrections, we illuminate a situation where disregarding correlations would lead to a substantial error in determining radiative correction: a radiative correction changes sign and even the order of magnitude when the presence of other electrons is accounted for.

In particular, we consider vacuum polarization (VP) corrections to energies of atomic states. To the leading order in $`\alpha Z`$ the VP may be accounted for by introducing the Uehling potential $`U_{\text{VP}}(r)`$ into the atomic Hamiltonian. This potential is attractive, and for a hydrogen-like ion the resulting VP corrections to the energies are always negative. For a complex atom, we find by contrast that, for orbitals with $`l>0`$, the total correction is positive. Briefly, the reason for such a counterintuitive effect is due to a readjustment of atomic orbitals when the $`U_{\text{VP}}(r)`$ potential is added to the self-consistent Dirac-Hartree-Fock (DHF) equations. The innermost $`1s`$ orbitals are “pulled in” by the short-ranged VP potential, leading to a decrease of the effective nuclear charge seen by the outer orbitals and thus to an increase of the electrostatic energy of these orbitals. Since for orbitals with $`l>0`$, overlap with $`U_{\text{VP}}(r)`$ and thus the lowest order correction are small, the resulting indirect “relaxation” contribution dominates the total VP correction to the energies. In the following we will present numerical results supporting this relaxation mechanism. Atomic units ($`\mathrm{}=|e|=m_e1`$) are used throughout.

Because of our interest in PNC in Cs, below we illustrate the relaxation effect with numerical results for this atom; however, the relaxation mechanism is also applicable in the cases of other many-electron atoms. We also notice that the relaxation mechanism described here is similar to that observed in calculations of the Breit corrections Lindroth et al. (1989); Derevianko (2002).

The conventional many-electron Hamiltonian may be represented as

$$H=\underset{i}{}h_0(i)+\frac{1}{2}\underset{ij}{}\frac{1}{r_{ij}},$$
(1)
where the single-particle Dirac Hamiltonian is

$$h_0(i)=c(𝜶_i𝒑_i)+\beta _ic^2+V_{\text{nuc}}(r_i).$$
(2)
The nuclear potential $`V_{\text{nuc}}(r)`$ is obtained from the nuclear charge distribution $`\rho _{\text{nuc}}(r)`$; which is we approximate by the Fermi distribution

$$\rho _{\text{nuc}}(r)=\frac{\rho _0}{1+\mathrm{exp}[(rc)/a]},$$
(3)
where $`\rho _0`$ is the normalization constant, $`c`$ and $`a`$ are the nuclear parameters. In our the numerical example for 133Cs, we use $`c=5.6748`$ fm and $`a=0.52`$ fm.

A common starting point for describing a multi-electron atom is the self-consistent field method. Here the many-body wave-function is approximated by a Slater determinant constructed from single-particle orbitals (bi-spinors) $`u_k(𝒓)`$. The orbitals are obtained by solving self-consistently the eigenvalue equations

$$\left(h_0+U_{\text{DHF}}\right)u_k(𝒓)=\epsilon _ku_k(𝒓),$$
(4)
where $`U_{\text{DHF}}`$ is the traditional DHF potential which depends on the orbitals occupied in the Slater determinant. The DHF energies for the core and several valence orbitals of Cs are listed in Table 1.

The polarization of the vacuum by the nucleus modifies the nuclear electric field seen by the electrons. To the leading order in $`\alpha Z`$, the VP may be conveniently described with the Uehling potential, which for a point-like nucleus of charge $`Z`$ reads

$$U_{\text{VP}}^{\text{p.c.}}\left(r\right)=\frac{2}{3\pi }\frac{\alpha Z}{r}_1^{\mathrm{}}𝑑t\sqrt{t^21}\left(\frac{1}{t^2}+\frac{1}{2t^4}\right)\mathrm{exp}\left[\frac{2r}{\alpha }t\right].$$
(5)
This potential must be folded with the nuclear charge distribution,

$$U_{\text{VP}}\left(r\right)=𝑑𝒓^{\mathbf{}}\rho _{\text{nuc}}(|𝒓𝒓^{}|)U_{\text{VP}}^{\text{p.c.}}\left(r^{}\right).$$

We approximated $`\rho _{\text{nuc}}(r)`$ with the Fermi distribution, Eq. (3). In the numerical evaluation of the extended-nucleus Uehling potential, we employed the routine from Ref. Hnizdo (1994). The Uehling potential $`U_{\text{VP}}\left(r\right)`$ generated by the Cs nucleus is shown in Fig. 1. Notice that the actual range of this potential is a few nuclear radii (instead of Compton wavelength $`\mathrm{\lambda ̄}_e384\text{fm}`$), because the potential for a point-like charge, Eq.(5), diverges logarithmically as $`r0`$; therefore the folded potential $`U_{\text{VP}}`$ is dominated by the contributions accumulated inside the nucleus.

How does one compute the VP corrections $`\delta \epsilon _k`$ to the energies of the atomic orbitals? Below we consider two possibilities: (i) lowest-order perturbative treatment,

$$\delta \epsilon _k^{(1)}=u_k|U_{\text{VP}}|u_k,$$
(6)
and (ii) the self-consistent approach. Indeed, as in Ref. Johnson et al. (2001), the VP potential may be introduced into the DHF equations,

$$\left(h_0+U_{\text{VP}}+U_{\text{DHF}}^{}\right)u_k^{}(𝒓)=\epsilon _k^{}u_k^{}(𝒓),$$
(7)
and a set of new energies $`\epsilon _k^{}`$ and orbitals $`u_k^{}(𝒓)`$ is obtained. Notice that the DHF potential is modified as well, since it depends on the new set of the occupied orbitals $`u_k^{}(𝒓)`$. The correlated VP correction to the energy of the orbital $`k`$ is simply

$$\delta \epsilon _k^{\text{DHF}}=\epsilon _k^{}\epsilon _k.$$
(8)
Additionally, we carried out an independent correlated calculation in the framework of the linearized Coupled DHF approximation Dalgarno (1966), which is equivalent to the random-phase approximation (RPA). This approximation describes a linear response of the atomic orbitals to the perturbing interaction, i.e. the VP potential. Numerical values obtained from the linearized coupled DHF calculations were in close agreement with the full DHF results.

The numerical results of our calculations are presented in Table 1. While analyzing this Table, we observe that the lowest order corrections, $`\delta \epsilon _k^{(1)}`$, are always negative, reflecting the fact that the Uehling potential is attractive (see Fig. 1). Owing to the short-ranged nature of VP, and the fact that only the $`s`$-orbitals have a significant overlap with the nucleus, the corrections to the energies of $`l=0`$ orbitals are much larger than those for $`l>0`$ orbitals. As to the correlated corrections, they differ quite substantially from the lowest order-corrections. A comparison of Eq. (7) and Eq. (4) reveals the origin of this discrepancy: the perturbation, in addition to the Uehling potential, contains a difference between the two DHF potentials

$$\delta U=U_{\text{VP}}+\left(U_{\text{DHF}}^{}U_{\text{DHF}}\right).$$
(9)
For orbitals with $`l>0`$, where the first term above is small, the modification of the DHF potential contributes significantly to the VP energy corrections.

The modification of the DHF potential induced by the vacuum polarization is clearly a many-body effect, not present in hydrogen-like system. Such an effect has been explored before, for example in calculations of the Breit corrections Lindroth et al. (1989); Derevianko (2002), and it is commonly referred to as a relaxation mechanism. Let us illustrate this relaxation mechanism. Denoting the correction to the occupied orbital wave functions as $`\chi _a(𝒓)=u_a^{}(𝒓)u_a(𝒓)`$, we write

$`\left(U_{\text{DHF}}^{}U_{\text{DHF}}\right)(𝒓)`$
$`{\displaystyle \underset{a}{}}{\displaystyle \chi _a^{}(𝒓^{\mathbf{}})\frac{1}{|rr^{}|}u_a(𝒓^{\mathbf{}})𝑑𝒓^{\mathbf{}}}+`$
$`{\displaystyle \underset{a}{}}{\displaystyle u_a^{}(𝒓^{\mathbf{}})\frac{1}{|rr^{}|}\chi _a(𝒓^{\mathbf{}})𝑑𝒓^{\mathbf{}}}\text{exchange},`$
where we discarded contributions non-linear in $`\chi _a(𝒓)`$, and “exchange” denotes non-local part of the perturbation. The first two (direct) terms can be interpreted as an electrostatic potential produced by a perturbation $`\delta \rho _{\text{el}}(r)`$ in the radial electronic density

$$\rho _{\text{el}}(r)=\frac{1}{4\pi r^2}\underset{a}{}u_a^{}(𝒓)u_a(𝒓).$$

We plot both the electronic density $`\rho _{\text{el}}(r)`$ and the VP-induced perturbation $`\delta \rho _{\text{el}}(r)`$ in Fig. 2. The minima of $`\rho _{\text{el}}(r)`$ correspond to positions of the electronic shells, marked on the plot by their values of principal quantum number $`n`$.

The figure 2 may be interpreted in the following way: the $`s`$ orbitals are “pulled in” by the attractive Uehling potential closer to the nucleus. As a result, screening of the nuclear charge by the inner orbitals becomes more efficient. For example, the modification of the effective charge felt by the $`n=2`$ electrons is simply the area under the $`\delta \rho _{\text{el}}(r)`$ curve, accumulated between $`r=0`$ and the radius of the shell ($`r0.08a_0`$); from Fig. 2 it is clear that the induced modification of the effective charge for the $`n=2`$ shell has a negative sign. Such an enhanced screening leads to a reduced attraction of the electrons by the nucleus and to the increase in the energy of the outer electrons. From Table 1, we see that this indirect relaxation contribution to the energy may be well comparable to the direct VP correction, $`\delta \epsilon _k^{(1)}`$. While for $`l=0`$ orbitals the direct correction gives a reasonable estimate, for all orbitals with $`l>0`$, the neglect of the relaxation would lead to even qualitatively incorrect result. Moreover, the higher the orbital angular momentum, the smaller is the direct correction, and the more important is the relaxation mechanism. For example, for $`4d`$ orbitals the VP correction in the lowest order is four orders of magnitude smaller than the correlated result.

To summarize, here we illuminated the importance of the self-consistent many-body treatment in calculations of vacuum polarization corrections. Including the VP Uehling potential into the atomic Hamiltonian causes re-adjustment (relaxation) of the electrostatic self-consistent field. The induced change in the electrostatic energies is substantial for states with the orbital angular momentum $`l>0`$. As illustrated in our numerical results for Cs, the relaxation mechanism determines the sign and even the order of magnitude of the total VP correction for orbitals with $`l>0`$.

###### Acknowledgements.

The work of A. D. and B. R. was supported in part National Science Foundation Grant No. PHY-00-99419, the work of W.R.J. was supported in part National Science Foundation Grant No. PHY-01-39928.


text:
# The last Gamma Ray Burst in our Galaxy? On the observed cosmic ray excess at particle energy 10¹⁸ eV

## 1 Introduction

For some time now, the detection of an excess in $`10^{18}`$ eV cosmic rays from the general direction of the Galactic Center (GC) by the air shower array AGASA Hayashida et al. (1999) has been a special riddle in Galactic cosmic ray research. The angular region of the excess observed by AGASA and SUGAR is illustrated in Fig. 1.

Gamma Ray Bursts (GRB) have long been argued to produce high energy cosmic rays. In this paper we will show what the observational signature should be and try to demonstrate that the AGASA excess can be attributed to the last one, or last few GRB in our inner Galaxy.

There is a significant excess around $`10^{18}`$ eV, seen by two different experiments, AGASA and SUGAR, of cosmic rays coming from the general direction of the GC region, that cannot be attributed to the expected gradient of cosmic ray sources and distribution. This excess has been supported at the time by the air fluorescence detector Fly’s Eye Bird (1999). However the sky coverage of AGASA does not include the GC, while that of SUGAR does Clay (2001). The SUGAR data suggest a point source to within their spatial resolution, while AGASA shows an extended source.

The excess starts to be significant around $`3\times 10^{17}`$ eV, peaks near $`10^{18}`$ eV, and cuts off sharply at about $`3\times 10^{18}`$ eV, see Teshima et al. (2001).

The flux of the excess particles can be turned into a luminosity of particles beyond $`10^{18}`$ eV of about $`4\times 10^{30}`$ erg/s. Since AGASA cannot observe the entire region, this inferred luminosity must be a lower limit, with the true luminosity possibly being a factor of 3 - 10 larger.

As suggested by the AGASA Collab. these events may arise from neutrons Hayashida et al. (1999); the possibility that they are due to photons was discussed and discarded by Bellido (2002) on the basis of SUGAR data. Firstly, neutrons are not deflected due to the Galactic magnetic fields. Secondly, a peak at exactly $`10^{18}`$ eV corresponds to the distance to the GC region (in distance from here between 5 and 11 kpc, using the region of highest star formation rate of about 3 kpc around the GC Güsten & Mezger (1983)), folding the neutron decay with an injected power law spectrum. Only with neutrons would there be a lower limit in energy, above which there can be significant flux. With photons there is no such threshold.

In this case the original proton must have had an energy about three times that observed, and so we require proton energies of at least up to $`6\times 10^{18}`$ eV.

There are three main mechanisms and respective sites to accelerate particles in the Galaxy: supernova explosions either in the interstellar medium, in young and hot star bubbles, or in massive star winds. In any of the three cases such an energy per nucleon cannot be reached for any reasonable parameter of shock velocity and/or magnetic field Lagage & Cesarsky (1983); Jokipii (1987); Biermann & Cassinelli (1993); Biermann & Strom (1993); Biermann (1994).

The only way to accelerate particles to such an energy per nucleon in a normal galaxy as ours is relativistic shocks. Such relativistic shocks are produced in GRBs Vietri (1998); Piran 1999 . GRBs are believed to occur in every galaxy in some small fraction of all supernovae, of order $`10^4`$, arising from the final evolutionary stages of very massive stars, either in a binary system or as a single star. Therefore we propose to ask: What are GRBs expected to produce in terms of energetic particles?

As shown in Rachen & Mészáros (1998), because of adiabatic losses, the highest energy particles that emerge from a GRB are mostly neutrons; protons are captive in the magnetic field and so suffer extensive adiabatic losses on the way out. These neutrons will decay after some distance, turning into protons, which are then caught by the magnetic field in the Galaxy, and rumble around with a rather short residence time scale. There is a small, but finite probability that they will produce a neutron again in interactions with the interstellar medium. These secondary neutrons then could travel undeflected to us to be observed.

We will try to follow the neutrons originally ejected from a GRB. Fig. 2 shows our concept in three key time steps.

## 2 The last GRB event in our Galaxy

In order to estimate the remaining traces of any activity of cosmic rays ejected and/or produced by GRBs, we discuss below various probabilities, the residence time for the protons resulting from the decay of the neutrons ejected, and present spectral and spatial arguments.

We will somewhat arbitrarily adopt here a GRB rate of 1 per $`10^6`$ years for the inner Galaxy (considering Piran (2001), this may seem pessimistic, but considering Frail (2001), even optimistic).

The residence time scale for energetic protons in the Solar neighborhood is around, Gaisser (1990); Ferrando (1994), $`\tau _{res}=\mathrm{\hspace{0.17em}2}\times 10^7\mathrm{yrs}(E/10^9\mathrm{eV})^{1/3}`$ and so is about $`2\times 10^4`$ yrs at $`10^{18}`$ eV . Obviously, near the GC this time scale could be different. Since the magnetic field strength is higher, and the path to the outer parts of the Galactic halo is also larger it is likely that the time is longer, and so we adopt $`10^5`$ years. After an injection event the number of cosmic rays at some energy will decay in a 3D-slab diffusion approximation with $`t^{3/2}`$ for times up to the diffusive reservoir time, and with $`t^{5/2}`$ thereafter. For the nominal injection time scale of every $`10^6`$ yrs, the expected diminution factor is about 1/300 at again the typical time to the next GRB. The decay length of a neutron at $`10^{18}`$ eV exceeds the thickness of the GC thin disk-like region of high density and high star formation rate. Using a 3 kpc radius of the region of high interstellar medium density and high star formation rate in the Galactic disk about the GC and a decay length of 10 kpc this leads to a suppression factor of about 1/3. Another important aspect of the model is the probability of a proton to produce again an energetic secondary neutron in a collision with a nucleus of the interstellar medium in the inner Galaxy. Using conservatively the typical grammage inferred from the boron to carbon ratio in cosmic rays, with the Solar neighborhood again as a reference, this probability is of the order of a few percent, estimated to 1/20 (simulated with cross sections from PDG (2002)). Next we compare in the following the flux of particles above $`10^{18}`$ eV with the flux at all cosmic ray energies combined. Assuming the injection spectrum to be around $`E^{2.2}`$ at first produces a diminution of the flux beyond $`10^{18}`$ eV by a factor of about 1/100 Rachen & Mészáros (1998); Bednarz & Ostrowski (1998). Photons are produced in decays of neutral pions, but the typical energy of secondary pions is lower than that of energetic secondary protons and neutrons. Thus the final photon energy is lower and so may become submerged in the higher cosmic ray flux at lower energies. In Fig. 3 we show a simulation for a power law injection spectrum with an arbitrary exponential cutoff in energy at $`E_{\mathrm{cut}}=10^{18}`$ eV: the differential flux is $`E^{2.2}\mathrm{for}E<E_{\mathrm{cut}}`$ and $`E^{2.2}\mathrm{exp}\left[(EE_{\mathrm{cut}})/E_{\mathrm{cut}}\right]\mathrm{for}EE_{\mathrm{cut}}`$. Down to an energy of about a factor ten below the cutoff energy the effect of the cutoff is clearly visible in the ratio of the various fluxes, and below that energy, the ratio approaches a constant. What is the probability that a GRB jet points more or less along the Galactic plane? As seen from the GC region the Galactic plane is quite thick, and so maybe half the sky has an appreciable column density – see, e.g., Zylka et al. (1990) – leading to a factor of 1/2. And, finally, what is the total CR production by a “typical” GRB? We adopt here the generous number of $`10^{51}`$ erg in pure cosmic rays from Piran 1999 ; Pugliese et al. (2000); Berezinsky (2002).

Therefore the expected flux today, from the last GRB occurring $`10^6`$ yrs ago is now $`10^{31}`$ erg/s. This is just above what is observed, and so allowing for uncertainties mainly due to the limited sky coverage of AGASA, a very plausible estimate to explain the data.

The observed spectrum would be completely dominated by the two step propagation of the secondary neutrons in such a picture. Therefore the spectrum is the folding of the production spectrum, with the decay probability inside the available space, so a hump from the minimum distance to get any neutrons, see Fig.3, to the maximum energy possible from GRB productions.

Last we wish to ask whether it could be that we observe the effect of a few GRBs. The most important parameters are 1) the time since the GRB event, and 2) the orientation with respect to the Galactic plane. So adopting a picket-fence model for the timing of subsequent GRBs, we note that the flux from a GRB that was twice as old is decreased by a factor of about 5 with respect to the younger one, and a GRB at three times the age shows a flux that is a factor of 16 down.

We conclude that the observed distribution is rather likely to be the result of several GRB events in the GC region.

## 3 Predictions and tests

Large numbers of photons, electrons and neutrinos are produced in the collisions that give rise to the second generation neutrons in such a picture. However, their mean energy is small relative to that of the neutrons. The fluxes of photons, electrons and neutrinos are expected to follow an energy spectrum similar to that of the interaction protons. Furthermore, photons and neutrinos map out the spatial and angular distribution of the proton interactions. It is convenient to express the secondary particle spectra in terms of the primary proton spectrum by multiplying it with appropriate reduction factors. Simulations using the Monte Carlo event generator SIBYLL 2.1 Fletcher et al. (1994); Engel et al. (1999) predict the following reduction factors: for secondary protons and antiprotons 0.27, for neutrons and antineutrons 0.09, for photons 0.11, for electron-positron pairs 0.05, and for neutrinos (all flavors) 0.13. All these numbers are normalized to a primary proton spectrum, using a power law of $`E^{2.2}`$, and the energy range $`10^{17}`$ \- $`10^{18}`$ eV. These numbers are the ratio of the fluxes far below the upper energy cutoff. Observable is the ratio of the uncharged components, e.g., the ratio of neutrons to gamma-rays, which is here close to 1; however, as Fig. 3 shows, near to the upper cutoff the photons drop off earlier than the neutrons. The curves can be shifted in energy for any other assumption of maximum energy, since they will look the same relative to maximum energy (here $`10^{18}`$ eV with a following exponential cutoff).

To see an appreciable flux of neutrons peaking near $`10^{18}`$ eV with a visible extension to $`2\times 10^{18}`$ eV requires that the primary proton/neutron flux extends to at least about $`6\times 10^{18}`$ eV. A measurement of the ratio of neutrons to photons, with a simultaneous determination of the injected powerlaw slope, would then allow to estimate the real cutoff energy of the injected proton/neutrons.

In Fig. 4 we show the 95% confidence level for the expected flux at an AGASA-like experiment at Earth, using 1000 random realizations of GRBs in the inner Galaxy. The horizontal line shows roughly the average value actually observed for the excess signal by AGASA. The density distributions are from Lepine (2003), and the magnetic field distributions from Stanev (1995); Beck et al. (1996); Harari et al. (1999); Beck (2001); Medina Tanco (2001). In these specific simulations the cross-sections were taken from PDG (2002).

The spatial distribution of the neutrons should mimic roughly a folding of the decay distance at some energy of the presumed original neutron with the matter distribution. Thus it should be similar to the gamma ray emission and the FIR emission along the Galactic plane. We note that the estimated residence time at $`10^{18}`$ eV energies corresponds to a real travel distance via diffusion of only about 2 to 3 kpc, while the propagation and interaction distance (i.e. the length along the meandering path) is 30 kpc.

It is interesting to consider the time evolution of such a neutron flux: During a first phase protons resulting from neutron decay interact to produce new neutrons, and so steadily increase the observable neutron flux. Second, we begin to use up the reservoir of protons, so the flux of neutrons decreases slowly, with $`t^{3/2}`$. And third, the flux begins to decay with $`t^{5/2}`$, as protons leak out from the probable interaction volume.

The spatial distribution strongly depends on the specific orientation and the exact location of the original double jet of the GRB - if that is the best model. Of course, it is very unlikely, that the GRB jet would point exactly at Earth - that might be damaging to us, since this is the same power in a few seconds as the Sun gives in minutes, and in a less benign form. Therefore, if we see several GRB events and their consequences superimposed on each other, then the discrepancy between AGASA and SUGAR in spatial distribution might be partially explainable.

## 4 Conclusions

We have shown that it is plausible that the observed AGASA excess of events near $`10^{18}`$ eV energies coming to us from the GC region is due to the last few GRB events in the Galaxy. We predict a corresponding flux in photons and neutrinos.

In fact, if the predicted details can be confirmed, we will have established i) that GRB cosmic ray signature can be detected, ii) the cosmic ray production of GRBs to be of order $`10^{51}`$ erg, iii) that their particle energy extends to at least $`6\times 10^{18}`$ eV, iv) that the maximum particle energy can be estimated with a measurement of both neutrons and photons, as well as the slope of the injection spectrum, and v) that their contribution to the overall energetics of Galactic cosmic rays is minor. To check this will be a major contribution of the Pierre Auger Observatory whose southern part is ideally located to observe the GC region Blümer et al. (2003). The combination of fluorescence and surface detectors of this experiment allow measurements in the energy region from several $`10^{17}`$ eV to the highest energies.

## 5 Acknowledgements

The authors appreciate comments from A. Watson and T. Kellmann. PLB is mainly supported through the AUGER theory and membership grant 05 CU1ERA/3 through DESY/BMBF (Germany); GMT is also supported by FAPESP and CNPq.


text:
# Nonclassical imaging for a quantum search of trapped ions

(June 4, 2004)

## Abstract

We discuss a simple search problem which can be pursued with different methods, either on a classical or on a quantum basis. The system is represented by a chain of trapped ions. The ion to be searched is a member of that chain, consists, however, of an isotopic species different to the others. It is shown that the classical imaging may lead as fast to the final result as the quantum imaging. However, for the discussed case the quantum method gives more flexibility and higher precision when the number of ions considered in the chain is increasing. In addition, interferences are observable even when the distances between the ions is smaller than half a wavelength of the incident light.

photon correlation, quantum search, linear ion trap

Quantum search algorithms grover  enable us to determine an object from a black box with $`N`$ elements by a number of measurements which is of the order of $`\sqrt{N}`$. A quantum search thus provides a polynomial speed-up in comparison with any known classical algorithm. Several methods have been proposed for implementing a quantum search Chuang98 . In this letter we discuss a very simple example where the search can be performed in the same system either on a classical or a quantum basis. Such a situation is useful for illustrating the particular aspects of a quantum search with respect to a classical search. We recall that fluorescence imaging, i.e., measurement of the mean intensity of radiation scattered by atoms upon excitation, provides information about the density profile. In fact, this is a common method used to image, for example, trapped ions or a Bose-Einstein condensate Walther92 ; Blatt98 ; Ketterle00 . We could, however, do more than just image the sources of fluorescence, e.g. by determining the spectrum of the scattered light. In this case, in addition to the position the motional state of the atoms is accessible. Instead of imaging the particles individually on a detector, one could also observe the fluorescence in the far field or Fourier plane of a lens eichmann . In this case, the intensity distribution corresponds to the interference pattern described by the first order correlation function of the fluorescence light. It results from the contribution of all particles at once, and their relative position can be deduced from the interference pattern. The intensity distribution thus contains more information than the direct imaging of the ions grover2 . A further step in sophistication would be to observe the scattered light using two detectors and to measure the second order correlation function. This corresponds to a nonclassical imaging technique where the associated spatial distribution is determined by the different paths the photons can take when reaching the two detectors; the corresponding interference pattern again relies on the contribution of all scatterers simultaneously but for certain excitation angles is purely due to quantum interferences skornia01 ; agarwal02 . In this type of experiment only the coincidence events at the two detectors are recorded. As the paths of the photons contributing to these coincidences change with the relative position of the detectors and/or the scatterers, the observed interference patterns again allow the positions of the individual particles to be retraced. However, owing to the larger variability of the parameters involved, a much richer interference structure is obtained. To illustrate this in more detail we consider a particular example, viz. a linear chain of ions of the same atomic species. In order to set up a search procedure we assume that one of the ions belongs to an isotopic species different to the others. By calculating the spatial photon - photon correlations we then show how the second order correlations can be used to reveal information on the position of the off-resonant, non-radiating isotope in the chain. We will also demonstrate why this information is more extensive than that obtained from the usual fluorescence imaging techniques or first order correlation function. In particular, it will be shown that the second order correlation function is able to provide data in parameter regions where the first order correlation function is not able to provide any information.

Let us consider a chain of ions with an energy level scheme as shown in Fig. 1a. The transition $`|g|e`$ is used for exciting the system, whereas the transition $`|e|f`$ serves for fluorescence detection. The excitation could be performed, by for example, a short pulse of radiation. The most basic search protocol would consist in resonantly exciting each ion one by one on the $`|g|e`$ transition using a focused laser beam and observing the fluorescence scattered by the ion on the $`|e|f`$ transition. With the isotope being off-resonant and remaining dark upon excitation, it will take at most $`N1`$ steps to localize the isotope, i.e. in the mean $`N/2`$ trials. One could next observe the fluorescence intensity in the far field or in the Fourier plane of a lens, i.e. without imaging the ions. The corresponding interference pattern - produced by the simultaneous superposition of the electromagnetic fields of all scatterers - can be used to extract information about the position of the isotope. To be explicit, let us calculate the far-field intensity produced by a chain of ions at some distance $`r`$. The positive frequency part of the field amplitude can be written as (see agarwalbook )

$`\stackrel{}{E}^{(+)}(\stackrel{}{r}_1,t)`$ $`=`$ $`\stackrel{}{E}_0^{(+)}(\stackrel{}{r}_1,t)`$
$``$ $`{\displaystyle \frac{k^2e^{ikr}}{r}}{\displaystyle \underset{j}{}}\mathrm{e}^{ik\stackrel{}{n}_1\stackrel{}{R}_j}\stackrel{}{n}_1\times (\stackrel{}{n}_1\times \stackrel{}{d}_{fe})A_{fe}^{(j)},`$
where $`\stackrel{}{n}_1`$ stands for a unit vector in the direction of the detector position, $`\stackrel{}{r}_1=\stackrel{}{n}_1r`$, $`k=\omega _0/c`$. $`\omega _0`$, $`\stackrel{}{d}_{fe}`$ and $`A_{fe}^{(j)}=|f_{jj}e|`$ are, respectively, the transition frequency, dipole matrix element and the dipole operator of ion $`j`$ for the transition $`|e|f`$ and the sum is over all ion positions $`\stackrel{}{R}_j`$ in the chain. The far-field intensity is thus determined by the first order correlation function $`G^{(1)}(\stackrel{}{r}_1)=A^{}(\stackrel{}{r}_1)A(\stackrel{}{r}_1)`$, where $`A(\stackrel{}{r}_1)=_jA_{fe}^{(j)}\mathrm{e}^{ik\stackrel{}{n}_1\stackrel{}{R}_j}`$. In the case of uncorrelated ions ($`A_{ef}^{(j)}A_{fe}^{(i)}=A_{ef}^{(j)}A_{fe}^{(i)}`$), $`G^{(1)}(\stackrel{}{r}_1)`$ can be simplified to agarwalbook

$$G^{(1)}(\stackrel{}{r}_1)=\underset{j}{}A_{ee}^{(j)}+\underset{ji}{}A_{ef}^{(j)}A_{fe}^{(i)}\mathrm{e}^{ik\stackrel{}{n}_1(\stackrel{}{R}_i\stackrel{}{R}_j)}.$$
(2)
According to (2), the information on the spatial structure of the ion chain is contained in $`G^{(1)}(\stackrel{}{r}_1)`$ as long as $`A_{ef}^{(j)}0`$. This is the case as long as, for example, no $`\pi `$-pulse is used for the excitation. The information about the localization of a non-radiating particle is also contained in (2), as $`G^{(1)}(\stackrel{}{r}_1)`$ will look different for different positions of the isotope. This will be discussed in more detail below. Note, however, that a result similar to (2) is obtained also in the case of classical dipoles or classical antennas, i.e. in the case of classical light sources grover2 .

As mentioned in the introduction, we could, however, do more than just measure the far-field intensity, e.g. by examining the quantum features of the emitted radiation. For that purpose let us study the spatial photon - photon correlations produced by the chain of ions by using two photodetectors (see Fig. 1b).

The photon - photon correlations are determined by the expression (see agarwalbook )

$$\stackrel{}{E}^{()}(\stackrel{}{r}_1,t)\stackrel{}{E}^{()}(\stackrel{}{r}_2,t)\stackrel{}{E}^{(+)}(\stackrel{}{r}_2,t)\stackrel{}{E}^{(+)}(\stackrel{}{r}_1,t),$$
(3)
where $`\stackrel{}{r}_2=\stackrel{}{n}_2r`$ defines the position of the second detector and $`\stackrel{}{E}^{()}(\stackrel{}{r}_i,t)`$ denotes the complex conjugate of $`\stackrel{}{E}^{(+)}(\stackrel{}{r}_i,t)`$, $`i=1,2`$. As can be seen from (3), the spatial photon - photon correlations are determined in terms of the atomic/ionic operators through

$$G^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)=A^{}(\stackrel{}{r}_1)A^{}(\stackrel{}{r}_2)A(\stackrel{}{r}_2)A(\stackrel{}{r}_1).$$
(4)
In order to demonstrate how the information on the spatial structure of the chain is contained in the quantity $`G^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$, let us first examine the case of two ions. If the ions are initially prepared in the state $`|e`$, we find (see agarwal02 )

$$G^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)=B(1+\mathrm{cos}(k(\stackrel{}{n}_2\stackrel{}{n}_1)(\stackrel{}{R}_A\stackrel{}{R}_B))),$$
(5)
where $`B`$ is a constant. Obviously, the information on the location of the two ions is contained in the nonclassical agarwal02  interference pattern $`G^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ via the atomic position variables $`\stackrel{}{R}_A`$ and $`\stackrel{}{R}_B`$. Let us consider next a chain of $`N`$ ions of which one is an off-resonant, non-radiating isotope. The explicit calculation in case of $`\pi `$-excitation leads to the following result:

$$G_{p,N}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)=\underset{i,jp(i<j)}{}|\gamma _{ij}(\stackrel{}{r}_1,\stackrel{}{r}_2)|^2,$$
(6)
where $`p=\overline{1,N}`$ stands for the position of the isotope, $`i=\overline{1,N1}`$ and $`j=\overline{2,N}`$ and $`\gamma _{lm}(\stackrel{}{r}_1,\stackrel{}{r}_2)=\alpha _l(\stackrel{}{r}_1)\beta _m(\stackrel{}{r}_2)+\alpha _m(\stackrel{}{r}_1)\beta _l(\stackrel{}{r}_2)`$ with $`\alpha _l(\stackrel{}{r}_1)=\mathrm{e}^{ik\stackrel{}{n}_1\stackrel{}{R}_l}`$ and $`\beta _m(\stackrel{}{r}_2)=\mathrm{e}^{ik\stackrel{}{n}_2\stackrel{}{R}_m}`$ ($`l,m=\overline{1,N}`$). The solution (6) can easily be interpreted as the sum of terms associated with all possible optical path differences between the photons when scattered by two different ions and recorded by the two detectors, on the assumption that the isotope at $`p`$ does not scatter at all. Equation (6) becomes even more transparent for equally spaced ions. In the case of four ions, for example, the system after initial excitation by a $`\pi `$-pulse will be in one of the pure states $`|feee`$, $`|efee`$, $`|eefe`$ or $`|eeef`$. In this case we get $`|\gamma _{12}(\stackrel{}{r}_1,\stackrel{}{r}_2)|^2|\gamma _{23}(\stackrel{}{r}_1,\stackrel{}{r}_2)|^2|\gamma _{34}(\stackrel{}{r}_1,\stackrel{}{r}_2)|^2=2+2c_1(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ and $`|\gamma _{13}(\stackrel{}{r}_1,\stackrel{}{r}_2)|^2|\gamma _{24}(\stackrel{}{r}_1,\stackrel{}{r}_2)|^2=2+2c_2(\stackrel{}{r}_1,\stackrel{}{r}_2)`$, where $`c_m(\stackrel{}{r}_1,\stackrel{}{r}_2)=\mathrm{cos}[k(\stackrel{}{n}_1\stackrel{}{n}_2)(\stackrel{}{R}_1\stackrel{}{R}_{m+1})]`$, so that

$`G_{p,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)=`$ (7)
$`=2(\delta _{p,1}+\delta _{p,4})(3+2c_1(\stackrel{}{r}_1,\stackrel{}{r}_2)+c_2(\stackrel{}{r}_1,\stackrel{}{r}_2))`$
$`+2(\delta _{p,2}+\delta _{p,3})(3+c_1(\stackrel{}{r}_1,\stackrel{}{r}_2)+c_2(\stackrel{}{r}_1,\stackrel{}{r}_2)+c_3(\stackrel{}{r}_1,\stackrel{}{r}_2)).`$
where $`\delta _{p,q}`$ is the Kronecker symbol.

According to (7), the information about the isotope position can be extracted from $`G_{p,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ due to the unique distribution of prefactors (or Fourier coefficients) $`c_j`$ for the different isotope localizations. From (7) it is seen, however, that the outcome for the isotope positions $`p=1`$ ($`2`$) is identical to the outcome for $`p=4`$ ($`3`$). This means that for equally spaced ions an additional measurement would be required to determine whether the isotope is on the left or right hand side with respect to the middle of the chain. A similar argument holds for the general outcome (6): for arbitrarily spaced ions there is a unique interference pattern for each isotope position $`p`$ since the $`\gamma _{ij}`$ will in general be different for different $`i,j`$, so that there is a unique combination of $`\gamma _{ij}`$ for each isotope position $`p`$. By measuring the spatial dependence of the photon-photon correlation function one could thus clearly distinguish the isotope position $`p`$ from any other position $`p^{}`$. Note in particular that the additional degree of freedom given by $`\stackrel{}{r}_2`$ due to the use of two detectors increases the parameter space available. This is exceedingly useful when the information from the first order correlation function $`G_{p,N}^{(1)}(\stackrel{}{r}_1)`$ is difficult to extract. This is the case for, for example, a large number of ions and/or a small ion spacing. There is in addition the particular situation where $`G_{p,N}^{(1)}(\stackrel{}{r}_1)`$ contains no information at all, as is the case, for example, for a $`\pi `$-pulse excitation. In this case the contrast of the interference pattern of $`G_{p,N}^{(1)}(\stackrel{}{r}_1)`$ vanishes, whereas that of $`G_{p,N}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ can still remain maximal skornia01 ; agarwal02 . In Fig. 2, we show the behavior in case of four equally spaced ions for $`G_{p,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ for $`p=1`$ (or $`4`$) and $`p=2`$ (or $`3`$) as a function of $`\varphi _i(\stackrel{}{r}_i)=\mathrm{arctan}(r_{i,y}/r_{i,x})`$, $`i=1,2`$. For the same excitation angle as used in this figure (i.e. $`\pi `$-pulse excitation) $`G_{p,4}^{(1)}(\stackrel{}{r}_1)`$ would correspond to a constant, independently of $`\varphi _1(\stackrel{}{r}_1)`$ and independent of the isotope position $`p`$. For a $`\pi /2`$-pulse excitation on the $`|f|e`$ transition, however, $`G_{p,4}^{(1)}(\stackrel{}{r}_1)`$ would - apart from a prefactor - show an interference pattern similar to $`G_{p,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ for $`\pi `$-pulse excitation and $`\varphi _2(\stackrel{}{r}_2)=0`$ (see horizontal lines in Fig. 2). More precisely, we obtain

$`G_{p,4}^{(1)}(\stackrel{}{r}_1)`$ $`=`$ $`1/2(\delta _{p,1}+\delta _{p,4})(3+2c_1(\stackrel{}{r}_1)+c_2(\stackrel{}{r}_1))`$
$`+`$ $`1/2(\delta _{p,2}+\delta _{p,3})(3+c_1(\stackrel{}{r}_1)+c_2(\stackrel{}{r}_1)+c_3(\stackrel{}{r}_1)).`$
where $`c_m(\stackrel{}{r}_1)=\mathrm{cos}[k\stackrel{}{n}_1(\stackrel{}{R}_1\stackrel{}{R}_{m+1})]`$ and $`p=\overline{1,4}`$. The result (up to a prefactor) is shown in Fig. 3a. Yet, it is here that the additional degree of freedom of $`G_{p,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ with respect to $`G_{p,4}^{(1)}(\stackrel{}{r}_1)`$ becomes important: the additional parameter $`\stackrel{}{r}_2`$ increases the available parameter space and in this manner allows a more flexible and precise search of the isotope. To show that, we plot in Fig. 3b $`G_{p,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ for $`p=1`$ (or $`4`$) and $`p=2`$ (or $`3`$), and with $`\stackrel{}{r}_1`$ and $`\stackrel{}{r}_2`$ such that $`||\varphi _1(\stackrel{}{r}_1)||\varphi _2(\stackrel{}{r}_2)||=1/\pi `$. The latter condition stands for a fixed distance between the two detectors and corresponds in Fig. 2 to a straight line in the $`(\varphi _1(\stackrel{}{r}_1),\varphi _2(\stackrel{}{r}_2))`$-plane (see tilted lines in Fig. 2). As is seen from Fig. 3b, the constraints for the angular resolving power are much more relaxed in comparison with that needed for $`G_{p,4}^{(1)}(\stackrel{}{r}_1)`$ (Fig. 3a).

The spatial second order correlation function also allows the isotope position to be determined in cases where the ions are separated by only $`\lambda /2`$ so that they cannot be individually resolved. This can be demonstrated by analyzing $`G_{1,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ for equally spaced ions with ion separation $`d=\lambda /2`$ (Fig. 4). If the detector positions are chosen such that $`\varphi _1(\stackrel{}{r}_1)=\varphi _2(\stackrel{}{r}_2)`$ \- corresponding in the $`(\varphi _1(\stackrel{}{r}_1),\varphi _2(\stackrel{}{r}_2))`$-plane to a straight line with slope 1 - we resolve for $`G_{1,4}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ a central maximum plus two side maxima in the central region $`\pi /4\varphi _1(\stackrel{}{r}_1)\pi /4`$, whereas for $`G_{1,4}^{(1)}(\stackrel{}{r}_1)`$ only the central maximum is obtained (Fig. 4a). However, it is from the position and amplitude of the side maxima that the information about the isotope position $`p`$ is derived, whereas the central maximum in this regard contains no information.

Let us finally consider nine equally spaced ions in the trap one of which is an isotope. There are $`5`$ different positions for the isotope which can be distinguished from $`G_{p,9}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ (see above). These five different possibilities are plotted in Fig. 5a on the assumption that $`||\varphi _1(\stackrel{}{r}_1)||\varphi _2(\stackrel{}{r}_2)||=1/\pi `$ and a $`\pi `$-pulse is used for excitation. It can be clearly seen from the figure that the requirement for the angular resolving power is far less demanding in comparison with Fig. 5b, where $`G_{p,9}^{(1)}(\stackrel{}{r}_1)`$ (apart from a prefactor) is plotted for a $`\pi /2`$-pulse excitation on the $`|f|e`$ transition. For a relative position of the two detectors such that $`||\mathrm{sin}\varphi _1(\stackrel{}{r}_1)||\mathrm{sin}\varphi _2(\stackrel{}{r}_2)||=0.378`$, one even obtains for $`G_{p,9}^{(2)}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ constant values, with amplitudes which depend on the isotope position $`p`$, but are independent from $`\varphi _1(\stackrel{}{r}_1)`$. These values are depicted in Fig. 5a by dashed lines. However, due to the nonlinear dependence of the two detector positions, this situation might be more difficult to implement experimentally. Nevertheless, such as special angular independent case could play an important role when one would be interested to normalize the correlation function, e.g. for the purpose of a better comparison between the different isotope positions.

In conclusion, it has been shown in a simple search problem using a chain of trapped ions that the search process can be strongly improved when interferences are observed than employing a one by one search. The speed up of the search process results from the fact that the interference patterns are produced by the light of all emitting ions of the chain at once so that there is a simultaneous contribution of all scatterers to the signal. This superposition of the signal is apparently sufficient for the speed up and no entanglement between the ions nor other quantum phenomena are required grover2 . However, quantum interferences allow to increase the parameter space available and in particular improve the precision when a larger number of ions is involved, since they relax the demands for the spatial resolving power of the detectors. Moreover, in the case of quantum interferences an interference pattern is obtained even at ion distances smaller than $`\lambda /2`$ which is the ultimate limit for a classical interference experiment.

Acknowledgement. One of the authors, (G.O.A.) gratefully acknowledges the support of the Alexander von Humboldt Foundation Fellowship.


text:
# Some aspects of the phenomenology of canonical noncommutative spacetimes1footnote 11footnote 1This work is based on Refs. [1] and [2], to which I refer the reader for all details.

## Abstract

I describe some phenomenological contexts in which it is possible to investigate effects induced by (string-motivated) canonical noncommutative spacetime. Due to the peculiar structure of the theory the usual criteria adopted for the choice of experimental contexts in which to test a theory may not be applicable here; care is required in taking into account the effects of IR/UV mixing. This invites one to consider contexts involving particles of relatively high energies, like high-energy cosmic rays and certain high-energy gamma rays observed from distant astrophysical sources.

Canonical noncommutative spacetime is characterized by coordinate noncommutativity of the form:

$$[x^\mu ,x^\nu ]=i\theta ^{\mu \nu }.$$
(1)
An increase of interest in this type of noncommutativity has been recently motivated by the observation that it may be relevant to the description of some aspects of string theories formulated in presence of a background $`B`$ field.

Since they reflect the properties of a background the $`\theta ^{\mu \nu }`$ parameters cannot be observer independent; thus, as a normal background field would, they break 1st ; 2nd ; mst ; carlson  some spacetime symmetries of the theory. This symmetry loss causes one-loop $`E(p)`$ dispersion relations to be deformed. Dispersion relation deformations in (non supersymmetric) $`U(1)`$ theory turn out to be as follows. For photons the modification to the self-energy is proportional to mst

$$\frac{\left(p\theta \right)^\mu \left(p\theta \right)^\nu }{\left(p\theta \right)^4},$$
(2)
(where $`\left(p\theta \right)^\mu `$ is shorthand for $`p_\nu \theta ^{\mu \nu }`$) and thus the deformed dispersion relation is

$$p_0^2=\stackrel{}{p}^2+\frac{C}{\left(p\theta \right)^2};$$
(3)
where $`C`$ is a constant. This dispersion relation deformation is polarization dependent, in the sense that only photons with polarization vector parallel to the direction of the vector $`\left(p\theta \right)^\mu `$ are affected by the modification. For neutral scalars the modification to the self-energy is proportional to $`\left(p\theta \right)^2`$ and thus the dispersion relation is modified in this way:

$$p_0^2=m^2+\stackrel{}{p}^2+\frac{C^{}}{\left(p\theta \right)^2},$$
(4)
with $`C^{}`$ a different constant. Finally, for Majorana fermions the dispersion relation receives only logarithmic corrections, mass suppressed.

A fundamental point are the energy scales suitable for phenomenology. Of course it is necessary to introduce an UV cutoff $`\mathrm{\Lambda }`$, since the theory is of course not viewed as fundamental (for example it does not include gravity). But through the IR/UV mixing mvrs  the lack of knowledge of the UV sector (above $`\mathrm{\Lambda }`$) also implies a lack of knowledge of the structure of the theory at scales lower than222For simplicity I am assuming that the matrix $`\theta ^{\mu \nu }`$ can be described in terms of a single characteristic scale $`\theta `$. $`\mu =\left(\mathrm{\Lambda }\theta \right)^1`$, which should therefore be treated as an IR cutoff in phenomenology.

In presence of canonical noncommutativity our lack of control at energy scales higher than a certain cutoff scale $`\mathrm{\Lambda }`$ implies, because of the IR/UV mixing, also a lack of control at energy scales below $`\mu `$.

This means that testing the theory at very low energies (below the energy scale $`\mu `$) does not yield reliable results. Of course in setting up a phenomenology that focuses on the regime $`\mu =\left(\mathrm{\Lambda }\theta \right)^1<E<\mathrm{\Lambda }`$ a key issue is the value of $`\mu `$. And for some plausible choices of cutoff $`\mathrm{\Lambda }`$ and noncommutativity parameters one finds 1st  that $`\mu `$ can be rather high, possibly as high as $`\mu TeV`$. This implies that to test the theory one must compare it with data corresponding to relatively high energies. In such cases our best chances are clearly found in astrophysics.

The deformation (3) of the photon dispersion relation brings about a deformation of the propagation of light; the speed of light acquires a dependence upon wavelength. A good opportunity to test this effect is provided by Blazars grbgac ; aus , which indeed have been observed billetal  at and above the TeV scale. For some choices of the cutoff $`\mathrm{\Lambda }`$ the analysis of recently-observed Blazars allows 1st  to set limits on $`\theta `$ in the neighborhood of $`\theta 10^{40}cm^2`$ The surprising fact that limit takes the form of an upper bound reflects some features of the IR/UV mixing. In fact, the observations only allow to determine that one of two options is realized: either the cutoff estimate is incorrect or a suitable upper bound on $`\theta `$ must be inforced. The mentioned polarization dependence of the effect implies that, for unpolarized sources, only a fraction of the flux can be used for the analysis.

In the context of observations of cosmic rays an opportunity for canonical noncommutativity to enter the analysis originates from (4), which suggests a modification of the dispersion relation for pions. This would affect the threshold energy requirements for photopion production,

$$p+\gamma =p+\pi ,$$
(5)
which is known to be relevant for the analysis of the interactions between high-energy cosmic rays and CMB (cosmic microwave background) photons. The usual special relativistic analysis of photopion production leads to the prediction that the flux of cosmic rays with energies above a cutoff energy of $`10^{19}`$ eV should be strongly suppressed. The value of this cutoff energy can be significantly higher 1st  if one takes into account the modified dispersion relation (4) for the pion, while keeping unmodified the law of energy-momentum conservation and the dispersion relations for the proton and the CMB photon. Support for an unmodified law of energy-momentum conservation is found 2nd ; douglas  in the analysis of field theories in canonical noncommutative spacetime333Note however that, while in canonical noncommutativity it appears appropriate to assume modified dispersion relations and unmodified laws of energy-momentum conservation, there are other proposals for quantum spacetime in which both a modification of the dispersion relation and a modification of energy-momentum conservation are present dsr ., the assumption of unmodified proton dispersion relation is justified by some general properties gaume  of charged spin-$`1/2`$ high-energy particles in canonical noncommutativity, and we assumed 1st  an unmodified dispersion relation for CMB photons because of their very low typical energies (likely below the mentioned low-energy cutoff of validity of our framework) and the strong experimental basis for an unmodified dispersion relation for low-energy photons. Forthcoming cosmic-ray observatories could be sensitive to the effects of canonical noncommutativity if $`\theta `$ takes a value which is roughly in the neighborhood 1st  of $`10^{38}cm^2`$.