text:
# Binary Tree Summation Monte Carlo Method for Potts Models

## 1 Introduction

Efforts to develop better and more efficient algorithms for Monte Carlo simulations have a long history, in which the Fortuin-Kasteleyn (FK) transformation  for Potts models has played a pivotal role. The most common use of this transformation has been to create algorithms in which clusters of spins in the Potts model are flipped simultaneously . In this paper, we present a new algorithm using the FK representation in the spirit of recent work by Newman and Ziff  on the percolation problem. The method produces independent samples and sums up a large number of configurations for each sweep. The partition function and thermodynamic averages for all values of the temperature $`T`$ can be computed from a single run.

Consider the Potts model with the Hamiltonian,

$$H(\sigma )=J\underset{i,j}{}\delta _{\sigma _i,\sigma _j},$$
(1)
where the summation is over nearest neighbor pair, $`\sigma _i=1,2,\mathrm{},q`$. The FK transformation allows us to write the partition function in the percolation representation as

$$Z=\underset{\mathrm{\Gamma }}{}p^b(1p)^{Mb}q^{N_c(\mathrm{\Gamma })},$$
(2)
where we sum over all configurations of bonds connecting nearest neighbor sites, $`p=1\mathrm{exp}\left(J/(kT)\right)`$, $`b`$ is the number of bonds present, $`M=dN`$ is the maximum number of possible bonds, and $`N_c(\mathrm{\Gamma })`$ is the total number of clusters for a given configuration of bonds.

In order to evaluate thermodynamic averages in this representation, we carry out the summation over all configurations in two steps. First, we sum over the number of bonds $`b`$, and then for each value of $`b`$, we sum over all configurations consistent with that number of bonds. Thus we can write

$$Q(p)=Z^1\underset{b}{}p^b(1p)^{Mb}c_bQ_b,$$
(3)
where $`c_b=_{\mathrm{\Gamma }_b}q^{N_c(\mathrm{\Gamma }_b)}`$ and

$$Q_b=\frac{1}{c_b}\underset{\mathrm{\Gamma }_b}{}q^{N_c(\mathrm{\Gamma }_b)}Q(\mathrm{\Gamma }_b),$$
(4)
where $`\mathrm{\Gamma }_b`$ is a configuration with exactly $`b`$ bonds. If we could compute $`c_b`$ and $`Q_b`$ for some observable $`Q(\mathrm{\Gamma })`$ for every $`b`$, then we could compute the function at any values of $`p`$ or $`T`$.

The problem then reduces to the computation of the normalization constants $`c_b`$ and the expectation values in an ensemble with probability distribution proportional to $`q^{N_c(\mathrm{\Gamma })}`$. For the special case of $`q=1`$ (bond percolation), the values of $`c_b`$ are the binomial coefficients and the configurations with fixed number of bonds are weighted uniformly. For general case, we describe two sampling methods. The first is very simple, and useful conceptually, but exponentially inefficient. The second turns out to be very efficient.

## 2 A Survival and Death Process

Starting with an empty lattice, one sweep consists of repeated application of the following steps until the process dies:

1. Pick an unoccupied neighbor pair at random for the next bond.
2. If inserting a bond

1. does not change the cluster number ($`\mathrm{\Delta }N_c=0`$), accept the configuration;
    2. merge two clusters, so that the cluster number decreases by 1, ($`\mathrm{\Delta }N_c=1`$), accept the configuration with probability $`1/q`$, or reject the configuration and terminate the process (and begin the next sweep from an empty lattice).
3. Take statistics of the survival configurations (with equal weights).

The probability distribution of the configurations generated in this manner with $`b`$ bonds is

$$P(\mathrm{\Gamma }_b)=\frac{b!(Mb)!}{M!}q^{N_c(\mathrm{\Gamma }_b)N_c(\mathrm{\Gamma }_0)}.$$
(5)
The survival probability after $`b`$ bonds is then

$$S_b=\underset{\mathrm{\Gamma }_b}{}P(\mathrm{\Gamma }_b)=\frac{b!(Mb)!}{M!}c_bq^N,$$
(6)
where $`N`$ is the total number of sites. As the survival probability decays exponentially with the number of bonds added, we expect very poor statistics for large values of $`b`$. However, the eventual survival probability, $`S_M`$, is equal to $`q^{N+1}`$, which is nonzero.

## 3 Binary Tree Summation Method

To deal with the exponentially decreasing survival probability, we have developed the following method for summing multiple bond sequences in a single sweep. For a configuration $`\mathrm{\Gamma }_b`$ we count the number of type-0 bonds $`n_0`$ of unoccupied pairs of sites that are on the same cluster, and the number of type-1 bonds $`n_1`$ of unoccupied pairs of sites that would connect two different clusters. Clearly $`n_0+n_1+b=M`$. At each step, we pick a type-1 bond with equal probability from among all current type-1 bonds. Starting with an empty lattice and merging clusters at each step, we continue until all sites are members of the same cluster. Along the way, we collect statistics $`Q(i)`$ for each of the $`N`$ configurations generated.

After a sweep has been carried out, we construct all possible paths that a full simulation of both type-0 and type-1 bonds would have taken if we had joined clusters in exactly the same sequence as in the actual simulation, but had also inserted a random number of type-0 bonds at each step. At each choice between types of bond, the type-1 bond is given a relative weight of $`n_1/q`$ and a type-0 bond $`n_0`$. The total weight for a path starting from an empty lattice to a particular configuration is the product of the factors $`n_0`$ or $`n_1/q`$ depending on the path taken.

To understand this algorithm, imagine that we had actually followed the branching process. At each step, a configuration may split into two configurations, one with a type-1 bond added (with probability $`1/n_1`$) and one with a type-0 bond added (with probability of $`1/n_0`$). The probability of appearance of a particular path is $`[n_{f(k)}(\mathrm{\Gamma }_k)]^1`$, where $`f(k)=0`$ or 1 depending on the choice of type-0 or 1 bond. In taking statistics, we have weighted with the inverse of the factor, multiplied by additional factors of $`q`$ proportional to the number of clusters $`N_c(\mathrm{\Gamma })`$. The net effect is the required sample average with overall weight $`q^{N_c(\mathrm{\Gamma })}`$. Now the key observation is that we do not need an explicit simulation for the type-0 bonds, since the type-0 bonds have no effect on the measured quantities.

The binary tree summation algorithm has several attractive features. First, the usual slowing down due to correlation between samples is absent. Each sweep is independent. Second, for each sweep, data for a very large number of samples are collected. Although they are highly correlated, an exponentially large number of paths can be summed efficiently with $`O(N^2)`$ operations per sweep. Third, unlike multicanonical simulations  or the flat histogram method , there are no unknown weighting factors to determine.

## 4 Implementation

Let $`w(b,i)`$ be the weight of the total contributions from all possible paths to the state specified by the number of bonds $`b`$ and merge sequence number $`i`$. These quantities can be calculated from the starting condition $`w(0,i)=\delta _{i,0}`$, the constraint that $`n_0(b,i)0`$, and the recursive equation

$$w(b+1,i)=w(b,i)n_0(b,i)+w(b,i1)n_1(b,i1)/q,$$
(7)
where $`n_1(b,i)=n_1(i)`$, $`n_0(b,i)=n_0(i)b+i`$. The value $`w(b,i)`$ is nonzero only for $`bi`$. The computation of the weights $`w(b,i)`$ is similar to that of binomial coefficients. The final contribution to the statistics at $`b`$ number of bonds is

$$Q_b=W_b^1\underset{i=0}{\overset{N1}{}}w(b,i)Q(i),$$
(8)
where the average is over simulation sweeps, and $`W_b=_iw(b,i)`$ is the total weight at a given $`b`$. We can then show that

$$W_b=\frac{M!}{(Mb)!}S_b,$$
(9)
which allows us to compute $`c_b`$. The result is numerically identical to computing the conditional survival probability $`S_{b+1}/S_b`$ from the expectation value of $`Q=(n_0+n_1/q)/(Mb)`$.

During the simulation, the values of $`n_0`$ and $`n_1`$ can be updated efficiently. For each cluster, we keep a list of unoccupied bonds with other clusters. When two clusters (A and B) are joined, we merge the smaller one with the larger one, and remove and count the number of bonds $`n_{\mathrm{AB}}`$ connecting the clusters. We update according to $`n_0n_0+n_{\mathrm{AB}}1`$, $`n_1n_1n_{\mathrm{AB}}`$. The timing of our program shows that this part of the algorithm scales nearly linearly with number of sites $`N`$, as expected.

## 5 Results and Discussions

We note that the coefficients $`c_b`$ are related to the density of states $`n(E)`$, which gives the coefficients of the partition function polynomial in the variable $`\mathrm{exp}\left(J/(kT)\right)`$. In the FK percolation representation, it becomes a polynomial in $`p/(1p)`$. By a proper change of variables, we can find exact results  of $`c_b`$ for the two-dimensional Ising model.

We define the following errors to test our method against exact results for the two-dimensional Ising model:

