ALMA OBSERVATIONS OF SPT-DISCOVERED, STRONGLY LENSED, DUSTY, STAR-FORMING GALAXIES

The primary observations for this work were obtained from ALMA under a Cycle 0 program (2011.0.00958.S; PI: D. Marrone) in which 47 sources identified in the SPT survey are each observed in both the compact and extended array configurations. The first data release for this program includes 20 sources observed in the compact array configuration; here, we focus on the four sources for which these low-resolution data were sufficient to resolve the targets into multiple components. Of the remaining 16 sources, at least 8 are not point-like at the resolution of these observations, but we defer lens modeling for these sources until the remaining data are in hand.

The sources were targeted for brief snapshot observations with the dual-polarization Band 7 (275–373 GHz) receivers on two dates, 2011 November 16 and 28. The first local oscillator was set to 343.8 GHz, with all four spectral windows configured in time domain mode with 128 channels of 15.625 MHz width centered at 5.125 and 7 GHz IF in each sideband. There were 16 and 14 antennas available on these days, respectively, arranged in a compact configuration. The total elapsed time (for all sources, including those not published here) in the observations was approximately 4.1 hr. The total integration time per source was 61 and 91 s in the first and second tracks, respectively, with 6.1-s sampling of the visibility data. The array alternately observed the science targets and gain calibrators (30 s), observing the calibrator every three to four minutes. Additional sources with known positions, precisely established against the International Celestial Reference Frame using very long baseline interferometry (Ma et al. 1998), were added to the tracks to verify astrometry and calibration and observed with the same cycle as the science targets.

The flux scale was set with observations of Callisto on the first day. On the second day, the flux scale was derived by setting the flux of quasar J0403-360 to 1.84 Jy, as reported by the ALMA staff from adjacent calibration observations. The absolute flux density scale is correct to within 15%. The antenna gains are equalized through gain calibration (amplitude and phase) on the main calibrators in each track. Short-timescale phase correction is achieved using the ALMA water vapor radiometry (WVR) system. Very little variation is observed in these gain amplitudes through the tracks, and there is no evidence for atmospheric decorrelation on the longest baselines in the phase scatter of the target or calibrator visibilities after WVR phase correction. The data were processed with the Common Astronomy Software Applications package (McMullin et al. 2007; Petry et al. 2012) using standard steps for a continuum observation.

The four well-resolved sources are listed in Table 1, and we refer to them throughout the paper with shortened versions of their coordinate-based names (e.g., SPT0346-52 for SPT-S J034640−5205.1). Observations on November 16 included baselines of 15–150 kλ, resulting in synthesized beams of 15 × 13 (FWHM) for the RA = 5^h sources. The uv coverage on November 28 was less uniform and spanned 15–240 kλ, providing a synthesized beam of 21 × 09 for the other two sources. Deconvolved source images and beam shapes are shown in Figure 1.

Zoom In Zoom Out Reset image size

Table 1. Source and Lens Parameters

							Source Intrinsic Properties
ID	zS	zL	rE	ML	L	μ	R1/2	LFIR	ΣFIR	S1.4mm	S860μm
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
SPT-S J034640−5205.1	5.656	...	1.124 ± 0.004	3.73 ± 0.04a	0.55 ± 0.01	5.4 ± 0.2	0.59 ± 0.03	37.3 ± 2.8	24	8.1	23.0
SPT-S J041840−4752.0	4.224	0.265	1.390 ± 0.012	2.39 ± 0.04	0.20 ± 0.03	21.0 ± 3.5	1.07 ± 0.17	3.8 ± 0.7	0.74	1.6	5.0
SPT-S J052903−5436.6	3.369	0.140	1.536 ± 0.017	1.64 ± 0.04	0.10 ± 0.03	9.4 ± 1.0	2.39 ± 0.24	3.8 ± 0.5	0.15	3.8	15.6
SPT-S J053816−5030.8	2.782	0.404	1.987 ± 0.009	7.15 ± 0.05	0.13 ± 0.02	20.5 ± 4.0b	...	4.5 ± 0.9	...	1.5	5.8
SPT0538-50 A						19.8 ± 4.6	0.52 ± 0.12	3.0c	2.4c
SPT0538-50 B						21.9 ± 3.7	1.61 ± 0.33	1.4c	0.12c

Notes. Column 1: SPT source name. Column 2: background source redshift. Column 3: lens redshift. Column 4: Einstein radius (arcsec). Column 5: lens mass, interior to r_E (10¹¹ M_☉). Column 6: lens ellipticity. Column 7: total magnification of background source. Column 8: source radius, determined as the half-width at half-maximum for the Gaussian source component in the model fit (kpc). Column 9: intrinsic far infrared luminosity (10¹² L_☉). Column 10: intrinsic source flux (luminosity per area; 10¹² L_☉/kpc²). Column 11: intrinsic 1.4 mm flux density, obtained as the SPT flux density divided by μ (mJy). Column 12: the same as (11), but for the ALMA 350 GHz (860 μm) flux density. The parameter uncertainties do not include a contribution from the cosmological parameters. An additional 4% uncertainty in mass is found by marginalizing over the WMAP7 parameter Markov chains. ^aAssuming that the lens is located at z = 0.8, see Section 4. ^bTotal magnification of the two source components, see Section 4. ^cDerived assuming that L_FIR is divided between components in the same ratio as the flux density in the model.

Download table as:  ASCIITypeset image

Table 1 shows source and lens redshifts, as available, for these four SPT sources. Redshifts for three of the dusty sources were obtained in another ALMA Cycle 0 project (2011.0.00957.S; PI: A. Weiss) through a spectral line survey conducted using the Band 3 (84–116 GHz) receivers. The complete results of this survey will be reported in Vieira et al. (2013) and Weiß et al. (2013). For each of these sources, multiple high-significance lines are detected in the Band 3 spectral scan, providing unambiguous redshifts. In the case of SPT0538-50, a redshift was determined from a combination of millimeter-wavelength and optical spectroscopy, as described in Greve et al. (2012), and the source was not included in the ALMA redshift search proposal.

The combination of NIR pre-imaging and submillimeter interferometric observations described below provided the basis for ground-based spectroscopic observations of the putative lens galaxies. The fields of SPT0346-52, SPT0418-47, and SPT0529-54 were targeted first with R-band pre-imaging, then with multi-object masks, using the Mask Exchange Unit of FORS2 (Appenzeller et al. 1998) on the Very Large Telescope (VLT) UT1, and exposing them for 3 × 900 s integrations each. The data were collected during 2012 February and March, in Service Mode (ESO program ID 088.A-0902) under an average seeing of ∼11. The galaxies were observed with slits of 1'' in width using the G300 V grism, yielding a velocity resolution of ∼650 km s⁻¹ or ∼13 Å, sampled at ∼3.3 Å pixel⁻¹. SPT0538-50 was observed with the XSHOOTER echelle spectrograph (Vernet et al. 2011) in longslit mode. For these observations, the slit widths were 1'' (UV-B), 09 (VIS-R), and 09 (NIR); the corresponding resolving powers are R = 5100, 8800, and 5600, sampled with 3.2, 3.0, and 4.0 wavelength bins, respectively (after on-chip binning by a factor of two of the UV-B and VIS-R detectors). 6 × 900 s integrations were obtained in 2010 February, however, half of these were taken under worse conditions and we use only the better three integrations in the spectrum presented here.

The observations were prepared and the data reduced using the standard ESO pipeline,³⁵ performing bias and flat corrections, background subtraction, fringe correction, registration and combination, wavelength calibration, and one-dimensional spectral extractions. Spectra are shown in Figure 2. With the exception of SPT0346-52, redshifts for the lensing galaxies were measured by fitting Gaussians to the Ca K+H absorption lines.

Zoom In Zoom Out Reset image size

The small primary beam size of the 12 m ALMA antennas at 350 GHz required that we improve upon the ∼10'' positional uncertainty of the SPT detections before proposing ALMA observations. These targets were initially followed up with the Submillimeter Array (SMA; Blundell 2004; Ho et al. 2004), or the Australia Telescope Compact Array³⁶ (ATCA; SPT0529-54 only). The SMA observations (M. Bothwell et al., in preparation) had a typical angular resolution of 3''× 12'' because of the low declination of the sources, but those data gave some indication that the submillimeter emission was resolved. The 100 GHz ATCA detection of SPT0529-54 was of low significance and did not reveal any ringlike structure.

Deep near-infrared (NIR) imaging data were acquired from the Southern Astrophysical Research (SOAR) Telescope (OSIRIS; Depoy et al. 1993), the Very Large Telescope (VLT ISAAC; Moorwood et al. 1998), and the Hubble Space Telescope (HST), as part of followup programs to examine lens properties and identify rest-frame optical emission from the background sources.

The luminosities in Table 1 are derived from photometric observations at millimeter to submillimeter wavelengths. In addition to SPT photometry (1.4 and 2 mm), we use 870 μm observations from LABOCA (Siringo et al. 2009) on the APEX telescope (in Max Planck time) and Herschel-SPIRE photometry at 250, 350, and 500 μm. The LABOCA data were analyzed according to the procedure described in Greve et al. (2012). The Herschel data are reduced as described in Weiß et al. (2013). All photometric data, including SPT measurements and the data previously reported in Greve et al. (2012) for SPT0529-54 and SPT0538-50, are provided in Weiß et al. (2013). The LABOCA measurements of the total 870 μm flux density agree with the total flux density in the modeled ALMA data (described below) to within calibration uncertainties, typically 10%. Luminosities are calculated using a graybody model like that described in Greve et al. (2012), though in the present case, the τ = 1 wavelength is a parameter of the fit. The effect of freeing this parameter is to broaden the peak of the spectral energy distribution (SED) as needed to match the short-wavelength Herschel data points.

