format: text

text:
# Rare Events and Scale–Invariant Dynamics of Perturbations in Delayed Dynamical Systems

## Generic surface growth model for DDSs.–

Any DDS can be interpreted as a spatially extended system by introducing a simple transformation, $`t=x+\theta \tau `$, where $`x[0,\tau ]`$ is the space variable, while $`\theta `$ is a discrete time variable arecchi2 . This is a powerful representation in which the formation and propagation of structures, defects and spatiotemporal intermittency can be more easily identified arecchi2 . Note that the time delay $`\tau `$ becomes the system size in such a way that the time dependence with the delayed variable is transformed into a short-range interaction within the horizontal space coordinate $`x`$ in the space-time representation.

The dynamics of propagation of localized disturbances, i.e. the Lyapunov vector, can be obtained by linearizing (1) which, after the transformation to space-time representation, becomes politispace

$$\varphi (x,\theta )=u\varphi (x,\theta )+v\varphi (x,\theta 1),$$
(2)
where $`=_x`$, $`\varphi (x,\theta )=\delta y(x,\theta )`$ is the Lyapunov vector, and $`u=/y`$ and $`v=/y_\tau `$.

Following Pikovsky et. al. pik1 ; pik2 , we can map this equation into a surface growth model by introducing the surface height as $`h(x,\theta )=\mathrm{log}|\varphi (x,\theta )|`$. Then, after defining $`_\theta \varphi (x,\theta )=\varphi (x,\theta +1)\varphi (x,\theta )`$, one easily obtains the corresponding growth equation for the surface. We finally arrive to the equation,

$$_\theta h=(1/v)h+(u/v)+1,$$
(3)
where $`h=_xh`$ and the boundary conditions $`h(\tau ,\theta )=h(0,\theta +1)`$.

In the following we restrict ourselves to the most studied case of DDSs in the chaotic regime, in which $`u=\mathrm{const}`$ and $`v`$ is a nonlinear function (see Eq.(6) below), so that $`1/v`$ will play the role of a noise source. As we shall see below, probability distribution and correlation of this noisy term will turn out to be essential for understanding the behavior of the growth equation. Note that, if $`P_v(v)`$ is the probability density function (PDF) of $`v`$, the corresponding PDF of $`\zeta =1/v`$ is given by $`P_\zeta (\zeta )=|\zeta |^2P_v(1/\zeta )`$. Therefore, if $`P_v`$ behaves as $`P_v(v)|v|^\nu `$ when $`v0`$, we have that the PDF of $`\zeta `$ is then given by $`P_\zeta (\zeta )|\zeta |^{(2+\nu )}`$ for large $`|\zeta |`$ (conversely for small $`|v|`$). This power-law decay distribution of the noise is very remarkable and leads to a dynamics dominated by rare fluctuations. The occurrence of large events of amplitude $`|\zeta |1`$ has a small but non-negligible probability (they take place whenever the function $`v`$ gets close to zero). As we will show below, these rare fluctuations drive the system towards a stationary state that displays multiscaling properties. Anticipating part of our numerical results, in Figure 1 we plot the PDF of the chaotic function $`\zeta `$ for the Mackey-Glass model, which, as expected from our theoretical arguments, decays as a power-law $`P_\zeta (\zeta )|\zeta |^{2.42}`$.

The important fact that the surface growth equation that describes the propagation of disturbances in DDSs, Eq.(3), is dominated by rare fluctuations strongly suggests that the universality class of growth to which the surface must belong is the well-know linear Zhang model:

$$_\theta h=D^2h+f_0+\eta (x,\theta ),$$
(4)
which is the representative model of the universality class of linear surface growth in the presence of power-law noise zhang1 ; krug ; baraba ; lopez ; baraba&stanley . This model is expected to describe the dynamics of the Lyapunov vector in DDSs for large enough systems (i.e. long enough delay times). KPZ-like nonlinearities such as $`(h)^2`$ are not expected to appear in the effective surface growth equation, since Eq.(3) is statistically invariant under the $`hh`$ transformation and this symmetry must be preserved. Moreover, for the class of systems we are interested in this Letter, $`u(x,\theta )=\mathrm{const}`$ and a simple scaling argument shows that the additive noise in Eq.(4) must be $`\eta \sqrt{\zeta }`$. So, we expect the effective noise term $`\eta `$ in Eq.(4) to be distributed according to a power-law, $`P(\eta )\eta ^{(1+\mu )}`$, where the noise index is $`\mu =2(1+\nu )`$. Further justification that Eq.(3) is in the universality class of the linear Zhang model can be given as follows. By making use of standard stochastic functional techniques, like Novikov’s theorem ojalvo&sancho ; novikov , spatial average of the stochastic drift term $`\zeta h`$ can be calculated explicitely for the simpler case of Gaussian noise to obtain $`D^2h+bh`$, where the constants $`D`$ and $`b`$ are related to space-time noise correlator integrals (see for instance Ref. ojalvo&sancho ). In the case of noise with zero mean value and an even spatially homogeneous correlation function $`\eta (x,t)\eta (x^{},t^{})`$ one obtains a vanishing drift term $`b=0`$.

In the context of surface kinetic roughening, Eq.(4) corresponds to the well-known linear Zhang model zhang1 ; krug ; baraba ; lopez ; baraba&stanley  that describes the surface height dynamics in the presence of an uncorrelated power-law distributed noise source. In the following, we shall briefly review the scaling properties of this growth model. The height $`h(x,t)`$ exhibits scale-invariant properties characterized by the scaling behavior of the surface fluctuations as measured by the global width, $`W(L,t)=\overline{[h(x,t)\overline{h(x,t)}]^2}^{1/2}`$, where the average is calculated over all $`x`$ in a system of horizontal size $`L`$ (overline) and noise (brackets). Space and time scale-invariance leads to the scaling behavior of the width,

