format: text

text:
# Astrometric Microlensing of Quasars

## 1 Introduction

Gravitational lensing acts on quasars in a number of ways. The best known effect is strong lensing which produces multiply imaged quasars. Only a small fraction of all quasars (about 1 in 500) are multiply imaged with galaxies acting as the main lenses (for an updated list, see the CASTLES web page http://cfa-www.harvard.edu/glensdata/). Most of these quasars – and possibly some “single” quasars as well – are also affected by microlensing: the coherent effects of stellar mass objects in the lensing galaxy. Microlensing is well known for affecting the apparent brightness of the quasar, which changes due to the relative motion between lens, source and observer (see, e.g., Paczyński 1986, Wambsganss, Paczyński & Schneider 1990, Lewis et al. 1998, Wyithe et al. 2000a, Wambsganss 2001, Woźniak et al. 2000a,b). However, stellar microlensing has yet another effect, which was first explored by Lewis and Ibata (1998): when the quasar crosses a (micro-)caustic, a new pair of microimages is created or destroyed. The relatively sudden (dis-)appearance of a highly magnified image pair does not just produce a strong fluctuation in apparent brightness, it can also shift the center-of-light of the quasar (the weighted sum of all the microimages) by tens of microarcseconds. Such a positional change is in principle detectable with current VLBI techniques in the radio regime, and will be measurable in the infrared and optical wavebands as well with (near) future technology (Keck-Interferometry, VLTI, SIM, GAIA).

A number of studies have explored astrometric microlensing in a different regime, namely stars in the Milky Way (or dark matter objects in its halo) acting on background stars either in the Magellanic Clouds or in the bulge of the Milky Way (e.g., Miyamoto & Yoshii 1995; Miralda-Escude 1996; Mao & Witt 1998; Boden, Shao & van Buren 1998; Goldberg & Woźniak 1998; Paczyński 1998; Han, Chun & Chang 1999; Han & Kim 1999; Han & Jeong 1999; Safizadeh, Dalal, Griest 1999; Dominik & Sahu 2000; Gould & Han 2000; Han & Kim 2000; Salim & Gould 2000; Delplancke et al. 2001; Belokurov & Evans 2002; Dalal & Lane 2003). The effects of Milky Way stars on background quasars were studied by Hosokawa, Ohnishi & Fukushima (1997), Sazhin et al. (1998) and Honma & Kurayama (2002). Williams & Saha (1995) had discussed large image shifts produced by substructure in the lensing galaxy. However, not much work has yet been done on cosmological astrometric microlensing – stars in lensing galaxies acting on even more distant quasars – beyond Lewis & Ibata (1998), with the recent exception of Salata & Zhdanov (2003). The present paper aims to explore this further.

We study quantitatively eight different cases with various values of the surface mass density, with and without external shear. We present example microlensing lightcurves with the corresponding center-of-light shifts and we investigate the correlation between high-magnification photometric events and large-offset astrometric events.

As first pointed out by Lewis & Ibata (1998), large offset events are typically correlated with high-magnification events, whereas the inverse is much less true. The reason for this is that the location of the newly appearing bright image pair during a caustic crossing is unrelated to the previous center-of-light of the microimages. In rare cases, the new image pair may appear at or close to the center-of-light of the pre-existing microimages. Such a situation would correspond to a large change in brightness with very little change in position. In most cases, however, the new bright image pair will appear at a location which is unrelated to the previous center-of-light, and hence produce a sudden jump. Since the positional offset is preferentially perpendicular to the direction of the external shear, the shear direction can be inferred this way.

After introducing the microlensing length and time scales (Section 2), we describe our simulations (Section 3). In Section 4, we illustrate the effect of astrometric and photometric fluctuations and present statistical correlations between positional offsets, magnitude fluctuations and time intervals between the measurements for four microlensing situations with different surface mass densities ($`\kappa =0.2,0.4,0.6,0.8`$), with and without external shear ($`\gamma =\kappa `$ or 0). Finally we discuss the possibilities of real detections of this phenomenon in the near future.

## 2 Microlensing basics: length, time and angular scales

### 2.1 Standard mass, length and time scales

The lensing effects of cosmologically distant compact objects in the mass range $`10^3M/M_{}10^3`$ on background sources is usually called “cosmological microlensing”. The source is typically a distant quasar, but in principle other objects can be microlensed as well, e.g. high-redshift supernovae (e.g., Rauch 1991) or gamma-ray sources/bursters (Torres et al. 2003; Koopmans & Wambsganss 2001). The only condition is that the source size be comparable to or smaller than the Einstein radius of the foreground lens.

The microlenses can be ordinary stars, brown dwarfs, planets, black holes, molecular clouds, or other compact mass concentrations. In most practical cases, the microlenses are part of a galaxy which acts as the main (macro-)lens. However, microlenses could also be located in, say, clusters of galaxies (Tadros, Warren & Hewett 1998, Totani 2003) or they could even be imagined “floating” freely and filling intergalactic space (Hawkins 1996, Hawkins & Taylor 1997).

The relevant length scale for microlensing is the Einstein radius in the source plane:

$$R_E=\sqrt{\frac{4GM}{c^2}\frac{D_SD_{LS}}{D_L}}3.2\times 10^{16}\sqrt{\frac{M}{M_{}}}\mathrm{h}_{75}^{0.5}\mathrm{cm}$$
(1)
where “typical” lens and source redshifts of $`z_L=0.5`$ and $`z_S=2.0`$ were assumed for the expression on the right hand side ($`G`$ and $`c`$ are the gravitational constant and the velocity of light, respectively; $`M`$ is the mass of the lens, $`D_L`$, $`D_S`$, and $`D_{LS}`$ are the angular diameter distances between observer – lens, observer – source, and lens – source, respectively, and a concordance cosmological model is assumed with $`\mathrm{\Omega }_{\mathrm{tot}}=1`$, $`\mathrm{\Omega }_{\mathrm{matter}}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$).

This length scale translates into an angular scale of:

$$\theta _E=\frac{R_E}{D_S}2.2\times 10^6\sqrt{\frac{M}{M_{}}}\mathrm{h}_{75}^{0.5}\mathrm{arcsec}.$$
(2)
It is obvious that image splittings on such small angular scales can not be observed directly. What makes microlensing observable in the first place is the fact that observer, lens(es) and source move relative to each other. Due to this relative motion, the microimage configuration changes with time, and so does the total magnification, i.e. the combined fluxes of all the microimages which make up the macro-image. This change in magnification over time (lightcurve) can be measured: microlensing is a dynamical phenomenon.

There are two time scales involved: the standard lensing time scale $`t_E`$ is the time it takes the lens (or the source) to cross a length equivalent to the Einstein radius, i.e.

$$t_E=(1+z_L)\frac{R_E}{v_{}}25\sqrt{\frac{M}{M_{}}}v_{600}^1\mathrm{h}_{75}^{0.5}\mathrm{years},$$
(3)
where the same typical assumptions are made as above, and the relative transverse velocity111For simplicity we assume here that all the transverse motion is done by the lensing galaxy $`v_{600}`$ is parametrized in units of $`v_{}=600`$km/sec and for a Hubble constant $`H_0=75h_{75}`$ km/sec/Mpc. This time scale results in relatively large values. However, microlensing fluctuations are expected (and observed, see Woźniak et al. 2000a,b) on much shorter time scales. The magnification distribution is highly non-linear with sharp caustic lines separating regions of low and high magnification. Often, the density of caustics is quite high, so that a source may encounter half a dozen or more caustic lines within the length of one Einstein radius. When a source crosses a caustic, we can expect a large change in magnification (and correspondingly in position) within the time it takes the source to cross its own radius. So the relevant time scale is the “crossing time”:

$`t_{cross}`$ $`=`$ $`(1+z_L){\displaystyle \frac{R_{source}}{v_{}(D_S/D_L)}}`$ (4)
$``$ $`0.58R_{15}v_{600}^1\mathrm{h}_{75}^{0.5}\mathrm{years}`$
$``$ $`210R_{15}v_{600}^1\mathrm{h}_{75}^{0.5}\mathrm{days}.`$
Here the quasar size $`R_{15}`$ is parametrized in units of $`10^{15}`$cm.

### 2.2 The special case of Q2237+0305

The quadruple quasar Q2237+0305 (Huchra et al. 1985; Irwin et al. 1989; Wambsganss et al. 1990; Wyithe et al. 2000a,b; Woźniak et al 2000a,b) is a very special and favorable case and of particular interest to microlensing studies. It was the first system in which microlensing was discovered (Irwin et al. 1989). Subsequently it received a lot of attention, both observational (Corrigan et al. 1991, Ostensen et al. 1996; Woźniak et al 2000a,b) and theoretical (Wambsganss et al. 1990; Wyithe et al. 2000a,b; Yonehara 2001). Due to the fact that the lensing galaxy is so close ($`z_G=0.039`$, Huchra et al. 1985), the physical and angular Einstein radii are considerably different from the standard case treated above:

$$r_{E,Q2237+0305}1.6\times 10^{16}\sqrt{\frac{M}{M_{}}}\mathrm{h}_{75}^{0.5}\mathrm{cm},$$
(5)
$$\theta _{E,Q2237+0305}7.3\times 10^6\sqrt{\frac{M}{M_{}}}\mathrm{h}_{75}^{0.5}\mathrm{arcsec},$$
(6)
The resulting time scales (Einstein time and crossing time) are much shorter than in almost all other multiple quasars:

$$t_{E,Q2237+0305}8.7\sqrt{\frac{M}{M_{}}}v_{600}^1\mathrm{h}_{75}^{0.5}\mathrm{years},$$
(7)
$`t_{cross,Q2237+0305}18R_{15}v_{600}^1\mathrm{h}_{75}^{0.5}\mathrm{days}.`$ (8)
For that reason, this quadruple system is ideally suited for microlensing studies. The length, time and angular scales for the “typical case” as well as for the special case of Q2237+0305 are summarized in Table 1.

## 3 The Simulations

In order to explore astrometric microlensing for a variety of realistic scenarios, we consider eight different cases, with the following values of the dimensionless surface mass density: $`\kappa =0.2,0.4,0.6`$ and 0.8. Each one of these we treated both without external shear ($`\gamma =0`$) and with external shear equal to the surface mass density ($`\gamma =\kappa `$), corresponding to an isothermal sphere model for the lensing galaxy. The shear always acts along the X-axis in the simulations. This means that the caustics are expanded in the X-direction and compressed in the Y-direction, as can be seen for example in the top left panel of Fig. 3. Horizontal bands appear on large scales (many Einstein radii). The effect is not always obvious on small scales however (sub-Einstein radii, bottom left panel of Fig. 3). Inversely, if we follow the light rays forward from source to observer (as opposed to backwards as in the simulations, which is explained below), the micro-image configuration of a source will appear compressed in the X-direction and expanded in the Y-direction (cf. the right hand panels of Fig. 3).

For each of these eight cases, we produced a magnification pattern in the source plane with a side length of $`L_A=20R_E`$, sampled on a $`1000^2`$ pixel grid, i.e. $`1R_E`$ is covered by 50 pixels, or $`1\mathrm{pix}=0.02R_E`$. This is our starting configuration A. The minimum source size we can consider is given by the pixel scale. In order to explore smaller source sizes as well (at the cost of smaller magnification patterns which may not be representative), we zoom in into the central part of the first configuration. We do higher resolution simulations (all on $`1000^2`$ pixel grids) in three steps of factors of two (configurations B, C and D, respectively), resulting in magnification patterns with side lengths of $`L_B=10R_E`$, $`L_C=5R_E`$, and $`L_D=2.5R_E`$, respectively. For the highest resolution simulation (configuration D), this means $`R_E=400`$ pixels, or 1 pix $`=0.0025R_E`$.

We used a modified version of the ray shooting code described in Wambsganss (1990, 1999). In the original version, rays are followed backwards from the observer through the lens plane (where all the deflections are determined) to the source plane (where the light rays are “collected” in pixels). The density of rays then indicates the magnification as a function of position in the source plane, often displayed as color coded magnification patterns. Lightcurves can be obtained by convolving a given source profile with this two-dimensional magnification map. However, all information about where the rays originated from in the lens/image plane is lost in this original algorithm. For the exploration of the centroid shift of the collection of microimages, it is exactly this information that is required. Hence we modified the code to record the positions of all the individual microimages brighter than a given magnification threshold. To do this, we defined a second, larger regular grid of $`4000^2`$ rays, which covers in the image plane a region of four times the angular side length of the magnification pattern (i.e. $`(80\theta _E)^2`$ for configuration A). For a set of $`1.6\times 10^7`$ test rays, we keep track of the positions in the image plane as well as in the source plane. In this way, we can find all the positions in the image plane that are mapped onto a certain area in the source plane, and hence identify all the microimages corresponding to a particular source position and size.

In order to get information on both the magnification and the microimage locations of a finite source at a given position, we convolved the magnification pattern with a luminosity profile. We used circular sources with Gaussian widths $`\sigma =`$ 2, 4, 8 and 16 pixels, corresponding to physical sizes between $`0.04R_E`$ and $`0.32R_E`$ in the starting configuration A, or $`1.5\times 10^{15}`$cm to $`1.2\times 10^{16}`$cm in the typical case described in Section 2.1. We did the same with the higher resolution magnification patterns, configurations B, C and D. However, it turned out that the high resolution cases, although well suited for studying individual caustic crossings for small sources, are not quite large enough for statistical investigations. For this reason, we restricted ourselves to cases A and B for the statistical evaluations below. The numerical values corresponding to the different source sizes we used are tabulated in Table 2.

Using these simulations, we determined the positions of the individual microimages, the center-of-light and the total magnification of the macro-image corresponding to a particular source position and profile. In Figs. 1 and 2, a microimage situation is shown for each of the eight cases considered: $`\kappa =0.2,0.4,0.6`$ and 0.8 with $`\gamma =0`$ (Fig. 1), and with $`\gamma =\kappa `$ acting along the X-axis direction (Fig. 2). In the left columns, a small part of the magnification pattern is shown, with the source position and profile superimposed. The panels on the right hand side show the particular microimage configurations, with the light centroid indicated by a plus sign (cf. Paczyński 1986).

In Figure 3, the method is illustrated for $`\kappa =0.8`$ and $`\gamma =0.8`$. The top left panel shows the magnification pattern in the source (quasar) plane with side length $`L=10\theta _E`$ (B configuration). The grey scale indicates regions of different magnifications in the source plane: the lighter the grey, the higher magnification. The large square indicates the region in the source plane for which all source positions were evaluated. The top right panel shows the corresponding microimage configuration in the image plane: the white regions are the parts of the image plane where light bundles appear which originated from within the square on the left (side length is 40 $`\theta _E`$). In other words, it shows what a large square shaped source would look like to the observer. The bottom panels show the same for the smaller square region indicated in the top left panel, zoomed eight times. The many isolated light patches indicate that microimages are spread over a very large area in the image plane.

Figure 4 shows a quasar microlensing scenario with microlensing parameters $`\kappa =0.6`$ and $`\gamma =0.6`$, and side length 10 $`\theta _E`$ (B configuration). The straight vertical white line marks the track of the quasar motion relative to the magnification pattern; the length of the path is 2.0 $`R_E`$. For this particular track (followed from the lowest part upwards), Figure 5 shows from top to bottom: the X- and Y-coordinates of the quasar relative to the starting position ($`\mathrm{\Delta }\theta _X`$ and $`\mathrm{\Delta }\theta _Y`$) as a function of time; the absolute value of the positional shift ($`|\mathrm{\Delta }\theta |`$) relative to the starting position as a function of time; and the corresponding light curve ($`\mathrm{\Delta }m`$) of the quasar. The solid and dotted lines correspond to two different values of the source size: $`\sigma =0.04R_E`$ (4 pixels in B configuration) and $`\sigma =0.16R_E`$ (16 pixels in B configuration), respectively. The track in Fig. 4 starts in a region of low magnification which is taken as the zero point of the magnitude scale on the lowest panel in 5.

The two panels in Figure 6 represent the centroid shift for the two different source sizes. The whole track has a length of 2.0 $`t_E`$ ($``$ 50 and 15 years for the ’typical’ lensing case and the Q2237+0305 case, respectively). The labels $`t_0`$, $`t_1`$ and $`t_2`$ correspond to the starting, middle and final positions, respectively. Snapshots of the ensemble of microimages at times $`t_0`$, $`t_1`$ and $`t_2`$ are shown in Figure 7 assuming $`\sigma =0.16R_E`$. All microimages are plotted up to a radius of 3 $`\sigma `$, without accounting for the declining source profile. On average, the brightness of the macroimage declines as the fourth power of distance to the center of light (Katz, Balbus, Paczyński 1986). On our plots, due to the finite resolution, some of the distant and faint microimages appear bigger than they are. This figure simply aims at illustrating the very large spatial spread of the microimages, covering roughly $`30\theta _E\times 40\theta _E`$. The lower panels show enlargements of the central parts with the path of the center-of-light superimposed.

## 4 Results and Discussion

### 4.1 Microlens-induced positional offsets and corresponding magnitude changes

In order to evaluate the correlations between the magnitude changes and the positional changes of a microlensed quasar, we simulated pairs of measurements separated by time intervals of 0.02 $`t_E`$ (Figures 8, 10) and 0.16 $`t_E`$ (Figures 9, 11), respectively. These values correspond to about half a year and four years in the “typical” case, and to two months and 1.5 years in the case of Q2237+0305 (cf. Table 1). We placed the source at a random position in the magnification pattern and chose a second position at a distance of $`0.02R_E`$ or $`0.16R_E`$ either parallel or perpendicular to the action of the external shear (X-axis). We determined differences in the magnifications $`|\mathrm{\Delta }m|`$ and the center-of-light positions $`|\mathrm{\Delta }\theta |`$ between these two source positions. Each pair of measurements ($`|\mathrm{\Delta }m|`$, $`|\mathrm{\Delta }\theta |`$) is represented as a point in the various panels of Figures 8 to 11. In each panel, the offset in position $`|\mathrm{\Delta }\theta |`$ is shown against the corresponding offset in magnitude $`|\mathrm{\Delta }m|`$ for about 20,000 such measurement pairs. In the left columns of all four figures, a small Gaussian source size of $`\sigma =0.02R_E`$ is assumed, in the right column the source size is $`\sigma =0.16R_E`$. In each panel, the thick (lower) line indicates the median of $`|\mathrm{\Delta }\theta |`$ as a function of $`|\mathrm{\Delta }m|`$, and the thin (upper) line shows the 95th-percentile: 5% of all simulations for a given $`|\mathrm{\Delta }m|`$ would result in offsets $`|\mathrm{\Delta }\theta |`$ which are above this thin line.

In Figure 8, the four lensing situations without external shear are considered: $`\kappa =0.2,0.4,0.6,0.8`$ (from top to bottom) with a small time step $`\mathrm{\Delta }t=0.02t_E`$. In all panels, the majority of the points cluster near the origin. This is easily understandable: the position of the source relative to the caustics has not changed by much during the short time interval, so the changes both in magnification and in the centroid position tend to be small. For small sources (left panels), however, $`|\mathrm{\Delta }m|`$ and $`|\mathrm{\Delta }\theta |`$ can reach relatively large values with the median lines indicating an almost linear statistical relation. The slope of these median lines slightly increases with increasing surface mass density $`\kappa `$. The time step corresponds roughly to the crossing time: $`\mathrm{\Delta }t=0.02t_Et_{\mathrm{cross}}`$, and there are indeed cases with easily measurable magnitude fluctuations ($`|\mathrm{\Delta }m|>0.2`$ mag), which result in center-of-light offsets of about 0.5 $`\theta _E`$. In the right panel, there are no significant changes in either magnification or center-of-light position, because the source has moved only a fraction of its own diameter ($`\mathrm{\Delta }t<t_{\mathrm{cross}}`$).

In Figure 9, the same is shown for a larger time step $`\mathrm{\Delta }t=0.16t_E`$. Many more points are now spread towards larger offsets and larger magnification changes. For the small source, median values of $`|\mathrm{\Delta }\theta |1\theta _E`$ are reached in the highest surface mass density cases. In fact, the 95th-percentile line for the $`\kappa =0.8`$ case indicates that for magnitude changes $`|\mathrm{\Delta }m|>0.4`$mag, 5% of the all cases result in center-of-light offsets larger than 4$`\theta _E`$. For the large source (right column) the expected offsets are still quite moderate, with median values of about 0.5 $`\theta _E`$. This is not too surprising, because the time interval corresponds to just about the crossing time for the source.

Figures 10 and 11 contain the same diagrams for the cases with external shear, $`\kappa =\gamma =0.2,0.4,0.6,0.8`$ (from top to bottom). Larger offsets are reached here than in the corresponding scenarios without shear. Especially the two middle rows with $`\kappa =\gamma =0.4`$ and 0.6 – which best represent the typical values of convergence and shear of a multiply imaged quasar – produce values of positional offset a factor of two higher than the corresponding cases without shear: For short time steps and small sources, medians of $`|\mathrm{\Delta }\theta |1\theta _E`$ and 95th-percentiles of $`|\mathrm{\Delta }\theta |2\theta _E`$ (small source, left column) are reached. For smaller values of the surface mass density, the caustic density and hence the number of microimages is not large enough to produce big positional offsets, whereas for higher values of $`\kappa `$, the density of caustics is so high that an additional microimage pair only produces a small fluctuation.

Results are even more dramatic when assuming a larger time step ($`\mathrm{\Delta }t=0.16t_E`$), as displayed in Figure 11: median values of $`|\mathrm{\Delta }\theta |`$ between 0.5$`\theta _E`$ and 3$`\theta _E`$ for the small source size (left column), and significant median values of $`|\mathrm{\Delta }\theta |0.5\theta _E`$ to $`1\theta _E`$ even for the large source (right column).

### 4.2 Microlens-induced positional offset as a function of time (with magnitude thresholds)

We evaluated the shifts in center-of-light positions for time intervals intermediate between $`0.02\mathrm{\Delta }t/t_E0.16`$ assuming a threshold for the magnitude fluctuations. In Figures 12 and 13, the median and 95th-percentile offsets are shown as a function of increasing time step $`\mathrm{\Delta }t`$, for $`|\mathrm{\Delta }m|0.2`$ mag (left panels) and $`|\mathrm{\Delta }m|0.5`$ mag (right panels). This was done with the following idea in mind: since very accurate measurements of quasar image positions are “expensive”, it is unlikely that all lensed quasar candidates can be astrometrically monitored. However, photometric monitoring is comparably cheaper. So ideally, one could determine the positions of the most promising multiple quasars with high accuracy once, then monitor them photometrically, and whenever a large microlens-induced magnitude change has been detected, a second astrometric measurement should be performed.

In each panel of Figure 12, the two sets of curves show the median (thick lines) and the 95th-percentile (thin lines) of the $`|\mathrm{\Delta }\theta |`$ distribution as a function of $`\mathrm{\Delta }t`$ for the small source ($`\sigma =0.02R_E`$, solid) and the large source ($`\sigma =0.16R_E`$, dotted). All curves show basically the same behaviour: the $`|\mathrm{\Delta }\theta |`$ values first increase with increasing $`\mathrm{\Delta }t`$, then flatten out. This behaviour shows that many or most jumps in magnitude and position are dominated by one fold-caustic crossing. The time scale is dominated by $`t_{\mathrm{cross}}`$. For a slightly larger time interval, the offset does not increase significantly any more. Only for much larger $`\mathrm{\Delta }t`$’s allowing for additional caustic crossings, would another increase in $`|\mathrm{\Delta }\theta |`$ be expected. Offsets of more than 100 microarcseconds could indeed be reached this way, but the characteristic time scale would be depressingly large (many decades!).

The qualitative behaviour of the cases with external shear displayed in Figure 13 is similar to those without shear. However, the expected offset values are significantly higher here, as already seen in Figs. 10 and 11: again, the cases $`\kappa =\gamma =0.4`$ and 0.6 appear most promising (second/third row), with median values of between 1 $`\theta _E`$ and 2 $`\theta _E`$, and 95th-percentiles of 4 $`\theta _E`$ or higher. These values translate into about 7 to 15 microarcseconds (median) and $`>`$ 30 microarcseconds (95th percentile) when applied to the quadruple quasar Q2237+0305 (cf. Table 1). In fact, brightness fluctuations of more than one magnitude have been measured in Q2237+0305 on time scales of a few months (Woźniak 2000a,b), so these events do occur and seem not to be very rare.

### 4.3 Discussion and comparison with previous work

We explored astrometric microlensing effects on short timescales (of order months to years) for a set of parameter values $`\kappa `$ and $`\gamma `$. We find that the relevant time scale for measuring relatively large jumps of occasionally many tens of microarcseconds can be as short as a few months. Such sudden changes of position – produced by caustic crossings – are statistically related to fluctuations in the apparent brightness of the quasar. Therefore, a good strategy for detecting this centroid shift would be to measure the positions of the most promising lensed quasars very accurately once, then to monitor the quasars in the optical, and – when a photometric microlensing event is detected – to perform one or a few more accurate astrometric measurements.

We found that the effect is most pronounced for values of the surface mass density and shear $`\kappa =\gamma =0.4`$ or 0.6. These parameters happen to be applicable to many of the lensed quasar images. The most favorable case is the quadruple lens Q2237+0305: because of the closeness of the lensing galaxy, the time scale is relatively short (cf. also Woźniak 2000a,b). We also investigated the positional shift of the image as a function of source size: whereas a quasar with a typical size of $`0.16\theta _E`$ produces median offsets of order $`\theta _E`$ (and 95th-percentiles of about $`2\theta _E`$), smaller sources ($`0.02\theta _E`$) reach median values of $`2.5\theta _E`$ and 95th-percentiles larger than $`5\theta _E`$ (where in the case of Q2237+0305, $`\theta _E7.3\times \sqrt{M/M_{}}\mathrm{h}_{75}^{0.5}\mu \mathrm{arcsec}`$, cf. Table 1).

Lewis and Ibata (1998) had investigated the astrometric microlensing effect specifically on Q2237+0305. They had found that substantial image shifts of $`100\mu `$arcsec are possible within months. We can confirm this in some rare cases. For typical caustic crossings with photometric fluctuations of about 0.5 mag, we find values between 20 and 40 $`\mu `$arcsec. Salata & Zhdanov (2003) have looked into the question of rms fluctuations of the quasar position; however, they only considered large sources ($`0.5R_E`$) and cases with $`\kappa +\gamma <0.8`$.

### 4.4 Prospects for observations

Ground-based differential astrometry in the near infrared is able to achieve measurement uncertainties of better than 10 microarcsec, as was reported very recently from the Palomar Testbed Interferometer (Lane & Muterspaugh 2003). This is only feasible for bright objects so far, but it is a very exciting result which opens up promising opportunities for the coming years. Another instrument promising extremely high astrometric precision in the near future (with planet detection as one of its main scientific drivers) is the PRIMA instrument at the ESO VLTI. It should become efficient in 2004 (Paresce et al. 2003; see also Delplancke et al. 2001). The current goal at ESO is to achieve 50 $`\mu `$arcsec accuracy with PRIMA in the $`H`$ and $`K`$ bands in 2005-2008 and 10 $`\mu `$arcsec accuracy in 2008-2010 (Th. Henning, private communication).

A number of space-based astrometric projects are also underway: The Space Interferometric Mission (SIM) is a five year mission scheduled for launch in 2009 (http://sim.jpl.nasa.gov). The mission’s goal is to reach an astrometric accuracy of about 1 microarcsecond for a predefined grid of objects brighter than 13 mag in the visible, which is quasi-inertially tied to a set of distant QSOs. For slightly fainter objects (the four images of the lensed quasar Q2237+0305 have magnitudes $``$ 17), SIM is expected to yield 4 microarcsecond absolute positions (Unwin et al. 2002). In order to detect center-of-light offsets for multiple quasars, in fact only relative astrometry between the quasar images is required. SIM will also be able to measure such positional shifts as a function of color, which will give us hints on the physical structure of the continuum emission region: presumably the cooler/redder part is more extended than the hotter/bluer part, which means that we expect larger changes in the center-of-light at shorter wavelengths. The GAIA satellite is an ESA mission currently scheduled for launch in June 2010 (see http://sci.esa.int/gaia; Perryman et al. 2001, Perryman 2002). With a nominal precision of a few microarcseconds for bright objects (about 10 microarcseconds for 15th mag objects), it will measure accurate positions of 500 000 quasars. So GAIA is expected to be an extremely useful instrument for astrometric microlensing purposes. It will provide many positional shifts of quasar images, along with their lightcurves in many filters.

## 5 Summary and Conclusions

We analyzed the shifts in the center-of-light positions of gravitationaly lensed quasar images and the corresponding flux variability due to the microlensing effect of stars in the lensing galaxy. We found the following results:

1. The center-of-light of multiply imaged quasars is expected to vary on time scales of months to years, due to astrometric microlensing effects.
2. This effect depends on the size of the lensed quasar, which we modelled with a Gaussian profile: it is larger for smaller sources. This means that measuring centroid shifts will allow us to constrain the size of the quasar.
3. We studied eight cases, four without external shear ($`\kappa =0.2,0.4,0.6,0.8`$; $`\gamma =0.0`$) and four with shear ($`\kappa =\gamma =0.2,0.4,0.6,0.8`$). The strongest effects are seen in the cases $`\kappa =\gamma =0.4`$ and 0.6, which in fact are very typical parameter values for many lensed quasar images.
4. The effect of center-of-light shift of quasars is not limited to multiply imaged quasars: single quasars may occasionally be microlensed too (as was suggested by Hawkins 1996). The offsets should be relatively small, however, due to the fact that the values of $`\kappa `$ and $`\gamma `$ will be lower for single quasars than for multiply imaged ones. Our simulation with $`\kappa =0.2`$ comes closest to describing these cases with very low convergence.
5. The centroid shift effect is statistically strongly correlated with changes in the apparent brightness of the quasar. Both occur during the crossing of a (micro-)caustic: photometric fluctuations $`0.5`$ mag produce offsets $`2\theta _E`$ in 50% of the cases, and $`5\theta _E`$ in 5% of the cases. These values correspond to about 15 and 35 microarcseconds, respectively, for the lensed quasar Q2237+0305.
6. The best strategy for measuring centroid shifts is therefore to monitor the apparent brightness of multiple quasars with ground based telescopes and to determine the relative positions of the various quasar images occasionally, ideally before, during and after a high magnification event.
7. In the cases with external shear, the center-of-light shift depends strongly on the direction of the shear, which can hence be determined in a statistical sense.

With the next generation of astrometric instruments providing an accuracy of order 10 microarcseconds, the astrometric microlensing effect of stars acting on background quasars will become detectable.


text:
# Upper limits on the strength of periodic gravitational waves from PSR J1939+2134

## 1 Introduction

The LIGO and GEO gravitational wave detectors performed their first science run, denoted S1, from 23 August to 9 September 2002 (Abbott et al., 2003a, for the LIGO Scientific Collaboration). Although these instruments were still to reach their design sensitivities, their performances during this period were sufficiently good to justify a serious test of our search algorithms on real interferometer data. Code to search for continuous gravitational waves from rapidly rotating neutron stars has been under development within the LIGO Scientific Collaboration (LSC) since the mid-to-late 1990s, and for S1 the LSC Pulsar Upper Limits Group (PULG) used two independent search methods to set upper limits on signals from the millisecond pulsar PSR J1939+2134. This object has the shortest spin period of any known pulsar and falls in a relatively quiet part of the detector bands. Its short period also enables us to place a relatively tight upper limit on the equatorial ellipticity of the neutron star. Here we present a brief outline of the results from this exercise. A more detailed presentation of the methods and discussion of the results can be found in Abbott et al., 2003b (for the LIGO Scientific Collaboration).

## 2 Detector sensitivities

In setting an upper limit on the gravitational wave amplitude from a pulsar we are saying something about the pulsar’s physical parameters as constrained by instrumental performance rather than simply making a statement about the strain sensitivity of our instrument. The strength of any continuous gravitational wave signal from PSR J1939+2134 seen by a detecting interferometer will depend not only on the pulsar’s distance, quadrupole moment and spin period, but also on the relative orientation between its axis of spin and the arms of the interferometer. The diurnal motion of the pulsar through the antenna pattern will also modulate the strength of the signal, so there is no simple relationship between the strain sensitivity curve of an interferometer and the upper limit that interferometer can place on the parameters of a particular pulsar. We can however define a *mean* performance, averaged over these unknowns, for a particular data set and detection method.

Figure 1 shows how our interferometers constrain monochromatic strain amplitudes at Earth as a function of signal frequency, averaged over the antenna pattern and relative orientation of the pulsar spin axis to the interferometer. The curves are derived from the ‘frequency domain method’ described below, and represent detection with a 1% false alarm rate and a 10% false dismissal rate using the actual observation times and data from S1. The lower solid curves show the design sensitivities of the interferometers for a 1 y observation period. The dotted lines indicate the strain amplitudes corresponding to hypothetical neutron stars at 8.5 kpc with equatorial ellipticities of $`10^3`$, $`10^4`$ and $`10^5`$. Finally the dots represent the upper limit on strain amplitude for a selection of pulsars. For this we have assumed a simple spindown model in which all the apparent loss in rotational kinetic energy of the neutron star goes into the emission of gravitational waves. It is clear from this figure that no detection was expected from the S1 data. Our primary aims in this exercise were to test the data reduction and analysis methods in the presence of realistic interferometer noise and set an upper limit on the strengths of such waves.

## 3 Search methods

Although the location, period and period derivative of PSR J1939+2134 are known to high precision the inclination of the neutron star’s spin axis to the line of sight ($`\iota `$) and the phase ($`\varphi `$), polarization ($`\psi `$) and, of course, amplitude ($`h_0`$) of the gravitational wave are unknown. All these factors affect the form of the signal and were searched over in S1 using two independent methods.

The ‘frequency domain method’ uses Fourier transform techniques to set a frequentist upper limit on the strain amplitude from the pulsar. The input time series is broken down into 60 s-duration short Fourier transforms (SFTs) which are high pass filtered at 100 Hz, Tukey windowed in the time domain and amplitude calibrated once per minute. These SFTs are then optimally combined to take account of the Doppler evolution and spindown of the signal using a Dirichlet kernel method (Jaranowski et al., 1998) and used to compute a detection statistic, denoted $``$, which is the likelihood ratio of the data, maximized over $`\iota `$, $`\varphi `$ and $`\psi `$. Monte Carlo injections into the data stream are then used to calculate the probability distribution function (pdf) of $``$. The upper limit is defined such that had a signal at the upper limit level or greater been present in an ensemble of similar experiments to the one performed, 95% of the experiments would yield values of $``$ that exceed the one seen in the real experiment. The computing requirement for this Monte Carlo work is extensive, and the searches ran off-line on the Medusa cluster at UWM, (296 single 1 GHz CPU nodes and 58 TB of disk space), and the Merlin cluster at AEI Potsdam (180 dual 1.4 GHz CPU nodes and 36 TB of disk space).

The ‘time domain method’ is an algorithm specifically developed to search for neutron stars with a complex, but known, rotational phase evolution. It uses heterodyne methods to reduce and filter the data to a rate of one complex value every 60 s, estimating the noise variance over the same period to take account of non-stationarity in the data. The algorithm proceeds to carry out the search as a standard Bayesian modelling problem, determining a posterior pdf for the unknown source parameters $`h_0`$, $`\iota `$, $`\varphi `$ and $`\psi `$. The prior for $`h_0`$ is chosen to be uniform over $`[0,\mathrm{}]`$ and 0 for $`h_0<0`$, and the three other priors are chosen to be least informative about the orientation of the neutron star. The time domain Bayesian upper limit is set via the posterior pdf for $`h_0`$, marginalized over the nuisance parameters $`\iota `$, $`\varphi `$ and $`\psi `$. We define the 95% upper credible limit for $`h_0`$ as the value at which the cumulative probability of $`h_0`$ (from zero) is 0.95.

Note that there is nothing special about the choice of domain (frequency or time) for the frequentist and Bayesian analyses – the above choice was influenced by the desire to have two very different approaches to solving the same problem.

## 4 Results

As expected the S1 data provided no evidence for continuous wave emission from PSR J1939+2134 at twice its rotation frequency (as determined by radio observations). However, we were able to improve the best upper limit for such signals by about a factor of 100, based on earlier work by Hough et al. (1983) and Hereld (1983).

The results from the four interferometers that took part in S1 are summarized in Table 1, the quoted uncertainties being largely due to calibration. The two analysis methods are very different, and the meaning of ‘upper limit’ is not the same for both, however we would hope that the concept of an upper limit is sufficiently robust that the two interpretations would deliver numerically similar results, as indeed they do.

In principle the Bayesian analysis allows the data from all four interferometers to be combined to give an improved joint upper limit, although the relative phase calibration of the instruments was not sufficiently reliable to do this for S1. It is however interesting to note that an upper limit for $`h_0`$ can, within a suitable model, be interpreted as an upper limit on the equatorial ellipticity of the neutron star. If we take the distance to the pulsar to be 3.6 kpc, its radius to be 10 km and its mass to be $`1.4M_{}`$ a value of $`h_0<1.4\times 10^{22}`$ corresponds to an ellipticity of less than $`2.7\times 10^4`$. Such a large upper limit to the ellipticity could be rejected a priori, as it implies a much greater spindown rate (via gravitational luminosity) than is seen in this pulsar. However as detector sensitivities improve we can expect the above methodology to deliver constraints on neutron star parameters not otherwise achievable.

## 5 Acknowledgements

The authors gratefully acknowledge the support of the United States National Science Foundation for the construction and operation of the LIGO Laboratory and the Particle Physics and Astronomy Research Council of the United Kingdom, the Max-Planck-Society and the State of Niedersachsen/Germany for support of the construction and operation of the GEO 600 detector. The authors also gratefully acknowledge the support of the research by these agencies and by the Australian Research Council, the Natural Sciences and Engineering Research Council of Canada, the Council of Scientific and Industrial Research of India, the Department of Science and Technology of India, the Spanish Ministerio de Ciencia y Tecnologia, the John Simon Guggenheim Foundation, the David and Lucile Packard Foundation, the Research Corporation, and the Alfred P. Sloan Foundation.

## References


text:
# Cooper-pair resonances and subgap Coulomb blockade in a superconducting single-electron transistor

## Abstract

We have fabricated and measured superconducting single-electron transistors with Al leads and Nb islands. At bias voltages below the gap of Nb we observe clear signatures of resonant tunneling of Cooper pairs, and of Coulomb blockade of the subgap currents due to linewidth broadening of the energy levels in the superconducting density of states of Nb. The experimental results are in good agreement with numerical simulations.

The single-electron transistorsingle  (SET) and its superconducting version is one of the most versatile tools in mesoscopic physics. It has been used for extremely sensitive charge measurements,rfset  for the construction of Cooper pair pumps and other adiabatic devices with applications in metrology,pumps  and more recently for building up superconducting quantum bits.chargequbits

The $`IV`$ characteristics of superconducting SETs present the usual features of quasiparticle tunneling (at voltages above $`2\mathrm{\Delta }_{\mathrm{Nb}}+2\mathrm{\Delta }_{\mathrm{Al}}`$), Josephson-quasiparticle tunneling (at half of these values), and Josephson effect (around zero bias). These features have been thoroughly investigated by now by many groups and the physics of a charge transport at these bias voltages is well understood. However, at low bias voltages also other transport processes could become important and can alter the performance of Josephson-based devices. In this paper, we study two such processes appearing in our Nb-based SET: Resonant tunneling of Cooper pairs, and transport through states inside the gap of Nb (subgap currents).

We have fabricated Al/AlOx/Nb/AlOx/Al single electron transistors using a lithographic technique described elsewhere.us  Measurements were done using a small dilution refrigerator equipped with well-thermalized and electrically filtered measuring lines. The superconducting gaps obtained for Nb and Al ($`\mathrm{\Delta }_{\mathrm{Nb}}=1.4`$ mV, $`\mathrm{\Delta }_{\mathrm{Al}}=0.2`$ meV), and also the measured critical temperatures for Nb ($`T_{\mathrm{C},\mathrm{Nb}}8.08.5`$ K), show that the films are indeed of good quality.

At voltages below the gap of Nb, a series of gate-dependent resonance peaks appears in the $`IV`$ characteristics of the SET (Fig. 1). We interpret this as resonant tunneling of Cooper pairs, a transport phenomenon first predicted theoretically and later observed in Al symmetrically-biased superconducting SETs.haviland  Below we describe the same process for our Nb-island SETs under the asymmetric bias shown in Fig. 1. We consider a generic process in which a charge $`\delta q_1`$ tunnels through the left junction and a charge $`\delta q_2`$ tunnels through the second junction, both into the island. During the process, a charge $`\delta q=\delta q_1+\delta q_2`$ is transferred into the island and a charge $`\delta Q=\delta q_1\delta q_2`$ is transferred through the external circuit in the forward direction. The change in the electrostatic free energy (including work done by the sources) associated with this process is

$`E(\delta q,\delta Q)`$ $`=`$ $`{\displaystyle \frac{(\delta q)^2}{2C_\mathrm{\Sigma }}}+{\displaystyle \frac{\delta q}{C_\mathrm{\Sigma }}}(q_0C_\mathrm{g}V_\mathrm{g})`$ (1)
$`+{\displaystyle \frac{\delta q}{C_\mathrm{\Sigma }}}\left(C_\mathrm{g}{\displaystyle \frac{C_2C_1}{2}}\right)V{\displaystyle \frac{\delta Q}{2}}V,`$
where $`q_0`$ is the initial charge of the island, $`C_1`$ and $`C_2`$ are the capacitances of the left and right junctions, $`C_\mathrm{g}`$ is the gate capacitance, and $`C_\mathrm{\Sigma }=C_1+C_2+C_\mathrm{g}`$. Resonant Cooper pair tunneling in superconducting SETs occurs when no energy is required for processes resulting in the transport of $`m`$ Cooper pairs $`\delta Q=2me`$ through the external circuit and the creation of an excess of $`n`$ Cooper pairs $`\delta q=2ne`$ on the island, i.e., $`E(2ne,2me)=0`$.

The dominant processes are those involving $`m=\pm 1`$ and $`n=\pm 1`$, corresponding to a single Cooper pair tunneling through either one of the junctions.haviland  For these processes, the resonant condition describes two lattices with a cell size of $`(4e/C_\mathrm{\Sigma })\times (2e/C_\mathrm{g})`$ (in $`V\times V_\mathrm{g}`$ plane) corresponding to odd and even values of $`q_0/e`$, and displaced with respect to each other by $`2e/C_\mathrm{\Sigma }`$ along the $`V`$ axis and by $`e/C_\mathrm{g}`$ along the the $`V_\mathrm{g}`$ axis. During the time of an experiment, $`q_0`$ can either be fixed, in which case only one lattice pattern will appear, or it can fluctuate between odd and even values, in which case what will be measured is the overlap between the two lattices, resulting in a checkerboard lattice with half the periodicity. This last situation occurs indeed in all of our three samples. To check this, we first measured the gate modulation at voltages above the quasiparticle threshold of Nb, where the transport is of single-electron type, and determined $`e/C_\mathrm{g}`$. This corresponds to the size of the checkerboard pattern along $`V_\mathrm{g}`$ that we see in our measurements at lower bias as well, e.g., $`e/C_\mathrm{g}=12.3`$ mV ($`C_\mathrm{g}=12.9`$ aF) for the sample of Fig. 1. (For clarity, we will concentrate our discussion on that particular sample from now on. The other two samples yielded similar results.) The observation that the checkerboard pattern is $`1e`$-periodic allows us to determine $`e/C_\mathrm{\Sigma }=68.3`$ $`\mu `$V and thus the charging energy $`E_\mathrm{C}=e^2/2C_\mathrm{\Sigma }=34.1`$ $`\mu `$eV.joyez  The SET’s junction asymmetry is reflected in the different absolute values of the positive-slope lines (for which $`n`$ and $`m`$ have different signs) versus the negative-slope lines ($`m`$ and $`n`$ have the same sign), yielding $`C_1=1.03`$ fF and $`C_2=1.29`$ fF.

To check for consistency, we have determined $`E_\mathrm{C}`$ also by standard Coulomb-blockade thermometry cbt  with the superconductivity both in Al and Nb suppressed. We find $`E_\mathrm{C}=36`$ $`\mu `$eV, which is in very close agreement with the results above.

The second intriguing feature in the low-bias $`IV`$ of our Nb-island SETs, is the presence of relatively large subgap currents: from figures 1 and 2(a) one can see that the Cooper pair resonances appear as being built upon a current which increases with bias voltage (at voltages larger than $`2\mathrm{\Delta }_{\mathrm{Al}}+2\mathrm{\Delta }_{\mathrm{Nb}}`$, it will merge into the quasiparticle current). The origin of the subgap currents in superconductors is still a matter of intense theoretical and experimental investigations.review  Here, we further explore these currents by using a magnetic field to suppress the gap of Al; the transition to the normal-metal state in the leads is clearly indicated by the disappearance of the Josephson effect and of the Cooper pair resonances as seen in Fig. 2(a).

While Al is in the normal state, a very interesting feature appears at low bias voltages \[Fig. 2(b)\], where we observe a dip in the conductance. In measurements on single Al-Nb junctions \[again with the gap of Al suppressed, Fig. 2(b) inset\] this feature is not present.us  Instead, below 0.2 mV the $`IV`$ is approximately ohmic. This shows that the dip is due to the Coulomb blockade of the subgap current, and that in the first approximation one could attempt to fit the graph with the standard analytical expressions for Coulomb-blockade thermometry cbt  (CBT) \[Fig. 2(b), main graph, red dashed line\]. Fitting yields $`T98`$ mK, which agrees well with the value measured using a calibrated resistor thermally anchored to the mixing chamber, and $`E_\mathrm{C}13`$ $`\mu `$eV, which is in reasonable agreement with the value determined from Cooper pair resonances, considering that we are not in the limit $`k_\mathrm{B}TE_\mathrm{C}`$ and therefore the standard approximations of Coulomb blockade thermometry are not accurate, especially for determining $`E_\mathrm{C}`$.cbt

A better model should take into account that we are at low temperatures $`k_\mathrm{B}T<E_\mathrm{C}`$ and also describe the nonlinear increase in the subgap current at higher voltages. To develop such a model, we introduce a lifetime broadening $`\mathrm{\Gamma }`$ of the quasiparticle energies, resulting in a density of states mi

$$\rho (E,\mathrm{\Delta }_{\mathrm{Nb}})=\left|\mathrm{Re}\left(\frac{Ei\mathrm{\Gamma }}{\sqrt{(Ei\mathrm{\Gamma })^2\mathrm{\Delta }_{\mathrm{Nb}}^2}}\right)\right|.$$
(2)
This effect could be caused by the proximity to the metal-insulator transition, due to the granular structure of the Nb films (with our films being still on the metalic side, $`\mathrm{\Gamma }\mathrm{\Delta }_{\mathrm{Nb}}`$).mi  Another possibility is the opening of conduction channels in the junctions,joyez  as suggested by the relatively large value of the room-temperature conductance \[for the sample presented in Fig. 1, Fig. 2, and Fig. 3(a), this was 30 $`\mu `$S\]. Irrespective to the microscopic origin, this density of states accounts well for the existence of subgap currents seen in Fig. 2(b),inset. At low energies $`E\mathrm{\Delta }_{\mathrm{Nb}}`$ the density of states can be approximated as $`\rho (E,\mathrm{\Delta }_{Nb})\mathrm{\Gamma }/\mathrm{\Delta }_{\mathrm{Nb}}+(3\mathrm{\Gamma }/2\mathrm{\Delta })(E/\mathrm{\Delta }_{\mathrm{Nb}})^2`$. In the case of a single junction of resistance $`R`$ with a normal-state metal this results in a subgap conductance $`R^1\rho (V,\mathrm{\Delta }_{\mathrm{Nb}})`$; therefore, ohmic for low enough voltages and increasing as $`V^2`$ for relatively larger voltages. In addition, the granular structure of the films and, possibly, impurities resulting from outgassing of the mask polymer is also making the film inhomogenous. As a result, the gap edge is smeared due to local fluctuations in the effective electron-electron interaction. Other fabrication techniques \[e.g. the use of stronger polymers such as PES zorin \] could improve the quality of the Nb films. For our samples, since the superconducting coherence length in bulk Nb is only $`\xi _{\mathrm{Nb}}=38`$ nm, we expect to be in the limit in which the size of the inhomogeneities is much larger than $`\xi _{\mathrm{Nb}}`$. In this case, it has been shown smeared  that the effective (average) density of states takes the form

$$\overline{\rho }(E)=_0^{\mathrm{}}𝑑\mathrm{\Delta }_{\mathrm{Nb}}𝒫(\mathrm{\Delta }_{\mathrm{Nb}})\rho (E,\mathrm{\Delta }_{\mathrm{Nb}}),$$
(3)
where $`𝒫`$ is the probability density associated with a certain value of the gap. For small levels of inhomogeneity, the function $`𝒫`$ is a Gaussian of standard distribution $`\sigma `$ centered around an average gap value (slightly smaller than the bulk value).

Finally, to complete our model, we calculate numerically the currents and conductances for the superconducting SET in the framework of the orthodox theory of sequential tunneling, by computing the tunneling probabilities into and out of the island (with the density of states above for Nb) and by solving the master equation for the charge on the island. In the case of suppressed Al gap, a very good fit is obtained for $`\sigma =0.38`$ meV, $`\mathrm{\Gamma }=37.5`$ $`\mu `$eV, $`T=95.7`$ mK, and charging energy, $`E_\mathrm{C}=34.1`$ $`\mu `$eV, as determined before \[Fig. 2(b), main graph, black continuous line\]. The temperature also agrees well with the value obtained by simple CBT fitting above and the value of the thermometer.

In the case when the leads are superconducting the comparison with the experimental data can be done only qualitatively: the main reason is the existence of Cooper pair resonances peaks, yielding an extra contribution to the current which cannot be eliminated in a straightforward way. However, a number of qualitative features can still be observed by averaging the resonance peaks over gate voltages. The conductances obtained by this procedure should reflect the voltage dependence of the subgap current in the region where the height of the resonance peaks is approximately constant. A reasonably fair qualitative agreement with the model presented above is obtained \[Fig. 3(a)\]. The main feature that we see both in the calculated conductivities and in the conductivities obtained from our data through the above procedure is the appearance of a peak at the onset of the Al quasiparticle threshold ($`2\mathrm{\Delta }_{\mathrm{Al}}`$). By using a magnetic field to partially suppress the gap of Al, we see that the peak moves to the left \[see Fig. 2(a) and Fig. 3(a)\]. To verify that this is indeed the case, we have measured single Al-Al junctions with the same fabrication parameters as the leads of the Nb-island SET, and determined the gap of Al at various magnetic fields \[Fig. 3(b), circles\]. The agreement with the gaps determined from the low-voltage features \[Fig. 3(a)\] of three Nb-island SET samples (square, triangle, star) is good.

In conclusion, we have fabricated and measured superconducting single-electron transistors and we have presented evidence for a number of transport processes occurring at bias voltages below the gap of the island: resonant tunneling of Cooper pairs, Coulomb blockade of the subgap current for normal leads, and the appearance of a step in the subgap current at $`2\mathrm{\Delta }_{\mathrm{Al}}`$ for the case of superconducting leads.

###### Acknowledgements.

This work was supported by the Academy of Finland (Acad. Res. Grant no. 00857 and Projects No. 7111994, 7118122, 7205476), EU (IST-1999-10673, HPMF-CT-2002-01893), and the Finnish National Graduate School in Materials Physics. J.K. acknowledges the support of the Department of Energy, Office of Basic Energy Sciences via the Chemical Sciences, Geosciences, and Biosciences Division. G. S. P. would like to thank D. V. Averin, Yu. V. Nazarov, and I. S. Beloborodov for useful discussions.


text:
# A Fast Algorithm for Cosmic Rays Removal from Single Images

## 1 INTRODUCTION

Cosmic ray hits cause defects in all astronomical images obtained with CCD detectors. A relatively efficient approach to removing traces of cosmic rays from such images is to use multiple frames of the same object and then combine them with some algorithm for rejection of the outlying data. Methods of this type have been presented among others by Shaw & Horne (1992) and Windhorst, Franklin, & Neuschaefer (1994). They may be found in most astronomical data processing packages.

In many cases we are not able to obtain multiple images of the same object or the required time resolution prevented us from using multiple image methods. A straightforward method that could be used in such cases would be sigma clipping. In practice it is almost impossible to find a good detection threshold for such an algorithm. In fact, it may leave some obvious cosmic rays untouched, while giving numerous spurious detections and rejecting valid data at the same time. Other types of single image algorithms may also be found in most astronomical data processing packages. Most of these methods rely on the sharpness of cosmic rays, relative to the atmospheric smoothing of real images. These methods are based on some interactive learning techniques (Murtagh, 1992; Salzberg et al., 1995) or involve a special fitting of a model representing a real image of a star with a superimposed cosmic ray (Rhoads, 2000). These methods were designed to correct direct images of the stellar objects and we find them not suitable for spectral images.

The exposure times of the spectra are usually longer than direct images and range from 600 seconds to 1800 seconds. Considerable amounts of cosmic ray events accumulate in the images during such exposures. We have constructed a simple and straightforward algorithm for detecting cosmic rays in single images. Our method does not need any model of the shape of the image features themselves. We analyze the histograms of pixel counts in small sub-frames in order to detect pixels deviating by some factor from the bulk of the pixels under consideration.

We would like to note here, that low value bad pixels tend to hide weak cosmic rays by artificially making the local variance larger than would be caused only by usual sources of noise in the image: readout noise, dark current and the Poisson noise associated with background and object signal. The cosmic rays have random distributions of their positions in the images. On the other hand, bad pixels may be found as deviating pixels at the same position in every frame taken with a given instrument. Therefore we suggest to make bad pixel correction as a separate step of the reduction process prior to cosmic rays correction.

We describe details of our algorithm in paragraph 2. In paragraph 3. we describe an implementation and the results of our tests. In paragraph 4. we present examples of application to real astronomical data. Paragraph 5 contains the summary.

## 2 The Algorithm

The proposed algorithm is based on the idea that cosmic rays deposit a portion of their energy in the pixels they hit, causing some extra signal in these pixels. The signal coming from cosmic rays does not have a Gaussian distribution. This should reflect in the distribution of counts in the image affected by cosmic rays. The image as a whole may have a large range of signal levels in different areas. We analyze relatively small sub-frames, to work with a more concentrated local distribution of counts. In most cases this distribution should be rather compact. The cosmic rays appear as points standing out in the high counts interval.

The proposed algorithm of cosmic rays detection generally consists of the following steps:

1. Select small size sub-frames, which cover the whole frame with substantial overlapping.

In each sub-frame:

1. Calculate standard deviation of the distribution of counts:

$$\sigma =\sqrt{\frac{c_i^2\frac{(c_i)^2}{n}}{n}}$$
(1)
2. Apply a single sigma clipping step to correct the estimate of standard deviation for outlying pixels.
3. Construct a histogram of the distribution of counts.
4. Find the mode of the distribution of counts, i.e. the peak of the histogram.
5. In the interval of counts higher than the mode, find gaps in the histogram i.e. bins with zero data points.
6. Find the first gap which is wider than a threshold, which is the standard deviation multiplied by an arbitrary number (usually 3.0).
7. If such gap exists, flag pixels with counts lying above the gap as affected by cosmic rays.

Figure 1 illustrates an example of a histogram with a few pixels identified as affected by cosmic ray events.

The next step consists of a replacement of the count numbers in the flagged pixels. If we consider a single image, the information about the real signal in the affected pixels is lost. In many cases however, the characteristic scale of spacial variations of signal are of the order of at least a few pixels. In such cases, one may use some interpolation to replace missing pixel counts. In our implementation, we decided to substitute the cosmic rays by the average of the counts in the neighboring pixels.

Cosmic rays are often multiple-pixel events. High signal pixels may screen neighboring pixels from detection. For this reason we run the procedure of detection and cleaning cosmic rays iteratively. The process is rapidly converging and usually there are no new detections after 2 or 3 iterations. We have also introduced a parameter called ”growing radius”, which tells the procedure to clean the pixels closer than this radius to the detected one, even when they do not possess enough signal to be flagged as affected by cosmic rays.

## 3 Implementation and Test

We have written a computer program, which implements the above described algorithm111The C-source of our program may be downloaded from the following Web-page: http://www.camk.edu.pl/{̃}pych/. This program has several input parameters which allow the user to control the process of detecting and cleaning cosmic rays. These parameters define:

* the size of the sub-frame box: X and Y dimensions of the rectangle,
* the threshold: the number by which the local variance is multiplied,
* the dispersion axis: the interpolation for bad pixel substitution is calculated along this axis; if no axis selected, it is calculated in the annulus around that pixel,
* the lower and upper radii of the annulus for the interpolation of data,
* the growing radius: the maximum distance of the pixels which are to be corrected in the neighborhood of the flagged pixel.

Our program also creates a map of detected cosmic rays. This map can be examined to check for any correlation with the real data. The existence of such correlation indicates that the detection threshold has been set too low and the real data has been modified.

To check the capabilities of our method, we have generated an artificial two-dimensional Echelle spectrum with the size of 2048 x 4096 pixels. For this purpose, we have used a procedure from IRAF package: noao.artdata.mkechelle, with its default parameters. We have also used the procedure noao.artdata.mknoise to obtain two noise images. The first one consists of readout noise (5 electrons) and a background of 500 counts with Poisson noise. The second image consists of 1000 cosmic rays with maximum energy 30000 electrons.

Our program finished with zero detections on the noiseless image of the Echelle spectrum. The result for the spectrum with the readout and Poisson noise added was one detection over the whole frame. Finally, we tested the procedure on a spectrum with Poisson noise and cosmic rays added to the original image of a spectrum. Our program has found all but one cosmic ray event of peak greater than 200 counts. The only event left, and with a peak of 430 counts, was located in the vicinity of a trace of the spectrum and was concealed by this signal.

Figure 2 presents histograms of the counts distribution. The upper-left panel (A) presents the original distribution of counts from cosmic rays. The Y axis is cut at a value of 10 pixels to highlight the distribution of high energy hits. The upper-right panel (B) presents the distribution of cosmic rays in the interval of small counts. The lower panels present the distribution of counts in the image obtained as a difference between the spectrum image after cosmic rays cleaning and the image without cosmic rays added. The Y-axis in the lower-left panel (C) is cut at a value of 10 pixels, to show the distribution of pixels with undetected cosmic rays. The features visible above 200 counts are the traces of a single undetected cosmic ray event. The lower-right panel (D) presents the distribution of the residual counts. It has the shape of a Gaussian with a small ’high-end’ tail, produced by undetected cosmic rays. The standard deviation of the residuals, $`\sigma =26.0`$, is similar to the variance of the original image without cosmic rays; $`\sigma =27.8`$. This similarity reflects the nature of the adopted substitution method for the pixels affected by cosmic rays, which is simply a mean value of the surrounding pixels.

Our tests were conducted on a PC with a 500 MHz Intel Celeron processor. The program was compiled using a GNU C Compiler (gcc version 3.2.2) under Red Hat Linux 9.0. The CPU time needed to process a single frame: 2048 x 4096, 32 bit pixels, was about 40 seconds. The CPU time required to process an image depends linearly on the number of the image pixels.

We have also run the IRAF procedure noao.imred.crutil.cosmicrays. The CPU time needed to process the same image was over 900 seconds. In comparison with the above, our algorithm may be classified as a fast one. The relative speed results from the fact that the whole detection process does not require extensive calculations. Our test demonstrated that the IRAF routine was not able to remove multiple pixel cosmic ray events. They were slightly modified, but most of them remained in the frame. At the same time, a large number of image pixels unaffected by cosmic rays were modified. We would like to stress here that our algorithm does remove cosmic rays effectively while leaving almost all of the image data untouched.

Figure 3 presents central fragments of our test images. The upper-left panel shows the artificial Echelle spectrum with the Poisson noise and cosmic rays added. The upper-right panel presents the residual signal remaining after subtraction of the original spectrum without cosmic rays from the image cleaned by our program. The lower-left panel presents the image after being processed by the IRAF procedure noao.imred.crutil.cosmicrays. Note that the cosmic rays were not removed. The lower-right panel shows the residual signal remaining after subtraction of the original spectrum without cosmic rays from the image obtained from the IRAF routine noao.imred.crutil.cosmicrays. The scales of the images of the residual signals are the same.

## 4 Examples of practical application

We have applied the above algorithm to several data sets. Our general experience is that the method works very well on images of both long-slit and Echelle spectra. The algorithm was originally designed to work with the spectroscopic data. Our method could be used to clean images of stellar fields with a very wide Point Spread Function (PSF). It works well in a sense that it does not produce too many spurious detections. The limitation of the presented method for these types of images comes from the fact that direct images usually have large count variations within small scales. This produces a large standard deviation of the counts and prevents cosmic rays in the neighborhood of bright objects from being detected. We suggest using the method proposed by Rhoads (2000) in such cases.

Figure 4 illustrates an example result of cosmic ray removal in a long-slit spectrum. A fragment of the spectrum of a 9.2 mag. variable star, obtained in a 1200 second exposure with the 1.88-m telescope at the David Dunlap Observatory, University of Toronto, is shown. The left panel presents the original image (after bias subtraction and flat-field correction), while the right panel shows the cleaned image. No evident cosmic rays are identified by eye in the cleaned frame.

Our method works best on low signal images. For the faint stars however, we encounter another problem. The night-sky emission lines, which are sharp features perpendicular to the dispersion axis, may sometimes be strong enough to be identified as cosmic rays. Our solution to this problem is to edit the map of cosmic rays, erase features identified as sky emission lines (replace the counts with zeros) and then subtract a modified cosmic ray map from the original image.

Figure 5 illustrates an example of the application on an Echelle spectrum of a faint star. The upper-left panel shows a fragment of the Echelle spectrum of a 17 mag. star, obtained in an 1800 second exposure with MIKE spectrograph (Magellan Inamori Kyocera Echelle) at the Las Campanas Observatory. The upper-right panel presents the map of cosmic rays detected in the image. A sky emission line is visible in the center. The lower-left panel presents the same fragment of the spectrum as shown in the upper-left panel, but after subtraction of the modified map of cosmic rays. The map was edited to remove the sky emission line. The lower-right panel shows the map of cosmic rays with the sky emission line removed.

## 5 Summary

We have presented a cosmic-ray rejection algorithm based on a simple analysis of the histogram of the image data. The most important advantage of our method is that it does not require modeling of the image data and may apriori be applied to any type of well sampled image data. We have checked that for the spectroscopic images, it is very effective in detecting cosmic rays, while avoiding numerous spurious detections. Our method does not require advanced and extensive computations, so it is relatively fast.

The weak point of this approach, is that bright objects may shield cosmic rays in their neighborhood from detection. This is caused by two factors. First: the sensitivity to cosmic ray events is reduced at the locations of bright objects because of the Poisson noise associated with the image photons. Second: we look for bright spots and the legitimate object may be brighter than a nearby comic ray trace.

Our suggestion is to use the algorithm presented above for spectroscopic data whenever multiple image methods cannot be employed for cosmic rays removal.

Future improvement to the presented method could be an introduction of better interpolation for the replacement of pixels affected by cosmic rays.

The author would like to thank Prof. S. Rucinski, J. R. Thomson and H. DeBond for help in preparing this article. This article has been prepared when W. Pych held the NATO Post-Doctoral Fellowship administered by the Natural Sciences and Engineering Council of Canada (NSERC). The author acknowledges also the support from the Polish Grant KBN 2 P03D 029 23 and the NSERC research grant to S. Rucinski.