Period-Luminosity Relations for Galactic classical Cepheids in the Sloan bands

We present the first period-luminosity (PL) and period-Wesenheit (PW) relations in the Sloan-Pans-STARRS gP1rP1iP1 bands for classical fundamental mode Cepheids in the Milky Way. We used a relatively modest number of 76 stars for the PL and 84-85 stars for the PW relations calibration. The data for the project were collected with the network of 40-cm telescopes of Las Cumbres Observatory, and Gaia Data Release 3 parallaxes were used for the calculations. These gri-band PL and PW relations calibrations will be a useful tool for distance determinations in the era of large sky surveys using the Sloan photometric system, especially with the near-future start of the Large Synoptic Survey of Space and Time (LSST).


INTRODUCTION
Since the discovery of a relationship connecting the periods and luminosities of classical Cepheids in the Small Magellanic Cloud (SMC) by Henrietta Swan Leavitt (Leavitt 1908;Leavitt & Pickering 1912), the periodluminosity (PL) relation (or Leavitt law) became one of the most important methods for estimating distances in the universe. The law serves as the first rung of the extragalactic distance ladder and is used to measure the distances to type Ia supernova host galaxies, which is a crucial step in the determination of the Hubble constant (e.g., Riess et al. 2022). For this reason classical Cepheids, together with other types of pulsating stars, such as type II Cepheids and RR Lyrae stars, became very important distance indicators. Astronomers dedicate considerable effort establishing ever more precise and accurate PL relations, as well as period-Wesenheit (PW) relations, where the Wesenheit index is a reddening-free magnitude by construction (Madore 1982).
PL and PW relations for classical Cepheids are available for a wide range of optical (e.g., Udalski et al. 1999;Tammann et al. 2003;Fouqué et al. 2007) and near-and mid-infrared bands (e.g., Macri et al. 2015; * Based on data from the Las Cumbres Observatory Ripepi et al. 2016;Wang et al. 2018), which have recently been extensively studied for their dependence on metallicity (e.g., Gieren et al. 2018;Ripepi et al. 2020;Breuval et al. 2021;Trentin et al. 2023;Molinaro et al. 2023) and also theoretically (e.g., Anderson et al. 2016;De Somma et al. 2022). Near-infrared PL relations have undeniable advantage over the optical range, mainly due to their decreased sensitivity to extinction, producing much smaller scatter. Nonetheless, the more symmetric shapes and smaller amplitudes of Cepheids at longer wavelengths, can cause the problem with their identification in the near-infrared domain (Pietrukowicz, Soszyński & Udalski 2021). This means that the optical range still plays an essential role for the discovery and classification of new variable stars, and precise PL relations are intensively studied for many optical filters. However, to our knowledge, the Sloan filters have not been widely used in the context of Cepheids so far, except for studies of variables in other galaxies, e.g., Hoffmann & Macri (2015) for NGC4258 or Adair & Lee (2023) for M33. The Sloan photometric system (ugriz;Fukugita et al. 1996) is a wide-band photometric system having wavelength coverage from 3000 up to 11000Å and was developed for the Sloan Digital Sky Survey (SDSS; Abazajian et al. 2003). Since its inception, this photometric system has gained wide-spread popularity for large-scale sky surveys. The most noteworthy are the Panoramic Survey Telescope And Rapid Response System (Pan-STARSS; Tonry et al. 2012), where the first part of the data release (Pan-STARRS1) of Andromeda (PAndromeda) provided PL relationships for pulsating stars in that galaxy (Kodric et al. 2018), and the Dark Energy Survey, using the DECam wide-field camera (Flaugher et al. 2015;Dark Energy Survey Collaboration et al. 2016). The Pan-STARSS and the Dark Energy Survey have become de facto standards for imaging in the Sloan bands for the northern and southern hemispheres, respectively.
Another ambitious upcoming project is the 10-year Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), which will be carried out at the Vera C. Rubin Observatory, currently under construction. The LSST is expected to provide an enormous amount of data and discover many new variable stars, among others, Cepheids in the Local Group galaxies and beyond. Hence, new possibilities are opening up in the field of distance measurements. However, in order to take full advantage of this project, precise calibration of PL relations for classical Cepheids (and other type of pulsating stars) has to be performed in the Sloan passbands, particularly in the Milky Way (MW), where good-quality geometric parallaxes for such objects are provided by the Gaia mission (Gaia Collaboration 2023). It is important to note, that unlike other photometric systems used before, generally, wide-field surveys in the Sloan system are not transforming their observations to the original reference system defined by Smith et al. (2002, where the magnitudes are often denoted as u ′ g ′ r ′ i ′ z ′ ), which is defined as the native system of the 1-m telescope of the US Naval Observatory Flagstaff Station. Instead, each survey releases their photometry in their own native systems, where the zero-points of the bands are calibrated using the transformed magnitudes from Smith et al. (2002). Therefore, it is a crucial detail to clarify exactly which system a given study is using.
Realizing the need for well-calibrated PL relationships in the Sloan bands, various PL and Wesenheit relationships have been established for different kinds of pulsating variable stars using data from the Zwicky Transient Facility (ZTF; Bellm et al. 2018Bellm et al. , 2019Dekany et al. 2020). These include RR Lyrae variables (Ngeow et al. 2022a), type II (Ngeow et al. 2022b) and anomalous Cepheids (Ngeow et al. 2022c), and even SX Phe variables (Ngeow et al. 2023). However, due to the relatively low number of calibrators for type II and anomalous Cepheids, as well as problems with blending due to the usage of stars in globular clusters and the large pixel size of the ZTF camera, significant systematic uncertainties remain in these calibrations. Furthermore, there is still no such calibration available for classical Cepheids, given that they are generally not located in clusters, and the vast majority of the ones with reliable Gaia parallaxes are brighter than the ZTF saturation limit.
The aim of this work is to derive the PL and PW relations for MW classical Cepheids in the three Sloan gri filters, and specifically calibrated to the Pan-STARRS system, for the first time. We expect that the relationships obtained here will become a useful tool for measuring distances in the universe in the upcoming era of large programs like the LSST.
The paper is organized as follows. In Sect. 2 we describe the sample of the selected classical Cepheids, the reduction of the data and further selection criteria for the derivation of the PL and PW relations, described later in Sect. 3. A short discussion is given in Sect. 4, and we conclude with a summary of our main results in Sect. 5.

