Identifying high-redshift GRBs with RATIR
O. M. Littlejohns, N. R. Butler, A. Cucchiara, A. M. Watson, A. S. Kutyrev, W. H. Lee, M. G. Richer, C. R. Klein, O. D. Fox, J. X. Prochaska, J. S. Bloom, E. Troja, E. Ramirez-Ruiz, J. A. de Diego, L. Georgiev, J. González, C. G. Román-Zúñiga, N. Gehrels, H. Moseley
aa r X i v : . [ a s t r o - ph . H E ] A p r Draft version September 16, 2018
Preprint typeset using L A TEX style emulateapj v. 5/2/11
IDENTIFYING HIGH-REDSHIFT GRBS WITH RATIR
O. M. Littlejohns , N. R. Butler , A. Cucchiara , A. M. Watson , A. S. Kutyrev , W. H. Lee , M. G. Richer ,C. R. Klein , O. D. Fox , J. X. Prochaska , J. S. Bloom , E. Troja , E. Ramirez-Ruiz , J. A. de Diego ,L. Georgiev , J. Gonz´alez , C. G. Rom´an-Z´u˜niga , N. Gehrels and H. Moseley Draft version September 16, 2018
ABSTRACTWe present a template fitting algorithm for determining photometric redshifts, z phot , of candidatehigh-redshift gamma-ray bursts (GRBs). Using afterglow photometry, obtained by the ReionizationAnd Transients InfraRed (RATIR) camera, this algorithm accounts for the intrinsic GRB afterglowspectral energy distribution (SED), host dust extinction and the effect of neutral hydrogen (localand cosmological) along the line of sight. We present the results obtained by this algorithm andRATIR photometry of GRB 130606A, finding a range of best fit solutions 5 . < z phot < . z phot in the ranges 4 < z phot . < z phot <
10 and can robustlydetermine when z phot >
4. Further testing highlights the required caution in cases of highly dustextincted host galaxies. These tests also show that our algorithm does not erroneously find z phot < z sim >
4, thereby minimizing false negatives and allowing us to rapidly identify all potentialhigh-redshift events.
Subject headings:
Gamma-ray bursts: general, gamma-ray bursts: specific (GRB 130606A), tech-niques: photometric INTRODUCTION
Gamma-ray bursts (GRBs) (Gehrels et al. 2009;M´esz´aros 2006) are bright objects that emit across theentire electromagnetic spectrum and therefore can beseen to high-redshift (Lamb & Reichart 2000). Thisallows observers to identify the locations of high-redshift, faint, host galaxy population otherwise missedby the current and future magnitude limited surveys(Tanvir et al. 2012). Such discoveries will allow us to un-veil the properties of early star formation (up to z ∼ < z <
10) have been identified in-cluding GRB 050904 (Totani et al. 2006), GRB 080913(Greiner et al. 2009a), GRB 090423 (Tanvir et al. 2009)and GRB 090429B (Cucchiara et al. 2011). However,it is only in a more recent example, GRB 130606A at z = 5 .
913 (Chornock et al. 2013), that high signal tonoise ratio spectrum has been obtained of a high-redshiftGRB afterglow. Using a method similar to that employedfor quasars (Fan et al. 2006; Becker et al. 2001), the de-tection of Ly β and Ly γ transmission in this spectrumallowed Chornock et al. (2013) to identify that the inter-galactic medium (IGM) is mostly ionized at this redshift, School of Earth & Space Exploration, Arizona State Univer-sity, AZ 85287, USA NASA Postdoctoral Program Fellow, Goddard Space FlightCenter, Greenbelt, MD 20771, USA Instituto de Astronom´ıa, Universidad Nacional Aut´onoma deM´exico, Apartado Postal 70-264, 04510 M´exico, D. F., M´exico NASA, Goddard Space Flight Center, Greenbelt, MD 20771,USA Astronomy Department, University of California, Berkeley,CA 94720-7450, USA Department of Astronomy and Astrophysics, UCO/Lick Ob-servatory, University of California, 1156 High Street, Santa Cruz,CA 95064, USA and that the epoch of reionization must have occurredearlier in the history of the Universe.Castro-Tirado et al. (2013) also present broad-band photometric and spectroscopic observations ofGRB 130606A, finding a lower HI column density of N HI = 7 . × cm − . Totani et al. (2013) consider apossible intervening damped Lyman- α (DLA) system at z DLA = 5 .
806 as a contributor to the absorption via HIgas. Considerations of the required silicon abundance insuch a DLA system lead Totani et al. (2013) to insteadattribute the residuals in their host-only model to diffuseIGM along the line of sight at 5 . < z IGM < . z phot , estimates for high-redshift dropout candidates. Whilst z phot = 9 . +0 . − . Littlejohns et al.(90% likelihood range) was later obtained, additionalrapid-response from GRB-dedicated NIR multiple filterinstruments would have assisted in providing a measureof z phot and allowing large aperture telescopes to beginobservations at the earliest possible epoch after the ini-tial trigger.The variable z phot is estimated by fitting templates tothe optical and NIR spectral energy distributions (SEDs;see for example Kr¨uhler et al. 2011 and Curran et al.2008) observed for the afterglow of the GRB, which hasan intrinsic synchrotron spectrum (Granot & Sari 2002;Sari et al. 1998). Through this template fitting methodthe redshift of a GRB, the host galaxy dust extinctionand the optical/NIR spectral index of the GRB after-glow can be estimated. Fitting a host extinction lawallows a study of the prevalence of dust in early galaxiesalong GRB sight lines (Zafar et al. 2011), whilst the lo-cal spectral index can be compared to that obtained inother energy ranges to determine key parameters of theunderlying synchrotron spectrum (Perley et al. 2013).To determine z phot , Kr¨uhler et al. (2011) make ex-tensive use of data from the Gamma-Ray BurstOptical/Near-Infrared Detector (GROND; Greiner et al.2008). GROND simultaneously images in seven fil-ters: g’ , r’ , i’ , z’ , J , H and K S , providing a broadoptical and NIR SED. Such an SED was used to de-termine the z phot = 4 . ± .
15 for GRB 080916C(Greiner et al. 2009b) and another indicated helped iden-tify that GRB 080913 was a high redshift event ( z =6 . ± . z phot isthe Reionization And Transients Infra-Red (RATIR;Butler et al. 2012) camera, which is mounted on the1.5-meter Harold L. Johnson telescope of the Ob-servatorio Astron´omico Nacional on Sierra San Pe-dro M´artir in Baja California, Mexico. This facil-ity, which became fully operational in December 2012,conducts autonomous observations (Watson et al. 2012;Klein et al. 2012) of GRB triggers from the Swift satellite(Gehrels et al. 2004). Since RATIR obtains simultaneousphotometry in r , i , Z , Y , J and H , it is an excellent in-strument for estimating z phot for high-redshift GRBs andas such allows us to optimize spectroscopic observationswith highly oversubscribed large telescopes.RATIR and GROND operate in similar manners, bothwith best responses of order a few minutes after the ini-tial Swift trigger time (Butler et al. 2013c; Updike et al.2010). The filters employed by both instruments aresimilar, but not identical. GROND has a broader spec-tral coverage, extending to both higher and lower wave-lengths. The lower wavelength coverage aids in the iden-tification of 2 . < z phot <
4, although RATIR is specifi-cally designed to target objects in the high-redshift Uni-verse. While GROND also has a K S -band filter, extend-ing coverage to higher redshifts, RATIR contains a Y -band filter that GROND lacks. With such functionality,RATIR is better able to provide more precise estimatesof z phot when a GRB occurs 7 < z <
8. RATIR isalso 100% time dedicated to GRB follow-up. Perhapsthe most important distinction between the two instru-ments is in their locations. RATIR is situated at latitude φ = +31 ◦ ′ ′′ , while GROND is mounted on the 2.2- meter MPI/ESO telescope at La Silla Observatory, Chile,with a latitude φ = − ◦ ′ ′′ , meaning the both forma complimentary pair of instruments routinely accessingdifferent parts of the sky.This work describes the development of our templatefitting routine. In § § Swift and RATIRobservations of GRB 130606A, before in § §
5. Finally, we presentour conclusions in § METHOD
To measure z phot we adopt a template fitting method-ology similar to that of Kr¨uhler et al. (2011), who usecombined data sets from GROND and both the Ultravi-olet/Optical Telescope (UVOT; Roming et al. 2005) andX-ray Telescope (XRT; Burrows et al. 2005) on boardthe Swift satellite. Another similar methodology is thatof Curran et al. (2008), although the algorithm we im-plement is significantly less computationally expensive.To this end we use models of the underlying physicalspectrum of the source and the mechanisms by whichthis spectrum is altered before reaching an observer atEarth. These processes are dust absorption from thehost galaxy, attenuation from the Intergalactic Medium(IGM) and the response of the RATIR filters to incomingphotons.Unlike Kr¨uhler et al. (2011) our algorithm does notconsider UVOT photometry, as only approximately 26%of
Swift
GRBs having UVOT detections, compared tothe ∼
60% of bursts detected by ground-based facilities(Roming et al. 2009). Following Schady et al. (2012), wealso favor the extinction laws of Pei (1992), instead ofthe more general “Drude” model (Liang & Li 2009) usedby Kr¨uhler et al. (2011). This was done in an attemptto limit the number of parameters being fitted as the“Drude” model has four parameters compared to the one( A V ) used in the standard Pei (1992) extinction laws. Intrinsic GRB spectrum
The emission of GRBs is attributed to synchrotronemission (Sari et al. 1998), which produces a spectrumconsisting of a series of broken power-laws. The loca-tion of the synchrotron cooling break, usually found be-tween the optical and X-ray regimes, evolves with timein a manner determined by the nature of the circum-burst medium (Granot & Sari 2002). The passage ofsuch a cooling break through the optical and NIR bandswould manifest itself through an achromatic temporalbreak in the GRB light curve, which we do not observefor GRB 130606A. Due to this, and to avoid the compli-cation of additional fitting parameters, we assume thatthe regime in which RATIR observes was entirely withinonly one of the power-law segments, thus: F ( λ ) = F (cid:18) λλ (cid:19) β opt , (1)where F ( λ ) is flux as a function of wavelength, λ , F is the normalization corresponding to the flux at a desig-nated wavelength, λ , and β opt is the local spectral indexdentifying high-redshift GRBs with RATIR 3of the power-law. Both λ and λ are in the rest frameof the radiating outflow. Traditionally, the spectrum ofa GRB is represented as a function of frequency ratherthat wavelength, however the expression in Equation 1 isentirely equivalent and easier to combine with the dustextinction and hydrogen absorption models detailed in § § Host dust absorption
Once emitted, the first external effect to act upon theGRB radiation is absorption by dust in the medium ofthe host galaxy along the GRB sight line. Most GRBhost galaxies detected in the optical and NIR regimes ex-hibit low amounts of dust extinction along the GRB lineof sight (Starling et al. 2007; Kann et al. 2006) with themajority of SED models for a sample covering 0 . < z < . A V .
2. NIR data have revealed a subset( ∼ A V & . z . a priori , the tem-plate algorithm implements models for the MW, LMCand Small Magellanic Clouds (SMC). We also consider acase with no host dust absorption.To calculate dust absorption as a function of wave-length, we use the empirical extinction curves presentedin Table 1 of Pei (1992) and the corresponding values of R V , where R V ≡ A V /E B − V and A V is the extinctionin magnitudes of the V -band. E B − V = A B − A V givesthe difference in extinction between the B - and V -bands.Both A V and E B − V are defined in terms of rest-frame B and V .We dust extinction at each wavelength given in Table1 of Pei (1992), using Equation 2, below. We then re-sample the extinction law on to a finer grid. A λ = A V (cid:18) R V E λ − V E B − V + 1 (cid:19) . (2) IGM attenuation
The algorithm then calculates attenuation from the in-tervening IGM due to absorption from clouds of neu-tral hydrogen along the line of sight (Gunn & Peterson1965). We adopt a similar methodology to that usedin the hyperz software (see also Bolzonella et al. 2000)by estimating the reduction of flux from Ly α and Ly β absorption using Equation 3, below, which is based onEquation 17 of Madau (1995): http://webast.ast.obs-mip.fr/hyperz/ D ( λ ) = , λ rest > λA × z ) R z ) exp " − A (cid:18) λ obs λα (cid:19) . dλ obs , < λ rest < λB × z ) R z ) exp − P j =3 , Aj (cid:18) λ obs λj (cid:19) . dλ obs , < λ rest < , λ rest <
912 ˚A. (3)
The second line of Equation 3 states the contributionfrom Ly α , whilst the third line totals that from Ly β ,Ly γ and Ly δ . ∆ λ A = 120 (1 + z ) ˚A, ∆ λ B = 95 (1 + z )˚A, A = 3 . × − is the coefficient of Ly α absorptionrelative to the Ly α forest contribution. A = 1 . × − , A = 1 . × − and A = 9 . × − are the coefficientsfor Ly β , Ly γ and Ly δ , respectively.To determine which filters require the application of D ( λ ), we calculate R IGM using Equation 4: R IGM = λ R λ λF ( λ ) A λ D ( λ ) T ( λ ) dλ λ R λ λF ( λ ) A λ T ( λ ) dλ (4)This ratio considers how many source counts passthrough the filter both with and without IGM absorp-tion, whilst weighting the flux at each wavelength by thefilter transmission curve, T ( λ ). A value of R IGM = 1indicates that there is no IGM absorption in the band,and so the flux evolves slowly across the bandpass inquestion. Similarly, when R IGM = 0 all flux within thefilter is absorbed by neutral hydrogen and is thereforenot evolving across the band. In these instances we canneglect the shape of the filter transmission curves, andtherefore use the first line of Equation 5. F filter = λ R λ F ( λ ) A λ D ( λ ) dλ λ R λ dλ , R IGM = 0 , R IGM λ R λ F ( λ ) A λ dλ λ R λ dλ , < R IGM <
1. (5)When the effects of neutral hydrogen absorption liewithin a particular filter, there is a sharp transitionwithin its wavelength coverage, 0 < R
IGM <
1, and thefilter requires a more detailed treatment. In such a filterthere is a fraction of wavelength coverage which expe-riences suppression of GRB flux. We need to calculatethis fraction so we can correctly predict the total aver-age magnitude across the filter. We therefore apply afactor of R IGM to the filter magnitude, which correctlyaccounts for the fraction of flux that successfully reachesthe telescope and passes through the filter in question.This is shown in the second line of Equation 5.
Fitting and prior probabilities
To implement the model described we used the amoeba fitting routine within idl as well as the idl
As- Littlejohns et al.tronomy Library (Landsman 1993). amoeba is a robustroutine that can avoid local minima in parameter space.To avoid regions of the parameter space containingunphysical values of β opt we imposed a prior probabil-ity distribution upon the spectral index in the opticaland NIR regime. This was implemented in the man-ner described in § β opt = β X ), (ii) there is a cooling break be-tween the two regimes in the fast cooling regime (imply-ing β opt = ), (iii) a cooling break in the slow coolingregime ( β opt = β X − ) (Granot & Sari 2002). The cool-ing break occurs at the synchrotron frequency emittedby electron with a Lorentz factor that causes the coolingtimescale from radiation, t cool , to equal the dynamicaltime of the system, t dyn (Granot et al. 2000). In thefast cooling regime t cool < t dyn and electrons cool via ra-diative losses faster than the dynamical timescale. Theslow cooling regime is the reverse, in which the majorityof electrons are unable to cool within t dyn .These three regimes for β opt lead to a prior probability p ( β opt | I ), given in full in Equation 6. I is the known in-formation about the parameter being fitted (in our priorthis is β opt ). p ( β opt | I ) = G ( β opt , β X , ∆ β X ) + G (cid:0) β opt , , ∆ β X (cid:1) + G (cid:0) β opt , β X − , ∆ β X (cid:1) . (6) Schady et al. (2012) discovered that approximately 50%(25/49) of their sample of X-ray and NIR detected longGRB SEDs can be well fitted assuming β opt = β X , sowe weighted p ( β opt | I ) accordingly. With no further in-formation on whether the burst is in the fast or slowcooling regime, these two options were equally weightedto be 25% of the total distribution each.In this instance our prior is the three expected valuesof β opt already discussed. Effectively, the prior proba-bility weights the parameter space and reduces the totalviable range a parameter can take. We include realisticestimates of the Gaussian width, σ , for each of the dis-tributions in p ( β opt | I ) and so ensure β opt is not forcedto take the exact value of one of three states discussed. Error estimation in z phot The primary output of the template fitting routines isa value for z phot . To obtain an error in this value we con-sider a uniform grid of values for redshift, z grid . This gridextends over the full range of potential redshifts that canbe fitted (see § z grid we fit β opt , A V (for the no, MW, LMC and SMC type dust models) andthe normalization to the intrinsic spectrum, N , using amoeba whilst holding z grid constant. To optimize thesearch of the parameter space using amoeba we invoke p ( β opt | I ). To do so we calculate p ( β opt | I ) within the fit-ting routine and include it as an additional term of ourfit statistic, χ , as shown in Equation 7. The first termof Equation 7 corresponds to the standard χ value. χ = χ − p ( β opt | I )) . (7)By calculating χ we force amoeba to find the bestfit solution at each point of z grid which accounts for our knowledge of β X . We can then find a measure of theprior weighted probability distribution, p ( D | θI ), where θ corresponds to the full set of parameters, directly from χ at each value of z grid using Equation 8. The value of z grid where p ( D | θI ) peaks is our estimate of z phot . p ( D | θI ) ∝ exp (cid:18) − χ (cid:19) . (8)To provide an error estimate on z phot we then find thenarrowest range of z grid containing 90% and 99.73% ofthe total weighted probability distribution. The formercorresponds to the 90% confidence interval, whilst thelatter would correspond to the 3 σ confidence interval ofa Gaussian distribution. Both were obtained for all fourextinction models discussed in § OBSERVATIONS OF GRB 130606A
The Burst Alert Telescope (BAT; Barthelmy et al.2005) on board
Swift triggered on GRB 130606A (
Swift trigger 557589) at 21:04:39 UT on 2013 June 6 th (Ukwatta et al. 2013). BAT identified prompt structureat the trigger time and later, brighter γ -ray emission atapproximately 150 seconds after the initial trigger time. Swift /XRT detected an uncatalogued fading X-raysource, providing a more accurate position for ground-based follow-up with narrow field instruments.The enhanced-XRT position for GRB 130606A wasfound to be RA(J2000) = 16 h m . s
12, Dec(J2000) =+29 ◦ ′ . ′′ . ′′ Swift /UVOT to astrometrically correct the positions ofsources in the XRT field of view.Chornock et al. (2013) observed the afterglow ofGRB 130606A with both the Blue Channel spectro-graph (Schmidt et al. 1989) on the Multiple Mirror Tele-scope (MMT) and the Gemini Multi-Object Spectro-graph (GMOS; Hook et al. 2004) on the Gemini Northtelescope. These observations had midpoints of 7.68and 13.1 hours after the initial BAT trigger, for theBlue Channel spectrograph and Gemini-N/GMOS re-spectively. From this rapid follow-up with large aperturefacilities, Chornock et al. (2013) obtained high signal tonoise ratio spectra allowing them to measure z = 5 . r - and i -bands and 1.85 hours in the Z -, Y -, J - and H -bands (Butler et al. 2013b). The observations for allsix photometric bands from this first night are shown inFigure 1.A second epoch of observations were taken the follow-ing night between 31.14 and 37.86 hours after the Swift trigger (Butler et al. 2013a). With 4.96 hours of expo-sure time in the r - and i -bands both yielding upper lim-its, showing that the source had faded in the i -band. Theburst remained bright enough in each of the NIR filtersto allow detections in all four bands with a total expo-sure time of 2.07 hours. In all four filters the burst hadfaded by more than two magnitudes, as shown in Table1.dentifying high-redshift GRBs with RATIR 5 Figure 1.
Photometry obtained for the field of GRB 130606A. Top left: r -band, top right: i -band, middle left: Z -band, middle right: Y -band, bottom left: J -band and bottom right: H -band. Each image is the sum of all observations taken on the night of the Swift /BATtrigger. The red circle in each panel is the refined
Swift /XRT error circle and the blue circle indicates the RATIR identified counterpart.All images are approximately 40 ′′ × ′′ in scale. Littlejohns et al.The RATIR camera consists of two optical and two in-frared cameras (Butler et al. 2012), allowing for simulta-neous image capture in four bands ( riZJ or riY H ). Splitfilters immediately above the infra-red detectors allow fornear-simultaneous images in six bands by dithering thesource across the infra-red detector focal plane. We em-ploy a dither pattern that allows for sequential exposurein ZJ followed by Y H . Dithering also allows for subtrac-tion of the sky background and the detector dark signalat the source position. We capture 80s exposure framesin ri and 67s exposure frames in ZY JH due to additionaloverhead. We apply the same image reduction pipeline –with twilight flat division and bias subtraction routineswritten in python and using astrometry.net (Lang et al.2010) for image alignment and swarp (Bertin 2010) forimage co-addition – to data taken from each camera.Photometry is calculated by running sextractor (Bertin & Arnouts 1996) over individual science framesand mosaics over a range of apertures between 2 and30 pixels diameter. A weighted average of the flux inthis set of apertures for all stars in a given field isthen used to construct an annular point-spread-function(PSF). This PSF is then fitted to the annular flux val-ues for each source to optimize signal-to-noise for pointsource photometry. We perform a field-photometric cal-ibration for the riZ bands using the Sloan Digital SkySurvey Date Release 9 (SDSS DR9;Ahn et al. 2012). TheRATIR riZ filters have been shown to be equivalent tothe SDSS filters ( riz ) at the .
3% level (Butler et al.in prep). We calibrate the
Y JH bands relative to theTwo Micron All Sky Survey (2MASS; Skrutskie et al.2006), employing the United Kingdom Infrared Telescope(UKIRT) Wide Field Camera (WFCAM; Hodgkin et al.2009; Casali et al. 2007) to estimate Y from JH for fieldstars. Photometric residuals for these bands are also ap-proximately 3%. MODELING GRB 130606A SED
The spectroscopic observations for GRB 130606A allowus to assess the accuracy of our template fitting routineusing this burst as a test case. The standard operationof our algorithm requires input RATIR photometry andan estimation of β X . The former is taken from the ini-tial RATIR GRB Coordinates Network (GCN) circular,whilst the latter is obtained from the UK Swift
ScienceData Centre (UKSSDC) automated analysis data prod-ucts (Evans et al. 2009).The algorithm then uses a broad 1-dimensional grid inredshift (0 z grid
12) to find z phot . This value pro-vides a rapid and robust estimate of z phot that aids indetermining whether follow-up with large aperture facil-ities is warranted.We then use 2-dimensional grids to refine our estimateof z phot . Two such grids are used for each extinction law.The first fits A V and model normalization, N , at fixedpoints in z grid β opt parameter space over the range 0 z grid
12. By determining the region of this parameterspace with the smallest total χ , we then produce asecond grid which focuses on the best fit solution with ahigher resolution in both z grid and β opt .In practice, the results from the finer, 2-dimensionalgrid provide a more precise estimate of z phot . However, when coordinating follow-up from large aperture facil-ities, those produced from the 1-dimensional grid z grid provide sufficiently robust estimates such that further ob-servations can be requested at the earliest possible epochafter the initial Swift /BAT trigger.For GRB 130606A we used RATIR photometry ob-tained between 7.38 and 7.79 hours after the
Swift /BATtrigger time (Butler et al. 2013b). As the RATIR imagereduction pipeline runs on data as it is available fromthe instrument, it was already possible at this epoch toidentify GRB 130606A as an r -band dropout candidate.We also obtained β X = 0 . +0 . − . from the UKSSDC.The best fit for each of the four extinction laws, thethree standard templates from Pei (1992) and a fit wherewe assumed no host extinction, are detailed in Table2. The corresponding templates are plotted with theRATIR photometry in Figure 2. It is important to notethat we only include dust extinction from the host galaxy.For GRB 130606A, there is a possibility of an interveningsystem, with its own unknown dust content (Totani et al.2013). If present, z DLA = 5 .
8, meaning that we would beunable to distinguish between the dust in the host andany potential dust in the DLA. Even should the dustcontent of the DLA be high, our value of z phot remainsrobust due to the small difference in redshift better thenhost galaxy and DLA.Table 2 suggests that either the MW or no-dust solu-tion best represents the data. Visual inspection of Fig-ure 2 suggests the LMC and SMC models are of at leastcomparable quality, however, the magnitude comparedto that measured by RATIR is an average across thewavelength coverage of each filter. For the LMC andSMC templates, the average values across the J - and H -bands are more discrepant than those of the MW or nodust models. When using a MW type extinction law, alarge quantity of dust is favored, with A V = 1 .
48. Conse-quently, z phot for the MW template is reduced, as a largequantity of dust contributes to the strong suppression ofintrinsic GRB flux in the r - and i -bands.Figure 3, which contains the probability maps over allallowed ranges of z grid , shows the best fit solutions arelocated in a narrow redshift range at z phot ∼ β opt ∼ β X and correspondsto z phot = 6 .
01. This is in good agreement with thecruder estimate made using the 1-dimensional grid inredshift, which is z phot = 5 . +0 . − . for an SED with nodust extinction (see Table 2).The LMC and SMC probability maps in Figure 4 con-tain similar structure. The first solution is at β opt ∼ β X ,as found when A V = 0. With the inclusion of a non-zero A V this solution is more extended in the z phot axis. Thisis because at marginally lower redshifts a higher A V al-lows absorption from the IGM to contribute to the fluxdeficit seen at shorter wavelengths. This elongation inthe z phot axis is more pronounced in the lower β opt solu-tion seen for both the LMC and SMC extinction laws.Some values of β opt do not exactly coincide with theexpected values resulting from a cooling break. Thesedentifying high-redshift GRBs with RATIR 7 Table 1
RATIR photometry for GRB 130606A. The first two columns denote the start and stop times ( T start and T stop respectively) of the observations on both nights in hours since the BAT trigger. The magnitudesfor each filter are given in the AB magnitude system. T start T stop r i Z Y J H (hours) (hours)7.38 14.19 24.48 ± ± ± ± ± ± > > ± ± ± ± µ m)110100 F λ ( µ J y ) NoneMWLMCSMCr i Z Y J H
Figure 2.
SED templates for each of the fits in Table 2. The colored points correspond to the measured RATIR photometry, with thefilter being marked above each measurement. The black lines indicate the best fits obtained using the specific host extinction laws detailedin the key.
Table 2
Details of the fitted solutions to the SED ofGRB 130606A, stating which host extinctionmodel was used, z phot , its 3 σ error, β opt , A V and χ for each model. These solutions wereobtained using the initial 1-dimensional grid of z grid values.Extinction z phot β opt A V χ /ν None 5.97 +0 . − . +0 . − . +0 . − . +0 . − . result from an enhancement in the prior probability dis-tribution between β opt = β X − . β opt = 0 . A V = 0 extinction law, the β opt ∼ β X isthe better of the two alternative solutions, and indeed afitted value of A V = 0 is retrieved for both the LMC andSMC extinction laws. The top right panel of Figure 4 shows the same pa-rameter space using the MW extinction law. In this in-stance there are four local maxima in the prior weightedprobability map. The template SED for each is shownin Figure 5. The two higher z phot solutions are equiv-alent to those found when A V = 0. There are two ad-ditional solutions, however, at lower redshift. It is theshape of the MW extinction law that allows for thesetwo solutions, as the 2175 ˚A bump can contribute to thesuppression of low wavelength emission. Even with thecontribution from this feature, the required A V is high,with the β opt ∼ β X solution requiring A V = 0 .
73 and thelower β opt solution requiring A V = 1 . χ is lowest for solution 2 ( z phot = 5 . A V =0 . β opt = 0 . β opt prior probability distribution, the effective χ , χ asdefined in Equation 7, is lower for solution 1 ( z phot =5 . A V = 1 . β opt = 0 . z phot β op t E X P (- χ - e ff / ) χ - /2) No dust z phot β op t E X P (- χ - e ff / ) χ - /2) MW dust z phot β op t E X P (- χ - e ff / ) χ - /2) LMC dust z phot β op t E X P (- χ - e ff / ) χ - /2) SMC dust Figure 3. e χ of the fits, where dark regions have the lowest χ . The narrow horizontal panels demonstrate the 1-dimensional probabilitydistributions for z phot , whilst the narrow vertical panels show the 1-dimensional probability distributions of β opt . Top left: no hostextinction, top right: Milky Way host extinction law, bottom left: LMC host extinction law, bottom right: SMC host extinction law. spectroscopic redshift of Chornock et al. (2013), we finda solution with β opt ∼ β X is most consistent with thereported value of z spec = 5 .
913 in all cases.As an independent test to determine which of thedust extinction laws best represents the data, we usedthe generic three component extinction law measuredin Massa & Fitzpatrick (1986) and Fitzpatrick & Massa(1990, 1988, 1986). Applying the prior probability dis-tribution detailed in Reichart (2001) yielded a fit where z phot = 6 . +0 . − . , β opt = 0 . R V = 3 .
21 and A V =0 .
00. Low dust extinction disfavors solutions 1 and 2found using a MW type extinction law (see Table 3),which both had the most discrepant values of z phot whencompared to z spec as reported in Chornock et al. (2013). FURTHER TESTS OF THE ALGORITHM WITHSIMULATED RATIR DATA
Table 3
Details of the multiple solutions to the SED ofGRB 130606A, stating which host extinction model wasused, z phot , β opt , A V and χ for each model. Thesesolutions were obtained using the initial 1-dimensional gridof z grid values.Extinction Solution z phot β % rmopt A V χ /ν None 1 6.01 0.86 0.00 2.55/3None 2 6.09 0.55 0.00 5.24/3MW 1 5.62 0.42 1.46 2.58/2MW 2 5.69 0.86 0.73 2.24/2MW 3 6.01 0.86 0.00 2.56/2MW 4 6.09 0.55 0.00 5.30/2LMC 1 6.01 0.86 0.00 2.55/2LMC 2 5.82 0.44 0.30 4.00/2SMC 1 6.00 0.86 0.00 2.58/2SMC 2 5.93 0.45 0.12 3.43/2 dentifying high-redshift GRBs with RATIR 9 z phot β op t E X P (- χ - e ff / ) χ - /2) No dust z phot β op t E X P (- χ - e ff / ) χ - /2) MW dust z phot β op t E X P (- χ - e ff / ) χ - /2) LMC dust z phot β op t E X P (- χ - e ff / ) χ - /2) SMC dust Figure 4. e χ of the fits, where dark regions have the lowest χ . The narrow horizontal panels demonstratethe 1-dimensional probability distributions for z phot , whilst the narrow vertical panels show the 1-dimensional probability distributions of β opt . Top left: no host extinction, top right: Milky Way host extinction law, bottom left: LMC host extinction law, bottom right: SMChost extinction law. We constructed simulated data to be processed withour algorithm in order to test the regions of the parame-ter space in which robust values of z phot could be recov-ered from RATIR photometry. Specifically, we wanted toestablish the range of GRB redshifts for which we couldestimate z phot , demonstrate whether the algorithm couldcorrectly recognize the input extinction law and considerthe effects of a high amount of dust extinction. The range of z phot The first tests conducted were simple, with simulatedinput parameters A V, sim = 0, β opt , sim = 1 .
02 and 0 02 was chosen as it correspondedto the peak of the Swift distribution of β X (obtained fromthe on-line repository at the UKSSDC, see Evans et al.2007). In this instance we have assumed that the op-tical and X-ray parts of the spectrum lie on the same power-law slope, such that β opt , sim = β X , sim . While oursimulated source does not have dust extinction, we didallow for the fitting routine to determine a value for A V .Figure 6 shows the fitted values of z phot compared to z sim .Figure 6 demonstrates some clear limitations of ouralgorithm when applied to RATIR data. First is the in-ability to measure z phot . 4. This is physically motivatedby the observed wavelength of Ly α . At z phot = 4, the ob-served wavelength of Ly α is ∼ r -band. The gray region 3 σ error region shows that thealgorithm cannot discern the precise value of z phot , how-ever the lack of an r -band dropout rules out solutionswhere z phot > 4. Similar analyses can be conducted us-ing GROND data sets, which routinely take simultane-0 Littlejohns et al. µ m)110100 F λ ( µ J y ) Solution 1: z phot = 5.62, β opt = 0.42, A V = 1.46Solution 2: z phot = 5.69, β opt = 0.86, A V = 0.73Solution 3: z phot = 6.01, β opt = 0.86, A V = 0.00Solution 4: z phot = 6.09, β opt = 0.55, A V = 0.00r i Z Y J H Figure 5. SED templates for the four possible solutions using a MW extinction law (see top right panel of Figure 4). The colored pointscorrespond to the measured RATIR photometry, with the filter being marked above each measurement. The black lines indicate the bestfits obtained using each of the four solutions detailed in the key and Table 3. z sim z pho t Figure 6. Fitted z phot as a function of the input redshift, z sim ,used to create simulated RATIR magnitudes. The dotted line de-notes equality and the gray region illustrates the 3 σ confidenceregion obtained at each value of z sim . Note A V, sim = 0 in eachsimulation. ous images in seven filters. Using the g’ -band, GRONDis able to determine z phot down to z ≈ . . . z sim . . 75. Thisis due to the dropout feature lying between the RATIR Y - and J -bands. In these instances, the fitting algorithmtends towards a single value of z phot . All templates where z grid is such that the SED features from neutral hydro-gen absorption occur between bands result in an other-wise identical set of fitted parameters. This results infits which are of identical statistical merit across a rangeof z grid . To calculate z phot the fitting routine takes anaverage of this range of z grid , rather than simply takingthe value of z grid that normally corresponds to a uniquebest fit.Despite not having a precise value for z phot when thedropout occurs between the Y - and J -bands, our al-gorithm does still indicate that it lies between the twobands, and thus robustly confirms it is of a high redshift.Such an effect is also present when the dropout occursbetween the lower wavelength filters, however, the gap inwavelength coverage is smaller and so the correspondingranges of z phot in which this occurs are narrower.By comparison, GROND does not include a Y -bandfilter, leaving a large gap in spectral coverage betweenthe z’ - and J -bands. Indeed, Kr¨uhler et al. (2011) alsoencounter a similar effect to that experienced when us-ing our fitting algorithm at 7 . . z sim . . 75. Fromsimulations of a comparable nature to those presentedhere, Kr¨uhler et al. (2011) conclude that the uncertaintyin their estimate of z phot rises from ∆ z phot ∼ . z phot = 6 . z phot ∼ . z phot = 8 . 0. The 3 σ confidence region for z phot as derived from RATIR pho-tometry similarly rises from ∆ z phot = 0 . 08 at z phot = 6 . z phot = 2 . 28 at z phot = 8 . 0. However, this increaseoccurs largely when 7 . . z phot . . 0. Thus, while bothRATIR and GROND can robustly infer when a GRBdentifying high-redshift GRBs with RATIR 11has 6 . . z phot . . 75, only RATIR can provide precisemeasurements when 6 . . z phot . . z sim & . z phot . Aside fromthe dropout feature once more falling between two pho-tometric bands, a single photometric detection in the H -band is insufficient for the number of fitted parametersin our model. Identifying an extinction law In cases where the dust content is high, knowing thetype of extinction law can provide valuable informationregarding the circumburst medium around the GRB.Thus we tested the ability of our algorithm to distin-guish between the three Pei (1992) extinction laws. Tooptimize the effects of the simulated dust, we chose a red-shift, z sim = 3 . 36. This value places the 2175 ˚A bumpat the midpoint between the RATIR Z - and Y -bands al-lowing the SED to most clearly define the shape of thebump. We simulated RATIR photometry at z sim withboth A V = 0 . 25 and A V = 0 . 5, using the three Pei (1992)extinction laws, in turn.As our algorithm is designed for use on early-timeGRB photometry, to enable rapid spectroscopic follow-up, we simulated the intrinsic GRB spectra to be bright,as might be expected during the first night of RATIRobservations. Each of the simulated RATIR magnitudesrequired a realistic error, representative of those typicallyreported by RATIR at these epochs, so we assigned anerror of ∆m AB = 0 . 03 to each of simulated magnitudes.In some cases this error may prove to be conservative,with RATIR capable of measuring magnitudes to accu-racies of ∆m AB ∼ . 02 within the first 12 hours afterthe initial Swift trigger (Littlejohns et al. 2013).With simulated RATIR photometry in hand, we thenused our fitting algorithm, which uses templates includ-ing all three extinction laws (and an A V = 0 model), tofind the most representative SED for the data. Table 4shows the full set of results obtained when fitting thesesimulated data.In all instances, our algorithm found the A V = 0 modelto be the worst fit to the simulated photometry. The χ value for each A V = 0 fit becomes larger with increased A V , in .For 5 of the 6 test cases, we found that the templatesolution with the correct extinction law yielded the fitwith the smallest χ . The margin by which this fit wasbetter than alternative extinction laws increased notablywhen increasing the amount of dust from A V = 0 . 25 to A V = 0 . 5, as expected.From the values shown in Table 4 it is clear that theMW type extinction law is the most easily identified byour algorithm, with the MW type template SED hav-ing the lowest χ and the recovered value of A V mostclosely matching A V , in in both MW tests. We demon-strate the SED templates from each of the extinction lawsimplemented by our algorithm in Figure 7. This morerobust identification is due to the larger prominence ofthe 2175 ˚A feature in comparison to the LMC and SMCextinction laws.The prominence of a strong 2175 ˚A bump in the MWextinction law has an additional implication for our fit-ting algorithm. At z sim = 3 . 36, the sharp drop in flux from neutral hydrogen absorption is not covered by theRATIR SED. In the absence of large quantities of dustin the host galaxy along the GRB sight line, our algo-rithm cannot determine the redshift of the GRB when z < 4. However, occurring at longer wavelengths thanneutral hydrogen absorption, a strong 2175 ˚A feature canprovide constraints on the galaxy host redshift.Table 4 demonstrates this, as the recovered value of z phot is in good agreement with z sim for a MW type in-put extinction law with A V, in = 0 . 5. Fitting the resultingsimulated photometry with a MW extinction law the 3 σ error bound states 3 . < z phot < . 73. Thus, in casesof strong MW type absorption, our algorithm can de-termine z phot to a slightly lower limit than indicated bythe results from § A V are best fitted with aMW type extinction law. Fitting dusty hosts We then considered the case where the intrinsic GRBemission is highly extinguished by dust in a host galaxyat a variety of simulated redshifts. With this in mind,we adopted A V, sim =0.5. Taking into consideration thecapabilities of the algorithm in the A V = 0 case, we onlyconsidered 4 < z sim < 10. As with the tests to identifydust extinction we used a bright intrinsic GRB spectrum,typical of observations within the first 12 hours after the Swift /BAT trigger. Accordingly, we assumed ∆m AB =0 . 03 for all bands in which we obtained detections.Figure 8 shows the results obtained when the same ex-tinction law was used in both the fitted template and pro-ducing the simulated photometry. In each case we plotall values where the template fitting algorithm could suc-cessfully recover a value of z phot . We have only includedplots for the MW and LMC type extinction laws, as dueto the high A λ for the SMC extinction law, the fitting al-gorithm was unsuccessful in producing a unique solutionfor most of the values of z sim with such a large quantityof SMC type dust in the host galaxy.Figure 8 show the recovered values of z phot remainsconsistent with z sim prior to z sim ∼ . z sim & . z phot that over predicts z sim . Thiseffect is largely due to the neutral hydrogen absorptionfeature falling between the RATIR Y - and J - bands, giv-ing large uncertainty in the precise wavelength at whichit occurs. Even with the larger uncertainty in the MWextinction law simulations, we can robustly determinethat z phot > . CONCLUSIONS We have presented a template fitting algorithm used todetermine photometric redshifts from RATIR data whena dropout candidate is present. This algorithm repre-sents the intrinsic GRB spectrum with a physically mo-tivated synchrotron model, includes dust extinction fromthe GRB host galaxy and absorption from interveningclouds of neutral hydrogen along the line of sight. Eachfitted SED therefore provides estimates of z phot , β opt and A V .2 Littlejohns et al. µ m)1000 F λ ( µ J y ) NoneMWLMCSMCr i Z Y J H Figure 7. Fitted SED templates to simulated RATIR photometry. The photometry was produced using a MW type extinction law with A V = 0 . 5. The black lines each describe a best fit template obtained by our fitting algorithm using a different extinction law, each beingdescribed in the key. The details for each fit are available in Table 4. z sim z pho t MW, A V =0.5 z sim z pho t LMC, A V =0.5 Figure 8. Fitted z phot as a function of the input redshift, z sim , used to create simulated RATIR magnitudes. The dotted line denotesequality. The left and right panels show results from using MW and LMC type extinction laws, respectively, to produce simulatedphotometry. The solid lines shows z phot found for each z sim , whilst the gray regions illustrates the 3 σ confidence region about each bestfit. dentifying high-redshift GRBs with RATIR 13 Table 4 Recovered values of A V , z phot and the corresponding χ values for testsconducted with moderate ( A V, sim = 0 . 25 and A V = 0 . 5) host dustextinction. In all cases z sim = 3 . 36. The type of extinction law and valueof A V used to produce the simulated RATIR photometry are shown in thefirst and second columns. Th extinction law used to fit z phot , the fittedvalues of A V and z phot and the associated χ value are also shown.Input Dust A V , in Fitted Dust A V , fit z phot β opt χ /ν MW 0.25 None 0.00 3.65 +0 . − . +0 . − . +0 . − . +3 . − . +0 . − . +0 . − . +0 . − . +3 . − . +0 . − . +3 . − . +1 . − . +3 . − . +0 . − . +0 . − . +0 . − . +3 . − . +0 . − . +2 . − . +3 . − . +3 . − . +0 . − . +0 . − . +2 . − . +0 . − . . < z phot < . z = 5 . β opt can help discern betweenthe multiple potential solutions, as demonstrated partic-ularly for the MW type extinction law.We also present the typical output of the algorithm,including SEDs and probability maps of the parameterspace for each extinction model. Analysis of the finelygridded probability maps, which focus on the most prob-ably region of z phot β opt parameter space, shows that foran LMC or SMC type extinction law the favored modelis of a host containing negligible quantities of dust. Thedegeneracy between several model solutions with a MWextinction law is partially lifted by using the Reichart(2001) prior probability distribution, which indicates alow value for A V best fits the data. With this priorprobability ruling out high A V solutions, the best fit MWsolution becomes one with negligible host dust content.With this solution, the three extinction laws all favor anear identical fit, which is consistent with a zero A V so-lution in both the obtained parameter values and χ fitstatistic.To ensure our algorithm is robust, we then create simu-lated RATIR photometry, typical of the quality expectedfrom the first night of observations after the Swift /BATtrigger. We firstly show that when there is no dustextinction in the host galaxy, we can successfully re-cover z phot using the template fitting, in the ranges4 . z sim . . . . z sim . . 5, for a varietyof β opt values. When 7 . . z sim . . 75 we remain ableto identify that z phot is within this range, and is thereforea target of high interest to larger facilities.Introducing a moderate quantity of dust extinctionto our simulations allows us to draw some conclusionsabout the dust present in the host. MW type dust is themost easily identifiable, due to the prominent feature at2175 ˚A. For weaker amounts of dust extinction it is dif-ficult to differentiate between LMC and SMC extinctionlaws. However, at A V ≈ . A V = 0 . z phot even when neutral hydrogen absorption occursslightly below the RATIR wavelength coverage. With alarge presence of MW type dust in a host galaxy our al-gorithm can determine z phot to an improved lower limitof z phot ∼ z sim . This isdue to the large amount of suppression of low wavelengthemission from the intrinsic GRB spectrum by the dustpopulation within the host galaxy. For both the MWand LMC extinction laws precise values of z phot wereobtained for z sim . . 5. After this point the accuracy of z phot reduces, although we remained able to robustlystate that the GRB was of high redshift. The associated3 σ errors were higher with the MW extinction law dueto the greater prominence of the 2175 ˚A feature in thistype of extinction law. This allows slightly lower redshiftsolutions to be found where the 2175 ˚A bump contributesto the suppression of lower wavelength emission.We also note that low mass stars in the Milky Wayprovide a population of interlopers that may occur neara Swift /XRT GRB position. To quantify the proba-bility that such a source is present in RATIR observa-tions we considered thick and thin disk components ofthe Milky Way with exponential vertical density pro-files (Juri´c et al. 2008). We determined the local den-sity of M-stars within a radius of 20 pc using the Hip-parcos Main Catalog (Perryman et al. 1997). This vol-ume was chosen to simultaneously maximize complete-ness and the total number of sources. We cross refer-enced these sources with the 2MASS All-Sky Catalog ofPoint Sources (Cutri et al. 2003), using the cross-matchservice provided by CDS, Strasbourg, to obtain their J band magnitude. We found the chance probability thatan M-star with J < 21 will be present in a 300 squarearcsecond region centered around the Swift /XRT posi-tion is less than 2.3 %. This corresponds to one (orfewer) chance M-star in observations of approximately43 different fields of view. Such events can be classifiedas non-GRB by their lack of temporal variability and theblackbody shape of the SED at low wavelengths.Since December 2012, RATIR has observed 56 GRBfields of view, with one instance of an M-star interloper.RATIR observations of the field of GRB 131127A founda red source near the Swift /XRT position, suggesting ahigh-redshift candidate ideal for follow-up (Butler et al.2013d). This source was observed at later epochs, reveal-ing no significant evidence for fading (Butler et al. 2013e;Im et al. 2013), leading to the conclusion that this sourcewas not a GRB.The presented methodology is aimed specifically atidentifying potential high-redshift GRBs, and providinga preliminary estimate of redshift. With z phot in hand,the justification of triggering larger spectroscopic facili-ties to measure a highly resolved spectroscopic redshiftis significantly higher. With the automated capabilitiesof RATIR it is therefore possible to request observationsfrom such large facilities at early epochs after the initial Swift trigger time, thereby obtaining high signal-to-noiseratio spectra required to answer fundamental questionsabout the high-redshift Universe.We thank Sandra Savaglio for a constructive refereereport on this work. We thank the RATIR projectteam and the staff of the Observatorio Astron´omico Na-cional on Sierra San Pedro M´artir. RATIR is a collab-oration between the University of California, the Uni-versidad Nacional Auton´oma de M´exico, NASA God-dard Space Flight Center, and Arizona State Uni-versity, benefiting from the loan of an H2RG detec-tor and hardware and software support from TeledyneScientific and Imaging. RATIR, the automation ofthe Harold L. Johnson Telescope of the Observato-rio Astron´omico Nacional on Sierra San Pedro M´artir,and the operation of both are funded through NASAdentifying high-redshift GRBs with RATIR 15grants NNX09AH71G, NNX09AT02G, NNX10AI27G,and NNX12AE66G, CONACyT grants INFR-2009-01-122785 and CB-2008-101958, UNAM PAPIIT grantIN113810, and UC MEXUS-CONACyT grant CN 09-283. Facilities: