J-PLUS: photometric calibration of large-area multi-filter surveys with stellar and white dwarf loci

July 31, 2019

J-PLUS: photometric calibration of large area multi-filter surveys

with stellar and white dwarf loci

C. López-Sanjuan

1

, J. Varela

1

, D. Cristóbal-Hornillos

1

, H. Vázquez Ramió

1

, J. M. Carrasco

2

, P. -E. Tremblay

3

,

D. D. Whitten

4

_{, V. M. Placco}

4

_{, A. Marín-Franch}

1

_{, A. J. Cenarro}

1

_{, A. Ederoclite}

5

_{, E. Alfaro}

6

_{, P. R. T. Coelho}

5

_,

F. M. Jiménez-Esteban

7, 8

, Y. Jiménez-Teja

9, 5

, J. Maíz Apellániz

7

, D. Sobral

10, 11

, J. M. Vílchez

6

, J. Alcaniz

9

,

R. E. Angulo

12, 13

, R. A. Dupke

9, 14, 15

, C. Hernández-Monteagudo

1

, C. L. Mendes de Oliveira

5

, M. Moles

1

, and

L. Sodré Jr.

5

1 _{Centro de Estudios de Física del Cosmos de Aragón, Unidad Asociada al CSIC, Plaza San Juan 1, 44001 Teruel, Spain}

e-mail: clsj@cefca.es

2 _{Institut de Ciències del Cosmos, Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, 08028 Barcelona, Spain} 3 _{Department of Physics, University of Warwick, Coventry, CV4 7AL, UK}

4 _{Department of Physics and JINA Center for the Evolution of the Elements, University of Notre Dame, Notre Dame, IN 46556,}

USA

5 _{Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de São Paulo, 05508-090 São Paulo, Brazil} 6 _{IAA-CSIC, Glorieta de la Astronomía s/n, 18008 Granada, Spain,}

7 _{Centro de Astrobiología (CSIC-INTA), ESAC Campus, Camino Bajo del Castillo s/n, 28692 Villanueva de la Cañada, Spain,} 8 _{Spanish Virtual Observatory, 28692 Villanueva de la Cañada, Spain}

9 _{Observatório Nacional, Rua General José Cristino, 77 - Bairro Imperial de São Cristóvão, 20921-400 Rio de Janeiro, Brazil} 10 _{Department of Physics, Lancaster University, Lancaster, LA1 4YB, UK}

11 _{Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands}

12 _{Donostia International Physics Centre (DIPC), Paseo Manuel de Lardizabal 4, 20018 Donostia-San Sebastián, Spain} 13 _{IKERBASQUE, Basque Foundation for Science, 48013, Bilbao, Spain}

14 _{University of Michigan, Department of Astronomy, 1085 South University Ave., Ann Arbor, MI 48109, USA} 15 _{University of Alabama, Department of Physics and Astronomy, Gallalee Hall, Tuscaloosa, AL 35401, USA}

Submitted July 2019

ABSTRACT

Aims.We present the photometric calibration of the twelve optical passbands observed by the Javalambre Photometric Local Universe Survey (J-PLUS).

Methods.The proposed calibration method has four steps: (i) definition of a high-quality set of calibration stars using Gaia informa-tion and available 3D dust maps; (ii) anchoring of the J-PLUS gri passbands to the Pan-STARRS photometric soluinforma-tion, accounting for the variation of the calibration with the position of the sources on the CCD; (iii) homogenization of the photometry in the other nine J-PLUS filters using the dust de-reddened instrumental stellar locus in (X − r) versus (g − i) colours, where X is the filter to calibrate. The zero point variation along the CCD in these filters was estimated with the distance to the stellar locus. Finally, (iv) the absolute colour calibration was obtained with the white dwarf locus. We performed a joint Bayesian modelling of eleven J-PLUS colour-colour diagrams using the theoretical white dwarf locus as reference. This provides the needed offsets to transform instrumental magnitudes to calibrated magnitudes outside the atmosphere.

Results. The uncertainty of the J-PLUS photometric calibration, estimated from duplicated objects observed in adjacent pointings and accounting for the absolute colour and flux calibration errors, are ∼ 19 mmag in u, J0378 and J0395, ∼ 11 mmag in J0410 and J0430, and ∼ 8 mmag in g, J0515, r, J0660, i, J0861, and z.

Conclusions. We present an optimized calibration method for the large area multi-filter J-PLUS project, reaching 1-2% accuracy within an area of 1 022 square degrees without the need for long observing calibration campaigns or constant atmospheric monitoring. The proposed method will be adapted for the photometric calibration of J-PAS, that will observe several thousand square degrees with 56 narrow optical filters.

Key words. methods: statistical – techniques: photometric

1. Introduction

The analysis of Milky Way (MW) stars and the understand-ing of extragalactic sources have greatly benefited from large (& 5 000 deg2) and systematic optical and near-infrared pho-tometric surveys, such as the second Palomar Observatory Sky Survey (POSS-II; Gal et al. 2004), the Sloan Digital Sky Survey (SDSS; Abazajian et al. 2009), the Two Micron All-Sky

Sur-vey (2MASS; Skrutskie et al. 2006), or the VISTA Hemisphere Survey (VHS; McMahon et al. 2013). These studies will move forward in the following decade with a bunch of on-going and planned next-generation surveys, some of them summarized in Table 1 for reference.

One fundamental step in the data processing of all the major surveys is the photometric calibration of the observations. The calibration process aims to translate the observed counts in

Table 1. Compilation of finished (F), on-going (O), and scheduled (S) optical and near-infrared large area (& 5 000 deg2_{) photometric surveys.}

Acronym Status Area Photometric system Reference

[deg2_]

POSS-II F 19 000 JFN Gal et al. (2004)

SDSS F 14 000 ugriz Abazajian et al. (2009)

Pan-STARRS O 31 000 grizy Chambers et al. (2016)

DES O 5 000 grizY Flaugher (2012)

Gaia O 41 253 G, GBP, GRP Gaia Collaboration et al. (2016)

DESI Legacy Imaging Surveys O 14 000 grz Dey et al. (2019)

SkyMapper O 20 000 uvgriz Wolf et al. (2018)

J-PLUS O 8 500 ugriz+ 7 medium bands Cenarro et al. (2019)

S-PLUS O 9 500 ugriz+ 7 medium bands Mendes de Oliveira et al. (2019)

LSST S 18 000 grizY Ivezic et al. (2008)

J-PAS S 8 500 56 bands (140Å) Benítez et al. (2014)

2MASS F 41 253 JHKs Skrutskie et al. (2006)

VHS F 19 000 JHKs McMahon et al. (2013)

UHS O 18 000 JKs Dye et al. (2018)

Euclid S 15 000 VIS+ YJH Laureijs et al. (2011)

tronomical images to a physical flux scale referred to the top of the atmosphere. Because accurate colours are needed to derive photometric redshifts for galaxies and atmospheric parameters for stars, and reliable absolute fluxes are involved in the estima-tion of the luminosity and the stellar mass of galaxies, current and future photometric surveys target a calibration uncertainty at 1% level and below to reach their ambitious scientific goals.

The traditional calibration approach relies in a network of standard stars with a well known flux across the wavelength range of interest. The monitoring of these standards with the survey photometric system permits to calibrate the observations. The calibration of large area multi-filter surveys has two main challenges that are not optimally tackled with this traditional method: (i) obtaining an homogeneous photometric calibration across areas of thousands of square degrees, and (ii) performing a consistent wavelength calibration for dozens of passbands.

Thanks to lessons learnt from SDSS, the repeated scan of calibration fields, and the constant monitoring of the sky condi-tions, methodologies such as ubercalibration, supercalibration, and hypercalibration; the estimation of photometric flat fields; or the forward photometric modelling have been successfully applied to reach 1% level precision in broad-band surveys (Pad-manabhan et al. 2008; Regnault et al. 2009; Wittman et al. 2012; Schlafly et al. 2012; Ofek et al. 2012; Burke et al. 2014, 2018; Scolnic et al. 2015; Magnier et al. 2016b; Finkbeiner et al. 2016; Zhou et al. 2018). These methodologies were envisioned to pro-vide an homogeneous calibration over large areas and can be also applied to multi-filter surveys, but their large number of pass-bands makes the calibration campaigns severely time consuming and the calibration observations can take as long as the scien-tific operations. To optimise the telescope time and speed up the survey progress, novel calibration strategies must be developed for projects such as the Javalambre Photometric Local Universe Survey (J-PLUS1; Cenarro et al. 2019), the Southern Photomet-ric Local Universe Survey (S-PLUS; Mendes de Oliveira et al. 2019), and the Javalambre Physics of the accelerating universe Astrophysical Survey (J-PAS2_{; Benítez et al. 2014).}