$$W(L,t)\{\begin{array}{ccc}t^\beta \hfill & \mathrm{if}& tL^z\hfill \\ L^\alpha \hfill & \mathrm{if}& tL^z\hfill \end{array}.$$
(5)
The critical exponents $`\alpha `$ and $`\beta `$ are the roughness and growth exponent, respectively, and characterize the global scaling behavior of the surface height. The dynamic exponent $`z`$ is related with the horizontal extent of the correlation length $`l_\times t^{1/z}`$ and fulfills the usual relation $`z=\alpha /\beta `$. A simple mean-field approximation krug ; baraba&stanley  provides estimates of the scaling exponents that depend on the index $`\mu `$ characterizing the power-law decay of the noise as follows. The dynamics is diffusive ($`z=2`$) for any value of $`\mu `$. However, the global roughness exponent depends on the noise index, $`\alpha =(3\mu )/\mu `$ for $`\mu 2`$, when the fluctuations are dominated by rare events, and $`\alpha =1/2`$ for $`\mu 2`$, where the noise distribution tail becomes irrelevant and standard Gaussian behavior is recovered krug ; baraba&stanley . Another interesting aspect of roughening dominated by power-law noise is multiscaling –namely, higher moments of the height-height correlation function do not scale with the same exponents– which indicates that the height distribution is far from Gaussian and dominated by rare fluctuations of the noise. In particular, multiscaling behavior in the linear Zhang model occurs for noise indexes in the range $`1<\mu <3`$ lopez2 . This behavior is analogous to that observed in the velocity fluctuations of turbulent fluids turbulence  as well as in some solid-on-solid models of epitaxial growth krug-mbe ; baraba&stanley . Multiscaling properties can be proven by calculating the $`q`$th order height-height correlation function $`G_q(x,t)=\overline{|h(x_0+x,t)h(x_0,t)|^q}^{1/q}`$, which scales as $`G_q(x,t)t^{\beta _q}`$ for short times and saturates, $`G_q(x)x^{\alpha _q}`$, for long enough time zhang1 ; krug ; baraba ; baraba&stanley . Exponents $`\alpha _q`$ and $`\beta _q`$ that depend on the index $`q`$ are the characteristic fingerprint of multiscaling behavior.

## Numerical results.–

We have carried out extensive numerical simulations of several chaotic delayed systems in order to compare with our theoretical findings. We have studied the well-known class of systems

$$\dot{y}=ay+F(y_\tau ),$$
(6)
which have been used in the past for a variety of applications ranging from biology mg  to lasers with feedback ikeda ; roy ; goedgebuer . For the particular choices $`F(\rho )=b\rho /(1+\rho ^{10})`$, $`F(\rho )=b\mathrm{sin}(\rho \rho _0)`$ and $`F(\rho )=b\mathrm{sin}^2(\rho \rho _0)`$, which correspond to the Mackey-Glass mg , Ikeda ikeda , and the delayed feedback model for secure optical communications goedgebuer , respectively. For this class of systems, the function $`u`$ is constant, $`u=a`$, and the additive noise in Eq.(4) is $`\eta (x,\theta )\sqrt{1/v}`$, where the nonlinear function $`v=dF/d\rho `$. In all our numerical simulations we have used the Adams-Bashforth-Moulton predictor-corrector scheme nr . From the numerical solution of the corresponding equation for the evolution of the perturbations of the leading Lyapunov exponent, Eq.(2), we have constructed the surface $`h(x,\theta )=\mathrm{log}|\varphi (x,\theta )|`$ following the procedure above discussed to get the space-time mapping. After that, the scaling properties of the resulting surface have been studied. For the sake of brevity we focus the discussion of numerical results on the Mackey-Glass model, but similar results were found for the other two systems studied. The parameters $`a=0.1`$ and $`b=0.2`$ are used in all the results we are presenting here, and simulation with a time delay varying from ten to a few thousand time units have been carried out. The region of interest here corresponds to delays $`\tau 17`$ for which the Mackey-Glass model is known to be hyper-chaotic farmer .

In Figure 1 (inset) we show a typical run of the simulation for the Mackey-Glass system with a time delay $`\tau =500`$ . The intermittent behavior of $`\zeta `$ is already apparent in the form of large peaks in the inset plot. We recorded the values taken by $`v`$ and collected statistics from long runs in order to obtain the PDF of the noise $`\zeta =1/v`$. In Fig. 1 we show our results for two instances of the time delay $`\tau =500,\mathrm{\hspace{0.17em}2000}`$. We obtain a power-law distribution $`P_\zeta (\zeta )|\zeta |^{2.42}`$ that is independent of the time delay used and corresponds to an exponent $`\nu =0.42\pm 0.02`$. The power spectrum of $`\zeta `$ (not shown) is flat indicating that $`\zeta `$ lacks long-range temporal correlations. We then expect the noise in the effective surface growth model (4) to be power-law distributed with $`\mu =2(1+\nu )=2.84\pm 0.04`$. This value of the noise index $`\mu >2`$ implies, according to the mean-field prediction for the linear Zhang model, that the roughness and growth exponent should be $`\alpha =1/2`$ and $`\beta =1/4`$, respectively. Very interestingly, it also implies that one should observe multiscaling behavior, since $`\mu (1,3]`$.

In Figure 2 we plot the time behavior of the surface width for different system sizes $`L`$ (i.e. time delays $`\tau `$) for the Mackey-Glass system. From that figure we obtain an estimation of the growth exponent $`\beta =0.25\pm 0.03`$ in the intermediate times regime, before a stationary state is reached. The scaling regime appears after an earlier short transient regime (of about one decade), and becomes longer for larger system sizes. For the largest system we used $`L=4\times 10^3`$ the scaling regime lasts for about two decades. After saturation, which occurs at times $`t_{sat}L^z`$, the surface height becomes scale-invariant and the width grows as $`W(tt_{sat},L)L^\alpha `$, with a roughness exponent $`\alpha =0.50\pm 0.05`$. From the inset of Fig. 2, it becomes clear that our estimation of the roughness exponent fits better in the limit of large system size, showing that long delays are required to reach the true scaling regime. These exponent values agree with the theoretical values predicted for the Zhang model for a noise index $`\mu =2.84>2`$.

Finally, we have also studied numerically the existence of multiscaling phenomena. In Figure 3, we plot the $`q`$th height-height correlation function for the Mackey-Glass system with time delay $`\tau =10^3`$. Multiscaling behavior is obtained for any time delay larger than $`\tau 17`$, which is the lower bound for chaos to appear in the system. However, multiscaling becomes clearer in larger systems (longer time delays), since the multiscaling region has a larger expand. As predicted by our theoretical arguments, multiscaling behavior is to be expected whenever $`\mu <3`$. Crossover to self-affine scaling is observed on large scales as generically expected for surface roughening driven by rare events (see baraba&stanley ; baraba  for details).

## Conclusions.–

Symmetry and scaling considerations have led us to conclude that the linear Zhang model for surface roughening driven by power-law noise describes the dynamics of the Lyapunov vector for long enough time delays. We then argue that DDSs generically display multiscaling of fluctuations due to the power-law nature of the driving noise.

The interpretation of DDSs as spatial systems has also been used in practical applications, in particular for the analysis of experimental data of $`CO_2`$ lasers with optical feedback arecchi2 . The technique allowed to identify many spacelike properties in experiments, including defects formation, transitions between weak and fully developed turbulence, etc. However, we have established a precise difference between high-dimensional chaos in delayed systems and in truly extended systems (like coupled-map lattices, partial differential equations, etc). Our results also suggest that different high-dimensional chaotic systems could be divided into a few universality classes based upon basic symmetries and conservation laws, akin to what occurs in scale-invariant surface growth.

###### Acknowledgements.

We thank L. Pesquera for useful comments and discussions. A.D.S. acknowledges a postdoctoral fellowship from CONICET (Argentina). Financial support from the MCyT (Spain) and FEDER under projects BFM2000-0628-C03-02, BFM2000-1108, BFM2001-0341-C02-02 as well as project OCCULT IST-2000-29683 from EU are acknowledged.


text:
# Lattice dilation near a single hydrogen molecule in an interstitial channel within a nanotube bundle

## I Introduction

The nature of adsorbed phases within a bundle of carbon nanotubes has been explored recently with many techniques hete ; oscdiff ; osche ; bob ; wang ; revaldo ; will ; kan&bill ; uptake ; boni ; rmp ; hecond ; dil ; interc . The varied motivations for this effort include interest in both basic scientific questions (e.g., novel phases of matter) and possible applications (gas storage, separation and sensing). Most of the theoretical and simulation research has assumed that the nanotubes are neither deformed nor displaced by the adsorbate. Such assumptions are convenient and might appear to be well justified by the small interaction between the gas and tubes (relative to the large cohesive energy of the tubes). However, our recent studies dil ; interc  found that the behavior of He, Ne and H2 gases within interstitial channels (ICs) contradicts that assumption. In particular, the ground state of the system corresponds to a slight expansion of the lattice of nanotubes in order to accommodate the gas (which is a moderately high density fluid at zero line pressure). The resulting expanded lattice provides a much larger binding energy for the particles because the interstitial potential of the gas is extremely sensitive to the tubes’ separation. For example, in the case of H2, the binding energy per molecule doubles when a 1% increase of lattice constant occurs. While the exact numbers found in that calculation are sensitive to the details of the potential, which is only approximately known, a large effect of dilation is expected to be present for all such potentials that have been proposed to describe these particles within ICs. The reason is that the particles experience strong forces in this environment, implying strong reaction forces acting on the neighboring tubes.

One interesting concomitant of this dilation was suggested in the previous study. The fact that the binding energy per molecule, $`|E_N|/N`$, of many molecules in a dilated lattice is much greater than that of a single molecule in an undilated lattice, $`|E_1^0|`$, suggests that there exists a strongly attractive, effective interparticle attraction, mediated by the lattice. Here, the superscript $`0`$ refers to the value in the absence of any dilation. If true, this “mediation hypothesis” implies that the molecules will condense to form an anisotropic liquid at and below a high critical temperature $`T_c`$ , an estimate of which is provided by the increase in binding energy due to the dilation:$`(|E_N||E_N^0|)/N300`$ K for H2. An analogy occurs in the BCS model of superconductivity, in which case the dynamical lattice deformation (phonons) mediates an electron-electron attraction, which drives the transition. The mediation hypothesis hinges on the implicit assumption that the strong binding found for the many-particle system is not present in the analogous one-body problem (including dilation). The heuristic idea underlying the hypothesis is simple and plausible: a single particle cannot deform its environment significantly, so only small dilation occurs in its vicinity. If several particles aggregate in the IC, instead, they can share the benefit of the lattice deformation without much additional cost in lattice energy. Thus, the large increase in binding is a cooperative, many-body effect and a high value of $`T_c`$ is expected.

The present paper is one contribution toward developing and testing this hypothesis. We compute the energy of a single H2 molecule in the presence of a deformed environment. The procedure employs a Born-Oppenheimer type of approximation, fixing the environment with a given deformation and finding the hydrogen eigenenergy in that environment. This is also analogous to the small-polaron theory employed to treat electrons in alkali-halide crystals kitt . It may therefore be appropriate to term the entity consisting of the hydrogen molecule and its surrounding dilation a “dilaton”. The net energy (tube lattice energy \+ H2 energy) is minimized with respect to possible deformations. We evaluate the ground state energy of the dilaton by solving the Schrödinger equation in the presence of an anisotropic lattice distortion, centered on the mean position of the impurity. We obtain a result that is consistent with the mediation hypothesis: the binding energy of a single molecule in the deformed lattice, $`|E_1|`$, is only slightly greater than that in the undeformed lattice, $`|E_1^0|`$. This means that the enhanced binding of the large $`N`$ system dil  is necessarily a consequence of the effective interaction between two or more dilatons. Hence, the hypothetical transition should indeed exhibit a high critical temperature.

The organization of this paper is the following. In section II, we describe the calculation of the elastic energy $`E_{latt}`$ of the deformed nanotube lattice. We compute the potential energy $`V`$ and the H2 energy $`E_1`$ inside the channel in Section III. Then, we find the ground state of the system (lattice + molecule) by minimizing the total energy $`E=E_{latt}+E_1`$ with respect to possible lattice distortions. Section IV summarizes our results and comments on its implications for the N-body problem.

## II Nanotube Bulge Energy

The interstitial channel between three parallel nanotubes will bulge beyond its equilibrium size to accommodate a molecule in the channel. Because of the rigidity of the tubes, the bulge will swell gradually along the channel axis and then contract again with some long length scale $`w`$ (Fig.1). The deformation will involve long wavelength acoustic modes of the nanotubes; it is energetically prohibitive for the unit cells of the tubes to undergo short wavelength, optical deformations. It is therefore appropriate to model the nanotube unit cells as rigid cylinders of length $`c=2.5`$ Å, so that each nanotube consists of a chain of cylinders secured end-to-end by spring forces Mizel . The forces within each chain of cylinders, as well as the forces between chains, can be obtained from the Born-von Karman force constants of reference  provided that nanotubes are assumed to be of the (10,10) variety.

Our ansatz for the shape of the bulge assumes that the three central nanotubes get pushed symmetrically outward a distance $`d(z)`$ by the molecule (Fig. 2). We adopt a Lorentzian form to describe the magnitude of the displacement away from equilibrium

$$d(z)=\frac{h}{1+(z/w)^2}.$$
(1)
Here, $`h`$ gives the maximum outward displacement of a unit cell while $`w`$ characterizes the extension of the bulge along the channel axis.

An internal energy $`E_{internal}`$ is required for the tubes to bend, which is determined by the frequency of transverse acoustic modes $`v_{TA}/c`$. Within the model of Ref. 16 adopted here, it has the value

$$E_{internal}=3\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}\frac{1}{2}M\frac{v_{TA}^2}{c^2}(d(l+1)d(l))^2$$
(2)
Here, $`v_{TA}2\times 10^6`$ cm/s is the velocity of the transverse acoustic modes that cause a tube to flex like a snake, while $`M=8\times 10^{22}`$ g is the mass of one nanotube unit cell, which contains 40 atoms Mizel . An inter-tube energy $`E_3`$ arises as the three tubes forming the channel spread apart. If they are members of a nanotube rope, the three tubes will push on $`9`$ neighoring tubes, contributing an additional inter-tube energy $`E_9`$, while the neighboring $`9`$ tubes themselves are assumed to remain in place. Utilizing the Lorentzian form (1) in the energy function of Ref. 16, one computes the deformation energy of the tubes as a function of the channel bulge shape:

