text:
# Intrinsic Decoherence Dynamics in Smooth Hamiltonian Systems: Quantum-classical Correspondence

## I Introduction

Quantum dynamics induces unitary transformations in Hilbert space, but most often it is only the dynamics projected onto a Hilbert subspace that is of interest. In general this reduced dynamics is nonunitary and therefore displays decoherence zurekreview . For example, if a system of interest is coupled to a bath, then averaging over the bath degrees of freedom introduces decoherence in the system dynamics. Likewise, in an isolated system, the reduced dynamics of a subsystem of this isolated system can display decoherence. We have termed decoherence in the latter case “intrinsic decoherence” since it does not involve an external bath batista .

Understanding decoherence is of crucial importance to a variety of modern fields such as quantum information processing qcbook  and quantum control of atomic and molecular processes brumerreview ; ricebook ; brumerbook . Our interest here is in the quantum-classical correspondence (QCC) between classical and quantum descriptions of the dynamics of decoherence. Specifically, we consider an initial quantum state subjected to either quantum or classical dynamics and compare the time evolution of the decoherence in both cases. We note that the formal theory of correspondence between quantum dynamics and classical Liouville dynamics wilkie  suggests that classical Liouville dynamics projected onto a subspace should also display decoherence. That is, as in the quantum case, the classical Liouville dynamics considered in the entire phase space is unitary and the classical Liouville dynamics projected onto a subspace is nonunitary. We therefore expect that the reduced classical Liouville dynamics propagated classically will show decoherence dynamics that is, at least qualitatively, parallel to that seen in the reduced quantum dynamics insofar as the loss of phase information, entropy production, etc. In the case of bath-induced decoherence we recently showed analytically that (a) one can indeed introduce a direct classical analog of quantum decoherence, and (b) examining the dynamics of decoherence classically gives new insights into both the dynamics of decoherence described quantum mechanically and into the conditions for QCC in decoherence dynamics gongprl .

Here we extend these considerations to intrinsic decoherence, both analytically and computationally. Specifically, in this paper we study QCC in the dynamics of intrinsic decoherence in smooth Hamiltonian systems, with an emphasis on the usefulness of classical dynamics in describing intrinsic decoherence. In particular, the classicality of early-time intrinsic decoherence dynamics is studied using a second-order perturbative treatment, and the interesting connection between decoherence rates at later times and the stability properties of classical trajectories is revealed by considering a simple approximate classical theory of intrinsic decoherence dynamics. The analytic and computational results shed new light on decoherence, dynamics of quantum entanglement, and quantum chaos. This study is also of interest to semiclassical decoherence studies bmiller , e.g., semiclassical descriptions of intrinsic decoherence dynamics in large molecular systems batista .

This paper is organized as follows. In Sec. II we introduce a second-order perturbation theory in an effort to understand QCC in early-time intrinsic decoherence dynamics. For simplicity we focus upon two degree-of-freedom systems, but the extension to larger systems is straightforward. Computational results of two sample cases in coupled-oscillator model systems, which strongly support the physical picture afforded by the perturbative treatment, are presented in Sec. III. Then, a classical theory of intrinsic decoherence dynamics for initially localized states is derived in Sec. IV. In the same section, detailed computational studies using this simple theory are carried out for the quartic oscillator model and one of its variants. Discussions and a summary comprise Sec. V.

## II Early-time Intrinsic Decoherence Dynamics

Consider a conservative system composed of two subsystems, with the total Hamiltonian given by

$`H(Q,P,q,p)={\displaystyle \frac{P^2}{2}}+{\displaystyle \frac{p^2}{2}}+V_1(Q)+V_2(q)+V_{12}(Q,q),`$ (1)
where $`(Q,P)`$ and ($`q,p`$) are dimensionless phase space conjugate variables, $`V_i`$ is the potential of the $`i`$-th subsystem, and $`V_{12}(Q,q)`$ describes arbitrary coupling between the two subsystems. As the system evolves the total system wavefunction $`|\psi (t)`$ becomes inseparable due to quantum entanglement, even if it is initially separable in $`Q`$ and $`q`$. As a result, measuring a subsystem would collapse the system wavefunction and therefore affect the properties of the other subsystem. Similarly, ignoring a subsystem decoheres the other one. The degree of intrinsic decoherence, which is induced by, and is a manifestation of, quantum entanglement between the two subsystems, can be measured by a well-known quantity: the quantum linear entropy linear-entropy

$`S_\mathrm{q}=1Tr_1(\widehat{\stackrel{~}{\rho }}^2),`$ (2)
where $`Tr_i`$ denotes a trace over the $`i`$-th subsystem, and $`\widehat{\stackrel{~}{\rho }}Tr_2\left(|\psi (t)\psi (t)|\right)`$ is the reduced density operator for the first subsystem. An increase in $`S_\mathrm{q}`$ suggests an increase of $`1/(1S_\mathrm{q})`$, which gives the number of orthogonal quantum states that are incoherently populated if the second subsystem is ignored. Below we choose $`q,p`$ as the “bath” variables and $`P,Q`$ as the system variables.

Let $`\rho _c(Q,P,q,p,t)`$ denote the phase space distribution function evolved classically, and $`\rho _W(Q,P,q,p,t)`$ denote the quantum (Wigner) phase space distribution function. Their time evolution equations are given by

$`{\displaystyle \frac{\rho _c}{t}}`$ $`=`$ $`\{H,\rho _c\},`$ (3)
$`{\displaystyle \frac{\rho _W}{t}}`$ $`=`$ $`\{H,\rho _W\}_M,`$ (4)
where $`\{\}`$ denotes classical Poisson bracket and $`\{\}_M`$ denotes quantum Moyal bracket moyal . We define classical and quantum reduced distribution functions as

$`\stackrel{~}{\rho }_c(Q,P,t)`$ $``$ $`{\displaystyle \rho _c(Q,P,q,p,t)𝑑q𝑑p},`$ (5)
$`\stackrel{~}{\rho }_W(Q,P,t)`$ $``$ $`{\displaystyle \rho _W(Q,P,q,p,t)𝑑q𝑑p}.`$ (6)
Since

$`S_\mathrm{q}(t)=12\pi \mathrm{}{\displaystyle \stackrel{~}{\rho }_W^2(Q,P,t)𝑑Q𝑑P},`$ (7)
where $`\mathrm{}`$ is the effective Planck constant, we can define a classical analog \[denoted $`S_c(t)`$\] to $`S_\mathrm{q}(t)`$ by replacing $`\stackrel{~}{\rho }_W`$ with $`\stackrel{~}{\rho }_c`$. That is,

$`S_c(t)12\pi \mathrm{}{\displaystyle \stackrel{~}{\rho }_c^2(Q,P,t)𝑑Q𝑑P}.`$ (8)
The main focus here is to compare $`S_c(t)`$ with $`S_\mathrm{q}(t)`$, i.e., the classical vs. quantum evolution of the intrinsic decoherence dynamics, as measured by the classical vs. quantum entropy.

Perturbative treatments have proved to very useful in understanding decoherence dynamics gongprl ; kim ; duan ; bacon . Here, to analytically examine classical vs. quantum intrinsic decoherence dynamics at early times, we apply the perturbative approach developed in our previous work gongprl  to the case of intrinsic decoherence dynamics. Specifically, consider a second-order perturbative expansion with respect to the time variable $`t`$ for both $`S_\mathrm{q}`$ and $`S_c`$, i.e.,

$`S_c(t)`$ $`=`$ $`S_c(0)+{\displaystyle \frac{t}{\tau _{c,1}}}+{\displaystyle \frac{t^2}{\tau _{c,2}^2}}+\mathrm{},`$
$`S_\mathrm{q}(t)`$ $`=`$ $`S_\mathrm{q}(0)+{\displaystyle \frac{t}{\tau _{\mathrm{q},1}}}+{\displaystyle \frac{t^2}{\tau _{\mathrm{q},2}^2}}+\mathrm{}.`$ (9)
Then, from the classical and quantum dynamics of the entire system one obtains

$`{\displaystyle \frac{1}{\tau _{c,1}}}=4\pi \mathrm{}{\displaystyle \stackrel{~}{\rho }_c(Q,P,0)\{H,\rho _c(Q,P,q,p,0)\}𝑑q𝑑p𝑑Q𝑑P},`$ (10)
$`{\displaystyle \frac{1}{\tau _{\mathrm{q},1}}}=4\pi \mathrm{}{\displaystyle \stackrel{~}{\rho }_W(Q,P,0)\{H,\rho _W(Q,P,q,p,0)\}_M𝑑q𝑑p𝑑Q𝑑P},`$ (11)
$`{\displaystyle \frac{1}{\tau _{c,2}^2}}`$ $`=`$ $`2\pi \mathrm{}{\displaystyle \stackrel{~}{\rho }_c(Q,P,0)\{H,\{H,\rho _c(Q,P,q,p,0)\}\}𝑑q𝑑p𝑑Q𝑑P}`$ (12)
$`2\pi \mathrm{}{\displaystyle \left[\{H,\rho _c(Q,P,q,p,0)\}𝑑q𝑑p\right]^2𝑑Q𝑑P},`$
and

$`{\displaystyle \frac{1}{\tau _{\mathrm{q},2}^2}}`$ $`=`$ $`2\pi \mathrm{}{\displaystyle \stackrel{~}{\rho }_W(Q,P,0)\{H,\{H,\rho _W(Q,P,q,p,0)\}_M\}_M𝑑q𝑑p𝑑Q𝑑P}`$ (13)
$`2\pi \mathrm{}{\displaystyle \left[\{H,\rho _W(Q,P,q,p,0)\}_M𝑑q𝑑p\right]^2𝑑Q𝑑P}.`$
Further, using the definitions of the classical Poisson and quantum Moyal brackets and assuming that initial classical and quantum distribution functions are identical and separable, i.e.,

$`\rho _c(Q,P,q,p,0)=\rho _W(Q,P,q,p,0)=\stackrel{~}{\rho }_1^0(Q,P)\stackrel{~}{\rho }_2^0(q,p),`$ (14)
we have

$$\frac{1}{\tau _{c,1}}=\frac{1}{\tau _{\mathrm{q},1}}=0,$$
(15)
$$\frac{1}{\tau _{c,2}^2}=2\pi \mathrm{}\left[\frac{\stackrel{~}{\rho }_1^0(Q,P)}{P}\right]^2C(0,0)𝑑Q𝑑P,$$
(16)
and

$`{\displaystyle \frac{1}{\tau _{\mathrm{q},2}^2}}`$ $`=`$ $`{\displaystyle \frac{1}{\tau _{c,2}^2}}+2\pi \mathrm{}{\displaystyle \underset{l_1l_20}{}\frac{[\mathrm{}/(2i)]^{(2l_1+2l_2)}}{(2l_1+1)!(2l_2+1)!}}`$ (17)
$`\times {\displaystyle \frac{^{(2l_1+1)}\stackrel{~}{\rho }_1^0(Q,P)}{P^{(2l_1+1)}}}{\displaystyle \frac{^{(2l_2+1)}\stackrel{~}{\rho }_1^0(Q,P)}{P^{(2l_2+1)}}}C(l_1,l_2)dQdP,`$
where $`C(l_1,l_2)`$ is a correlation function given by

$`C(l_1,l_2)`$ $``$ $`{\displaystyle \frac{^{(2l_1+1)}V(Q,q)}{Q^{(2l_1+1)}}}{\displaystyle \frac{^{(2l_2+1)}V(Q,q)}{Q^{(2l_2+1)}}}_{\stackrel{~}{\rho }_2^0}`$ (18)
$``$ $`{\displaystyle \frac{^{(2l_1+1)}V(Q,q)}{Q^{(2l_1+1)}}}_{\stackrel{~}{\rho }_2^0}{\displaystyle \frac{^{(2l_2+1)}V(Q,q)}{Q^{(2l_2+1)}}}_{\stackrel{~}{\rho }_2^0}.`$
Here $`_{\stackrel{~}{\rho }_2^0}`$ denotes the ensemble average over the zero-time “bath distribution function” $`\stackrel{~}{\rho }_2^0(q,p)`$. It is worth emphasizing that in our derivations we have used the same initial state for the classical and quantum dynamics.