The present paper summarizes the efforts in the quest for an optimised photometric calibration procedure for J-PLUS. The survey started in November 2015 and in the last four years

sev-1 _j-plus.es 2 _j-pas.org

eral calibration methods have been implemented and tested. The growing amount of data, the improved knowledge of the tele-scope optics and the filter system, and the efforts of the commu-nity to produce other high-quality legacy datasets (Table 1) have permitted the fine tuning of the calibration method to achieve the 1% precision goal in most of the J-PLUS filters. As reference, we provide a brief description of the previous calibration pro-cedures applied to J-PLUS data in Sect. 3, and the instructions to update public J-PLUS photometry with the new calibration method presented along this work in Sect. 6.

This paper is organised as follows. In Sect. 2, we present the J-PLUS data and the ancillary datasets used in the calibra-tion process. A summary of the previous calibracalibra-tion methods is presented in Sect. 3. The current concordance photometric cal-ibration methodology is detailed in Sect. 4, and the calcal-ibration precision is presented in Sect. 5. The recipes to apply the new calibration to J-PLUS data are outlined in Sect. 6. We present our conclusions in Sect. 7. Magnitudes are given in the AB sys-tem (Oke & Gunn 1983).

2. J-PLUS photometric data

J-PLUS is a photometric survey of several thousand square de-grees that is being conducted from the Observatorio Astrofísico de Javalambre (OAJ, Teruel, Spain; Cenarro et al. 2014) using the 83 cm Javalambre Auxiliary Survey Telescope (JAST/T80) and the T80Cam, a panoramic camera of 9.2k × 9.2k pixels that provides a 2 deg2field of view (FoV) with a pixel scale of 0.5500pix−1 (Marín-Franch et al. 2015). The J-PLUS filter sys-tem, composed of twelve bands, is summarized in Table 2. The J-PLUS observational strategy, image reduction, and main sci-entific goals are presented in Cenarro et al. (2019).

Table 2. J-PLUS photometric system, extinction coefficients, and limiting magnitudes (5σ, 300

aperture) of J-PLUS DR1 (Cenarro et al. 2019).

Passband (X) Central Wavelength FWHM mDR1

lim kX=

AX

E(B−V) Comments

[nm] [nm] [AB]

u 348.5 50.8 20.8 4.916 In common with J-PAS

J0378 378.5 16.8 20.7 4.637 [OII]; in common with J-PAS

J0395 395.0 10.0 20.7 4.467 Ca H+K J0410 410.0 20.0 20.9 4.289 Hδ J0430 430.0 20.0 20.9 4.091 G band g 480.3 140.9 21.7 3.629 SDSS J0515 515.0 20.0 20.9 3.325 Mgb Triplet r 625.4 138.8 21.6 2.527 SDSS

J0660 660.0 13.8 20.9 2.317 Hα; in common with J-PAS

i 766.8 153.5 21.1 1.825 SDSS

J0861 861.0 40.0 20.2 1.470 Ca Triplet

z 911.4 140.9 20.3 1.363 SDSS

by optical artefacts were masked. The DR1 is publicly available at the J-PLUS website3_.

The new calibration process presented in Sect. 4 uses J-PLUS DR1 in combination with ancillary data from Gaia and the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS), so we describe these datasets in the following.

2.1. Pan-STARRS DR1

The Pan-STARRS 1 is a 1.8 m optical and near-infrared tele-scope located on Mount Haleakala, Hawaii. The teletele-scope is equipped with the Gigapixel Camera #1 (GPC1), consisting of an array of 60 CCD detectors, each 4 800 pixels on a side (Cham-bers et al. 2016).

The 3π Stereoradian Survey (referred as PS1 hereafter; Chambers et al. 2016) covers the sky north of declination δ = −30◦_{in four SDSS-like passbands, griz, with an additional} pass-band in the near-infrared, y. The entire filter set spans the range 400 − 1 000 nm (Tonry et al. 2012).

Astrometry and photometry were extracted by the Pan-STARRS 1 Image Processing Pipeline (Magnier et al. 2016a,b,c; Waters et al. 2016). PS1 photometry features a uniform flux cal-ibration, achieving better than 1% precision over the sky (Mag-nier et al. 2016b; Chambers et al. 2016). In single-epoch pho-tometry, PS1 reaches typical 5σ depths of 22.0, 21.8, 21.5, 20.9, and 19.7 in grizy, respectively (Chambers et al. 2016). The PS1 DR1 occurred in December 2016, and provided a static-sky cat-alogue, stacked images from the 3π Stereoradian Survey, and other data products (Flewelling et al. 2016).

Because of its large footprint, homogeneous depth, and ex-cellent internal calibration, PS1 photometry provides an ideal reference for the calibration of the gri J-PLUS broad-bands.

2.2. Gaia DR2

The Gaia spacecraft is mapping the 3D positions and kinemat-ics of a representative fraction of MW stars (Gaia Collaboration et al. 2016). The mission will eventually provide astrometry (po-sitions, proper motions, and parallaxes) and optical spectropho-tometry for over a billion stars, as well as radial velocity mea-surements of more than 100 million stars.

In the present paper, we used the Gaia DR2 (Gaia Collab-oration et al. 2018b). It contains five-parameter astrometric

de-3 _{www.j-plus.es/datareleases/data_release_dr1}

terminations and provides integrated photometry in three broad-bands G, GBP(330 − 680 nm), and GRP(630 − 1 050 nm) for 1.4 billion sources with G < 21. The typical uncertainties in Gaia DR2 measurements at G= 17 are ∼ 0.1 marcsec in parallax, ∼ 2 mmag in G−band photometry, and ∼ 10 mmag in GBPand GRP magnitudes (Gaia Collaboration et al. 2018b).

3. Previous calibration methods applied to J-PLUS data

The different procedures implemented to perform the photomet-ric calibration of the J-PLUS DR1 observations have provided precious knowledge to reach the optimised method presented in Sect. 4. Thus, a proper presentation of these methods is manda-tory to understand the strengths and weaknesses of each proce-dure, and motivate the need for a new methodology.

The ultimate goal of any calibration strategy is to obtain the zero point (ZP) of the observation, that relates the magnitude of the sources in passband X on top of the atmosphere with the magnitudes obtained from the analogue to digital unit (ADU) counts of the reduced images. We simplify the notation in the following using the passband name as the magnitude in such fil-ter. Thus,

X= −2.5 log₁₀(ADUX)+ ZPX. (1)

In the estimation of the J-PLUS DR1 raw catalogues, the reduced images were normalized to a one-second exposure and ZPX= 25 was used. This defined the instrumental magnitudes Xins. 3.1. Spectro-photometric standard stars

The main sources for the spectro-photometric standard stars (SSSs) are the spectral libraries CALSPEC4, the Next Genera-tion Spectral Library5_{, and STELIB (Le Borgne et al. 2003).} Fol-lowing the calibration procedure based on Bouguer fitting lines, each SSS is observed at different airmasses throughout the night to derive the atmospheric extinction coefficient and the photo-metric zero point of the system. The synthetic magnitudes of the SSSs were estimated by convolving the reference spectra with the J-PLUS photometric system, and the instrumental magni-tudes were estimated from the Moffat (1969) profile fitting to the observed light distribution of the SSSs. For this procedure to

be accurate, the atmospheric conditions must be stable along the night.

• Strengths: Consistent absolute flux calibration of the 12 J-PLUS filters.

• Weaknesses: The calibration observations consume signifi-cant fraction of telescope time. The typical magnitudes of the SSSs can produce saturated images in the broad-bands. Can be only applied in full photometric nights.

3.2. Comparison with broad-band photometry

The significant (∼ 80%) overlap between J-PLUS and SDSS footprints allows calibrating the J-PLUS broad-band observa-tions against the corresponding ones in SDSS. This technique was used to calibrate the ugriz bands by comparing J-PLUS 600 aperture instrumental magnitudes and SDSS PSF magnitudes. Because of differences in the effective transmission curves be-tween SDSS and J-PLUS photometric systems, colour-term cor-rections need to be applied to the SDSS magnitudes to obtain the corresponding J-PLUS photometry. These corrections are of particular importance in the case of the u band, where filters are known to be significantly different.

The same procedure was applied using the PS1 photometry as reference. In this case, any J-PLUS observation is covered by PS1, but only the griz bands are available. The colour-term corrections are significant in the case of the g band.

• Strengths: Reliable and accurate flux calibration of the J-PLUS broad-bands. High density of sources to perform the calibration. Can be also applied in non photometric nights. • Weaknesses: The calibration of the seven J-PLUS medium

bands is missing. Only a fraction of the J-PLUS area is cov-ered by SDSS, while we have no access to the u band with PS1. We inherit any flux calibration bias affecting the ref-erence photometry (see Lorenzo-Gutiérrez et al. 2019, for caveats about SDSS photometry at g. 15).