$$E_{latt}(h,w)=E_{internal}+E_3+E_9$$
(3)
## III Hydrogen dilaton energy

As in most studies of hydrogen near nanotubes, we ignore the internal degrees of freedom of the molecule and focus on the motion of its center of mass. We compute the ground state energy $`E_1(h,w)`$ of the molecule in the deformed channel $`E_1(h,w)`$ by solving the Schrödinger equation

$$\left[\frac{\mathrm{}^2}{2m}^2+V(r,z)\right]\psi (r,z)=E_1\psi (r,z)$$
(4)
using the diffusion method. Here, $`V(r,z)`$ is the potential energy due to the deformed tubes. Since the deformation of the tubes is expected to vary over a long distance, we evaluate the potential using a local approximation (confirmed a posteriori) in which the potential is that of a uniform system having the local value of the function $`d(z)`$: $`V(r,z)U_{d(z)}(r)`$, where $`r`$ is the radial distance from the center of the channel. The potential energy $`U_{d(z)}(r)`$ is obtained by summing the contribution from the three tubes surrounding the channel and according to the nanotube potential model described in Ref.. In this way, we obtain the total energy as:

$$E(h,w)=E_{latt}(h,w)+E_1(h,w)$$
(5)
The ground state of the system is obtained by minimizing this expression with respect to the deformation parameters $`h`$ and $`w`$. Figure 3 shows the function $`E(h,w)`$ and the presence of a minimum $`E_{min}=282.5`$ K at $`h=0.005`$ Å and $`w=8`$ Å. This represents a binding energy increase of 1 K with respect to the case of non-deformed tubes. The deformation is extremely small in the transverse direction ($`h/a`$ 0.0003, where $`a=17`$ Å is the unperturbed lattice constant) and it extends over a much larger distance ($`16`$ Å) along the axis of the channel.

The probability density $`|\psi |^2`$ and the external potential $`V`$ for the minimum energy configuration are plotted in Figure 4. The molecule’s localization in the transverse direction is very similar to that found in the undilated case, with a root mean square transverse displacement $`\sqrt{r^2}0.2`$ Å.The localization in the direction parallel to the axis has a much larger spread of the wave function, $`\sqrt{z^2}1.6`$ Å. Note that the longitudinal spread of the wave function is much smaller than the spread of the deformation that is responsible for the localization. This is evident in Figure 4: the equipotential lines in Fig. 4(a) show more anisotropy than the constant probability lines of Fig. 4(b). The line shapes are consistent with a simple model calculation in which one assumes an anisotropic harmonic oscillator potential with radial and axial frequencies in the ratio $`\omega _r/\omega _z=150.`$

## IV Summary

We have investigated the ground state properties of the hydrogen dilaton, consisting of a localized molecular wave function plus the accompanying lattice deformation. As assumed from the start, the rigidity of the tubes implies that the longitudinal extent of the deformation is large compared to the other length scales in the problem. The ground state energy obtained by minimizing the total energy is just 1 K below the ground state energy of a molecule interacting with the rigid nanotube lattice; this value is the net result of a 10 K increase in lattice energy and an 11 K decrease in energy of the molecule due to the deformation. While the exact values of these energies are sensitive to the assumptions in the potential energy calculation, their small magnitude (compared to the binding energy of the many-body system) should be a robust conclusion of our study. This leaves open the question posed in the introduction: how does one account for the huge cohesive energy found in our previous study of the many-body system?

That large cohesive energy (about 200 K) is necessarily a consequence of inter-dilaton couplings: 2 body, 3 body,… One may speculate, for example, about whether a dilaton dimer, consisting of two molecules plus their deformation, has a much lower energy than two isolated dilatons. This might well be the case because a second molecule can exploit the rather extended deformation we have found for a single molecule. This argument would imply that the binding energy of the dimer is at least 10 K because the dimer need not produce a deformation that is longer than that shown in Fig. 4 for one molecule.

We intend in future work to explore this dimer problem by a technique similar to that used here, which is essentially a variational, strong-coupling method. If we find greatly enhanced cohesion in that problem, it will imply that the two-dilaton interaction is strongly attractive, with a resulting high temperature critical transition.

###### Acknowledgements.

We are grateful to Gerry Mahan for many helpful discussions. A. M. acknowledges the support of a 2001 Packard Fellowship for Science and Engineering. M.M.C and M.W.C. acknowledge support from the National Science Foundation.


text:
# Asymptotic Capture-Number and Island-Size Distributions for One-Dimensional Irreversible Submonolayer Growth

## Abstract

Using a set of evolution equations \[J.G. Amar et al, Phys. Rev. Lett. 86, 3092 (2001)\] for the average gap-size between islands, we calculate analytically the asymptotic scaled capture-number distribution (CND) for one-dimensional irreversible submonolayer growth of point islands. The predicted asymptotic CND is in reasonably good agreement with kinetic Monte-Carlo (KMC) results and leads to a non-divergent asymptotic scaled island-size distribution (ISD). We then show that a slight modification of our analytical form leads to an analytic expression for the asymptotic CND and a resulting asymptotic ISD which are in excellent agreement with KMC simulations. We also show that in the asymptotic limit the self-averaging property of the capture zones holds exactly while the asymptotic scaled gap distribution is equal to the scaled CND.

Recently, considerable theoretical effort has been carried out towards a better understanding of the scaling properties of the island-size distribution in submonolayer epitaxial growth. scaling ; BErapid ; re\_sc\_1d ; JMF\_recent ; APFsurf ; Vvedensky ; recent ; BE\_recent  For example, in the pre-coalescence regime the island size distribution $`N_s(\theta )`$ (where $`N_s`$ is the number of islands of size $`s`$ at coverage $`\theta `$) satisfies the scaling form, scaling

$$N_s(\theta )=\frac{\theta }{S^2}f\left(\frac{s}{S}\right),$$
(1)
where $`S`$ is the average island size, and the scaling function $`f(u)`$ depends on the critical island size and on the island morphology.scaling

One of the standard tools used in these studies is the rate-equation (RE) approach which involves a set of deterministic coupled reaction-diffusion equations describing the coverage-dependence of $`N_s(\theta )`$ through a set of rate-coefficients usually called capture numbers. Smo ; Venables  For the irreversible growth of point islands, rate-equations valid in the pre-coalescence regime may be written in the form

$$\frac{dN_1}{d\theta }=12R\sigma _1N_1^2RN_1\underset{s2}{}\sigma _sN_s$$
(2a)
$$\frac{dN_s}{d\theta }=RN_1\left(\sigma _{s1}N_{s1}\sigma _sN_s\right),\text{for}s2$$
(2b)
where the capture numbers $`\sigma _s`$ ($`\sigma _1`$) correspond to the average capture rate of diffusing monomers by islands of size $`s`$ (monomers) and $`R=D/F`$ is the ratio of the monomer diffusion rate to the deposition rate. Accordingly, the central problem in using the RE approach is the determination of the average capture numbers $`\sigma _s(\theta )`$ and the corresponding capture number distribution (CND).

Recently, we have developed a self-consistent rate-equation approach to irreversible submonolayer growth in which correlations between the size of an island and the corresponding average capture zone are explicitly taken into account in order to accurately predict the scaled island size and capture number distributions.JMF\_recent ; APFsurf  Our method involves numerical integration of the island-density RE’s (2) along with the analytical solution of a set of approximate evolution equations for the average Voronoi-area surrounding an island of size $`s`$, and is based on the following two assumptions: (1) the average capture area per freshly nucleated dimer is proportional (before rescaling) to the average area per island and (2) the combined effects of the preferred nucleation of dimers in large capture zones and the preferred “break-up” of large capture zones due to nucleation may be approximated by a uniform rescaling of the average capture zone of each island. JMF\_recent ; APFsurf  Using this approach, we have obtained numerical results for the scaled capture number and island-size distributions which agree well with kinetic Monte Carlo (KMC) simulations in both one and two dimensions over a wide range of experimentally relevant values of $`D/F`$ ($`D/F=10^510^9`$).

Recently, it has been argued Coment03  that because our method does not explicitly take into account spatial fluctuations in the nucleation, then it must lead to a diverging ISD in the asymptotic limit corresponding to infinite $`D/F`$. However, as shown by Bartelt and Evans,BErapid  in the asymptotic limit the scaled ISD is related to the scaled CND as,

$$f(u)=f(0)\mathrm{exp}\left[_0^u𝑑x\frac{2z1C^{}(x)}{C(x)zx}\right],$$
(3)
where $`C(s/S)=\sigma _s/\sigma _{av}`$ is the scaled CND, $`z`$ is the dynamical exponent describing the dependence of the average island size on coverage ($`S\theta ^z`$), and $`f(0)`$ is determined by the normalization condition,

$$_0^{\mathrm{}}𝑑uf(u)=1.$$
(4)
As pointed out in Ref. BErapid, , Eq. 3 implies that if $`C(u)>zu`$ then no divergence will occur. However, if $`C(u)`$ crosses $`zu`$ at some value $`u_c`$ then the ISD will be cut off at $`u_c`$ (i.e. $`f(u)=0`$ for $`uu_c`$) if $`C^{}(u_c)>2z1`$ while a divergence in the ISD will occur if $`C^{}(u_c)<2z1`$. An example of a divergent asymptotic ISD is the usual mean-field theory with $`C(u)=1`$. Thus the question of whether or not our method leads to a divergence in the asymptotic limit is entirely determined by the asymptotic scaled CND.