Sample of stars
For the purpose of deriving the PL relations we observed 96 classical Cepheids pulsating in the fundamental mode. The chosen stars are bright, with the great majority having Gaia G-magnitudes in the range from about 6 to 12 mag with a few fainter stars up to 15.7 mag, which were mostly observed serendipitously, as they lie in the same field as a bright (primary) target. They are located at distances between 0.57 and 12.48 kpc (with a median distance of 2.3 kpc), which assure good-quality Gaia DR3 parallaxes (Gaia Collaboration 2023). Their periods taken from Pietrukowicz, Soszyński & Udalski (2021), ranging from about 3 to 69 days. The locations of the selected Cepheids are marked in Fig. 1 and their main physical parameters are summarized in Table 1.

Data and reduction
The data were collected between May 2020 and July 2022 within various observing programs 1 . The images of target stars were obtained in three filters modeled after the original Sloan g ′ r ′ i ′ filters, using 16 robotic 40-cm telescopes, which are part of the Las Cumbres Observatory (LCO) Global Telescope Network 2 . They are equipped with 3K × 2K SBIG STL-6303 cameras, 0°7 5°b 2 4 6 8 10 12 Distance [kpc] with a field of view of 29.2 × 19.5 arcmin 2 and pixel size of 0.571 arcsec pixel −1 without binning. Observations were taken in the air mass range of 1.00−1.87, and measured average seeing of about 2.35, 2.26 and 2.30 arcsec in the Sloan g, r and i bands, respectively. We aimed to obtain well-covered light curves for each of our targets. Therefore, during the course of the project, their observations were continuously monitored, and new observing requests were sent to LCO to fill gaps in the light curve phases. The raw images were reduced and processed with the LCO BANZAI pipeline 3 and we accessed them from the LCO Archive 4 . We performed aperture photometry on these images with a fixed aperture of 14 pixels, corresponding to about 8 arcsec (which for the vast majority of images contain virtually all light from the stars) using the standard DAOPHOT package (Stetson 1987). Next, we cross-matched the list of stars from each image with the ATLAS All-Sky Stellar Reference Catalog version 2 (ATLAS-REFCAT2; Tonry et al. 2018), which provides a compilation of magnitudes for stars between magnitudes 6 and 19 in the Pan-STARRS system (for the bands used here, generally denoted as g P 1 r P 1 i P 1 ). The average DAOPHOT photometric errors were 0.04 mag for Sloan gr and 0.05 mag in Sloan i for stars with Sloan g < 14.0 mag, which are the magnitudes of most of our local reference stars. The average DAOPHOT photometric errors for the Cepheids themselves are 0.002 mag in all three filters.
In order to obtain the light curves of the Cepheids in our fields, first we identified any possible variable stars around our target Cepheids, and removed them from the list of potential comparison stars. From the remaining ones, we selected stars with photometric uncertainties better than 0.1 mag and up to a few magnitudes fainter than the target Cepheid in the field. We then accepted them as the reference points for standardization that we used to calculate the coefficients in the following transformation equations: g inst = a 1 · g cat + a 2 + a 3 · (g − r) cat (1) where g inst , r inst and i inst are our Sloan instrumental magnitudes normalized to one second, while g cat , r cat and i cat are the catalog magnitudes from Tonry et al. (2018) in the Pan-STARRS system. The colors in the equations were calculated based on three subsequent images in each filter, which was possible as the observation strategy was to perform consecutive exposures in all three bands for all epochs. We calculated the a 3 , b 3 , and c 3 coefficients based on the constantmagnitude stars from all fields for a given telescope, and then fixed them for that telescope, in order to reduce the number of derived coefficients. The a 2 , b 2 , and c 2 coefficients are the zero-points of the transformations. The a 1 , b 1 , and c 1 coefficients were added in order to remove the trend that appeared in the residuals when they were not present in the equations. The reason for the appearance of this trend is the presence of nonlinearity affecting the SBIG 6303 cameras, as noted before by the LCO staff 5 . While a procedure was proposed to minimize it, the authors of the report pointed out themselves that the form of the proposed correction is not optimal. We noticed that adding the a 1 , b 1 , and c 1 coefficients in our transformation equations removed the trend well, so we decided to stick to this solution. The number of used comparison stars in a given image ranged from a few to several dozens. The fitting procedure was done iteratively by removing deviating stars in each of three iterations (first by using the Random Sample Consensus algorithm from the scikit-learn Python package, and then two 3σ clipping steps). In the end, the residuals were inspected by eye, and single outliers were discarded if necessary. The derived coefficients from Equations (1)-(3) were then used to transform the brightnesses of the Cepheids in a given frame to the Pan-STARSS photometric system. The typical rms value of the transformations was about 0.035 mag, but depended on the quality of images and the number of comparison stars used. We also performed a bootstraping procedure with 1000 iterations for Equations (1)-(3) to estimate the errors of the standardized magnitudes of the Cepheids in each frame, which we utilized in the next step. The mean shift between our calibrated and catalog magnitudes in the range of our Cepheid magnitudes in the three Sloan filters is less than 0.01 mag for all but one camera. Thus, this value can be adopted as the maximum error of the accuracy of the photometric zero-point, which propagates to the final PL/PW relations as a systematic error.
For each star in each of the three Sloan filters we phased the data, and subsequently calculated intensityaveraged mean apparent magnitudes by fitting the Fourier series to the Cepheid light curve with the curve fit Python algorithm. The applied Fourier order depended on the properties of a given light curve and was chosen by visual inspection. The minimum chosen order was equal to two and the maximum used was 14. High orders were allowed to fit well the maximum brightness in specific cases. During the fitting procedure we applied weights to each point, which we calculated from the errors determined from bootstraping in the previous step. Additionally, we performed Monte Carlo simulations, but errors estimated with this method were of an order smaller than the ones returned from Fourier series fit, so we conservatively adopted the latter as the uncertenties of obtained mean apparent magnitudes. The final light curves with fitted Fourier series are presented in Fig. A.1.