The ALMA interferometer measures visibilities, Fourier components of the sky intensity distribution across a two-dimensional range of spatial frequencies, rather than directly imaging the emission. To properly compare these data with a source model, we must perform our analysis in the visibility plane, where the measurement and its noise are well understood (e.g., Bussman et al. 2012). Past techniques for modeling interferometric lens observations (e.g., Wucknitz et al. 2004) have generally operated on reconstructed images that are subject to difficult-to-model biases and noise properties. Furthermore, residual errors in calibration in the visibility data, such as those arising from imperfect knowledge of the antenna positions or uncompensated atmospheric delay, are often corrected as part of the imaging process through an iterative clean/self-calibration technique (Cornwell et al. 1999). However, the inclusion of the cleaning step in the determination of these corrections, which are then applied to the visibility data themselves, changes the data in ways that are not easily included in the modeling uncertainties. Here, we describe a visibility-based lens modeling technique that simultaneously determines these self-calibration phases so that we incorporate the full range of uncertainty present in the measurements.

We model the lenses by generating model lensed images that are subjected to simulated observations and compared to the data. The source is assumed to have a symmetric Gaussian light profile with four free parameters: flux density F_S, radius R_S, and positional offsets from the lens center X_S, Y_S. The lens is modeled as a Singular Isothermal Ellipsoid (SIE) with five free parameters: mass inside the Einstein radius M_L, ellipticity _L, orientation angle θ_L (east of north), and position X_L, Y_L. The SIE profile has been shown to be a good approximation to galaxy-scale density profiles (Treu & Koopmans 2004; Koopmans et al. 2006, 2009). Since lensed image positions provide very precise measurements of the projected mass interior to the images, we report this robustly measured quantity in Table 1, the total mass (M_L) inside the Einstein radius (r_E) of the lensing galaxies. However, the total halo masses associated with these galaxies will be much higher. Magnification, μ, is calculated as the ratio of the total lensed to unlensed flux.

Given a set of lens and source parameters, we make a high-resolution image of the lensed source, which we pad with zeroes and Fourier transform to generate model visibilities. For each ALMA visibility, we interpolate the Fourier transformed image to the ALMA u, v coordinates. We also correct for the antenna primary beam attenuation by multiplying the sky model images by the primary beam pattern before sampling the Fourier modes. We use a symmetric Gaussian with FWHM of 18'' for the primary beam, which leads to minimal attenuation for these well-centered sources.

The agreement between data and model visibilities is determined by calculating the χ² between them, which requires an estimate of the visibility noise. Because we observe strong sources, the visibility scatter has a contribution from the sky signal and cannot be used to determine the noise level. Instead, we derive noise levels by measuring visibility scatter after differencing visibilities that are adjacent in time for the same baseline/polarization/IF, which has the effect of removing all sky signal. This gives results that are identical to those found by scaling the visibility noise to obtain a reduced χ² of unity for the best-fit models. We explore the model parameter space using a Markov Chain Monte Carlo (MCMC) method with Metropolis–Hastings sampling. An example of the parameter degeneracies in these fits is shown in Figure 3.

Zoom In Zoom Out Reset image size

An additional complication of the modeling procedure is the presence of uncorrected antenna-based phase errors in the data. Observations of test sources, quasars with positions referenced to the International Celestial Reference Frame (ICRF), show significant residual phase errors after the primary phase calibration step. For simple source structures, self-calibration using CLEAN components as input for the phase correction is a standard procedure to improve image fidelity (e.g., Taylor et al. 1999). However, in this application, it would add significant complication to model the clean/self-cal process as part of the fitting, and leaving the phase errors uncorrected can significantly bias the model parameters because of the strong sensitivity to image flux ratios in the model. We have therefore developed a procedure to determine the self-cal phases as part of the model fitting. We optimize the χ² for each step in the Markov chain by adjusting the N − 1 antenna phases. We find that the resulting phases vary little over the chain and closely resemble those found for nearby point sources added to the tracks to test the calibration and astrometry, giving us confidence in this method. Additional details of the method and simulations of its effectiveness are provided in the Appendix.

The model parameter space is complex, with the possibility of multiple isolated minima separated by high barriers in χ². To decrease the possibility of missing important minima in the posterior, we search the space more broadly by "tempering," which is similar to the simulated annealing method (e.g., Press et al. 2007). A control parameter, T (the analog of temperature), is introduced to flatten the posterior surface and to make the minima more accessible. This is achieved by raising the posterior to the power 1/T (with T > 1).