3.3. Comparison with SDSS spectroscopy

This method starts by convolving the SDSS stellar spectra with the spectral response for each J-PLUS passband, yielding syn-thetic magnitudes. Their comparison to the observed magnitudes in 600aperture provides estimates for the zero points. Although the sky coverage of the SDSS spectra is smaller and sparser than J-PLUS photometry, it can be used to calibrate those J-PLUS passbands that have no photometric counterpart. In particular, given the spectral coverage of the SDSS spectra, these are used to calibrate the J-PLUS passbands from J0395 to J0861, includ-ing gri broad-bands. With the installation of the BOSS trograph (Smee et al. 2013), the wavelength range of the spec-tra was extended to the blue, thus allowing the calibration of the J0378 band in areas of the sky for which BOSS spectra are available. However, u and z bands still fall out of the covered range by SDSS spectroscopy. Given the large FoV of T80Cam at JAST/T80, dozens of high-quality SDSS stellar spectra in a single J-PLUS pointing are frequent.

• Strengths: Consistent flux calibration of the seven medium-band J-PLUS filters. Can be applied in non photometric nights.

• Weaknesses: The calibration of u and z is missing. SDSS spectroscopy does not cover all the J-PLUS area. Source density is low with respect to the photometric case. We in-herit any flux calibration bias affecting SDSS spectra.

3.4. Stellar locus regression

The previous procedures were designed to be applied to any sin-gle exposure or any combination of exposures in a given filter, independently of the observations in any other band. However, by combining the information from different bands, it is possible to apply methods that enable anchoring the calibration through-out the spectral range. One particular approach is the use of the stellar locus (Covey et al. 2007; High et al. 2009; Kelly et al. 2014; Kuijken et al. 2019). This procedure takes advantage of the way stars with different stellar parameters populate colour-colour diagrams, defining a well-limited region (stellar locus) whose shape depends on the specific colours used.

The implemented stellar locus regression (SLR) method first constructs the median stellar locus in all the 2145 possible colour-colour combinations in J-PLUS. The initial photometry used to estimate the median locus relies on the previous cali-bration procedures: SDSS photometry for u and z, and SDSS spectroscopy for the rest of the J-PLUS passbands. The SLR works with relative colours, so a reference filter is needed. In this case, the i band provided the best results and was anchored to the available broad-band photometric reference: SDSS or PS1 in those pointings outside the SDSS footprint. Then, the distance of the 2145 stellar loci in each pointing to the median ones was minimized in an iterative process, leading to eleven offsets per pointing.

• Strengths: Consistent relative flux calibration of the 12 J-PLUS filters in all the surveyed area. Can be applied in non photometric nights.

• Weaknesses: Needs a minimum density of stars to populate the stellar locus. Can not be used for standalone calibration of one image. Needs at least one reference band with ex-ternal calibration. The current version does not include the effect of MW dust reddening in the estimation of the median stellar locus.

3.5. Summary

A detailed description of the previous methods can be found in Varela & Cristóbal-Hornillos (2017). The tests performed re-veal that the best current option is the SLR, because it provides a consistent calibration in all the J-PLUS filters, pointings, and atmospheric conditions. The SLR has been therefore the refer-ence calibration method in J-PLUS DR1, and all the available calibrations for a given filter and pointing are accessible in the J-PLUS database ADQL table jplus.CalibTileImage.

Fig. 1. Flowchart of the calibration method presented in this work. Ar-rows that originate in small dots indicate that the preceding data prod-uct is an input to the subsequent analysis. Datasets are shown with their project logo, and external codes or models with grey boxes. The rounded-shape boxes show the calibration steps. The asterisk marks those steps based on dust de-reddened magnitudes. White boxes show intermediate data products, and ovals highlight publicly available data products of the calibration process.

4. J-PLUS photometric calibration with stellar and white dwarf loci

In this section, we provide the details about the proposed methodology for the photometric calibration of the multi-filter J-PLUS project. We started by gathering the needed informa-tion to define a high-quality set of stars for calibrainforma-tion (Sect.4.1). Then, we anchored the J-PLUS photometry in the gri broad-bands to the PS1 photometry (Sect.4.2). Next, we homogenize the photometric solution along the J-PLUS area in the other nine passbands with the instrumental stellar locus (Sect. 4.3). Finally, we estimate the absolute colour calibration of the J-PLUS pass-bands with the white dwarf locus (Sect. 4.4). The performance and the error budget of the obtained photometric calibration are presented in Sect. 5. To guide the reader, a flowchart of the cali-bration process is presented in Fig. 1.

The strengths of the new method are that it permits a consis-tent flux calibration of the 12 J-PLUS filters in all the surveyed area, can be applied in non photometric nights, no previous cal-ibration of the medium bands is needed, and includes the effect of MW dust in the stellar locus estimation. The weaknesses are similar to the SLR ones, mainly the need of a minimum den-sity of stars to populate the stellar and white dwarf loci, and the need of one reference band with external absolute calibra-tion. The former issue is circumvented thanks to the large FoV of T80Cam at JAST/T80, that always provides a few hundred high-quality stars for calibration, and to the large area already covered by J-PLUS DR1, that provides enough numbers of the sparse white dwarfs to take advantage of their locus. The latter issue is mitigated thanks to the excellent external photometric reference provided by PS1 in the gri passbands.

The J-PLUS instrumental magnitudes used for calibration were measured in a 600diameter aperture. This aperture ensures a low flux contamination from neighbouring sources and is not dominated by background noise, but it is not large enough to capture the total flux of the stars. Thus, we applied an aperture correction Caperthat depends on the pointing and the passband. The aperture correction was computed from the growth curves of bright, non-saturated stars in the pointing. The typical num-ber of stars used is 50 and the median aperture correction varies from Caper = −0.09 mag in the u band to Caper = −0.11 mag in the z band, with a median value of Caper= −0.1 mag for all the filters. We assumed that the J-PLUS 600magnitudes corrected by aperture effects provided the total flux of stars.

4.1. Step 1: definition of the high-quality stellar set for calibration

The initial stage of our methodology aims to define a high-quality sample of stars to perform the photometric calibration. We started by cross-matching the J-PLUS sources with signal-to-noise S/N > 10 and SExtractor photometric flag equal to zero (i.e. with neither close detections nor image problems) in all the 12 passbands against the Gaia DR2 catalogue using a 1.500 radius6. We discarded those J-PLUS sources with more that one Gaiacounterpart, and those with either S/N < 3 in Gaia par-allax, noted $ [arcsec], or without a photometric measurement in any G, GBP, or GRPpassband. We obtained 496 798 unique high-quality stars for calibration.

We applied the correction to the G photometry and the Vega to AB conversions presented in Maíz Apellániz & Weiler (2018). The median G magnitude of the calibration sample is G= 15.7 mag, with 99% of the sources having G. 17.5 mag.

We worked with dust de-reddened magnitudes and colours in several stages of the calibration process. We computed the extinction coefficients kXof each J-PLUS passband X using the extinction law presented in Schlafly et al. (2016, S16 hereafter)7. These coefficients, presented in Table. 2, assume RV = 3.1. The de-reddened J-PLUS photometry, either instrumental or cal-ibrated, is noted with the subscript 0 and was obtained as

X₀= X − k_XE(B − V). (2)

We estimated the colour excess E(B − V) [mag] of each J-PLUS + Gaia matched source from the 3D dust maps provided by Bayestar178 (Green et al. 2018). As stated by the authors,

6 _{The full J-PLUS versus Gaia catalogue can be found in the ADQL}

table jplus.xmatch_gaia_dr2 at J-PLUS database

−1.0 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0

(G

BP

− G

RP

)

0 −5 0 5 10 15

M

G0

White Dwarf (WD) : 293 Main sequence (MS) : 465 583 Giant branch (GB) : 30 922

Fig. 2. Absolute magnitude in the G band versus GBP− GRPcolour diagram,

corrected by dust reddening, of the 496 798 high-quality sources in com-mon between J-PLUS DR1 and Gaia DR2. The black dots are individual measurements. The coloured solid con-tours show density of objects, starting on 25 mag−2 _{and increasing by a}

fac-tor of ten in each step. We define three populations on this diagram: main se-quence stars (465 583; white area), gi-ant branch stars (30 922; grey area), and white dwarfs (293; blue area).

the colour excess EB17 retrieved by Bayestar17 is not directly E(B − V), so we scaled the output from Bayestar17 to ensure the same (r − i) colour excess both in J-PLUS and PS1. This im-plies E(B − V)= 0.92 × EB17(see Green et al. 2018, for further details). We will study the impact of the assumed extinction law in our results in Sect. 5.5.