Reddening
The reddening of Cepheids is a complicated topic and variously treated in the literature (e.g., Harris 1981;Turner 1989;Kovtyukh et al. 2008; see Turner 2016, and references therein for a review; Madore, Freedman & Moak 2017). Multiple extinction maps have been published for our Galaxy, e.g., Schlegel, Finkbeiner & Davis (1998, hereafter SFD map) or the more recent Gaia/2MASS 3D Dust Extinction Map (Lallement et al. 2019) and Bayestar2019 reddening map (Green et al. 2019) available through the dustmap code 6 , all having their advantages and disadvantages. For the Cepheids in our sample, however, those maps turned out to be inadequate. Even though the SFD map provides reddening values for all our target stars, it turned out to be overestimated in most of the cases even after accounting for the simple model of dust in the Galaxy in the line of sight. As for the Gaia/2MASS 3D map, it is not available for about a half of our stars, as they are located at distances larger than 3 kpc which is its limit. Similarly, the Bayestar2019 reddening map did not provide reddening values for about a half of the objects, due to a lack of coverage of specific Galactic regions. Despite the abovementioned problems, we have tested the maps but they generated an artificially large scatter in the resulting PL relations. For this reason we have decided to follow the approach from Breuval et al. (2021) and adopt reddening values from Fernie et al. (1995) 7 , after application of the 0.94 scaling factor suggested by Groenewegen (2018). Also, we have adopted an uncertainty of 0.05 mag if not provided. Reddening information was only available for 86 stars from our sample (see column 6 in Table 1), but as we still had far enough objects for the construction of precise PL relations, we decided against the idea of mixing reddenings from different sources to keep consistency.
And finally, we adopted extinction vectors (R λ ) for the Pan-STARRS1 g P 1 r P 1 i P 1 filters from Green et al.
(2019, their Table 1), equaling: R g = 3.518, R r = 2.617 and R i = 1.971. We also calculated R λ directly from the Fitzpatrick (1999) reddening law using the extinction 8 Python package, assuming R V = 3.1, and effective wavelengths for the Pan-STARSS filters from Table 4 of Tonry et al. (2012), and obtained R g = 3.666, R r = 2.583 and R i = 1.901.

Distances
In order to calibrate PL relations, one needs to recalculate the derived apparent mean magnitudes of stars into their absolute magnitudes, which requires knowing the distances to the stars. The best source of geometric parallaxes known today is the third Gaia Data Release (hereafter Gaia DR3) catalog (Gaia Collaboration 2023). Lindegren et al. (2021a) proposed a correction for the zero-point (ZP) offset of Gaia parallaxes, which can be calculated with the dedicated Python code 9 they provided. This correction takes into account the ecliptic latitude, magnitude, and color of a given star. The parallaxes of our Cepheids range between about 0.08 to 1.76 mas (see column 3 in Table 1), and the corrections from about −9 to −52 µas (with a mean of about −29 µas). Following the recommendation of Lindegren et al. (2021a) to include an uncertainty of a few µas in the ZP, we adopted 5 µas as a systematic uncertainty. This value is propagated into the distance modulus, and hence the PL relation determination, as 0.025 mag for the median parallax for our sample, equal to 0.43 mas.
The quality of the parallaxes can be judged based on two parameters: the Renormalized Unit Weight Error (RUWE) and the Goodness-of-fit (GOF). The RUWE indicator is sensitive to the photocentric motions of unresolved objects, such as astrometric binaries. Ideally, it should be equal to 1.0 (Lindegren et al. 2021a). The GOF parameter was used by, e.g., Riess et al. (2021) as an indicator of the level of asymmetry of a source, and was adopted by the authors to be less than 12.5. Following Breuval et al. (2021) and Wielgórski et al. 8 https://extinction.readthedocs.io/en/latest/ 9 https://gitlab.com/icc-ub/public/gaiadr3 zeropoint (2022) we rejected stars with RUWE > 1.4 (seven stars) or GOF > 12.5 (five stars) from our sample of Cepheids. All stars that met the second condition also met the first, which means that RUWE is more constraining than GOF in our case. The rejected stars were FN Vel, RZ CMa, SU Cru, U Aql, UW Car, V496 Aql, and WW Car (marked in Table 1 with the "a" footnote in their IDs). Moreover, we have rejected from the fitting two additional stars (BB Sgr and IT Car) because of the poor quality of their light curves (marked in Table 1 with the "b" footnote, see also Fig. A.1).
Five of our Cepheids fall into the problematic range of transition of Gaia window classes where the parallax ZP could be affected (see Fig. 1 in Lindegren et al. 2021b). The three of them fall into G = 11.0 ± 0.2 mag and two into G = 12.0 ± 0.2 mag. Following Breuval et al.
(2021), we quadratically added 10 µas to the parallax uncertainties of the affected stars. Finally, following the suggestion of Riess et al. (2021) we increased all Gaia DR3 parallax uncertainties by 10% to account for a possible excess uncertainty. The mean of the calculated parallax uncertainties (given in column 3 in Table 1) is 21 µas.