$`ϵ_0`$ $`=`$ $`\left|q^{N1}c_M/c_01\right|,`$ (10)
$`ϵ_1`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{b=0}{\overset{M}{}}}\left|c_b/c_b^{\mathrm{exact}}1\right|,`$ (11)
$`ϵ_{\mathrm{MAX}}^Q`$ $`=`$ $`\underset{T}{\mathrm{max}}\left|Q(T)Q^{\mathrm{exact}}(T)\right|,`$ (12)
$`ϵ_{\mathrm{AVE}}^Q`$ $`=`$ $`{\displaystyle _0^1}𝑑x\left|Q\left(T(x)\right)Q^{\mathrm{exact}}\left(T(x)\right)\right|,T(x)={\displaystyle \frac{x}{1x}}.`$ (13)
Due to a special cancellation for this algorithm, $`ϵ_0`$ is exactly zero. The error $`ϵ_1`$, along with the maximum and average errors in the energy and specific heat, are listed in Table 1. Since the algorithm asymptotically takes $`O(N^2)`$ operations per sweep, a fair comparison with other methods should compare the total CPU times. A comparable N-fold way transition matrix Monte Carlo (TMMC) run took 1.9 microsecond per sweep per site on the same machine. Thus, the present method is superior for small lattices with linear size $`L16`$. For $`L=32`$ it is comparable to TMMC . For much larger lattices, it becomes less favorable mainly due to the $`O(N^2)`$ nature of the algorithm. It is quite likely that we can speed up the computation using special properties of the weights.

Our method is applicable for Potts models with any number of states, including fractional or negative values. Work is currently in progress to apply this algorithm to a number of problems of interest.


format: text

text:
# Non-adiabatic coupling and adiabatic population transfer in quantum molecular systems

## Abstract

We show that a counter-intuitive pulse sequence leads to adiabatic passage between the vibrational levels of three harmonic potentials through parallel dark states in adiabatic approximation. However, the adiabatic assumptions break down for very intense pulses and non-adiabatic couplings result in the population transfer by light-induced potential shaping.

Population transfer by shaping light induced potentials (LIP) is a very robust scheme to transfer vibrational population between different electronic potentials in diatomic molecules newPRL ; SolPRA00 ; us ; MalCPH01 ; GarPRL98 ; RodPRA00 . One of the main advantages of the scheme is that the vibrational quantum number is conserved during the overall process. Since the selectivity of the transfer is guaranteed by a symmetry rule, the method is not very sensitive to the exact position of the energy resonance, i.e. to the topological or energetic features of the potentials involved.

The simplest scenario of population transfer by shaping LIP can be described in terms of three electronic states (potential curves) sequentially coupled by two laser pulses, and the vibrational population is driven by two photon off resonant absorption. Although there are several possible scenarios for selective population transfer newPRL ; BanPRA01 , the first scheme proposed, called APLIP, involves a sequence of two very strong (TW/cm2) picosecond pulses applied in counter-intuitive order GarPRL98 . The pulse that couples the intermediate electronic state with the final excited electronic state must precede the pulse that couples the ground electronic state with the intermediate one. Then the system dynamics follows a characteristic pattern in the vibrational basis representation, which is illustrated in Fig. 1 for a test system of three harmonic oscillators with the same force constant. (More details about the model are presented below.) According to Fig. 1(a), the overall population on the ground electronic state is rapidly transferred to the final electronic state, while very few population is temporally excited to the intermediate electronic state. So, in principal it is possible to suppose there is a dark state consisting of the initial and final potentials, similar to the population transfer through the dark state of the three-level system BeRMP98 . Looking into the vibrational populations in detail, it can be seen that the passage is mediated by substantial excitation of high energy vibrational levels in both initial and final potentials, Fig.1(b).

The scenario previously considered is very similar to the well known Stimulated Raman Adiabatic Passage (STIRAP) scheme BeRMP98 . In STIRAP there are three levels sequentially coupled by two fields working in counter-intuitive order. This arrangement of pulses prepares the system initially in a dark state, $`|\mathrm{\Phi }_0`$, which is an eigenstate of the Hamiltonian dressed by the field in the adiabatic approximation. The dark state correlates at late times with the final level $`|3`$ and never overlaps with the intermediate level. The population is transferred from $`|1`$ to $`|3`$ in a completely similar way to the overall electronic populations of Fig. 1(a).

Normally, the STIRAP scheme involves nanosecond (or longer) pulses with intensities in the GW/cm2 (or less). Therefore it is natural to ask if APLIP is a short time, strong field version of STIRAP (and applied to electronic states instead of vibrational or atomic levels). Indeed the requirement of stronger laser fields can be understood under the general assumptions of STIRAP. Since the adiabaticity condition in STIRAP is usually expressed by the pulse area relation, $`\mathrm{\Omega }\tau 1`$, where $`\mathrm{\Omega }`$ is the effective Rabi frequency of the pulses and $`\tau `$ their time widths, any reduction in the pulse width must imply a corresponding increase in the pulse amplitude. Nevertheless, in this note we show that the population transfer in APLIP can not be reduced to the generalized description of STIRAP in multi-level system. In order to do so we show that the natural extension of STIRAP to the three potential curves scenario in the strong field limit predicts a different dynamical evolution than the one shown in Fig. 1.

Let us consider in detail the Hamiltonian of the system. The Hamiltonian for three electronic potentials ($`V_i(x)`$, $`i=1,2,3`$) coupled by two laser fields ($`E_1(t)`$ and $`E_2(t)`$) in the Born Oppenheimer and rotating wave approximation (RWA) reads

$$H^{RWA}=\left(\begin{array}{ccc}T+U_1(x)& \frac{1}{2}\mu _{12}(x)E_1(t)& 0\\ \frac{1}{2}\mu _{12}(x)E_1(t)& T+U_2(x)& \frac{1}{2}\mu _{23}(x)E_2(t)\\ 0& \frac{1}{2}\mu _{23}(x)E_2(t)& T+U_3(x)\end{array}\right),$$
(1)
where $`T`$ is the kinetic energy operator and $`U_i(x)`$ are the diabatic dressed potentials ($`U_1(x)=V_1(x)`$, $`U_2(x)=V_2(x)\omega _1`$, $`U_3(x)=V_3(x)\omega _1\omega _2`$, $`\omega _i`$ being the carrier frequencies of the lasers), $`\mu _{12}(x)`$ and $`\mu _{23}(x)`$ are the dipole moments (atomic units are used throughout). In this representation the global wave function of the system is expanded in terms of electronic wave functions, $`\{\mathrm{\Xi }_i(q;x)\}`$, and vibrational wave packets, $`\{\psi _i(x,t)\}`$, where $`q`$ is an index representing the collective electron coordinates and $`x`$ is the vibrational coordinate.

After expanding the vibrational wave packet in the set of eigenfunctions of the vibrational Hamiltonian for each potential,

$$\psi _\alpha (x,t)=\underset{j}{}d_j^{(\alpha )}(t)\varphi _j^{(\alpha )}(x),$$
(2)
where $`\varphi _j^{(\alpha )}(x)`$ is the j’th vibrational level in the electronic state $`\alpha `$ (in the following we use Greek letters to designate electronic states and Roman letters to designate vibrational levels), we arrive at the following time dependent Schrödinger equation (TDSE):

$$\{\begin{array}{cc}& i\dot{d}_j^{(1)}=(\omega _j^{(1)}+\omega _1+D_1^0D_2^0)d_j^{(1)}_k\frac{\mathrm{\Omega }_{jk}^{(1)}}{2}d_k^{(2)}\hfill \\ & i\dot{d}_j^{(2)}=_k\frac{\mathrm{\Omega }_{jk}^{(1)}}{2}d_k^{(1)}+(\omega _j^{(2)}D_2^0)d_j^{(2)}_k\frac{\mathrm{\Omega }_{jk}^{(2)}}{2}d_k^{(3)}\hfill \\ & i\dot{d}_j^{(3)}=_k\frac{\mathrm{\Omega }_{jk}^{(2)}}{2}d_k^{(2)}+(\omega _j^{(3)}\omega _2+D_3^0D_2^0)d_j^{(3)},\hfill \end{array}$$
(3)
where $`\omega _j^{(\alpha )}+D_{(\alpha )}^0`$ is the eigenvalue corresponding to the $`|\varphi _j^{(\alpha )}`$ eigenfunction, $`\mathrm{\Omega }_{ij}^{(1)}=E_1(t)\varphi _i^{(1)}\mathrm{\Xi }_1|\mu _1|\varphi _j^{(2)}\mathrm{\Xi }_2`$ and $`\mathrm{\Omega }_{ij}^{(2)}=E_2(t)\varphi _i^{(3)}\mathrm{\Xi }_3|\mu _2|\varphi _j^{(2)}\mathrm{\Xi }_2`$ are the Rabi frequencies, and $`D_\alpha ^0`$ are the potential zero energies.

In order to compare the dynamics of the general system with that of STIRAP, we establish a correspondence between the $`3\times N`$ equations of motion (Eq. (3)) and the equations for $`N`$ $`3`$-level ladder systems, where both diabatic and adiabatic states are known, corresponding to that of the STIRAP Hamiltonian. In doing this connection we neglect all the contributions from the continuum wave functions in all the electronic states. To simplify the notation we use a model of three symmetrically displaced harmonic oscillators (Fig. 2) with the same force constant (SDHO model), so that the separation between the minima of the potentials is constant, $`r_0^{(3)}r_0^{(2)}=r_0^{(2)}r_0^{(1)}=R`$. We assume that the equilibrium configurations of the excited potentials are displaced to larger inter nuclei distances. Furthermore we consider only processes in two photon resonance. This allows us to define a constant energy splitting, $`\delta \omega =\omega _{n+1}^{(1)}\omega _n^{(1)}=\omega _{n+1}^{(3)}\omega _n^{(3)}`$, and a one photon detuning, $`\mathrm{\Delta }=D_2^0\omega _0^{(1)}\omega _1=(\omega _2+D_2^0)\omega _0^{(3)}D_3^0`$. The TDSE is:

$$\{\begin{array}{cc}& i\dot{d}_j^{(1)}=(j\delta \omega \mathrm{\Delta })d_j^{(1)}_k\frac{\mathrm{\Omega }_{jk}^{(1)}}{2}d_k^{(2)}\hfill \\ & i\dot{d}_j^{(2)}=_k\frac{\mathrm{\Omega }_{jk}^{(1)}}{2}d_k^{(1)}+(j\delta \omega )d_j^{(2)}_k\frac{\mathrm{\Omega }_{jk}^{(2)}}{2}d_k^{(3)}\hfill \\ & i\dot{d}_j^{(3)}=_k\frac{\mathrm{\Omega }_{jk}^{(2)}}{2}d_k^{(2)}+(j\delta \omega \mathrm{\Delta })d_j^{(3)}.\hfill \end{array}$$
(4)
Finally we invoke the Condon approximation. Therefore, the state to state Rabi frequencies can be expressed as $`\mathrm{\Omega }_{ij}^{(1)}=p_{ij}^{12}\mu _1E_1(t)`$ and $`\mathrm{\Omega }_{ij}^{(2)}=p_{ij}^{23}\mu _2E_2(t)`$, where we have defined the geometric Franck-Condon parameters $`p_{ij}^{\alpha \beta }=\varphi _i^{(\alpha )}|\varphi _j^{(\beta )}`$. This is the model used to obtain the results shown in Fig. 1, where we have chosen $`\delta \omega =710^4`$ a.u., $`\mathrm{\Delta }=0.015`$ a.u., $`\mu _1=\mu _2=1`$, $`E_i(t)=E_0S(t)`$, with $`E_0=0.08`$ a.u. and $`S(t)`$ being a Gaussian envelope function. Both pulses have the same envelope function, with width $`\sigma =2.5`$ ps and they are time delayed, so that the second pulse precedes in $`2.5`$ ps the first pulse. To obtain the numerical results shown in Fig. 1 we have numerically integrated the TDSE based on a grid discretization of the Hamiltonian in Eq. (1) and not on the discrete basis representation of Eq. (4). Details of the numerical propagator are given in reference SolPRA00 .

In order to go further in the comparison between the set of Eqs. (4) and those of STIRAP, we change the representation from the diabatic basis $`\{\varphi _n^{(1)},\varphi _n^{(2)},\varphi _n^{(3)}\}`$ to the adiabatic basis of the 3-levels ladder system for each $`n`$ sub-system, $`\{\mathrm{\Phi }_+^{(n)},\mathrm{\Phi }_0^{(n)},\mathrm{\Phi }_{}^{(n)}\}`$, using the block diagonal rotation matrix $`\widehat{R}=\widehat{R}_1\mathrm{}\widehat{R}_n\mathrm{}\widehat{R}_N`$ with

$$\widehat{R}_n=\left(\begin{array}{ccc}\mathrm{sin}\phi _n\mathrm{sin}\theta _n& \mathrm{cos}\theta _n& \mathrm{cos}\phi _n\mathrm{sin}\theta _n\\ & & \\ \mathrm{cos}\phi _n& 0& \mathrm{sin}\phi _n\\ & & \\ \mathrm{sin}\phi _n\mathrm{cos}\theta _n& \mathrm{sin}\theta _n& \mathrm{cos}\phi _n\mathrm{cos}\theta _n\end{array}\right),$$
(5)
where the angles are defined by $`\mathrm{tan}\theta _n=\mathrm{\Omega }_{nn}^{(1)}/\mathrm{\Omega }_{nn}^{(2)}`$ and $`\mathrm{tan}(2\phi _n)=2\left(\sqrt{(\mathrm{\Omega }_{nn}^{(1)})^2+(\mathrm{\Omega }_{nn}^{(2)})^2}\right)/\mathrm{\Delta }`$.

In the new representation, the dynamics of the system is followed by the amplitude coefficients $`\{a_+^{(n)},a_0^{(n)},a_{}^{(n)}\}`$, which are delocalized and follow the electronic transitions between the potentials at different moments of time. The new basis is quasi-adiabatic since the transformation matrix $`\widehat{R}`$ diagonalizes only each $`3\times 3`$ sub-system matrices, while there remain couplings between different sub-systems.

To illustrate the nature of these couplings, we detail the analysis for a system of 6-levels obtained by truncating the expansion of the wave function (Eq. (2)) to only the first two vibrational levels. We consider these levels as belonging to two coupled 3-level ladder systems, whose Hamiltonian is

$$H=\left(\begin{array}{cc}\widehat{H}_{11}& \widehat{H}_{12}\\ & \\ \widehat{H}_{21}& \widehat{H}_{22}\end{array}\right)=\frac{1}{2}\left(\begin{array}{ccccccc}0& \mathrm{\Omega }_{11}^{(1)}& 0& 0& \mathrm{\Omega }_{12}^{(1)}& 0& \\ & & & & & & \\ \mathrm{\Omega }_{11}^{(1)}& 2\mathrm{\Delta }& \mathrm{\Omega }_{11}^{(2)}& \mathrm{\Omega }_{12}^{(1)}& 0& \mathrm{\Omega }_{12}^{(2)}& \\ & & & & & & \\ 0& \mathrm{\Omega }_{11}^{(2)}& 0& 0& \mathrm{\Omega }_{12}^{(2)}& 0& \\ & & & & & & \\ 0& \mathrm{\Omega }_{21}^{(1)}& 0& 2\delta \omega & \mathrm{\Omega }_{22}^{(1)}& 0& \\ & & & & & & \\ \mathrm{\Omega }_{21}^{(1)}& 0& \mathrm{\Omega }_{21}^{(2)}& \mathrm{\Omega }_{22}^{(1)}& 2(\delta \omega \mathrm{\Delta })& \mathrm{\Omega }_{22}^{(2)}& \\ & & & & & & \\ 0& \mathrm{\Omega }_{21}^{(2)}& 0& 0& \mathrm{\Omega }_{22}^{(2)}& 2\delta \omega & \end{array}\right),$$
(6)
and we next change the representation using $`\widehat{R}=\widehat{R}_1\widehat{R}_2`$ ($`\widehat{R}_n`$ is given by Eq. (5)). We obtain the quasi-adiabatic Hamiltonian

$$\begin{array}{ccc}\widehat{H}^{CD}\hfill & =\left(\begin{array}{cc}\widehat{R}_1^1& 0\\ 0& \widehat{R}_2^1\end{array}\right)\left(\begin{array}{cc}\widehat{H}_{11}& \widehat{V}_{12}\\ \widehat{V}_{21}& \widehat{H}_{22}+\delta \omega \widehat{I}\end{array}\right)\left(\begin{array}{cc}\widehat{R}_1& 0\\ 0& \widehat{R}_2\end{array}\right)\hfill & \\ & & \\ & & \\ & =\left(\begin{array}{cc}\widehat{H}^{(1)}& \widehat{R}_1^1\widehat{V}_{12}\widehat{R}_2\\ \widehat{R}_2^1\widehat{H}_{21}\widehat{R}_1& \widehat{H}^{(2)}+\delta \omega \widehat{I}\end{array}\right).\hfill & \end{array}$$
(7)
In this representation each sub-system Hamiltonian $`\widehat{H}^{(n)}`$ is diagonal. Neglecting the coupling between sub-systems, $`\widehat{H}_{12}^{QA}=\widehat{R}_1^1\widehat{V}_{12}\widehat{R}_20`$, and using a counter-intuitive sequence of pulses, the Hamiltonian has two parallel (independent) dark states as in STIRAP. Therefore, the differences between APLIP and STIRAP must come from the inter sub-system couplings. Substituting $`\mathrm{cos}(\theta _n)=\mathrm{\Omega }_{nn}^{(2)}/\mathrm{\Omega }_{nn}^{(e)}`$ and $`\mathrm{sin}(\theta _n)=\mathrm{\Omega }_{nn}^{(1)}/\mathrm{\Omega }_{nn}^{(e)}`$ where $`\mathrm{\Omega }_{nn}^{(e)}=\left((\mathrm{\Omega }_{nn}^{(1)})^2+(\mathrm{\Omega }_{nn}^{(2)})^2\right)^{\frac{1}{2}}`$, we obtain, for our reduced 6-level Hamiltonian:

$$\widehat{H}_{12}^{QA}=(\widehat{H}_{21}^{QA})^{}=\left(\begin{array}{ccc}\begin{array}{c}\frac{\mathrm{sin}\phi _1\mathrm{cos}\phi _2}{\mathrm{\Omega }_{11}^{\left(e\right)}}\left(\mathrm{\Omega }_{11}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}+\mathrm{\Omega }_{11}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \\ +\frac{\mathrm{cos}\phi _1\mathrm{sin}\phi _2}{\mathrm{\Omega }_{22}^{\left(e\right)}}\left(\mathrm{\Omega }_{22}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}+\mathrm{\Omega }_{22}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \end{array}& \frac{\mathrm{cos}\phi _1}{\mathrm{\Omega }_{22}^{\left(e\right)}}\left(\mathrm{\Omega }_{22}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\mathrm{\Omega }_{22}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(1\right)}\right)& \begin{array}{c}\frac{\mathrm{sin}\phi _1\mathrm{sin}\phi _2}{\mathrm{\Omega }_{11}^{\left(e\right)}}\left(\mathrm{\Omega }_{11}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}+\mathrm{\Omega }_{11}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \\ +\frac{\mathrm{cos}\phi _1\mathrm{cos}\phi _2}{\mathrm{\Omega }_{22}^{\left(e\right)}}\left(\mathrm{\Omega }_{22}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}\mathrm{\Omega }_{22}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \end{array}\\ & & \\ \frac{\mathrm{cos}\phi _2}{\mathrm{\Omega }_{11}^{\left(e\right)}}\left(\mathrm{\Omega }_{11}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\mathrm{\Omega }_{11}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(1\right)}\right)& 0& \frac{\mathrm{sin}\phi _2}{\mathrm{\Omega }_{11}^{\left(e\right)}}\left(\mathrm{\Omega }_{11}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\mathrm{\Omega }_{11}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(1\right)}\right)\\ & & \\ \begin{array}{c}\frac{\mathrm{cos}\phi _1\mathrm{cos}\phi _2}{\mathrm{\Omega }_{11}^{\left(e\right)}}\left(\mathrm{\Omega }_{11}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}+\mathrm{\Omega }_{11}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \\ \frac{\mathrm{sin}\phi _1\mathrm{sin}\phi _2}{\mathrm{\Omega }_{22}^{\left(e\right)}}\left(\mathrm{\Omega }_{22}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}+\mathrm{\Omega }_{22}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \end{array}& \frac{\mathrm{sin}\phi _2}{\mathrm{\Omega }_{22}^{\left(e\right)}}\left(\mathrm{\Omega }_{22}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}\mathrm{\Omega }_{22}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(1\right)}\right)& \begin{array}{c}\frac{\mathrm{cos}\phi _1\mathrm{sin}\phi _2}{\mathrm{\Omega }_{11}^{\left(e\right)}}\left(\mathrm{\Omega }_{11}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}+\mathrm{\Omega }_{11}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \\ \frac{\mathrm{sin}\phi _1\mathrm{cos}\phi _2}{\mathrm{\Omega }_{22}^{\left(e\right)}}\left(\mathrm{\Omega }_{22}^{\left(1\right)}\mathrm{\Omega }_{12}^{\left(1\right)}+\mathrm{\Omega }_{22}^{\left(2\right)}\mathrm{\Omega }_{12}^{\left(2\right)}\right)\hfill \end{array}\end{array}\right).$$
(8)
Now let us concentrate in the case when both Rabi frequencies are equal, $`\mu _{12}A_1=\mu _{23}A_2=\mathrm{\Omega }_0`$. We consider the initial state to be the ground vibrational level, $`|\varphi _0^{(1)}`$ which initially correlates with $`|\mathrm{\Phi }_0^{(1)}`$. Then the only non zero couplings between this state and any other quasi-adiabatic state are

$$\mathrm{\Phi }_+^{(2)}|H_{21}^{QA}|\mathrm{\Phi }_0^{(1)}=\mathrm{cos}\phi _2(p_{11}^{12}p_{21}^{23}p_{21}^{12}p_{11}^{23})\frac{\mathrm{\Omega }_0S_1(t)S_2(t)}{\sqrt{(\mathrm{\Omega }_{11}^{(1)}(t))^2+(\mathrm{\Omega }_{11}^{(2)}(t))^2}},$$
(9)
and

$$\mathrm{\Phi }_{}^{(2)}|H_{21}^{QA}|\mathrm{\Phi }_0^{(1)}=\mathrm{sin}\phi _2(p_{11}^{12}p_{21}^{23}p_{21}^{12}p_{11}^{23})\frac{\mathrm{\Omega }_0S_1(t)S_2(t)}{\sqrt{(\mathrm{\Omega }_{11}^{(1)}(t))^2+(\mathrm{\Omega }_{11}^{(2)}(t))^2}}.$$
(10)
Equations (9), and (10) separate the molecular contribution (geometrical factors) from the pulse contribution (depending on the pulse amplitude). However it can be seen that the APLIP dynamic behavior cannot stem from these terms. For instance, in the model of symmetrically displaced harmonic oscillators the symmetry of the system imposes that $`p_{12}^{12}=p_{21}^{23}`$ and $`p_{11}^{23}=p_{11}^{12}`$. Therefore the quasi-adiabatic initially populated state, $`|\mathrm{\Phi }_0^{(1)}`$, is not coupled with the rest of the system. So, this state is indeed the adiabatic state whose evolution reproduces exactly the STIRAP behavior independently from the detuning or intensity of the pulses.

Increasing the number of vibrational levels in the expansion (Eq. (2)) does not introduce any different type of couplings and the same conclusions apply. It is very simple to generalize the expression for the general $`3N`$-level system, in this case, the Hamiltonian matrix terms are:

$$\begin{array}{cc}& \mathrm{\Phi }_0^{(n)}|H^{QA}|\mathrm{\Phi }_0^{(m)}=0\hfill \\ & \mathrm{\Phi }_0^{(n)}|H^{QA}|\mathrm{\Phi }_+^{(m)}=\mathrm{\Omega }_0\mathrm{cos}\phi _m(p_{nn}^{12}p_{nm}^{23}p_{nm}^{12}p_{nn}^{23})\frac{S_1(t)S_2(t)}{\sqrt{(p_{nn}^{12}S_1(t))^2+(p_{nn}^{23}S_2(t))^2}}\hfill \\ & \mathrm{\Phi }_0^{(n)}|H^{QA}|\mathrm{\Phi }_{}^{(m)}=\mathrm{\Omega }_0\mathrm{sin}\phi _m(p_{nn}^{12}p_{nm}^{23}p_{nm}^{12}p_{nn}^{23})\frac{S_1(t)S_2(t)}{\sqrt{(p_{nn}^{12}S_1(t))^2+(p_{nn}^{23}S_2(t))^2}}.\hfill \end{array}$$
(11)
For the SDHO model, the geometrical symmetry requires that: $`p_{11}^{12}=\varphi _1^{(1)}|\varphi _1^{(2)}=\varphi _1^{(2)}|\varphi _1^{(3)}=p_{11}^{23}`$, and $`p_{12}^{12}=\varphi _1^{(1)}|\varphi _2^{(2)}=\varphi _2^{(2)}|\varphi _1^{(3)}=p_{21}^{23}`$. Once again, this makes any couplings involving the $`|\mathrm{\Phi }_0^{(n)}`$ states identically zero. The system exhibits $`N`$ parallel (independent) trapped states. Therefore, if the initial wave function is a coherent superposition of vibrational levels of the ground electronic state, $`|\mathrm{\Psi }(x,0)=_{n=0}^Na_n(0)|\varphi _n^{(1)}(x)`$, a sequence of two coherent counter-intuitive pulses will read this wave function in the adiabatic representation as $`|\mathrm{\Psi }(x,0)=_{n=0}^Na_n(0)|\mathrm{\Phi }_0^{(n)}(x,0)`$, that will adiabatically passage at final times to $`|\mathrm{\Psi }(x,T)=_{n=0}^Na_n(0)\mathrm{exp}\left(i(n1)\delta \omega T\right)|\varphi _n^{(3)}(x)`$ (where the phase factor only depends on the eigenvalues of the dressed states $`|\mathrm{\Phi }_0^{(n)}`$, i.e. the geometrical or Berry phase is zero). At each instant of time only the initial levels $`|\varphi _n^{(1)}`$ and the corresponding $`|\varphi _n^{(3)}`$ levels are populated, the final population in $`|\varphi _n^{(3)}`$ being completely determined by the initial population, $`|a_n(0)|^2`$. The same conclusion applies for the final probabilities if the initial state is an incoherent sum of vibrational levels. The dynamics follows as $`N`$ independent STIRAP type systems and does not reproduce the characteristic features of APLIP.

This result is in flagrant contradiction with the numerical evidence (direct solution of the TDSE with the Hamiltonian of Eq.(1)), in which many of high energy vibrational levels of the ground and final potentials are considerably populated (Fig. 1(b)). Although the analytic proof derived in this paper is valid for a specific “ideal” system, the SDHO model, of course any small asymmetry in the model cannot explain the clear difference in behavior of STIRAP and APLIP. The difference, therefore, must stem from a different source. The validity of the quasi-adiabatic Hamiltonian rests upon the adiabatic approximation, that is, the neglect of the contribution from all the terms coming from $`\widehat{R}^1\dot{\widehat{R}}=\widehat{R}_1^1\dot{\widehat{R}}_1\widehat{R}_2^1\dot{\widehat{R}}_2\mathrm{}\widehat{R}_N^1\dot{\widehat{R}}_N`$, that couple states $`\mathrm{\Phi }_0^{(n)}`$ and $`\mathrm{\Phi }_\pm ^{(n)}`$ belonging to the same sub-system. In STIRAP, the adiabatic condition, $`\mathrm{\Omega }_0\tau 1`$ guarantees that these terms can be neglected. Surprisingly enough, for very large $`\mathrm{\Omega }_0`$ and detuning, the APLIP dynamics shows that the contribution of these terms cannot be neglected. This is in agreement with the numerical results obtained for population transfer by shaping LIPs in more general scenarios us ; RodPRA00 .

I.R.S. gratefully acknowledges support from the Secretaría de Estado de Educación y Universidades (Spanish Government).


text:
# QUANTUM NONLOCALITY AND QUANTUM DYNAMICS

## .1 Introduction

It was recently claimed in ref.Gisin  that the usual ”quantum statics” supplemented by the no-signaling condition implies a completely positive, linear quantum dynamical law. The proof offered therein relies on the central notion that a mixed quantum state described by a density matrix $`\widehat{\rho }`$ is always equivalent to, or can be prepared as, a classical probabilistic mixture of pure states $`\{p_i,|\psi _i\}`$, where $`p_i`$ denotes the probability of occurrence of the pure state $`|\psi _i`$. This assumption was criticized in ref.Bona  as leading directly to linearity, regardless of the no-signaling condition. We wish to give an expounded version of this criticism, and reexamine the role played by quantum nonlocality in enforcing dynamical linearity. We also review the issue of seemingly viable nonlinear alternatives, particularly the interesting class of nonlinear dynamics that retains the quantum theory of pure states in its linear form. For these cases we point out that the no-signaling condition leads to a simple criterion of dynamical (non)linearity: Quantum dynamics is (non)linear if and only if the entangled pure states of noninteracting systems evolve (non)linearly .

## .2 No-signaling and the linearity of quantum dynamics

Let us begin by recalling the argument of ref.Gisin . Consider an arbitrary quantum system, prepared in a mixed quantum state represented by the probabilistic mixture $`\{p_i,|\psi _i\}`$. If the density matrix corresponding to this mixture reads

$$\widehat{\rho }=\underset{i}{}p_i|\psi _i\psi _i|,$$
(1)
a dynamical map $`\widehat{\rho }𝐠(\widehat{\rho })`$, not necessarily linear, decomposes as a superposition of individual evolutions of its pure state components $`|\psi _i`$. In other words,

$$𝐠(\widehat{\rho })g\left(\underset{i}{}p_i|\psi _i\psi _i|\right)=\underset{i}{}p_ig\left(|\psi _i\psi _i|\right),$$
(2)
where $`𝐠\left(|\psi _i\psi _i|\right)`$ need not be a pure state. Since this identity must apply to every probabilistic mixture described by the density matrix $`\widehat{\rho }`$, it is concluded that the dynamical map $`𝐠`$ must be linear.

The no-signaling condition enters the proof via an argument regarding the reduced state of a quantum system in an entangled pure state with a space-like separated counterpart. Since such a state cannot be prepared locally, one may question whether it represents an essentially indecomposable quantum state, i.e., an elementary mixture Bona-1 ; Gyftopoulos  \[although this terminology is not employed in ref.Gisin \]. On the other hand, the assumption of a potentially nonlinear dynamics necessarily implies the possibility that elementary and probabilistic mixtures may evolve differently. Suppose the latter claim is true, and \[elementary\] mixtures generated as reduced local states corresponding to entangled pure states may evolve differently than probabilistic mixtures with identical density matrices. In ref.Gisin  it is shown that such a reduced state can be converted into any genuine probabilistic mixture with an identical density matrix by a remote maximal measurement on the entangled counterpart. If the no-signaling condition must apply, then the reduced state must necessarily evolve in the same linear manner as a probabilistic mixture, or the remote measurement could be detected locally. An additional argument concerning the complete positivity of the evolution concludes the proof.

Nevertheless, the above recourse to the no-signaling condition is not necessary once one accepts the interpretation of Eq.(2) that probabilistic mixtures must evolve linearly. It is sufficient to note instead that a similar statement must apply to probabilistic mixtures of pure entangled states \[entangled mixtures\] of two noninteracting systems, say A and B. In a limiting case, the same can be said about all uncorrelated states of the form $`\widehat{\rho }_{A+B}=\widehat{\rho }_A\widehat{\rho }_B`$. If we require now that the uncorrelated states of noninteracting systems must evolve into uncorrelated states, it follows necessarily that the overall dynamical map $`g_{A+B}`$ \[which is linear\] must be a direct product of local dynamical maps, $`𝐠_{A+B}=𝐠_A𝐠_B`$. Moreover, each of the local maps $`𝐠_A`$ and $`𝐠_B`$ must be linear, by the same line of reasoning derived from Eq.(2). Since the linear map $`𝐠_A𝐠_B`$ must apply to entangled pure states as well, one must conclude that the corresponding reduced states also evolve linearly, whether they are probabilistic mixtures or not.

At this point one is led to reconsider the motivation behind the claim of ref.Gisin . The original observation was that any attempt ”to modify quantum physics, e.g., by introducing nonlinear evolution laws for pure states, …easily leads to the possibility of superluminal communication” \[Gisin , pg.1\]. From this it was reasoned that, perhaps, the converse may be true, and linearity is secured by the relativistic impossibility of superluminal signaling \[the no-signaling condition\]. As pointed out in ref.Bona , and emphasized in the preceding paragraph, this is not the case. Instead, the crucial conjecture leading to linearity is seen to be the postulated identification of mixed quantum states as probabilistic mixtures in the dynamical sense implicit in Eq.(2), or alternatively, the dynamical equivalence of probabilistic and elementary quantum mixtures. This equivalence proves sufficient to support the no-signaling condition, but the latter does not seem to be a prerequisite for linearity.

The term ”dynamical equivalence” is used here to emphasize a distinction from a mere ”static equivalence”. We reserve the latter term to mean that any maximal measurement on a given mixed state can be reproduced identically by a similar measurement on a unique ”probabilistic mixture” of orthogonal pure states, as well as on any member of an \[infinite\] family of statistically equivalent ”probabilistic mixtures” of nonorthogonal states. This is a purely statistical statement, which in itself does not imply a ”wavefunction collapse” during measurement, and thus need not involve the ”projection postulate”. In contrast, a ”dynamical equivalence” states that a ”statically equivalent” probabilistic mixture at time $`t`$ is equivalent to the evolved by pure state propagation of the ”statically equivalent” probabilistic mixture from an earlier time $`t_0<t`$. Note that a ”static equivalence” does not require a priori a ”dynamical equivalence”. The reason is simply that the physical behavior of a system is not defined only by its statistics, but also by interactions and, perhaps, self-interactions. In particular, the ”static equivalence” of mixed states to probabilistic mixtures cannot assert the absence or presence of state-dependent self-interactions, and cannot rule out dynamical nonlinearity on its own. On the other hand, the ”dynamical equivalence” expressed by Eq.(2) testifies precisely to the absence of any such self-interactions.

If we accept this distinction between static and dynamic equivalence, we may resolve that a corrected version of the theorem proposed in ref.Gisin  can read:

Usual quantum statics and the dynamical equivalence of mixed quantum states to probabilistic mixtures implies a linear, completely positive quantum evolution law, which necessarily complies with the no-signaling condition.

This leaves us with the question why exactly nonlinear quantum evolution laws appear incompatible with the no-signaling condition nlin-ftl-1 ; nlin-ftl-2 ; nlin-ftl-3 ; Polchinski ; Mielnik . The answer remains essentially with the remote preparation procedure discussed in ref.Gisin . It is not necessarily nonlinearity itself that is incompatible with no-signaling, but the induced distinguishability of probabilistic and elementary mixtures set against intrinsic quantum nonlocality. That is, if elementary mixtures can be converted by action-at-a-distance into probabilistic mixtures, then any local means to distinguish between the two types of mixtures, such as a nonlinear dynamics, or the nonlinear observables discussed in refs.Polchinski ; Mielnik , necessarily opens the possibility of superluminal communication.

## .3 Nonlocality and nonlinear quantum dynamics

One may object that the concept of remote preparation leading to the above conclusion involves an implicit acceptance of the ”projection postulate” Bona , contrary to the claim of ref.Gisin  that a strict statistical interpretation of nonlocal correlations should suffice. The reason lies in the intrinsic interconnection between the ”probabilistic mixture” interpretation of mixed quantum states, in particular of the output states of quantum measurements, and the ”projection postulate”. For suppose again that two noninteracting systems $`A`$ and $`B`$, at rest relative to each other, are initially in an entangled pure state $`|\mathrm{\Psi }_{AB}`$. Their respective local states are necessarily ”elementary mixtures” that cannot be prepared locally. Let a maximal preparatory measurement be effected and recorded on system $`B`$ at time $`t_0`$ as observed in the common rest frame, and let system $`A`$ be subject to an observation, also immediately recorded, at an arbitrarily close time $`t_0+\delta t`$. Then if the state of system $`A`$ is to be identified at $`t_0+\delta t`$ as a genuine ”probabilistic mixture”, amenable at least in principle to local preparation, it follows that the remote preparation at $`t_0`$ must transform the joint state of $`A`$ and $`B`$ into a separable mixture, realizable as a statistical superposition of products of local states \[otherwise the claim that the state of $`A`$ is no longer an ”elementary mixture” cannot be upheld\]. In other words, the total state must undergo a disentanglement process. We note in passing that such a disentanglement transformation will leave in general the measured system $`B`$ entangled with a measuring device, in which case the total state of $`(A+B+\mathrm{measuring}\mathrm{device})`$ should become a separable state of $`(B+\mathrm{measuring}\mathrm{device})`$ and $`A`$. Now, if at $`t_0+\delta t`$ this total state of separable form is indeed a ”probabilistic mixture”, then by the very definition of the latter every pair of $`A`$ and $`B`$ \[$`+\mathrm{measuring}\mathrm{device}`$\] in an ensemble representative of this mixture must actually exist in one of the contributing product states, say $`\widehat{\rho }_A^{out}\widehat{\rho }_B^{out}`$. Let the joint recorded observations of $`A`$ and $`B`$ corroborate one such product state in each run of the experiment. This evidently implies that each particular sample of $`A`$ and $`B`$ must undergo an evolution from the initial, pre-measurement state $`|\mathrm{\Psi }_{AB}`$ at $`t_0`$ to some final, post-measurement state $`\widehat{\rho }_A^{out}\widehat{\rho }_B^{out}`$ at $`t_0+\delta t`$, occurring with a certain probability among a set of possible product states. System $`A`$, in particular, must evolve from a local state $`\widehat{\rho }_A^{in}=Tr_B\left[|\mathrm{\Psi }_{AB}\mathrm{\Psi }_{AB}|\right]`$ at $`t_0`$ into the local state $`\widehat{\rho }_A^{out}`$ at $`t_0+\delta t`$. But this is exactly the ”projection at-a-distance” process that was to be avoided. As usual with the ”projection at-a-distance” paradox, one can only notice that for sufficiently small $`\delta t`$ the recorded observations of $`A`$ and $`B`$ become space-like separated events, and in suitably chosen referentials the observation of $`A`$ can precede the preparatory measurement on $`B`$. Hence the claim of a causal relationship between the two cannot be given any empirical support, and the interpretation that a local measurement can produce an ”instantaneous disentanglement” of an entangled state into a ”probabilistic mixture” is rendered meaningless. It is, however, fair to say that a local measurement can produce a ”disentanglement” of the total density matrix into a separable form, even if the corresponding mixture cannot be characterized as ”probabilistic”.

It is worth pointing out that the above line of reasoning applies as well if we begin with the assumption that local measurements produce local ”probabilistic mixtures”. Suppose that the preparatory measurement by a suitable measuring device leaves system $`B`$ in a ”probabilistic mixture” of local states. Then the total state must be a separable mixture of states of $`B`$ and joint states of $`A`$ and of the measuring device. But since system $`A`$ and the measuring device at the location of $`B`$ were initially in mutually uncorrelated states, and a local operation cannot produce nonlocal entanglement, it follows that the total post-measurement state must be separable with respect to all contributing parties. That is, it must be a statistical superposition of products of local states of $`A`$, $`B`$ and the measuring device, respectively. Since this separable form corresponds to a ”probabilistic mixture” for $`B`$, it must necessarily correspond also to a ”probabilistic mixture” for both $`A`$ and the measuring device. Obviously, we have just recovered the remote preparation process under discussion, and the rest of the above argument applies identically.

It becomes evident now that the statement ’ local measurements produce local and/or remote ”probabilistic mixtures”’ is equivalent to the ”projection postulate”. According to this logic, we are presented with the following orthogonal alternatives:

(L) We accept the ”probabilistic mixture” interpretation of output states of quantum measurements, and thus we implicitly accept the ”projection postulate”, with all inherent difficulties. On the other hand, the ”projection postulate” validates the phenomenon of remote preparation or ”projection at-a-distance”, which in turn can be invoked to argue that quantum dynamics is necessarily linear.

(NL) We circumvent the interpretation of quantum measurement outputs as ”probabilistic mixtures”, to the effect that mixed quantum states, including said measurement outputs, are to be treated as ”elementary mixtures”. In this case the concept of ”remote preparation” is rendered untenable, and although ”elementary mixtures” are statically equivalent to corresponding ”probabilistic mixtures”, there no longer is definite support for their dynamical equivalence. Hence a linear quantum dynamics is no longer indispensable.

The first option (L) has been reviewed in the preceding Section. For the sake of a closed argument, we find it instructive to explore here the question whether the second point of view (NL) can somehow yield viable nonlinear dynamical alternatives. Aside from a redefinition of mixed quantum states as elementary mixtures, we have in mind only an extension of the dynamical law, while all other ”static” notions \[including observables\] are retained in their usual form.

In this case, we must acknowledge that there already exists considerable evidence, both experimental tests-1 ; tests-2 ; tests-3 ; tests-4 ; tests-5 ; tests-6  and theoretical Mielnik , against a nonlinear dynamics of pure states in closed \[isolated\] systems. In truth, the past decade has brought us exceptional experimental support tests-1 ; tests-2 ; tests-3 ; tests-4 ; tests-5 ; tests-6  in favor of a linear \[and unitary\] propagation of pure states. Let us accept this fact as an empirical constraint on all acceptable quantum dynamics, which we cast in the form of a pure state condition: all pure states in the Hilbert space of states of a closed \[isolated\] quantum system evolve in time according to a linear dynamics.

The immediate question we face now is: once accepted the pure state condition, does it make any sense whatsoever to consider yet the possibility of a nonlinear dynamics, or is the problem already closed?

Counter-intuitive as may seem to the contemporary physicist heavily trained in linear dynamics, the fact that pure states must propagate linearly does not preclude, in itself, an intrinsically nonlinear dynamics for mixed states, and so does not imply Eq.(2) once decided that mixed quantum states must be regarded as ”elementary mixtures”. In fact, the literature already offers a number of positive, trace preserving nonlinear dynamics that reduce to a simple unitary evolution on all pure states of the state space Jordan ; Czachor-1 ; Czachor-2 ; Czachor-3 ; Czachor-4 ; Czachor-5 ; Czachor-6 ; Czachor-7 ; Beretta-1 ; Beretta-2 ; EQD .

A simple such example is provided by the nonlinear von Neumann equation Czachor-3

$$i\mathrm{}\dot{\widehat{\rho }}=[H\left[\frac{\widehat{\rho }}{Tr\widehat{\rho }}\right]^q+\left[\frac{\widehat{\rho }}{Tr\widehat{\rho }}\right]^qH,\widehat{\rho }],$$
(3)
with $`q>0`$ a real scalar, and $`H`$ a self-adjoint operator. This dynamics produces an evolution of the form

$$\widehat{\rho }(t)=S_{\widehat{\rho }}(t)\widehat{\rho }(0)\left[S_{\widehat{\rho }}(t)\right]^{},$$
(4)
where the propagator $`S_{\widehat{\rho }}(t)`$ is a unitary operator with a nonlinear dependence on $`\widehat{\rho }`$, and satisfies the equation of motion

$$i\mathrm{}\dot{S}_{\widehat{\rho }}(t)=\left[H\left[\frac{\widehat{\rho }}{Tr\widehat{\rho }}\right]^q+\left[\frac{\widehat{\rho }}{Tr\widehat{\rho }}\right]^qH\right]S_{\widehat{\rho }}(t).$$
(5)
One can easily check that pure states evolve into pure states, since $`(d/dt)(\widehat{\rho }^2\widehat{\rho })\dot{\widehat{\rho }}\widehat{\rho }+\widehat{\rho }\dot{\widehat{\rho }}\dot{\widehat{\rho }}=0`$. Moreover, in this case \[$`\widehat{\rho }^2=\widehat{\rho }`$\] the nonlinearity disappears, and the von Neumann equation reduces to the familiar linear form

$$i\mathrm{}\dot{\widehat{\rho }}=[H,\widehat{\rho }].$$

The dynamics described by Eq.(4) becomes physically meaningful if, in addition, it complies with the no-signaling condition. In particular, the uncorrelated states of two noninteracting systems, $`A`$ and $`B`$, must propagate into uncorrelated states \[separability\], and entangled mixed states that produce identical local initial conditions must generate identical local evolutions \[locality\]. To reconcile these requirements, one may adapt Polchinski’s conjecture Polchinski  \[see also Czachor-1 ; Czachor-2 ; Czachor-Doebner \], and postulate the total propagator as the product of local propagators for the local reduced states, i.e.,

$$S_{\widehat{\rho }_{AB}}^{(AB)}(t)=S_{Tr_B\widehat{\rho }_{AB}}^{(A)}(t)S_{Tr_A\widehat{\rho }_{AB}}^{(B)}(t).$$

Accordingly, the total generator must read

$$G^{(AB)}\left(\widehat{\rho }_{AB}\right)=G^{(A)}\left(Tr_B\widehat{\rho }_{AB}\right)+G^{(B)}\left(Tr_A\widehat{\rho }_{AB}\right).$$
(6)
Despite its artificial aspect, this definition finds a natural justification if the density matrix in the original equation of motion (3) is interpreted from the beginning as a local state obtained by averaging out possible entanglement with a noninteracting environment \[$`\widehat{\rho }\widehat{\rho }_{loc}=Tr_{env}[\widehat{\rho }_{loc+env}]`$\]. It was argued Czachor-Doebner  that such well-behaved nonlinear dynamics may indeed be formulated in a manner that avoids both the ”wavefunction collapse” paradigm of quantum measurement theory, and a conflict with quantum nonlocality. In a separate paper nlin-th  we address the problem of an extended framework for the formulation and analysis of these theories .

Here we remark that as a direct corollary of the above observation the problem of linearity in quantum dynamics cannot be considered closed in absence of a test of linear propagation for mixed states. Although this task may seem considerably more demanding than a test of pure state dynamics, we wish to point out that it actually amounts to a test of linearity for the dynamics of entangled pure states. Indeed, the following simple argument shows that no nonlinear theory may account for the linear propagation of both the pure states of isolated systems and the entangled pure states of noninteracting systems without a conflict with the no-signaling condition. Consider again two noninteracting systems $`A`$ and $`B`$, and suppose that all pure states of the joint system $`(A+B)`$ evolve linearly, under a common propagator. By the same reasoning used previously for the case of probabilistic mixtures, uncorrelated pure states must evolve into uncorrelated pure states and the overall linear propagator must be a direct product of individual linear propagators for $`A`$ and $`B`$. If one considers now an arbitrary entangled state, it follows that the corresponding reduced states must also evolve linearly. On the other hand, there are an infinity of entangled mixed states that produce the same reduced state for either one of the systems. But if entanglement is not to be distilled locally \[the no-signaling condition\], a given reduced state must evolve in the same manner, regardless of any entanglement with a remote counterpart. Hence the reduced states generated by entangled mixed states must also evolve linearly, which means in effect that the overall dynamics must be linear.

Conversely, any nonlinear dynamics that is well-behaved with respect to the no-signaling condition, and also conforms to the pure state condition, must necessarily predict that pure entangled states of noninteracting systems must evolve, and perhaps decohere into mixed states, in a nonlinear manner. We are led in this way to the linearity criterion stated in the opening paragraph, and which we recall here for convenience:

Linearity criterion: If all pure states of nonseparable, isolated systems evolve linearly, then the overall quantum dynamics is linear if and only if the entangled pure states of noninteracting systems evolve linearly.

This is our basic result concerning the previously defined alternative (NL), under the constraint of the pure state condition.

## .4 Summary

We showed that, contrary to the assertion of ref.Gisin , the relativistic ”no-signaling condition” is not necessary for a rationalization of linear quantum dynamics from ”quantum statics”. It turns out that the ”probabilistic mixture” interpretation of mixed quantum states in its dynamical form is already sufficient for this purpose. This is also the essential obstacle faced by nonlinear quantum theories, and not the incompatibility between a nonlinear quantum dynamics and relativity, as advocated lately. Perhaps it may be argued that the Hilbert space structure of the quantum state space is determined by the relativistic structure of space-time Svetlichny . But once this structure is accepted as the fundamental component of quantum statics, and is supplemented by the ”probabilistic mixture” point of view, there is no need for a recurring reference to relativity in order to justify the linearity of quantum dynamics.

On the other hand, the ”probabilistic mixture” interpretation of mixed states implies a covert acceptance of the ”projection postulate”. If one chooses to dismiss ”projections-at-a-distance”, e.g., as incompatible with relativity, then mixed quantum states must be regarded exclusively as ”elementary mixtures”. A suitable physical argument for this point of view may be, for instance, that there always exists a finite chance of nonvanishing, nondistillable entanglement with noninteracting components of the environment. In this case previous arguments Gisin ; nlin-ftl-1 ; nlin-ftl-2 ; nlin-ftl-3 ; Polchinski  that rule out dynamical nonlinearity by conflict with the dynamical ”probabilistic mixture” interpretation under the ”no-signaling condition” loose their ground, and well-behaved, relativistically consistent nonlinear dynamics become possible Jordan ; Czachor-1 ; Czachor-2 ; Czachor-3 ; Czachor-4 ; Czachor-5 ; Czachor-6 ; Czachor-7 .

Certainly, the viability of these theoretical alternatives is essentially limited by empirical evidence. We have given particular attention to the very likely empirical constraint tests-1 ; tests-2 ; tests-3 ; tests-4 ; tests-5 ; tests-6  that all pure states of isolated systems must evolve unitarily into pure states, as in linear quantum theory. An important result concerning this possibility is a simple linearity criterion, which suggests that a decisive empirical answer to the (non)linearity problem in quantum dynamics can be provided in two steps:

i) a test of linear \[and unitary\] propagation of pure states in closed systems, which is already available tests-1 ; tests-2 ; tests-3 ; tests-4 ; tests-5 ; tests-6 ;

ii) a test of linear \[and unitary\] propagation of entangled pure states of noninteracting systems.


text:
# Design of UHECR telescope with 1 arcmin resolution and 50 degree field of view

## 1 Introduction

The observation of ultra high-energy cosmic rays (UHECRs) and high-energy neutrinos (HE-$`\nu `$’s) with directional accuracy of $`1`$ arcmin (’) is expected to lead to revolutionary progress in particle astrophysics through “multi-particle” observations in cooperation with multi-wavelength astronomical observations. An excellent example is the recent prompt observation of X-rays from a gamma-ray burst (GRB) with directional precision of 1’ by Beppo-SAX, an event that marked the start of multi-wavelength observations using many advanced astronomical telescopes . The directional accuracy of 1’ is several tens times better than that of $`1^{}`$ currently operating, constructing, or planning UHECR detectors  aim for.

Upon entering the atmosphere, UHECRs produce an air shower (AS) through interaction with air nuclei. The AS axis coincides with the UHECR incident angle, and the intrinsic angular resolution of the shower is 0.3’ or smaller for a UHECR of 1017 eV. This narrow dispersion results from the high Lorenz factor of secondary particles ($`10^4`$) in the direction of the parent particle, and the very low transverse momenta of secondary particles (rarely exceeding 1 GeV/c) as reported by the ATLAS group . For UHECRs, more precise reconstruction of the primary direction can be expected due to the higher Lorenz factor.

ASs can be detected via the isotropic fluorescence emission of the air nitrogen that it excites. An optical system with imaging devices can be used to track the longitudinal development of an AS. This technique was first employed in the Fly’s Eye detector  and will be part of several future projects on UHECR detection. The key advantage of fluorescence detection is the capacity for achieving fine angular resolution. The signal to noise ratio ($`S/N`$) for each pixel of the imaging devices is inversely proportional to the square root of the angular diameter of the field of view (FOV). A better $`S/N`$ will allow the device to trigger charge levels, making it possible to analyze less energetic or more distant events.

Galactic and extragalactic magnetic fields are considered to have intensities on the order of $`\mu `$G, nG, and below . Thus, UHECRs are expected to experience only slight deflection due to galactic and extragalactic magnetic fields, and cannot be confined by the galactic magnetic field. The existence of an “ankle” at $`E10^{19}`$ eV is often interpreted as a crossover from a steeper galactic component to a harder component of extragalactic origin.

As the possible sources of astrophysical accelerators capable of producing UHECRs are very limited, UHECRs that have been deflected only slightly can be traced back to their sources. Thus, high-resolution measurements of the arrival direction will lead to significant progress in cosmic-ray astronomy. The magnetic deflection of UHECRs during propagation induces charge separation. One of the important implications of 1’ angular resolution is that this separation of charge can be resolved, allowing unambiguous particle identification to be performed based on the deflection from sources. Such observations may also allow us to study the structure of extragalactic magnetic fields.

Based on the above considerations, we propose a telescope design that features a FOV of 50  $`\times `$ 50  and 1’ resolution for the arrival direction. This is realized through a combination of a large-area image intensifier (II) connected to a CCD camera  and a customized Baker-Nunn optical system. In this paper, we describe the conceptual design of the modified Baker-Nunn optical system. The expected performance of the telescopes for UHECR detection is also presented assuming that the image resolution of the CCD camera is better than that of the optical system.

## 2 General Concept of Optical System

Our design requirements can be summarized as follows:

1. a large FOV of 50 $`\times `$ 50 to cover a large solid angle for UHECR detection;
2. an image resolution of $``$0.02 to achieve 1’ resolution for the arrival direction;
3. a large caliber with a diameter of $``$2000 mm to distinguish faint fluorescence.

In addition, cost performance will be a critical issue with respect to construction of a large-aperture array of telescopes.

Among several approaches to meet our requirements, spherical-mirror optics with an aperture diaphragm on the mirror’s curvature center appears to offer some advantages for the proposed telescope. This simple optics produces homogeneous images within the entire FOV. Typical off-axis aberrations, principally coma aberration, are eliminated, although spherical aberration is a factor. Thus, from a review of the variants of this optics, the Schmidt, Makstov, Baker Super-Schmidt, and Baker-Nunn optics were selected as candidates. The Schmidt optics has a one-element corrector lens at the diaphragm, and provides higher resolution than optics without the corrector lens. However, we found that it was difficult to achieve the required image resolution over the wide FOV. The Makstov and Baker Super-Schmidt optical systems require thick corrector lenses, and as such are not suited to large-scale telescopes.

The Baker-Nunn optics has a three-element corrector lens consisting of convex and concave lenses with spherical or aspherical surfaces . The basic premise of this optics is that each element has easy-to-make spherical surfaces, residual spherical aberration after assembly can be canceled by adjusting the spacing between lens elements, and chromatic aberration can be eliminated by selecting appropriate lens materials. Figure 1 shows a schematic of the Baker-Nunn optics. The second lens element is an aperture diaphragm.

Although the advantages of the Baker-Nunn optics are not related directly to the specific area of improvement pursued in the proposed design, the three-element corrector lens in this optics provides a large degree of freedom for optimization. As a result, it is considered the most suitable for our purposes. In our design, the first and third lens elements are 2400 mm in diameter, and the second element is 2000 mm in diameter, and all elements are 10 mm thick. Considering the cost, and complexity of fabrication processes, we designed the corrector lens as an acrylic Fresnel lens as appropriate for a large-aperture telescope. As described in the following section, the surface shape of each lens element was determined through an optimization process.

## 3 Optimization by Ray Tracing

### 3.1 Ray Tracing

We developed ray-tracing program written in C++ in order to optimize the optical system and examine its performance. The program models the refraction and reflection of incident rays on each surface, and allows a bundle of on- and off-axis rays to be traced numerically. As is common in designing an optical system, the shape of the front surface of each optical component is expressed as

$$Z=\frac{h^2/R}{1+\sqrt{1h^2/R^2}}+Ah^4+Bh^6+Ch^8,$$
(1)
where $`Z`$, $`h`$, and $`R`$ represent the surface sag, radial height, and curvature of the surface. $`Z`$ is positive in the direction of on-axis incident rays. $`A`$, $`B`$, and $`C`$ represent the 4th, 6th, and 8th order aspheric coefficients. The surfaces of each lens element and the mirror are numbered 1 to 7 as shown in Fig. 1.

The focal surface is defined as the surface minimizing the image spot for a bundle of parallel rays. The spot size is defined as the root mean square (RMS) of positions of parallel rays at the focal surface. The fineness of the bundle of rays analyzed is determined according to the purpose of ray tracing: a rather rough interval between rays is selected for optimization, and a fine interval is selected for performance evaluation. The spot size ($`s`$) is a function of the incident angle ($`\theta `$) and the wavelength ($`\lambda `$) of parallel rays.

### 3.2 Optimization

We used the modified Powell’s method as the algorithm for optimization based on the following evaluation function:

$$F=\underset{\theta }{}(\underset{\lambda }{}s(\theta ,\lambda )^2)^2,$$

where $`\theta `$ varies from 0 to 25, and $`\lambda `$ varies from 330 to 410 nm, taking into account the air-fluorescence emission spectrum. We adopted the fourth power of spot size as the evaluation function in order to obtain flat resolution over the entire FOV.

In principle, the curvature, aspheric coefficients, thickness, and position of each lens element are free parameters in optimization. In order to reduce the number of free parameters, we imposed the symmetry described in Ref.  in which the corrector lens is symmetric with respect to the plane including the center of the second element and which is perpendicular to the axis of the ray. In addition, the thickness and position of each surface were fixed. These constraints allowed much more efficient optimization. The optimized parameters of the optical system are given in Table 1. The axis and origin of the $`Z`$ coordinate are coincident with the axis of the rays and the center of the corrector lens. Note that we took the effect of the region obscured by the focal surface into account in the optimization.

Figure 2 shows the surface sag of lenses 2, and 3 as a function of radial height. The surface sag remained well below 0.5 mm when the width of Fresnel grooves was chosen to be 10 mm. Thus, the 10-mm thickness of each lens element is sufficient, although the thickness may be changed according to the support method employed for the corrector lens. We confirmed that the same level of spot-size resolution can be obtained by another optimization using a different lens thickness.

Figure 3 shows the spot diagrams of the optimized optical system. The shape and size of spot diagrams changed with the wavelength and the incident angle of the incoming rays. The four strongest fluorescence peaks, $`\lambda =`$ 336, 358, 391, and 425 nm, which dominate the filtered fluorescence spectrum after propagation through air, were selected as wavelengths for the diagrams. The spot diagrams for incident angles of 5, 15, and 25 are shown.

Figure 4 shows the spot size distribution as a function of incident angle, taking the characteristic wavelengths in the fluorescence spectrum into account. Solid and dashed curves represent the spot size optimized to cover a FOV of $`\pm `$25 and $`\pm `$20, respectively. The lens parameters shown in Table 1 are the result of the former optimization, and give a spot size of smaller than 0.02 over a FOV of $`\pm `$25. The results were confirmed by cross checks using Mathematica Optica  with correction for the refractive index of the atmosphere.

In parallel with the development of this optical system, the present authors have also developed a high-density CCD imaging detector with integrated large-area image intensifier, which has been demonstrated to have a resolution of $``$ 0.02 . This imaging system is mounted at the focal surface of the optimized optics to give the proposed telescope design.

## 4 Performance

A solid angle of nearly 2$`\pi `$ sr can be observed using an array of 8 telescopes (i.e., a station), with 7 telescopes arranged azimuthally and 1 telescope directed upward. It should be noted that the telescope directed upward is primarily important in HE-$`\nu `$ detection rather than UHECR detection .

A Monte Carlo simulation was performed to examine the resolution obtained for the arrival direction of UHECRs. The event generator calculated the signal strength produced by atmospheric fluorescence emission from ASs with a given geometry, and built simulation events by calculating the number of photoelectrons detected by each pixel of a CCD camera at each mirror. The simulation includes the longitudinal and lateral profiles of the AS, the model of the atmosphere, the air fluorescence yield, and the detector optics. We also took into account fluctuation of the first interaction depth, impact point, and directional angles of shower cores assuming appropriate distributions.

Current and planned fluorescence detectors  have a pixel resolution of $``$1. The angular resolution of the reconstructed track was evaluated by analyzing event samples that triggered the current 1-pixel detectors. Such signals are triggered by having at least 6 pixels ($`1^{}`$ resolution) fire at 4$`\sigma `$ greater than the night sky background, with a corresponding reconstructed track extending more than 5 in the FOV . A night sky background level of 100 photoelectrons/$`\mu `$s/($`1^{}\times 1^{}`$) was assumed. Figure 5 shows the comparison of event displays between the 0.02 and 1 systems.

Applying a stereo reconstruction of the AS, the arrival direction is determined by the geometrical intersection of shower-detector planes (SDPs). Since the resolution of the normal vector of a SDP is the primary factor affecting the resolution of arrival direction, we discuss here for simplicity the accuracy of SDP measurement. The angular resolution of SDP is plotted in Fig. 6 for 1019 and 1020 eV protons. The angular resolution is determined from the integral of the distribution encompassing 68 % of events. Figure 7 shows the variation in the angular resolution of the SDP with proton energy. Resolution of $`<`$ 1’ is obtained at above 1018 eV for the SDP. In the case of stereo reconstruction with two or more stations, the pointing resolution is also expected to be $``$ 1’ after reconstruction of the geometrical intersection of SDPs.

The $`S/N`$ of each pixel is proportional to $`(A/\mathrm{\Delta }\theta )^{\frac{1}{2}}`$ , where $`A`$ represents the light collection area, corresponding to the reflection area for simple mirror optics and the area of the aperture diaphragm for the Baker-Nunn optics, and $`\mathrm{\Delta }\theta `$ denotes the angular size of each pixel. Note that this proportionality holds only when the effects of lateral distribution are neglected. Assuming the same size of $`A`$, the $`S/N`$ value for the 0.02-pixel imaging device is expected to be about 7 times higher than that from the 1-pixel one. However, this is expected to be affected by the quantized number of photoelectrons at these small dimensions. To examine this, the AS profile was simulated by the Monte Calro method and the quantitative improvement in the $`S/N`$ was examined for the 0.02 pixels. The AS events of protons with energy of 1020 eV were generated isotropically for a given impact parameter ($`R_\mathrm{P}`$). The maximum value of $`S/N`$ (obtained event-by-event) was averaged for a given $`R_\mathrm{P}`$, and the maximum values are shown in Fig. 8. Note that the event selection applied here differs from that applied to compare the angular resolution as it is not necessary to apply the 6 pixel condition to compare the $`S/N`$ of the $`1^{}`$\- and 0.02-pixel detectors for distant ASs. The solid curve shows the expected $`S/N`$ value for 1 pixels calculated from that of 0.02 pixels by the theoretical proportionality. The expected and obtained $`S/N`$ values are in good agreement for the distant events ($`R_\mathrm{P}>`$ 70 km). For the closer events, the lateral distribution of the shower-core image becomes significant compared with that for the 1 pixel image (see the bottom panel of Fig. 5). As a result, the expected $`S/N`$ ratio is appreciably lower than the simulation result.

We also define another signal to noise ratio ($`S/N_{\mathrm{TRK}}`$) appropriate for the track captured by the imaging devices. The signal ($`S`$) is the sum of charges detected by all pixels associated with the track. The background ($`BG`$) is the sum of background for all pixels associated with the track. Then, $`S/N_{\mathrm{TRK}}`$ is defined as $`S/N_{\mathrm{TRK}}=S/\sqrt{BG}`$. Figure 9 shows the distributions of $`S/N_{\mathrm{TRK}}`$ for protons of 1019 and 1020 eV. The 0.02 device exhibited excellent $`S/N_{\mathrm{TRK}}`$ values. The significant improvement in $`S/N_{\mathrm{TRK}}`$ was obtained because the lateral extent of the AS is smaller than the pixel resolution of 1. Thus, we have demonstrated the improved resolution of the arrival direction offered by the proposed system. The detailed performance of the detector, including aperture, energy resolution, and shower-max resolution, will be discussed elsewhere.

## 5 Conclusion

We have presented the conceptual design of a new type of telescope for UHECR detection based on the Baker-Nunn optical system. The telescope has a spot size resolution of smaller than 0.02 within the entire FOV of 50, and has been shown to be capable of 1’ angular resolution for the arrival direction of UHECRs. This design has the potential to open new possibilities in the field of UHECR and HE-$`\nu `$ physics through high-resolution observation.

The authors express their gratitude to Y. Higashi for helpful advice and fruitful discussion, and to T. Aoki for continuous support. We also thank N. Manago and M. Jobashi for their useful help in the development of the Monte Carlo simulation.


text:
# Lyman-break galaxies: are they young spheroids?

## 1 Introduction

Since the discovery of a large population of actively star-forming galaxies at $`3z3.5`$, identified by their red-shifted Lyman continuum breaks at 912Å in the rest-frame (hereafter LBGs, see Steidel et al. 1996a,b), it has became clear that studying their properties as a whole, as well as their individual spectra, could allow us to obtain information on star formation histories, ages and metallicities of such high redshift objects. This information can then be used to compare the LBGs to present-day galaxies. ¿From Hubble Deep Field and ground-based far-UV spectra, Steidel et al. (1996a,b) found that these galaxies are very similar to nearby star-forming galaxies. The observation of space density, star formation rate, masses, physical sizes and morphologies, led them to the conclusion that the LBGs they discovered could be the progenitors of the spheroidal components of present-day luminous galaxies (i.e. elliptical galaxies or spiral bulges). An argument in favor of ellipticals in clusters as the descendants of LBGs is also given by the strong clustering observed for these objects (Adelberger et al. 1998; Steidel et al. 1998). On the other hand, other authors (Lowenthal et al. 1997; Sawicki and Yee, 1998; Sommerville, Primack & Faber, 2001) favor the idea that LBGs can be sub-$`L^{}`$ galaxies or that they are low mass starbursting objects which later merge to form more luminous galaxies. A crucial test for the nature of LBGs is to measure their masses. The inferred stellar masses for LBGs are in the range $`M_{}=10^910^{11}M_{}`$ with an average mass of $`10^{10}M_{}`$ (Papovich et al. 2001). However, these are only lower limits since they measure the young starburst component.The inferred ages for the starbursting populations range from 30 Myr to 1 Gyr. Again, these are lower limits to the ages of LBGs, since it cannot be excluded that star formation was occurring even before the observed starburst.

By means of near-IR observations, Shapley et al. (2001) studied 81, $`z3`$ LBGs and found a strong correlation between the inferred star-formation rates and dust extinction: younger galaxies have more dust and higher star-formation rates than older galaxies. Dustier galaxies are also more luminous. To explain these correlations they proposed an “evolutionary sequence” for LBGs, in which the younger galaxies evolve into the older, less reddened and quiescent ones. They derived ages $`t_{sf}`$ for LBGs ranging from several Myr to 1 Gyr, and best-fit star-formation rates in the range $`\psi (t_{sf})5940\mathrm{h}^2\mathrm{M}_{}\mathrm{yr}^1`$ with a median $`\psi (t_{sf})45\mathrm{h}^2\mathrm{M}_{}\mathrm{yr}^1`$. The total masses in stars derived by integrating the average star formation rate over time are $`10^{10}M_{}`$. They detected also evidences of winds and outflows driven by SNe and massive stars, which could enrich the intergalactic medium in metals and energy.

In the following we focus on MS 1512-cB58, the brightest LBG known so far owing to its gravitationally lensed nature. This object is at $`z=2.7276`$ and has a luminous mass of $`10^{10}M_{}`$, a star formation rate of $`\psi (t_{sf})40\mathrm{M}_{}\mathrm{yr}^1`$ (Pettini et al. 2001) and an effective radius of $`r_L2\mathrm{kpc}`$ (Seitz et al. 1998), for a $`\mathrm{\Omega }_m=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,h=0.70`$ cosmology. From the study of MS 1512-cB58 Pettini et al. (2001) derived also the interstellar medium (ISM) chemical abundances. In particular, they concluded that the abundances of O, Mg, Si, P, S are $`2/5`$ of their solar values, whereas N, Mn, Fe and Ni are underabundant by a factor of $`3`$. Depletion into dust, which is known to be present in MS 1512-cB58, probably accounts for some of the Fe-peak element underabundances, but this is not likely to be an important effect for nitrogen. They indicate that MS 1512-cB58 is young with an age, less than the typical timescale for the release of N from intermediate mass stars and Fe from type Ia SNe. They suggest $`t_{sf}300\mathrm{M}\mathrm{y}\mathrm{r}`$ while Ellingson et al. (1996) derived an age younger than $`35`$ Myr for the starburst in MS 1512-cB58 by means of a fit of the spectral energy distribution (SED). Finally, an outflow involving the bulk of the ISM and proceeding at a speed of $`255\mathrm{Kmsec}^1`$ seems to be present in MS 1512-cB58. This velocity probably exceeds the escape velocity of the galaxy, thus suggesting that MS 1512-cB58 is probably suffering a galactic wind, although firm conclusions on this point are not possible.

Models for the chemical evolution of elliptical galaxies, which are able to reproduce the observed properties of local ellipticals, assume that most of the star formation in these objects took place several billion years ago and lasted for a period no longer than 1 Gyr, at the end of which a strong galactic outflow occurred. These kinds of models are known as “supernova-driven galactic wind models” for ellipticals and were first introduced by Larson (1974). Supernova-driven wind models including detailed chemical evolution and supernova energetics have been later proposed by Arimoto & Yoshii (1987) and Matteucci & Tornambè (1987) and they were aimed at explaining the mass-metallicity relation and the color-magnitude diagram of ellipticals. A multi-zone, more-refined version of these models was then presented by Martinelli, Matteucci & Colafrancesco (1988). These models are also called “monolithic models” as opposed to “hierarchical models” (see e.g. Kauffmann, 1996), the two main competing pictures for galaxy formation. The main difference between these two scenarios of galaxy formation is that in the first case ellipticals formed in a short timescale (1-2 Gyr) and at very high redshifts (z $``$ 3) and evolved passively afterwards, whereas in the second case galaxies assembled in a hierarchical fashion. In particular, the hierarchical process took place over a large time interval with the consequence of having active star formation until recently in large ellipticals.

The aim of this letter is to test whether the observed chemical properties of LBGs, and in particular of MS 1512-cB58, can be explained under the assumption that LBGs are normal elliptical galaxies forming at high redshift, and to set independent constraints on the age of this galaxy. In particular, abundance ratios, especially those between elements formed on different timescales, can set important constraints on the age and the nature of high-redshift objects. In Section 2 we describe the chemical evolution model we adopt, in Section 3 we present the model results and their comparison with the observed abundances. Finally in Section 4 some conclusions are drawn.

## 2 The Model

The elliptical galaxy model we adopted here is a multi-zone model (Martinelli, Matteucci & Colafrancesco, 1998; Pipino et al. 2002), in which the gas collapses into the potential well of a dark matter halo at high redshifts ($`z3`$) on free-fall timescales. The star formation rate is very efficient and lasts until the energy injected into the ISM by stellar winds and supernovae (Ia,b and II) triggers a galactic wind which develops first in the most external regions and later in the more internal ones. As a consequence, while the outer regions are suffering an outflow the inner ones are still actively forming stars. This kind of model can account for the abundance and color gradients observed in ellipticals (Menanteau, Jimenez & Matteucci, 2001). Generally, after the onset of the galactic wind no star formation is assumed to take place in that particular region, thus supernovae (SNe) of type II and Ib, which originate from massive stars, no longer occur whereas SN Ia continue to explode until the present time. In fact, type Ia SNe are believed to originate from white dwarfs in binary systems. Therefore, the evolution of the galaxy, after the onset of the wind, is dominated by type Ia SNe which inject into the ISM heavy elements (mostly iron) and thermal energy.

For the equations of chemical evolution and related details, we direct the reader to the above mentioned papers, while here we recall only the main assumptions of the model. The initial mass function (IMF) is expressed by the classic Salpeter (1955) law: $`\phi (m)m^{(1+1.35)}`$, defined in the mass range $`0.1100M_{}`$. Given the uncertainties in the knowledge of the mass of this object we assumed several values for the initial luminous mass (gas plus stars) of the galaxy. In this letter we present three cases, with initial luminous masses $`10^9`$, $`210^{10}`$ and $`10^{11}M_{}`$ and effective radii of 0.5, 2 and 3 kpc (see Matteucci, 1992), respectively. Each galaxy is divided in several concentric shell roughly 0.5 kpc wide.

The star formation rate (SFR) $`\psi (t)`$ is related to the gas mass by the following equation:

$$\psi (t)=\nu M_{gas}(t)$$
(1)
with the efficiency of star formation $`\nu `$ being in the range 15–12 $`Gyr^1`$ (see Pipino et al. 2002). In this formulation the star formation rate is maximum at the beginning and then decreases with time, as shown in Figure 1 (where the average star formation rate inside the optical radius of each galaxy is plotted), and strongly depends upon the galaxy mass. The potential energy of the gas is computed by assuming that the galaxy is surrounded by a massive and extended dark matter halo with $`M_{dark}=10M_{lum}`$ and $`\frac{r_L}{r_D}=0.1`$, where $`r_D`$ represents the radius of the dark matter core (Bertin et al., 1992). The onset of the galactic wind in each zone depends on the time at which the energy, restored by SNe, exceeds the gas binding energy. At this time, $`t_{GW}`$, all the gas present in the galaxy starts to be lost. Clearly $`t_{GW}`$ is extremely sensitive to the assumptions made about the feedback between SNe and ISM. The detailed calculation of the ISM thermal energy due to stellar winds and SNe can be found in Pipino et al. (2002). In that paper we adopted new energetic prescriptions for SNe taking into account the most recent results concerning the evolution of a SN remnant in the ISM. In particular, the most crucial parameter in order to compute the amount of energy transferred by a single SN into the ISM is the cooling time of the remnant. This is the time at which radiative losses cannot be neglected and most of the energy deposited by the shock wave into the ISM, during the adiabatic phase, is lost radiatively. In particular, for type II supernova remnants we adopted the Cioffi, Mc Kee & Bertschinger (1988) cooling time, which is a function of density and metallicity of the ISM. On the other hand, we assumed that SN Ia can transfer all of their initial blast wave energy into the ISM. In fact, since type Ia SNe originate from long living systems, their explosions occur in a medium already heated by SNII and therefore the cooling is neglible (see Recchi et al. 2001). Under these assumptions the $`210^{10}M_{}`$ model develops a galactic wind in the most external region (roughly 3 kpc from the center) at $`31\mathrm{M}\mathrm{y}\mathrm{r}`$ and subsequently in the inner region (0-1 kpc) at $`200\mathrm{M}\mathrm{y}\mathrm{r}`$ after the major episode of star formation has started. This situation is certainly compatible with what observed in MS 1512-cB58, where the outflowing gas is probably located in front of the star forming region at a distance of a few kpc (Pettini et al. 2001). The amount of stars formed at the time of the onset of the external wind is $`1.110^{10}M_{}`$ all over the galaxy with the stars being more concentrated towards the center. The models for $`10^9`$ and $`10^{11}M_{}`$ do not seem to reproduce the properties of MS 1512-cB58 since they predict too low and too high star formation rates at any galactic time, respectively, as evident in Figure 1. They also develop a wind in the most external regions too early and too late, respectively.

By means of our model we can predict the evolution of the abundances in the gas of H, He, N, O, Mg, Si, Fe versus time and metallicity. We adopted different sets of stellar yields in different stellar mass ranges: in particular, for massive stars ($`M>8M_{}`$) the yields of Woosley & Weaver (1995), for low and intermediate mass stars ($`0.8M/M_{}8`$) the yields of Renzini & Voli (1981) and for type Ia SNe the yields of Nomoto et al. (1984). It is worth noting that the model described above has been already compared with the observational properties of local ellipticals (see Pipino et al. 2002 and references therein), and they reproduce both colors and stellar chemical abundances of normal ellipticals.

## 3 Model results

In Figure 1 we show the predicted and observed star formation rates for MS 1512-cB58. The value quoted by Pettini et al. (2001) is 40$`M_{}yr^1`$ but we can consider a range of possible values of 20–80 $`M_{}yr^1`$, given the uncertainties connected with the derivation of the star formation rate. From the comparison between the predictions of the $`210^{10}M_{}`$ model (the best model) and the observed value, we infer an age for the galaxy inside one $`r_L`$ of $`100\pm 60`$Myr. In Figures 2 and 3 we show the predicted \[N/Fe\], \[Mg/Fe\] and \[Si/Fe\] versus \[Fe/H\] and time, respectively, and compare them with the observed ones in MS 1512-cB58. The observed values have been corrected for possible dust depletion of Fe, although the underabundance of Fe-peak elements (relative to the undepleted S) cannot be entirely due to this effect. In fact, if that would be the case we should expect the abundances of Mg and Si to be also depleted, although to a lesser extent, relative to S. This is not what is observed since the ratios of Mg and Si relative to S are solar, as discussed by Pettini et al. (2001). We assumed that roughly half of the total Fe mass is hidden in grains (Pettini, private communication). Therefore, the observed Fe abundance has been increased by roughly a factor of 2. As Figure 2 shows, the agreement between model predictions and observations is good for the \[el/Fe\] vs. \[Fe/H\] relationships, except for Mg which is predicted to be lower than observed. However, this is expected since all the available yields from massive stars underestimate the solar Mg abundance (e.g. Thomas, Greggio & Bender, 1998). Unlike the star formation rate, the abundance ratios in the gas are quite insensitive to the galaxy mass, whereas they depend strongly upon the nucleosynthesis, The IMF and timescales for the production of the different elements.

In Figure 3, where the abundance ratios are plotted as functions of galactic age and redshift, we indicate the observational values by dotted lines and see where these lines intercept the predicted curves. The interception gives an indication of the possible age of MS 1512-cB58. For \[Mg/Fe\] and \[Si/Fe\] the estimated galactic age can be as low as $`10`$ Myr, whereas for the \[N/Fe\] ratio the inferred age is $`35`$ Myr, in excellent agreement with other independent estimates (Ellingson et al. 1996). From the \[$`\alpha `$/Fe\] ratios we cannot exclude ages larger than 10 Myr but we can exclude ages larger than 300 Myr since the ratios at these times are considerably lower than observed. The high values of the \[$`\alpha `$/Fe\] ratios indicate, in fact, that the object is younger than the typical timescale for the enrichment from SNe Ia (0.3-0.5 Gyr for ellipticals, Matteucci & Recchi, 2001). Nitrogen is also a good cosmic clock since it is mainly produced by low and intermediate mass stars, although massive stars contribute as well for a small amount of this element. As a consequence, N starts to be produced in a non neglible way only after 35 Myr (the lifetime of a 8 $`M_{}`$ star) while its bulk appears only after 300 Myr (the lifetime of a 2.7 $`M_{}`$). The ratio of N relative to any $`\alpha `$-element is even more sensitive to age then the \[$`\alpha `$/Fe\] ratio and this is because N is mainly a secondary element and therefore sensitive to the overall metallicity. The age-\[Fe/H\] relationship predicted by the model for the $`210^{10}M_{}`$ galaxy suggests that \[Fe/H\]=-0.85$`\pm 0.1`$ dex (the dust-corrected iron abundance) is reached at an age between 20 and 35 Myr. It is worth noting that we obtain a good agreement with observations also for the predicted \[N/O\] ratio which is independent of dust depletion. In fact, we considered the oxygen abundance measured for MS 1512-cB58 by Teplitz et al. (2000) (\[O/H\]=-0.35 dex) together with the N abundance by Pettini et al. (2001) and compared the resulting \[N/O\] vs. \[O/H\] and time. In Figure 4 we show the two plots for the model galaxies which show that the observed \[N/O\] ratio is compatible with an age between 25 and 30 Myr, in good agreement with that estimated from the \[N/Fe\] ratio.

## 4 Discussion and conclusions

In the previous section we have shown that a SN-driven wind multi-zone model for the formation and evolution of an elliptical galaxy with initial luminous mass of $`210^{10}M_{}`$, which develops a wind in the most external region ($`3`$ kpc from the center) after $``$ 31 Myr, can well reproduce the dust-corrected abundances measured by Pettini et al. (2001) in the Lyman-break galaxy MS 1512-cB58. The fit to the chemical abundances basically suggests that MS 1512-cB58 could be a very young elliptical galaxy with an age $`20t_{sf}35`$ Myr, experiencing its main starburst and starting to develop an outflow. Independent age estimates (Ellingson et al., 1996), based on the SED of this galaxy, suggest also an age of $``$ 35 Myr, in good agreement with our independent estimate. The abundance ratios represent a strong constraint since they depend mainly upon the timescales of the production of the different chemical elements. On the other hand, the star formation rate depends strongly upon the galaxy mass and suggests a stellar mass of $`1.110^{10}M_{}`$ for MS 1512-cB58 and an age of $`100\pm 60`$ Myr, but the observed value is affected by the assumed Hubble constant whereas chemical abundances are not, and therefore this age estimate should be considered less certain than that derived from the abundances. One of the uncertainties involved in this comparison is that the abundance of Fe is very probably affected by dust depletion but we do not know exactly how much depletion occurs. We assumed that the observed Fe abundance should be increased by a factor of two and this assumption produces values for the \[$`\alpha `$/Fe\] ratios which are consistent with the \[N/Fe\] ratio. On the other hand, if we took the observed abundances without correction the higher \[$`\alpha `$/Fe\] ratios would suggest a younger age than what we have already inferred whereas from the higher \[N/Fe\] ratio a much older one, thus leading to a less consistent picture. Moreover, our predicted \[N/O\] vs. \[O/H\] is also in good agreement with the observed values and neither N or O should be substantially affected by dust depletion. In any case, dust depletion can only be part of the cause for the measured Fe underabundance, which is mostly due to the young age of the galaxy as compared to the typical timescale for Fe production by type Ia SNe in ellipticals ($`0.30.5`$ Gyr, as shown by Matteucci & Recchi, 2001).

In conclusion, this letter supports the idea that LBGs are young galaxies which could be the progenitors of present day spheroids (ellipticals and bulges have very similar properties), although we cannot predict the precise fate of cB58. We also like to stress the importance of using abundance ratios as cosmic clocks for dating high redshift objects.

## 5 Acknowledgments

We are grateful to Max Pettini for his many important suggestions and to an anonymous referee for constructive comments.