Equation (15) shows that the zero first-order linear entropy increase rate, i.e., $`1/\tau _{\mathrm{q},1}`$=0, has a strict classical analog. Further, Eq. (16) indicates that classical Liouville dynamics also predicts a second-order entropy production rate $`1/\tau _{c,2}^2`$ that is the analog of the second-order quantum decoherence rate $`1/\tau _{\mathrm{q},2}^2`$. Thus, we can identify two categories of early-time intrinsic decoherence dynamics: classical if $`\tau _{c,2}\tau _{\mathrm{q},2}`$, and nonclassical if $`\tau _{\mathrm{q},2}`$ appreciably differs from $`\tau _{c,2}`$.

To simplify Eqs. (16) and (17) we introduce the Fourier transform \[denoted $`F(Q_1,Q_2)`$\] of $`\stackrel{~}{\rho }_1^0(Q,P)`$, i.e.,

$`F(Q_1,Q_2)`$ $``$ $`{\displaystyle \stackrel{~}{\rho }_1^0(\overline{Q},P)\mathrm{exp}\left[\frac{i\mathrm{\Delta }QP}{\mathrm{}}\right]𝑑P},`$ (19)
where $`\overline{Q}(Q_1+Q_2)/2`$ and $`\mathrm{\Delta }Q=Q_1Q_2`$. We then obtain

$$\frac{1}{\tau _{c,2}^2}=\frac{1}{\mathrm{}^2}|F(Q_1,Q_2)|^2\mathrm{\Delta }Q^2C(0,0)𝑑Q_1𝑑Q_2,$$
(20)
and

$`{\displaystyle \frac{1}{\tau _{\mathrm{q},2}^2}}={\displaystyle \frac{1}{\tau _{c,2}^2}}+{\displaystyle \frac{1}{\mathrm{}^2}}{\displaystyle |F(Q_1,Q_2)|^2\underset{l_1l_20}{}\frac{\mathrm{\Delta }Q^{(2l_1+2l_2+2)}}{(2l_1+1)!(2l_2+1)!}\frac{1}{2^{(2l_1+2l_2)}}C(l_1,l_2)dQ_1dQ_2}.`$ (21)
Equations (20) and (21) are general results. For the special case of $`V_{12}(Q,q)=f(Q)g(q)`$, Eqs. (20) and (21) can be rewritten in a simple and more enlightening form:

$`{\displaystyle \frac{1}{\tau _{c,2}^2}}={\displaystyle \frac{g^2(q)_{\stackrel{~}{\rho }_2^0}g(q)_{\stackrel{~}{\rho }_2^0}^2}{\mathrm{}^2}}{\displaystyle |F(Q_1,Q_2)|^2\mathrm{\Delta }Q^2\left[\frac{df(\overline{Q})}{dQ}\right]^2𝑑Q_1𝑑Q_2},`$ (22)
and

$`{\displaystyle \frac{1}{\tau _{\mathrm{q},2}^2}}={\displaystyle \frac{g^2(q)_{\stackrel{~}{\rho }_2^0}g(q)_{\stackrel{~}{\rho }_2^0}^2}{\mathrm{}^2}}{\displaystyle |F(Q_1,Q_2)|^2\mathrm{\Delta }Q^2\left[\frac{\mathrm{\Delta }f(\overline{Q})}{\mathrm{\Delta }Q}\right]^2𝑑Q_1𝑑Q_2},`$ (23)
where $`\mathrm{\Delta }f(\overline{Q})/\mathrm{\Delta }Q`$ is the finite-difference function

$`{\displaystyle \frac{\mathrm{\Delta }f(\overline{Q})}{\mathrm{\Delta }Q}}`$ $``$ $`{\displaystyle \frac{f(\overline{Q}+\mathrm{\Delta }Q/2)f(\overline{Q}\mathrm{\Delta }Q/2)}{\mathrm{\Delta }Q}}`$ (24)
$`=`$ $`{\displaystyle \frac{f(Q_1)f(Q_2)}{Q_1Q_2}}.`$
As a result: (1) If $`f(Q)`$ depends only linearly or quadratically upon the coupling coordinate $`Q`$, a common approximation, then $`(1/\tau _{\mathrm{q},2}^21/\tau _{c,2}^2)=0`$ for any initial state. That is, in this case there exists perfect QCC in early-time dynamics of intrinsic decoherence, regardless of $`\mathrm{}`$, and irrespective of the potentials $`V_1(Q)`$ and $`V_2(q)`$. (2) Even in the case of highly nonlinear $`f(Q)`$, as long as $`F(Q_1,Q_2)`$ decays fast enough with $`|Q_1Q_2|`$ such that $`\mathrm{\Delta }f/\mathrm{\Delta }Qdf/dQ`$, QCC would still be excellent. The smaller the $`\mathrm{}`$, the more rigorous is this requirement. (3) If $`\mathrm{\Delta }f/\mathrm{\Delta }Q`$ differs significantly from $`df/dQ`$ over the $`Q`$-coordinate scale of the initial state, quantum entropy production can be totally unrelated to classical entropy production. Such cases of poor QCC are of fundamental interest, but are not the focus of this paper.

The second-order perturbative treatment is most reliable for early-time dynamics and for relatively weak decoherence. The above results are particularly significant for studies on the control of intrinsic decoherence, where early-time dynamics of weak decoherence is important. In these circumstances it is useful to understand the extent to which (quantum) intrinsic decoherence is equivalent to classical entropy production, i.e. to increasing $`S_c(t)`$. In particular, if there exists good correspondence between classical and quantum decoherence dynamics, then the essence of decoherence control is equivalent to the suppression of classical entropy production, and various classical tools may be considered to achieve decoherence control. If not, then fully quantum tools are required.

The above perturbation results clearly demonstrate that quantum dynamics of intrinsic decoherence has a direct analog in classical Liouville dynamics. This rather intriguing result motivates us to computationally examine QCC in the dynamics of intrinsic decoherence over all time scales.

## III Computational Results: Two Sample Cases

To computationally examine QCC in the dynamics of intrinsic decoherence, we consider coupled-oscillator model systems with smooth Hamiltonians. In all the model systems studied below, we choose

$`V_1(Q)+V_2(q)={\displaystyle \frac{\beta }{4}}(Q^4+q^4),`$ (25)
where $`\beta =0.01`$. Since $`V_1(Q)`$ and $`V_2(q)`$ have no simple harmonic terms, any observed agreement between classical and quantum behavior cannot be attributed to the similarity between classical and quantum harmonic oscillator dynamics. If the coupling potential $`V_{12}(Q,q)`$ is quadratic in both $`Q`$ and $`q`$, i.e., $`V_{12}(Q,q)=\alpha Q^2q^2/2`$, then the resultant coupled-oscillator system is the well-known quartic oscillator model eckhardt ; gongpre ; channel . Because this model is well-studied and can display strongly chaotic (e.g., $`\alpha =1.0`$, $`\beta =0.01`$) or integrable (e.g., $`\alpha =0.03`$, $`\beta =0.01`$) dynamics, it is used in Sec. IV as an ideal model to study QCC in intrinsic decoherence dynamics for both integrable and chaotic cases.

Our perturbation theory approach predicts good classical-quantum agreement at short times for some potentials and initial conditions and poor agreement for others. We examine both these cases computationally.

It suffices to consider one case of poor agreement, since poor QCC at early times invariably translates to similar behavior at later times. Consider then $`V_{12}(Q,q)`$ to be some highly nonlinear potential. Computations of the quantum dynamics and thus the time dependence of $`S_\mathrm{q}(t)`$ are straightforward gongpre . $`S_c(t)`$ is computed directly using Monte-Carlo simulations with an importance sampling technique (where the Monte-Carlo simulations are based upon Eq. (27) below). From the analytical results above we see that $`V_{12}(Q,q)`$ and the scale of the initial state play decisive roles in QCC in early-time intrinsic decoherence dynamics. In particular, we expect poor QCC if $`V_{12}(Q,q)=f(Q)g(q)`$ differs significantly from a linear or quadratic function of $`Q`$ such that $`\mathrm{\Delta }f/\mathrm{\Delta }Q`$ differs significantly from $`(df/dQ)`$ over the $`Q`$-coordinate scale (i.e., the support) of the initial state. To confirm this computationally we consider $`f(Q)g(q)=\mathrm{sin}^2(10Q)q^2`$, with the initial distribution functions of the two subsystems given by

$`\stackrel{~}{\rho }_1^0(Q,P)`$ $`=`$ $`{\displaystyle \frac{1}{\pi \mathrm{}}}\mathrm{exp}\left[{\displaystyle \frac{(QQ_0)^2}{2\sigma _Q^2}}{\displaystyle \frac{(PP_0)^2}{2\sigma _P^2}}\right],`$
$`\stackrel{~}{\rho }_2^0(q,p)`$ $`=`$ $`{\displaystyle \frac{1}{\pi \mathrm{}}}\mathrm{exp}\left[{\displaystyle \frac{(qq_0)^2}{2\sigma _q^2}}{\displaystyle \frac{(pp_0)^2}{2\sigma _p^2}}\right].`$ (26)
Here the dimensionless effective Planck constant is chosen to be $`\mathrm{}=0.005`$ throughout, except for one case in Sec. V, and $`\sigma _Q/25=25\sigma _P=\sqrt{\mathrm{}/2}`$, $`\sigma _q=\sigma _p=\sqrt{\mathrm{}/2}`$, $`Q_0=0.5`$, $`P_0=0.5`$, $`q_0=0`$, with $`H(Q_0,P_0,q_0,p_0)=0.24`$. Note that $`\stackrel{~}{\rho }_1^0(Q,P)`$ is strongly squeezed in $`P`$ and that this initial distribution function is considerably delocalized in $`Q`$. Further, since $`|f(Q)|=|\mathrm{sin}^2(10Q)|1.0`$, $`|\mathrm{\Delta }Q||df/dQ|`$ can be much larger than $`|\mathrm{\Delta }f(Q)|`$. Thus, for this case the perturbation result predicts that at early times there can be substantial classical entropy production with insignificant quantum decoherence. As shown in Fig. 1, this is nicely confirmed by the numerical results of $`S_\mathrm{q}(t)`$ and $`S_c(t)`$. In particular, Fig. 1 shows that, at $`t=1.0`$ $`S_c(t)`$ (discrete points) is $``$ 0.9 while $`S_\mathrm{q}(t)`$ (solid line) is still less than 0.2. Evidently, QCC in this case is indeed very poor from the very beginning.

There remains then the important question of the quantitative degree of QCC in circumstances where our perturbative analysis predicts good short-time QCC. In particular, it is important to investigate whether or not good QCC predicted perturbatively remains for a considerable amount of time. If so, then the perturbative treatment provides a useful guide to our understanding of QCC in intrinsic decoherence dynamics. If not, then our perturbative results make sense only for extremely weak decoherence. Dramatically, our computational studies strongly support our analytical perturbation results, even in the presence of significant decoherence. For example, consider the case, where the parameters for the initial state are the same as in the previous case (therefore the initial state is also much delocalized), but the coupling potential is given by $`V_{12}(Q,q)=Q^2\mathrm{sin}^2(q)`$. This coupling potential is highly nonlinear in $`q`$ but still quadratic in $`Q`$. In accord with the second-order perturbation results, such a coupling potential should still give rise to good early-time QCC in the intrinsic decoherence dynamics of the first subsystem. This is confirmed by the quantitative comparison between $`S_\mathrm{q}(t)`$ and $`S_c(t)`$ shown in Fig. 2. More importantly, Fig. 2 shows that outstanding QCC remains even when both $`S_\mathrm{q}(t)`$ and $`S_c(t)`$ have increased to close to their saturation value of unity. Numerous other computational results (not shown) are consistent with the two cases shown here.

These results show the usefulness of the second-order perturbation theory in understanding QCC in intrinsic decoherence dynamics emanating from squeezed initial states. Also of interest is intrinsic decoherence dynamics associated with sufficiently localized initial states, which, in accord with the previous perturbation results, should display excellent early-time QCC for any coupling potential $`V_{12}(Q,q)`$. We now computationally examine QCC at later times for localized states as initial conditions and explain the results in terms of a simple classical theory of intrinsic decoherence dynamics.

## IV Localized Initial States

### IV.1 Simple Classical Approach