Period-Luminosity Relations
To calibrate the PL relations for a sample of selected classical MW Cepheids in the three Sloan filters, g P 1 r P 1 i P 1 , first we dereddened the apparent mean magnitudes (m λ ) derived from the Fourier series fitting as described in Sec. 2.3, and then calculated absolute magnitude (M λ ) of each Cepheid as: where ϖ is the parallax in arcseconds. Next, we plotted the derived absolute magnitudes against the decimal logarithm of the periods related by the following relationship: where a λ and b λ are the slope and intercept we are looking for, respectively. The pivot period logP 0 = 1 was added to the original equation in order to minimize correlation between the slope and the intercept. It was adopted to be equal to 1, because we wanted to use the value close to the median of our Cepheid sample (∼ 0.98 days) but simultaneously being a round number. We used the curve fit function from the scipy Python library for the fitting in each passband with 3σ-clipping and obtained the slope and intercept from Equation (5) together with their corresponding errors, which we adopted as the final statistical uncertainties of the calculated coefficients. The obtained values are given in Table 2, where 76 stars were used for the final fit from the range of distances from 0.6 to 7.8 kpc (with a median of about 2.2 kpc). As a control test, we also performed a bootstraping procedure on the residuals with 10, 000 iterations in each passband and got very similar values for the uncertainties compared to the ones from the least-square method, which assures us that the adopted slope and intercept errors are reliable. The resulting rms of the PL relation is 0.25, 0.23, and 0.22 in the Sloan -Pan-STARRS g P 1 r P 1 i P 1 bands, respectively, for Green et al. (2019)  To avoid any bias in the absolute magnitude due to its nonlinear relationship with parallax, we also calculated the Astrometry-Based Luminosity (ABL), following the approach recommended by Feast & Catchpole (1997) and Arenou & Luri (1999), where: Then the PL relation given by Equation (5) will take the form: We performed the fitting of Equation (7) as described above and the resulting coefficients with their corresponding uncertainties are given in Table 2.
The PL relations for the three filters are presented in Fig. 2. The error bars in the plot are derived from the error propagation, taking into account the statistical errors on the mean magnitudes and parallaxes, where the latter ones are dominant. Solid and dashed lines represent the fit from the classical and ABL methods, respectively. The coefficients from these two approaches differ slightly, and the difference is higher for slopes (within 2σ), while the intercepts are within 1σ.

Period-Wesenheit Relations
PL relations based on reddening-free Wesenheit indices (or sometimes called Wesenheit magnitudes, originally introduced by Madore 1982) have proved to be particularly useful for distance determination due to their independence of reddening. Following Ngeow et al. (2022b) we defined three Wesenheit indices, W ri r , W gr r , and W gi g , for which, however, we calculated coefficients based on R λ from two sources: Green et al. (2019) and Fitzpatrick (1999); see Section 2.3. The resulting equations are given in column 2 of Table 3.
The mean apparent magnitudes used for deriving the Wesenheit indices were not corrected for reddening. We then fitted the relations analogous to Equation (5) and Equation (7): where again a and b are the slope and intercept sought, respectively, and logP 0 is the pivot period, chosen to be 1, and we applied the 3σ-clipping procedure to remove outliers. Because PW relations do not require reddening values, the number of stars used for calibration was higher than for the case of absolute magnitudes, and ranged between 84 and 85 stars (from a distance range between 0.6 and 12.5 kpc, with a median of 2.4 kpc). The smallest rms we found was for the R λ calculated based on Fitzpatrick (1999, rms = 0.21 − 0.27) and the largest was based on Green et al.
(2019, rms = 0.27 − 0.30). We performed the fitting for the Wesenheit coefficient with the ABL method, as well. The coefficients with their corresponding uncertainties returned by the least-squares algorithm are given in Table 3. Figure 3 presents an example of the fit for parameter R λ from the Wesenheit index calculated based on the extinction coefficients from Green et al. (2019). Also in this case the slope values from the classical and ABL methods differ noticeably, while the intercepts agree very well within the errors.

