text:
# Investigation of the possible third image and mass models of the gravitational lens PMN J1632–0033

## 1 Introduction

The central density profile $`\rho (r)`$ of galaxies is poorly known, but is important for understanding galaxy structure and the nature of dark matter. Theoretical cold-dark-matter (CDM) models predict a steep central density cusp, $`\rho r^\beta `$ with $`\beta `$ 1–1.5 (see, e.g., Navarro, Frenk, & White, 1997; Moore, 1999). However, measurements of spiral galaxy rotation curves suggest the dark matter has a constant-density core. This has been argued most strongly for low-surface-brightness galaxies (de Blok et al., 2001) but may hold for all spiral galaxies (Debattista & Sellwood, 2000; Salucci & Burkert, 2000), including our own Galaxy (Binney & Evans, 2001). The conflict may be due to observational effects such as beam smearing, at least in part (van den Bosch & Swaters, 2001), but the apparent conflict has encouraged some theorists to invent mechanisms that avoid producing cusps (e.g., bar-induced evolution, Weinberg & Katz 2002; warm dark matter, Bode, Ostriker, & Turok 2001; self-interacting dark matter, Spergel & Steinhardt 2000).

For elliptical galaxies, where rotation curves cannot be measured, the central density profile is even less well understood. Observations of early-type galaxies with the Hubble Space Telescope (HST) show that the light distribution is cuspy, with cusp parameters that may be correlated with global galaxy properties (Lauer et al., 1995; Faber et al., 1997). But such studies do not probe mass directly, and are limited to nearby galaxies due to the requirements on angular resolution.

Multiple-image gravitational lenses can be used to measure the central surface density of galaxies at significant redshift, by searching for faint central images of the background object. In most cases, lenses should produce an odd number of images (Dyer & Roeder, 1980; Burke, 1981), one of which is faint and appears near the center of the lens galaxy. This “central image” or “odd image” corresponds to the central maximum in the time-delay surface. However, of the $``$80 known systems in which a galaxy produces multiple images of a quasar or radio source, almost all consist of either 2 or 4 images. One system, APM 08279+5255, definitely has three quasar images (Lewis et al., 2002b), but the third image may be due to a “naked cusp” configuration rather than a central time-delay maximum (Lewis et al., 2002a).

No central image has been identified definitively, even among the radio lenses, where the central image could in principle be seen through the light and dust of the lens galaxy. The absence of central images has been attributed to the demagnification of the image by the high central surface density of the lens galaxy (see, e.g., Narasimha, Subramanian, & Chitre, 1986; Wallington & Narayan, 1993; Rusin & Ma, 2001; Evans & Hunter, 2002; Keeton, 2002). An “extra” radio component has been detected in several radio lenses, and in three of those cases it is still unclear whether the extra component is an additional quasar image or faint radio emission from an active galactic nucleus (AGN). Those three systems are 0957+561 (Harvanek et al., 1997), MG J1131+0456 (Chen & Hewitt, 1993), and the subject of this paper: PMN J1632–0033 (Winn et al., 2002).

Radio components A and B of J1632–0033 are two images of a $`z=3.42`$ quasar, but the third, faint component C is of unknown origin (see Fig 1). In an HST image, the lens galaxy appears to be a fairly circular early-type galaxy with an optical effective radius of $`0\stackrel{}{\mathrm{.}}2`$, but it is too faint to be characterized in detail. The position of component C is near the optical lens galaxy position, as expected for either an AGN or a central odd image. If A, B, and C are all images of the same quasar, then they should have similar radio spectra. However, Winn et al. (2002) only presented measurements of the flux density of C at one frequency.

Here we present new VLBA data that, in combination with the previous data, allow us to measure the flux densities of all three components at four widely spaced frequencies. Our analysis shows that C has a significantly different radio spectrum than A and B. We argue that there are two alternative explanations for this difference. Either C is not a third quasar image, and instead is the active nucleus of the lens galaxy; or, C is a third image, but its spectral index is inverted by free-free absorption in the lens galaxy.

Whichever of these possibilities is true, the additional information from the VLBA data provides more constraints on mass models of the lens galaxy than are usually available for a two-image lens. If the central component is a third image, then its position and flux are valuable constraints. If the central component is the lens galaxy, then our data provide an unusually precise position for the lens galaxy, and a strong upper limit on the flux of a third image. The VLBA maps also reveal the relative orientations of the radio jets associated with A and B. As with all gravitational lensing studies, we still do not have enough information to determine the complete density profile of the lens galaxy, but we can use parameterized models to explore the implications of the data.

Historically, two classes of parameterized models have been used to study the central-image problem: models with constant-density cores, and models with central cusps. Early studies used softened power laws, $`\rho (r^2+r_c^2)^{\beta /2}`$, and found that central images are suppressed if the core radius $`r_c`$ is sufficiently small (see, e.g., Wallington & Narayan, 1993). Since the demagnification of the central image is determined by the central surface density, steeper density profiles (larger $`\beta `$) allow for larger core radii for the same flux limit on the central image. Power laws shallower than isothermal ($`\beta <2`$) produce central images even in the limit $`r_c0`$.

With the gradual realization that theoretical dark-matter profiles and observed surface-brightness profiles of early-type galaxies have cusps rather than cores, attention has shifted to the properties of density distributions with central cusps. Rusin & Ma (2001) used the absence of odd images in the class sample, presently the largest homogeneous sample of radio lenses, to argue that the global mass profiles of the lens galaxies cannot be much shallower than isothermal. Muñoz, Kochanek, & Keeton (2001) reached a similar conclusion about the particular radio lens class B1933+503. Keeton (2002) found that the distribution of stars observed in HST images of nearby early-type galaxies is often sufficiently concentrated to suppress the central image without recourse to dark matter.

We describe our observations of J1632–0033 in § 2, and our interpretation of the results in § 3. In § 4, we compute lens models for this system under the two different hypotheses for the nature of the central component. We consider density distributions that are scale-free power laws, and also determine the required properties of a constant-density core or an inner cusp. Finally, in § 5, we summarize our conclusions and discuss future observations that could test more definitively whether PMN J1632–0033 is a three-image quasar.

## 2 Observations

We observed J1632–0033 with the VLBA on 2002 March 14 and 15, for one eight-hour session each day. On the first day, the array included one VLA antenna in place of the VLBA antenna at Pie Town, but on the second day, the array consisted of the usual ten antennas. During both sessions we alternated between observations at the standard 8.4 GHz (3.6 cm) band and the standard 1.7 GHz (18 cm) band. For both bands, the observing bandwidth of 32 MHz per polarization was divided into 4 sub-bands. Both senses of polarization were recorded with 2-bit sampling. The data were correlated in Socorro, New Mexico, producing 16 channels of width 500 kHz from each sub-band, with an integration time of one second.

Calibration was performed with aips555The Astronomical Image Processing System (aips) is developed and distributed by the National Radio Astronomy Observatory. using standard procedures. We used the observations of J1632–0033 to solve for residual delays, rates, and phases directly (rather than using phase-referencing), with a fringe-fitting solution interval of 2 minutes. We used the multiple-field clean deconvolution algorithm available in aips to deconvolve a field centered on A and, simultaneously, a field of the same size centered on the mid-point between B and C. For the 8.4 GHz data, the fields were $`1024\times 1024`$ with a scale of 0.2 mas pixel-1. For the 1.7 GHz data, the fields were $`512\times 512`$ with a scale of 1.0 mas pixel-1.

At first, the data from each day were analyzed separately. Once we verified that the results of the two sessions were consistent, we combined the visibility data from each frequency to produce final maps. The data were self-calibrated, based on the clean model derived from the preliminary maps, with a 15-second solution interval. The final 1.7 GHz maps are shown in Fig. 2, and the final 8.4 GHz maps are shown in Fig. 3. In the vicinity of components B and C, the noise level in the maps is close to the theoretical limit. However, the noise level is three times larger in the vicinity of component A due to large sidelobes. Further improvement is difficult because of the near-equatorial location of the object, which results in poor coverage of the visibility plane even for an eight-hour Earth-rotation synthesis observation.

The flux densities of A, B, and C at each frequency are reported in Table 1. Also included in this table are results based on the previous data of Winn et al. (2002), for those cases in which all three components were detected (or could have been detected, given the angular resolution). These previous data were re-analyzed using natural weighting rather than uniform weighting, to maximize the sensitivity to the faint component C and to match the present analysis. The results were consistent with previous results in all cases. Of particular interest was the improved 15 GHz map, in which component C was detected with flux density $`0.84\pm 0.19`$ mJy. This is consistent with the upper limit of 1.5 mJy given by Winn et al. (2002), but the positive detection of C in the lower-noise map provides another valuable point of comparison.

