text:
# Probing the Ionization State of the Universe at 𝑧>6 Based on observations obtained with the W. M. Keck Observatory, which is jointly operated by the California Institute of Technology and the University of California.

## 1 Introduction

The discovery of quasars at redshifts greater than 6 from the Sloan Digital Sky Survey (SDSS; Fan et al. 2001, Fan et al. 2003) has opened a window for spectroscopic exploration of the high redshift intergalactic medium. The detection of broad, black, Ly$`\alpha `$ absorption troughs, as predicted by Gunn & Peterson (1965; hereafter GP), indicates that we are beginning to probe the era when hydrogen was reionized by an early generation of stars. Becker et al. (2001) detected the first complete GP trough in the spectrum of SDSSp J103027.10$`+`$052455.0 ($`z=6.28`$; hereafter SDSS J1030$`+`$0524). They compared the absorption to that seen in the spectra of lower redshift quasars and concluded that the fraction of neutral hydrogen in the IGM increases sharply at $`z6`$. Fan et al. (2003) report that the most distant known quasar, SDSS J114816.64$`+`$525150.3 ($`z=6.37`$111This is a new redshift, different from that given in Fan et al. (2003), based on the spectrum presented in this paper. See §5.4 for a discussion of the $`z=6.41`$ redshift derived by Willott, McLure & Jarvis (2003) for this object.; hereafter SDSS J1148$`+`$5251), also shows a complete GP trough in its spectrum. Djorgovski, Castro, Stern & Mahabal (2001) argued that the spectra of even lower redshift quasars show indications of neutral regions in the IGM.

There ensued a spate of theoretical papers discussing the details of how reionization takes place. But the most exciting new observational result is undoubtedly the discovery by the Wilkinson Microwave Anisotropy Probe (WMAP) that the microwave background has traversed a moderate optical depth ($`\tau =0.17\pm 0.04`$) electron scattering medium, implying that the universe must have been ionized at $`z10`$–20 (Kogut et al. 2003). At first glance this appears inconsistent with reionization at $`z6`$, but since a relatively small neutral fraction ($`1`$%) would suffice to produce black GP troughs, the results are not contradictory. There is mounting evidence that the ionization history of the intergalactic medium (IGM) was complex. Wyithe & Loeb (2003) suggested that the hydrogen in the IGM could have been reionized twice, and Cen (2003) argued that two phases of reionization are likely under a wide range of conditions. In this scenario, the first reionization occurs at $`z15`$, driven by the formation of massive, zero-metallicity Population III stars (Venkatesan, Tumlinson & Shull 2003; Mackey, Bromm & Hernquist 2003). The ionizing radiation from massive stars is dramatically reduced, however, when Pop III supernovae seed the gas with metals (due to a reduction in the formation of massive stars, increased line-blanketing in hot star atmospheres, and increased mass loss rates from both hot stars and red giants.) Consequently the IGM recombines and remains relatively neutral, $`n_{\mathrm{HI}}/n_\mathrm{H}>0.1`$, until an increasing Pop II massive star population and a declining IGM density allow the second reionization to occur at $`z6`$. This neutral fraction is sufficient to produce the observed high optical depths in the Ly$`\alpha `$ GP troughs of high redshift quasars.

To advance our understanding of the ionization history of the universe, we need high resolution, high signal-to-noise ratio spectra of the IGM absorption at $`z>6`$. Recently several groups have discovered Ly$`\alpha `$-emitting galaxies at $`z6.5`$ (Hu et al. 2002, Kodaira et al. 2003). These objects are interesting and relevant as direct probes of star and galaxy formation at $`z>6`$, and the mere detectability of their Ly$`\alpha `$ emission is an important clue to the ionization state of the IGM (Haiman 2002). However, their continua are far too faint (even redward of Ly$`\alpha `$) to make them useful for studies of the intervening IGM. Only high-redshift quasars offer the possibility of acquiring high quality information on the state of the IGM at high redshift.

In this paper we present high quality spectra of the only two quasars known to have GP absorption troughs222SDSS J1048$`+`$4637 at $`z=6.23`$ (Fan et al. 2003) could also have a GP trough, but our recent spectra show that it has strong, intrinsic, broad absorption lines that make it less useful for studies of the IGM at $`z>6`$ (Fan et al. 2003b)., SDSS J1030$`+`$0524 ($`z=6.28`$) and SDSS J1148$`+`$5251 ($`z=6.37`$), taken with the ESI spectrograph on the Keck Telescope. The following sections describe the observations and data reduction (§2), discuss the analysis methods used (§3), present the properties of the IGM absorption in the spectra (§4), and discuss the results (§5). Appendix A explores the effect of light travel time on the observed expansion of ionization fronts around high-redshift quasars.

## 2 Observations

Spectra of SDSS J1030$`+`$0524 and SDSS J1148$`+`$5251 were taken with the Keck Echellette Spectrograph and Imager (ESI; Sheinis et al. 2002) between January 2002 and February 2003. ESI covers the wavelength range 4000–11000 Å with 10 Echelle orders having a constant resolution of $`R=\lambda /\mathrm{\Delta }\lambda 10000`$ ($`30`$ km/s). The data were taken with the slit at the parallactic angle in mostly good but variable seeing (from 0.6″ to $`>1`$″, with mean seeing around 0.8″). The sky transparency ranged from photometric to light cirrus. A log of the observations is given in Table 1.

The most challenging data reduction problem for these observations is subtracting the sky background, including the bright OH lines that dominate the near infrared atmospheric emission. The sky is considerably brighter than the object, and small errors in the sky subtraction can easily produce offsets in the emission that either create the appearance of light in the absorption trough or remove light that should be present there. To assist in the sky subtraction, we shifted the object along the slit for the various integrations, which helps to average over any small-scale flat-fielding, bias or dark current variations not corrected by the standard CCD calibrations.

The sky brightness varies as a function of time in a complex way in the near infrared, with the atmospheric OH emission varying rapidly and independently of other sky emission. In addition the ESI spectrograph, though a highly stable instrument, does undergo small changes over time, with the wavelengths shifting slightly from one observation to the next. Presumably these changes (which are typically at the level of a few tenths of a pixel or less) are due to minor thermal and gravity-induced flexures.

We extract a sky-subtracted spectrum from each independent observation. The separate spectra are then combined using weights to optimize the signal-to-noise ratio of the result. It is necessary to be very careful with the weights for combining to avoid any bias in the results. A simple unweighted average is safely free of bias but leads to a much noisier spectrum than an optimally weighted sum. On the other hand, using weights $`1/\sigma ^2`$ with $`\sigma `$ derived from the noise estimates of the individual spectra leads to a very noticeable bias toward lower flux values, because the estimated noise is larger for brighter pixels. In the black absorption troughs this leads to fluxes that are systematically negative.

The dangers of using a weighted combination for Poisson data are well-documented (for a nice analysis, see Wheaton et al. 1995). The cure for the bias is to choose weights that are completely independent of the data. We have combined the spectra using weights that are derived from the sky noise in adjacent pixels (i.e., along the slit). Note that using the sky noise from the same pixel would still lead to a small bias, because pixels that have positive sky fluctuations (making the residual counts a little low) would be given slightly lower weights, thus leading to a positive bias in the sum. Using the sky noise in neighboring pixels produces a completely unbiased result that has a nearly optimal signal-to-noise ratio. An overall scale factor is determined for each spectrum to match the fluxes, thereby correcting for variable atmospheric transparency.

Our extraction algorithm is highly customized for ESI data, which have a number of unusual characteristics (e.g., highly curved orders, a strongly non-linear wavelength dispersion, and some odd notches in the sensitivity in the bluer orders.) The images are first geometrically rectified to produce orders with orthogonal wavelength and spatial coordinates lying along rows and columns. This removes the very large curvature of the ESI Echelle orders and also corrects small variations in the wavelength dispersion along the slit that manifest themselves as tilted or curved sky lines. The order curvature removal also automatically corrects the spectral trace for differential atmospheric refraction using the airmass of the observation; this approach is more effective than tracing the orders directly because high-redshift quasar spectra have long stretches with no detectable light, foiling standard tracing algorithms.

The resampling algorithm used in the rectification is a flux-conserving interpolation scheme that is accurate to second order in the pixel shifts (as opposed to the more commonly used schemes that are only first order, making them considerably more dispersive.) The interpolation method is based on the high-order, monotonic, flux-conserving advection schemes that are widely used in modern hydrodynamic simulations (van Leer 1977).

Extraction of spectra from the rectified orders is straightforward. We use an optimal extraction algorithm with a second-order fit to the sky along the slit and a spatial profile that is assumed constant for each order. For orders where the profile is not detected, the trace position and width are predicted based on those orders that are detected. The wavelength scale derived from calibration lamps is adjusted using the positions of sky lines for each observation. The final spectra for the quasars are shown in Figure 1.

## 3 Analysis of the Spectra

High signal-to-noise ratio spectra promise two improvements in our understanding of the GP troughs. By reducing the noise in the black regions of the spectrum, they improve our lower limits on the absorbing optical depth and may allow us to detect faint light in the dark troughs. But since the optical depth limits increase only as the logarithm of the signal-to-noise ratio, which itself increases only as the square root of the integration time, it is difficult to make dramatic improvements in the limits. The main benefit of high signal-to-noise ratio spectra is that they allow us to set strong detection limits over narrower bands in the spectra. The optical depth distribution is very unlikely to be a smooth function of redshift; instead, at the transition between a neutral and an ionized IGM, there will be ionized, partially transparent bubbles sprinkled along the line of sight. To detect these bubbles we need to be able to recognize narrow emission spikes in the trough where the quasar light leaks through the IGM. The new spectra presented in the paper are far better for that purpose than our earlier observations (which were 30 min exposures taken with ESI.)

We begin with an approach similar to that used by Becker et al. (2001). The residual fluxes in the GP troughs of Ly$`\alpha `$ and Ly$`\beta `$ are measured directly from the spectra; the fluxes are weighted using the sky noise from adjacent pixels (as described above for the extraction) to improve the signal-to-noise ratio while avoiding bias in the sums. Then the mean transmission and optical depth of the IGM are determined by comparing the residual flux to the (unabsorbed) quasar continuum at the positions of the troughs.

Because very little neutral hydrogen is required to create a GP trough in Ly$`\alpha `$ (a 1% neutral fraction will suffice), the presence of a GP trough due to Ly$`\alpha `$ provides only a lower limit to the amount of neutral hydrogen. The GP trough from the Ly$`\beta `$ transition can be a more sensitive test because the Ly$`\beta `$ cross-section is a factor of 5.27 lower than Ly$`\alpha `$. Interpretation of the Ly$`\beta `$ trough is complicated by the fact that it is overlaid by the Ly$`\alpha `$ forest at $`z5`$, which is unpredictable along a given line of sight due to cosmic variance. In spite of that, using fairly secure assumptions, the Ly$`\beta `$ trough can give much stronger limits on the neutral hydrogen fraction. We consequently measure residual fluxes in both the Ly$`\alpha `$ and Ly$`\beta `$ absorption regions.