Below we show that $`S_\mathrm{q}(t)`$ and $`S_c(t)`$ are often in excellent agreement, over large time scales, for initially localized states, significantly extending the perturbation theory result. In doing so we compare the quantum $`S_\mathrm{q}(t)`$ with full classical mechanics as well as with a simple classical theory derived in this section. The latter provides further insight into the origins of increasing $`S_c(t)`$.

To derive the simplified classical result we first use Liouville’s theorem to reexpress $`S_c(t)`$ as

$`S_c(t)`$ $`=`$ $`12\pi \mathrm{}{\displaystyle \rho _c(Q(t),P(t),q(t),p(t),t)}`$ (27)
$`\times \rho _c(Q(t),P(t),q^{},p^{},t)dq^{}dp^{}dQ(t)dP(t)dq(t)dp(t)`$
$`=`$ $`12\pi \mathrm{}{\displaystyle \rho _c(Q,P,q,p,0)\rho _c(Q(t),P(t),q^{},p^{},t)𝑑q^{}𝑑p^{}𝑑Q𝑑P𝑑q𝑑p}`$
$`=`$ $`12\pi \mathrm{}{\displaystyle \rho _c(Q,P,q,p,0)\rho _c(Q^{\prime \prime },P^{\prime \prime },q^{\prime \prime },p^{\prime \prime },0)𝑑q^{}𝑑p^{}𝑑Q𝑑P𝑑q𝑑p},`$
where $`(Q(t),P(t),q(t),p(t))`$ is the phase space location of the trajectory emanating from $`(Q,P,q,p)`$ at $`t=0`$, and $`(Q^{\prime \prime },P^{\prime \prime },q^{\prime \prime },p^{\prime \prime })`$ is the phase space location of the trajectory at time zero if the classical trajectory is propagated backwards from $`(Q(t),P(t),q^{},p^{})`$. Because the initial state $`\rho _c(Q,P,q,p,0)`$ is assumed highly localized in phase space, $`(Q^{\prime \prime },P^{\prime \prime },q^{\prime \prime },p^{\prime \prime })`$ must be very close to $`(Q,P,q,p)`$ in order for the term $`\rho _c(Q,P,q,p,0)\rho _c(Q^{\prime \prime },P^{\prime \prime },q^{\prime \prime },p^{\prime \prime },0)`$ in Eq. (27) to be appreciable and thus to contribute to $`S_c(t)`$. Hence a convenient approximation can be made: we assume that, at time $`t`$, only those backward trajectories near $`(Q(t),P(t),q(t),p(t))`$ need be taken into account. This means that we treat $`Q^{\prime \prime }Q`$, $`P^{\prime \prime }P`$, $`q^{\prime \prime }q`$, $`p^{\prime \prime }p`$, $`\delta q^{}[q^{}q(t)]`$, and $`\delta p^{}[p^{}p(t)]`$ as sufficiently small such that

$`Q^{\prime \prime }`$ $``$ $`Q+M_{13}(t)\delta q^{}+M_{14}(t)\delta p^{},`$
$`P^{\prime \prime }`$ $``$ $`P+M_{23}(t)\delta q^{}+M_{24}(t)\delta p^{},`$
$`q^{\prime \prime }`$ $``$ $`q+M_{33}(t)\delta q^{}+M_{34}(t)\delta p^{},`$
$`p^{\prime \prime }`$ $``$ $`p+M_{43}(t)\delta q^{}+M_{44}(t)\delta p^{},`$ (28)
where $`M_{ij}`$ ($`i,j=1,2,3,4`$) is the stability matrix associated with the backward trajectories emanating from $`(Q(t),P(t),q(t),p(t))`$:

$`M_{ij}={\displaystyle \frac{(Q,P,q,p)}{(Q(t),P(t),q(t),p(t))}}.`$ (29)
Although this approximation should be less reliable for chaotic systems, we demonstrate below that, it is, nonetheless, computationally useful in both integrable and chaotic cases.

For simplicity we consider below a specific example in which $`\stackrel{~}{\rho }_1^0(Q,P)`$ and $`\stackrel{~}{\rho }_2^0(q,p)`$ are symmetric Gaussian states given by Eq. (26) with $`\sigma _Q=\sigma _P=`$$`\sigma _q=\sigma _p\sigma =\sqrt{\mathrm{}/2}`$. Other distributions can be considered in an analogous fashion. Substituting Eqs. (26) and (28) into Eq. (27) and evaluating the integrals, we obtain

$$S_c(t)=1\frac{1}{2}\frac{\mathrm{exp}\left[\frac{U^2X^2+V^2Y^22UVZ}{2\sigma ^2(X^2Y^2Z^2)}\right]}{\sqrt{X^2Y^2Z^2}}_{\rho _c^{}},$$
(30)
where

$`X`$ $`=`$ $`M_{13}^2+M_{23}^2+M_{33}^2+M_{43}^2,`$
$`Y`$ $`=`$ $`M_{14}^2+M_{24}^2+M_{34}^2+M_{44}^2,`$
$`Z`$ $`=`$ $`M_{13}M_{14}+M_{23}M_{24}+M_{33}M_{34}+M_{43}M_{44},`$
$`U`$ $`=`$ $`(QQ_0)M_{14}+(PP_0)M_{24}+(qq_0)M_{34}+(pp_0)M_{44},`$
$`V`$ $`=`$ $`(QQ_0)M_{13}+(PP_0)M_{23}+(qq_0)M_{33}+(pp_0)M_{43},`$ (31)
and

$`\rho _c^{}(Q,P,q,p)=4\pi ^2\mathrm{}^2[\stackrel{~}{\rho }_1^0(Q,P)]^2[\stackrel{~}{\rho }_2^0(q,p)]^2.`$ (32)
Equations (30) and (31) indicate that the classical dynamics of intrinsic decoherence is closely related to the classical stability matrix elements averaged over a rescaled initial distribution function. This interesting connection provides insight into a variety of interesting aspects of quantum intrinsic decoherence dynamics furuya ; miller ; arul ; fujisaki ; prosen . For example, for chaotic dynamics in which classical trajectories are highly unstable and therefore in which $`|M_{ij}|`$ increases rapidly, $`S_c(t)`$ should increase much faster than for the case of integrable dynamics. This observation can, with the assumption that there is fairly good QCC in intrinsic decoherence dynamics, directly explain previous results on quantum signatures of classical chaos in the dynamics of quantum entanglement furuya . Further, because Eqs. (30) and (31) are expressed in terms of classical stability matrices, characteristics of the time dependence of $`S_c(t)`$ and therefore of $`S_\mathrm{q}(t)`$ can be easily related to the time and space fluctuations in the instability of classical trajectories.

### IV.2 Computational Results: Localized Initial States

Consider then QCC over large time scales for localized initial states for both integrable and chaotic cases. To do so we examine the quartic oscillator model as well as results where the coupling potential is replaced by the nonlinear potential $`V_{12}(Q,q)=0.5Q^2q^2+Q^4q^2`$. The initial states are chosen to be localized initial states, and both the full classical dynamics and the approximate time dependence of $`S_c(t)`$ in Eq. (30) are compared to the quantum result. Specifically, we realize the ensemble average in Eq. (30) by Monte-Carlo simulations, using only $`2\times 10^4`$ sampling classical trajectories from which the stability matrix elements $`M_{ij}`$ are evaluated. The initial Gaussian states are chosen symmetric with $`\sigma =\sqrt{\mathrm{}/2}=0.05`$ and are sufficiently localized so that Eq. (28) should be a valid approximation. Full classical results for $`S_c(t)`$ (represented again by discrete circular points), which are much more demanding computationally, are also provided below.

Figures 3 and 4 compare results for $`S_c(t)`$ obtained from Eq. (30) and from exact classical results with $`S_\mathrm{q}(t)`$ (dashed line) for integrable and chaotic dynamics in the quartic oscillator model, respectively. A number of observations are in order. First, it is clear that in both cases the approximate $`S_c(t)`$ are in excellent agreement with the full classical results, confirming the utility of the simple model \[Eq. (30)\]. Second, both $`S_c(t)`$ and $`S_\mathrm{q}(t)`$ are seen, in the chaotic case, to relax faster towards $`1.0`$ than they do in the integrable case. Third, the oscillation amplitudes of $`S_\mathrm{q}(t)`$ in the chaotic case are much smaller than that in the integrable case. Hence, the fast relaxation and small-amplitude oscillations of $`S_\mathrm{q}(t)`$ shown in Fig. 4 may be regarded as fingerprints of the underlying classical chaos. Finally, and most importantly, the entire time dependence, including oscillations in Fig. 3 and Fig. 4 of $`S_\mathrm{q}(t)`$ are beautifully captured by both the exact and the approximate $`S_c(t)`$. It should also be noted that the QCC time scale shown in Fig. 4 is appreciably longer than is the QCC break time $`t_b5.0`$ for the same $`\mathrm{}`$, obtained by quantitatively comparing the structure of the classical and quantum distribution functions gongthesis . This can be understood by the fact that $`S_c(t)`$ \[or $`S_\mathrm{q}(t)`$\] describes the reduced distribution functions $`\stackrel{~}{\rho }_c(Q,P,t)`$ \[or $`\stackrel{~}{\rho }_W(Q,P,t)`$\] which is insensitive to the fine structure of $`\rho _c(Q,P,q,p,t)`$ \[or $`\rho _W(Q,P,q,p,t)`$\].

Calculations for many other initial states confirm that the QCC results shown in Figs. 3 and 4 are typical, indicating that (a) QCC is essentially exact over large time scales and (b) the simple classical theory of intrinsic decoherence dynamics introduced above provides a useful approximation to the exact results.

Figure 5 shows one case, however, where the approximate $`S_c(t)`$ and exact classical or quantum results differ quantitatively. Here the system is still the quartic oscillator model with $`\alpha =1.0`$, $`\beta =0.01`$, but with an initial state of special type. In particular, both the initial average position $`Q_0`$ and the initial average momentum $`P_0`$ are set to zero. Initial states of this type are called channel states channel , and effectively give rise to very weak coupling between the two subsystems over a considerably large time scale. Indeed, since at short times $`df/d\overline{Q}\mathrm{\Delta }f(\overline{Q})/\mathrm{\Delta }Q0`$ for $`\overline{Q}Q_0=0`$, one obtains from Eqs. (22) and (23) that $`1/\tau _{c,2}^2=1/\tau _{2,\mathrm{q}}^20`$. Therefore the early-time intrinsic decoherence rate should be small, as seen in Fig. 5, although the underlying classical dynamics is strongly chaotic. For this reason one expects that the dynamical behavior of $`S_c(t)`$ and $`S_\mathrm{q}(t)`$ should differ from previous cases. As shown in Fig. 5, in this case both $`S_\mathrm{q}(t)`$ and $`S_c(t)`$ increase in a step-wise fashion, distinctly different from that in Figs. 3 and 4. Agreement between them remains excellent. However, the approximate $`S_c(t)`$ misses some of the important structure.

In Fig. 6 we show the QCC result for a simple variant of the quartic oscillator model, i.e., $`V_{12}(Q,q)=0.5Q^2q^2+Q^4q^2`$, a coupling potential that is neither linear nor quadratic. As seen in Fig. 6, even with such nonlinear coupling $`S_c(t)`$ and $`S_\mathrm{q}(t)`$ are in excellent agreement over large time scales. This emphasizes the fact that the good QCC results observed in the quartic oscillator model are not due to the fact that the coupling potential therein is quadratic. Hence, we conclude that for initially localized states, our simple classical theory of intrinsic decoherence dynamics \[see Eq. (30)\] is generally useful in describing intrinsic decoherence dynamics in smooth Hamiltonian systems.

## V discussion and summary

Quantum entanglement between individual subsystems has no classical analog. Nevertheless, as shown in this work, the quantum dynamics of quantum entanglement, as manifest in the quantum dynamics of intrinsic decoherence, does have a classical analog in classical Liouville dynamics describing classical correlations between classical subensembles. Hence it is useful to isolate the conditions under which there is good QCC in intrinsic decoherence dynamics. This is done analytically, for early-time dynamics and for weak decoherence, by a second-order perturbative theory. Interestingly, as demonstrated by our computational studies in Sec. III and Sec. IV, the physical picture of QCC afforded by the perturbative treatment can be still very useful even when the time scale under investigation is relatively long and the degree of intrinsic decoherence is significant. In particular, under the circumstances where there is good early-time QCC, classical Liouville dynamics can provide a simple means of understanding different aspects of intrinsic decoherence dynamics, for relatively large time scales and for both integrable and chaotic dynamics. Further, we have derived an approximate but very simple classical theory of linear entropy production of intrinsic decoherence dynamics associated with localized initial states, and shown that the rate of entropy production is closely related to the stability properties of classical trajectories.