For the VLA and merlin666The Multi-Element Radio Linked Interferometry Network (merlin) is a U.K. national facility operated by the University of Manchester at Jodrell Bank Observatory on behalf of PPARC. data, the flux densities were measured by fitting a model consisting of 3 point sources to the visibility data, using difmap (Shepherd, Pearson, & Taylor, 1994). The VLBA data sets were too large for model-fitting in visibility space. For these data, the flux densities were measured in image space with aips, using an aperture surrounding each component. The quoted uncertainties are the larger of the rms noise and the spread in values obtained using different choices for the aperture.

With these measurements of the flux densities of all three components at four widely spaced frequencies, we have the opportunity to test whether C has the same radio continuum spectrum as A and B. Assuming that the flux density of each component obeys $`S_i(\nu )\nu _i^\alpha `$ (where $`\alpha _i`$ is the spectral index of component $`i`$), a linear least-squares fit to the data in Table 1 gives $`\alpha _\mathrm{A}=0.27\pm 0.07`$, $`\alpha _\mathrm{B}=0.25\pm 0.07`$, and $`\alpha _\mathrm{C}=+0.29\pm 0.18`$. The uncertainties in the spectral indices include a contribution from the uncertainties in the absolute flux density scale, which is at least 5% at all frequencies.

The comparison between $`\alpha _\mathrm{A}`$ and $`\alpha _\mathrm{C}`$ can be made more precise by examining the flux density ratios, because the ratios are independent of the uncertainties in the absolute flux density scale. The difference in spectral indices, $`\alpha _\mathrm{C}\alpha _\mathrm{A}`$, is the logarithmic slope of $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$. The flux density ratios $`S_\mathrm{B}(\nu )/S_\mathrm{A}(\nu )`$ and $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$ are plotted in Fig. 4. The solid line represents the best-fit power law to $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$, which has a slope $`\alpha _\mathrm{C}\alpha _\mathrm{A}=+0.51\pm 0.17`$, implying the spectral indices of C and A are discrepant at the $`3\sigma `$ level. The dashed line is a fit to the data under the hypothesis that C is a third quasar image whose spectrum had been modified by free-free absorption (see § 3). Also plotted are upper limits on $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$ from VLA data at the two highest frequencies.

The position angles of the jets in A and B were measured by fitting elliptical components to the VLBA maps with aips, using a procedure (JMFIT) that corrects for the ellipticity of the beam. For modeling purposes, we adopted the average of the 1.7 GHz and 8.4 GHz results, giving $`\varphi _\mathrm{A}=8\mathrm{°}`$ and $`\varphi _\mathrm{B}=69\mathrm{°}`$, measured east of north. We estimate the uncertainty in each position angle to be $`\pm 15\mathrm{°}`$, an error range that is approximately twice the difference between the 1.7 GHz and 8.4 GHz results. These ranges of position angle are illustrated by dotted lines in Figs. 2 and 3, although it is difficult to account for the elliptical beam in a visual inspection.

## 3 The nature of component C

The flux density ratio $`S_\mathrm{B}(\nu )/S_\mathrm{A}(\nu )`$ is nearly independent of observing frequency, as should be the case for gravitationally lensed images of a single compact source. By contrast, the flux density ratio $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$ increases significantly with frequency. This immediately suggests that C is not a third image of the same background source. Given that its position is consistent with the optical position of the lens galaxy, it would be natural to conclude that C represents faint radio emission from the lens galaxy, as discussed further in § 3.1. However, there are several mechanisms that might cause a third quasar image to appear to have a different spectral index than the two brighter images. These are discussed in § 3.2–§ 3.5. The only mechanism we cannot rule out is free-free absorption (§ 3.5).

### 3.1 Lens galaxy emission

How does the flux density of C compare with the flux density expected from a typical galaxy at the lens redshift? Given the radio flux density and spectral index of component C that were presented in § 2, and the lens redshift estimate777The lens redshift $`z_l=1.0\pm 0.1`$ was estimated by Winn et al. (2002) by requiring the photometric properties of the lens galaxy to conform to the fundamental plane of elliptical galaxies, a method invented by Kochanek et al. (2000). of $`z1`$, the implied 1.4 GHz luminosity of the lens galaxy is $`1.7\times 10^{24}h^2`$ W Hz-1.

This is significantly more powerful than the radio emission observed from ordinary galaxies, which is thought to arise from star formation and is limited to $`\stackrel{<}{}10^{23}`$ W Hz-1. However, it is within the observed range of luminosities from galaxies harboring active galactic nuclei (AGN). For example, Sadler et al. (2002) identified a sample of 420 AGN in the 2dF Galaxy Redshift Survey, of which 24% had a 1.4 GHz radio luminosity at least as large as the implied luminosity of component C (after correcting for their different choice of cosmology).

We conclude that the identification of component C as the active galactic nucleus of the lens galaxy is a reasonable one. This would make J1632–0033 unusual among the known lenses, but not unique; the lens galaxy of B2045+265 has a bright central component that is almost certainly an active nucleus (Fassnacht et al., 1999).

### 3.2 Variability

If the source is variable, then the instantaneous flux density ratios between the lensed images will fluctuate, due to the differential time delays between the images. If, in addition, the degree of variability depends on frequency, then the flux density ratios will vary with frequency. However, the differential time delay between components B and C is expected to be very small ($`<`$0.5 day). This is much smaller than typical time scales for quasar variability ($``$months). Indeed, our measurements, which were spaced apart by one day, showed no sign of variability. We therefore discount the possibility that variability is the reason for the discrepancy of spectral indices.

### 3.3 Radio substructure

If the background radio source is extended and has sub-components with different spectral indices, then these sub-components may each be magnified by a different factor (see, e.g., Patnaik & Porcas, 1999). In that case, the sum of flux densities of the sub-components in a lensed image would not necessarily have the same spectral index as the corresponding sum in a different image.

However, in the case of J1632–0033, the source is very compact. From Fig. 3 we can limit the separation of sub-components in image A to $`<`$5 mas, corresponding to a change in magnification of only $`\mathrm{\Delta }\mu \stackrel{<}{}0.01`$ in an isothermal model. Furthermore, the agreement of $`\alpha _\mathrm{A}`$ and $`\alpha _\mathrm{B}`$ argues against the presence of a significant effect from frequency-dependent radio substructure.

### 3.4 Resolved-out radio structures

If a radio source has structure on large angular scales, with Fourier components corresponding to smaller spatial frequencies than those measured by an interferometer, then that portion of the radio source is invisible to the interferometer. The invisible radio structure is said to be “resolved out.” We must consider the possibility that the spectral indices of A and C are different only because the components are being resolved out by different amounts at each frequency.

We reject this hypothesis based on three facts. First, there cannot be much resolved-out structure in components A and B. During a monitoring program, we measured their flux densities at 8.4 GHz with the VLA on 2002 March 14, simultaneous with the VLBA measurement described in § 2. The VLBA flux densities were $`83\pm 7\%`$ of the VLA flux densities, despite being measured on angular scales $``$200 times smaller. Second, if C is a third quasar image, it is demagnified by a factor of $``$102 and would have an even smaller fraction of resolved-out structure than A or B. Third, even if there is resolved-out structure in A, it would be expected to have a steep radio spectrum, because the extended jets of quasars generally have steeper spectra than the compact cores. This would only worsen the discrepancy between the spectral indices of A and C, by causing the true spectral index of A to be steeper than observed.

### 3.5 Propagation effects

Propagation effects that are chromatic, and that affect each lensed image by a different amount, will cause the spectra of the images to differ. For example, reddening by dust causes the optical colors of lensed images to differ (see, e.g., Motta et al., 2002; Falco et al., 1999). At radio wavelengths, the propagation effects are due to interstellar plasma. The polarizations of lensed images may differ due to Faraday rotation (see, e.g., Patnaik et al., 1993). Or, in a phenomenon known as “scatter-broadening,” the angular size of a lensed image may be increased due to scattering by electron-density fluctuations along the line of sight (see, e.g. Jones et al., 1996). However, these two effects do not alter the spectral indices of lensed images. A possible exception might arise from scatter-broadening, if it causes a frequency-dependent amount of flux to be resolved out (see § 3.4 above), but our data do not show the characteristic $`\nu ^2`$ dependence of angular size due to scatter-broadening.

One propagation effect that would alter the spectral indices of lensed images is free-free absorption by electrons in the lens galaxy. Free-free absorption decreases strongly with frequency, which would flatten or invert the spectral index of the central component. We must therefore consider the possibility that component C is a third quasar image whose spectral index has been inverted by free-free absorption.

Assuming that C is affected by free-free absorption, and A is not, we would expect the flux density ratio to obey $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$ $`=\mu _{\mathrm{CA}}\mathrm{exp}(\tau _\nu )`$, where $`\mu _{\mathrm{CA}}`$ is the lensing magnification ratio and $`\tau _\nu `$ is the optical depth to free-free absorption. Here $`\nu =(1+z)\nu _{\mathrm{obs}}`$ is the frequency at the lens redshift. A useful approximation for $`\tau _\nu `$ has been given by Altenhoff et al. (1960) and Mezger & Henderson (1967):

$$\tau _\nu =0.08235\times \left(\frac{T_e}{{}_{}{}^{}K}\right)^{1.35}\left(\frac{\nu }{\mathrm{GHz}}\right)^{2.1}\left(\frac{N_e^2𝑑s}{\mathrm{pc}\mathrm{cm}^6}\right)=\left(\frac{\nu }{\nu _c}\right)^{2.1},$$
(1)
where $`T_e`$ and $`N_e`$ are the electron temperature and number density, and $`\nu _c`$ is defined implicitly. This formula predicts a sharp cut-off at low frequencies, and the asymptotic limit $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$ $`=\mu _{\mathrm{CA}}`$ as $`\nu \mathrm{}`$.

We solved for the values of the two parameters $`\nu _c`$ and $`\mu _{\mathrm{CA}}`$ that provide the best fit to the data. The results were $`\nu _c=3.2`$ GHz ($`\nu _{c,\mathrm{obs}}=1.6`$ GHz) and $`\mu _{\mathrm{CA}}=0.0045`$. The dotted line in Fig 4 shows the best-fitting curve, which is statistically consistent with the data ($`\chi ^2=1.1`$ with 3 degrees of freedom). There are no published cases of free-free absorption by a lens galaxy for comparison. However, the best-fitting value $`\nu _c=3.2`$ GHz is comparable to the values of $`\nu _c=`$ 1–10 GHz that have been inferred for the central regions of some active galaxies, using multi-frequency VLBI measurements of radio jets (for recent examples, see Marr, Taylor, & Crawford 2001, Jones et al. 2001, Kameno et al. 2000, Taylor 1996, and Levinson, Laor, & Vermuelen 1995).

Because the data are quantitatively consistent with both the third-image hypothesis and the free-free absorption hypothesis, we have only our prior prejudice on the relative likelihood of these two possibilities to guide us. One might object that the third-image hypothesis is contrived, because neither a central image nor free-free absorption has previously been observed in a radio lens. However, we do not consider this “double novelty” as a serious objection because a central image might be expected to suffer more absorption than other images, due to the potentially large electron column density near the center of the lens galaxy. For this reason, we take a neutral position, and consider the consequences of both hypotheses in the section that follows.

## 4 Updated lens models

In this section we use information from the VLBA maps to update the lens models for J1632–0033 that were presented by Winn et al. (2002). The goal of our analysis is to investigate the improvement that results when two extra pieces of information are available besides the usual image positions and fluxes: the relative orientations of the radio jets, and a strong constraint (either a measurement or an upper limit) on the flux of a central image. In § 4.1, we assume that component C is due to the lens galaxy. In § 4.2, we explore the consequences if component C is a third quasar image.

### 4.1 Case 1: Component C is the lens galaxy

In this case, we use the model constraints given in Table 2. The position of C is taken to be the center of the lens galaxy, which is known more precisely than the lens galaxy position in most other systems. The upper limit on the flux density of a third quasar image, expressed as the magnification ratio $`\mu _3/\mu _\mathrm{A}`$, represents the $`5\sigma `$ limit in the 1.7 GHz VLBA map. The magnification ratio $`\mu _{\mathrm{BA}}\mu _\mathrm{B}/\mu _\mathrm{A}`$ is the average of the 15 measurements of $`S_\mathrm{B}(\nu )/S_\mathrm{A}(\nu )`$ presented in this paper and by Winn et al. (2002), and the uncertainty is the variance in those measurements. In theory, our error estimates for the positions and magnification ratio should be enlarged (to $``$1-2 mas and 10-20%, respectively) to account for possible mass substructure in the lens galaxy (see, e.g., Metcalf & Madau, 2001; Dalal & Kochanek, 2002). However, in all the cases described below, the error statistic is dominated by either the third-image flux or the jet orientations, causing the precise values of the uncertainties on the positions and magnification ratio to be unimportant.

We define the goodness-of-fit parameter as

$$\chi ^2=\underset{i=\mathrm{A},\mathrm{B}}{}\left[\frac{(𝐑_i^{}𝐑_i)^2}{(\mathrm{\Delta }𝐑_i)^2}+\frac{(\varphi _i^{}\varphi _i)^2}{(\mathrm{\Delta }\varphi _i)^2}\right]+\frac{(\mu _{\mathrm{BA}}^{}\mu _{\mathrm{BA}})^2}{(\mathrm{\Delta }\mu _{\mathrm{BA}})^2},$$
(2)
in which primed quantities are those determined by the model, and unprimed quantities are those observed. In addition, we reject models that predict a third image brighter than the upper limit. This is computationally faster than the more orthodox procedure of adding another error term to Eq. 2 for the third-image flux, but we have verified that the results are very nearly equivalent. This is because the flux of the third image depends sensitively on the parameters in our models, causing $`\chi ^2`$ to increase very rapidly across the boundary in parameter space where the third-image constraint is violated.

Although in principle there could be a complicated angular structure in the lens model, due to the ellipticity of the lens galaxy and due to external perturbations, we do not have enough constraints to explore such models. Instead, we use mass models with circular symmetry, to which we add an external shear field (of strength $`\gamma `$ and position angle $`\theta _\gamma `$) to simulate the combined effect of galaxy ellipticity and tidal fields from neighboring mass concentrations.

We begin with global power-law models, for which $`\rho (r)r^\beta `$. The model has 7 free parameters: $`\beta `$, $`\gamma `$, $`\theta _\gamma `$, the mass normalization, the source location $`(x_s,y_s)`$, and the position angle of the source jet $`\varphi _\mathrm{s}`$. We applied the 8 constraints of Table 2 using a $`\chi ^2`$-minimization code. Because the constraint on the third-image flux is one-sided, and because of the way in which we have implemented that constraint, the allowed models fit the constraints perfectly ($`\chi ^2=0`$). The optimal value is $`\beta =2.05`$ and the $`2\sigma `$ bounds are $`1.95<\beta <2.28`$. This range brackets the isothermal value $`\beta =2`$. The upper bound is enforced by the observed jet orientations. The lower bound is enforced by the upper limit on the flux of the third image.

The lower bound is robust, even if our upper limit on $`\mu _3/\mu _\mathrm{A}`$ is too strong due to the possible extinction of a third image by free-free absorption, or the lowering of its peak flux density by scatter-broadening. For example, if we weaken the limit on the third-image flux density by a factor of 10 ($`\mu _3/\mu _\mathrm{A}<0.01`$), the lower bound changes only slightly: $`1.91<\beta <2.28`$.

Interestingly, there is a different (weaker) lower bound on $`\beta `$ that results from the relative positions of A and B with respect to the lens center. For $`\beta <1.83`$, the radial critical curve is sufficiently large to encompass B. In such models, the parity of B is reversed, and a brighter third image is predicted in a location outside the radial critical curve; image B becomes the central odd image in a three-image system. Such models can be ruled out.

Next, we investigate the possibility of a constant-density core. Given the preceding results, we assume that the density profile at large radii is isothermal and derive an upper bound on the size of a hypothetical constant-density core. We adopt a two-dimensional surface density of the form

$$\kappa (R)=\frac{b}{2\sqrt{R^2+R_c^2}},$$
(3)
where $`R`$ is the two-dimensional coordinate. The constraint on the third-image flux requires $`R_c<1.6`$ mas, corresponding to 0.8% of the optical radius of the lens galaxy, or $`9h^1`$ parsecs at $`z=1`$.111Here and elsewhere, we compute angular-diameter distances assuming a flat $`\mathrm{\Omega }_\mathrm{m}=0.3`$ cosmology, giving $`5.61h^1`$ kpc per arc second.

Finally, we explore models with central cusps. We use the cusped models introduced by Muñoz, Kochanek, & Keeton (2001), in which the mass density is

$$\rho (r)=\frac{\rho _0}{(r/r_b)^\beta [1+(r/r_b)^2]^{(\eta \beta )/2}}.$$
(4)
The density falls off with a logarithmic slope of $`\eta `$ for $`rr_b`$, and $`\beta `$ for $`rr_b`$, where $`r_b`$ is the “break radius.” This model is convenient because the lensing deflection can be computed analytically.