The parallax measured by Gaia can be used to estimate the distance to the calibration stars. As discussed in Luri et al. (2018) and Bailer-Jones et al. (2018), such estimation should account by the inherent asymmetry in the parallax to distance transformation. To properly account for the uncertainties in the 3D dust maps and in the distances estimated from Gaia paral-laxes, we extracted 10 000 random points $randfrom a Gaussian distribution $ ± σ$ in the parallax of each source. Then, we imposed a positive parallax value and computed the attenuation at the corresponding distance drand = 1/$randand sky position using each time a random dust map solution from Bayestar17. Then, the median and the ±34% of the attenuation distribution were recorded as the value of the colour excess E(B − V) and its error. We checked that the colour excess distribution is Gaussian in most cases, providing a proper description of E(B − V) for each calibration star. This procedure naturally accounts for the asymmetry in the distances and applies a d > 0 prior (i.e. no Galaxy model has been assumed in the computation of the dis-tances; see Luri et al. 2018 and Bailer-Jones et al. 2018 for an extensive discussion).

The extinction coefficients of G, GBP, and GRPwere obtained as for the J-PLUS passbands; kG = 2.600, kGBP = 3.410, and

kGRP = 1.807. We note that this provides first order de-reddened

magnitudes and colours, since the proper extinction correction of Gaiaphotometry is colour and dust-column dependent (Daniel-ski et al. 2018; Gaia Collaboration et al. 2018a). However, the low extinction at the J-PLUS pointings makes this simple correc-tion sufficient for our goal, i.e. to define a sample of calibration stars.

We estimated the G−band absolute magnitude of the J-PLUS + Gaia sources as

MG0= G − kGE(B − V)+ 5 log10($)+ 5. (3) This estimation assumes a dust de-reddening using the Bayestar17 colour excess with the simplified extinction

coef-ficients aforementioned, and the inverse of the parallax as a dis-tance proxy. We note that the latter is a crude approximation to the Bayesian distance provided by Bailer-Jones et al. (2018). Because we aim to define general populations to calibrate the J-PLUS photometry, all these simplifications fulfil our require-ments.

The absolute magnitude - colour diagram of the J-PLUS+ Gaiasample of high-quality stars is presented in Fig. 2. We se-lected three populations on this diagram, named main sequence (MS) stars, giant branch (GB) stars, and white dwarfs (WDs). Formally, WD = [ MG0 > 7 ] ∩ [ (GBP− GRP)0< 0.35 ] ∩ [ MG0 > 10.5 + 7 × (GBP− GRP)0], (4) GB = [ MG0 < 4.1 ] ∩ [ (GBP− GRP)0> 0.35 ] ∩ [ MG0 < −1.5 + 8 × (GBP− GRP)0], (5) and MS= ( GB ∪ WD )c, (6)

where the superindex c denotes the absolute complement set. These broad classes can contain other types of objects, such as hot sub-dwarfs or unresolved binaries in the case of the main sequence area. We note that main sequence and giant branch stars could be used together in the next calibration steps, but we preferred to split them to minimize secondary branches in those colour-colour diagrams that include J-PLUS filters sensitive to gravity (i.e. J0515).

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

r

[mmag]

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

r

−

P

r

[mmag]

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

r

[mmag]

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

r

−

P

r

[mmag]

Fig. 3. Residuals of the comparison between PS1 and J-PLUS photometry in the r band,∆r = ∆r − ∆ratm, as a function of the (X, Y) position of

the source on the CCD. Upper left panel: Stacked residual map of all the J-PLUS DR1 pointings. Upper right panel: Stacked residual map after applying the plane correction estimated pointing-by-pointing. Lower left panel: Residual map of the pointing pid= 00315. The gradient in the

residuals of the individual sources (coloured circles) is fitted with a plane (coloured squares). The direction of maximum variation is shown with the arrow. Lower right panel: Residual map after applying the plane correction.

4.2. Step 2: anchoring gri broad-bands with PS1 data The next step in our calibration process aims to anchor the J-PLUS photometry to the PS1 photometric solution in the shared gribroad-band filters. The PS1 photometry is currently the ref-erence of other broad-band photometric surveys such as SDSS (Finkbeiner et al. 2016), the Dark Energy Spectroscopic Instru-ment (DESI) Legacy Surveys (Dey et al. 2019), or the Hyper Suprime-Cam Subaru Strategic Program (HSC, Aihara et al. 2019). Moreover, PS1 observations cover all the sky visible from OAJ, providing a consistent reference for any J-PLUS observa-tion.

We cross-matched our MS calibration set with the PS1 DR1 catalogue using a 1.500radius9. As in the Gaia case, we discarded those sources with more than one counterpart in the PS1 cata-logue or without a photometric measurement in any PS1 pass-band. We used the PS1 PSF magnitudes as reference. As stated by Magnier et al. (2016c), the PSF magnitudes in PS1 were op-timised to minimize the difference with respect to aperture cor-rected magnitudes, and thus are a good proxy for the total flux of stars.

To account for the differences between the J-PLUS and PS1 photometric systems, we applied the following transformation equations, TX= XPS1− XJ−PLUS, where X is the passband under 9 _{The full J-PLUS versus PS1 catalogue can be found in the ADQL}

table jplus.xmatch_panstarrs_dr1 at J-PLUS database

study CPS1 = gPS1− iPS1, (7) Tg = 0.8 − 88.6 × CPS1+ 22.5 × C2PS1[mmag], (8) Tr = 4.9 − 3.2 × CPS1+ 8.2 × C2PS1[mmag], (9) Ti = −2.2+ 3.9 × CPS1+ 7.6 × C2PS1[mmag], (10) Tz = −13.0 + 24.4 × CPS1+ 6.2 × C2PS1[mmag]. (11) These equations were estimated in two steps. First, we obtained an initial transformation by convolving the Pickles (1998) stel-lar library with both PS1 and J-PLUS photometric systems. We applied these initial transformations to the full MS calibration set with PS1 counterpart and accounted by residual correlations with (g − i)PS1colour in the range 0.4 < (g − i)PS1 < 1.4. This is the validity range of the reported transformation equations, so we only kept sources in this colour interval when comparing J-PLUS and PS1 photometry. The median residuals with colour between both photometric systems are below 2 mmag, but we can not trace the presence of absolute systematic differences. This issue will be explored in more detail in Sect. 5.3.

In the following, we use the r band as example, but the methodology was the same for the other broad-bands. We es-timated the difference between the transformed PS1 PSF cali-brated magnitudes, rPS1, and the J-PLUS instrumental magni-tudes, rins, as

We fitted the∆r distribution in each pointing pidwith a Gaussian function of median µr and dispersion σr. The zero point offset accounting for the atmosphere transparency of the observations was estimated as

∆ratm(pid)= µr. (13)

One important issue regarding large FoV instruments is the possible variation of the zero point with the position of the sources on the CCD. This can be due to the differential varia-tion of the airmass across the observavaria-tion, the presence of scat-tered light in the focal plane, or the change of the effective filter curves with position (see Regnault et al. 2009; Ibata et al. 2014; Starkenburg et al. 2017, for further details). We explore the pres-ence of such position-dependent effect by studying the residual difference

∆r = ∆r − ∆ratm (14)

as a function of the (X, Y) pixel position of the sources on the CCD. For this exercise, we combined the information of all the J-PLUS pointings. We find a clear gradient across the CCD in the difference between PS1 and J-PLUS photometry (Fig. 3, up-per left panel). This position-dependent effect impacts the pho-tometry of the sources at 2% level. Moreover, this gradient is not universal and depends on the pointing. The origin of such gradi-ent is still unclear and is under investigation. From the practical point of view, we performed a fit of the ∆r residuals in each pointing to a plane,

Pr(X, Y)= A × X + B × Y + C, (15)

where A, B, and C are the parameters that define the plane. This provided a position-dependent zero point for each source in the pointing pid, estimated as

ZPr(pid, X, Y) = ∆ratm(pid)+ Pr(pid, X, Y) + 25. (16) We present an example of this procedure for the J-PLUS pointing pid = 00315 in the lower panels of Fig. 3. The global residual, after applying the pointing-by-pointing plane correction, reduces to 0.5% level (Fig. 3, upper right panel). The improvement in the photometric precision of the J-PLUS calibration thanks to the plane correction is demonstrated in Sect. 5.1, where the common sources from adjacent pointings are used to estimate the uncer-tainties in the calibration process.

At the end of this step, the calibration of the J-PLUS gri passbands is anchored to the PS1 photometric solution. We also calibrated the z band, and it will be used as a control check (Sect. 5.3) of the calibration procedure.

4.3. Step 3: homogenization with the instrumental stellar locus