We also considered using the Ly$`\gamma `$ absorption since Ly$`\gamma `$ has an even lower intrinsic optical depth. However, the Ly$`\gamma `$ GP trough suffers from overlying absorption by both Ly$`\alpha `$ and Ly$`\beta `$ at lower redshifts and consequently does not give additional useful information.

It is necessary to estimate the brightness of the quasar continuum at the Ly$`\alpha `$ and Ly$`\beta `$ troughs in order to measure the transmission of the IGM. For this purpose we use the far UV quasar power law spectrum of Telfer, Zheng, Kriss, & Davidsen (2002) and the composite spectrum from the Large Bright Quasar Survey (LBQS; Francis et al. 1991333We use the improved LBQS template created by S. Morris, described in Brotherton et al. (2001) and available at http://sundog.stsci.edu/first/QSOComposites. Note that the LBQS composite is very similar to the FIRST Bright Quasar Survey composite derived by Brotherton et al.). The templates are matched to the observed spectra at a rest wavelength of $`1290`$ Å. The Telfer et al. power law has a break near 1280 Å and rises less steeply into the UV than a simple extrapolation of the continuum redward of Ly$`\alpha `$. In the near UV ($`\lambda >1280`$ Å), the spectral slope is $`\alpha _{NUV}=0.72`$ ($`f_\nu \nu ^\alpha `$), while in the far UV ($`\lambda <1280`$ Å), $`\alpha _{EUV}=1.57`$ for radio-quiet objects and $`1.96`$ for radio-loud objects. While there is no evidence that either of our quasars is radio-loud, we compute continuum levels using both EUV slopes in order to represent the range of continua seen in quasars.

Broad emission lines can also contribute to the continuum. Our spectra are displayed in Figures 2 and 3 along with the LBQS composite spectrum. For these figures, our spectra have been denoised using a wavelet shrinkage algorithm (Donoho & Johnstone 1994) that is optimized for the case where the noise amplitude is variable but known (heteroscedastic data, in the parlance of statisticians). We use a Haar (1910) wavelet transform, which allows error propagation and has the great benefit of keeping the noise in transform space uncorrelated at each scale before the filtering is applied. Cycle-spinning (Coifman & Donoho 1995) is utilized to make the result shift-invariant and to remove the artifacts generated by the square-wave Haar basis functions.

Denoising thus uses a variable smoothing length to suppress the noise as much as possible on all scales while remaining consistent with the spectrum given the estimated noise array. The resulting spectrum is heavily smoothed in flat regions (such as the dark troughs), which suppresses the noise. (Values in the smoothed areas are consequently highly correlated.) The smoothing threshold is set to $`3\sigma `$, so only structures more significant than that remain in the spectrum. Note that high resolution structures are preserved (compare the absorption on the blue edge of the Ly$`\alpha `$ emission with Fig. 1) while the noise in both the continuum and the troughs is dramatically reduced.

For SDSS J1030$`+`$0524 (Fig. 2), the LBQS composite is clearly a good match to the spectrum, so we estimate the continuum in the Ly$`\alpha `$ and Ly$`\beta `$ absorption troughs using the mean flux in the relevant wavelength window from the LBQS spectrum. On the other hand, SDSS J1148$`+`$5251 shows considerably weaker emission lines than the composite (Fig. 3); we therefore use the Telfer et al. power law continuum alone as the estimate of the brightness at the absorption troughs. The continua estimated by using these different templates differ by about 30%.

## 4 Results

### 4.1 SDSS J1030$`+`$0524

For SDSS J1030$`+`$0524, our results can be described very simply: the GP absorption troughs of both Ly$`\alpha `$ and Ly$`\beta `$ are very black. There is no evidence of any residual light in the troughs to the limit of our observations. The Ly$`\alpha `$ and Ly$`\beta `$ absorption trough details are shown in Figure 4.

The results are summarized in Table 2. The band limits were chosen to exclude the first bright peaks on either end of the Ly$`\beta `$ GP trough and are identical in redshift for the two troughs. The residual fluxes for both the Ly$`\alpha `$ and Ly$`\beta `$ troughs are consistent with zero.

For the Ly$`\beta `$ trough, the continuum estimate includes a reduction factor of 0.12 for the overlying Ly$`\alpha `$ forest absorption at $`z=5`$ (as determined by Becker et al. 2001) The $`\tau (\text{Ly}\beta )`$ values are converted to equivalent Ly$`\alpha `$ optical depths by multiplying by 5.27, the ratio of the oscillator strengths of the two lines. Thus the Ly$`\beta `$ GP trough measurements for SDSS J1030$`+`$0524 imply $`\tau (\text{Ly}\alpha )>22.8`$ ($`1\sigma `$). The errors quoted in the table reflect only the statistical uncertainties. As found by Becker et al., the strongest limit on the optical depth comes from the Ly$`\beta `$ trough, despite the additional uncertainty that comes from the overlying Ly$`\alpha `$ absorption.

These results are consistent with the measurements of Becker et al. (2001) (confirmed by Pentericci et al. 2002) but are naturally stronger given the much improved signal-to-noise ratio of the data. Our limits on the residual flux in the GP troughs and the Ly$`\alpha `$ optical depth of the IGM are lower by the expected factors given the longer integrations reported here. The limits remain consistent with zero flux in the GP troughs of SDSS J1030$`+`$0524.

There is no evidence for any significant peaks in either trough. Figure 5 displays the signal-to-noise ratio of the spectra in the troughs. The distribution is consistent with zero-mean noise in the troughs; the reduced $`\chi _\nu ^2`$ for the combined Ly$`\alpha `$ and Ly$`\beta `$ troughs is 0.97.

### 4.2 SDSS J1148$`+`$5251

The GP troughs for SDSS J1148$`+`$5251 are also shown in Figures 4 and 5. The mean transmission properties of the troughs are given in Table 2, which includes the transmission measured for both a high redshift range ($`6.1<z<6.32`$) and a lower redshift interval ($`6.0<z<6.1`$). Note that we now quote a measurement, not a lower limit, on $`\tau `$, as flux has been detected in the troughs. The correction factor for Ly$`\alpha `$ forest absorption overlying the Ly$`\beta `$ GP trough was taken to be 0.11, from Figure 2 of Becker et al. (2001). The $`1\sigma `$ errors on the transmission are smaller than for SDSS J1030$`+`$0524, both because the seeing and photometric conditions were slightly better for SDSS J1148$`+`$5251 and because its continuum is brighter.

In contrast to SDSS J1030$`+`$0524, there is light detected in the GP troughs of SDSS J1148$`+`$5251 both in narrow spikes and in broader bands. The mean transmission nonetheless remains low. Figure 6 summarizes the Ly$`\alpha `$ optical depth as a function of redshift, including additional high-redshift quasars from Fan et al. (2001). There is a noticeable steepening of the $`\tau (z)`$ relationship at $`z5.8`$.

A closer examination of the troughs reveals some interesting features. The Ly$`\alpha `$ trough extends from $`z=5.93`$ to 6.32 and looks dark at high resolution except for a weak emission spike at $`z=6.08`$. The reality of that peak is supported by its coincidence with a Ly$`\beta `$ emission peak at the same redshift (Fig. 4). We believe this is likely to be a real detection of an ionized bubble in the IGM that permits quasar light to leak through at the redshift of the bubble; this feature is discussed in detail below (§5.1).

Several emission spikes are detected in the Ly$`\beta `$ trough as well. Two especially strong peaks are seen at 7205 Å and 7236 Å (Ly$`\beta `$ redshifts of $`z=6.03`$ and 6.06). There also is significant non-zero flux in the troughs even between the peaks (best seen in Fig. 5). If this flux is quasar light transmitted by the IGM, it certainly suggests strongly that the IGM is not very neutral at $`z6`$. That interpretation is not so clear cut, however. We argue below (§5.2) for the alternative interpretation that both the Ly$`\alpha `$ and Ly$`\beta `$ troughs are contaminated by light from an intervening galaxy at $`z=5`$. If that suggestion is correct, the optical depths derived from the troughs of SDSS J1148$`+`$5251 are only lower limits.

## 5 Discussion

### 5.1 The $`z=6.08`$ emission feature in SDSS J1148$`+`$5251

The universe at $`z6`$ was very likely ionized by galaxies rather than by quasars, since there do not appear to be enough quasars to do the job (Fan et al. 2001). There must have been regions in the preionized universe where a concentration of unusually luminous galaxies (possibly the precursor to a massive cluster) was capable of ionizing the nearby IGM, creating an ionized bubble. It is therefore not surprising to find weak emission spikes embedded in the mostly black GP troughs (Figs. 4 and 5) where the line-of-sight to the quasar passes through small “holes” in the IGM Ly$`\alpha `$ and Ly$`\beta `$ transmission created by ionized bubbles. That can also explain the detection of Ly$`\alpha `$-emitting galaxies at $`z6.5`$ (Hu et al. 2002, Kodaira et al. 2003), which in this picture ought to be found in concentrations with other galaxies.

According to the scenario presented by Cen (2003), the universe is first ionized at $`z15`$ by Pop III stars and then recombines to become substantially neutral, with neutral fractions remaining above 0.1 until the second reionization occurs at $`z6`$. At $`z=6`$ to 7 the neutral fraction $`n_{\mathrm{HI}}/n_\mathrm{H}`$ is typically $`0.15`$. As Cen points out, that makes it easier to detect high redshift galaxies; it also makes it easier for lower luminosity sources to create ionized bubbles.

The peak seen in the Ly$`\alpha `$ trough of SDSS J1148$`+`$5251 at 8600 Å ($`z=6.08`$) has an intensity about 8% of the extrapolated (pre-absorption) continuum, implying an IGM optical depth of $`\tau (\text{Ly}\alpha )=2.5`$. The optical depth of this window at Ly$`\beta `$ would be only about 0.5, which would produce a very strong peak in the Ly$`\beta `$ trough; but since the overlying Ly$`\alpha `$ forest absorption is expected to reduce that by a factor of 0.11 (corresponding to an overlying optical depth $`\tau =2.2`$), the expected amplitude of the Ly$`\beta `$ trough peak is in fact very similar to that observed, albeit with a rather large uncertainty in the Ly$`\alpha `$ forest attenuation.

How large would an ionized region need to be to produce this feature? The damping wings from neutral hydrogen on either side of the ionized bubble create a high optical depth through the region unless the bubble is relatively large. This led Barkana (2002) to question Djorgovski et al.’s (2001) argument that dark regions in the spectrum of a lower redshift were created by neutral regions of the IGM. Miralda-Escudé (1998; see also Barkana 2002) derived the Ly$`\alpha `$ damping wing optical depth for hydrogen extending from $`z_1`$ to $`z_2`$ when observed at a wavelength $`\lambda _\alpha (1+z)`$:

$$\tau _{\mathrm{damp}}=6.43\times 10^9\tau _{\mathrm{GP}}\left(\frac{n_{\mathrm{HI}}}{n_\mathrm{H}}\right)\times \left[I\left(\frac{1+z_2}{1+z}\right)I\left(\frac{1+z_1}{1+z}\right)\right],$$
(1)
where $`\tau _{\mathrm{GP}}`$ is the Gunn-Peterson optical depth,

$$\tau _{\mathrm{GP}}=3.42\times 10^5\left(\frac{1+z}{7.08}\right)^{3/2},$$
(2)
and we are using a $`\mathrm{\Lambda }`$ cosmology with $`H_0=70`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_m=0.3`$, and $`\mathrm{\Omega }_b=0.04`$. The function $`I(x)`$ is defined to be

$$I(x)=x^{1/2}\left(\frac{x^4}{1x}+\frac{9}{7}x^3+\frac{9}{5}x^2+3x+9\right)\frac{9}{2}\mathrm{ln}\left|\frac{1+x^{1/2}}{1x^{1/2}}\right|.$$
(3)
Given an ionized region at redshift $`z`$, Eqn. (1) allows us compute the combined damping wing optical depth from the neutral regions on either side of the bubble. Although the constant factor multiplying $`\tau _{\mathrm{GP}}`$ in Eqn. (1) is small, the leading term in $`I(x)`$ can be very large near $`z=z_1`$ or $`z_2`$, so $`\tau _{\mathrm{damp}}`$ can be substantial.

Figure 7 displays the resulting optical depth for a bubble located at $`z=6.08`$. The neutral IGM is assumed to extend from $`z=6`$ to $`z=6.32`$, but the results are relatively insensitive to those limits. The bubble diameter required to produce a window with $`\tau =2.5`$ is 1.33 Mpc if the IGM is completely neutral outside the bubble. If the IGM is already partly ionized, the required bubble size is much smaller (0.22 Mpc if $`n_{\mathrm{HI}}/n_\mathrm{H}=0.15`$).

Barkana (2002) computed the diameter of an H II bubble as a function of the halo mass $`M`$ of the ionizing source (including the escape probability for ionizing photons in a low metallicity galaxy) to be

$$D_S=0.150\mathrm{Mpc}\left(\frac{7.08}{1+z}\right)\left(\frac{N_{\mathrm{ion}}}{20}\frac{M}{10^9\mathrm{M}_{\mathrm{}}}\right)^{1/3},$$
(4)
where $`N_{\mathrm{ion}}`$ is the number of ionizing photons per baryon that escape the galaxy. We can invert this to get the mass required to produce a bubble of a given diameter:

$$M=7\times 10^{11}\mathrm{M}_{\mathrm{}}\left(\frac{D_S}{1.33\mathrm{Mpc}}\frac{1+z}{7.08}\right)^3\frac{20}{N_{\mathrm{ion}}}.$$
(5)
We conclude that a plausible star-forming galaxy (or group of galaxies) along the line of sight can produce the observed peak at $`z=6.08`$ in the Ly$`\alpha `$ and Ly$`\beta `$ GP troughs of SDSS J1148$`+`$5251, even if the surrounding IGM is completely neutral. Obviously if the IGM is partly ionized, the required bubble diameter would be much smaller (Fig. 7) and so the required galaxy mass would also be far smaller ($`<10^9\mathrm{M}_{\mathrm{}}`$). With future instrumentation it will be interesting to search the vicinity of the quasar for evidence of star-forming galaxies at $`z=6.08`$ that are the source of the ionizing continuum radiation required to make this bubble.

### 5.2 Evidence for an intervening $`z=5`$ galaxy

The strong spikes seen in the Ly$`\beta `$ trough of SDSS J1148$`+`$5251 at 7205 Å and 7236 Å are more difficult to explain as the result of ionizing bubbles in the IGM. The strength of these peaks is not easily reconciled with the absence of corresponding Ly$`\alpha `$ peaks, given the expected overlying Ly$`\alpha `$ forest absorption. The peak at 7205 Å has an amplitude $`8.7\pm 0.3`$ ($`\times 10^{18}\text{erg}\text{s}^1\text{cm}^2\text{Å}^1`$), while the brightest plausible continuum at that wavelength is 24 (from the LBQS template spectrum in Fig. 3). If we assume that the 7205 Å peak is quasar light leaking through the $`z=6`$ Ly$`\beta `$ and $`z=5`$ Ly$`\alpha `$ absorption, the implied combined optical depth is $`\tau (\text{Ly}\alpha \text{@}z=5)+\tau (\text{Ly}\beta \text{@}z=6)=1.0`$. But there is no hint of a Ly$`\alpha `$ peak at 8539 Å: the flux limit there is $`0.1\pm 0.2`$ with an estimated continuum level of 17, implying $`3\sigma `$ upper limits $`\tau (\text{Ly}\alpha \text{@}z=6)>3.3`$ and $`\tau (\text{Ly}\beta \text{@}z=6)=\tau (\text{Ly}\alpha \text{@}z=6)/5.27>0.6`$. The overlying Ly$`\alpha `$ forest absorption must then be $`\tau (\text{Ly}\alpha \text{@}z=5)<0.4`$.

We conclude that if the 7205 Å peak is the result of an ionized bubble at $`z=6.03`$, it must fall in a highly transparent window in the $`z=5`$ Ly$`\alpha `$ forest that barely attenuates the continuum at all, allowing $`>65`$% of the light through. High resolution spectra of the Ly$`\alpha `$ forest at $`z5`$ (Songaila, Hu, Cowie & McMahon 1999; Djorgovski et al. 2001; Songaila & Cowie 2002) do show small transparent windows, but they are rare. The spectra presented in this paper are typical. There is one transparent window seen at 7800 Å ($`z=5.4`$) in our spectrum of SDSS J1030$`+`$0524 (Fig. 2), and no windows at all are seen in SDSS J1148$`+`$5251. In the $`5<z<5.5`$ Ly$`\alpha `$ forest, only 0.5% of the bandpass has optical depths less than 0.4. Songaila & Cowie (2002) discuss the distribution of the transmitted fraction as a function of redshift; using their analytical model, at $`z=5`$ only 1.7% of the spectrum is expected to have a transmission $`>0.65`$. This fraction drops to only 0.6% if we use the $`1\sigma `$ upper limit, $`\tau (\text{Ly}\alpha \text{@}z=5)<0.2`$.

We conclude that it is problematic to reconcile the strong peak at $`z6.03`$ in the Ly$`\beta `$ trough with the absence of a corresponding Ly$`\alpha `$ peak. Unless the Ly$`\alpha `$ forest is extraordinarily transparent at the position of the peak, we would expect a strong, easily detected peak in the Ly$`\alpha `$ GP trough.

There is, however, another possible explanation for the strong peaks in the Ly$`\beta `$ trough of SDSS J1148$`+`$5251. A careful examination of the quasar spectrum (Fig. 8) reveals a strong C IV absorption feature at 9200A ($`z=4.943`$). That suggests a second interpretation for the spikes, namely that they are Ly$`\alpha `$ emission from a foreground galaxy at $`z=4.943`$. Usually one does not see emission associated with C IV absorbers; but ordinarily such emission lines would be lost in the Ly$`\alpha `$ forest and not noticed. This particular emission line, falling as it does in the Ly$`\beta `$ trough, stands out like a sore thumb. Moreover there are multiple C IV absorption systems and multiple (putative) Ly$`\alpha `$ emission components seen, suggesting the presence of a fairly massive mass concentration.

Is it plausible that the emission in the Ly$`\beta `$ trough could be from an intervening $`z=4.9`$ galaxy? The implied Ly$`\alpha `$ luminosities are similar to those of Ly$`\alpha `$ emitters detected in the Subaru deep field at $`z=4.86`$ (Ouchi et al. 2003). Moreover, the velocity offset between the Ly$`\alpha `$ emission and the C IV absorption seen in Fig. 8 is in the range typically observed for systems that show both Ly$`\alpha `$ emission and metal line absorption ($`500`$ km/s, Shapley et al. 2002).

The main argument against the intervening emission model is that Ly$`\alpha `$ emitters are sufficiently rare on the sky that it is unlikely to find one along a random line-of-sight. But perhaps this is not a random line-of-sight. Figure 8 shows that there are both multiple Ly$`\alpha `$ emission components and multiple C IV absorption components present in the spectrum (note particularly the structure in both components of the strong C IV doublet at 9200 Å.) The two presumed Ly$`\alpha `$ lines are split by 1200 km/s, and the C IV structure extends over at least 250 km/s. It is possible that there is a massive structure (a protocluster?) along the line-of-sight that is amplifying the quasar through gravitational lensing. If that is correct, then the likelihood of finding a Ly$`\alpha `$ emitter in front of the $`z=6.37`$ quasar is in fact not so small, since the foreground object assisted in the discovery of the quasar by boosting its apparent brightness. Even lensing by a factor of two might create a strong discovery bias since the luminosity function is steep for bright quasars (Wyithe & Loeb 2002; Comerford, Haiman & Schaye 2002).

For the more familiar case where a nearby ($`z<1`$) galaxy is lensing a high-redshift quasar, the lensing magnification is limited to less than $`2`$ because stronger lensing produces multiple images with a noticeably wide splitting. But when the lens redshift is $`z5`$, as in the scenario proposed here, the Einstein ring radius is much smaller and so the splitting could not have been detected from the ground (in agreement with the imaging observations by Fan et al. 2003.) Figure 9 shows the magnification and observed image positions as a function of impact parameter for a singular isothermal sphere galaxy model with a velocity dispersion of 250 km/s. This velocity dispersion is consistent with the velocity structure seen in the strong C IV absorption system and is slightly larger than the dispersion expected for an $`L_{}`$ galaxy ($`v=220`$ km/s). In this geometry, the multiple image separation is only 0.3 arcsec, which would be difficult to detect even if the components were of equal brightness.

Note also that since there are multiple C IV absorption and Ly$`\alpha `$ emission components, it is possible that there are multiple lensing galaxies along the line-of-sight to this source. That could possibly increase the amplification factor further, although it would also lead to larger splitting and shifts. We coadded the two-dimensional spectra of SDSS J1148$`+`$5251 for the Echelle order that includes the Ly$`\beta `$ trough to search for evidence that the peaks are either broader than the quasar trace or are shifted relative to the center of the trace. Neither effect was seen. The width of the peaks agrees with the width of the rest of the order, and the central positions agree to better than 0.1″ (which is the approximate uncertainty in the trace position.)

This appears to be the most problematic aspect of interpreting these peaks as Ly$`\alpha `$ emission from a foreground lensing galaxy: we see no direct evidence any any images (including the spectrum) for lensing. That does not contradict the lensing hypothesis given the small lens scale at $`z=5`$, but neither does it provide any additional support. And there is certainly reason to think it unlikely that we would stumble across a lens with such a small Einstein ring radius, since the area on the sky that is lensed by such objects is rather small. But as we discuss below, the additional supporting evidence in favor of the hypothesis is sufficient to make a reasonably strong case.

### 5.3 Continuum emission in the troughs of SDSS J1148$`+`$5251

A puzzling aspect of the Ly$`\alpha `$ and Ly$`\beta `$ GP trough transmissions for SDSS J1148$`+`$5251 is that the Ly$`\alpha `$ residual flux is far higher than expected from the Ly$`\beta `$ flux. This is reflected in Table 2 by the fact that $`\tau (\text{Ly}\alpha )`$ inferred from the Ly$`\beta `$ trough is far higher than the optical depth actually measured in Ly$`\alpha `$. Either there is too much light in Ly$`\alpha `$ trough or there is too little in the Ly$`\beta `$ trough (which presumably would require a great deal of extra absorption in the $`z5`$ Ly$`\alpha `$ forest.) Neither explanation appears particularly palatable based on a conventional model of a quasar with IGM absorption.

An intervening Ly$`\alpha `$-emitting galaxy at $`z=5`$ offers a more natural explanation: the Ly$`\alpha `$ and Ly$`\beta `$ troughs are contaminated by continuum emission from the galaxy. We find that the observed brightness of the continuum is fully consistent with the hypothesis of an intervening galaxy, which strongly supports our suggestion. Wyithe & Loeb (2002) pointed out that a lensing galaxy could contaminate the black GP trough, although they were concerned mainly with lensing by low redshift galaxies (which is a priori considerably more likely.)

Between 8420 Å and 8900 Å ($`z=5.93`$ to 6.32), the only obvious emission in the Ly$`\alpha `$ trough is the peak at $`z=6.08`$ discussed above (Fig. 4). If we sum the light in that wavelength range, excluding the peak ($`8600<\lambda <8620`$), we get an estimate for the continuum level of $`3.9\pm 0.6\times 10^{20}\text{erg}\text{s}^1\text{cm}^2\text{Å}^1`$. We take that as an estimate of the brightness of the galaxy continuum. The AB magnitude of the flux detected in the Ly$`\alpha `$ trough is 26.4, which corresponds to an absolute AB magnitude of $`21.9`$ at a rest wavelength of 1450 Å if the light comes from a galaxy at $`z=5`$.

The observed equivalent width (EW) of the peak at 7205 Å is then 930 Å, and the rest-frame EW is 160. This EW is perfectly consistent with the EW distribution determined by Malhotra & Rhoads (2002) for Ly$`\alpha `$-emitting galaxies at $`z=4.5`$. It is slightly lower than the median EW for their sample, as is expected for a more massive galaxy (Malhotra & Rhoads, private communication).

### 5.4 The extent of the Strömgren sphere

Many authors have discussed the possibility that high-redshift quasars are selected because they are strongly lensed. Wyithe & Loeb (2002) and Comerford, Haiman, & Schaye (2002) determined how lensing might bias the quasar luminosity function. Haiman & Cen (2002; hereafter HC) considered the specific case of lensing of SDSS J1030$`+`$0524; they argue that if the quasar is embedded in a largely neutral IGM, the large extent of its H II region implies that it must in fact have a large ionizing luminosity and so cannot have been too strongly lensed. If the ionizing luminosity were too low, the IGM near the quasar would have remained neutral and so the Ly$`\alpha `$ emission line would have been much more heavily absorbed on the blue side. This argument also sets a minimum age for the quasar, since it must have been in existence long enough to ionize a substantial Strömgren sphere.

The HC result is directly relevant to the question of whether SDSS J1148$`+`$5251 could be strongly gravitationally lensed. On the blue side of Ly$`\alpha `$, emission is detectable to $`z=6.32`$ (Fig. 4), a comoving distance of 2.8 Mpc from the quasar. (Note that this distance is somewhat uncertain due to the difficulty in determining an accurate redshift for the quasar; this is discussed further below.) According to HC — see also Madau & Rees (2000), Cen & Haiman (2000) and references therein — the radius of a Strömgren sphere ionized by a high-redshift quasar as a function of time is approximately

$$R_S=7.0\text{Mpc}\left(\frac{\dot{N}_{ph}}{10^{58}\text{s}^1}\frac{t_Q}{10^7\text{yr}}\right)^{1/3}\left(\frac{1+z_Q}{7.37}\right)^1,$$
(6)
where $`\dot{N}_{ph}`$ is the ionizing photon luminosity, $`t_Q`$ is the age of the quasar, and as above we use cosmological parameters $`H_0=70`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_m=0.3`$, and $`\mathrm{\Omega }_b=0.04`$. This formula results from matching the total number of ionizing photons emitted to the number of hydrogen atoms within a sphere, and it ignores both recombinations (which reduce the radius) and the Hubble expansion (which increases it.) It uses the mean IGM density assuming that the IGM contains nearly all the baryons; if the universe is locally overdense by a factor $`D`$, the radius is reduced by $`D^{1/3}`$.

Inverting Eqn. (6) to solve for the quasar parameters, we find

$$\frac{\dot{N}_{ph}}{10^{58}\text{s}^1}\frac{t_Q}{10^7\text{yr}}=0.063\left(\frac{R_S}{2.8\text{Mpc}}\frac{1+z_Q}{7.37}\right)^3.$$
(7)
The ionizing luminosity for SDSS J1148$`+`$5251 (assuming no lensing) is rather uncertain. We have estimated it using two different quasar spectral templates. The Telfer et al. (2002) quasar template, which we have used for the Ly$`\alpha `$ and Ly$`\beta `$ continuum normalization, gives an ionizing photon luminosity of $`\dot{N}_{ph}=0.9\text{}1.3\times 10^{58}\text{s}^1`$, where the smaller value assumes the EUV spectral index for radio-loud quasars ($`\alpha _{EUV}=1.96`$) and the larger assumes the radio-quiet quasar index ($`\alpha _{EUV}=1.57`$). The Elvis et al. (1994) template, which was used by HC, gives the substantially smaller value $`\dot{N}_{ph}=0.2\times 10^{58}`$. We think the Telfer et al. spectrum is preferable but report results using both spectra below. The corresponding range of values for $`t_Q`$ is 5–$`7\times 10^5`$ yrs or $`3\times 10^6`$ yrs (for the Telfer and Elvis templates, respectively.)

Such a short lifetime for SDSS J1148$`+`$5251 would be surprising, since its high luminosity implies a large black-hole mass which should take a considerable length of time to assemble. The $`e`$-folding timescale for an accreting black hole is

$$t_{acc}=4\times 10^7\text{yr}\left(\frac{ϵ}{0.1}\right)\eta ^1,$$
(8)
where $`ϵ`$ is the radiative efficiency and $`\eta `$ is the ratio of the quasar’s luminosity to the Eddington luminosity (Haiman & Loeb 2001). This long natural timescale makes it unlikely that we will observe quasars with ages as short as a few hundred thousand years unless the radiative efficiency $`ϵ`$ is very low, allowing the quasar to accrete quietly and rapidly.

On the other hand, following HC’s argument, if the quasar is lensed then its true luminosity could be substantially smaller. With a plausible age of $`t_Q10^7`$ yrs, SDSS J1148$`+`$5251 could be easily lensed by a factor of 3 or more and still be able to ionize the IGM out to the observed distance of 2.8 Mpc.

We note that for the sizes and ages derived for this object, the Strömgren sphere is still in its early phase of rapid relativistic expansion. Even in that stage of evolution, however, Eqn. (6) accurately describes the observed size of the H II region as measured by absorption along the line-of-sight. See the Appendix for further details.

One more point that is worth making is that these discrepancies become even greater if the IGM is not completely neutral at $`z6`$. The expected value of $`\dot{N}_{ph}t_Q`$ decreases in direct proportion to the neutral fraction of the gas, so if the IGM at $`z=6`$ has $`n_{\mathrm{HI}}/n_\mathrm{H}=0.15`$ as suggested by Cen (2003), the quasar would have to be very highly magnified through lensing, in an extremely overdense region of the IGM, or extraordinarily young in order to produce an H II region as small as that we observe.

What other explanations could there be for the small H II region observed in front of SDSS J1148$`+`$5251? One possibility is that the redshift of SDSS J1148$`+`$5251 is underestimated, which would increase the observed size of the ionized region. If the quasar redshift is 6.41 instead of 6.37 (Willott, McLure & Jarvis 2003), the Strömgren radius increases to 4.7 Mpc. That would increase the $`\dot{N}_{ph}t_Q`$ product in Eqn. (7) by a factor of 4.8 to 0.31, which reduces the discrepancy with the observed luminosity and expected lifetime. However, a redshift as high as 6.41 produces a very poor match between the LBQS template and the observed spectrum (Fig. 3), so we consider a redshift that large unlikely. Matching the LBQS composite spectrum to the Mg II emission line detected at $`2\mu `$m by Willott et al. (2003) reveals that $`z=6.37`$ matches the data about as well as $`z=6.41`$; in fact, from the LBQS spectrum, which does include the offsets typically seen between the redshifts of different emission lines, a redshift of 6.41 would appear to be an upper limit for SDSS J1148$`+`$5251. It will probably require higher SNR infrared spectra of this object to resolve the question of its true redshift.

Another possibility is that clumpy structure in the IGM conspires to limit the size of the Strömgren sphere along the line of sight. Certainly there is substantial IGM structure near the quasar, as it can be seen in the narrow absorption lines on the blue edge of the Ly$`\alpha `$ emission line for wavelengths approaching the GP trough. A cloudy IGM can also modify the radiative transfer for ionizing radiation so that the Strömgren sphere is not fully ionized (Cen & Haiman 2000; Fan et al. 2002), so that gas within the ionized region contributes substantially to the observed Ly$`\alpha `$ optical depth.

However, these effects ought to be relatively small for H II regions that are still in their very early phase of relativistic expansion (see Appendix A). In that phase there is a surfeit of ionizing photons, with far more photons present than are required to ionize the neutral hydrogen. The expansion rate of the ionization front is limited by the light-travel time; recombinations are completely negligible, and even dense gas clumps are fully ionized by the copious radiation.

Stopping the ionization front with a dense cloud along the line of sight is difficult for the same reason. For example, suppose there is a dense cloud in the IGM at $`z=6.32`$. If the cloud is optically thick to Lyman continuum radiation, it will stop the Strömgren sphere expansion, thus producing an anomalously small ionized region. A cloud dense enough to stop the ionization front must have a density greater than

$$n_H=\left(\frac{\dot{N}_{ph}}{4\pi R_S^2\alpha _B\mathrm{}}\right)^{1/2}=0.12\mathrm{cm}^3\left(\frac{\dot{N}_{ph}}{10^{58}\mathrm{s}^1}\frac{\mathrm{kpc}}{\mathrm{}}\right)^{1/2}\left(\frac{R_S}{2.8\mathrm{Mpc}}\right)^1,$$
(9)
where $`\alpha _B`$ is the case B recombination coefficient and $`\mathrm{}`$ is the line-of-sight distance through the cloud. If we assume the cloud is spherical, its mass is

$$M=2\times 10^6\mathrm{M}_{\mathrm{}}\left(\frac{\mathrm{}}{\mathrm{kpc}}\right)^{5/2}\left(\frac{\dot{N}_{ph}}{10^{58}\mathrm{s}}\right)^{1/2}\left(\frac{R_S}{2.8\mathrm{Mpc}}\right)^1.$$
(10)
The cloud would have to be either dense or quite massive to stop the quasar’s ionization front.

Another possible explanation for the small size of the ionized region is that the edge of the observed emission at $`z=6.32`$ does not indicate the edge of the ionized gas. Only a trace of neutral hydrogen is needed to generate a large optical depth to Ly$`\alpha `$ (eq. 2). Consequently, if the ionization fraction in the H II region drops below $`10^4`$, the IGM becomes opaque. This can occur if radiation transfer in the clumpy IGM reduces the ionizing flux seen by the gas. This also appears unlikely in the Strömgren sphere’s relativistic expansion phase, but more detailed calculations are required to confirm that $`z=6.32`$ really marks the position of the quasar’s ionization front.

### 5.5 Could SDSS J1030$`+`$0524 also be lensed?

One might wonder whether the arguments presented in favor of a foreground object lensing SDSS J1148$`+`$5251 could also be applied to SDSS J1030$`+`$0524. SDSS J1030$`+`$0524 does show some absorption lines from intervening metal line systems (Figs. 1 and 2), but none are as deep as the strongest C IV absorption line system seen in SDSS J1148$`+`$5251, nor do we see any evidence for associated Ly$`\alpha `$ emission. The absorption spectrum of SDSS J1030$`+`$0524 will be discussed further by Madau & Bolte (in preparation.)

There is no sign that continuum emission from an intervening galaxy contaminates the Ly$`\alpha `$ or Ly$`\beta `$ GP troughs in SDSS J1030$`+`$0524. Indeed, those troughs are very black, with residual levels even below those detected in SDSS J1148$`+`$5251 (see Table 2). From the flux limit in the Ly$`\alpha `$ trough, we place a limit of $`m_{\mathrm{AB}}>27.4`$ on any intervening galaxy. We can rule out the presence of a lensing galaxy at $`z<3.5`$ with an absolute AB magnitude brighter than $`20`$ (at a rest wavelength of 2000 Å).

Finally, Haiman & Cen (2002) discussed the limits that can be placed on any lensing of SDSS J1030$`+`$0524 from the size of its Strömgren sphere. The quasar has an observed $`R_S=4.5`$ Mpc at $`z=6.28`$, so we find $`(\dot{N}_{ph}/10^{58})(t_Q/10^7)=0.25`$. The ionizing photon luminosity is $`0.13\times 10^{58}\text{s}^1`$ or 0.61–$`0.87\times 10^{58}\text{s}^1`$ (from the Elvis and Telfer templates), leading to age estimates of $`2\times 10^7`$ yrs or 3–$`4\times 10^6`$ yrs. These age estimates appear plausible, although the higher ionizing luminosities using the Telfer et al. template do leave more room for lensing than was found by HC.

### 5.6 Summary

Our proposal gives a consistent picture for SDSS J1148$`+`$5251: the GP troughs of both Ly$`\alpha `$ and Ly$`\beta `$ are likely to be contaminated with light from an intervening galaxy at $`z5`$. The existence of this galaxy is supported by the detection of Ly$`\alpha `$ emission, C IV absorption, and continuum emission. There could be several such intervening systems, since both the emission and absorption are complex with multiple components spread over $`1000`$ km/s. The intervening system amplifies the quasar’s light via gravitational lensing, thereby enhancing the likelihood of discovery for the quasar. The possibility of strong lensing is supported by the relatively small H II region created by the quasar’s ionizing radiation, which indicates that the quasar is probably considerably less luminous than its apparent brightness would indicate.

There may be some way to reconcile the observations of SDSS J1148$`+`$5251 presented here with a model based purely on IGM absorption. No such explanation is obvious to us, though there are so many excellent theoreticians who are interested in this problem that there will doubtless be many good ideas proposed shortly after the publication of these results! From the observational side, clearly an HST image of SDSS J1148$`+`$5251 would be very interesting. It could reveal the presence of multiple images (or set strong limits on their absence.) An HST image taken using a narrow-band filter in the Ly$`\beta `$ trough would directly determine whether the peaks seen at 7205 and 7236 Å are extended and offset (as expected if they are Ly$`\alpha `$ emission from an intervening galaxy) or are point-like and coincident with the quasar (as expected if they are the result of “leaks” in the IGM absorption.)

If the intervening galaxy suggestion turns out to be correct, then our observations of SDSS J1030$`+`$0524 remain the single best measurement of the transmission of the IGM at $`z>6`$. There is no light detected in the Ly$`\alpha `$ or Ly$`\beta `$ troughs for this object at a very low level. On the other hand, if SDSS J1148$`+`$5251 is shown not to have an intervening system contaminating its GP troughs, then the IGM toward it is in fact quite transparent, with both numerous high-ionization holes and a significant transmission across the whole Ly$`\alpha `$ and Ly$`\beta `$ troughs. The resolution of these questions will have to wait for the discovery of additional $`z>6`$ quasars.

Many thanks to Michael Bolte and Piero Madau for sharing their ESI exposures of SDSS J1030$`+`$0524 with us. We appreciate a helpful conversation with James Rhoads and Sangeeta Malhotra. RLW thanks Holland Ford and the Johns Hopkins University Astronomy & Physics Department for their hospitality during a sabbatical, when much of this work was carried out, and appreciates continuing support from the Space Telescope Science Institute. RHB acknowledges support from NSF grant AST-00-98355 and the Institute of Geophysics and Planetary Physics (operated under the auspices of the U. S. Department of Energy by Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48). XF acknowledges support from the University of Arizona and an Alfred P. Sloan Research Fellowship. MAS acknowledges support from NSF grant AST-00-71091.

## Appendix A Appendix: Early Evolution of a Strömgren Sphere

The early evolution of the ionized region around a quasar is marked by a period of very rapid expansion, with the ionizing front moving out at nearly the speed of light. In that case the finite light travel time across the Strömgren sphere cannot be ignored. However, the size of the sphere as inferred from observations of absorption along the line-of-sight to the quasar turns out to have exactly the same evolution with time as one derives under the assumption of an infinite speed of light. This appendix derives the evolution of the size both as seen in the frame of the ionizing source and as seen through the line-of-sight absorption.

We assume that the mean free path for ionizing photons in the neutral gas is very short and that recombinations and the expansion of the universe can be ignored. These are good approximations: the mean free path for 1 Rydberg photons is only 1 kpc at $`z=6`$ and declines as $`(1+z)^3`$. Cosmological simulations indicate that the gas clumping factor $`C=n_H^2/n_H^2`$ is likely to be well below the value $`C100`$ required to make recombinations a significant factor (Gnedin & Ostriker 1997; Madau & Rees 2000; Cen & Haiman 2000).

With these assumptions, ionizing photons propagate freely to the edge of the ionized region where they are immediately absorbed by the neutral gas. Figure 10 shows the relevant geometry. The photons that would have traveled beyond the Strömgren radius $`R_S`$ (the shaded region) have been absorbed by neutral hydrogen atoms within $`R_S`$. The size of the H II region is determined by a balance between the number of photons in the outer shell and the number of atoms within the ionized sphere:

$$\dot{N}_{ph}\left(t\frac{R_S}{c}\right)=\frac{4\pi R_S^3n_H}{3},$$
(A1)
where $`\dot{N}_{ph}`$ is the ionizing photon luminosity and $`t`$ is the age of the source. This cubic equation can be analytically solved to give

$$R_S(t)=\frac{3y}{y^2+y+1}ct,$$
(A2)
where

$$y=\left(\frac{3^{3/2}}{2}\frac{t}{t_c}+\sqrt{\frac{27}{4}\left(\frac{t}{t_c}\right)^2+1}\right)^{2/3}$$
(A3)
and the time $`t_c`$, which roughly marks the end of the relativistic expansion period, is

$$t_c=\left(\frac{3\dot{N}_{ph}}{4\pi n_Hc^3}\right)^{1/2}=1.2\times 10^7\text{yr}\left(\frac{\dot{N}_{ph}}{10^{57}\text{s}^1}\right)^{1/2}\left(\frac{1+z_Q}{7}\right)^{3/2}\left(\frac{n_{\mathrm{HI}}}{n_\mathrm{H}}\frac{\mathrm{\Omega }_bh^2}{0.0196}\right)^{1/2},$$
(A4)
where $`z_Q`$ is the quasar redshift. Here we have explicitly included the dependence on the neutral fraction $`n_{\mathrm{HI}}/n_\mathrm{H}`$ and the cosmological parameters.

The expansion law in Eqn. (A2) is shown in Figure 11. After an initial period where the ionization front moves out at nearly the speed of light, the expansion slows and approaches the standard $`t^{1/3}`$ law (Madau & Rees 2000, Cen & Haiman 2000).

This is not, however, the expansion law that is observed via absorption along the line-of-sight. When we observe Ly$`\alpha `$ absorption from the neutral gas just beyond the edge of the H II region, the Ly$`\alpha `$ photons being detected were emitted only a short time after the first ionizing photons. Specifically, if the Ly$`\alpha `$ photons cross the Strömgren boundary at a radius $`R_S`$, they were emitted when the age of the quasar was only $`t^{}=tR_S/c`$ (which is the time for light to cross the outer, shaded shell in Fig. 10). The observed Strömgren evolution, as inferred from absorption, is not $`R_S(t)`$ versus $`t`$ but rather $`R_S(t)`$ versus $`t^{}`$. With some algebraic manipulation of the analytical expressions for $`R_S`$ and $`t^{}`$, we can express $`R_S`$ directly in terms of the observed (apparent) time $`t^{}`$ as

$$R_S=R_c\left(\frac{t^{}}{t_c}\right)^{1/3},$$
(A5)
where $`R_c=ct_c`$.

This is a surprising result (at least, it surprised us!) The expansion law for the observed radius is exactly the same as the expansion law derived if one completely ignores light-travel time effects (eq. 6). In the frame of the quasar the ionization front’s expansion velocity is limited by the speed of light, but for the observer there is an initial period of superluminal expansion. The speed of light drops out of the solution along the line-of-sight because the delay required to allow light to travel from the source to the edge of the H II region is exactly compensated by the speedup that results from that edge being closer to the observer, allowing light originating there to reach the observer sooner than light from the quasar core.

Note that the size of the observed H II region in Eqn. (A5) may be larger than the ionized region that exists at the end of the quasar’s lifetime. For example, suppose a quasar shines for a period $`t_c`$ and then shuts down. According to Eqn. (A2), at $`t_c`$ the radius of the H II region is $`0.7R_c`$. But the observer measures a radius $`R_c`$ (eq. A5). What is happening in this case is that the H II region continues to expand after the quasar shuts off until it reaches a maximum size $`R_c`$ at $`t=2t_c`$. At that point the inner sphere in Figure 10 is completely empty of ionizing photons and so the expansion stops. The observer is detecting Ly$`\alpha `$ photons that were emitted just at the end of the quasar’s lifetime, but by the time those photons cross the ionization front it has expanded to its maximum size.

Our conclusion is that when all the light travel time effects are taken into account, Eqn. (A5) is the correct formula for computing the age of the quasar given the luminosity. Ages determined using that formula are, in fact, the age of the quasar at the time it emitted the photons being detected, and so such ages are direct measurements of the minimum quasar lifetime even during the superluminal expansion phase.


format: text

text:
# Unitary Gate Synthesis for Continuous Variable Systems

## I Introduction

One of the central problems of the quantum information theory is to establish what resources are sufficient for universal quantum computation. In this context, the question whether a given Hamiltonian $`H`$ can simulate  another one has attracted considerable attention recently Dur01 ; Bennett02 ; Dodd02 ; Nielsen02 ; Wocjan02 ; Vidal02 . In its simplest form, this problem may be formulated as follows. Consider two parties, traditionally referred to as Alice and Bob, possessing a single qubit each. The interaction between those two qubits is governed by a fixed Hamiltonian $`H`$, that is determined by the physical properties of the systems that represent the qubits. In addition to the interaction $`H`$, Alice and Bob may attach (local) ancillas to their qubits and perform arbitrary local unitary operations on their subsystems. It is usually assumed that these local operations are very fast compared to the evolution induced by the Hamiltonian $`H`$. The task for Alice and Bob is to simulate the evolution due to different Hamiltonian $`H^{}`$. Two kinds of simulations should be distinguished. The infinitesimal time simulation Dur01 ; Bennett02 ; Dodd02 ; Nielsen02 ; Wocjan02 ; Vidal02  consists of simulating the action of the Hamiltonian $`H^{}`$ for an infinitesimally short time interval $`\mathrm{\Delta }t`$. The gate synthesis Khaneja01 ; Vidal02b ; Hammerer02 ; Masanes02 ; Bremner02 ; Zhang02  requires the implementation of the unitary transformation $`U^{}=\mathrm{exp}(iH^{}t)`$ for some finite time $`t`$.

It turns out that in the two-party setting all nonlocal Hamiltonians $`H`$ are qualitatively equivalent. Given enough time $`\tau `$, Alice and Bob can, with the help of local ancillas, simulate the evolution $`\mathrm{exp}(iH^{}t)`$ for any $`H^{}`$ Bennett02 . The central question, then, is what is the optimal simulation. The latter may be defined as a simulation that requires the shortest time. For the two-qubit case, this problem has been completely solved and the optimal protocols for Hamiltonian Bennett02 ; Vidal02  and unitary gate Vidal02b  simulations have been determined. The situation becomes much more complicated for higher dimensional systems and for higher number of involved parties. Simulation protocols suggested for these generic settings are rather involved and it is not known which protocols are optimal.

It should be stressed that most of the work focused on discrete variable systems: qubits or, more generally, qudits. Recently, however, Kraus et al. extended the notion of Hamiltonian simulation to continuous variable systems Kraus02 . They assumed that Alice and Bob possess a single-mode system each and these two modes are coupled via quadratic Hamiltonian (we assume $`\mathrm{}=1`$ throughout this paper):

$$H=c_{11}x_Ax_B+c_{12}x_Ap_B+c_{21}p_Ax_B+c_{22}p_Ap_B,$$
(1)
where $`x_j`$ and $`p_j`$ are two conjugate quadratures of the $`j`$th mode. Kraus et al. showed that almost every Hamiltonian (1) is capable to simulate any other Hamiltonian of the form (1) provided that Alice and Bob can apply fast local phase shifts described by single-mode Hamiltonians $`H_A=x_A^2+p_A^2`$ and $`H_B=x_B^2+p_B^2`$.

These results are interesting both from the theoretical and experimental points of view. In particular, the off-resonant interaction of light with the collective atomic spin Kuzmich98 ; Kuzmich00 ; Kuzmich00theory ; Duan00 ; Julsgaard01 ; Schori02 ; Kuzmichbook  can be described by the effective unitary transformation

$$U=\mathrm{exp}(itH_{AL}),$$
(2)
where the Hamiltonian

$$H_{AL}=\kappa x_Ap_B$$
(3)
is a special instance of (1). The typical geometry of the experiments is such that a light with strong coherent field polarized along the $`x`$ axis propagates along the $`z`$ axis through the atomic sample, whose spin is also polarized along the $`x`$ axis. The $`x`$ and $`p`$ quadratures are defined as the properly normalized $`y`$\- and $`z`$-components of the collective spin operators describing the polarization state of light and atomic ensemble, respectively Kuzmich98 ; Kuzmich00theory ; Duan00 . In recent beautiful experiments it was demonstrated that the interaction (2) can be employed to squeeze the atomic spin Kuzmich00 , entangle two distant atomic ensembles Julsgaard01 , and transfer the quantum state of light into the atoms Schori02 . Schemes for teleportation and swapping of the quantum state of collective atomic spin have been suggested Kuzmich00theory ; Duan00 . These experiments and proposals in fact rely on the quantum non-demolition (QND) measurement of the atomic quadrature, possibly accompanied by a suitable feedback.

As showed by Kraus et al., the Hamiltonian (3) can simulate any Hamiltonian (1). In particular, $`H_{AL}`$ can be used to implement a beam splitter and a two-mode squeezer. This is very appealing since it suggests that, for instance, the storage of the quantum state of light into atoms and the subsequent readout of the quantum memory – the transfer of quantum state of atoms onto light – can be implemented in a unitary way if the Hamiltonian (3) is used to simulate a beam splitter.

However, currently there are technical difficulties that will complicate the actual practical realization of this procedure. The effective unitary transformation (2) describes the modification of the polarization state of the light pulse after the passage through the atomic sample due to its coupling with atoms. This means that the Hamiltonian simulation requires several passages of the light pulse through the atomic sample (c.f. the detailed description of the simulation protocol in Sec. II). In currently envisaged experiments, the pulse width must be at least $`1`$ $`\mu \mathrm{s}`$ Polzikprivate  which corresponds to the length $`300`$ m. This long pulse would have to be stored somewhere (e.g. in an optical fiber) until its tail leaves the atomic sample. Only then can the pulse (with properly applied phase shifts) be fed to the atomic sample again.

These practical considerations imply that the approach relying on the infinitesimal time simulation is not very convenient from the experimental point of view. It is possible to simulate a gate by concatenating a large sequence of short-time Hamiltonian simulations but this would require a large number of manipulations and passages of the light pulse through the sample. Since, in practice, every round of the gate synthesis procedure is necessarily accompanied by some losses and other errors, the accumulation of the errors would negatively influence the simulation.

In this paper, we show how to implement several important two-mode interactions with the Hamiltonian (1) such that the number of the applications of the Hamiltonian (1) is minimized. We demonstrate that only three sequences of evolution governed by Hamiltonian $`H`$, interspaced by (fast) local phase shifts on both subsystems, suffice to implement a two-mode squeezer and a beam splitter. For the specific Hamiltonian (3) we also design a single-mode squeezing gate that involves four evolution steps. These results illustrate that several important quantum information processing tasks such as entangling the light and collective atomic spin, or a transfer of the quantum state of light into atomic clouds and vice versa, can be carried out with a small number of repeated passages of the light pulse through the atomic sample.

This paper is structured as follows. In Sec. II we introduce the notation, the canonical form of the interaction Hamiltonian (1) and we describe the gate simulation protocol. In Sec. III, we consider the simple interaction Hamiltonian (3) and we show how to implement the two-mode squeezing operation, beam splitter transformation and also single-mode squeezing as a sequence of three (or four) intervals of evolution governed by Hamiltonian (3) combined with local phase shift operations. In Sec. IV we extend this analysis to the generic interaction Hamiltonians (1). Finally, Sec. V contains conclusions.

## II Description of the simulation protocol

In this section we describe the simulation protocol. The gate synthesis consists of a sequence of $`N`$ intervals of evolution governed by the Hamiltonian $`H`$ followed by local unitary phase shift transformations. The resulting unitary gate $`G`$ is given by

$$G=V_N^{}e^{iHt_N}V_N\mathrm{}V_2^{}e^{iHt_2}V_2V_1^{}e^{iHt_1}V_1.$$
(4)
The local phase shift operation applied to modes A and B reads

$$V_j=e^{i\varphi _{jA}a^{}a}e^{i\varphi _{jB}b^{}b},$$
(5)
where $`a`$ and $`b`$ are the annihilation operators of modes $`A`$ and $`B`$, respectively. Note that $`V_{j+1}^{}V_j`$ is still of the form (5), with $`\varphi _A=\varphi _{jA}\varphi _{j+1,A}`$ and $`\varphi _B=\varphi _{jB}\varphi _{j+1,B}`$. With the help of the useful identity

$$U^{}\mathrm{exp}(iHt)U=\mathrm{exp}(iU^{}HUt)$$
(6)
we can rewrite Eq. (4) as

$$G=e^{iH_Nt_N}\mathrm{}e^{iH_2t_2}e^{iH_1t_1},$$
(7)
where $`H_j=V_j^{}HV_j`$.

The Hamiltonian (1) is characterized by four parameters. However, by means of local rotations, we can always transform this Hamiltonian to a diagonal form

$$H_c=c_1x_Ap_B+c_2p_Ax_B,$$
(8)
where $`c_1=\sigma _1`$ and $`c_2=\sigma _2det[C]/|det[C]|`$, and $`\sigma _1`$ and $`\sigma _2`$ are the singular values of the matrix $`C`$ defined as $`(C)_{ij}=c_{ij}`$ Kraus02 . In close analogy to the qubit case Bennett02 , we may refer to $`H_c`$ as the canonical form of $`H`$. Mathematically, we have

$$\mathrm{exp}(iH_ct)=\stackrel{~}{V}^{}\mathrm{exp}(iHt)\stackrel{~}{V},$$
(9)
where $`\stackrel{~}{V}`$ is a local rotation (5). This shows that, without loss of generality, we may assume that $`H`$ has the canonical form (8). In particular, it follows that $`H`$ is able to simulate an arbitrary $`H^{}`$ (1) if and only if $`H_c`$ is able to simulate an arbitrary canonical Hamiltonian (8).

In Eq. (5) the phase shifts $`\varphi _j`$ may be arbitrary. In what follows, we focus on the phase shifts that preserve the canonical form of $`H`$. There are four inequivalent possibilities:
(a) $`\varphi _A=0`$, $`\varphi _B=0`$,

$$H_1=c_1x_Ap_B+c_2p_Ax_B.$$
(10)
(b) $`\varphi _A=\pi /2`$, $`\varphi _B=3\pi /2`$,

$$H_2=c_2x_Ap_B+c_1p_Ax_B.$$
(11)
(c) $`\varphi _A=\pi `$, $`\varphi _B=0`$,

$$H_3=c_1x_Ap_Bc_2p_Ax_B.$$
(12)
(d) $`\varphi _A=\pi /2`$, $`\varphi _B=\pi /2`$,

$$H_4=c_2x_Ap_Bc_1p_Ax_B.$$
(13)
From the structure of these Hamiltonians we can deduce that two different noncommuting canonical Hamiltonians $`H_1`$ and $`H_2`$ are available. Furthermore, we can see that $`H_3=H_1`$ and $`H_4=H_2`$, hence we can implement any transformations of the form $`\mathrm{exp}(iH_1t)`$ and $`\mathrm{exp}(iH_2t)`$ where $`t`$ is an arbitrary real number, positive or negative. The two specific cases $`c_1=c_2`$ and $`c_1=c_2`$ when $`H_1=\pm H_2`$ and the simulation is not possible correspond to the Hamiltonians of a two-mode squeezer and a beam splitter, respectively.

## III XP coupling

Having established the notation and described the simulation protocol, we may proceed to the unitary gate synthesis. Namely, we would like to decompose the unitary transformation $`G`$ that we want to simulate into a sequence of unitary evolutions governed by Hamiltonians $`H_1`$and $`H_2`$ that were defined in the previous section,

$$G=e^{iH_2t_N}e^{iH_1t_{N1}}\mathrm{}e^{iH_2t_2}e^{iH_1t_1}.$$
(14)
We are particularly interested in the simulations that involve the lowest possible number of steps $`N`$, because such simulations require low number of local manipulations in the eventual experimental implementation.

We note here that Eq. (14) is an example of a decomposition of a group element into a product of $`N`$ other group elements. In the present case, the underlying group is the symplectic group $`\text{Sp}(4,R)`$ of all linear canonical transformations of the quadratures of the two modes $`A`$ and $`B`$ Simon94 ; Arvind95 . It is worth mentioning here that the related problem of a decomposition of the symplectic transformation into a sequence of simple evolutions associated with the common passive and active linear optical elements has been studied recently. Braunstein has shown that any $`N`$ mode symplectic transformation can be implemented as a sequence of an $`N`$-mode passive linear interferometer followed by $`N`$ single-mode squeezers and another passive interferometer — the so-called Bloch-Messiah decomposition Braunstein99 . The decompositions of this kind have also been applied to investigate the properties of nonlinear optical couplers Fiurasek00 ; Rehacek02 .

In this section, we shall consider the simplest and also the experimentally relevant coupling between the two systems described by the interaction Hamiltonian (3). Without loss of generality, we may assume that the coupling constant is equal to unity, hence the two relevant Hamiltonians read

$$H_1=x_Ap_B,H_2=p_Ax_B.$$
(15)
All the canonical Hamiltonians (8) have the important property that the $`x`$ and $`p`$ quadratures are not mutually coupled when we write down the Heisenberg equations of motion for $`x_j`$ and $`p_j`$. This means that the evolution of the operators $`𝒙=(x_A,x_B)^T`$ and $`𝒑=(p_A,p_B)^T`$ is governed by the following linear canonical transformations:

$$𝒙_{\mathrm{out}}=𝑺𝒙_{\mathrm{in}},𝒑_{\mathrm{out}}=𝑹𝒑_{\mathrm{in}}.$$
(16)
This decoupling of $`x`$ and $`p`$ quadratures greatly simplifies the analysis. The transformation (16) must preserve the canonical commutation relations $`[x_j,p_k]=i\delta _{jk}`$. From these conditions we can express the matrix $`𝑹`$ in terms of $`𝑺`$,

$$𝑹=(𝑺^1)^T,$$
(17)
hence the evolution of $`p`$ quadratures is uniquely determined by the evolution of the $`x`$ quadratures.

Our task is to implement two-mode unitary gates (symplectic transformations) as a sequence of a small number of unitary transformations generated by the Hamiltonians (15). The matrices $`𝑺_1`$ and $`𝑺_2`$ associated with the unitary evolutions $`U_1=\mathrm{exp}(iH_1t)`$ and $`U_2=\mathrm{exp}(iH_2t)`$ read

$$𝑺_1(t)=\left(\begin{array}{cc}1& 0\\ t& 1\end{array}\right),𝑺_2(t)=\left(\begin{array}{cc}1& t\\ 0& 1\end{array}\right).$$
(18)
The factorization (14) can be rewritten in terms of the matrices $`S_j`$ as follows,

$$𝑺=𝑺_2(t_N)𝑺_1(t_{N1})\mathrm{}𝑺_2(t_2)𝑺_1(t_1),$$
(19)
where $`𝑺`$ is the matrix associated with the gate $`G`$. Since $`det𝑺_1=det𝑺_2=1`$, we are restricted to a three parametric subgroup of transformations $`𝑺`$ such that $`det𝑺=1`$. In what follows, we will discuss the implementation of three important gates: a beam splitter, a two-mode squeezer and a single-mode squeezer.

### III.1 Beam splitter

The beam splitter operation is described by the matrix

$$𝑺_{BS}(\varphi )=\left(\begin{array}{cc}\mathrm{cos}\varphi & \mathrm{sin}\varphi \\ \mathrm{sin}\varphi & \mathrm{cos}\varphi \end{array}\right).$$
(20)
We show that this transformation can be implemented as a sequence of three evolutions (18),

$$𝑺_{BS}(\varphi )=𝑺_1(\gamma )𝑺_2(\beta )𝑺_1(\alpha ).$$
(21)
The explicit multiplication yields

$$𝑺_{BS}(\varphi )=\left(\begin{array}{cc}1+\alpha \beta & \beta \\ \alpha +\gamma (1+\alpha \beta )& 1+\gamma \beta \end{array}\right).$$
(22)
If we compare the elements of the matrices on left- and right-hand sides of Eq. (22), we obtain a set of equations for the parameters $`\alpha `$, $`\beta `$ and $`\gamma `$ whose solution yields:

$$\alpha =\mathrm{tan}\frac{\varphi }{2},\beta =\mathrm{sin}\varphi ,\gamma =\alpha .$$
(23)
The parameters $`\alpha `$ and $`\beta `$ can be negative but this is not an obstacle as explained in the previous section since we can change the sign of the Hamiltonian $`H_1`$ or $`H_2`$ by $`\pi `$ rotation of one of the systems. Two cases of particular importance are (i) the balanced beam splitter ($`\varphi =\pi /4`$), that requires $`\alpha =1\sqrt{2}`$ and $`\beta =\sqrt{2}/2`$, and (ii) the swap ($`\varphi =\pi /2`$) that exchanges the quantum states of the two systems, $`\alpha =1`$, $`\beta =1`$.

The swap gate is closely related to the two-step protocol for mapping the state of collective atomic spin on light that was suggested by Kuzmich and Polzik Kuzmichbook . In fact, their two-step protocol can be obtained by simply removing the last step of the present three-step swap gate. The removal of the third step means that the mapping adds some noise and the procedure is thus only approximate. A possible way how to improve its performance is to use squeezed light. For details, see Kuzmichbook .

### III.2 Two-mode squeezer

Let us now turn our attention to the two-mode squeezer, described by the following matrix,

$$𝑺_{\mathrm{TMS}}(r)=\left(\begin{array}{cc}\mathrm{cosh}r& \mathrm{sinh}r\\ \mathrm{sinh}r& \mathrm{cosh}r\end{array}\right).$$
(24)
Similarly as in the case of the beam splitter, we attempt to implement this transformation as a sequence of three evolutions, c.f. Eq. (21). By comparison of the right-hand side of Eq. (22) with the matrix (24) we again obtain a system of nonlinear equations for the parameters $`\alpha `$, $`\beta `$, and $`\gamma `$ having the solution

$$\alpha =\mathrm{tanh}\frac{r}{2},\beta =\mathrm{sinh}r,\gamma =\alpha .$$
(25)
Note the similarity with the results for the beam splitter, the difference is essentially that the goniometric functions have been replaced by their hyperbolic counterparts. The parameters are finite for any finite $`r`$. However, $`\beta `$ grows exponentially with $`r`$ and for large $`r`$ we have $`\beta e^r`$. On the other hand, for small $`r`$ we get $`\beta r`$. This implies that we may reduce the synthesis time if we implement the two-mode squeezing transformation as a sequence of $`n`$ two-mode squeezers with $`r^{}=r/n`$. The reduction of the time is achieved at the expense of a higher number of steps of the gate synthesis protocol. For modest values of squeezing $`r`$, the sequence of three evolutions is advantageous since it involves the minimum necessary number of manipulations of the systems $`A`$ and $`B`$.

### III.3 Single-mode squeezer

After dealing with two-mode gates, let us now focus on the single-mode gate, namely the single-mode squeezer. Since $`det𝑺=1`$, the squeezing of quadrature $`x_A`$ will necessarily by accompanied by anti-squeezing of $`x_B`$, and vice versa,

$$𝑺_{\mathrm{SMS}}(r)=\left(\begin{array}{cc}e^r& 0\\ 0& e^r\end{array}\right).$$
(26)
It turns out that in this case a sequence of three transformations (18) is insufficient and we must consider a sequence of four basic evolutions:

$$𝑺_{\mathrm{SMS}}(r)=𝑺_2(\delta )𝑺_1(\gamma )𝑺_2(\beta )𝑺_1(\alpha ).$$
(27)
We proceed as before and derive equations for the four parameters appearing in (27):

$`e^r`$ $`=`$ $`1+\alpha \beta +\delta (\alpha +\gamma +\alpha \beta \gamma ),`$
$`e^r`$ $`=`$ $`1+\gamma \beta ,`$
$`0`$ $`=`$ $`\beta +\delta (1+\gamma \beta ),`$
$`0`$ $`=`$ $`\alpha +\gamma (1+\alpha \beta ).`$
This system of equations has a one-parametric class of solutions, given by

$$\beta =\frac{e^r1}{\alpha },\gamma =\alpha e^r,\delta =\frac{e^r(1e^r)}{\alpha },$$
(28)
and $`\alpha `$ is arbitrary but nonzero. We may choose the optimal value of $`\alpha `$ that minimizes the total time $`T`$ needed for the implementation of the squeezing operation,

$$T=|\alpha |+|\beta |+|\gamma |+|\delta |.$$
(29)
Assuming that $`r>0`$ we obtain by solving $`dT/d\alpha =0`$ the optimal value,

$$\alpha =\sqrt{\frac{e^{2r}1}{1+e^r}}.$$
(30)
In the limit of small $`r`$, all four parameters $`\alpha `$, $`\beta `$, $`\gamma `$, and $`\delta `$ are proportional to $`\sqrt{r}`$. This stems from the fact that the single-mode squeezing Hamiltonian $`H_{SMS}=x_Ap_A`$ cannot be obtained as a linear combination of the two-mode Hamiltonians $`H_1`$ and $`H_2`$ and only the terms of the order of $`O(t^2)`$ or higher in (14) may give rise to the contribution proportional to $`H_{SMS}`$.

## IV Generic quadratic coupling

In this section we shall assume that the interaction Hamiltonian has the generic canonical form. Although the mathematical analysis will be more involved we shall still be able to derive analytical formulas for the interaction times characterizing the gate synthesis. Without loss of generality, we may assume that $`c_1=1`$ in (8). We have to distinguish two classes of Hamiltonians giving rise to qualitatively different evolutions of the quadratures in the Heisenberg picture. For $`c_2>0`$ the dynamics resembles an amplifier while for $`c_2<0`$ we obtain oscillatory dynamics reminiscent that of a beam splitter. We shall discuss these two cases separately.

### IV.1 Amplifier-like Hamiltonians

Suppose first that $`c_2>0`$ and introduce a more convenient notation $`c_2=s^2`$, $`s>0`$, hence

$$H_1=x_Ap_B+s^2p_Ax_B,H_2=s^2x_Ap_B+p_Ax_B.$$
(31)
It is an easy exercise to derive matrices $`𝑺_1`$ and $`𝑺_2`$ corresponding to the unitary evolutions governed by $`H_1`$ and $`H_2`$, respectively,

$$𝑺_1^+(t)=\left(\begin{array}{cc}\mathrm{cosh}(st)& s\mathrm{sinh}(st)\\ \frac{1}{s}\mathrm{sinh}(st)& \mathrm{cosh}(st)\end{array}\right),$$
(32)
$$𝑺_2^+(t)=\left(\begin{array}{cc}\mathrm{cosh}(st)& \frac{1}{s}\mathrm{sinh}(st)\\ s\mathrm{sinh}(st)& \mathrm{cosh}(st)\end{array}\right).$$
(33)
In what follows we will focus on the implementation of the beam splitter and two-mode squeezer transformations. We have seen in the previous section that these transformations could be implemented as a sequence of three basic evolutions governed by Hamiltonians $`H_1`$ or $`H_2`$ that were of the form (15). Moreover, there was an inherent symmetry in this gate synthesis; we have found that $`\gamma =\alpha `$. It turns out that these basic symmetry properties remain valid also for the generic Hamiltonians (31), and we can thus decompose the two-mode squeezing transformation as

$$𝑺_{TMS}(r)=𝑺_1^+(\alpha /s)𝑺_2^+(\beta /s)𝑺_1^+(\alpha /s).$$
(34)
Here the parameters $`\alpha =st_1`$ and $`\beta =st_2`$ are the rescaled interaction times. The nonlinear equations for the parameters $`\alpha `$ and $`\beta `$ are much more complicated than before. Nevertheless, analytical results can be obtained. From the condition $`S_{12}=S_{21}`$ we get

$$\mathrm{tanh}\beta =\frac{2\mathrm{sinh}\alpha \mathrm{cosh}\alpha }{\mathrm{cosh}^2\alpha (s^2+1+s^2)\mathrm{sinh}^2\alpha }.$$
(35)
The condition $`S_{11}=S_{22}`$ is satisfied due to the symmetry ($`\gamma =\alpha `$), and the parameter $`\alpha `$ can be determined from the last independent equation $`S_{12}/S_{11}=\mathrm{tanh}r`$, which yields

$`{\displaystyle \frac{2y(s+s^1)}{1+y^2(s^2+1+s^2)}}=\mathrm{tanh}r,`$
where $`y=\mathrm{tanh}\alpha `$. This is a quadratic equation for $`y`$ whose solution reads

$$\mathrm{tanh}\alpha =\frac{s+s^1\sqrt{1+(s^2+1+s^2)(\mathrm{cosh}r)^2}}{(s^2+1+s^2)\mathrm{tanh}r}.$$
(36)
We have selected the root that yields the correct limit $`\mathrm{tanh}\alpha =0`$ when $`r0`$. In the opposite limit $`r\mathrm{}`$ we obtain

$$\mathrm{tanh}\alpha _{\mathrm{}}=\frac{1}{s^1+1+s}.$$
(37)
On inserting this back into Eq. (35) we find that

$$\underset{r\mathrm{}}{lim}\mathrm{tanh}\beta =1.$$
(38)
It is easy to check that the equations for $`\alpha `$ and $`\beta `$ have finite solutions for any finite $`r`$. In the limit $`r\mathrm{}`$, $`\alpha `$ approaches a finite asymptotic value, cf. Eq. (37), while $`\beta `$ grows to infinity.

Suppose now that we want to implement the beam splitter transformation (20). The calculations of the parameters $`\alpha `$ and $`\beta `$ parallel those for the two-mode squeezer. From the condition $`S_{12}=S_{21}`$ we express $`\beta `$ in terms of $`\alpha `$:

$$\mathrm{tanh}\beta =\frac{2\mathrm{sinh}\alpha \mathrm{cosh}\alpha }{\mathrm{cosh}^2\alpha +(s^21+s^2)\mathrm{sinh}^2\alpha }.$$
(39)
Since $`(s^21+s^2)1`$, it follows that $`|\mathrm{tanh}\beta ||\mathrm{tanh}(2\alpha )|`$. From the condition $`S_{12}/S_{11}=\mathrm{tan}\varphi `$ we obtain quadratic equation for $`\mathrm{tanh}\alpha `$ leading to,

$$\mathrm{tanh}\alpha =\frac{(s^1s)\mathrm{sign}(s^1s)\sqrt{\frac{s^21+s^2}{\mathrm{cos}^2\varphi }1}}{(s^21+s^2)\mathrm{tan}\varphi }.$$
(40)
The $`\mathrm{sign}`$ function in the above formula selects the root that yields the correct limit $`\alpha =0`$ when $`\varphi 0`$. The formula (40) is applicable in the interval $`\varphi [0,\pi /2]`$. In the limiting case $`\varphi =\pi /2`$ we have

$$\mathrm{tanh}\alpha _{\pi /2}=\frac{\mathrm{sign}(ss^1)}{\sqrt{s^21+s^2}},$$
(41)
which implies that $`|\mathrm{tanh}\alpha _{\pi /2}|<1`$ iff $`s1`$. Furthermore, it can be shown that $`\alpha `$ is a monotonic function of $`\varphi `$ in the interval $`[0,\pi /2]`$. We can thus conclude that with the interaction Hamiltonian of the amplifier type (31) we can implement any beam splitter transformation (20) with the mixing angles in the interval $`[0,\pi /2]`$, which includes the two important cases of a balanced beam splitter ($`\varphi =\pi /4`$) and the swap ($`\varphi =\pi /2`$). It follows from the expressions (39) and (40) that the simulation becomes more and more time consuming for Hamiltonians close to the two-mode squeezing Hamiltonian $`H_{\mathrm{TMS}}=x_Ap_B+p_Ax_B`$, i.e., when $`s1`$.

### IV.2 Beam splitter-like Hamiltonians

Having derived the gate synthesis parameters for the Hamiltonians (31), we proceed to the interaction Hamiltonians leading to oscillatory dynamics,

$$H_1=x_Ap_Bs^2p_Ax_B,H_2=s^2x_Ap_B+p_Ax_B.$$
(42)
The S matrices associated with these Hamiltonians read

$$𝑺_1^{}(t)=\left(\begin{array}{cc}\mathrm{cos}(st)& s\mathrm{sin}(st)\\ \frac{1}{s}\mathrm{sin}(st)& \mathrm{cos}(st)\end{array}\right),$$
(43)
$$𝑺_2^{}(t)=\left(\begin{array}{cc}\mathrm{cos}(st)& \frac{1}{s}\mathrm{sin}(st)\\ s\mathrm{sin}(st)& \mathrm{cos}(st)\end{array}\right).$$
(44)
We shall not repeat the details of the derivations of the parameters $`\alpha `$ and $`\beta `$ and we only summarize the results here. The beam splitter transformation can be accomplished by the following choice:

$$\mathrm{tan}\beta =\frac{2\mathrm{sin}\alpha \mathrm{cos}\alpha }{\mathrm{cos}^2\alpha +(s^2+1+s^2)\mathrm{sin}^2\alpha }$$
(45)
and

$$\mathrm{tan}\alpha =\frac{s^1+s\sqrt{(s^2+1+s^2)(\mathrm{cos}\varphi )^2+1}}{(s^2+1+s^2)\mathrm{tan}\varphi }.$$
(46)
This reveals that simulation of any beam splitter with $`\varphi [0,\pi /2]`$ is possible and the parameters satisfy $`|\alpha |\pi /2`$ and $`|\beta |\pi /2`$.

Consider now the two-mode squeezing operation (24). After some algebra, one obtains

$$\mathrm{tan}\beta =\frac{2\mathrm{sin}\alpha \mathrm{cos}\alpha }{\mathrm{cos}^2\alpha (s^21+s^2)\mathrm{sin}^2\alpha }$$
(47)
and

$$\mathrm{tan}\alpha =\frac{s^1s\mathrm{sign}(s^1s)\sqrt{\frac{s^21+s^2}{\mathrm{cosh}^2r}1}}{(s^21+s^2)\mathrm{tanh}r}.$$
(48)
The function $`\mathrm{tan}\alpha `$ must be real which implies that the term under the square root must be non-negative. This constraint, in turn, limits the amount of two-mode squeezing that can be produced via a three-step protocol (34). It holds that $`rr_{\mathrm{th}}`$ where

$$\mathrm{cosh}r_{\mathrm{th}}=\sqrt{s^21+s^2}.$$
(49)
The origin of this bound lies in the fact that the dynamics governed by the Hamiltonians (42) and captured by the matrices (43) and (44) is oscillatory and fully periodic with period $`2\pi /s`$. Squeezing above the threshold $`r_{\mathrm{th}}`$ can be achieved only if we concatenate several three-step protocols. It thus appears that the amplifier-like Hamiltonians (31) are in certain sense more versatile than the beam-splitter like Hamiltonians (42), because the former allow to implement any two-mode squeezing gate and also any beam splitter with $`\varphi [0,\pi /2]`$ via a three-step protocol (34).

## V Conclusions

In this paper we have addressed the problem of gate synthesis for continuous variable systems. We have assumed that two single-mode systems $`A`$ and $`B`$ interact via a quadratic Hamiltonian. We have studied how to implement a unitary symplectic gate $`G`$ with the use of the interaction Hamiltonian $`H`$ as a resource. The gate synthesis protocol consists of a sequence of evolutions governed by $`H`$ and followed by fast local phase shifts applied to the systems $`A`$ and $`B`$. We have focused on the gate simulation protocols that involve the minimal necessary number of steps, because these protocols require a low number of local control operations which is important from the experimental point of view. We have shown that a three-step protocol suffices for simulation of the two-mode squeezer as well as a beam splitter. For the specific case of the Hamiltonian (3) we have also established a four step implementation of a single-mode squeezer. Our results are applicable to any physical systems coupled via quadratic Hamiltonians. In particular, the gate synthesis protocols proposed in the present paper may find applications in the experiments where light interacts with atomic ensembles via a Kerr-like coupling Kuzmich00 ; Julsgaard01 ; Schori02 .

###### Acknowledgements.

I would like to thank S. Massar, E.S. Polzik, N.J. Cerf and F. Grosshans for helpful discussions. I acknowledge financial support from the Communauté Française de Belgique under grant ARC 00/05-251, from the IUAP programme of the Belgian governement under grant V-18, from the EU under projects RESQ (IST-2001-35759) and CHIC (IST-2001-32150) and from the grant LN00A015 of the Czech Ministry of Education.