In this Rapid Communication we rigorously address the question of the asymptotic behavior obtained using our method by analytically deriving the asymptotic scaled CND along with the corresponding asymptotic scaled ISD for the case of irreversible growth of point islands in one-dimension. We find that, contrary to the claims in Ref. Coment03, , our method leads to a non-divergent, asymptotic ISD as well as to an asymptotic scaled CND which is close to that obtained in simulations. We then show that by slightly modifying our original analytical form for $`C(u)`$, an improved analytical expression for the asymptotic scaled CND may be obtained which is in excellent agreement with KMC simulations. The resulting analytic expression leads to a scaled ISD which is also in excellent agreement with KMC simulations. We also demonstrate that the asymptotic scaled gap distribution is identical to the scaled CND and as a result the self-averaging property of the capture zones holds exactly.

For clarity, we briefly review our method and its application to the case of irreversible growth of point islands in one-dimension.APFsurf  In this case, we have shown that the local capture-number $`\stackrel{~}{\sigma }(y)`$ for an island with gap-size $`y`$ (corresponding to the distance to the nearest island) is

$$\stackrel{~}{\sigma }(y;\theta )=\frac{2\xi _1}{\xi ^2}\mathrm{tanh}\left(\frac{y}{2\xi _1}\right),$$
(5)
where the monomer capture length $`\xi `$ and the nucleation length $`\xi _1`$ are defined as

$$2\sigma _1N_1=1/\xi _1^2,1/\xi _1^2+\underset{s2}{}\sigma _sN_s=1/\xi ^2.$$
(6)
From the evolution equations for the distributions of gap lengths, the average “gap-length” $`\stackrel{~}{y}_s`$ (before rescaling) corresponding to an island of size $`s`$ may be obtained as the solution of the equation,

$$s2=_{\theta _y}^\theta RN_1(\varphi )\stackrel{~}{\sigma }(\stackrel{~}{y}_s;\varphi )𝑑\varphi .$$
(7)
The coverage $`\theta _y`$ is defined by $`\stackrel{~}{y}_s=bY(\theta _y)`$ where $`b`$ is the proportionality factor (before rescaling) between the average gap-length of a freshly nucleated dimer and the average gap-length of all islands, $`Y(\theta )=1/N(\theta )`$ is the average gap-length at coverage $`\theta `$, and $`N(\theta )=_{s2}N_s(\theta )`$ is the average island density. Physically, the integral in Eq. 7 may be interpreted as corresponding to the average number of particles added to the dimer since it was formed at coverage $`\theta _y`$, neglecting any change in the capture zone due to nucleation and break-up. To include the effects of break-up, the gap lengths are rescaled by a rescaling factor $`a=1/_{s2}N_s\stackrel{~}{y}_s`$ so that the average gap length for an island of size $`s`$ is given by $`y_s=a\stackrel{~}{y}_s`$ while the capture number is given by $`\sigma _s=\stackrel{~}{\sigma }(y_s)`$.

We now consider the asymptotic limit corresponding to infinite $`D/F`$. In this limit, and assuming that $`\theta >>\theta _x`$ (where $`\theta _xR^{1/3}`$ corresponds to the coverage at which the island-density equals the monomer density), Eq. 5 may be rewritten as,

$$\stackrel{~}{\sigma }(y;\theta )y/\xi ^2.$$
(8)
This implies that in the asymptotic limit the scaled gap-distribution $`B(s/S)=y_s/Y`$ is identical to the corresponding scaled capture number distribution $`C(s/S)`$. This result also demonstrates the “self-averaging property” of the capture zones, i.e. in the asymptotic limit the average capture number $`\sigma _s`$ is exactly equal to the local  capture number evaluated at the average capture zone for an island of size $`s`$.

In the asymptotic limit that $`\theta >>\theta _x`$ one also has $`dN_1/d\theta =1RN_1/\xi ^20`$, which implies that $`RN_1/\xi ^2=1`$. Using these results in Eq. 7 and assuming that also $`\theta _y>>\theta _x`$,largey  we obtain

$$s2=\stackrel{~}{y}_s(\theta \theta _y).$$
(9)
Since in the asymptotic limit $`S1`$, this may be rewritten as

$$u=s/S=(\stackrel{~}{y}_s\theta /S)(1\theta _y/\theta ).$$
(10)
In the asymptotic limit, one has $`\theta /S=N=1/Y`$. Using $`\stackrel{~}{y}_s=b/N(\theta _y)`$ and $`N(\theta _y)\theta _y^{1/4}`$, $`N(\theta )\theta ^{1/4}`$ (since $`\theta _y,\theta \theta _x`$ and $`z=3/4`$ for irreversible growth of point-islands in one-dimension), one has $`(\theta _y/\theta )=(bY/\stackrel{~}{y}_s)^4`$. Thus, the asymptotic gap-distribution before rescaling will satisfy,

$$u=\stackrel{~}{B}(u)b^4/\stackrel{~}{B}(u)^3,$$
(11)
where $`\stackrel{~}{B}(u)=\stackrel{~}{y}_s/Y`$. Rescaling, we obtain for the asymptotic scaled gap distribution $`B(u)=a\stackrel{~}{B}(u)`$ the equation

$$u=B(u)/ab^4[B(u)/a]^3,$$
(12)
where the rescaling factor $`a`$ and the proportionality factor $`b`$ must satisfy the requirement that

$$_0^{\mathrm{}}B(u)f(u)𝑑u=1$$
(13)
It is easy to see that for any positive constants $`a`$ and $`b`$, Eq. 12 has a unique positive solution for every positive $`u`$ while for $`az`$ one has $`B(u)>zu`$. Similarly, one may show that for $`a<2/3`$ the ISD predicted by Eq. 3 will diverge while for $`2/3<a<z`$ there will be a cut-off at $`u_c`$, i.e. $`f(u>u_c)=0`$.

We now consider the case $`b=1`$ as was assumed in our numerical calculations for finite $`D/F`$.JMF\_recent ; APFsurf  Eq. 3 along with the requirement that $`f(u)0`$ rapidly for large $`u`$, indicates that the asymptotic $`B(u)`$ should satisfy $`B(u)zu`$ for large $`u`$. BE\_recent  Applied to Eq. 12 this implies that $`a=z`$ and thus $`C(0)=B(0)=3/4`$, in very good agreement with our KMC simulation results. This leads to the following expression,

$$u=\widehat{B}(u)1/\widehat{B}(u)^3$$
(14)
where $`\widehat{B}(u)B(u)/z`$ and $`\widehat{B}(0)=1`$. Unfortunately, integrating Eq. 3 numerically using Eq. 14 to determine $`f(u)`$, we find that the resulting “cross-integral” $`_0^{\mathrm{}}B(u)f(u)𝑑u1.15`$ is somewhat larger than $`1`$. This implies that a smaller value of the rescaling factor $`a`$ must be used. By carrying out the appropriate numerical integrals we find that $`a0.70`$ satisfies the “cross-integral” normalization condition (13). The corresponding scaled CND and ISD with $`b=1`$ are shown in Fig. 1. As can be seen, the calculated asymptotic ISD is shifted to the right from the simulation result and is cut off to zero at $`u_c1.85`$ where $`C(u)`$ crosses $`zu`$. Thus we find that, contrary to the claims in Ref. Coment03, , our method leads to a non-divergent ISD in the asymptotic limit of infinite $`D/F`$.

It is also interesting to consider the general case $`b1`$, i.e. the capture-area of a freshly-nucleated dimer is proportional to the average area per island with some unknown proportionality constant $`b`$. In this case, we assume that the value of the rescaling constant $`a`$ is fixed to the value $`a=z`$ by the requirement that the asymptotic behavior of $`B(u)`$ is given by $`B(u)zu`$ for large $`u`$. We then search for a value of $`b`$ ($`b0.87`$) such that the sum-rule (13) is satisfied. As can be seen in Fig. 1, the resulting scaled gap-distribution $`B(u)`$ is now very close to the simulation results for all $`u`$ except for small $`u`$ ($`u<0.7`$) where it is now significantly lower than for the case $`b=1`$. Accordingly, the peak of the corresponding island-size distribution (Fig. 1(b)) is somewhat lower than the peak of the simulated distribution, and the distribution itself is somewhat wider.

We now present an analytical form for the capture-number distribution $`C(u)=B(u)`$ which provides excellent agreement with simulations and thus strongly supports our conjecture that $`B(0)=z`$ and that $`B(u)zu`$ for large $`u`$. Since Eq. 12 with $`b=1`$ and $`a=z`$ leads to an expression for $`B(u)`$ (14) which satisfies $`\widehat{B}(0)=1`$ and $`B(u)zu`$ for large $`u`$, but a cross-integral (13) which is just slightly larger than $`1`$, we propose that the correct analytic expression for $`B(u)`$ has the same form as (14) but with additional terms corresponding to higher order powers of $`[\widehat{B}(u)]^1`$. As an example, we chose the following form,

$$u=\widehat{B}(u)(b_2/\widehat{B}(u)^3)\mathrm{exp}[b_1/\widehat{B}(u)^n],$$
(15)
where $`b_2=e^{b_1}`$ in order to satisfy $`\widehat{B}(0)=1`$, and $`n`$ is a free parameter for which we will choose $`n=5`$. We then use Eq. 3 to calculate $`f(u)`$ and search for a value of $`b_1`$ ($`b_10.993`$) such that the sum-rule (13) is satisfied. Our results for the corresponding $`B(u)`$ and $`f(u)`$ obtained in this manner are shown by the solid curves in Fig. 1. As can be seen, the agreement between the calculated scaled gap-distribution $`B(u)`$ and simulations is now excellent. Accordingly, the agreement between the corresponding calculated ISD and simulations is also very good. This result also confirms that $`C(u)=B(u)`$ in the asymptotic limit.

We have also tried several other values of $`n`$ ($`n=3,4,6`$) and obtained similarly good agreement with simulations.agreement  Finite expansions in powers of $`[\widehat{B}(u)]^1`$ satisfying $`\widehat{B}(0)=1`$ also lead to reasonably good agreement. In particular, the expression

$$u=\widehat{B}(u)\gamma \widehat{B}(u)^3\frac{1\gamma }{2}[\widehat{B}(u)^9+\widehat{B}(u)^{12}],$$
(16)
(with $`\gamma 0.39`$ to satisfy (13)) leads to results which are essentially identical to those obtained using Eq. 15. Thus, it would appear that almost any form similar to (14), but with higher order corrections, which satisfies the condition $`\widehat{B}(0)=1`$ as well as the sum-rule (13), leads to reasonably accurate results for the asymptotic scaled ISD and CND.

In conclusion, we have derived analytical expressions for the asymptotic scaled capture-number, gap-, and island-size distributions for the case of one-dimensional irreversible growth of point islands. In particular, we have shown that our method leads to an asymptotic scaled CND which is close to that obtained in simulations and as a result to a non-divergent asymptotic ISD. Our analytical results also indicate that in the asymptotic limit the scaled gap-distribution $`B(u)`$ is identical to the scaled capture-number distribution $`C(u)`$. Finally, by slightly modifying our analytical expression (14) for the scaled CND and solving self-consistently, we have obtained analytic expressions for the scaled CND and corresponding ISD which are in excellent agreement with simulations.

We note that our analytical results also demonstrated that the “self-averaging property” of the capture zones (i.e. the average capture number $`\sigma _s`$ is exactly equal to the local capture number evaluated at the average capture zone of islands of size $`s`$) is exact in the asymptotic limit for the case of irreversible growth of point-islands in one-dimension. If this “self-averaging property” holds generally (as suggested by the good agreement obtained in our numerical results for finite $`D/F`$ in both one- and two-dimensions APFsurf ; JMF\_recent ) then it may be possible in general to obtain accurate asymptotic ISDs without having to know the exact distribution of areas or the redistribution of areas due to nucleation events. This represents a significant simplification of the calculation of ISDs since the determination of the full distribution of capture zones is a very difficult problem for which no exact solution seems to be possible even in the most simple cases.BE\_recent ; recent

Finally, we consider the extension of these results to other cases of interest such as growth in two-dimensions. The present analytical work was essentially dependent on the asymptotically valid relation $`B(u)=C(u)`$. It is not clear that such an explicit relation may be derived for the case of two-dimensional growth, and additional mathematical difficulties occur due to complicated expressions for the local capture numbers. As a result, it is not clear if a straightforward extension of the present work to this case is possible, and further work is required to answer this question.

J. G. Amar would like to acknowledge support from the Petroleum Research Fund and from the NSF (DMR-0219328).


text:
# Superwind-driven Intense H2 Emission in NGC 6240 II: Detailed Comparison of Kinematical and Morphological Structures of the Warm and Cold Molecular Gas

## 1 INTRODUCTION

NGC 6240 is a nearby111In this paper, the distance is assumed to be 98 Mpc, following, for consistency, Heckman, Armus, & Miley (1987) and Ohyama et al. (2000). luminous infrared galaxy (e.g., Wright, Joseph & Meikle 1984; Joseph & Wright 1984), comprising two merging nuclei (north and south nuclei: hereafter, N and S nuclei, respectively) separated by 1.9″ (Condon et al. 1982; Fried & Schulz 1983; Eales et al. 1990; Thronson et al. 1990; Herbst et al. 1990; Colbert, Wilson, & Bland-Hawthorn 1994; Sugai et al. 1997; Tacconi et al. 1999; Ohyama et al. 2000; Scoville et al. 2000; Tecza et al. 2000; Beswick et al. 2001), and is well known for its unusually luminous H2 emission in the near infrared (e.g., Rieke et al. 1985; DePoy, Becklin, & Wynn-Williams 1986). Although many excitation mechanisms for the H2 have been proposed to date \[shock excitation (Rieke et al. 1985; DePoy et al. 1986; Lester, Harvey, & Carr 1988; Elston & Maloney 1990; Herbst et al. 1990; van der Werf et al. 1993; Sugai et al. 1997); UV fluorescence (Tanaka, Hasegawa, & Gatley 1991); X-ray heating (Draine & Woods 1990); and formation pumping (Mouri & Taniguchi 1995)\], H2 line ratio analyses have revealed that shock excitation is responsible for almost all the H2 emission (Sugai et al. 1997; Ohyama et al. 2000; Tecza et al. 2000). However, because of its huge H2 luminosity and exceptionally large H2/Br $`\gamma `$ intensity ratio, the origin of NGC 6240’s H2 emission has been the subject of considerable debate; star formation activity is not strong enough to power all the observed H2 luminosity (e.g., Rieke et al. 1985; Draine & Woods 1990). Interestingly, most of the H2 emission comes from the region between the nuclei, although both Br $`\gamma `$ and \[Fe II\] (originating from OB type stars and supernovae in star-forming regions) are detected from each nucleus (Herbst et al. 1990; van der Werf et al. 1993; Sugai et al. 1997; Ohyama et al. 2000; Tecza et al. 2000). Two possible models have been proposed to explain this intense, off-nuclear, shock-excited H2 emission. In the first, the emission is caused by shock heating at the interface between the two colliding galaxies (Tecza et al. 2000 and references therein), while in the second it is due to the interaction between the hot gas associated with the superwind activity of the S nucleus and the tidally-produced molecular gas concentration between the two nuclei (Ohyama al. 2000).

Recent high-resolution, sensitive CO observations of NGC 6240 at radio wavelengths have revealed that a significant fraction of the system’s molecular gas is concentrated not at the two nuclei but rather occupies the region between them, forming a massive thick disk with highly turbulent motion (Tacconi et al. 1999). Hereafter we refer to the CO gas between the nuclei as the off-nuclear molecular gas or molecular gas concentration. A violent tidal disturbance during the merging process of the two nuclei could be responsible for the formation of such an off-nuclear molecular gas disk, although the creation of such a massive concentration of gas away from the merging nuclei has no theoretical explanation (e.g., Barnes 2002 and references therein; see also Tecza et al. 2002). The peak of the H2 emission is located between the S nucleus and the CO peak (Ohyama et al. 2000; Tecza et al. 2000).

The first step in exploring the reason for the source’s unusually intense H2 emission was the analysis of the velocity field and the excitation conditions of the H2 lines (Ohyama et al. 2000; Tecza et al. 2000 and references therein). In this paper we take the second step by conducting a detailed comparison of the morphological and kinematical conditions of H2 and CO lines. This work is motivated by the fact that H2 and CO are molecular species with different excitation conditions; comparisons of their emission lines should give us some new insights on the conditions under which high H2 emitting efficiency (H2 flux per molecular gas content) is likely. Such a comparison has only recently become possible, as a result of the availability of sub-arcsecond, sensitive mapping observations of both CO and H2. Tecza et al. (2000) conducted three-dimensional H2 mapping observations of this galaxy under sub-arcsec seeing conditions, while Tacconi et al. (1999) carried out high-resolution CO ($`J=21`$) mapping observations with a sub-arcsec beam. Both data have been re-sampled to make a direct comparison between the cold and warm molecular gas at the same scale (Tecza et al. 2000). Based on these data, we discuss how H2 emission can be excited efficiently, focusing on the H2/CO intensity ratio as an indicator of the H2 emitting efficiency.

## 2 Data and Analysis

All data that we analyze in this paper were reproduced from Figure 3 of Tecza et al. (2000), in which grids of both the CO ($`J=21`$) and H2 velocity profiles were presented. These profiles were reproduced from the original data of both Tecza et al. (2000) and Tacconi et al. (1999) for H2 and CO, respectively. The velocity resolution of the original H2 data was 150 km s-1, while the spatial resolution was 0.8″ - 1.0″ , the seeing conditions at the time of observation (Tecza et al. 2000). In the CO observations, the full width at half maximum (FWHM) of the elliptical Gaussian was beam 0.9″$`\times `$ 0.5″ , while the velocity resolution was 6.5 km s-1. These data were re-sampled onto $`0.6\times 0.6`$ arcsec2 sub-regions ($`290\times 290`$ pc2) around the double nuclei over $`3.6\times 4.2`$ arcsec2 ($`1.7\times 2.0`$ kpc2) regions, and re-sampled at 140 km s-1 velocity resolution (Tecza et al. 2000). These data enabled us to compare the CO and H2 velocity profiles directly, with a resolution that was sufficient to separate the two merging nuclei and the surrounding regions.

The velocity profiles were integrated over blue and red velocity regions ($`dVVV_{\mathrm{sys}}=3000`$ km s-1 and $`dV=0300`$ km s-1 for blue and red, respectively, where $`V_{\mathrm{sys}}`$ is the systemic velocity) in each sub-region to measure the blue and red intensities of CO and H2 \[$`I_{\mathrm{blue}}`$(CO), $`I_{\mathrm{red}}`$(CO), $`I_{\mathrm{blue}}`$(H2), and $`I_{\mathrm{red}}`$(H2)\]. The velocity cut-offs at $`dV=\pm 300`$ km s-1 were applied in the integration because it is necessary to measure the intensities in the same velocity range for all sub-regions in order to make a fair comparison of their spatial distribution. We found that some sub-regions had a weak CO intensity at larger velocity offsets from the systemic velocity, so that the intensity could not be accurately measured out to the velocity cut-off. We eliminated the integrated intensity data of the sub-regions from the following analyses if the normalized intensity of CO (normalized to be unity at the CO peak position: see Figure 3 of Tecza et al. (2000) and Figure 2 of this paper) was 0.05 (the typical noise level) or less at the cut-off velocity. The cut-off velocity ($`dV=\pm 300`$ km s-1) was chosen as a compromise: although a larger cut-off velocity is desirable for studying the properties of the high-velocity gas components, applying a larger ($`dV>300`$ or $`<300`$ km s-1) velocity cut-off would discard a larger number of sub-regions due to a lack of useful integrated intensity data. Note that since the observed line widths ($`1000`$ km s-1 and $`1600`$ km s-1 for CO and H2 at full width at zero intensity, respectively: Tacconi et al. 1999; Tecza et al. 2000) are wider than the applied velocity cut-off (at $`dV=\pm 300`$ km s-1), about 19% and 33% of the total CO and H2 fluxes fall outside the velocity range for the integration, and were not counted in the analysis that follows. Single Gaussian profile fitting was used to measure the velocities of both lines \[$`V`$(CO) and $`V`$(H2)\], although some areas show more complicated emission line profiles, such as a double-peaked profile (Tecza et al. 2000). Figure 1 shows the total (blue+red) integrated intensity and the velocity maps of both CO and H2 \[$`I_{\mathrm{blue}}`$(CO)$`+I_{\mathrm{red}}`$(CO), $`I_{\mathrm{blue}}`$(H2)$`+I_{\mathrm{red}}`$(H2), $`V`$(CO), and $`V`$(H2)\]. The peaks of both lines are located between the two nuclei and slightly east of the S nucleus ($`0.6`$″) 222 Note that the higher resolution H2 map presented by Tecza et al. (2000) and other previous H2 maps (e.g., Sugai et al. 1997) showed the H2 peak between the double nuclei, although our Figure 1 shows the peak at a slightly different position. This is an artefact that was probably caused by the re-sampling of the data onto the $`0.6\mathrm{}\times 0.6\mathrm{}`$ sub-regions adopted in producing Figure 3 of Tecza et al. (2000) and our Figure 1. A similar artefact is also seen in the map of CO. The adopted velocity range for integrating the profiles (from $`dV=300`$ km s-1 to $`dV=300`$ km s-1) may also affect the positions of the CO and H2 peaks in our maps. . The velocities of both lines change along PA$`45`$ deg, in the sense that they are blueshifted around the S nucleus and are redshifted in both the NE and SW regions. Simple models of the velocity field, such as a rotation or a pure expansion/contraction, have difficulty reproducing these observations.