Clearly, the linear entropy is just one of many possible representation-independent measures of intrinsic decoherence, and $`S_\mathrm{q}(t)S_c(t)`$ does not mean that the quantum dynamics is equivalent to the corresponding classical Liouville dynamics. For example, if the saturation value of the linear entropy in the long time limit is of particular interest, then the measures $`1/[1S_c(t)]`$ and $`1/[1S_\mathrm{q}(t)]`$ (which gives the number of orthogonal states that are incoherently populated) should be more useful in describing QCC. Indeed, our results in Figs. 4 and 6 suggest that as time increases one has $`1/[1S_c(t)]>>1/[1S_\mathrm{q}(t)]`$. This is consistent with our previous observation gongprl  that, decoherence can dramatically improve QCC, but even strong decoherence does not necessarily suffice to ensure that quantum entropy production is the same as classical entropy production.

It should also be pointed out that the model quantum systems studied in this paper are still far from the semiclassical regime. This is indicated, in the chaotic case of the quartic oscillator model for example, by the fact that the QCC break time is relatively short compared to the time scale that we examined. Correspondence will worsen quantitatively with increasing $`\mathrm{}`$, although, as discussed above, $`\mathrm{}`$ is far from the only factor influencing the quality of the QCC. However, we note that for localized initial states the qualitative features of the time dependence of classical and quantum linear entropy may remain similar to one another with much larger effective Planck constants. For example, Fig. 7 displays fairly good QCC between $`S_\mathrm{q}(t)`$ and $`S_c(t)`$, in the chaotic case of the quartic oscillator model, with $`\mathrm{}=0.05`$ and with an initial symmetric Gaussian state.

A number of interesting extensions of this work are under consideration. First, it seems straightforward but necessary to consider cases in which the coupling potential depends upon both position and momentum. Second, we propose to further investigate the role of the dynamics of the subsystems in addition to that of the coupling potential (e.g., the dynamics in coupled Morse oscillator systems). Third, it is interesting to study QCC in intrinsic decoherence dynamics in terms of the decay of off-diagonal density matrix elements. Such studies are ongoing, with preliminary studies han  indicating that comparing the time dependence of off-diagonal density matrix elements to its direct classical analog gongprl  will provide deeper insights into QCC in intrinsic decoherence dynamics.

To summarize, we have shown that classical dynamics can be very useful in describing intrinsic decoherence dynamics in smooth Hamiltonian systems. In particular, we have identified conditions under which excellent quantum-classical correspondence in the early-time dynamics of intrinsic decoherence is possible via a second-order perturbative treatment, have presented a simple classical theory of intrinsic decoherence dynamics emanating from localized initial states, and have provided supporting computational results. The hope is that by extending this study to high-dimensional Hamiltonian systems, we may use purely classical approaches to describe (at least qualitatively) the dynamics of quantum entanglement or intrinsic decoherence in polyatomic molecular systems.

Acknowledgments: This work was supported by the U.S. Office of Naval Research and by the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# Java Physics Generator and Analysis Modules

## I INTRODUCTION

In performing Linear Collider physics and detector studies, one encounters a number of existing Fortran77, and new C++ and Fortran95 software modules for event generation, detector simulation and analyis. Stand-alone packages exist for Pandora(C++), Pythia(F77) and Whizard(F95) physics generation, and in some cases have been integrated into separate simulation packages for TESLA SimDet(F77), Asian JLC QuickSim(C++) and American Java Fast Monte Carlo detector studies. Individual physics analyses proceed by adopting one of the generator and detector simulation packages and often are hard to compare to similar analyses done in different environments. A Java software framework lcstudies  allows the use of different event generators, detector simulation packages and analyis modules in the same high level Object-Oriented (OO) environment. An environment which provides modern development tools for Java physics and detector analysis modules using prototype HEP class libraries hep.lcd . The model allows networked graphical interactive applications, or batch single-pass processing of complete physics generation, parallel detector simulations and event analysis for full 500 fb-1 simulated data samples. Java histogram objects files can be saved for final presentation using the Java Analysis Studio (JAS) jas .

In this paper, Java Native Interfaces (JNI) for Pandora, Pythia and Whizard Linear Collider physics event generators are described and used in defining a Java framework with a common generator interface package, Sec. II. Sections III and IV list the generator and detector simulation modules, respectively. A simple comparison of the different physics generators is made for a LC bench mark process in Sec. V. Detector simulations using the U.S. Fast Monte Carlo, TESLA SimDet and Asian QuickSim packages are compared in Sec. VI. Details of interface implementations for the generators, and for executing SimDet and QuickSim detector simulations and accessing simulated quantities are outlined in Sec. VII. Java software package organization and a brief discussion of the documentation is left for Sec. VIII.

## II JAVA ANALYSIS FRAMEWORK

Several prototype Java interfaces, classes and libraries have been writtenlcstudies  to provide a high level object-oriented physics generation, detector simulation and analysis environment for $`e^+e^{}`$ Linear Collider studies.

### II.1 Generator interface

A protype Java object model has been developed to generalize the interface to physics generators such as Pandora, Pythia and Whizard. The following methods are implemented through Java Native Interfaces (JNI) to the underlying C++, Fortran77 or Fortran95 package:

* setup() - Set up generator configuration using default parameters.
* initialize() - Initialize for selected beams, energy and physics process with given parameters.
* generateEvent() - Generate an event.
* terminate() - Summarize generated data sample.

One additional method is required to provide information about the generator.

* getName() - Return generator name and version number.

Generator specific methods have been added such as

* list(int) - List particles in the event.

The useage of the Java generator interface is outlined in the following excerpt:

```
// Import the Pythia generator classes.
import hep.generator.pythia.*;

// Create a Pythia process.
  pythia = new Pythia();

// Set up a "User-defined" Pythia process
// with a String of parameters settings.
  pythia.give(parameters);

// Initialize the e+e- Center-of-Mass system.
  pythia.init("CMS","e+","e-",Ecm);

// Run Pythia to generate and list events.
  for (int n=1; n<=NEvents; n++) {
    pythia.generateEvent();
    pythia.list(1);
  }

```

Generator processes interfaces have been defined for classes which dynamically load predefined process settings, such as PythiaProcess for initializing Pythia. A number of reference process classes have been written, such as eetoZH for Higgsstrahlung and eetottbar for Top-pair production processes.

### II.2 EventGenerator and analysis modules

The generator interface described above is used to provide “EventGenerator” modules for the American LCD frameworkhep.lcd . A Java package hep.lcd.generator is abstracted to introduce different generator modules as alternative “EventSource’s”:

* Generator - Defines a standard HEP physics event generator module.
* GeneratedEvent - Defines a generated event.
* HEPEvent - Provides a standard /HEPEvt/ implementation of a generated event.

The useage of a LCD EventGenerator module is illustrated by

```
// Import the LCD generator classes for Pythia.
import hep.lcd.generator.pythia.*;

// Create a Pythia generator for given process.
  generator = new PythiaGenerator(processName);

// Run to generate and dump the first event.
  for (int n=1; n<=NEvents; n++) {
    event = generator.generateEvent();
    if (n==1) event.dump();
  }

```

Within the LCD Java framework, event processing and analysis modules are added as “Driver’s” or “Processor’s” defined in a hep.lcd.util.driver package. The following classes outline the use of detector simulation and analysis modules:

```
public class AnlPythiaZHJets implements Driver

  public AnlPythiaZHJets()

    add(detectorSimulation);
    add(analysis);
 or
    add(new MCFast());
    add(new SimDetModule());
    add(new QuickSimModule());
    ...
    add(new PythiaZHJetAnalysis());

class PythiaZHJetAnalysis implements Processor
  void process(LCDEvent event)
    particles = event.get("MCParticles");
    tracks  = event.get("SimDet Tracks");
    objects = event.get("Energy Flow Objects");
    ...

```

Note that the Java Virtual Machine (JVM) allows multiple detector simulations to run in independent computer memory space.

### II.3 Analysis classes and tools

The utility of Java physics and analysis class libraries is illustrated in the following excerpts for the Particle physics class

```
  int PDGID = particle.getType().getPDGID();
  double[] PV = particle.getMomentum();
    double pT = Math.sqrt(PV[0]*PV[0]+PV[1]*PV[1]);
    double tanTh = pT/PV[2];

```

and the Histogram analysis class

```
  histogram("Mass").fill(particle.getMass());

```

Java framework applications can be built using standard HEP analysis tools and run through convenient Java packages so that

* gmake - Compiles Java and native code then creates a shared object library.
* java -mx64M GenWhizard … - Loads Java and all shared object libraries, executes Whizard and writes a xxx.javahist unbinned histogram file.

After completing several jobs, the Java Analysis Studio jas  can be used to combine histograms and make presentation plots.

```
  jas &
    Open files: xxx.javahist yyy.javahist ...
    Overlay and format plots
    Save as PostScript files

```

## III MONTE CARLO EVENT GENERATORS

Three representative LC physics generators have been used in this development: Pandora V2.2, Pythia V6.2 and Whizard V1.22. Pythia which has been extensively tested by many experiments was used to provide checks on specific processes. Pandora and Whizard are new developments which also offer complete event generation capability.

### III.1 Pandora Monte Carlo

The Pandora V2.2 Pandora  Monte Carlo provides simulations of selected two- and three-body processes and some illustrative beyond the Standard Model processes. Initial state radiation, beam polarization, and final state spin correlations are included for all processes. In stand-alone applications Pandora is used with the Pandora-Pythia V3.2 Pandora-Pythia  interface package to generate final state events. A “straight-forward” interface from the Java language to underlying Pandora C++ methods, Sec. VII, is provided in a hep.generator.pythia Java package, Sec. VIII.

### III.2 Pythia Monte Carlo

The Pythia V6.2 Pythia  Monte Carlo program allows generation of high-energy physics events including hard and soft interactions, parton distributions, initial and final state parton showers, multiple interactions, fragmentation and decay. A Java-C-Fortran77 interface described below provides full Pythia functionality within the Java environment.

### III.3 Whizard Monte Carlo

The Whizard V1.22 Whizard  Monte Carlo system is designed for the efficient calculation of multi-particle scattering cross sections and simulated event samples. Tree-level matrix elements are generated automatically for arbitrary partonic processes by calling external programs (O’MEGA, MADGRAPH and CompHEP). In this development, Java-C-Fortran95 interfaces to new Whizard procedures for initialization, phase space integration, event generation and summary was written to allow easy portability and compiler independence.

## IV DETECTOR SIMULATION MODELS

Models of the expected LC detector response to 500 GeV $`e^+e^{}`$ annihilation events are being developed. Preliminary versions used in this study are described below, and compared in Sec. VI. The detailed Java interface implementations are outlined in Sec. VII.

### IV.1 LCD Fast Detector Simulation

The American LCD V1.4 detector simulation includes charged particle momentum smearing based on detailed error estimates, gaussian energy smearing for photons and neutral hadrons, and acceptance and energy threshold requirements. It doesn’t account for confusion in reconstructing calorimeter energy deposits resulting in perfect energy flow separation of charged and neutral components of jets.

### IV.2 TESLA SimDet Detector Simulation

The European SimDet V4.0 detector simulation includes parameterized charged and neutral energy smearing based on full (Brahms) Monte Carlo simulations, acceptance requirements, and a realistic energy flow algorithm.

### IV.3 JLC QuickSim Detector Simulation

The Asian QuickSim detector simulation includes the generation of tracker hits followed by track and vertex reconstruction. Simulated calorimeter cell energy depositions are formed based on detailed detector simulations. Detailed cluster finding and track-cluster matching are performed before writing out combined track information..

## V PHYSICS GENERATOR COMPARISONS

The Pandora, Pythia and Whizard Linear Collider physics event generators are compared in Fig 1 for LCD Fast Monte Carlo detector simulations of Higgsstrahlung, $`e^+e^{}`$$``$$`ZHiggs`$, events. In these studies all generators were run at ECM = 500 GeV, and the Higgs mass was set to 115 GeV. Pythia hadronization was used following all parton-level event generation processes. For Pandora, the NLC500 machine parameters were selected in simulating the beamstrahlung energy spectrum. While both Pythia and Whizard event generators used Circe beamstrahlung simulation of the TESLA machine. These different machine settings are not significant in this comparison.

Overall, the generators were found to be in good agreement. In this study, Whizard V1.22 generated events did not include the natural Z width as seen in the jet-jet mass distributions shown in Fig 1d; a Whizard deficiency that has been corrected following this work.

In related studies higgs , the Java framework was used to compare the simulation of background processes for $`WW`$$`ZZ`$ and $`q\overline{q}`$ production from Pandora with explicit final state simulations in Whizard. These comparisons showed discrepancies due to the more precise treatment of polarization in Pandora for exclusive processes. The simulation of inclusive, interferring processes in Whizard was found to be the most complete, and the Whizard MadGraph amplitude calculation provided handelling of parton color flow information.

## VI DETECTOR COMPARISONS

In a $`e^+e^{}`$ Linear Collider the Higgs can be reconstructed through hadronic decayshiggs  even in multijet final states. The jet-jet mass distributions for the American LCD Fast Monte Carlo (FMC), TESLA SimDet and JLC QuickSim detector simulations of Higgstrahlung signal-only events with the Z and Higgs decaying hadronically are shown in Fig 2. In these current simulations, the LCD FMC jet-jet mass resolution is significantly better since it assumes “perfect” energy flow. More-realistic TESLA SimDet and JLC QuickSim detector simulations give comparable jet energy resolutions but different mean reconstructed jet-jet masses. Both these simulation models may be improved with the development of new techniques to resolve charged and neutral calorimeter clusters, and to correct for the loss of low energy tracks and clusters.

The overall performance of the Java framework, event generators, detector simulations and jet analyis was quite acceptable for full 500 fb-1 production runs. Each generator was run separately with all three detector simulations being performed for the same generated event. The Java framework itself and event generation and fragmentation typically took much less than 100 msec per event, with the longest times for multi-body processes. The detailed native SimDet and QuickSim detector simulations each took about 300-500 msec per event, depending on the type of events. The Java jet-finding which included multiple passes to merge jets and the subsequent analysis took roughly 100 msec per event for each type of particle jet being found, e.g. for Monte Carlo particles, for jets found from reconstructable tracks and clusters, or for Energy Flow object jets. The typical overall event generation, multiple detector simulation, and jet analysis took about 1-1.5 seconds per event.

Batch jobs were labelled by the generator name and parameters that were set, and histogram folders were created for the different detector simulations. In the analysis of results, JAS would open separate output histogram files in different folders and use the folder names to automatically generate plot legends.

## VII JAVA NATIVE INTERFACE (JNI) IMPLEMENTATIONS

The physics generator, physics event, detector simulation and software framework implementations111http://www-lc.lbl.gov/software/docs/ for all developed software modules are outlined in Table 1.

In generating events through the Pandora-JNI, a Java-C++ interface is used to allow the setting of different beam configurations, and to create physics processes to be added to Pandora’s Box. Pandora processes parameters, such as the Higgs mass, are required for the Java interface construction. Once a Pandora process has been created no Java object to C++ object interface is needed. Pandora generates primary partons which are accessed through the HEPEvt-JNI and passed to Pythia for fragmentation.

The Pythia-JNI allows for independent event generation and fragmentation processing. The interface can be used to pass configurations through pygive or to use Java classes such as “eetoZH” to configure Pythia for simulating the Higgsstrahlung process. The Pythia interface will be quickly replaced when new C++ versions of Pythia are made available and tested.

The Whizard interfacing to Fortran95 was simplified by creating a few simple Fortran95 modules. Whizard is configured by reading a “whizard.in” data file. Future developments of the Whizard Java interface are planned.

Java interfaces to the HEPEvt data structure and StdHEP output were written for multi-purpose use. After a generator fills HEPEvt, framework modules read the information and create event objects. These objects are later downloaded into HEPEvt structures that are accessed by fragmentation processes or by the different detector simulations. The design allows any combination of physics generator, fragmentation process or detector simulation. User analysis interfaces have also been written to allow one to use existing analysis packages.

An interface to Circe Circe  is used to handle the initialization of beamstrahlung simulation parameters from the Java framework for the Pythia and Whizard generators. The Tauola Tauola  $`\tau `$-decay package was used to handle polarization effects. A simple Java interface was written to allow testing of individual Tau decay modes.

The detector simulation packages have simple interfaces to set up, initialize and execute the simulations. Care was taken to pass events into the packages through separate HEPEvt interfaces. The simulated detector responses for SimDet and QuickSim are accessed from underlying storage into Java track and cluster classes to form reconstructed objects that are added to the LCD event for subsequent analysis. The SimDet and QuickSim interfaces follow a similar design model which could be generalized for other applications.

## VIII JAVA SOFTWARE PACKAGE AND API ORGANIZATION

The interface software and LCD modules have been organized into several Java packages or class libraries. The interfaces to the imported physics generator and detector simulation modules are in the HEP packages listed in Table 2. The LCD framework modules are organized in the packages listed in Table 3. The LCD modules provide reference implementations using the more generic interfaces to the HEP packages.

The LCD event generator and detector simulation modules are designed for simple use, see Sec. II.2. The Generator, Fragmentation and GeneratedEvent interfaces allow the design of abstract models for sourcing events to be based to the different detector simulation modules. Simulated reconstructed objects implement a basic Particle class allowing the development of abstract analysis methods to reconstruct event quantities from different detector models.

Special utility packages, pandorapythia and whizardpythia, based on the Pandora-Pythia Pandora-Pythia  software model for handelling the physics interfacing from Pandora and Whizard generated parton color flow information to Pythia hadronization, are located in the hep.lcd.util package.

Software Application Programming Interface (API) documentation on all Java packages is available on the web222http://www-lc.lbl.gov/software/java/api/. Class inheritance and method signiture documentation is automatically generated. Some documentation exists for most classes while comments on many self-explanatory methods are missing.

In building stand-alone event generation, detector simulation and physics analysis packages, the Java classes are compiled according to the following rule

```
# For a collection of Java source files
 JAVAFILES := $(wildcard hep/.../*.java)

# the corresponding class files are obtained
 CLASSFILES:= $(patsubst %.java, \
                %.class,$(JAVAFILES))

# and compiled with
%.class: %.java
   @echo "Compiling" $(*F)
   @javac $(JAVACFLAGS) $<

```

Header files for the JNI implementations are automatically generated from the Java source code using a utility

```
$(HeaderFile): $(JNISource)
   @if [ -f $(TargetFile) ]; \
      then rm $(TargetFile); fi;
   @javah -jni -d $(TargetDir) $(Source)

```

Shared object libraries and native code are compiled following

```
   gcc -shared -I$(JavaJDK)/include \
     -I$(JavaJDK)/include/linux \
     ...
     -L/usr/lib/gcc-lib/i386-redhat-linux/2.96 -lg2c \
     -L/usr/lib -lc \
     -o $(SharedLib)

```

where JavaJDK locates the Java development kit (JDK), the include files are used by the JNI implementations, and the library list usually includes Linux specific libraries and general libraries that normally would be loaded into an executable image. The shared object library is located by the LD\_LIBRARY\_PATH environment variable and loaded into the JavaVM following

```
public class Pythia
{
  ...

  /** Static initializer loads Pythia library. */
  static {
    String libname = "pythia";
    System.out.println(" Load lib"+libname+".so");
    System.loadLibrary(libname);
  }
}

```

###### Acknowledgements.

I thank Tony Johnson for many useful discussions on Java object modeling and for presenting this paper at the conference. Tony provided an initial interface to the HEPEvt data structure. Michael Peskin and Tim Barklow provided guidance on interfacing to Pandora and on passing the parton color flow information to Pythia for fragmenting high-energy $`e^+e^{}`$ jets. The Pythia generator interface was extended to include pygive & pyinit in consultation w/ Stephen Mrenna. Tim Barklow, Wolfgang Kilian and Thorsten Ohl helped in defining the Whizard interface and in implementing simple Fortran95 procedures. I thank Heinz-Jürgen Schreiber and Akiya Miyamoto for reviewing the detailed interfacing to the TESLA SimDet and JLC QuickSim detector simulation packages. The Java physics generator and detector simulation modules described in this paper are used within the Java based software framework under development by the American LCD software grouplcstudies ; hep.lcd ; jas . Work supported by the U.S. Department of Energy contract DE-AC03-76SF00098.


text:
# Unexpected residual regularities in liquid drop mass calculations

## I Introduction

The prediction of nuclear masses is of fundamental importance for a complete understanding of the nuclear processes that power the Sun, and the synthesis and relative abundances of the elements Rol88 . Möller at al. Moll95 , Duflo and Zuker Duf94 , Goriely et al. Gor01 , among many others, have developed mass formulae that calculate and predict the masses (and often other properties) of as many as 8979 nucleotides. There is a permanent search for better parameterizations that decrease the difference with the experimental masses and produce reliable predictions for unstable nuclei.

The atomic mass excesses have been tabulated for 8979 nuclei ranging from 16O to A=339 in Moll95 , calculated with a Finite Range Droplet macroscopic Model and a folded-Yukawa single-particle microscopic model, FRDM for short. The microscopic sector includes a Lipkin-Nogami calculation of pairing gaps besides Strutinsky shell-corrections. Ground-state energies are minimized with respect to shape degrees of freedom. Only 9 constants are adjusted to the ground-state masses of 1654 nuclei, with a mass model error of 0.669 MeV in the entire region of nuclei ranging from 16O to 263106.

The average mass error is found to gradually increase in a systematic way as the lighter region is approached. This behavior was related in Moll95  with the reliability of the Strutinsky method for the folded-Yukawa single-particle method, which is expected to be less accurate for light nuclei because the smooth, average quantities are less accurately determined from the relatively few levels occurring in light nuclei.

The main motivation for the present work arose from the detailed study of the nuclear mass error distribution shown in Fig. 1, which displays the number of nuclei within a given interval of FRDM mass errors Moll95 , for the 1654 isotopes with measured masses in 1995. The symmetrical distribution with two peaks around the origin was completely unexpected. In this article we study the regularities in the differences between measured and calculated masses, showing that they are closely related to the presence of this double peak.

Fig. 1 shows the FRDM distribution of nuclei Moll95  for different intervals of mass errors. The intermediate curve (diamonds) corresponds to a mass error interval of $`0.1`$ MeV where the double peak is quite evident.

Enlarging the mass error interval to $`0.2`$ MeV (triangles) softens the curvature but the presence of the double peak is still very clear. In the opposite direction, an interval of $`0.05`$ MeV (squares) produces larger fluctuations and some apparent interference, but it is evident that the distribution still displays the peaks. We can thus conclude that the effect is real and not an artifact of a particular distribution interval.

We carried out additional tests to verify the persistence of the double peak, including only nuclei with A $`>`$ 40 and A $`>`$ 60, normalizing the distribution according to the rule $`A^{1/3}`$, and including the experimental errors. In all circumstances the double peak remained Vel03 .

Figure 2 contains the distribution of FRDM mass differences Moll95

$$\mathrm{\Delta }M(N,Z)M_{exp}(N,Z)M_{th}(N,Z),$$
(1)
in the proton number (N) - neutron number (Z) space. We can see large domains with a similar error (each tone is associated to the magnitude of the error). It is remarkable that very well defined correlated areas of the same gray tone exist, which is a clear indication of remaining systematics and correlation.

Both the distribution of mass errors in the N-Z plane and the presence of the double peak in the error distribution suggest the presence of important correlations which are not taken into account in the FRDM mass formula. In the next sections these correlations are quantified and parametrized.

## II Regularities in the mass errors

In order to visualize the oscillatory patterns suggested in Fig. 2, different cuts were performed along selected directions on the N-Z plane. Given the large number of chains which can be studied, we have selected those with the largest number of nuclei with measured mass. For each cut a Fourier analysis was performed, and the squared amplitudes are plotted as a function of the frequencies on the right hand side of each figure. The discrete Fourier transforms $`F_k`$ are calculated as