In the previous section, we calibrated the J-PLUS gri broad-bands thanks to the PS1 photometry. However, we have no ac-cess to a high-quality photometric reference in the seven J-PLUS medium-bands. To perform the calibration of these passbands, and of the u and z broad-bands, we used the stellar and the white dwarf loci (Sect. 4.4).

The stellar locus technique assumes that the intrinsic distri-bution of stars defines a narrow region in colours space, and that such locus is independent of the position on the sky. Thus, we can calibrate a set of filters by matching the observed locus to a reference one (e.g. High et al. 2009; Kelly et al. 2014; Kui-jken et al. 2019). First, we tested the photometric calibration of

−1

0

1

2

3

(g

_{− i)}

0

−1.0

−0.5

0.0

0.5

1.0

1.5

2.0

(g

−

r)

0 Stellar locus

−1

0

1

2

3

(g

_{− i)}

0

−2.0

−1.5

−1.0

−0.5

0.0

0.5

1.0

(i

−

r)

0

Fig. 4. Dust de-reddened colour-colour diagrams of the J-PLUS pass-bands anchored to the PS1 photometric solution. Dots are individual MS calibration stars. Upper panel: (g − r)0versus (g − i)0stellar locus.

Lower panel: (i − r)0versus (g − i)0stellar locus. The red line in both

panels shows the median stellar locus in the range −0.5 < (g − i)0< 2.4.

the gri bands performed in the previous section by estimating the MS stellar locus in the (X−r)0vs. (g−i)0colour-colour diagram, where X= {g, i}. We present these diagrams in Fig. 4. We found a clearly defined stellar locus, that is parametrized with a linear interpolation from the median of the (X − r)0colour distribution at different (g − i)0values in the range −0.5 < (g − i)0< 2.4. The dispersion of MS stars with respect to the parametrized stellar locus is 13 mmag for the g band and 12 mmag for the i band. This exercise demonstrates that the anchoring to the PS1 pho-tometry provides a well calibrated gri J-PLUS magnitudes, and thus can be used to study the stellar locus in the other J-PLUS passbands.

For each of the J-PLUS filters X that we aim to calibrate, we have to construct the colour-colour diagram (X − r)0vs. (g − i)0 and define the stellar locus. Because the gri filters were already calibrated and the impact of the MW interstellar extinction had been removed, any pointing-by-pointing discrepancy with re-spect to the stellar locus can be attributed to the effect of the atmosphere in X at the moment of the observation. For a given (g−i)0colour, the stellar locus defines an intrinsic (X−r)0colour that satisfies the following equation in each J-PLUS pointing pid, (X − r)SL₀ = h [Xins,i− 25+ ZPX(pid, Xi, Yi) − kXE(B − V)i] −

−1

0

1

2

3

(g

_{− i)}

0

3.2

3.4

3.6

3.8

4.0

4.2

4.4

(J

0660

ins

−

r)

0

Initial stellar locus J0660ins

−1

0

1

2

3

(g

_{− i)}

0

3.2

3.4

3.6

3.8

4.0

4.2

4.4

(J

0660

ins

−

r)

0

Initial stellar locus J0660ins+ ∆J0660atm

−1

0

1

2

3

(g

_{− i)}

0

3.2

3.4

3.6

3.8

4.0

4.2

4.4

(J

0660

ins

−

r)

0

Final stellar locus J0660ins+ ∆J0660atm

Fig. 5. Dust de-reddened (J0660ins− r)0versus (g − i)0colour-colour

diagram. Dots are individual MS calibration stars. Upper panel: Initial instrumental J0660 photometry. The solid line shows the linear fit to the data in the range 0.20 < (g − i)0 < 1.25. Middle panel: J0660

photometry corrected with the offsets ∆J0660atm. The red line shows

the linear fit in the upper panel. Bottom panel: Final instrumental stellar locus (red line) estimated as the median of the colour distribution in the range −0.5 < (g − i)0 < 2.4. The final ∆J0660atmwere estimated with

respect to this locus.

In the case of the J-PLUS medium bands, we have no access to the intrinsic stellar locus (X − r)SL

0 . We circumvent this issue including an extra term in the zero point,

ZPX(pid, X, Y) = ∆Xatm(pid)+ PX(pid, X, Y) + ∆XWD+ 25, (18) where∆XWD is a new offset that provides the absolute calibra-tion of the passband outside the atmosphere. In this seccalibra-tion, we detail the estimation of∆Xatmand PX, and we deal with∆XWDin Sect. 4.4. We note that the calibration of gri against PS1 implies ∆griWD∼ 0.

We use the filter J0660 as example in the following, but the procedure was the same in the other J-PLUS filters. We de-fined the initial version of the instrumental stellar locus (ISL), noted (X − r)ISL₀ , with a linear fit to the dust-corrected colour-colour data of those MS calibration stars in J-PLUS with 0.20 < (g − i)0< 1.25. In this process, the magnitudes outside the atmo-sphere in the gri bands and the instrumental magnitudes in the J0660 band were used (upper panel in Fig. 5). Formally, (X − r)ISL₀ = h [Xins,i− kXE(B − V)i] − [ri− krE(B − V)i] i, (19)

where the index i runs over all J-PLUS MS calibration stars at a given (g − i)0colour. We estimated the offsets ∆J0660atmas the median difference between the MS calibration stars with 0.20 < (g − i)0 < 1.25 in a given pointing and the initial ISL. Thanks to these initial offsets, the pointing-by-pointing differences are largely suppressed (middle panel in Fig. 5). Then, we estimated the final ISL with a linear interpolation from the median of the (J0660ins−r)0+∆J0660atmcolour distribution at different (g−i)0 in the range −0.5 < (g − i)0< 2.4, and computed the final offsets in each pointing as the difference with respect to this final locus (bottom panel in Fig. 5). We checked that extra iterations do not improve the results. After this process, the dispersion of all MS calibration stars with respect to the J0660 instrumental stellar locus had decreased from 57 mmag to 12 mmag.

The next step is to estimate the plane correction outlined in Sect. 4.2. Because we did not have access to external photome-try, we used the colour distance to the final ISL as reference to define the residual

∆J0660 = (J0660ins− r)0+ ∆J0660atm− (J0660 − r)ISL0 . (20)

The stacked residuals along the CCD position with respect to the final ISL are presented in the upper left panel of Fig. 6. As in the r-band case, a clear gradient emerges, but with a different direc-tion. Because the stellar locus has its own physical dispersion, we enhanced the signal in each individual pointing by splitting the CCD on 16 regions (4 × 4 grid) and computing the median ∆J0660 in each of these regions (Fig. 6, lower panels). In this process, we assumed that any measured trend is due to varia-tions in the J0660 photometry alone. Then, we fitted a plane to these median differences to obtain PJ0660(pid, X, Y). The stacked residuals after applying the plane correction are at 0.5% level (upper right panel in Fig. 6). The inclusion of the plane correc-tion further decreases the dispersion with respect to the ISL to 9 mmag. As in the broad-band case, the improvement in the pho-tometry thanks to the plane correction is evaluated in Sect. 5.1.

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

J

0660

[mmag]

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

J

0660

−

P

J 0660

[mmag]

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

J

0660

[mmag]

0

2000 4000 6000 8000

X [pixels]

0

2000

4000

6000

8000

Y

[pixels]

−15

−10

−5

0

5

10

15

∆

J

0660

−

P

J 0660

[mmag]

Fig. 6. Residuals of the comparison between the final ISL and the (J0660ins− r)0+ ∆J0660atm colour as a function of the (X, Y) position of

the source on the CCD. Upper left panel: Stacked residual map of all the J-PLUS DR1 pointings. Upper right panel: Stacked residual map after applying the plane correction estimated pointing-by-pointing. Lower left panel: Residual map of the pointing pid= 00315. The median residuals

with respect to the instrumental stellar locus in 16 regions (4 × 4 grid, dotted lines) covering the CCD (coloured circles) are fitted with a plane (coloured squares). The direction of maximum variation is shown with the arrow. Lower right panel: Median differences after applying the plane correction.

4.4. Step 4: absolute colour calibration with the white dwarf locus

The properties of white dwarfs make them excellent standard sources for calibration (Holberg & Bergeron 2006). The model atmospheres of WDs can be specified at ∼ 1% flux level with an effective temperature (Te_ff) and a surface gravity (log g). These parameters can be accurately estimated by spectroscopic analy-sis of the Balmer line profiles, providing a reference flux model for calibration. They are also mostly photometrically stable and statistically present lower levels of interstellar reddening than main sequence stars (WDs are intrinsically faint, so we only de-tect the nearby ones). Because of these properties, a significant observational and theoretical effort is still on-going to provide the best possible WD network to ensure a high-quality calibra-tion of deep photometric surveys (e.g. Bohlin 2000; Holberg & Bergeron 2006; Narayan et al. 2016, 2019, and references therein).