To compare the properties of CO and H2 emissions, we here introduce a new index, the H2 emitting efficiency. The efficiency is defined as the intensity ratio of H2 and CO \[$`I`$(H2)/$`I`$(CO)\], and was calculated over the entire velocity range on each sub-region (Figure 2). Note, however, that these ratios are shown only at velocities where the normalized CO intensity is larger than 0.05 (the typical noise level) to avoid obtaining erroneous ratios in regions or velocities where CO is weak. Next, the integrated H2 emitting efficiencies for blue \[$`R_{\mathrm{blue}}I(\mathrm{H}_2)/I(\mathrm{CO})𝑑V/\mathrm{\Delta }V`$, integrating over the range $`dV=300`$ to $`0`$ km s-1\], for red ($`R_{\mathrm{red}}`$, integrating over the range $`dV=0`$ to $`+300`$ km s-1), and for blue plus red (total) ($`R_{\mathrm{total}}`$, integrating over the range $`dV=300`$ to $`+300`$ km s-1) were calculated, in order to examine the spatial and velocity distributions of the efficiencies. Here $`\mathrm{\Delta }V`$ is the velocity range for the integration, i.e., $`\mathrm{\Delta }V=300`$ km s-1 for blue- and red-, and $`\mathrm{\Delta }V=600`$ km s-1 for total-integrated H2 emitting efficiencies. A velocity cut-off (at $`dV=\pm 300`$ km s-1) was applied in integrating the H2 emitting efficiency, in the same way as for the calculation of the integrated intensities. Also, the integrated H2 emitting efficiencies in sub-regions where the normalized intensity of CO at the cut-off velocity was less than 0.05 were not used in the following analysis. The blue-to-red ratio of the H2 emitting efficiency ($`R_{\mathrm{blue}}`$/$`R_{\mathrm{red}}`$) could then be calculated. These values ($`R_{\mathrm{total}}`$ and $`R_{\mathrm{blue}}`$/$`R_{\mathrm{red}}`$) are mapped together with the velocity difference between CO and H2 \[$`\mathrm{\Delta }V_{\mathrm{H}_2\mathrm{CO}}V`$(H2)$`V`$(CO)\] in Figure 3.

Figure 3 illustrates our main results: (1) The total H2 emitting efficiency ($`R_{\mathrm{total}}`$) is smaller around the CO peak, and becomes higher with distance from the peak; (2) the blue-to-red ratio of the H2 emitting efficiency ($`R_{\mathrm{blue}}`$/$`R_{\mathrm{red}}`$) is smaller around the CO peak, and becomes larger with distance from the peak; and (3) the velocity difference ($`\mathrm{\Delta }V_{\mathrm{H}_2\mathrm{CO}}`$) is smaller around the CO peak, while the H2 emission becomes more blueshifted with respect to CO with distance from the peak. All these trends are seen almost systematically as a function of position, and are evident well over our resolution ($`0.6\mathrm{}`$ in space and 140 km s-1 in velocity). It seems noteworthy that the peak of H2 emitting efficiency corresponds neither to the peak position of CO nor to that of H2.

We searched for correlations between the measured quantities ($`R_{\mathrm{total}}`$ and $`R_{\mathrm{blue}}`$/$`R_{\mathrm{red}}`$); the results are shown in Figures 4 and 5. In these plots, data are marked with different symbols for different groups of sub-regions, which were divided according to their spatial positions around the CO intensity peak: the central region with higher H2 and CO intensities between the two nuclei (nuclear region), the region northeast of the central region (NE region), and the southwestern region (SW region). Two other sub-regions, the most southwesterly regions where H2 intensity increases with distance from the nucleus (the SW sub-group) and the most northeasterly regions with double-peaked profiles (the NE sub-group), are grouped separately. These sub-regions are also indicated in Figure 2. We found a tight negative correlation between $`I`$(CO) and $`R_{\mathrm{total}}`$ (Figure 4). A Spearman rank-order test for the correlation between log\[$`R_{\mathrm{total}}`$\] and log\[$`I`$(CO)\] showed that the probability $`P`$ of consistency with the null hypothesis (i.e., no correlation between the two quantities) is less than 0.001 for any sub-grouping (all the regions, all the regions except the SW sub-group, and all regions except for both SW and NE sub-groups). Therefore we consider that we have definitely shown that the total H2 emitting efficiency is enhanced in the outer regions with weaker CO emission. Two other positive correlations were found, between $`\mathrm{\Delta }V_{\mathrm{H}_2\mathrm{CO}}`$, $`R_{\mathrm{blue}}`$/$`R_{\mathrm{red}}`$, and $`R_{\mathrm{total}}`$ (Figure 5). Although these correlations look less tight, a Spearman rank order test showed that the probability $`P`$ of the null hypothesis is 0.20, 0.02, and 0.04 for the above three sub-groupings, respectively, in $`\mathrm{\Delta }V_{\mathrm{H}_2\mathrm{CO}}`$ vs. $`R_{\mathrm{blue}}`$/$`R_{\mathrm{red}}`$ diagram, and 0.02, 0.001, and $`<`$0.001 for the same three sub-groupings, respectively, in $`\mathrm{\Delta }V_{\mathrm{H}_2\mathrm{CO}}`$ vs. $`R_{\mathrm{total}}`$ diagram, indicating the presence of significant correlations between these quantities at greater than the 98% significance level if we exclude the SW sub-group.

## 3 DISCUSSION

### 3.1 H2 Excitation Mechanism

Our analysis has shown that the H2 emitting efficiency and the kinematical conditions of the warm gas are closely related. Such a correlation is a good indication of shock heating, in which bulk kinetic energy is converted to energy radiated through the H2 line. Although our analyses were carried out only for fluxes within the velocity range of $`dV=300300`$ km s-1, it seems reasonable to interpret the correlation as an indication of the major mechanism of the excitation of H2, because this velocity range contains 81% and 67% of the total fluxes of CO and H2, respectively (see Section 2 above). In some inner sub-regions where we can trace the H2 emitting efficiency at even bluer velocities ($`dV<300`$ km s-1), we found that the efficiency is generally higher at bluer velocity ($`dV<300`$ km s-1) than at lower velocity ($`dV=3000`$ km s-1), and this fact supports a model in which shock excitation of H2 is a major excitation mechanism around NGC 6240. Although other H2 excitation mechanisms (e.g., UV fluorescence, X-ray heating, and formation pumping) cannot be rejected, their contribution to the total must be less important. We note that the same conclusion has previously been reached based on a line-ratio analysis of the near-infrared spectra (Sugai et al. 1997; Tecza et al. 2000; and Ohyama et al. 2000).

### 3.2 The Superwind Model

Since correlations between the H2 emitting efficiency and the kinematical conditions of the molecular gas were found not only on the SW but also on the NE side of the CO peak, on scales of $`2\mathrm{}`$, the molecular gas concentration seems likely to be surrounded by gas that has a simple global motion. The motion is probably directed predominantly toward the observer, and the interaction between the gas that has this blueshifted motion and the molecular gas concentration would increase the H2 emitting efficiency on the bluer velocity side of the CO emission. The question then arises: what is the origin of this blueshifted gas motion? In the following, we discuss two possible explanations.

One possibility is that the merging nuclei are interacting with the off-nuclear molecular gas as they rotate around each other. Although there are no good theoretical explanations for the formation of such an off-nuclear molecular gas concentration in merging galaxies (e.g., Barnes 2002 and references therein; see also Tecza et al. 2002), here we assume that such concentration has formed as a result of tidal effects during the merging process. Since the H2 emitting efficiency is relatively higher at the SW portion of the concentration, a collision between the S galaxy and the concentration might contribute to the excitation of H2 emission. The extinction maps (Scoville et al. 2000; Tecza et al. 2000) have revealed the presence of a dusty region at 0″– 1.0″ NNE of the S nucleus, and the position of the maximum extinction ($``$ 0.2″– 0.3″ NNE of the S nucleus) coincides spatially with the CO peak. This indicates that the S galaxy is behind the dusty molecular gas concentration from our point of view. Since the S nucleus is known to be moving toward us around the concentration (Fried & Ulrich 1985; Tecza et al. 2000), a violent collision between the two would enhance the blueshifted H2 emission at the SW side of the CO peak. If this is the case, however, the redshifted H2 emission, rather than the blueshifted component, would be enhanced on the opposite side of the CO peak, because the N galaxy is rotating away from us on the opposite side of the concentration. This model therefore cannot explain the fact that blue-enhanced H2 emitting efficiency is observed on both sides of the CO peak.

The alternative model is that an outflow associated with the superwind activity of the S nucleus interacts with the off-nuclear molecular gas concentration. The superwind activity of this galaxy has been recognized as a result of its large-scale ($`20`$ kpc) shock-excited filamentary nebulae (e.g., Armus, Heckman, & Miley 1990; Keel et al. 1990) and the expanding motion of the H2 emitting clouds around the nuclear region (Ohyama et al. 2000). Since both of the individual galaxies are known to show flat rotation curves in their stellar velocity fields (Tecza et al. 2000), disk structures are likely to be present around each nucleus, even though the system is an advanced merger. Therefore we expect any superwind outflow to exhibit a bi-directional conical shape pointing toward the disk polar direction (e.g., Tomisaka & Ikeuchi 1988; Heckman et al. 1990; Suchkov et al. 1994; Strickland & Stevens 2000), although the shape of the outflow might be distorted due to the inhomogeneous distribution of the medium around the nucleus. An $`H`$band image of the S nucleus reveals an elongated structure in the NW-SE direction, suggesting the presence of a disk whose minor axis is directed NE-SW at an inclination angle of $`i60`$ deg (Scoville et al. 2000; Tecza et al. 2000). Also, as mentioned earlier, the S nucleus seems to be located behind the molecular gas concentration. Assuming that the off-nuclear molecular gas concentration formed during the merging process of the two nuclei, before the superwind started to blow, we would expect that, if a bi-conical superwind outflow emanated from the S nucleus, it would interact with the gas concentration from the rear. This interaction would produce blueshifted H2 emission around the CO peak (i.e., not only at the SW but also at the NE sides of the CO peak) if the opening angle of the outflow were wide enough to cover the whole molecular gas concentration. The oppositely directed cone of the superwind, pointing away from us, would not contribute significantly to the H2 luminosity, as it would not interact with the off-nuclear molecular gas. Since this model has no major difficulties in reproducing the observations, we regard it as more plausible than the model in which the S nucleus is colliding with the molecular gas concentration. A schematic picture of the superwind model is shown in Figure 6. Note that the real situation would be much more complicated due to the complicated medium distribution around the merging nuclei.

Following Heckman et al. (1993), clouds accelerated by the superwind outflow will have a terminal velocity of $`v_{\mathrm{cloud}}=400L_{\mathrm{bol},11}^{1/2}r_{\mathrm{kpc}}^{1/2}N_{21}^{1/2}`$ km s-1, where $`L_{\mathrm{bol},11}`$ is the bolometric luminosity in units of $`10^{11}`$ $`L_{\mathrm{}}`$, $`r_{\mathrm{kpc}}`$ is the initial distance of the cloud from the nucleus, and $`N_{21}`$ is the column density of the interstellar clouds in units of $`10^{21}`$ cm-2. Adopting $`L_{\mathrm{bol},11}=4.6`$ (Sanders & Mirabel 1995), $`r1.0`$ kpc (the projected distance between the S nucleus and the CO peak multiplied by a de-projection factor of $`2`$), and $`N_{21}=120`$333 Here we adopt $`N_{21}=1`$ as a minimum acceptable value corresponding to the typical value for the Milky Way. We also adopt $`N_{21}20`$ as a maximum acceptable value corresponding to one tenth of the total column density inferred from molecular gas observations of $`N_{21}200`$ (Tacconi et al. 1999).  , this equation implies that the clouds would be blown off at a speed of between a few and several hundreds of km s-1, which is almost comparable to the observed line width of H2 ($`500`$ km s-1 in FWHM: e.g., Ohyama et al. 2000; Tecza et al. 2000) and CO ($`260400`$ km s-1 in FWHM: Tacconi et al. 1999). Since the column density is likely to decline outward around the CO peak, falling off in a similar way to the CO surface brightness, (Figure 1a) 444The total column density of the molecular gas integrated over the line of sight is proportional to the observed integrated CO intensity, so the total CO intensity map is similar to the map of the total column density. Although the column density of each cloud may not be proportional to the total one, they are likely to be almost proportional to each other if we assume a simple homogeneous distribution of small clouds within the molecular gas concentration., i.e., the column density of each cloud is larger near the CO peak and is smaller at the outer region, so that clouds in the outer part of the molecular gas concentration would be accelerated up to faster speeds than those near to the CO peak. Therefore the observed trend in the relatively blueshifted H2 emission with respect to the CO emission in the outer part of the molecular gas concentration can be explained naturally by the model of superwind-molecular gas concentration interaction. This fact further suggests that kinematical disturbance due to the superwind also raises the H2 emitting efficiency within the molecular gas concentration, given the positive correlation between the H2 emitting efficiency and $`\mathrm{\Delta }V_{\mathrm{H}_2\mathrm{CO}}`$.

It has been pointed out that a C-type shock ($`v_{\mathrm{shock}}40`$ km s-1), rather than a J-type shock with faster shock velocity, is responsible for exciting most of the H2 emission, although most of the emission comes from clouds within a faster global velocity field ($`dV>100`$ km s-1) (Sugai et al. 1997). One way to resolve this apparent contradiction is to suggest that internal cloud-cloud collisions within the clumpy medium, at speeds of $`30`$$`50`$ km s-1, excite the H2 emission, while the clouds move along the global velocity field as a whole (Tecza et al. 2000). Since the molecular gas concentration is known to form a turbulent thick rotating disk (Tacconi et al. 1999), it is likely that the concentration is actually composed of numerous smaller clouds (Tecza et al. 2000). A high-resolution velocity dispersion map of the CO (Figure 4 of Tacconi et al. 1999) reveals the presence of regions with higher velocity dispersion just around the CO peak (especially to the SE and NW sides of the CO peak). Therefore it seems very likely that the interaction of small clouds within the molecular gas concentration with the external superwind provides an enhanced random velocity at the outer part of the CO concentration, as well as the blueshift of the H2 with respect to the CO.

### 3.3 The Cloud-Crushing Model

One of the merits of utilizing cloud-cloud collisions within the molecular gas concentration to give rise to the H2 emission is that a higher efficiency of energy deposition into the molecular gas clouds from the supernova explosions can be expected as a result of the cloud-crushing mechanism (Cowie, McKee, & Ostriker 1981), and the observed huge H2 luminosity can be explained (Draine & Woods 1990; Elston & Maloney 1990). This therefore raises the question of whether the observed trends in H2 emitting efficiency, as well as the total H2 luminosity, can be explained in terms of the cloud-crushing mechanism. In this model, we may consider, as a simple case, that some fraction of the kinetic energy of the small clouds, randomly moving within the bulk motion of the rotating molecular gas concentration, is dissipated through internal shocks as a result of cloud-cloud collisions, and gives rise to H2 emission. Then, the H2 surface brightness at a velocity of $`v`$ can be expressed as $`I_{\mathrm{H}_2}(v)=0.5\mathrm{\Sigma }_{\mathrm{H}_2}(v)\zeta (v)<v_{\mathrm{internal}}^2>ϵ/\mathrm{\Delta }t`$, where $`\mathrm{\Sigma }_{\mathrm{H}_2}(v)`$ is the surface mass density of the molecular gas concentration whose velocity is $`v`$, $`\zeta (v)`$ is the mass fraction of the small clouds that give rise to H2 emission through the cloud-crushing mechanism at a velocity $`v`$, $`<v_{\mathrm{internal}}^2>`$ is the mean of the squared internal random velocity of the small clouds, $`ϵ`$ is the fraction of the energy dissipated in the shock which is radiated in the H2 1-0 $`S`$(1) line, and $`\mathrm{\Delta }t`$ is the timescale of the cloud-cloud collisions. Since the surface mass density of the clouds $`\mathrm{\Sigma }_{\mathrm{H}_2}(v)`$ is proportional to the observed CO surface brightness \[$`I_{\mathrm{CO}}(v)`$\] at the velocity of $`v`$, the H2 emitting efficiency can be expressed as $`I_{\mathrm{H}_2}(v)/I_{\mathrm{CO}}(v)\zeta (v)<v_{\mathrm{internal}}^2>ϵ/\mathrm{\Delta }t`$, by definition. Since we are considering a situation in which the cloud-crushing mechanism is responsible for all the H2 emission, we may substitute a constant $`ϵ`$ ($`0.3`$) and a constant internal velocity ($`<v_{\mathrm{internal}}^2>^{1/2}=4050`$ km s-1) into the equation, representing an extreme case of the cloud-crushing mechanism working at its maximum efficiency (Draine & Woods 1990). Then, the H2 emitting efficiency can be expressed as $`I_{\mathrm{H}_2}(v)/I_{\mathrm{CO}}(v)\zeta (v)/\mathrm{\Delta }t`$. This expression indicates that the H2 emitting efficiency is higher when (1) a larger mass fraction of the clouds is involved in the production of H2 emission through the cloud-crushing mechanism, and/or (2) the collision frequency is higher (i.e., the collision timescale is shorter). The observed close relationship between the H2 emitting efficiency and the kinematical conditions in the molecular gas clouds can be explained by this expression, if the cloud-crushing mechanism is efficient: The outer part of the off-nuclear molecular-gas concentration is likely to be more strongly influenced kinematically by the interaction with the external superwind (see Section 3.2 above), and the H2 emitting efficiency is expected to be higher, since a larger mass fraction of the clouds will produce H2 as a result of internal cloud-cloud shocks and a higher frequency of the cloud-cloud collisions. In addition, the efficiency is expected to be higher at bluer velocities, because these are likely to be caused as a result of a more violent kinematical influence from the superwind. Therefore, the cloud-crushing mechanism within a molecular gas concentration disturbed by the interaction with the external superwind outflow seems to be able to reproduce the observed trends of the H2 emitting efficiency in both space and velocity, as well as the higher-than-normal H2 luminosity. Although the mechanism by which the off-nuclear molecular gas concentration between the two merging nuclei is created is not clear, the presence of such a gas concentration is necessary if the model is to reproduce the observed intense and efficient H2 emission. We look forward to future theoretical studies of the origin of such a gas concentration.