Stars rejected from the PL/PW relation fit
As described in Section 2.4, seven stars from the original sample were rejected from the fitting because they have not met the criteria of the Gaia RUWE and GOF quality indicators. We marked these stars in Figures 2 and 3. Six of them (FN Vel, RZ CMa, U Aql, UW Car, V496 Aql, and WW Car) lie very close to the found relations, but the obvious outlier is SU Cru, which deviates significantly toward brighter absolute magnitudes. This variable is defined in the Simbad Astronomical Database 10 as a double or multiple star. Szabados (1996) listed this star as one of the spectroscopic binaries found among classical Cepheids, which has a red photometric companion. Additionally, we have also rejected BB Sgr and IT Car because of the bad quality of their light curves (see Figure A.1).
After 3σ-clipping, AD Cam was rejected from the fit of the PL relations for being too faint. It is also deviating in the PW relations, although not so significantly.
The light curve of that star presented in Figure A.1 is not very well covered, which could affect the quality of the applied Fourier series and consequently the mean magnitude. Supplementing the data for this star will possibly eliminate this problem in the future.
Apart from AD Cam, EW Car is also rejected after 3σ-clipping in all PW relations. There is no information about the reddening for that Cepheid in Fernie et al. (1995), so it was not used for fitting the PL relation. Although the fit of the Fourier series for EW Car (Figure A.1) looks satisfactory, the quality of its light curve leaves much to be desired, as it was observed incidentally in the field containing the much brighter AQ Car, meaning the exposure times were kept short to not saturate that target, resulting in lower-quality light curves for EW Car. Ideally, a bigger telescope should be employed to provide data of higher quality in the future. Nevertheless, the current data are good enough to obtain good values of the Wesenheit indices, which might suggest some reddening issue for EW Car.
In the case of W gi g calculated based on the extinction coefficients from Fitzpatrick (1999) besides AD Cam and EW Car, SV Vul also has been rejected, for being too bright. This long-period Cepheid is another case of a spectroscopic binary listed by Szabados (1996). They also list as spectroscopic binaries six additional stars common with our sample: AV Sgr, LS Pup, RT Mus, RU Sct, VW Cen and YZ Sgr. If we remove them from the sample then the slopes and intercepts of the PL and PW relations would still agree within the 1σ uncertainties. The identified binary Cepheids constitute about 7% to 9% of the sample of stars used for the PL and PW relations determination, respectively. As Karczmarek et al. (2023) show, the influence of binary Cepheids on the distance modulus (i.e. the difference between the zero-points of the PL relations in two galaxies) for this binarity fraction is very small for galaxies having similar percentages of binary Cepheids (see their Figure 7), and can be neglected. This means that there is no need to reject these stars from our sample.

Test for first overtone contamination
Although all Cepheid variables used by us in this analysis have been classified as fundamental mode pulsators (Pietrukowicz, Soszyński & Udalski 2021), we performed a test to confirm this. To avoid any significant contamination with first overtone Cepheids, we have adopted a period cutoff of 5 days, which corresponds to logP = 0.7 days, and repeated the PL/PW relation fitting. The new slopes and zero-points agree within the 1σ uncertainties with the values from Ta-bles 2 and 3 which assures us about the purity of our sample of stars regarding their pulsation modes.

The influence of the parallax ZP on the PL/PW relations
To evaluate the influence of the ZP offset calculated based on Lindegren et al. (2021a) on the found PL/PW relations we repeated the fitting without introducing parallax corrections as a test. The slopes without a ZP correction are steeper by slightly more than 1σ, while the ZPs differ much more, and they are, on average, smaller by about 0.12 mag for the PL relations and about 0.14 mag in the case of the PW relations. The difference might be caused by the fact that the latter includes Cepheids from larger distances than the former.
The same year when Lindegren et al. (2021a, hereafter L21) introduced their ZP offset corrections of Gaia parallaxes, also Groenewegen (2021, hereafter G21) published his own, independent recipe. Both works use quasars and wide binaries as sources with known parallaxes for their calibrations. The main difference between the two approaches are the selection criteria: L21 use the G-magnitude, pseudocolor, and ecliptic latitude of a given star, while G21 decided to use G-magnitude, (G BP − G RP ) color, and for the spatial component HEALPix formalism was used to transform the sky coordinates into a sky pixel (where the HEALPix level is a parameter defining its resolution). We have calculated new parallax corrections according to the G21 instructions, where the mean new correction was about −20 µas, and have repeated the fitting of the PL/PW relations for our Cepheids. We obtained quite good agreement in the resulting coefficients with the ones calculated with the L21 corrections. The new slopes are slightly steeper, but they are well within the 1σ uncertainties, particularly for the absolute magnitudes and W ri r . For W gr r and W gi g , however, the differences are higher (within about 2σ). The intercepts for all relations are slightly smaller for G21, with the ZP shift being about 0.02 mag, which agrees within 1σ with L21.

Dependence on metallicity?
To check for any possible trend with metallicity, we adopted the metallicity values for the Cepheids in our sample from the summary provided by Breuval et al. (2021) in their Table 3. We then plotted the residuals of our PL and PW relations from Figures 2 and 3, respectively, against metallicity values (see Figure 4(a) and (b)). For the PL relations, metallicity information was available for 65 stars, and for the PW relations for 68 stars. For both cases, the metallicity range used is between −0.30 and 0.55 dex, with an average and standard deviation of 0.09 ± 0.16 dex. The comparison of the PL/PW relations residuals with metallicity does not show any clear trends in any band, but it should be borne in mind that the range of metallicity used is rather narrow, and trends could possibly be revealed for a wider range of metallicities. Because of the small range of the metallicities of our Cepheids, we also decided not to investigate a Period-Luminosity-Metallicity (PLZ) relation, however, such an analysis will be the subject of future projects.