A set of well characterised WDs can be used to obtain global offsets in the calibration of photometric systems (Holberg & Bergeron 2006). This procedure implies two steps: first, we have to obtain the properties of the calibration WDs (i.e. ef-fective temperature and gravity) to derive their theoretical fluxes. Second, the observed fluxes are compared against those obtained from the convolution of the WD modelled spectra with the tar-geted photometric system. The difference between both

mea-surements provides offsets that corrects the initial calibration of the studied passbands.

The application of the scheme above to multi-filter surveys is severely time consuming, implying repeated observations of the sparse population of reference WDs. As an example, only one WD in the calibration network from Narayan et al. (2019) has been observed in J-PLUS DR1. Instead of the one-to-one comparison, we statistically analysed the distribution of WDs in eleven J-PLUS colour-colour diagrams (Figs. 7, 8, and 9) to obtain the offsets ∆XWD. These offsets translate the ISL outside the atmosphere (Eq. 18) and complete the calibration process.

The observational WD locus is well described by the the-ory and presents two branches, corresponding to hydrogen (DA) and helium (DB+ DC) white dwarfs (e.g. Holberg & Bergeron 2006; Ivezi´c et al. 2007; Ibata et al. 2017; Gentile Fusillo et al. 2019; Bergeron et al. 2019). Such populations are evident for X = {u, J0378, J0395, J0660}, where the hydrogen lines are more prominent. We aim to match the WD locus estimated with J-PLUS instrumental magnitudes to the expected from theory, obtaining the absolute colour calibration of the J-PLUS pass-bands. We used the r band as reference in this analysis, and thus ∆rWD = 0 by construction.

now an homogeneous instrumental photometry, so we can use the whole WD population present in the J-PLUS DR1 area to transport the final ISL outside the atmosphere.

The statistical modelling of the WD locus is described in Sect. 4.4.1. The WDs selected using the Gaia absolute mag-nitude - colour diagram (Sect. 4.1 and Fig. 2) were first cleaned from outliers (Sect. 4.4.2). Then, a joint analysis of the eleven possible (X − r)0 vs. (g − i)0 colour-colour diagrams was per-formed to obtain the∆XWD values (Sect. 4.4.3). We present the results of this analysis in Sect. 4.4.4. The uncertainties related with the absolute colour calibration performed in this section are discussed in Sect. 5.3.

4.4.1. Modelling the WD colour-colour diagrams

We have developed several tools to use PRObability Functions for Unbiased Statistical Estimations (PROFUSE10) of galaxy distributions across cosmic time. This includes galaxy luminos-ity functions (López-Sanjuan et al. 2017; Viironen et al. 2018), galaxy merger fractions (López-Sanjuan et al. 2015), mass-to-light ratio vs. colour relations (López-Sanjuan et al. 2019a), Hα emission-line fluxes (Vilella-Rojo et al. 2015; Logroño-García et al. 2019), stellar populations (Díaz-García et al. 2015), or star/galaxy classification (López-Sanjuan et al. 2019b). In this case, we applied our previous knowledge to perform a Bayesian modelling of the white dwarf locus.

The intrinsic distribution of interest is noted D, and provides the real values of our measurements for a set of parameters θ,

D(Creal₁ , Creal₂ |θ), (21)

where Creal 1 and C

real

2 are the real values of the colours unaffected by both observational errors and systematic offsets. In our case, Creal

1 = (g − i)0 and Creal2 = (X − r)0. We derived the posterior of the parameters θ that define the intrinsic distribution D with a Bayesian model. Formally,

P(θ | Cobs₁ , Cobs₂ , σC1, σC2) ∝ L (C obs 1 , C obs 2 |θ, σC1, σC2) P(θ), (22) where σC1 and σC2 are the uncertainties in the observed (g −

i)0 and (X − r)0colours, respectively, L is the likelihood of the data given θ, and P(θ) the prior in the parameters. The posterior probability is normalised to one.

The likelihood function associated with our problem is L (Cobs₁ , Cobs₂ |θ, σC1, σC2)= Y k Pk(Cobs1,k, C obs 2,k|θ, σC1,k, σC2,k), (23) where the index k spans the WDs in the sample, and Pktraces the probability of the measurement k for a set of parameters θ. This probability can be expressed as

Pk(Cobs1,k, C obs

2,k|θ, σC1,k, σC2,k)=

Z

D(Creal₁ , Creal₂ |θ) PG(Cobs1 | C real 1 , σC1,k) × PG(Cobs2 | C real 2 , σC2,k) dC real 1 dC real 2 , (24)

where the real values Creal₁ and Creal₂ derived from the model D are affected by Gaussian observational errors,

PG(x | x0, σx)= 1 √ 2πσexp  −(x − x0) 2 2σ2  , (25) 10 _{profuse.cefca.es}

−1.0 −0.5

0.0

0.5

(g

_{− i)}

0

2.5

3.0

3.5

4.0

4.5

5.0

(u

ins

−

r)

0

0

1

2

Probabilit

y

0

1

Probability

WDs Outlier WDs Total model DA DB + DC ∆uWD =−3.864 σint= 0.031

Fig. 7. Dust de-reddened (uins− r)0vs. (g − i)0colour - colour diagram

of the 293 high-quality white dwarfs in J-PLUS DR1 (clean sample, cyan dots; outliers, red dots). The solid lines show the theoretical locus for DA (orange) and DB+DC WDs (magenta). The grey scale shows the most probable model that describes the observations. The upper and right blue histograms show the (g − i)0 and (uins− r)0 projections of

the data, respectively. The projections of the total, DA, and DB+DC models are represented with the black, orange, and magenta lines. The model in all the J-PLUS colour-colour diagrams shares the parameters µ = −0.808, s = 0.400, α = 2.65, fDA = 0.85, log g = 8.01, and

∆C1 = 0.007 (see text for details). The values of the filter-dependent

parameters σintand∆XWDare labelled in the panel.

providing the likelihood of an observed colour given its real value and uncertainty. We have no access to the real values of the colours, so we marginalised over them in Eq. (24) and the likelihood is expressed therefore with known quantities.

We explore the parameters posterior distribution with the emcee code (Foreman-Mackey et al. 2013), a Python imple-mentation of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) proposed by Goodman & Weare (2010). The emcee code provides a collection of solutions in the parameter space, denoted θMC, with the density of solutions being proportional to the posterior probability of the parameters. We obtained the central values of the parameters and their un-certainties from a Gaussian fit to the θMCdistribution.

We define in the following the intrinsic distribution assumed for the WD locus, and the prior imposed to their parameters. The WD population was described as

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

3.5

4.0

4.5

5.0

5.5

(J

0378

ins

−

r)

0

0

1

Probability

∆J0378WD =−4.480 σint = 0.029

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

3.5

4.0

4.5

5.0

5.5

(J

0395

ins

−

r)

0

0

1

2

Probability

∆J0395WD =−4.597 σint = 0.029

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

3.0

3.5

4.0

4.5

(J

0410

ins

−

r)

0

0

1

Probability

∆J0410WD =−3.663 σint = 0.014

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

2.8

3.0

3.2

3.4

3.6

3.8

4.0

4.2

(J

0430

ins

−

r)

0

0

1

Probability

∆J0430WD =−3.602 σint = 0.013

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

−0.6

−0.4

−0.2

0.0

0.2

0.4

(g

−

r)

0

0

1

2

3

Probability

∆gWD =−0.003 σint = 0.005

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

3.0

3.2

3.4

3.6

3.8

(J

0515

ins

−

r)

0

0 1 2 3 4 5 6

Probability

∆J0515WD =−3.443 σint = 0.015

Fig. 8. Similar to Fig. 7, but for X= J0378, J0395, J0410, J0430, g, and J0515 passbands. We omit the (g − i)0projection because it is shared

by all the panels.

sample; and σintis the intrinsic dispersion (i.e. related to physi-cal properties) of the WD locus.

The theoretical loci for DA and DB+DC WDs were ob-tained from the 3D model atmospheres presented in Tremblay et al. (2013) and Cukanovaite et al. (2018), respectively. The high-resolution spectral models at different gravities (log g = 7, 7.5, 8, 8.5, and 9) were convolved with the J-PLUS filter sys-tem to obtain the theoretical WD locus. We performed a linear interpolation in the provided colours to access other gravity val-ues during the modelling. The colour variations due to the vari-ety of gravities in the WD population under study are absorbed by the σintparameter.

We included at this stage two systematic offsets in the mod-elling, ∆C1 and∆C2. These offsets affect the theoretical WD locus, displacing it to match the observed distribution. The gri broad-bands are anchored to the PS1 photometric solution,