$$F_k=\frac{1}{\sqrt{N}}\underset{j}{}\frac{\mathrm{\Delta }M(j)}{\gamma }\text{exp}\left(\frac{2\pi ijk}{N}\right),$$
(2)
where N is the number of mass differences $`\mathrm{\Delta }M`$ in a given series. The parameter $`\gamma `$ makes $`F_k`$ dimensionless. Given that it only affects the global scale of the Fourier amplitudes, we made the simple selection $`\gamma =1`$ MeV. The Fourier amplitudes are plotted as functions of the frequency $`f=k/N`$. The presence of a few frequencies with notorious large components underlines the existence of an oscillatory behavior of the mass errors.

### II.1 Fixed N or Z

We start our analysis for fixed N or Z, i.e. we selected different chains of isotopes and isotones. Those isotopic chains with 20 or more nuclei with measured masses are presented in Fig. 3 and 4. Until 1995, the element which had most isotopes with measured masses was Cs (Z=55) with 34.

Fig. 3 displays the mass errors for the isotope chains Z = 46 to 56, and their Fourier analysis, nearly all exhibiting a clear peak for the low frequency $`f1/20=0.05`$.

Fig. 4 displays the mass errors for the isotope chains Z = 30,36, 37, 38, 40, 87, 89, 91, and their Fourier analysis.

A similar analysis was performed for fixed N, i.e. isotonic chains, which have up to 30 nuclei with measured masses. Fig. 5 displays the mass errors for the isotone chains N=50 to 61, and their Fourier analysis.

Fig. 6 displays the mass errors for the isotone chains N=78 to 89, and their Fourier analysis.

Many isotopic chains exhibit similar oscillatory patterns in the mass errors. However, a Fourier decomposition for these cases is crude because there are only a few points, typically between 15 and 30, in each chain.

### II.2 Fixed A or Tz

To map the mass error data in term of variables with the maximum possible number of nuclei along each chain, the following transformation is employed

$$\stackrel{~}{A}=Int[\sqrt{2}(Nsin\theta +Zcos\theta )],\stackrel{~}{T}_z=Int[\sqrt{2}(Ncos\theta Zsin\theta )].$$
(3)
Both $`\stackrel{~}{A}`$ and $`\stackrel{~}{T}_z`$ are, by construction, integer numbers. To avoid introducing artificial noise, the data are softened by the interpolation of mass errors for unphysical values of $`\stackrel{~}{T}_z,\stackrel{~}{A}`$, i.e. those with $`\stackrel{~}{T}_z`$ even and $`\stackrel{~}{A}`$ odd, or viceversa. This process is necessary to eliminate the large number of zeros which are induced by the transformation, which create artificial high frequency noise in the data.

Definition (3) has the advantage that, for $`\theta =45^{}`$, $`\stackrel{~}{A}=A,\stackrel{~}{T}_z=2T_z`$. There are seven values of $`T_z`$ which have more than 40 nuclei with measured masses. They are $`2T_z=5,6,7,8,16,17,18`$. Their mass differences $`\mathrm{\Delta }M`$ are plotted as a function of the mass number $`A`$ in the seven inserts shown in Fig. 7. The clustering of negative and positive errors is evident, again exhibiting the presence of residual correlations between the mass errors and the atomic and neutron numbers (Z, N), or, equivalently, with the isospin projection and the mass numbers.

The squared amplitudes of the Fourier transforms of the mass differences are shown in Fig. 8. For the first four isotopic chains ($`2T_z=5,6,7,8`$) there is a prominent peak at $`f0.033`$, i.e. a period $`\mathrm{\Delta }A30`$, which is about one third the size of each set of nuclei. The remaining three chains exhibit also a peak at low frequencies, but lower and wider.

In order to underline the relevance of the lower frequencies in the mass error distribution, in Fig. 9 the mass errors are plotted as a function of A, for $`2T_z=6`$ (left), and 8 (right), in the top panel. The middle panel exhibits the distribution obtained including only the three low frequencies: $`f=0.011,0.022,0.033`$. The second is the most important, with a period $`\mathrm{\Delta }A46`$.

### II.3 Rotated $`\stackrel{~}{A}`$ and $`\stackrel{~}{T}_z`$

We found that the best orientation, in order to have as many isotopes as possible with the same $`\stackrel{~}{T}_z`$, is $`\theta =56^{}`$. With this transformation, there are 174 isotopes with $`\stackrel{~}{T}_z=0`$.

Figure 10 displays the mass errors for 7 values of $`\stackrel{~}{T}_z`$, from $`\stackrel{~}{T}_z=4`$ to 2. The regularities seen in Fig. 2 as regions with the same gray tone are seen here in the different plots, as groupings of nuclei with similar positive or negative mass differences, for the same $`\stackrel{~}{A}`$ region. Besides the two large groups with positive and negative mass errors below $`\stackrel{~}{A}=50`$, there are evident regions with negative errors close to $`\stackrel{~}{A}=100`$, and with positive mass differences for $`150<\stackrel{~}{A}<200`$.

A Fourier analysis of the previous results is presented in Fig. 11. It is clear that for all the chains a few low frequencies again dominate the spectrum. In some chains there are also some higher frequencies which seem to be relevant, while for $`\stackrel{~}{T}_z=0`$ their contribution is small. These frequencies, close to $`f=0.5`$, are associated to oscillations with period 2, i.e. strong fluctuations between one nucleus and its closest neighbors.

In Fig. 12 the contribution of very few frequencies to the mass differences, for the chains with $`\stackrel{~}{T}_z=2,1,0`$ and $`1`$ are presented. The upper panels show the full sequences of mass differences, the medium panels the contribution of the four frequencies with the largest amplitudes, and the lower panels the remaining mass errors. In all cases it is apparent that the regularities are quite significant, and that the substraction of these systematic errors would lead to important improvements in mass calculations.

### II.4 Bustrofedon: all data aligned

Plotting the mass differences for different Z, Fig. 13 top, and for different N, Fig. 13 bottom, is very common in mass calculations. Both plots exhibit some degree of structure. In this way we obtain a plot of mass differences as a function of Z, with all the isotopes plotted along the same vertical line, see Fig. 13. The difficulty in quantifying these regularities lies in the simple fact that the are many nuclei with a given N or Z. For this reason in the previous subsections we have analyzed the data using different cuts.

Another way to organize the FRDM mass errors for the 1654 nuclei with measured masses is to order them in a single list, numbered in increasing order. To avoid jumps, we have ordered the isotopes along a $`\beta o\upsilon \tau \rho o\varphi \eta \delta \stackrel{´}{o}\nu `$ (bustrofedon) line, which literally means ”in the way the ox ploughs”. Nuclei were ordered in increasing mass order. For a given even A, they were accommodated following the increase in N-Z, and those nuclei with odd A starting from the largest value of N-Z, and going on in decreasing order. The middle panel exhibits the same mass differences plotted against the order number, from 1 to 1654. It provides a univalued function, whose Fourier transform can be calculated. The squared amplitudes are presented in Fig. 14.

As it was the case for the cuts for given N, Z, A ot Tz, the fully ordered data Fourier amplitudes are completely dominated by a few low frequencies. It is also remarkable that for ln(f) $`>`$ -3 the frequency distribution clearly resembles a power law.

Removing the five frequencies with the largest Fourier components ($`|F_k|>3`$) generates a noisier pattern. Figure 15 shows the mass differences plotted as function of the order index, and the results of leaving only the frequencies with amplitudes $`|F_k|`$ larger (or smaller) than 3. There only five frequencies with Fourier amplitudes whose absolute value is larger than 3. From them the dominant frequency is the one with period $`\mathrm{\Delta }i=207`$. The regular pattern generated by these five frequencies in displayed in the middle panel of Fig. 15. The remnant, shown in the bottom panel, is clearly closer to white noise, while some bumps remain.

From the analysis presented in this section, we can conclude that there are conspicuous correlations in the FRDM mass differences, whose periodic character is clearly exhibited by the Fourier analysis of different cuts. It is worthwhile to study the slope of the squared amplitudes against the frequencies (i.e. the power spectra) in a log-log plot, which in a random matrix context was found to offer a signature of quantum chaos Rel02 . This is of great interest in itself and will be analyzed in a separate publication Vel03b .

## III Removing the regularities

Having established that there are patent regularities in the differences between the masses calculated using the FRDM Moll95  and the measured ones, we will proceed to eliminate them in a simple way, by removing the two frequencies which contribute the most to the mass errors as functions of N and Z. We introduce an amplitude and a phase for each frequency, having a total of six parameters for protons and six for neutrons. The functions which minimize the errors separately for protons and neutrons are:

$`\mathrm{\Delta }_1(N)=.30\text{sin}(2\pi .012N+.284)+.20\text{sin}(2\pi .047N+1.06)`$ (4)
$`\mathrm{\Delta }_1(Z)=.23\text{sin}(2\pi .020Z0.83)+.10\text{sin}(2\pi .053Z+3.74)`$ (5)
In a fit including all nuclei, we found

$`\mathrm{\Delta }_2(N)=.21\text{sin}(2\pi .011N+0.86)+.34\text{sin}(2\pi .049N+0.17)`$
$`\mathrm{\Delta }_2(Z)=.14\text{sin}(2\pi .025Z0.74).19\text{sin}(2\pi .075Z5.53)`$
$`\mathrm{\Delta }_2(N,Z)=\mathrm{\Delta }_2(N)+\mathrm{\Delta }_2(Z),`$ (6)
while including only nuclei with A $``$ 65 the best fit is

$`\mathrm{\Delta }_3(N)=.65\text{sin}(2\pi .009N+2.03)+.19\text{sin}(2\pi .047N+1.00)`$
$`\mathrm{\Delta }_3(Z)=2.90\text{sin}(2\pi .018Z3.51)+2.61\text{sin}(2\pi .019Z0.74)`$
$`\mathrm{\Delta }_3(N,Z)=\mathrm{\Delta }_3(N)+\mathrm{\Delta }_3(Z).`$ (7)
The $`\mathrm{\Delta }_1`$ are functions only of N or Z, and were adjusted using six parameters, optimized for the 1350 nuclei with A $``$ 65. $`\mathrm{\Delta }_2`$ and $`\mathrm{\Delta }_3`$ are obtained by including at the same time the corrections in Z and N, for all the nuclei in the first case, and for those with A $``$ 65 in the second. The fitted frequencies are all between 0.01 and 0.05, in agreement with the dominant Fourier frequencies found in the previous section.

The effect of removing these sinusoidal components as functions of N and Z is shown in Fig. 16, for those nuclei with masses larger than 65.

While the regular pattern is more apparent as a function of N, on the right hand side panels, in both cases the ‘corrected’ masses are more compact and compressed at smaller errors.

Obtaining a better fit of the known nuclear masses by the inclusion of six or twelve extra variables is by no means surprising. However, the use of sinusoidal functions of N and Z is strongly motivated by the data themselves. The fit does not only reduce the overall error, but makes the double peak in the error distribution disappear, as shown in Fig. 17 for all the nuclei, and in Fig. 18 for nuclei with A $``$ 65.

In both figures the error width reduction is evident, as well as the concentration of the mass differences around zero, with a single peak. The continuous curve represents a Gaussian curve with the same width and normalization. After removing the oscillatory components, the r.m.s. mass errors is reduced from 0.681 MeV to 0.597 MeV for all nuclei, and from 0.534 MeV to 0.462 MeV for nuclei with A $``$ 65.

## IV Final remarks

In the present study we have shown that the differences between the masses calculated using the FRDM of Möller at al Moll95  and the measured ones have a well defined oscillatory component as function of N and Z, which can be removed with an appropriate fit, significantly reducing the error width, and concentrating the error distribution on a single peak around zero.

The remaining correlations can only originate in the microscopic terms in the mass formula, which in the FRDM are evaluated using the Strutinsky method Bra72 . Having shown that these correlations have a simple and clear dependence in the proton and neutron numbers, we are studying the possible removal of these effects by a refinement of the Strutinsky method, whose results will be reported elsewhere Vel03b .

## Acknowledgements

Acknowledgements: Relevant comments by R. Bijker, O. Bohigas, J. Dukelsky, J. Flores, J.M. Gomez, P. Leboeuf, S. Pittel, A. Raga, P. van Isacker, H. Vucetich and A. Zuker are gratefully acknowledged. This work was supported in part by Conacyt, México.