Interferometric phase calibration procedures generally leave some residual phase errors due to imperfect baseline solutions, uncompensated atmospheric delays, or other effects. The magnitude of these errors depends on many factors, such as calibration interval and calibrator-source separation, and their importance depends on the signal to noise ratio in the data and the complexity of the imaging task. In the present application, small phase errors can redistribute flux between lensed images, and without properly accounting for such effects in our lens modeling, the derived model parameter distributions may be significantly in error.

A standard procedure to correct antenna-based phase errors is self-calibration. Using a source model derived from images of the corrupted visibilities, phase corrections are derived, the source is re-imaged, and another iteration can be made using a new source model produced from the phase-corrected visibilities. For our present purposes, the most significant disadvantage of this phase correction scheme is that the uncertainties associated with the source model against which the data are self-calibrated are not included in the lens modeling, and any structural errors introduced by noise or phase errors may become a permanent part of the lens model. Here, we propose a method to incorporate the phase correction into the lens modeling procedure, using the lensed structure as the source model for self calibration.

To implement this self calibration technique, we use a perturbative approach in which we assume that the current data is equal to the model plus an antenna based phase corruption. We can write the expression for χ² for the visibility phases between the data and model as

$$\begin{equation} \chi ^2 = \left[\delta \Phi _i +\delta \phi _k \frac{\partial \Phi _i}{\partial \phi _k}\right] C^{-1} \left[\delta \Phi _i + \frac{\partial \Phi _i}{\partial \phi _k} \delta \phi _k \right] , \end{equation} \tag{ A1 }$$

where Φ_i is the phase of the i'th visibility, δΦ_i is the phase difference between the model and the data, and δϕ_k is the phase delay in the k'th antenna. ∂Φ_i/∂ϕ_k is a Jacobian matrix containing the gradients of the observed visibilities with respect to changes in antenna phases. For N antennas and M visibilities, ∂Φ/∂ϕ is an N × M matrix whose ik'th element is 1 if the first antenna of i'th visibility is k, −1 if the second antenna of i'th visibility is k, and zero otherwise. C is a diagonal covariance matrix whose nonzero elements are approximately (σ_xy/|V|)², where σ_xy is the rms error on the real or imaginary part of each visibility and |V| is the visibility amplitude. Note that this approximation to the covariance matrix is only valid in the limit of high signal-to-noise data (Wrobel & Walker 1999).

To minimize χ², we set its derivative with respect to antenna phases to zero (∂χ²/∂ϕ_i = 0). This allows us to write the antenna phase offsets as

$$\begin{equation} \delta \phi _l = -(F^{-1})_{lk} \frac{\partial \Phi }{\partial \phi } \, \, (C^{-1})_{ji} \, \, \delta \Phi _i, \end{equation} \tag{ A2 }$$

where F is the Fisher matrix calculated as

$$\begin{equation} F_{ij}=\frac{\partial \Phi _k}{\partial \phi _i} (C^{-1})_{kl} \frac{\partial \Phi _l}{\partial \phi _j}. \end{equation} \tag{ A3 }$$

At every iteration of the MCMC code, the value of χ² is minimized for the postulated model by deriving calibration phases using Equation (A2). The resulting χ² is used to evaluate the likelihood and determine the next link in the chain, thereby incorporating the uncertain phase correction in the parameter exploration. In the limit of an intrinsically Gaussian distribution for these phase parameters, this is equivalent to marginalizing over the phases. In Figure 7, we compare the results of this procedure with the standard self-calibration based on CLEAN components, and the raw data. The improvement in source subtraction is significant, even compared to the standard CLEAN procedure.

Zoom In Zoom Out Reset image size

A second test of the simultaneous fitting of the lensed emission and antenna phases is shown in Figure 8. Simulated observations of a typical lens model were created with realistic noise levels, and antenna phase errors were added to these visibilities. The MCMC fitting algorithm was applied to these multiple realizations of the same source to verify that the antenna phase errors are recovered. The excellent agreement between input and recovered phase demonstrates that the simultaneous fitting of the lens model and antenna phases does not bias the antenna phase measurement, despite the complicated source structure of the lens models.

Zoom In Zoom Out Reset image size

The models in Figure 4 were derived using the compact configuration data that were available at the time of submission. Higher resolution data were delivered later, and permit a direct comparison of the model based on the low-resolution data with the higher resolution observations of SPT0346-52. Figure 9 shows the predictions of the best fit model (from fitting to the compact data, presented earlier in this work) for the uv-coverage of the extended data (red contours). The black contours show the extended data. The contours demonstrate a high degree of agreement between the predictions and the new observations. A lens model for the extended configuration data shows a consistent model, resulting in magnification (from modeling the extended data alone) of 5.26 ± 0.12 in agreement with the magnification derived from the compact data.

Zoom In Zoom Out Reset image size