im-plying that∆C1∼ 0 and∆C2= −∆XWD. We have assumed the r band as absolute reference for the J-PLUS colours, so the offset in C2is the needed one to transform X instrumental magnitudes into calibrated magnitudes outside the atmosphere.

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

3.7

3.8

3.9

4.0

4.1

4.2

4.3

4.4

4.5

(J

0660

ins

−

r)

0

0 1 2 3 4 5 6

Probability

∆J0660WD=−3.900 σint= 0.014

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

−0.2

0.0

0.2

0.4

0.6

(i

−

r)

0

0 1 2 3 4 5

Probability

∆iWD= 0.004 σint= 0.006

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

3.0

3.2

3.4

3.6

3.8

4.0

4.2

4.4

(J

0861

ins

−

r)

0

0

1

2

3

Probability

∆J0861WD=−3.369 σint= 0.015

−1.0 −0.5 0.0

0.5

(g

_{− i)}

0

1.8

2.0

2.2

2.4

2.6

2.8

3.0

3.2

(z

ins

−

r)

0

0

1

2

3

Probability

∆zWD =−2.229 σint= 0.012

Fig. 9. Similar to Fig. 7, but for X= J0660, i, J0861, and z passbands. We omit the (g − i)0projection because it is shared by all the panels.

4.4.2. Removal of outlier WDs

One of the main advantages of WDs as calibration sources is their well-known physics. However, our initial WD sample de-rived from Gaia data in Sect. 4.1 can be contaminated by unre-solved WD+M binaries, foreground galaxies and neighbouring stars in the 600aperture, calcium and magnetic white dwarfs that are not reproduced by our assumed theoretical tracks, variable ZZ Ceti stars, etc. All these kind of sources could have colours far from the model expectations, biasing our analysis. Thus, our first goal is to clean up the initial WD sample from these physical outliers.

One way to define the clean WD sample is to visually in-spect the stamps, the J-PLUS photo-in-spectra, and the ancillary data of the 293 initial sources to select and remove the outliers. This process is subjective and time-consuming, so we decided to apply a statistical and automatic procedure to define the outlier WDs. We included a new component in the distribution of the J-PLUS WDs to account for the presence of outlier sources. This component was defined as a uniform density and is regulated with a new parameter called fout. Formally,

D(Creal₁ , Creal₂ |θ_WD, f_out)= (1 − fout) × DWD(Creal1 , C real

2 |θWD)

+ fout× U, (27)

where the function U provides a uniform probability density in colour-colour space. To minimize the degeneracies be-tween parameters in those colour-colour diagrams where DA and DB+DC white dwarfs are not well separated, we fixed µ = −0.8, s = 0.4, α = 2.8, fDA = 0.85, log g = 8, and ∆C1= 0. Thus, we only had three free parameters in this analy-sis, θWD = {σint, ∆C2} and fout. We obtained the most probable values for these parameters and computed the probability of each

white dwarf to be part of the desired WD locus. We only retained those sources with a probability larger than 97.5%, and the rest of the WDs were marked as outliers.

The selection of the outlier WDs was done in sequence, start-ing in the z filter and movstart-ing to the bluer passbands. We started in the reddest band because WD+M binaries, which dominate the outlier WDs, are easily detected in this colour-colour dia-gram. Those WDs marked as outliers in one band were not used in the subsequent analysis. We found 28 outlier WDs with this procedure, 10% of the initial sample. We repeated the full pro-cess, starting again with the z band, and no additional WD was marked as outlier. The remaining 265 WDs were used in the joint study presented in the next section.

4.4.3. Joint modelling of the WD locus

Table 3. Estimated offsets to transport the instrumental stellar locus out-side the atmosphere, intrinsic dispersion of the WD locus, and the final median zero points in J-PLUS DR1. The r band was used as reference in the estimation of the colour offsets.

Passband (X) ∆XWD σint hZPXi

[mag] [mag] [mag]

u −3.864 ± 0.005 0.031 ± 0.005 21.15 J0378 −4.480 ± 0.005 0.029 ± 0.005 20.53 J0395 −4.597 ± 0.005 0.029 ± 0.005 20.40 J0410 −3.663 ± 0.004 0.014 ± 0.005 21.34 J0430 −3.602 ± 0.004 0.013 ± 0.004 21.39 g −0.003 ± 0.002 0.005 ± 0.003 23.60 J0515 −3.443 ± 0.003 0.015 ± 0.003 21.57 r · · · 23.66 J0660 −3.900 ± 0.003 0.014 ± 0.004 21.12 i 0.004 ± 0.002 0.006 ± 0.003 23.35 J0861 −3.369 ± 0.004 0.015 ± 0.006 21.65 z −2.229 ± 0.004 0.012 ± 0.005 22.79

4.4.4. Absolute colour calibration from WD locus modelling We present in this section the final results of the WD locus mod-elling and the estimation of the offsets ∆XWDneeded to transport the J-PLUS instrumental magnitudes to calibrated magnitudes outside the atmosphere. The analysed data and the final mod-elling of the WD locus are presented in Figs. 7, 8, and 9. We summarize the obtained∆XWDand σintin Table 3.

We find a good agreement between the J-PLUS photomet-ric data and the WD locus model. The (g − i)0 distribution is parametrised with µ = −0.808 ± 0.007, s = 0.400 ± 0.007, and α = 2.65 ± 0.13. This distribution has a clear tail towards red colours (Fig. 7), that is properly described thanks to the skew-ness parameter α.

We obtain a DA fraction of fDA= 0.85 ± 0.01 and a median gravity of the WD population of log g = 8.01±0.03. This grav-ity value is similar to previous analysis (see Jiménez-Esteban et al. 2018; Gentile Fusillo et al. 2019; Tremblay et al. 2019; Bergeron et al. 2019, and references therein). We were able to constraint these two parameters using instrumental magnitudes in most of the J-PLUS filters thanks to the presence of the DA and DB+DC branches, and to the variation of the locus curvature with the gravity.

The offsets ∆XWD permit to obtain the final zero points of the 511 J-PLUS DR1 pointings. The typical error in these o ff-sets is ∼ 5 mmag, and thus a few hundred WDs are enough to provide robust results. We report the median J-PLUS zero points in Table 3. We found that the offset in the (g − i)0colour is not zero, with∆C1 = 7 ± 2 mmag. Thanks to the joint WD locus modelling, we are able to find residual differences in the cali-bration of the J-PLUS g and i passbands with respect to the PS1 photometric system (Sect. 4.2). We further discuss this issue in Sect. 5.3.

Finally, we comment on the values obtained for the intrinsic dispersion of the WD locus. We find that the bluer bands have a dispersion of ∼ 0.03 mag, larger than the ∼ 0.015 mag ex-hibited by the other filters. This highlights the larger impact of gravity variations in the photometry of the bluer J-PLUS filters and their importance for the study of individual WD properties. Such analysis is beyond the scope of the present paper, and will be addressed in future works of the J-PLUS collaboration.

5. Calibration performance and error budget

The methodology presented in the previous section aims to pro-vide the photometric calibration of the multi-filter J-PLUS ob-servations. In this section, we test the performance of the calibra-tion process by studying the photometric differences of sources observed by two adjacent pointings (Sects. 5.1 and 5.2). We also discuss the absolute colour (Sect. 5.3) and flux (Sect. 5.4) un-certainties in our calibration. The impact of the assumed MW extinction is explored in Sect. 5.5. We compare the new J-PLUS calibration with the previous ones in Sect. 5.6. Finally, the calibrated stellar locus is compared against stellar libraries in Sect. 5.7. We summarize the error budget of the calibration process in Table. 4.

5.1. Internal precision from overlapping areas

We measured the relative uncertainty (i.e. the precision) in the calibration by comparing the photometry of those MS stars ob-served independently in the overlapping areas between adjacent pointings. We computed the differences in the calibrated magni-tudes and estimated the median of those sources shared by every pair of overlapping pointings. We have 1173 unique pair point-ings in J-PLUS DR1. Then, the distribution of these median differences was used to estimate the relative uncertainty in the calibration. The distributions are well described by Gaussian functions and the desired precision is obtained as σ/

√

2, where σ is the measured dispersion. We used the pointing-by-pointing median instead of the total distribution for individual sources be-cause (i) the calibration was performed pointing-by-pointing, so this is the natural reference unit; (ii) we minimize the larger sta-tistical weight of the densest pointings; and (iii) we minimize the broadening of the distribution due to the uncertainties in the magnitude measurements.