First, we assume that the mass distribution at large radii is isothermal ($`\eta =2`$). We use the constraints in Table 2 to find the allowed region in the two-dimensional space spanned by $`\beta `$ and $`r_b`$. As before, the allowed models fit the constraints perfectly ($`\chi ^2=0`$). The results are shown in the left panel of Fig. 5. The $`1\sigma `$ allowed region is shaded. For large $`r_b`$, the overall mass distribution must be nearly isothermal, as expected from the preceding results. Shallower inner cusps are allowed only when $`r_b\stackrel{<}{}0\stackrel{}{\mathrm{.}}02`$, corresponding to 10% of the optical radius, or $`100h^1`$ parsecs at $`z=1`$. As $`r_b`$ approaches zero, the constraint on $`\beta `$ widens rapidly. At $`r_b1.6`$ mas the range of allowed exponents encompasses $`\beta =0`$, which is the limit of a constant-density core derived in the previous section.

The allowed region of cusp parameters is qualitatively different if we allow the mass profile of the galaxy at large radii to be steeper than isothermal. The upper and lower right panels of Fig. 5 show the results for the cases $`\eta =3`$ and $`\eta =4`$. Note that both NFW (Navarro, Frenk, & White, 1997) and Hernquist (1990) profiles, with $`\beta =1`$ central cusps and asymptotic slopes of $`\eta =3`$ and $`\eta =4`$ respectively, are not permitted as global models of the mass distribution, because their inner cusps are too shallow. Unlike the case of $`\eta =2`$, the range permitted for the central cusp exponent $`\beta `$ does not widen as $`r_b`$ shrinks. The central image and jet orientation constraints instead pinch off the permitted parameter space to set a minimum break radius: $`r_b>0\stackrel{}{\mathrm{.}}12`$ ($`r_b>0\stackrel{}{\mathrm{.}}2`$) for $`\eta =3`$ ($`\eta =4`$). Furthermore, for both $`\eta =3`$ and $`\eta =4`$, the allowed values of $`\beta `$ bracket the isothermal value $`\beta =2`$, regardless of the break radius. The non-isothermal profile at large radii must be compensated by a nearly isothermal profile within a fairly large break radius. Again, we see that if the mass profile obeys a power law over the range of radii between images A and B, then the favored distribution is close to isothermal.

### 4.2 Case 2: Component C is a third image

For the three-image scenario, we use the model constraints given in Table 3. In this case, the lens galaxy position is taken from the optical measurement of Winn et al. (2002), with error bars enlarged to 23 mas (0.5 pixel) to allow for systematic errors.222The smaller error bars of 9 mas quoted by Winn et al. (2002) assume a particular form for the point-spread function and the galaxy profile. Ordinarily, this would not introduce much additional error, but in this case component B and the lens galaxy are both faint and merged in the HST image. The value adopted for $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$ is based on the fit to the free-free absorption model described in § 3.5. We define the goodness-of-fit parameter as

$$\chi ^2=\underset{i=\mathrm{B},\mathrm{C},\mathrm{Gal}}{}\left[\frac{(𝐑_i^{}𝐑_i)^2}{(\mathrm{\Delta }𝐑_i)^2}\right]+\underset{i=\mathrm{A},\mathrm{B}}{}\left[\frac{(\varphi _i^{}\varphi _i)^2}{(\mathrm{\Delta }\varphi _i)^2}\right]+\frac{(\mu _{\mathrm{BA}}^{}\mu _{\mathrm{BA}})^2}{(\mathrm{\Delta }\mu _{\mathrm{BA}})^2}+\frac{(\mu _{\mathrm{CA}}^{}\mu _{\mathrm{CA}})^2}{(\mathrm{\Delta }\mu _{\mathrm{CA}})^2},$$
(5)
in which primed quantities are those determined by the model, and unprimed quantities are those observed.

For a scale-free power-law profile in an external shear field, we have three degrees of freedom. The model does a good job of reproducing the configuration, with $`\chi ^2=2.2`$. The logarithmic slope of the density profile in the optimal model is $`\beta =1.91\pm 0.02`$ ($`2\sigma `$). If we take the radial density profile to be isothermal with a constant-density core, as in Eq. 3, we can also achieve a reasonable fit ($`\chi ^2=3.0`$), with the optimal value $`r_c=5.5`$ milliarcseconds, or $`31h^1`$ parsecs, and a $`2\sigma `$ confidence region of 3.6–7.8 mas.

When we apply the cuspy lens model of Eq. 4, we find that the allowed region of cusp parameters is much more tightly constrained than they are in the two-image lens models. Figure 6 shows the results for models in which the radial density distribution at large radii is isothermal (left panel), or obeys a power law with logarithmic slope $`\eta =3`$ (upper right), or $`\eta =4`$ (lower right). Both boundaries of the allowed region (shaded) are enforced by the constraint on the third-image flux. In particular, for $`\beta >2`$, no third image is produced. For large $`r_b`$, the cusp must conform to the value $`\beta =1.91`$ derived above. As $`r_b`$ becomes smaller, the inner mass distribution must become shallower to compensate for the steeper mass distribution at large radii.

## 5 Discussion and summary of conclusions

For the purpose of finding central odd images, the best systems are radio-loud two-image quasars with a large magnification ratio between the two bright images. Radio loudness is important because the radio components can be observed through the lens galaxy, whereas optical components are likely to be hidden by the lens galaxy. Free-free absorption is potentially significant but can be overcome by observing at high frequencies. Two-image systems are better than four-image systems, as emphasized recently by Keeton (2002). This is because the projected source position in four-image systems is usually closer to the lens center, where the central image magnification is smallest. In addition, four-image systems are subject to a larger magnification bias than two-image systems (Wallington & Narayan, 1993; King & Browne, 1996; Keeton, Kochanek, & Seljak, 1997; Rusin & Tegmark, 2001; Finch et al., 2002), leading to a smaller intrinsic source flux and a fainter odd image. Among two-image systems, the asymmetric doubles are favorable because the source position is close to the radial caustic, producing the brightest possible central image for a given mass distribution.

The subject of this paper, PMN J1632–0033, is exactly this type of system: an asymmetric two-image radio lens. For this reason, it seemed reasonable that the third component found by Winn et al. (2002) could be an example of the elusive and long-sought central odd images. However, we have shown that the third component has a spectral index that differs by $`3\sigma `$ from those of the two bright quasar images.

We have argued that the two most plausible interpretations are that the central component is emission from the lens galaxy, or else that it is a third image with the extra complication that its spectral index has been modified by free-free absorption. The data are quantitatively consistent with either scenario. Lens galaxies are rarely radio-loud, but the flux density of the component is consistent with typical radio powers of active galactic nuclei. No central image has ever been securely identified for a radio lens, nor has free-free absorption by a lens galaxy ever been reported (with one possible exception, class B0128+437; I. Browne 2002, private communication), but one might expect the two phenomena to be related due to the hot and dense conditions near galaxy centers.

Given the new VLBI data, we have investigated models for the mass distribution of the lens galaxy under various assumptions. Whether or not component C is a third image, we showed that if the mass distribution obeys a power law $`\rho r^\beta `$, then it must be nearly isothermal. If C is the lens galaxy then $`1.95<\beta <2.28`$. If C is a third image then the profile is tightly constrained to be only slightly shallower: $`1.89<\beta <1.93`$. We also derived limits on a constant-density core radius, or on the break radius and inner exponent of a central cusp.

These results can be added to the growing collection of evidence that early-type galaxies have nearly isothermal mass profiles. In the few other gravitational lens systems where it is possible to measure the radial density profile, the results favor isothermal profiles (Kochanek, 1995; Cohn et al., 2001; Muñoz, Kochanek, & Keeton, 2001; Rusin et al., 2002; Koopmans & Treu, 2002), although there are some possible counter-examples (Treu & Koopmans, 2002; Chen, Kochanek, & Hewitt, 1995). Moreover, the findings are consistent with studies of early-type galaxies based on stellar dynamics (see, e.g., Gerhard et al., 2001; Rix et al., 1997) and X-ray halos (see, e.g., Fabbiano et al., 1989).

To further elucidate the nature of the central component and the central density profile of the lens galaxy, it would help to have a better optical image of the lens galaxy. The current HST image provides only a weak detection of the lens galaxy. A better image could be used to test more precisely whether radio component C is coincident with the center of the lens galaxy, and to study its surface brightness profile in detail. For example, with a better measurement of the galaxy light profile, one might be able to explore the balance between the luminous and dark matter in this system.

Even more helpful would be sensitive observations at high radio frequencies ($`\stackrel{>}{}`$30 GHz). As shown in Fig. 4, the power-law and the free-free absorption fits diverge significantly at high frequencies, where free-free absorption becomes negligible. For free-free absorption, the ratio $`S_\mathrm{C}(\nu )/S_\mathrm{A}(\nu )`$ must converge to a constant value. For lens galaxy emission, the spectrum need not follow a power law at the high frequencies, but it would require an unlikely coincidence for the intrinsic radio spectrum of the lens galaxy to turn over at high frequencies and become identical to the radio spectra of A and B. Observations at low radio frequencies ($`\stackrel{<}{}`$1 GHz) could also distinguish the hypotheses, in principle, but it would be harder to achieve the necessary angular resolution and signal-to-noise ratio. High frequency observations present their own challenges, including high noise levels and susceptibility to poor atmospheric conditions, but the reward would be a definitive discovery of one of the long-sought but never-seen central images of a gravitational lens.

We thank the anonymous referee for constructive criticism. We also thank I. Browne and L. Koopmans for interesting discussions about radio propagation effects. J.N.W. is supported by an Astronomy & Astrophysics Postdoctoral Fellowship, under NSF grant AST-0104347. C.S.K. is supported by NASA ATP grant NAG5-9265.


format: text

text:
# Order parameter node removal in the d-wave superconductor 𝑌𝐵𝑎₂𝐶𝑢₃𝑂_{7-𝑥} under magnetic field

## I

###### Acknowledgements.

We are indebted to Malcom Beasley for pointing out to one of us (G.D.) reference b21  on the effects of Bean Livingston barrier as seen in a tunneling experiment, to Dr. H. Castro and A. Kohen for helpful conversations. This work was supported by the Heinrich Herz-Minerva Center for High Temperature Superconductivity, by the Israel Science Foundation and by the Oren Family Chair of Experimental Solid State Physics.


text:
# Quantum Andreev map: A paradigm of quantum chaos in superconductivity

(March 2003)

## Abstract

We introduce quantum maps with particle-hole conversion (Andreev reflection) and particle-hole symmetry, which exhibit the same excitation gap as quantum dots in the proximity to a superconductor. Computationally, the Andreev maps are much more efficient than billiard models of quantum dots. This makes it possible to test analytical predictions of random-matrix theory and semiclassical chaos that were previously out of reach of computer simulations. We have observed the universal distribution of the excitation gap for large Lyapunov exponent and the logarithmic reduction of the gap when the Ehrenfest time becomes comparable to the quasiparticle dwell time.

One-dimensional (1D) quantum mechanical models with a chaotic classical limit were first studied by Casati, Chirikov, Ford, and Izrailev in 1979 Cas79 . These models have since developed into one of two paradigms of quantum chaos Izr90 ; Haa92 . The other paradigm is the 2D billiard of irregular shape Boh84 . Two is the lowest number of dimensions for non-integrable (chaotic) dynamics in autonomous systems, since a single constant of motion is sufficient for integrability in 1D. The 1D models get around this constraint through a periodically time-dependent external force (“kick”), which eliminates the energy as constant of motion — but still conserves the quasi-energy (analogously to quasi-momentum conservation in a periodic lattice). The two paradigms share a common set of phenomena in the fields of quantum chaos and localization Fis82 ; Alt96 ; Bor96 ; Fra97 .

The combination of chaos and superconductivity produces an entirely new phenomenology, notably the appearance of an excitation gap as a signature of quantum chaos Mel96 . The paradigm common to most of the literature is the 2D billiard connected to a superconductor Gut02 , introduced under the name “Andreev billiard” in Ref. Kos95 . The name refers to the Andreev reflection which occurs at the interface with the superconductor, where an electron at energy $`\epsilon `$ above the Fermi level is converted into a hole at energy $`\epsilon `$ below it.

From the point of view of computational efficiency, compact quantum maps such as the kicked rotator Cas79  (a particle confined to a circle and driven periodically in time with a strength that depends on its position) are much more powerful than 2D models such as billiards. Indeed, there exists a highly efficient diagonalization technique that works only for maps Ket99 . The lack of a 1D map for quantum chaos with superconductivity has hindered the numerical test of a variety of analytical predictions Vav01 ; Ost01 ; Lod98 ; Lam01 ; Tar01 ; Ada02 ; Sil02 ; Vav02 . Most notably, numerical efforts have not been able to distinguish the conflicting predictions Mel96  of random-matrix theory (RMT) and the semiclassical Bohr-Sommerfeld (BS) quantization: RMT predicts an excitation gap at the Thouless energy while BS gives an exponentially vanishing density of states without a true gap. Recent analytical work Tar01 ; Ada02 ; Sil02 ; Vav02 ; Lod98  has predicted that diffraction creates a gap in the BS density of states at the inverse Ehrenfest time. This has never been seen in computer simulations, because the Ehrenfest time scales logarithmically with the system size and is usually far too small to play a role. For these reasons there is a real need for something like a “quantum Andreev map”. Does it exist? And if it does, can it be simulated more efficiently than the Andreev billiard? These are the issues addressed in this paper.

We show how to construct quantum Andreev maps out of any conventional quantum map (not necessarily chaotic), in much the same way as any normal billiard can be turned into an Andreev billiard by coupling it to a superconductor. The construction is guided by the classical electron and hole dynamics on the Poincaré surface of section of an Andreev billiard. The Andreev kicked rotator is a particular example of such an Andreev map. We certify that it possesses the phenomenology of the Andreev billiards and search for predictions of RMT and semiclassics. We leave for future investigations the application of the Andreev map to other kicked models (possibly with a different phenomenology), such as the kicked top Haa92  and the Fermi-Ulam model Jain93 .

A quantum map is represented by the Floquet operator $`F`$, which gives the stroboscopic time evolution $`u(p\tau _0)=F^pu(0)`$ of an initial wave function $`u(0)`$. (We set the stroboscopic period $`\tau _0=1`$ in most equations.) The unitary $`M\times M`$ matrix $`F`$ has eigenvalues $`\mathrm{exp}(i\epsilon _m)`$, with the quasi-energies $`\epsilon _m(\pi ,\pi )`$ (measured in units of $`\mathrm{}/\tau _0`$). This describes particle excitations in a normal metal. We also need hole excitations. A particle excitation with energy $`\epsilon _m`$ (measured relatively to the Fermi level) is identical to a hole excitation with energy $`\epsilon _m`$. This means that hole excitations in a normal metal have Floquet operator $`F^{}`$ and wave function $`v(p)=(F^{})^pv(0)`$.

Particles and holes are coupled by connecting the normal metal via a lead to a superconducting reservoir. On a Poincaré surface of section, the lead is represented by a spatially localized region in which Andreev reflection converts electrons into holes and vice versa, with phase shift $`i`$. (A weak energy dependence of this phase shift is ignored for simplicity, but can be accounted for straightforwardly.) Analogously, for the quantum Andreev map we assume that Andreev reflections occur whenever an excitation ends up in a certain subspace of Hilbert space. This subspace $`n_1,n_2,\mathrm{}n_N`$ consists of $`N`$ out of $`M`$ states in a chosen basis and corresponds to a lead with $`N`$ propagating channels. The $`N\times M`$ matrix $`P`$ projects onto the lead. Its elements are $`P_{nm}=1`$ if $`m=n\{n_1,n_2,\mathrm{}n_N\}`$ and $`P_{nm}=0`$ otherwise. The dwell time of a quasiparticle excitation in the normal metal is $`\tau _{\mathrm{dwell}}=M/N`$, equal to the mean time between Andreev reflections. The fact that Andreev reflections only occur at multiples of the stroboscopic time $`\tau _0`$ is technically convenient, and should be physically irrelevant for $`\tau _{\mathrm{dwell}}\tau _0`$.

Putting all this together we construct the quantum Andreev map from the matrix product

$$=𝒫\left(\begin{array}{cc}F& 0\\ 0& F^{}\end{array}\right),𝒫=\left(\begin{array}{cc}1P^\mathrm{T}P& iP^\mathrm{T}P\\ iP^\mathrm{T}P& 1P^\mathrm{T}P\end{array}\right).$$
(1)
(The superscript “T” indicates the transpose of a matrix.) The particle-hole wave function $`\mathrm{\Psi }=(u,v)`$ evolves in time as $`\mathrm{\Psi }(p)=^p\mathrm{\Psi }(0)`$. The Floquet operator can be symmetrized (without changing its eigenvalues) by the unitary transformation $`𝒫^{1/2}𝒫^{1/2}`$, with

$$𝒫^{1/2}=\left(\begin{array}{cc}1(1\frac{1}{2}\sqrt{2})P^\mathrm{T}P& i\frac{1}{2}\sqrt{2}P^\mathrm{T}P\\ i\frac{1}{2}\sqrt{2}P^\mathrm{T}P& 1(1\frac{1}{2}\sqrt{2})P^\mathrm{T}P\end{array}\right).$$
(2)
In order to establish the correspondence of the 1D quantum Andreev maps to 2D Andreev billiards, we examine the spectral properties of the map. The Floquet operator $``$ possesses a particle-hole symmetry which entails that its $`2M`$ eigenvalues $`\mathrm{exp}(i\epsilon _m)`$ come in inverse pairs. This symmetry is the analogue of the particle-hole symmetry in Andreev billiards, in which excitation energies $`\pm \epsilon `$ occur symmetrically around the Fermi level. To avoid double-counting of levels, we restrict the quasi-energy to the interval $`(0,\pi )`$. The excitation spectrum of particles and holes consists of the $`M`$ quasi-energies in this interval, and the mean level spacing $`\pi /M`$ is twice as small as the level spacing $`\delta =2\pi /M`$ for particles and holes separately. The energy scale for the proximity-induced excitation gap is the Thouless energy $`E_\mathrm{T}=N\delta /4\pi =N/2M=1/(2\tau _{\mathrm{dwell}})`$.

