An integral-field spectroscopic strong lens survey

We present the observational results of a survey for strong gravitational lens systems consisting of extended emission-line galaxies lensed by intervening early-type galaxies, conducted using integral field units (IFUs) of the Magellan IMACS and Gemini GMOS-N spectrographs. These data are highly valuable for corroborating the lensing interpretation of Hubble Space Telescope imaging data. We show that in many cases, ground-based IFU spectroscopy is in fact competitive with space-based imaging for the measurement of the mass model parameters of the lensing galaxy. We demonstrate a novel technique of three-dimensional gravitational lens modeling for a single lens system with a resolved lensed rotation curve. We also describe the details of our custom IFU data analysis software, which performs optimal multi-fiber extraction, relative and absolute wavelength calibration to a few hundredths of a pixel RMS and nearly Poisson-limited sky subtraction.


Introduction
In the pursuit of strong lensing science, spatial resolution is of paramount importance, and therefore high-resolution imaging with the Hubble Space Telescope (HST) currently sets the gold standard at optical wavelengths. This paper presents a different observational approach: a spatially resolved spectroscopic survey for strong galaxy-galaxy gravitational lenses using the optical integral-field units (IFUs) of the Gemini GMOS-N and Magellan IMACS spectrographs. A significant goal of this paper is to demonstrate the potential of this technique for the confirmation and study of strong lenses. We will show that, in comparison with HST imaging, ground-based IFU spectroscopy is always complementary, often competitive and occasionally superior. These favorable comparisons are all due to the combination of large collecting aperture with simultaneous spatial and spectral resolution, currently only available at ground-based observatories. Specifically, we show that IFU observations can be valuable for corroborating the lensing interpretation of high-resolution imaging data by using the third (wavelength) dimension to cleanly separate the narrow-band images of two separate galaxies along a single line of sight. We show that in many cases, ground-based IFU data alone can provide competitive constraints on the mass model of the lensing galaxy as compared to HST data. Finally, we show that IFU data can be superior to space-based imaging in cases where the lensed emission-line flux gives a clearer picture than the lensed continuum, or when the presence of resolved velocity structure in the lensed galaxy allows a three-dimensional (3D) approach to lens modeling.
Numerous previous studies have used the unique 3D capability of high spatial sampling IFUs for galaxy-scale strong lensing science. The majority of these works have used IFU spectroscopy to measure differential substructure lensing effects on the continuum and emission-line regions of lensed quasars (Lamer et al 2006, Mediavilla et al 1998, Metcalf et al 2004, Motta et al 2004, Sugai et al 2007, Wayth et al 2005, Wisotzki et al 2003. IFUs have also been used to resolve the kinematic structure of high-redshift galaxies lensed into magnified 3 arcs by intervening galaxy clusters (Swinbank et al 2003(Swinbank et al , 2006(Swinbank et al , 2007. The results we present here represent the first application of IFU spectroscopy to the confirmation and modeling of spectroscopically selected strong gravitational lens candidates consisting of a superposition of two galaxies along a single line of sight. A further goal of this paper is to present the details of our custom IFU data calibration and analysis software. This software (described in section 3) performs optimal extraction, relative and absolute wavelength calibration to a few hundredths of a pixel RMS and nearly Poissonlimited sky subtraction.

Sample and observations
The lens candidates observed for this study were selected from the spatially unresolved spectroscopic database of the Sloan Digital Sky Survey (SDSS; York et al 2000) using the technique described in Bolton et al (2004). They are identified through the presence of an emission-line redshift significantly higher than the absorption-line redshift within the same spectrum, as observed through the 3 -diameter SDSS spectroscopic fiber aperture. By design of the selection process, the candidates consist of a close angular superposition of a bright foreground early-type galaxy with a faint background star-forming galaxy. If the (initially unknown) impact parameter is small enough, the background galaxy will be multiply imaged by the gravity of the intervening early-type, furnishing a powerful probe of the distribution of mass in the latter. Many such candidates, including all the confirmed lenses presented in this paper, have been successfully observed with the Advanced Camera for Surveys (ACS) aboard the HST by the Sloan Lens ACS (SLACS) Survey (Bolton et al 2005. The observations for this work were obtained with IFUs built by the University of Durham Astronomical Instrumentation Group for the Inamori Magellan Areal Camera and Spectrograph (IMACS; Bigelow and Dressler 2003) and the Gemini-North Multi-Object Spectrograph (GMOS-N;Hook et al 2003). The data were collected during March-May 2004 (GMOS-N) and September 2004 (IMACS). IMACS data were also collected during November 2003 and March 2004, but observing and instrument conditions were too poor for scientific use. Table 1 gives information on the candidate systems for which data are presented in this paper.
The IMACS IFU (Schmoll et al 2004) observes two 5 × 7 fields of view in the focal plane, separated by roughly one arc-minute: one FOV (field of view) for the object and one from which to estimate the sky background. The fields are sampled by a close-packed hexagonal array of lenslets which subtend 0.2 from side to side, for a total of 2000 lenslets between the two fields. The lenslets feed the light to optical fibers, which reformat the fields via a defined fieldmapping into a 1D array of output lenslets (a 'pseudo-slit') for dispersion. (See figure 1 for a schematic diagram of this mapping.) All of the IFU optics are confined to the space of a narrow cartridge whose thickness occupies three adjacent mask slots in the slit-mask server. The mask server inserts and removes the IFU in the same manner as a simple mask, and thus the role of the IFU in the spectrograph can be understood as that of a special slit-mask, with two 5 × 7 'slits' on the input side and one long 'slit' on the output. The IMACS IFU can be used with either of the two IMACS cameras: f /4 ('long') or f /2 ('short'). All IMACS data presented in this paper were obtained using the f /2 camera, which uses grisms for dispersion. The GMOS-N IFU (Allington-Smith et al 2002;Murray et al 2003) operates in an identical manner to the IMACS IFU, with a few notable exceptions. Firstly, while the object field is 5 × 7 (1000 lenslets) as in Table 1. Information for double-redshift systems with IFU data presented in this paper. SDSS spectroscopic ID consists of plate number, modified Julian date and fiber number. The columns z FG and z BG respectively contain foreground and background redshifts measured for each system from SDSS spectroscopy. The magnitudes r are extinction-corrected (Schlegel et al 1998)  IMACS, the background field is half this size: 5 × 3.5 (500 lenslets). Secondly, the 1500 total fibers are reformatted not into one single pseudo-slit, but rather into two parallel pseudo-slits separated by approximately 3200 pixels on the detector. In the 'two-slit' mode employed in the current work, broad-band filters must be used to limit the wavelength domain of the individual pseudo-slits so as to prevent overlapping of spectra. Table 2 lists the various unique spectrograph configurations used in the current work, along with their general characteristics. Information on which configurations were used for which targets is given in table 3 (further below).

Data calibration and analysis
The format of the data delivered by the GMOS-N and IMACS IFUs to their respective CCD-mosaic detectors is sufficiently complicated to justify analysis with specially developed software. Furthermore, the scientific application of strong lens modeling requires highly accurate wavelength calibration, flat-fielding and sky subtraction. For these reasons, we have developed our own set of IFU data-analysis software tools written in the IDL language.
Although there is an officially maintained and distributed package for the reduction of GMOS IFU data under IRAF, there is no such package for IMACS. Having developed our own IMACS routines, we preferred adapting them to GMOS usage over using the official GMOS package. Schematic representation of the focal-plane to detector-plane mapping implemented by the combination of the IFU with a dispersive element (grism or grating). The 5 × 7 field of view in the focal plane is fully sampled by 1000 hexagonal lenslets with side-to-side diameter of 0.2 . The lenslets feed optical fibers, which reformat the lenslet rows in blocks of 50 as shown. The blocks are arranged end to end with some inter-block spacing along a 'pseudo-slit' and dispersed in wavelength to produce 1000 spectra on the CCD. A separate and similarly sampled field of view is used for sky estimation, and the 50-fiber blocks from this field alternate in sequence with the object-field blocks on the detector mosaic. Table 2. Unique Gemini GMOS-N and Magellan-IMACS IFU observing modes used to collect the data presented in this paper. Three additional targets were observed using GMOS-N in the g band, but the data were not successfully reduced due to problems with tracing, flat-fielding and wavelength calibration. All GMOS-N IFU observations were binned 2× in the dispersion direction. We refer to our software as 'kungifu' (kung eye eff you), and in this section we describe its function. Figure 2 shows a flow-chart outlining the operation of the software. The kungifu package can be obtained for use by other investigators by contacting the authors, though the authors are unable to offer user support beyond the distributed documentation. The kungifu software is distributed along with template scripts and associated demonstration data frames, which allow interested users to reproduce all of the operations described in this section and associated inverse-variance  to examine the results graphically. The template scripts can then be adapted to suit the user's particular analysis needs.

Bias subtraction and data formatting
The IMACS detector consists of a mosaic of eight 2048 × 4096-pixel CCDs. The bias level of raw frames varies between the CCDs, within each CCD and from one exposure to the next. For 7 any given exposure, the bias is well characterized as the sum of row-overscan and columnoverscan terms. Our IMACS bias-subtraction routine first fits a smooth b-spline model (de Boor 1977) to the overscan region at the end of each row for a specified breakpoint spacing. This model is evaluated for all pixels in all rows of the image and subtracted. The IMACS bias level also has significant column-dependence, and a similar b-spline fit to the overscan region at the end of each column is subtracted next. The bias-subtracted and overscan-trimmed individual CCD images are then stored to a single multi-extension FITS (MEF) file. The detailed GMOS-N bias pattern is more stable than the IMACS bias pattern, but cannot be estimated from the overscan alone. Thus we perform GMOS-N bias subtraction using mean bias images and also subtract an overall average bias value for each of the three GMOS-N CCDs from the overscan region to account for slight variations. In two-slit IFU mode, the central CCD records spectra from both pseudo-slits. We break the central CCD into two logically separate images to separate these two regions and store them as separate MEF extensions along with the first and third CCDs, after bias subtraction and overscan trimming.

Flat-field modeling and tracing
Relative flux calibrations of the IMACS and GMOS-N IFUs are accomplished using dispersed IFU exposures of uniform spatial illumination by (approximately) flat-spectrum incandescent lamps. Note that highly uniform spatial illumination with facility calibration equipment over such small IFU fields of view is much more easily achieved than similarly uniform illumination over wide direct-imaging fields of view. These raw spectroscopic flat images measure the product of two non-uniform responses: the relative sensitivity of individual CCD pixels, and the fiber response image (including both individual fiber profiles and relative throughput differences between fibers). In an ideal spectrograph the illumination pattern of the fibers would be fixed relative to the detector, and these two effects would never need to be distinguished. In actuality both GMOS-N and IMACS exhibit limited flexure between successive exposures that causes the IFU fiber spectra to shift their position in CCD coordinates. Thus, we perform a factorization of the pixel-response and fiber-response calibrations by assuming an approximate scale separation between them. Since the IFU fiber spectra for all data obtained for this work run approximately along CCD rows, horizontal cross-sections through the IFU spectroscopic flat frames follow the smooth variation of the flat-lamp spectrum, modulated by gradual transitions from one fiber to the next. We generate smooth models of the spectroscopic flat-field images by fitting b-spline models to these cross-sections, with a breakpoint spacing chosen ideally to be smaller than the typical scale of flat-lamp features but greater than the scale of pixel-to-pixel defects. The resulting model-flat images are an approximation to the fiber response, and the ratio of rawto model-flat images give the approximate pixel response ('pixel flat'). Figure 3 demonstrates this flat-field factorization graphically. We derive a master pixel flat from the median-image of many individual pixel flats in a given spectrograph configuration, for application to all science frames in that configuration. We choose not to use imaging flats to calibrate the pixel response, since this calibration can in general be wavelength-dependent. Using pixel flats derived from spectroscopic flat-field frames ensures that the correction is derived with illumination at the particular wavelength appropriate to each pixel. We describe our use of the model-flat images for the calibration of fiber response in section 3.4 below.
We also use the model-flat images to determine the location of individual fiber traces on the CCD. For both the IMACS and GMOS-N IFUs, fibers along the output pseudo-slit are  Cross-section through an IMACS-f /2 model flat for two 50-fiber blocks. Note the substantial variation in fiber-to-fiber throughput and the significant overlap of the wings of adjacent fibers. Also note the nonzero level in between blocks, indicating a scattered-light contribution to the count levels.
grouped into blocks of 50, with median inter-fiber spacings on the detector of 3.5 and 5.7 pixels respectively for IMACS (in the f /2 or 'short' camera mode) and GMOS-N. Figure 4 shows a cross-section through an IMACS-f /2 model flat for two 50-fiber blocks. Fibers are approximately equally-spaced within the blocks, and we use this fact to our advantage to locate and trace all 50 fibers in a block at once. For the purpose of locating block positions relative to one another, our automatic tracing routine also makes use of a table of inter-block separations in units of the approximate local inter-fiber spacing (which will be independent of pixel scale), which is determined once for each IFU by a careful analysis of a model-flat cross-section.

Scattered-light subtraction
All IMACS and GMOS-N spectroscopic exposures exhibit a non-negligible scattered-light background not directly associated with the flux through the fibers, as evidenced by nonzero count levels in the inter-block regions where the flux from the fibers drops essentially to zero. For all IFU frames used in this work-calibrations and science exposures-we subtract a model of this background after pixel-flat correction. We estimate this scattered-light image from the observed levels between the fiber blocks. We use the flat-field-derived trace solution to define bands in the inter-block regions running parallel to the spectra, with a reasonable buffer to avoid the wings of the fiber cross-sections. Each band is fit with a b-spline as a function of column, and this fit is subsequently evaluated for each column. We then interpolate this fit across the fiber blocks by fitting a b-spline in each column, taking the band centers in that column as dependent 9 variables and the previous-fit evaluations in that column as independent variables. The b-spline breakpoint spacing in each fit may be adjusted according to the signal-to-noise in the scattered light levels and the degree of structure that one wishes to model.

Extraction
IMACS and GMOS-N IFU observations distribute the photons from a few dozen square arcseconds of the sky over a few tens of millions of CCD pixels. Thus for all but the brightest objects, low signal-to-noise presents a challenge. This fact combined with the well-behaved profile of the IFU fibers on the CCD suggests optimal spectrum extraction as a natural approach (e.g. Hewett et al 1985, Horne 1986). In an optimal extraction, the specific flux of a spectrum at a given wavelength is determined not from a simple sum over pixels within a defined aperture, but rather from the amplitude of a maximum-likelihood fit of a model cross-section to the observed spectral cross-section at that wavelength. This gives the maximum signal-to-noise in the extracted spectrum, and is unbiased to the extent that the model cross-section matches the actual cross-section. The most apparent obstacle to optimal extraction of IMACS and GMOS-N IFU data is the significant overlap between neighboring spectra. Fortunately, the situation is less dire for fiber-fed IFUs such as those of IMACS and GMOS-N than for multi-object multi-fiber spectrographs, since adjacent fibers on the detector are also adjacent on the sky (an explicit design feature). Since the 0.2 -diameter IFU lenslets will critically sample all but the very best ground-based seeing, the blending of neighboring fiber spectra on the detector leads to no significant loss of information (Allington-Smith et al 2002). Furthermore, modeling of the fiber profile is unnecessary, since the model flats described in section 3.2 provide us with a high signal-to-noise determination of the spectrum cross-sectional shape for all fibers and at all wavelengths. Before using the model flats for the extraction of spectra from the pixel-flat corrected and scattered-light subtracted object frames, we first normalize the model flats (also scattered-light subtracted) by dividing out a crude approximation to the flat-field lamp spectrum as a function of wavelength. This step is not crucial (particularly if one eventually performs an absolute flux calibration), but it prevents the flat-lamp spectrum from being imprinted on the data before flux calibration. We note that it is not important to use an exceedingly accurate model of the lamp spectrum, but only to divide all pixels with the same wavelength by the same value. Next, we shift the model flats perpendicular to the dispersion direction with a flux-conserving dampedsinc kernel so as to maximize the cross-correlation between the model-flat image and the objectframe image to be extracted. This shifting of the model flat accounts for the slight (typically of order 1 pixel or less) flexure that can occur between the object frames and the flat frames taken immediately following.
Our approach to extraction is described mathematically as follows. We define boundaries between fiber spectra by lines exactly half-way between the fiber traces, and in each CCD column i we associate with fiber j all pixels (i.e. rows) k falling between the boundaries on either side. Pixels split by the boundary are associated fractionally (see figure 5). Let w (i, j) k express this weighting: w (i, j) k = 1 for rows k wholly within the boundaries for fiber j in column i, 0 for rows wholly outside, and between 0 and 1 for rows fractionally included. Let d ik be the data frame to be extracted, σ 2 ik be the statistical variance of d ik , and f ik be the aligned, normalized model flat-field image of the fiber response. The optimally extracted specific flux in fiber j at column i (corresponding to a particular wavelength by the dispersion solution for that fiber), which we denote I i j , will be given by the value that minimizes which is A simple adjustment of this expression suggests a more succinct conceptual and operational approach: This instructs us to obtain the optimally extracted specific flux by dividing the data image by the model-flat image, then computing a weighted average over the appropriate fiber/column window, with the statistical weight given by the product of the squared model-flat image and the inverse-variance image. We implement the extraction algorithm in this manner. The flat-fielding procedure thus yields calibrated specific intensity measurements (flux per unit wavelength per unit area) in a relative sense at any given wavelength across the entire spatial field, though the absolute extracted flux values are not meaningful unless and until an absolute flux calibration is applied.
The optimal-extraction technique also provides a natural means for rejecting cosmic-ray (CR) hits in single exposures, because CRs generally will not have the same shape as the fiber cross-section. We flag pixels with highly statistically significant positive deviations between the data frame and the optimal-extraction model frame as CR pixels and grow the resulting CR mask to laterally-adjacent pixels, then repeat the extraction with CR pixels given zero weight. Figure 6 illustrates various elements of the extraction process for a small subregion of one IMACS CCD.
We store the extracted IFU object spectra for each CCD in the IMACS and GMOS-N mosaics as a separate image extension in a single MEF file for a single exposure. These spectrum images have a horizontal dimension equal to that of the CCD and a vertical dimension equal to the number of fibers with spectra falling on the CCD, and are not rectified in wavelength but rather have a unique wavelength sampling for each spectrum. We note here that storing the data in this 2D spectrum-image form instead of in the 'data-cube' form of two spatial dimensions by one spectral dimension affords a distinct advantage even for wavelength-rectified data (see section 3.7) in that it preserves the integrity of the detector frame.

Relative and absolute wavelength calibration
We establish wavelength calibration using exposures of He-Ne-Ar (for IMACS) and Cu-Ar (for GMOS-N) arc lamps. We first process arc images by subtracting the bias level, dividing by a pixel flat, and subtracting a scattered-light image model. We then use the most closelyassociated (in time) model spectroscopic flat, globally sinc-shifted perpendicular to the dispersion direction if needed for alignment, to perform an optimal extraction as described in section 3.4. The model flat image is not normalized in wavelength since the wavelength calibration is as yet unknown. Within the resulting set of extracted arc-spectra, we identify as many individual bright lines as possible. The centroids of these lines are found in each spectrum through an iterative linearized Gaussian centroiding algorithm that also measures the width of each line in the dispersion direction. The median line-spread width for each fiber is saved for later use.
Due to global curvature of the pseudo-slit image on the detector mosaic, each fiber spectrum has its own unique dispersion solution, as would be the case for the individual rows in a single long-slit spectrum. Furthermore, discrete offsets in the dispersion direction of 50-fiber blocks relative to one another, as well as offsets of individual fiber spot positions from the mean within their blocks, make a global 2D dispersion solution impractical. Rather, we use the set of 12 arc-line positions in each fiber to derive the individual dispersion solution in that fiber relative to a baseline defined by the average arc-line position across all fibers, which we refer to as the 'pixel-wavelength' baseline. We derive this solution in the following iterative fashion. First, for each arc line, we compute the simple average centroid position (CCD column) across all fibers, which we take to be the pixel-wavelength of that line. We then fit the position of all lines in each fiber spectrum with a low-order (linear or quadratic) polynomial function of this single pixelwavelength baseline. We then use the average residuals of these fits across all fibers to correct the adopted pixel-wavelength values for each line and refit. This process converges rapidly after a few iterations, and the resulting polynomial solutions are accurate to within a few hundredths of a pixel RMS difference between arc-line centroids as measured and as predicted by the fit (equivalent to a few hundredths of an Angstrom RMS at the dispersions used in this work). This process-which amounts to the determination, but not the application, of a wavelengthrectifying solution-has several advantages. Firstly, one need not know the identifications of the arc lines used. Secondly, one need not exclude blended lines, since any bias introduced by blended lines into the independent variable (pixel) will be exactly canceled in the dependent variable (pixel-wavelength). Finally, the strategy reduces absolute wavelength calibration to a 1D problem.
Once the solution for pixel as a function of pixel-wavelength has been determined for each individual fiber, we determine the single solution for physical wavelength as a function of pixel wavelength. We first use the (inverse) pixel-wavelength solution to determine the central pixel-wavelength of each pixel in the extracted arc frame. We then perform a single b-spline fit to the extracted arc flux as a function of pixel wavelength across all fibers. From a signalto-noise point of view we are effectively stacking all the individual arc spectra, but rather than re-binning/rectifying in (pixel-)wavelength, we fit to the data in their native pixel sampling. We identify peaks in the fitted b-spline model, and fit a polynomial solution for the absolute wavelength calibration. A further benefit of our decomposition of the wavelength calibration problem into relative and absolute steps is that even arc lines detected at signal-to-noise ratio (SNR) 1 in individual fiber spectra can be robustly centroided in the global b-spline model and used to constrain the absolute solution. Note that the separation into relative and absolute wavelength calibration is entirely analogous to the usual separation of flux calibration into the relative step of flat-fielding and a subsequent step of absolute calibration using a flux standard.
When applying the wavelength calibration to the science data, we also correct for slight flexure between the object frames and the dispersion solution from the calibrating arc frame by fitting for a low-order 2D polynomial transformation as defined by the positions of known night-sky emission lines.

Sky subtraction
Both the GMOS-N and IMACS IFUs incorporate dedicated fiber bundles for the observation of blank sky, from which to estimate the foreground level to be subtracted from the science data. As discussed by Kelson (2003), there is a distinct advantage in the estimation and subtraction of the night-sky spectrum before performing any rectification in wavelength. The multiple fibers of the IFU fields of view each have a slightly different wavelength sampling on the detector, and hence the discretely sampled line-spread function (LSF) observed for night-sky emission lines depends upon the sub-pixel location of the line's central wavelength. The native binning of the CCD, when considered for all blank-sky fibers together, provides a finely-sampled observation Figure 7. GMOS-N i-band sky subtraction for a single 900 s exposure: before (above) and after (below). A total of 100 fiber spectra are displayed in each case, arranged vertically. Wavelength direction is horizontal, and covers approximately the same range as is shown in figure 8. of the night sky spectrum. We thus fit a b-spline model to this data as a function of the central pixel-wavelength of each native spectral pixel, which we then evaluate for all fibers (object and sky) and subtract as our sky model. The extracted 1D LSF of the IMACS and GMOS-N IFUs exhibits some variation due to global distortions and fiber-optic heterogeneity, which we characterize via the LSF width measured for each fiber from arc frames as described in section 3.5 above. This LSF width is treated as a second independent variable in the b-spline model, fitted with linear or quadratic dependence. Figure 7 shows the results of i-band sky subtraction for a single 900 s GMOS-N frame, and figure 8 characterizes the typical statistical quality of sky subtraction across the observational sample.

Rectification and combination
To facilitate the combination of multiple exposures (and to make analysis more straightforward), we rebin our sky-subtracted IFU spectra on to a uniform wavelength baseline. We use a constantwavelength binning across the spectrum, with bin size slightly larger than the largest native pixel-width, so that only nearest-neighbor correlations will be present in the rebinned frames. The wavelength-bin boundaries are specified in heliocentric vacuum wavelengths, corrected with a heliocentric velocity shift appropriate to each observation, converted to air wavelengths, and mapped into the extracted IFU frames using the arc-frame dispersion solution, as adjusted to match the observed night sky lines. Multiple exposures are then combined after re-binning, with further CR rejection. Finally, the rebinned data from individual CCDs are combined onto a single mosaic image that maintains the native orientation of the spectra on the detector mosaic. As a data product, we prefer this mosaic to a 3D 'data cube' (2 spatial plus 1 spectral dimension) because it allows the reduced data to be displayed all at once in the frame of the detector. Data cubes may always be constructed using the IFU field-mapping. We do not rebin our data spatially, since all observations were made with single undithered telescope pointings.

Absolute flux calibration
Flux calibration is currently not implemented as an integral aspect of the kungifu tasks, and is not necessary for gravitational-lens modeling using emission-line images at a single wavelength. Wavelength-dependent IFU sky subtraction performance for IMACS (top) and GMOS (bottom) in an i-band wavelength range containing strong OH rotational transition lines. Shown in gray is the median across fibers of the ratio of sky counts to noise counts. The RMS across fibers of the ratio of sky-subtracted residual counts to noise counts is shown in black. When the noise is estimated correctly and the sky is subtracted perfectly, this second ratio will be 1. Pixel size is approximately 1.3 Å for IMACS and 0.9 Å for GMOS-N (including a binning factor of 2). Curves are computed for each exposure separately using all sky fibers; curves shown are a median across all individual exposures used in this paper. No adjustment is made for differing exposure times, airmasses or moon phases between the individual exposures.
In order to judge the depth of our observations and the scale of emission-line luminosities observed, we implement a somewhat crude flux calibration for our gravitational-lens candidate sample through a bootstrap connection to the flux-calibrated SDSS spectra, which were obtained through a single 3 -diameter fiber aperture that can be synthesized from within the IFU field of view.

Narrow-band image construction
The 3D data provided by IFUs afford the opportunity to construct narrow-band images of any bandwidth and at any wavelength within the data cube. Furthermore, emission-line regions of the spectrum may be decomposed into continuum and emission-line components through suitable profile modeling. To a very good approximation-largely as a consequence of our initial selection method-continuum emission seen in the IFU spectra of our lens candidates can be attributed to the foreground galaxy. Similarly, high-redshift emission lines are of course entirely due to the background galaxies. Thus decomposition into continuum and emission line components is approximately equivalent to decomposition into foreground and background galaxy images.
Our strategy for continuum/emission-line decomposition is as follows. First, we select a small (± ∼ 10 Å) wavelength range about the wavelength of a background emission line detected in the SDSS spectroscopy. We parameterize the emission-line component as a Gaussian function of wavelength, and the continuum as either a constant or linear function of wavelength, which are sufficient for our purposes over these small wavelength ranges. The central wavelength and width of the Gaussian term, as well as the slope of the continuum (parameterized as a fractional change over the wavelength range considered), are treated as global model parameters applying to all fiber spectra. For a given trial set of these parameters, a basis is constructed and the amplitudes of the continuum and emission-line components are fitted with a linear least-squares fit to the IFU data over the selected wavelength range. This linear fit is wrapped within a nonlinear optimization of the global emission-line and continuum shape parameters that minimizes the overall χ 2 using the MPFIT IDL implementation of the Levenberg-Marquardt algorithm. For the case of [O ii] 3727 emission, a double Gaussian line profile is used, with the line spacing fixed to its known rest-frame value of 2.78 Å and the relative strength of the two line components fitted as an additional global parameter. In several systems, slight velocity shifts are detected in the background emission lines, and these are modeled with an additional Gaussian-derivative amplitude that approximates the small velocity shifts. One system shows an appreciably resolved rotation curve, which is treated separately as described below.
The linearly fitted coefficients of the final parameterized basis functions are our best decomposition of the spectrum into background (emission-line) and foreground (continuum) components. Using these components, we can form emission-line and continuum images in the focal plane of the telescope using the known IFU field mapping. Note that the modeling process yields optimal signal-to-noise ratio for these images, better than simple flux-averaging within defined continuum and emission-line wavelength aperture. Figures 9, 10 and 12 show reconstructed narrow-band images of those systems for which the original SDSS emission-line detection was confirmed by IFU observation, with related information given in table 3. IFU observations of several additional candidate systems described in table 4 failed to confirm the original SDSS line detection, and the associated data are not presented in this paper.

Strong lens morphology and modeling
The emission-line images of the systems shown in figure 9, along with a number of systems in figure 10 (specifically, J0956, J1259, J1416, J1521, J1630 and J1702), show multiple features that suggest strong gravitational lensing. Bolton et al (2006) show how many of the same narrow-band IFU emission-line images confirm the lensing interpretation of HST imaging data by showing the spatial coincidence of the detected high-redshift emission-line features with the putative lensed images seen with HST. Confirmation of a lens candidate as a gravitational lens solely through IFU observation, however, relies on the successful fitting of a lensing mass model to the putative strongly lensed images. Here, we describe the IFU lens modeling strategy that we apply to all possible lens images, with successful results for the 7 systems shown in figure 9.
Where the signal-to-noise ratio and spatial resolution allow, we analyze all candidate lens systems uniformly by fitting a singular isothermal ellipsoid mass model (SIE: Kassiola and Table 3. Observation and data display details for double-redshift systems with IFU data presented in this paper. Observation modes are described in table 2. 'IFU line' column specifies the particular background emission line used to generate the narrow-band images and models shown in figures 9, 10 and 12 (with observed wavelengths given by the background redshifts in table 1). 'Figure levels' columns give the limits of the gray-scaling used in figures 9, 10 and 12. Units are 10 −17 erg cm −2 s −1 Å −1 lenslet −1 for continuum images and 10 −17 erg cm −2 s −1 lenslet −1 for emission-line images. 'Lens model?' column specifies whether or not a singular isothermal ellipsoid lens model was successfully fit to the narrow-band emission-line data. The higher success rate for lens modeling with IMACS data over GMOS data is a consequence of the selection of more promising candidates for IMACS observation: i.e. those with higher line fluxes and larger anticipated Einstein radii. The IMACS targeting choices were greatly informed by experience with the earlier GMOS observations.  Figure 9. Lens systems that admit successful IFU gravitational-lens modeling. There are three panels for each system: narrow-band continuum level (left), narrow-band emission-line level (center), and model emission-line level as reproduced by the best-fit gravitational lens model (right). All panels are 5 × 7 . Gray-scale is linear, and varies from system to system and between continuum and emission-line panels for visual ease. Individual gray-scale limits are given in table 3. Black corresponds to values at the limit and above, and white corresponds to values at −0.25 times the limit and below. In reading order from upper left, the individual systems are: J0037, J0044, J0737, J1402, J2238, J2302 and J2321.
where r q = q x 2 + y 2 /q, and x and y are angular coordinates aligned along the major and minor axes of the iso-density contours. The dependence of the model on x and y only through r q means that all iso-density contours are aligned, similar (not confocal), concentric ellipses. The parameter b expresses the strength of the lens and q gives the minor-to-major axis ratio of the iso-density contours (hence 0 < q 1 by convention). When q = 1, the SIE reduces to the singular isothermal sphere (SIS)  and b is equal to the angular radius of ring images of on-axis sources. This 'Einstein radius' is in turn related to the velocity dispersion σ of the lensing distribution through (D LS and D S are angular-diameter distances from lens to source and observer to source.) For purposes of comparison between models with differing q values, we adopt the same intermediate-axis normalization as Kormann et al (1994), whereby the mass interior to a given iso-density contour at fixed b is constant with changing q. A great advantage of the SIE is that its projected gravitational potential can be expressed analytically. For our purposes we need only the derivatives of this potential, which we use in the simple form given by Keeton and Kochanek (1998). We perform all lens modeling with our own IDL routines. As is always the case when fitting gravitational lens models, we must take care not to fit for more parameters than are constrained by the data. Fortunately the SIE is generally free from fundamental degeneracies among its parameters with regard to the constraints furnished by real lenses. In this paper, we are concerned with spatially extended lensed star-forming galaxies which will generally be of sufficient physical size to average over any microlensing effects, and thus we use image surface brightness rather than image positions to constrain our lens models. We begin the model-fitting procedure for each individual system with a best-guess choice for the parameters of the lens model and source galaxy that gives a reasonable qualitative approximation to the observed emission-line image morphology. For systems that show spatially resolved emissionline brightness distributions (most cases), we fit for lens and source parameters by generating an unlensed source-galaxy image (either Gaussian or Sérsic, with ellipticity if necessary) and viewing it through the potential of the parameterized lens model. This image is then smeared by a Gaussian PSF model (whose FWHM is fixed to 0.7 for all systems in the current analysis), integrated over the hexagonal IFU lenslets, and used to calculate χ 2 directly relative to the narrow-band IFU data. The model parameters are optimized nonlinearly to minimize the χ 2 statistic. We test the sensitivity to the fixed PSF FWHM value for the systems J2238 and J2302 by increasing the fixed value in steps from 0.7 to 1.6 re-optimizing the model parameters at each step. Over this full range, the best-fit Einstein radius increases by 2.8% for J2238 and 2.2% for J2302. We also present the correlations between the Einstein radius b and other select model parameters for these two systems in table 6. For some lenses, we constrain the center of the lensing potential to be coincident with the center of the lensing galaxy, which we determine by fitting a Sérsic model to the IFU continuum image. For one lens consisting of unresolved images (J0037), we fit the images with hexagonally-sampled Gaussians (constrained to have the same width as one another). The image positions and fluxes from these fits are then used to constrain the lens model, with χ 2 computed from the image position and flux residuals. Table 5 gives the measured angular Einstein radii that result from our lens model fitting procedure. Also shown are the Einstein radii measured from a similar analysis of HST-ACS imaging of these same systems , see also Koopmans et al 2006).

3D lens modeling
With one exception, the internal velocity gradients in the lensed emission line galaxies observed by this program are unresolved spectroscopically. The exception is J1029, for which the rotation curve of the lensed galaxy is clearly visible. Figure 11 shows a section of the extracted IFU spectrum array for this system. The appearance of resolved velocity structure in the lensed source affords a unique opportunity for 3D strong lens modeling: two spatial dimensions plus one velocity dimension. Since multiple points in the lensed image plane coming from the same point in the source plane not only have the same surface brightness but also the Table 5. Lens strength b (SIE Einstein radius) as measured by IFU lens modeling. Quoted errors on the Einstein radii are square-root-diagonal entries from the covariance matrix of the nonlinear fit. Also shown are HST-ACS b-values from Bolton et al (2007), where available. Formal statistical errors on b HST are at the level of 1 to 2 milli-arcsec; independent analyses of HST lens data using the pixelized source-plane techniques of Koopmans et al (2006) give b values that agree within 0.03 RMS of the values obtained using the methods of Bolton et al (2007), with no evidence for relative systematic bias. IFU Einstein radii are converted to velocity dispersions σ SIE using equation (6), which can be compared with the stellar velocity dispersion σ ,SDSS measured from SDSS spectroscopy.  J2238 J2302 Mass minor-to-major axis ratio 0.06 0.27 0.02 0.27 Mass position angle (degrees) 1.9 −0.21 1.1 0.03 Mass x-center (arcsec) 0.03 0.13 0.01 −0.29 Mass y-center (arcsec) 0.05 −0.34 0.02 0.02 Lensed source x-center (arcsec) 0.03 0.14 0.01 −0.23 Lensed source y-center (arcsec) 0.02 −0.40 0.01 −0.31 same velocity, this technique can in principle constrain the lensing model in more detail than techniques using only velocity-unresolved data. In practice, the spatial resolution of our J1029 data is not sufficient to enable a high-quality analysis. Nevertheless, the wavelength-dependent (i.e. velocity-dependent) lensed image structure, seen in iso-wavelength slices through the continuum-subtracted data-cube in the first and third rows of figure 12, reveals intriguing behavior that invites a plausible explanation in terms of 3D lens modeling. We describe our modeling technique here, which predicts the data-cube slices shown in the second and fourth rows of figure 12. Though the data and model disagree in detail, the qualitative features of the data are reproduced by the model. We adopt the same SIE model as above to describe the mass distribution of the lensing galaxy. We model the lensed background galaxy as an infinitesimally thin, inclined exponential disk with a constant circular orbital velocity as a function of galactocentric radius and a single overall systemic average velocity along the line of sight. For any given set of parameters, we compute the model image and velocity field of the background galaxy as seen through the lens potential. We then convolve the image and the intensity-weighted velocity and squared-velocity images with the spatial PSF and the spectral resolution, and assign to each IFU fiber in the model data cube a Gaussian line profile with the appropriate line intensity, central wavelength, and broadening. The χ 2 statistic is then computed with respect to the continuum-subtracted datacube and the parameters are optimized nonlinearly. The model shown in figure 12 represents the best result obtainable. The mismatches in detail between data and model are likely attributable to the inadequacy of the parameterized model employed. Unfortunately the limited spatial resolution of the data prevent us from exploring more detailed models for the source galaxy, which might in turn allow us to constrain more detailed models for the mass distribution of the lens. (In the simpler 2D modeling context above, we experienced similar difficulty with J1630.)

Summary and conclusions
We have shown that the high spatial resolution IFUs of the Gemini-N and Magellan telescopes can be used successfully to confirm the spatially extended nature of higher-redshift line emission detected along the line of sight behind lower redshift elliptical galaxies. The third (wavelength) dimension permits exceptionally accurate decomposition of the reconstructed narrow-band images into continuum and emission line components. The emission-line images can then be analyzed in a lensing framework virtually free from contamination by the continuum of the foreground galaxy. In a number of cases, the IFU narrow-band emission-line images have sufficient resolution and signal-to-noise ratio to constrain gravitational-lens models that agree with models from high-resolution HST imaging to within 0.05 RMS (i.e. 1/4 of a lenslet). This result indicates that much-though certainly not all-of the science that can be done with Figure 12. Narrow-band IFU imaging of SDSSJ1029+6115 centered on the redshifted Hα line of the background galaxy. Leftmost panels in top two rows show continuum data and Gaussian-convolved Sérsic-model fitted to those data, respectively. Remaining panels show 1 Å slices through the continuumsubtracted data cube from 8210 to 8220 Å, with data above and lens model (described in section 6) below. The model provides a reasonable qualitative description of the data, though shortcomings are apparent.
space-based imaging of these lenses could perhaps be done with ground-based telescopes, provided sufficient aperture, time and image quality. Space-based imaging enjoys a great advantage in the interpretation and analysis of strong galaxy-galaxy lenses with complex and irregular spatial morphology in the lensed images. This is the case for J1630, shown in figure 10, which can be seen at HST resolution to consist of five distinct and irregularly oriented knots of emission in the source plane Koopmans et al 2006). However, even in cases where space-based imaging is clearly superior for lens modeling, IFU spectroscopy is of great value as an independent confirmation of the lensing hypothesis.
We have also demonstrated the technique of 3D strong gravitational lens modeling, which is currently only feasible using large ground-based telescopes such as Gemini and Magellan. Although the demonstration data presented here are of insufficient resolution for high-quality lens modeling, we consider this application to hold great promise, particularly if it can be deployed in combination with adaptive optics.