text:
# The CAPMAP Instrument and its First Season

## 1 INTRODUCTION

The status of measurements of CMB polarization has evolved rapidly. Only a year separates the publication of tight upper limits Hedman et al. (2001, 2002); O’dell et al. (2002) and the detection of polarization Kovac et al. (2002). In the summer of 2001, we started designing the CAPMAP instrument. CAPMAP is a 16-element array comprised of twelve 90 GHz and four 40 GHz correlation receivers (four 90 GHz elements have been deployed for initial observations since February 2003). A new feed system has been designed to couple CAPMAP to the 7-meter Crawford Hill antenna to produce a 4 beam. The goal is to concentrate sensitivity on a $`1^{}`$ diameter cap around the North celestial pole (NCP). As indicated in Figure 7, CAPMAP will probe higher angular scales than those probed by MAP and DASI, and will be sensitive to the predicted peak in polarization. I will discuss in more details specifics of the telescope, the optics, and the receivers.

## 2 OPTICS

### 2.1 Telescope

In order to achieve high angular resolution, CAPMAP uses the Lucent Technologies 7-meter offset Cassegrain antenna located in Holmdel, NJ at the Crawford Hill Observatory. Built in 1977, it was originally designed to test polarization multiplexing of signals to geosynchronous satellites at K band and for radio astronomy up to 300 GHz. As a result, it demonstrates many of the requirements necessary for CMB measurements:

* low sidelobe level.
* low cross polarization level.
* good pointing accuracy and tracking capabilities on an alt-az mount.
* large f/D ratio, large useable focal plane.
* temperature-controlled Cassegrain focus room.

A detailed description of the antenna can be found in the Bell Labs Technical Journal Chu et al. (1978). Prior to the summer of 2001, when we started using the antenna, it had primarly been used during the last eight years to train winter-over operators for the south pole AST/RO telescope. We invested large amounts of work in the following two years to repair and upgrade the various mechanical, electrical, and software systems.

### 2.2 Feed system

If the full aperture were used, the 7-meter antenna would produce a $`0.02^{}`$ beam in the sky. This resolution is smaller than needed to probe the scales where the polarization peaks, at 10 or $`\mathrm{}`$ 1000. Thus, the feed system is designed to produce a FWHM beam size of 4 and to minimize pickup of stray radiation from the ground. This is achieved by aggressively under-illuminating both mirrors. The secondary and primary edge illuminations are respectively -56dB and -30dB. The small ($`3^{}`$) beam pattern required to produce such low edge tapers is the main price that must be paid for the advantage of the large effective focal ratio of $`f/D=5.6`$. A standard small-angle corrugated horn had to be rejected as a feed because its size would make it prohibitive to cryogenically cool and because its sidelobe pattern could not be controlled to the desired level.

The solution for the feed optics is a 3 cm aperture corrugated feedhorn followed by a 12 cm diameter HDPE meniscus lens shown in Figure 2. The lens is held in front of the horn with a metal shroud. The shroud is coated with a $`1/4^{\prime \prime }`$ thick microwave absorbing honeycomb to prevent stray radiation from illuminating the horn. This feed system produces the necessary beam size ($`2.6^{}`$ FWHM at 90 GHz) to achieve the designed secondary illumination and was optimized to keep its sidelobes below $``$ -40 dB (see Figure 3).

We note that because this antenna in not surrounded by ground screens, the extreme under-illumination and the obsessive control of sidelobes are necessary design features to keep the spillover below levels that would interfere with the observation of CMB polarization.

### 2.3 Array configuration

Because of its slow optics, the Crawford Hill telescope is well suited for an array of receivers at its focus. To determine the size of distortions for receivers displaced from the focal point, we performed extensive studies using both measurements from receivers at different positions in the focal plane and physical optics simulations Grasp8 Student Edition . Figure 4 (left) shows a simulated beam contour plot for a receiver located 39 cm above the focal point and pointed towards the center of the secondary. Distortions remain small even this far from the focal point.

The plate scale for the focal plane is such that a 15 cm shift in the focal plane corresponds to a $`0.25^{}`$ shift of the beam in the sky. For the initial deployment of four radiometers, we arranged the receivers in a 15 cm radius diamond pattern, obtaining the composite beam map in Figure 4 (right). This configuration allows us to observe a $`1^{}`$ diameter cap around the NCP by scanning the telescope back and forth in azimuth. This scan strategy naturally provide many systematic checks.

## 3 RECEIVERS

The receivers are heterodyne correlation polarimeters. With approximately 100K system noise and 14 GHz bandwidth, the receiver sensitivity is $``$1 mK$`\sqrt{s}`$. The signal is amplified by two back-to-back MMIC LNAs222Monolithic microwave integrated circuit low noise amplifiers(see section 3.1). After being down-converted to IF (2-18 GHz), filtered, and further amplified, the signal is detected with a correlation multiplier. The signal is phase-switched in the LO at 4 kHz to avoid the $`1/f`$ noise of the LNA, which have a $`1/f`$ knee before phase switching of $``$1 kHz. The IF is divided into three sub-bands to allow frequency discrimination. A single LO signal is generated and then distributed to each receiver after amplification in the MMIC power amplifiers.

### 3.1 MMIC LNAs

We take advantage of recent microwave technology developments by using a new generation of HEMT amplifier. We use MMIC LNAs and power amplifers developed, integrated, and tested at JPL (by Todd Gaier). We cascade two MMIC LNAs together to obtain the necessary gain. Figure 6 (right) shows the noise temperature and gain versus frequency for a typical MMIC cooled to 20K. Figure 3.1 (left) shows the distribution of noise and gains for the devices we have acquired so far. Note that 20dB of gain and less than 55K noise temperature are not uncommon.

### 3.2 Data Acquisition

The CAPMAP support electronics and data acquisition modules are all housed in the focal point cabin one meter away from the dewar. The vertex cab, support electronics, data acquisition module, and both dewar stages are temperature controlled. All such critical environmental and cryogenic temperatures are monitored and recorded at 1 Hz. All radiometer channels are synchronously sampled at 96 kHz using an ICS-610333http://www.ics-ltd.com/ data acquisition board. These data are then averaged down to 4 kHz, demodulated in software, and recorded at 100 Hz. For this season, the azimuth scan rate was one beam per 0.3 seconds.

## 4 CONCLUSION

CAPMAP complements MAXIPOL, B2K, CBI, DASI, and MAP in frequency coverage, experimental technique and angular resolution. An example of the sensitivity of CAPMAP to CMB polarization anisotropy is given in Figure 7. For this season, four receivers were deployed and observed during fair weather between January 16 and April 30, 2003 with sensitivities near the predicted ones based on 100 K system noise and 14 GHz bandwidth. As indicated in Figure 3, the beams are symmetric with acceptable sidelobes. Results from this first campaign will be detailed in a forthcoming publication. We expect to deploy the full 16 element array in early November 2003.

### 4.1 Acknowledgements

CAPMAP is a collaborative effort involving researchers at various institutions: Jeff McMahon, Phil Farese, Suzanne Staggs, and myself at Princeton University; Keith Vanderlinde, Colin Bischoff, Matthew Hedman, Dorothea Samtleben, and Bruce Winstein at University of Chicago; Eugenia Stefanescu and Josh Gundersen at University of Miami; and Todd Gaier at JPL.

This work was supported by the NSF grants #PHY-0099493 , #PHY-9984440, #PHY-0114422.


text:
# Discovery of a Clustered Quasar Pair at z ≈ 5: Biased Peaks in Early Structure Formation

## 1. Introduction

Quasars at large redshifts represent a powerful probe of structure formation in the early universe. Dissipative mergers and tidal interactions during the early stages of galaxy assembly may be fueling both bursts of star formation and early AGN activity (see, e.g., Silk & Rees (1998), Franceschini et al.  (1999), Monaco et al.  (1999), Kauffmann & Haehnelt (2000), Granato et al.  (2001), etc.). This fundamental connection is supported by the remarkable correlations between the masses of central black holes in nearby galaxies, and the velocity dispersions and luminosities ($``$ masses) of their old, metal-rich stellar populations (Magorrian et al.  (1998), Ferrarese & Merritt (2000), Gebhardt et al.  (2000)), and by the high metallicities observed in the high-$`z`$ quasar spectra (Hamann & Ferland (1999), and references therein). Quasars can also be used directly to probe evolution of large-scale structure out to high redshifts, as demonstrated clearly by the 2QZ survey (Croom et al.  (2001), Hoyle et al.  (2002)).

There should also be a fundamental connection between the formation of galaxies and the large-scale density field in the early universe. The highest density peaks, where presumably the first luminous objects formed, should be strongly clustered, due to biasing (Kaiser (1984)). This is a generic and robust prediction for essentially every model of structure formation, independent of any astrophysical details of galaxy formation.

Luminous high-$`z`$ quasars are likely situated in massive hosts (Turner (1991)). Such massive halos should be rare, and may be associated with $`4`$ to $`5`$-$`\sigma `$ peaks of the primordial density field (Efstathiou & Rees (1988), Cole & Kaiser (1989), Nusser & Silk (1993)). High-$`z`$ quasars can thus be used as biased tracers of the early large-scale structure, possibly marking the cores of future rich clusters.

A search for protoclusters around known high-$`z`$ objects such as quasars thus provides an important test of our basic ideas about the biased galaxy formation. We have conducted a search for clustered protogalaxies and AGN in the fields of selected quasars at $`z>4`$. Preliminary results from our program have been described by Djorgovski (1999) and Djorgovski et al.  (1999), and a complete account will be presented elsewhere.

Here we present the discovery of a clustered quasar pair at $`z5`$, which we interpret as a signature of a primordial large-scale structure, possibly a core of a forming rich cluster. This is the most distant large-scale structure currently known. One of the quasars, SDSS 0338+0021, was discovered by Fan et al.  (1999); the other was found in our much deeper search in its vicinity.

## 2. Observations and Data Reductions

Deep images of the field were obtained at the Palomar 200-inch Hale telescope using the prime-focus $`Cosmic`$ imager (Kells et al.  (1998)). The instrument FOV is 9.7 arcmin square, with 0.286 arcsec pixels. Multiple dithered exposures totalling 1800 s, 3000 s, and 7200 s were obtained in Gunn $`r`$, $`i`$, and $`z`$ filters respectively on 29 and 30 November 2000 UT, and of 3600 s, 1800 s, and 2400 s in Gunn $`g`$, $`r`$, and $`z`$ filters respectively on 31 December 2000 UT, all in good conditions. The data were reduced using standard procedures. Limiting magnitudes (3-$`\sigma `$ in a 3 arcsec apertures) are $`r_{lim}25.9`$, $`i_{lim}25.5`$, and $`z_{lim}23.1`$ mag in the SDSS ($`AB_\nu `$) photometric system. Multi-color selection was used to identify candidates for objects at $`z\text{ }>\text{ }5`$, as illustrated in Fig. 1. Details of the data reduction and candidate selection procedure will be presented elsewhere.

The first set of imaging data was used to select candidates for multi-slit spectroscopy, and spectra of several objects, including the known QSO (SDSS 0338+0021), were obtained at the WMKO Keck-I 10-m telescope on 29 December 2000 UT in good conditions. We used the LRIS instrument (Oke et al.  (1995)) with a 400 lines mm-1 grating ($`\lambda _{\mathrm{blaze}}=8500`$Å) and a GG495 order sorting filter, through 1.2 arcsec wide slitlets, with a mean dispersion in the wavelength region of interest of $`1.86`$ Å pixel-1. Two exposures of 1800 s each were obtained through a single slitmask at a PA = 340, close to the parallactic angle at the time. Exposures of arc lamps were used for wavelength calibration, with the flexure shifts corrected using the measurements of selected night sky lines. An average of archival response curves for this grating was used for the flux calibration.

Only one high-priority color-selected candidate, which we designated RD 657 (for “red dropout” and a serial number in our CCD object catalog), could be accommodated on this slitmask. It turned out to be a quasar at $`z5`$. The position of this object (J2000) is:

$$\alpha =03^h38^m30.03^s\delta =+00^{}18^{}40.4^{\prime \prime },$$