The quantization condition $`\mathrm{det}(e^{i\epsilon })=0`$ can be written equivalently as

$$\mathrm{det}[1+S(\epsilon )S^{}(\epsilon )]=0,$$
(3)
in terms of the $`N\times N`$ scattering matrix Oss02

$$S(\epsilon )=P[e^{i\epsilon }F(1P^\mathrm{T}P)]^1FP^\mathrm{T}.$$
(4)
Eq. (3) for the Andreev map has the same form as for the Andreev billiard Fra96 , but there $`S`$ is given in terms of a Hamiltonian $`H_0`$ instead of a Floquet operator $`F`$. The two approaches become entirely equivalent in the context of RMT, when $`H_0`$ is chosen from one of the Gaussian ensembles and $`F`$ is chosen from one of the circular ensembles proof . They are also equivalent in the semiclassical limit, when the billiard can be represented by a Poincaré map which can be quantized approximately Bog92 .

In the mean-field limit $`MN1`$, RMT predicts a hard gap in the excitation spectrum of size $`E_{\mathrm{RMT}}=\gamma E_\mathrm{T}`$ (with $`\gamma =2^{3/2}(\sqrt{5}1)^{5/2}=0.60`$), and above the gap a square-root singularity in the density of states $`\rho (\epsilon )=\pi ^1\mathrm{\Delta }^{3/2}(\epsilon E_{\mathrm{RMT}})^{1/2}`$ (with $`\mathrm{\Delta }=0.068N^{1/3}\delta `$) Mel96 . Sample-to-sample fluctuations of the lowest excitation energy $`\epsilon _0`$ around the mean-field gap have been calculated in Refs. Vav01 ; Ost01 . A universal probability distribution was predicted for the rescaled energy $`x=(\epsilon _0E_{\mathrm{RMT}})/\mathrm{\Delta }`$. While the mean-field prediction of RMT has been tested numerically in an Andreev billiard Mel96 , the numerical error bars are too large to extract the predicted universal gap fluctuations.

To demonstrate the efficiency of the quantum Andreev maps, we specialize to the quantum kicked rotator. The Floquet operator is Izr90

$`F`$ $`=`$ $`\mathrm{exp}\left(i{\displaystyle \frac{\mathrm{}\tau _0}{4I_0}}{\displaystyle \frac{^2}{\theta ^2}}\right)\mathrm{exp}\left(i{\displaystyle \frac{KI_0}{\mathrm{}\tau _0}}\mathrm{cos}\theta \right)`$ (5)
$`\times \mathrm{exp}\left(i{\displaystyle \frac{\mathrm{}\tau _0}{4I_0}}{\displaystyle \frac{^2}{\theta ^2}}\right),`$
with $`I_0`$ the moment of inertia of the particle and $`K`$ the (dimensionless) kicking strength. The particle moves freely along the circle for half a period, is then kicked with a strength $`K\mathrm{cos}\theta `$, and proceeds freely for another half period. The transition from classical to quantum behavior is governed by the effective Planck constant $`\mathrm{}_{\mathrm{eff}}\mathrm{}\tau _0/I_0`$. Since we would like to compare the kicked rotator to a chaotic billiard, without localization, we follow the usual procedure of quantizing phase space on the torus $`\theta ,p(0,2\pi )`$, rather than on a cylinder, with $`p=i\mathrm{}_{\mathrm{eff}}/\theta `$, the dimensionless angular momentum Izr90 . For $`\mathrm{}_{\mathrm{eff}}=2\pi /M`$, with integer $`M`$, the Floquet operator is an $`M\times M`$ unitary symmetric matrix. In angular momentum representation it has elements

$`F_{kk^{}}=e^{(i\pi /2M)(k^2+k^2)}(UQU^{})_{kk^{}},`$ (6a)
$`U_{kk^{}}=M^{1/2}e^{(2\pi i/M)kk^{}},`$ (6b)
$`Q_{kk^{}}=\delta _{kk^{}}e^{(iMK/2\pi )\mathrm{cos}(2\pi k/M)}.`$ (6c)
Upon increasing $`K`$ the classical dynamics varies from fully integrable ($`K=0`$) to fully chaotic \[$`K7`$, with Lyapunov exponent $`\lambda \mathrm{ln}(K/2)`$\]. For $`K<7`$ stable and unstable motion coexist (a so-called mixed phase space).

To introduce the Andreev reflection we use a projection operator which is diagonal in $`p`$-representation,