We summarize our finding in Table 4 and Fig. 10. The rel-ative uncertainty is ∼ 18 mmag in u, J0378, and J0395; ∼ 9 mmag in J0410 and J0430; and ∼ 5 mmag in the other filters. In Table 4 and Fig. 10, we also present the relative uncertain-ties derived with the stellar locus regression method and with our methodology when the plane correction is neglected. We found that the SLR calibration is clearly improved by the new procedure even without the plane correction at filters bluer than J0515. This is due to the inclusion of the MW extinction in our methodology, that is more prominent in the bluer bands. A great improvement in the redder bands (factor of 2-3) is feasible as a consequence of the plane correction, where this improvement is mild (∼30%) in the three bluer bands. This is due to the intrinsic properties of the stellar locus in these passbands, that is broader because of metallicity differences in the stars.

We conclude that the photometric precision of J-PLUS DR1 has been improved by a factor of two with respect to previous calibration processes without the need of time consuming cali-bration observations or constant atmospheric monitoring.

5.2. Photometric precision from giant branch stars

−150 −100 −50 0 50 100 150 hu1− u2i [mmag] 0.0 0.5 1.0 1.5 2.0 100 × Probabilit y σISL+P= 17 mmag −150 −100 −50 0 50 100 150 hJ03781− J03782i [mmag] 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 100 × Probabilit y σISL+P= 19 mmag −150 −100 −50 0 50 100 150 hJ03951− J03952i [mmag] 0.5 1.0 1.5 100 × Probabilit y σISL+P= 17 mmag −100 −50 0 50 100 hJ04101− J04102i [mmag] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 100 × Probabilit y σISL+P= 9 mmag −100 −50 0 50 100 hJ04301− J04302i [mmag] 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 100 × Probabilit y σISL+P= 9 mmag −50 0 50 hg1− g2i [mmag] 0 1 2 3 4 5 6 7 100 × Probabilit y σISL+P= 4 mmag −60 −40 −20 0 20 40 60 hJ05151− J05152i [mmag] 0 1 2 3 4 5 100 × Probabilit y σISL+P= 6 mmag −40 −20 0 20 40 hr1− r2i [mmag] 0 1 2 3 4 5 6 7 8 100 × Probabilit y σISL+P= 4 mmag −60 −40 −20 0 20 40 60 hJ06601− J06602i [mmag] 0 1 2 3 4 5 6 7 100 × Probabilit y σISL+P= 4 mmag −40 −20 0 20 40 hi1− i2i [mmag] 0 1 2 3 4 5 6 7 100 × Probabilit y σISL+P= 4 mmag −60 −40 −20 0 20 40 60 hJ08611− J08612i [mmag] 0 1 2 3 4 5 6 100 × Probabilit y σISL+P= 5 mmag −60 −40 −20 0 20 40 60 hz1− z2i [mmag] 0 1 2 3 4 5 6 100 × Probabilit y σISL+P= 5 mmag

Fig. 10. Distribution of median differences in the photometry of MS stars independently observed by two adjacent pointings. In all the panels the black histogram shows the results from the SLR (reference photometry in J-PLUS DR1), the grey filled histogram shows the results after applying ∆Xatm, and the coloured histogram after applying∆Xatmand PX(X, Y). The solid line is the best Gaussian fit to the latest case. The uncertainty in

the calibration is labelled in the panels and was estimated as the dispersion of the fitted Gaussian divided by the square root of two. We present, from top to bottom and from left to right, the filters u, J0378, J0395, J0410, J0430, g, J0515, r, J0660, i, J0861, and z.

With the above caveats in mind, the results summarized in Table 5 present the same trends and lead to the same conclusions than in Section 5.1. The typical dispersion in the u and J0378 bands is ∼ 40% larger than in the MS case. This reflects the in-herent difficulties in the calibration of these passbands and their larger photometric errors. The final dispersion in the rest of the passbands is mildly larger by ∼ 10% with respect to the MS case in Section 5.1. We conclude that the zero points obtained with the MS stars also provide a good calibration for the photometry of the independent GB population. Thus, a proper calibration of any other astrophysical source in the images is expected, as is also demonstrated with the WD locus analysis presented in Sect. 4.4.

5.3. Colour uncertainties

The uncertainties in the colour calibration of the J-PLUS DR1 photometry are presented in Table 4. The modelling process in

Sect. 4.4 provides the best solutions for∆XWD and also their dispersions, typically ∼ 5 mmag. These errors must be added to the uncertainties in Sect. 5.1 to have the error in the calibration when X − r colours are analysed.

We further study the offsets implied by the WD modelling in the common PS1 filters giz. The r reference filter is dis-cussed in the next section. We found ∆gWD = −3 ± 2 mmag and∆iWD= 4 ± 2 mmag. In the case of the z band, we compared the final calibration zero point at each pointing estimated from the instrumental stellar and white dwarf loci (Sects 4.3 and 4.4), and by direct comparison with PS1 photometry (Sect. 4.2). We found a difference of 0 ± 5 mmag between both procedures.

Table 4. Error budget of the J-PLUS photometric calibration.

Passband σSLR σISL σISL_+P σWD σISL_+P+WD σcal [mmag]a [mmag]b [mmag]c [mmag]d [mmag]e [mmag]f

u 37 23 17 5 18 18 J0378 38 23 19 5 20 20 J0395 37 22 17 5 18 18 J0410 28 16 9 4 10 11 J0430 27 18 9 4 10 11 g 21 13 4 2 4 7 J0515 18 11 6 3 7 8 r 14 12 4 0 4 6 J0660 19 17 4 3 5 7 i 12 12 4 2 4 7 J0861 14 12 5 4 6 8 z 15 12 5 4 6 8

Notes.(a)_{Stellar locus regression (SLR) was used as calibration method. Uncertainty from duplicated MS stars in overlapping pointings.} (b)_{Instrumental stellar locus (ISL) or PS1 was used to homogenize the photometry. Uncertainty from duplicated MS stars in overlapping pointings.} (c)_{ISL (or PS1) and the plane correction were used to homogenize the photometry. Uncertainty from duplicated MS stars in overlapping pointings.} (d)_{Uncertainty in the colour calibration from the Bayesian analysis of the white dwarf locus.}

(e)_{Final uncertainty in the J-PLUS (X − r) colours, σ}2

ISL+P+WD= σ2ISL+P+ σ2WD. (f)_{Final uncertainty in the J-PLUS flux calibration, σ}2

cal= σ 2

ISL+P+WD+ σ2r, where σr = 5 mmag (Sect. 5.4).

4000 5000 6000 7000 8000 9000

Wavelength [ ˚A]

0

5

10

15

20

25

30

35

40

σ

cal

[mmag]

Fig. 11. Final calibration uncertainty in J-PLUS DR1. The black pentagons show the accuracy achieved with the calibration procedure presented in this work. The red dots show the accuracy of the Stellar Locus Regresion methodology used as reference in J-PLUS DR1. The dark (light) grey area marks a precision of 10 mmag (20 mmag).

photometry of the Pickles stellar library had colour residuals at ∼ 10 mmag level. We corrected the colour dependence of these residuals, but global offsets at such level can not be discarded. Thanks to the white dwarf locus, we have been able to estimate these global offsets.

5.4. Absolute flux uncertainty

The last source of error in our analysis is related with the abso-lute flux calibration, that is determined by the reference r band. We stress that any change in the r band calibration will modify accordingly the offsets ∆XWD to keep anchored the white dwarf locus. The colour offsets derived in giz from the PS1 initial cal-ibration are at 5 mmag level (Sect. 5.3), and we can assume a similar precision for the r band. Moreover, Narayan et al. (2019) found a 4 mmag offset between the PS1 photometry and their network of 19 WDs defined for calibration purposes. Thus, we assume a σr = 5 mmag uncertainty in the absolute flux

calibra-Table 5. Precision of the J-PLUS photometric calibration from giant branch stars.

Passband σSLR σISL σISL+P [mmag]a _[mmag]b _[mmag]c

u 45 26 22 J0378 40 26 26 J0395 44 25 19 J0410 32 17 12 J0430 29 18 10 g 22 13 5 J0515 20 11 8 r 15 13 5 J0660 19 17 6 i 14 12 5 J0861 16 12 6 z 15 12 6

Notes. (a) _{Stellar locus regression (SLR) was used as calibration}

method. Uncertainty from duplicated GB stars in overlapping tiles.

(b)_{Instrumental stellar locus (ISL) or PS1 was used to homogenize the}

photometry. Uncertainty from duplicated GB stars in overlapping tiles.

(c)_{ISL (or PS1) and the plane correction were used to homogenize the}

photometry. Uncertainty from duplicated GB stars in overlapping tiles.

tion of the reference r band. We present our total error budget for absolute flux photometry in the last column of Table 4 and in Fig. 11. When compared with the SLR uncertainty, reported in the first column of Table 4, a factor of two improvement is reached.

5.5. Impact of the assumed Milky Way extinction