giving the projected separation from SDSS 0338+0021 of 196 arcsec. Setting the photometric zero-points on the magnitudes of SDSS 0338+0021 published by Fan et al.  (1999), i.e., $`r=21.70`$, $`i=19.96`$, and $`z=19.74`$, the magnitudes of the new QSO in the same SDSS system are $`r=23.01`$, $`i=21.24`$, and $`z=21.16`$ mag, with estimated uncertainties of $`0.1`$ mag, making it about 3.5 times fainter than the SDSS QSO. Fig. 2 shows the finding chart for the field, with both quasars marked.

Additional long-slit spectra of both objects were obtained using LRIS on 31 December 2000 UT in good conditions. We used a 600 lines mm-1 grating ($`\lambda _{\mathrm{blaze}}=7500`$Å) through a 1.0 arcsec wide slit, covering the wavelength range $`\lambda 49007450`$ Å, with a mean dispersion in the wavelength region of interest of $`1.255`$ Å pixel-1. Two exposures of 1200 s were obtained for SDSS 0338+0021, and four exposures of 1200 s for RD 657, with the object dithered on the slit between the exposures. The slit PA was always close to the parallactic angle. Arc lamps were used for wavelength calibration, with the flexure shifts corrected using the measurements of selected night sky lines, and exposures of standard star Hiltner 600 were used for the flux calibration.

We rebinned both grating data sets to a common sampling grid of 2 Å pixel-1, smoothed with Gaussians with $`\sigma =2`$ Å (lower than the optimal smoothing filter for these data, thus resulting in no loss of information). Fig. 3 shows the combined spectra of the two objects. The redshifts determined from the emission lines (taking into account the absorption of the blue side of Ly$`\alpha `$) are $`z_e=5.02\pm 0.02`$ for SDSS 0338+0021, and $`z_e=4.96\pm 0.03`$ for RD 657, with the large uncertainties due to the intrinsic difficulty of centering of broad emission lines. We also note the presence of a Lyman limit system at $`z_{LLS}=4.99`$ in the spectrum of the brighter QSO.

While absolute values of the redshifts cannot be measured very precisely, we used a simple cross-correlation to evaluate the redshift difference between the two objects. We excluded the portions of the spectra blueward of the centroids of the Ly$`\alpha `$ lines, since the differences in the IGM absorption between the two lines of sight could significantly affect the results. We obtain $`\mathrm{\Delta }z=0.063\pm 0.008`$, which corresponds to the restframe velocity difference of $`\mathrm{\Delta }v=3150\pm 400`$ km s-1 (these are the formal errors; possible systematic errors due to the mismatch of the spectra are hard to estimate precisely, but could be of the same order). Looking for correlated absorption systems in the two spectra would require data with a higher S/N and a higher resolution.

We obtained preliminary spectra of a number of other, fainter, color-selected candidates in the field. While the results are still inconclusive, none of them are luminous AGN.

## 3. Discussion and Conclusions

In what follows, for the sake of consistency with the previous work we will use the “standard quasar cosmology” with $`H_0=50`$ km s-1 Mpc-1, $`\mathrm{\Omega }_0=1`$, and $`\mathrm{\Lambda }_0=0`$.

In this cosmology, the observed angular separation of 196 arcsec corresponds to a projected spatial separation of $`1.126`$ proper Mpc, or $`6.744`$ comoving Mpc at $`z=4.99`$. If we assume that the observed $`\mathrm{\Delta }z=0.063\pm 0.008`$ is due entirely to cosmological expansion, the radial separation is $`4.27\pm 0.56`$ proper Mpc ($`25.6\pm 3.6`$ comoving). Quadratic sum of the two suggests a spatial separation of the two quasars of $`4.4\pm 0.6`$ proper Mpc ($`26.5\pm 3.6`$ comoving). However, this interpretation of the observed redshift difference as being due only to the Hubble expansion is uncertain. If we did not have any radial separation information, we could estimate the spatial separation from the projected separation alone: statistically, for a pair randomly oriented in space, the ratio of the two is $`\pi /2`$, so that the most probable spatial separation corresponding to our projected separation is $`1.77`$ proper Mpc ($`10.6`$ comoving). In the discussion below (which is similar to that by Schneider et al.  (2000)) we will consider both of these possibilities, i.e., comoving separations of $`26.5`$ Mpc (the cosmological expansion model) and $`10.6`$ Mpc (the deprojection model).

Individual quasars represent rare events in the general population of galaxies at any redshift; how likely is to have two of them so close together? In order to estimate the probability of finding such a close quasar pair by chance at this redshift, we use the evolving QSO luminosity function (QLF) by Fan et al.  (2001). For SDSS 0338+0021, we adopt the absolute blue magnitude $`M_B=26.56`$ from Fan et al.  (1999). Using the mean $`i`$ and $`z`$ magnitude difference of $`1.35`$ mag, we estimate $`M_B=25.2`$ for RD 657. The Fan et al.  (2001) QLF gives a number density of $`4.35\times 10^8`$ Mpc-3 for quasars with $`M_B25.2`$ at this redshift. The comoving volumes enclosed by spheres with radii equivalent to the physical separations of the quasars of 10.6 and 25.6 Mpc are $`4.97\times 10^3`$ Mpc3 and $`7.01\times 10^4`$ Mpc3, respectively.

The first question we can ask is, what is the a priori probability of finding such a close pair of QSOs at this redshift, regardless of the particular survey strategy? Assuming a Poissonian distribution of quasars, the probabilities of finding two QSOs at these luminosities in these volumes are $`2.3\times 10^8`$ and $`4.6\times 10^6`$, respectively. A similar reasoning was applied to the two serendipitously discovered quasar pairs at $`z>4`$ (Schneider et al. 1994a , Schneider et al.  (2000)).

However, the volume in which we found this pair was not selected at random: it is centered on a previously known QSO. We can thus ask an alternative question, which is specific for our experiment, namely, what is the probability that another QSO is found at random in this particular volume? (We note that the same answer would apply whether or not there is an already known QSO in its center.) The probabilities then become $`2.2\times 10^4`$ and $`3.0\times 10^3`$, respectively. Thus, it is still unlikely that this pair represents a chance occurence, suggesting that there is some physical clustering present.

We note that as of early 2003, we observed a total of 14 fields of QSOs at $`z\text{ }>\text{ }4.8`$ covering the FOV $`25`$ arcmin diameter ($`52`$ comoving Mpc in the cosmology used here) to a comparable depth. Our spectroscopic follow-up is still incomplete, and thus it would be premature to include this additional volume in the present computation, but to date no other cases of comparably bright QSO pairs have been found, suggesting that this system must be a relatively rare event. A full analysis and estimates of the QSO clustering and bias will be presented in a future paper, once the survey is complete.

We also note that in a deeper Keck survey of $`20`$ QSOs at $`z44.7`$, but covering a smaller FOV, $`6\times 8`$ arcmin2 ($`12\times 16`$ comoving Mpc2) we found at least 2 AGN companions to the known, bright QSOs, with sub-Mpc separations (not gravitational lenses), as well as a number of clustered faint galaxy companions (Djorgovski (1999), Djorgovski et al.  (1999), and in prep.).

The clustering strength cannot be meaningfully measured from a single pair of objects in a survey of as yet poorly defined coverage. With this caveat in mind, the small probability of a random occurrence of such a pair implies an effective $`r_0`$ which could be considerably greater than the observed pair separation, i.e., $`r_0>>10`$ comoving Mpc. At low redshifts, there is some spread or results between different groups (see, e.g., Boyle & Mo (1993), Croom & Shanks (1996), or Sabbey et al.  (2001)), but most authors find that the observed clustering length of quasars is comparable to that of galaxy groups, $`r_01020`$ Mpc (Bahcall & Chokshi (1991); see Hartwick & Schade (1990) for a review and references). A standard parameterization of the evolution of clustering in comoving coordinates is given by the formula

$$\xi (r,z)=\left(\frac{r}{r_0}\right)^\gamma (1+z)^{(3\gamma +ϵ)}$$

where $`\gamma 1.8`$ and $`ϵ`$ is the evolution parameter. For our chosen cosmology, the expected value for the CDM cosmogony is $`ϵ0.8`$, and this is consistent with observations of the evolution of galaxy clustering at $`z<1`$ (see, e.g., Carlberg et al.  (2000)). Thus, one expects a strong $`decrease`$ in the clustering strength at higher redshifts, and in any model gravitational clustering is always expected to grow in time. How do we then explain the apparent $`increase`$ in the strength of quasar clustering at high redshifts?

The most natural explanation is that quasars represent highly biased peaks of the density field, and that the bias itself evolves in time. Ever since the first detections of QSO clustering (e.g., Shaver (1984), Shanks et al.  (1987), Iovino & Shaver (1988), Mo & Fang (1993), etc.) it was considered possible that QSOs represent biased tracers of the density field, but the evolution of such bias was not clear. La Franca et al.  (1998) found a turn-up in the clustering strength of quasars even at redshifts as low as $`z2`$, but this was not confirmed in a much larger sample by Croom et al.  (2001). Outram et al.  (2003) find no evidence for an increase in the QSO clustering power spectrum amplitude out to $`z2.2`$.

The first hints of such an effect at high redshifts were provided by the three few-Mpc quasar pairs at $`z>3`$, found in the statistically complete survey by Schneider et al. 1994b , as first pointed out by Djorgovski et al.  (1993), and subsequently confirmed by more detailed analysis (Kundic (1997), Stephens et al.  (1997)). Intriguingly, the frequency of the few-Mpc separation quasar pairs at lower redshifts is roughly what may be expected from normal galaxy clustering (Djorgovski (1991); see also Zhdanov & Surdej (2001)). There is even a hint of a possible superclustering of quasars at $`z>4`$, on scales $`100h^1`$ comoving Mpc (Djorgovski (1998)), comparable to the scales of the first Doppler peak in CMBR fluctuations. Observations of large numbers of field galaxies at $`z33.5`$ also show a relatively strong clustering, with $`r_0510`$ Mpc, comparable to the galaxy clustering at $`z0`$ (Steidel et al.  (1998), Adelberger et al.  (1998)); this is also almost certainly a manifestation of biasing.

However, the bias should be even stronger at higher redshifts, and what is observed at $`z3`$ should be even more pronounced at $`z5`$: the earliest massive galaxies, including quasar hosts, should be strongly clustered. An example may be the possible grouping of Ly$`\alpha `$ emitters at $`z4.86`$ in the Subary Deep Field (Ouchi et al.  (2003), Shimasaku et al.  (2003)). Strong increase in biasing at high redshifts is also indicated in numerous theoretical studies, e.g., by Brainerd & Villumsen (1994), Matarrese et al.  (1997), Moscardini et al.  (1998), Blanton et al.  (2000), Magliocchetti et al.  (2000), Valageas et al.  (2001), Basilakos & Plionis (2001), etc. What these studies show is that a simple $`(r_0,ϵ)`$ parameterization of the clustering evolution is inadequate, and that the evolution of the bias factor, $`b`$, plays a key role. The effective bias factor generally increases with the redshift and the object mass (e.g., especially for the more luminous Lyman-break galaxies or the quasar hosts). For example, Croom et al.  (2002) find a marginally stronger clustering for the brighter QSOs, which might be residing in more massive hosts, and thus be more biased.

The chief uncertainty in our current understanding and interpretation of the structure evolution at high redshift, as indicated by luminous objects we can observe, is the evolution of bias. Observations of the clustering of quasars and galaxies around them at $`z45`$ and beyond can provide valuable empirical constraints in this endeavor. A better understanding of the primordial clustering of luminous sources is also important for models of the cosmic reionization (see Djorgovski et al.  (2001) and references therein), and thus for the interpretation of CMBR fluctuations at high angular frequencies. The quasar pair described here may be indicative of the biased clustering at $`z5`$, and more extensive and deeper surveys will provide additional observational input for the models of galaxy and large-scale structure formation.

We thank the anonymous referee for the insightful comments and suggestions. We wish to acknowledge the expert assistance of the staff of the Palomar and W. M. Keck Observatories. The LRIS data were obtained in the course of a collaborative project with F. Harrison and P. Mao. SGD acknowledges partial funding from the Ajax Foundation. The work of DS was carried out at the Jet Propulsion Laboratory, Caltech, under a contract with NASA.