Slope versus wavelength
The comparison of the slopes obtained by us using the two methods (linear and ABL) with other slope values presented in the literature, given in Figure 5, shows general agreement with both theoretical and observational results. The general trend is that the slopes become steeper for longer wavelengths. This is the behavior predicted by theoretical models (e.g., Bono et al. 2010), also observed in the data for any system. Di Criscienzo et al. (2013) provided a theoretical scenario for the properties of classical Cepheids in SDSS filters. They used sets of nonlinear convective pulsation models computed for the typical chemical compositions of the Galaxy, transformed subsequently into the the SDSS photometric system. They noted the discrepancy of their resulting slopes with semiempirical and empirical ones from the literature, with the latter being systematically shallower. This is opposite in our case, as our results are systematically steeper than the theoretical ones. Simultaneously, the slopes obtained by us for the Pan-STARSS bands are shallower than those for the Johnson-Cousins bands (see Fouqué et al. 2007;Breuval et al. 2021Breuval et al. , 2022, which is, however, predicted by the aforementioned theoretical studies. The comparison of our slopes with the slopes for all fundamental mode classical Cepheids from M31 derived by Kodric et al. (2018) in the Pan-STARSS photometric system yields very satisfying agreement. Their results match our trends well, particularly the linear method, which was applied in that study.
Worth noting is also the fact that the slopes obtained by us with the linear fit are steeper than those obtained with the ABL method, and the difference is more significant for longer wavelengths. The differences in the slope values (of either method) in the r P 1 and i P 1 filters calculated based on the reddening vectors from Green et al. (2019) and Fitzpatrick (1999) are statistically insignificant, but they are largest for the g P 1 band. This behavior, however, is not surprising, as reddening laws differ the most from each other in the blue part of the spectrum.

Comparison of the PL relations with the literature
For the sake of comparison, we used the theoretical slopes from Di Criscienzo et al. (2013, see their Table 2) predicted for all stars with log P ≤ 2.0 in the SDSS photometric system (marked also with black triangles in Figure 5), applied our zero-points, and plotted them with black, solid lines in Figure 6 over our PL relations from Table 2. We also used the PL relations derived for M31 by Kodric et al. (2018) for fundamental mode Cepheids from the full range of periods given in their Table 3, shifted them by the distance modulus from Li et al. (2021) equal to µ 0 M 31 = 24.407 mag, and then overplotted with the magenta, solid line in Figure 6(a) and (b). Analogously, we plotted the PW relation for M31 available for W ri r and used the theoretical slope for W gi g over the corresponding relations derived by us in Table 3 (see Figure 6(c) and (d)).
Such a comparison shows very good agreement of our resulting PL relations with the theoretical predictions for the SDSS system (similar to Pan-STARSS), particularly in the r P 1 and i P 1 bands, while the g P 1 band shows the worst compliance, as our slopes in the case of the PL relations are noticeably steeper, while the W gi g slopes are shallower than those reported by Di Criscienzo et al. (2013). Very positive and encouraging is the comparison with the empirical relations for M31 obtained in the same Pan-STARSS photometric system, although the conclusion about the PL relations in the g P 1 band remains unchanged. This band shows the worst agreement, particularly for the reddening vectors obtained directly from the Fitzpatrick (1999) reddening law. Also the W ri r values calculated based on this law differ much from the W ri r values from Kodric et al. (2018), potentially resulting in smaller measured distances.
As a consequence, we recommend using the g P 1 band (and Wesenheit indices based on it) for distance scale purposes with caution. It is particularly sensitive to the applied reddening law, as well as characterized by the largest resulting uncertainties, among others, from the larger intrinsic spread of the instability strip in the blue bands. The r P 1 and i P 1 bands, however, are very promising tools for future applications, especially in the future LSST survey.

SUMMARY
In this work we presented a calibration of the PL relations in the Sloan -Pans-STARRS g P 1 r P 1 i P 1 filters and W ri r , W gr r , and W gi g Wesenheit indices for classical Cepheids from the MW pulsating in the fundamental mode. According to our knowledge, this was performed using MW Cepheids for the first time here. We obtained the data for our analysis using the 40-cm telescopes of the LCO Telescope Network in the three g ′ r ′ i ′ Sloan filters, used parallaxes from the Gaia DR3 catalog, and ATLAS-REFCAT2 (Tonry et al. 2018) to calibrate the photometry into the Pan-STARRS version of the Sloan photometric system. The obtained light curves, which we present in Figure A.1 in the Appendix A, are generally well covered in pulsation phase, which allowed us to determine reliable mean magnitudes for a sample of 96 Cepheids used for this analysis. The homogeneous reddening values were adopted from Fernie et al. (1995, which, however, were not available for all Cepheids from our sample) and the extinction coefficients from Green et al. (2019) and calculated directly from the Fitzpatrick (1999) reddening law. After application of the selection criteria on the Gaia RUWE and GOF parameters and rejecting stars with poor-quality light curves, the number of stars used for the PL relations was 76 (see Table 2) and for the PW relation between 84 and 85 (see Table 3). The Wesenheit magnitudes were calculated based on the total to selective extinction ratios from two sources: Green et al. (2019) and Fitzpatrick (1999).
The predominant contribution to the absolute magnitude uncertainties and subsequently the statistical errors of our PL/PW relations comes from the parallax uncertainties, while the mean magnitude errors are minor. The systematic uncertainties of the ZPs of our relations come from the Gaia parallax ZPs and were estimated to be 0.025 mag, and the photometric ZP uncertainties, which we adopted as 0.01 mag.
We have checked the behavior of the PL/PW relations in case of removing the confirmed binary Cepheids from our sample, and we did not notice any significant change.
This test led us to the conclusion that there is no need for rejecting those stars from the fit.
We performed a comparison of the residuals of our PL/PW relations with the metallicities of those Cepheids, for which they were provided in Table 3 in Breuval et al. (2021), and did not notice any trends.
Because of the small range of available metallicities, we have decided not to perform a fit to obtain the PLZ relation. To do such a fit properly, a much larger range of metallicities is required, preferably using Cepheids also from nearby galaxies, which will be the aim of projects in the future.
We compared our resulting slopes for the Galactic PL relations in the Pan-STARSS photometric system with theoretical slopes in the SDSS system (Di Criscienzo et al. 2013), as well as empirical ones in other passbands. We observe overall good agreement with the general trend that the slopes become steeper for longer wavelengths. Moreover, our PL relations, especially in the r P 1 and i P 1 bands, agree well with the empirical PL relations for M31 derived by Kodric et al. (2018), which are also in the Pan-STARSS photometric system, while the agreement for the g P 1 band is the poorest.
We believe the relations derived here will be a useful tool for future distance determinations in the era of the upcoming LSST project, although we keep in mind that improvements are still easily possible. This could be done by extending the sample of used Cepheids, particularly by long-period Cepheids; by collecting data using homogeneous observations; by improving the coverage of the light curves; and by improving the extinction estimates of Cepheids. Future Gaia data releases are also expected to provide parallaxes of even better quality. We thank the anonymous referee for valuable comments which improved this paper.

2
The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 695099).
We also acknowledge support from the National Science Center, Poland grants MAESTRO UMO-2017/26/A/ST9/00446, BEETHOVEN UMO-2018/31/G/ST9/03050, and DIR/WK/2018/09 grants of the Polish Ministry of Science and Higher Education. We also acknowledge financial support from UniverScale grant financed by the European Union's Horizon 2020 research and innovation programme under the grant agreement number 951549. We gratefully acknowledge financial support for this work from the BASAL Centro de Astrofisica y Tecnologias Afines (CATA) AFB-170002 and the Millenium Institute of Astrophysics (MAS) of the Iniciativa Cientifica Milenio del Ministerio de Economia, Fomento y Turismo de Chile, project IC120009. W.G. also gratefully acknowledges support from the ANID BASAL project ACE210002. P.W. gratefully acknowledges financial support from the Polish National Science Center grant PRELUDIUM 2018/31/N/ST9/02742. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/ gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.  (Tody 1986(Tody , 1993, DAOPHOT (Stetson 1987), Astropy (Astropy Collaboration 2013), Sklearn (Pedregosa et al. 2011), NumPy (van der Walt et al. 2011Harris et al. 2020), SciPy (Virtanen et al. 2020), and Matplotlib (Hunter et al. 2007). Note-Source: source of extinction vector (R λ ); Band: the Pan-STARRS g P 1 r P 1 i P 1 bands; a λ : slope of the fit; b λ : zero-point of the fit; rms: rms of the derived relations; N : number of stars used for fitting.