We thank Matthias Tecza who kindly allowed us to see his doctoral thesis and gave us useful comments on the “Data and Analysis” section of this paper.


text:
# Universal Description of Spherical Aberration Free Lenses Composed of Positive or Negative Index Media

## Abstract

We find that the function that describes the surface of spherical aberration free lenses can be used for both positive and negative refractive index media. With the inclusion of negative index, this function assumes the form of all the conic sections and expands the theory of aplanatic optical surfaces. There are two different symmetry centers with respect to the index that create an asymmetric relationship between positive and negative index lens profiles. In the thin lens limit the familiar formulas for image position and magnification hold for any index.

It was known to René Descartes as early as 1637 that aspherical surfaces are required to convert plane waves into spherical waves and vice versa. In the short wavelength limit of geometric optics, the right conic surface can bring all parallel incident light rays to a single common focus point. Most commercial lenses have spherical profiles and only perform this function approximately; they are said to have *spherical aberration*. It is only somewhat well known that a hyperboloid is the ideal lens shape for this purpose when using media with a relative refractive index greater than one. Less well known, is the fact that an ellipsoid is the ideal shape for a lens composed of media with relative index between zero and one hecht .

Much of the fundamental behavior of media with negative refractive index was predicted in 1968 by Veselago veselago . Such media are not naturally occurring, so their properties were not extensively pursued at that time. Recently, negative index media were demonstrated using engineered composites smithPRL ; shelby , and the same design paradigm used for those composites also enables the construction of media with index less than one and even somewhat close to zero. This makes the question of optical element design with these media pertinent. In this paper we will discuss ideal (i.e. spherical aberration free) shapes for lenses composed of media with both positive and negative index. We will derive universal formulas for lens design that apply to media of any index, and discuss the symmetry between positive and negative index lenses. We will also demonstrate that, with the inclusion of negative index, a more compete theory of aplanatic optical surfaces is obtained.

The lenses discussed here need to be clearly distinguished from the “perfect” lens of Pendry perfectLens . These lenses can focus incident plane waves, i.e. waves from a source object at infinite distance. They are said to possess an aplanatic point at infinity bornWolf . Plane waves are *not* focused by the “perfect” lens; it has aplanatic points only at distances within one lens thickness. The importance of the “perfect” lens lies in its near field focusing capability. This work describes lenses that can focus or collect energy from distant sources where the lenses are composed of materials from an extended parameter space, namely negative refractive index.

Aspherical lenses are uncommon because they are difficult to manufacture. The requirement of $`\lambda /4`$ to $`\lambda /10`$ surface roughness for good optics is quite difficult to achieve in the visible range for non spherical surfaces. However, this requirement is automatically satisfied for surfaces on the composite materials that have been used to implement negative index media. These media are composed of unit cells that must be significantly smaller than the operational wavelength in order for them to function as a homogenous media xlowFreqPlasmons ; smithPRL . The manufacturing process that enables construction on the required sub unit cell length scales can by default control the surface profile with one unit cell accuracy. Though the current technology may never extend to visible light, whatever the operational wavelength of these media, sub-wavelength surface accuracy for any desired shape is available.

The results discussed here were first obtained by ray tracing hecht . The ray tracing algorithm employed only fundamental assumptions about the media dispersion , boundary matching, and the conservation of energy. For isotropic media, our algorithm is equivalent to Snell’s law, which applies to interfaces with relative index of either sign shelby ; parazzoliSnell ; pachecoDispersion ; FoteinopoulouFDTD . The conservation of energy was ensured using the Poynting vector, with each ray carrying power into a surface being required to carry power out.

Here, however, we show an analytical argument using optical path length (OPL) following Hecht hecht . We wish to find a surface that refracts parallel rays to a single point (Fig. 1). Parallel rays have planar wave fronts on which the phase is constant. Rays converging to a point have spherical wave fronts on which the phase is constant. The desired lens surface must be contoured such that, the optical path length in the transition region is the same for all rays. The optical path length or phase advance (retardation) is hecht ; bornWolf

$$OPL=_Cn𝑑s$$
(1)
The equivalence of the *OPL* in the transition region for the two rays shown in Fig.1 yields,

$$\sqrt{\left(fz\right)^2+\rho ^2}f=nz$$
(2)
From this we find

$$\frac{\left(za\right)^2}{a^2}+\frac{\rho ^2}{b^2}=1$$
(3)
which is an equation for conic sections where the $`a,b`$ parameters are given in terms of the index, $`n`$, and focal distance, $`f`$.

$`a`$ $`=f{\displaystyle \frac{1}{1+n}}`$ (4a)
$`b^2`$ $`=f^2{\displaystyle \frac{1n}{1+n}}`$ (4b)
Since $`b^2`$ can be negative, Eq. (3) can and does generate all the conic sections for different values of the index. Further, nothing in this derivation precludes $`n`$ from being negative, as long as the optical path length as given by Eq. (1) is applicable. We can see this as follows. Geometric rays follow the energy direction indicated by the Poynting vector. It is well established that, in isotropic negative index media, the phase advance direction is opposite to the Poynting vector veselago ; modulate ; pachecoDispersion ; FoteinopoulouFDTD , so it is appropriate to use negative index in Eq. (1).
The focal length, $`f`$, can also be negative to find lens surfaces that generate diverging spherical waves from plane waves (i.e. a virtual focus). The following functional form gives a surface centered at $`z=r=0`$ with the correct concavity for plane waves on the negative $`z`$ side and spherical waves on the positive $`z`$ side.

$$z(\rho )=a\left(1\sqrt{1\frac{\rho ^2}{b^2}}\right)$$
(5)
This equation is valid for all values of $`f`$ and $`n`$, except $`n=1`$, where it is singular. This limiting case is handled by letting $`n=1+\delta `$ with $`\delta 1`$, and substituting into Eq. (5) with Eq. (4). We then find

$$z(\rho )\frac{1}{4f}\rho ^2$$
(6)
which is a paraboloid.

From Eqs. (4),(5) and (6), we find the following behavior for the ideal lens surface. These results are also confirmed by ray tracing, Fig. 2.

| index | conic | converge/diverge | focus |
| --- | --- | --- | --- |
| $`1<n`$ | hyperboloid | convex/concave | far |
| $`n=1`$ | no solution | - | - |
| $`0<n<1`$ | ellipsoid | concave/convex | far |
| $`n=0`$ | sphere | concave/convex | far = near |
| $`1<n<0`$ | ellipsoid | concave/convex | near |
| $`n=1`$ | paraboloid | concave/convex | near |
| $`n<1`$ | hyperboloid | concave/convex | near |

The conic sections characterized by the $`a,b`$ parameters above, can be equally well characterized by another set of parameters, the eccentricity, $`\epsilon `$ and the major axis radius of curvature, $`R`$, (Fig. 3 inset). This pair of parameters is useful in that it separates the shape, given by $`\epsilon ,`$ from the length scale, given by $`R`$. Using the usual definition, eccentricities less than one indicate ellipses and greater than one indicate hyperbolas. Applying Eq. (4) leads to the simplification

$$\epsilon \frac{\sqrt{a^2b^2}}{\left|a\right|}=\left|n\right|$$
(7)
Thus the eccentricity is symmetric in $`n`$, (Fig. 3). If this were the only parameter, we could take any positive index lens design and use it for lenses composed of negative index media.

The radius of curvature can be calculated from Eq. (5).

$$R\frac{1}{z^{\prime \prime }\left(0\right)}=\frac{b^2}{a}=f(1n)$$
(8)
This is recognized as the familiar lens makers formula usually derived using spherical optics with the paraxial approximation. This is not surprising since any of the conic sections approximate a sphere when examined with small enough aperture; the two derivations must agree on the lens central radius of curvature. Ray tracing confirms this formula applies for both positive and negative index and correctly gives a change in the concavity at $`n=1`$. Because this formula is antisymmetric about $`n=1`$, (Fig. 3), we cannot use the same lens design for positive and negative index. If we do, the positive and negative index lenses will have different focal lengths. In fact, for $`\left|n\right|>1`$, one lens will be converging and the other diverging. Notably, the positive index analog of the $`n=1`$ parabolic lens is absent, since the $`n=1`$ lens requires zero radius of curvature.

We have not yet shown how, if at all, the position of the optical focus of one of these lenses relates to the geometric foci of its conic surface. It can be easily shown algebraically or by ray tracing, that the optical focus lies exactly at the position of one of the conic section’s foci. Ray tracing also confirms a simple rule that specifies which focus. For positive index the correct focus is the one furthest from the branch of the conic that is used for the lens surface. For negative index, the opposite is true. For the lenses with parabolic ($`n=1`$) and spherical ($`n=0`$) surfaces that posses just one geometric focus, it is also the optical focus.

As pointed out in Hecht hecht , two conic surfaces can be combined to construct lenses that convert diverging spherical waves to converging spherical waves. Ray tracing confirms that in the thin lens limit, the usual formulas apply for image distance

$$\frac{1}{f}=\frac{1}{o}+\frac{1}{i}$$
(9)
and magnification

$$M=\frac{i}{o}=\frac{f}{fo}$$
(10)
regardless of the sign of the index. In Fig. 4, a double concave parabolic lens is shown to function as a converging lens when composed of $`n=1`$ material.

The extension of spherical aberration free lens design to negative index media yields a more complete theory, but the practical limitations of current negative index media, notably strong chromatic aberration (due to media dispersion) may preclude technological interest in frequency ranges where good positive index materials are available. However, in frequency ranges where this is not the case, (e.g. the millimeter wave range), and composite media are already in use, negative index lenses can be advantageous. Where larger radius of curvature is desirable, it is worth noting that an $`n=2`$ lens has three times the radius of curvature of an $`n=2`$ lens of the same focal length, (Fig. 2). Even more significant is the fact that an $`n=1`$ lens is impedance matched to the surrounding media perfectLens ; *none* of the incident radiation is reflected. This is impossible with a positive index lens.

###### Acknowledgements.

This work was supported by DARPA Contract No. MDA972-01-2-0016 and by DARPA through a grant from ONR, Contract No. N00014-00-1-0632.