$$(P^\mathrm{T}P)_{kk^{}}=\delta _{kk^{}}\times \{\begin{array}{c}1\mathrm{if}LkL+N1,\hfill \\ 0\mathrm{otherwise}.\hfill \end{array}$$
(7)
(We checked that similar results are obtained when $`P`$ is diagonal in $`\theta `$-representation.) The position $`L`$ of the lead to the superconductor is arbitrary. The Floquet operator $``$ of the “Andreev kicked rotator” is then obtained by inserting Eqs. (6) and (7) into Eq. (1). We apply the symmetrization (2), so $``$ is a unitary symmetric matrix. The real symmetric matrix $`\frac{1}{2}(+^{})`$ can be diagonalized efficiently with $`𝒪(M^2\mathrm{ln}M)`$ operations (and not $`𝒪(M^3)`$ as with standard methods) by means of the Lanczos technique, if the multiplication with the matrix $`U`$ is carried out with the help of the Fast-Fourier-Transform algorithm Ket99 . The eigenvalues $`\mathrm{cos}\epsilon _m`$ uniquely determine the quasi-energy $`\epsilon _m(0,\pi )`$.

In the inset to Fig. 1 we show the density of states $`\rho (\epsilon )`$ for system size $`M=8192`$, kicking strength $`K=45`$ (strongly chaotic dynamics), and several widths $`N`$ of the lead to the superconductor. The density of states has been averaged over 250 different positions of the lead. The data points fall on top of the RMT prediction Mel96  without any adjustable parameter. Reducing the kicking strength down to $`K=1.2`$, one enters the regime of mixed classical dynamics. We see that the gap disappears, as predicted in Ref. Sch99 .

To test RMT beyond the mean-field limit, we study the statistical fluctuations of the gap. The main panel of Fig. 1 shows the probability distribution of the smallest eigenvalue $`\epsilon _0`$ in the chaotic regime. To improve statistics we sampled 6000 different positions of the lead. We rescaled the energy $`x=(\epsilon _0E_{\mathrm{RMT}})/\mathrm{\Delta }`$, as prescribed by Ref. Vav01 . Good agreement is observed with the universal scaling distribution Tra94 , again without any adjustable parameters.

It is predicted theoretically that deviations from RMT should appear if the Ehrenfest time $`\tau _\mathrm{E}=\lambda ^1\mathrm{ln}M`$ is no longer small compared to $`\tau _{\mathrm{dwell}}`$. For $`\tau _\mathrm{E}\tau _{\mathrm{dwell}}`$ the semiclassical Bohr-Sommerfeld approximation Mel96 ; Sch99  should be valid, with a diffraction induced gap of order $`\mathrm{}/\tau _\mathrm{E}`$ Lod98 . To search for these deviations from RMT we consider rotators with smaller kicking strengths (but still in the fully chaotic regime), thus smaller Lyapunov exponent, and much larger $`M`$.

In Fig. 2 we show the density of states for $`M=131072`$ and $`K=14`$. Strong deviations from the RMT prediction are clearly visible. Also plotted is the result of a BS calculation Mel96 , in which we slightly smoothed the singular delta functions. This approximation agrees better with the exact result. Most remarkably, it reproduces the three distinct peaks in the density of states, which now can be identified with trajectories of certain lengths. All trajectories with lengths that are odd multiples of $`\tau _{\mathrm{dwell}}=5`$ contribute to the peak at $`\epsilon /E_\mathrm{T}=\pi `$, odd multiples of $`4`$ contribute at $`\epsilon /E_\mathrm{T}=5\pi /4`$, and odd multiples of $`3`$ contribute at $`\epsilon /E_\mathrm{T}=5\pi /3`$.

A systematic reduction of the excitation gap is observed upon increasing the ratio $`\tau _\mathrm{E}/\tau _{\mathrm{dwell}}`$, as shown in Fig. 3. The main panel is a semi-logarithmic plot of $`\epsilon _0/E_T`$ as a function of $`M[2^9,2^{19}]`$, for $`M/N=\tau _{\mathrm{dwell}}=5`$ and $`K=14`$. Existing theories Vav02 ; Sil02  predict a linear initial decrease of $`\epsilon _0`$ with $`\mathrm{ln}M`$ at fixed $`\tau _{\mathrm{dwell}}=M/N`$. We fit the data to the prediction of Vavilov and Larkin Vav02 ,

$$\frac{\epsilon _0}{E_{\mathrm{RMT}}}=1\frac{\alpha }{2\lambda \tau _{\mathrm{dwell}}}\left(\mathrm{ln}M2\mathrm{ln}\frac{M}{N}\alpha ^{}\right).$$
(8)
We find $`\alpha =0.59`$ and $`\alpha ^{}=3.95`$. Once $`\alpha `$ and $`\alpha ^{}`$ are extracted, no free parameter is left, and the resulting curve, shown with a solid line in the inset to Fig. 3, correctly reproduces the parametric dependence on $`\tau _{\mathrm{dwell}}=M/N`$ at fixed $`M`$. As a further check, we tried the slightly different expression

$$\frac{\epsilon _0}{E_{\mathrm{RMT}}}=1\frac{\alpha }{2\lambda \tau _{\mathrm{dwell}}}\left(\mathrm{ln}M\alpha ^{\prime \prime }\right),$$
(9)
with $`\alpha ^{\prime \prime }=\alpha ^{}+2\mathrm{ln}5`$. The resulting curve (dotted line in the inset to Fig. 3) shows significant deviations from the numerical data. We conclude that Eq. (8) gives the correct parametric dependence of the Andreev gap.

A discrepancy remains in the value of the numerical coefficients. While the coefficient $`\alpha ^{}`$ is model dependent, the prefactor $`\alpha `$ is expected to be universal. Our numerics gives $`\alpha =0.59\pm 0.08`$, in between the two competing predictions $`\alpha =0.23`$ of Ref. Vav02  and $`\alpha =2`$ of Ref. Sil02 .

In conclusion, we have constructed a quantum map that accounts for the presence of superconductivity. The “Andreev kicked rotator” introduced above, has been shown to be equivalent to the Andreev billiards studied so far. Owing to the fact that it is one-dimensional rather than two-dimensional, it is much more efficient computationally, which permits to observe two theoretical predictions that are currently out of reach of billiard simulations: The universal gap fluctuations for large Lyapunov exponent and the logarithmic reduction of the gap for small Lyapunov exponent. We foresee that the Andreev kicked rotator on a cylinder (instead of on a torus) can be an equally effective tool to study the interplay of superconductivity and localization.

We have benefitted from discussions with İ. Adagideli and J. Tworzydło. This work was supported by the Dutch Science Foundation NWO/FOM.


text:
# Secure communication using mesoscopic coherent states

(Revised version: 04/17/03.)

## Abstract

We demonstrate theoretically and experimentally that secure communication using intermediate-energy (mesoscopic) coherent states is possible. Our scheme is different from previous quantum cryptographic schemes in that a short secret key is explicitly used and in which quantum noise hides both the bit and the key. This encryption scheme can be optically amplified. New avenues are open to secure communications at high speeds in fiber-optic or free-space channels.

For the encryption of data with perfect secrecy shannon  that cannot be broken with any advance in technology, one may in principle employ one-time pad with secret key obtained by the BB84 bennett-brassard  quantum cryptographic technique for key expansion. Such an approach is possible alamos , however, it is slow and inefficient because the key length needs to be as long as the data, and it also requires a nearly ideal quantum communication line that is difficult to obtain in long distance commercial systems such as the Internet core. On the other hand, for both military and commercial applications, there are great demands for secret communications that are fast and secure but not necessarily perfectly secure. (There are many practical issues, human as well machine based, that would make theoretical perfect security in specific models not so important in real life schneier ). In the following, a new scheme based on ideas similar to those of Ref. yuenAKE  is described for secure data encryption that can be operated at optical speeds with conventional optical technology, and a prototype experimental implementation is presented. In this scheme, a short secret key is classically extended and then used to encrypt data in a way that the quantum noise of the coherent states protects both the data and the key.

The following line of reasoning describes the ideas yuen\_tobe  that led to the development of our present kind of quantum cryptographic schemes. One crucial element for obtaining security in BB84 involves the detection of small intrusion on weak signals, which is difficult to achieve in a network environment. This problem would be alleviated if quantum signal sets of higher energy are selected for different bit values by a secret key shared between the sender (Alice) and the receiver (Bob). It is important to remember that some shared secret key is needed in BB84 for message authentication during protocol execution. The resulting scheme is acceptable as key expansion if the new key is secure even if the shared secret key is known to the attacker after the user communications are completed. When a secret key is used to identify the signal set, it would be a secret CDMA (Code Division Multiple Access) scheme classically, which does not allow key expansion because the user and the attacker have the same observation. We would discuss elsewhere how a corresponding KCQ (Keyed CDMA in Quantum Noise) scheme can be used to obtain key expansion in the quantum case. In this paper, we are concerned with the use of KCQ for data encryption.

There are two basic problems with classical encryption that does not employ the inefficient one-time pads. One is that the total data uncertainty $`H(𝐗)`$ given observation $`𝐘`$ is bounded above by the key uncertainty, $`H(𝐗)H(𝐊)`$ shannon . The other is that the key $`𝐊`$ may be found by a known-plaintext attack when the eavesdropper (Eve) knows the output-input pairs $`(𝐘,𝐗)`$ for some data length. In our scheme, $`H(𝐗)`$ is not bounded by $`H(𝐊)`$ because Eve cannot have the observation $`𝐘`$ that Bob obtained via the optimal quantum measurement utilizing the key $`𝐊`$. To extract information from even a full copy of the quantum signal without knowing $`𝐊`$, Eve has to make a sub-optimal measurement that would yield information on all possible signal sets for the purpose of either estimating $`𝐗`$ or finding $`𝐊`$ from a known-plaintext attack. As a result, Bob has a better channel/observation than Eve. Also, in contrast to classical cryptography, one can prove the security of our scheme against ciphertext-only attacks, although only individual attacks are described in this paper. One can show that, in a properly designed system, even an exponentially powerful search with known-plaintext attacks cannot succeed because Eve does not have $`𝐘`$ as above. Practically, our scheme can run at high speeds because there is no need for a long key $`𝐊`$.

Consider the following scheme in which each data bit is encoded into a coherent state of an infinite-dimensional space, referred to as a “qumode”. As in yuenAKE , there are $`M`$ possible states

$`|\alpha _0(\mathrm{cos}\theta _l+i\mathrm{sin}\theta _l),\theta _l=2\pi l/M`$ (1)
for a real $`\alpha _0`$ and $`l\{0,\mathrm{},M1\}`$, forming $`M/2`$ pairs $`\{l,l+M/2\}`$. A seed key $`𝐊`$ is used to drive an encryption mechanism whose output is a much longer running key $`𝐊^{}`$ that is used to determine, for each qumode carrying bit $`b`$ $`(=0,1)`$, which pair of signals (signal set) is to be used. Each pair may be macroscopically distinguishable since the inner product of any two basis states is $`\mathrm{exp}(2|\alpha _0|^2)`$. For large $`M`$, a lower bound on the obtained mean-square error $`(\delta \theta )^2`$ yuen2  that goes as $`1/|\alpha _0|^2`$ shows that asymptotically when $`M|\alpha _0|`$ the attacker’s error probability $`P_e^E`$ tends to $`1/2`$, the guessing level, in an individual attack on the data bit $`b`$. That this result holds in the limit $`M\mathrm{}`$ for fixed $`|\alpha _0|`$ is intuitively obvious. A two-mode coherent-state realization similar to Eq. (1), with $`|\alpha _0\mathrm{cos}\theta _l|\alpha _0\mathrm{sin}\theta _l`$, can also be used and the modes can be interpreted as ones of polarization, time, frequency, or whatever.

Numerical calculation of the optimal POVM for individual attack on bit discrimination for the $`M`$-ry case has shown nu-pavia  that the minimum probability of error $`P_e^E`$ for an eavesdropper can be made arbitrarily close to $`1/2`$ for a given coherent-state amplitude $`\alpha `$. The value $`P_e^E1/2`$ for a fixed average number of photons $`|\alpha |^2`$ is achieved by increasing the number of levels $`M`$. As shown in Figure 1, $`P_e^E`$ goes very fast to the asymptotic pure-guessing limit of $`1/2`$ as $`M`$ increases.

The above POVM calculation demonstrates that in this scheme an eavesdropper cannot obtain the bits sent regardless of the precision of her devices. The optimal POVM gives the maximum amount of information she could obtain from the sequence of physical signals sent without knowing the key. This uncertainty is due to the quantum noise of light and cannot be overcome with one’s precision capabilities. Bob, on the other hand, by knowing the key can extract information with greater precision. His decision has to be made only between two nearly orthogonal states in the same basis defined by a given $`𝐊^{}`$. His probability of error is helstrom  $`P_e^B=\frac{1}{2}\left(1\sqrt{1e^{2T|\alpha |^2}}\right)`$, where $`T`$ is the transmissivity of the channel. For large values of $`|\alpha |`$ the minimum probability of error $`P_e^B`$ is negligible, which makes possible an excellent signal recovery by the legitimate receiver. The case of collective attacks is more complicated and cannot be discussed here, in large part because there is no meaningful approach for evaluating the optimal bit error, even in the classical case. However, the entropy bound (Holevo’s theorem) could be used to show the ideal nature of this scheme for the criterion of data entropy.

The attacker can also try to find the key $`𝐊`$ based on her copy of the quantum signals, with or without some known-plaintext (data) corresponding to the signals. Even in a known-plaintext attack, the signal quantum fluctuation would yield, from the number of possibilities in each qumode, an exponential number of possible $`𝐊^{}`$ in a sequence of data bits. To identify $`𝐊`$ from such noisy observations of possible $`𝐊^{}`$ would involve an exponential search, which can always be launched against a key in known-plaintext attacks but which is currently believed to be proven impossible computationally even with a quantum computer. Note that the attacker has a much more difficult job of estimating the signal pair from $`M/2`$ possible ones, than the user who tries to discriminate 2 possible known states. A detailed quantitative treatment would be given elsewhere. It is important to note that, in the case of classical cryptography, if the running key $`𝐊^{}`$ is used directly as one-time pad on the data, the result is well known menezes  to be insecure against known-plaintext attacks. In our case, this attack is thwarted by quantum effects as explained above.

Consider an implementation of the above scheme as depicted in Fig. 2. This particular type of KCQ scheme will be refered to as $`\alpha \eta `$ (for coherent states and efficiency). In the encryption scheme for this new protocol, sketched on the left side in Fig. 2, Alice uses an explicit short secret key $`𝐊`$, extended to a longer key $`𝐊^{}`$ by another encryption mechanism such as a stream cipher, to modulate the parameters of a multimode coherent state. For the free-space implementation to be presented, the qumodes are the two orthogonal modes of polarization. In this case, Alice uses the running key $`𝐊^{}`$ to specify a polarization basis from a set of $`M/2`$ uniformly spaced two-mode bases spanning a great circle on the Poincaré sphere, as shown on the right in Fig. 2. Each basis consists of a polarization state and its antipodal state at an angle $`\pi `$ from it, representing the 0 and 1 bit value for that basis.

The message $`𝐗`$ is encoded as $`𝐘(𝐗,𝐊^{})`$. This mapping of the stream of bits onto points on the surface of the Poincaré sphere is the information to be shared by Alice (A) and Bob (B). Because of his knowledge of $`𝐊^{}`$, Bob is able to make a precise demodulation operation, producing the plaintext $`𝐗`$. He uses $`𝐊^{}`$ to apply the requisite polarization transformation to the received sequence of polarization states to return them to the linearly orthogonally polarized condition, representing the two original bits of the message $`𝐗`$. Bob’s demodulation is the inverse mapping transformation that was utilized by Alice.

Figure 3 shows a table-top experimental setup we have implemented as a proof-of-principle demonstration of this scheme. The modulation systems utilized both at the transmitter and receiver ends are electro-optic modulators (EOMs)(New Focus, model 4104), the laser (Toshiba, model TOLD9225M) operates at 670nm, and the detectors are single photon counting modules (Perkin-Elmer, model SPCM-AQR-16) with interference filters of 10nm bandwidth in front. A polarization beam splitter (PBS) is used at the receiver to discriminate between the orthogonal linear polarization states. Lenses (L) are used to optimize the beam Rayleigh range within the EOMs. A personal computer containing an interface card (National Instruments, model PCI-6111E) is used to control the EOM’s (digital-to-analog operation). The same card is also used for counting the output pulses from the detectors. In this configuration, a horizontally (H) or vertically (V) polarized light pulse representing bit 0 or 1 is generated and transformed into an elliptically polarized light state by application of the voltage $`V_k(k𝐊^{})`$. This voltage introduces a phase difference $`\mathrm{\Delta }\varphi _k`$ between the physical axes of the EOM, where $`\mathrm{\Delta }\varphi _k=\frac{\pi }{2}\frac{V_k}{V_\pi }+\varphi _0`$, and $`V_\pi `$ and $`\varphi _0`$ are specific to each modulator. This system operates, with bulk optics, at 200kHz rate for demonstration purposes and faster fiber based systems ($``$1GHz) are being implemented for free space as well as fiber channels.

Figure 4 shows a sequence of bits as received by Bob. The clear separation of the 0 and 1 histograms allow him to make bit decisions with no error. The same sequence of bits as seen by Eve are shown in Fig. 5, giving her a very high probability of error ($`P_e^E1/2`$) in bit decisions because of her lacking the key.

As an illustration of how the quantum noise of light can be utilized by the legitimate receiver on his behalf, in Fig. 6 we show the uncertainty in the polarization angle produced by a two-mode measurement for an average total photon number $`n38`$. Regions of low-variance values are seen around 0 and $`\pi /2`$ settings of the input polarization state \[cf. Fig. 6(c)\]. Such determination of the polarization angle with angle-dependent uncertainty shows that a higher degree of precision in angle determination can be achieved by an observer with a prior information on how to set the measuring apparatus than another who does not have this knowledge. From Fig. 6(c) one can observe that for PBS-axis orientations close to the incoming field polarization direction the uncertainty in determination of the angle is small. The polarization angle is obtained by averaging $`\mathrm{arctan}\sqrt{n_V/n_H}`$ over the occurrence of $`n_H`$ and $`n_V`$, the photon numbers reaching the detectors \[cf. Fig. 6(b)\]. The approximation that $`\theta `$ can be obtained through $`\mathrm{arctan}\sqrt{n_V/n_H}`$ is adequate for mesoscopic signals but becomes inadequate as $`n0`$ (quantum phase domain). Setting PBS without knowledge of this preferred orientation leads, in average, to a larger error in determination of the polarization angle. On the other hand, an observer setting his analyzer system close to the field polarization direction (0 or $`\pi /2`$) is only limited by the optical precision of the analyzer used (which can be made arbitrarily small), and the noise in the detectors.

The cryptography system should be designed so that the uncertainties caused by the quantum noise of light in the measurement of the polarization angles is large. It can be shown that the number of bases $`N_\sigma `$ within a standard deviation of the measured angle is $`N_\sigma =M/(\pi |\alpha |)`$. Fig. 7 shows experimental results that confirm this dependence.

The effect of noise on signal recovery by an eavesdropper in an opaque attack can be simulated by sending repeatedly the same bit, but varying $`𝐊^{}`$, from A to B when B (playing Eve) does not apply the key to demodulate the signals. In the following sequence $`S`$, the first number in a brace is the basis set by Alice (selected from $`M=200`$ possible positions with $`N_\sigma 6`$, and correctly recovered by Bob with the use of the key) while the second number is the basis extracted by Eve through a single measurement of $`n_H`$ and $`n_V`$. $`S=\{110,117\},\{84,78\},\{108,99\},\{90,91\},\{100,107\},\{102,97\},\{84,84\},\{110,105\},\{110,111\},\{114,105\},\{82,86\},\{100,95\},\{92,72\},\{108,108\},\{102,90\},\{108,97\},\{96,93\},\{110,103\},\{112,121\},\{86,86\},\{102,100\},\{88,91\},\{102,94\},\{106,98\},\{118,135\}.`$  Clearly, Eve makes a large number of errors in determining the transmitted bases.

In conclusion, we have demonstrated that under individual attacks Yuen’s encryption protocol is secure with an adequate number of bases $`M`$ and without the need for intrusion detection. Key expansion is also possible under this scheme, due to the better observation available to Bob with knowledge of the key. Significantly, the encryption system allows for signal amplification as long as the security is guaranteed at the source for the following reason. While Eve has to resolve $`M`$ levels to tell the bit or $`𝐊^{}`$ without knowing the key, Bob has to resolve only two levels with the key. Thus, amplifier noise would hinder Eve’s attack while it will not disrupt Bob’s decryption. Furthermore, there is no need to decrypt/re-encrypt at the nodes of a properly designed communication line with a moderate number of amplifiers.

Acknowledgements – This work is supported by the DARPA grant F30602-01-2-0528. We thank Adam Rybaltovski for some help in the early stages of this work and Paul Voss for advice on data acquisition programming. E-mail: barbosa@northwestern.edu.