Dust Extinction Law in Nearby Star-Resolved Galaxies. II. M33 Traced by Supergiants

The dust extinction curves toward individual sight lines in M33 are derived for the first time with a sample of reddened O-type and B-type supergiants obtained from the LGGS. The observed photometric data are obtained from the LGGS, PS1 Survey, UKIRT, PHATTER Survey, GALEX, Swift/UVOT and XMM-SUSS. We combine the intrinsic spectral energy distributions (SEDs) obtained from the ATLAS9 and Tlusty stellar model atmosphere extinguished by the model extinction curves from the silicate-graphite dust model to construct model SEDs. The extinction traces are distributed along the arms in M33, and the derived extinction curves cover a wide range of shapes ($R_V \approx 2-6$), indicating the complexity of the interstellar environment and the inhomogeneous distribution of interstellar dust in M33. The average extinction curve with $R_V \approx 3.39$ and dust size distribution $dn/da \sim a^{-3.45}{\rm exp}(-a/0.25)$ is similar to that of the MW but with a weaker 2175 Ang bump and a slightly steeper rise in the far-UV band. The extinction in the $V$ band of M33 is up to 2 mag, with a median value of $ A_V \approx 0.43$ mag. The multiband extinction values from the UV to IR bands are also predicted for M33, which will provide extinction corrections for future works. The method adopted in this work is also applied to other star-resolved galaxies (NGC 6822 and WLM), but only a few extinction curves can be derived because of the limited observations.


INTRODUCTION
Interstellar dust efficiently absorbs and scatters starlight, affecting observations and physical processes. Dust extinction or dust attenuation is of vital importance to recover the intrinsic spectral energy distributions (SEDs) of celestial objects and infer the properties of dust. Extinction represents the amount of light lost due to absorption and scattering of dust along a sight line. The extinction at a given wavelength depends on the grain size distribution and the optical properties of the grains (Salim & Narayanan 2020). In contrast to extinction, attenuation depends on both extinction and the complexity of star-dust geometry in galaxies, including scattering back into the sight line, varying column densities or optical depths and the contribution by unobscured stars (Salim & Narayanan 2020). Cardelli et al. (1989, CCM hereafter) found that the dust extinction law in the Milky Way (MW) from ultraviolet (UV) to near-infrared (IR) bands could be characterized by one parameter named the total-to-selective extinction ratio R V [= A V /E(B−V )], which depends on the interstellar environment along the sight line. However, the CCM extinction law is limited to only a set of sight lines in the MW, and it is not generally applied to external galaxies (Clayton et al. 2015). The properties of dust extinction curves or dust attenuation curves in galaxies and the physical mechanisms that shape them are fundamental extragalactic astrophysics questions and are important for deriving the physical properties of galaxies (Salim & Narayanan 2020). On the one hand, external galaxies allow us to study dust in diverse interstellar environments, which is a necessary intermediate step to understanding distant galaxies. On the other hand, whereas interpretation can sometimes be difficult in the MW disk because we see the projected material of the entire disk, high-latitude observation of face-on galaxies can provide clearer sight lines (Galliano et al. 2018).
Although the average extinctions in the LMC and M31 are similar to that in the MW, the extinction curve in the bar region of the SMC rises steeply in the UV bands and lacks 2175Å.
For the late-type spiral M33 (Sc, Nilson 1973) (≈ 840 kpc Freedman et al. 1991), which is the third largest member in the Local Group Galaxies, the latest study on attenuation was carried out by Moeller & Calzetti (2022). Moeller & Calzetti (2022) combined archival images from UV to IR to derive the ages, masses, and the values of E(B − V ) for the young star cluster population in M33 and found that all the star clusters have moderate-to-small internal extinction [E(B − V ) < 0.6 mag]. Hagen (2017) imaged the galaxy from FUV to NIR and measured the spatial variation of the dust attenuation law in M33 for the first time. They found that the attenuation curves tend to be steeper and with an MW-like 2175Å bump between the arms in M33, while along the arms, the curves seem to be shallower with a weak 2175Å bump. The median attenuation curve derived in Hagen (2017) is quite steep with a 2175Å bump and is somewhat different from the fairly shallow attenuation curve with a strong 2175Å bump obtained by Gordon et al. (1999) in the M33 nucleus study. Hagen (2017) found a median value of extinction in V band A V = 0.53 mag, which is twice the value of the fairly small mean amount of dust extinction (A V ≈ 0.25 mag) derived from the star formation study in M33 by Verley et al. (2009) because of the different assumed stellar models and the lack of FIR observations in Hagen (2017). The dust attenuation laws derived in Hagen (2017)

and Gordon
With the improvement of the observation resolution, individual stars in M33 could be distinguished and their photometry information and spectral types could be obtained, providing us with a completely new prospect for exploring the extinction law toward individual sight lines in M33. Wang et al. (2022, Paper I hereafter) derived dozens of extinction curves toward individual sight lines in M31 with the combination of the intrinsic SEDs from the stellar model atmospheres and model extinction curves from the dust model. In this work, the method adopted in Paper I is also applied to calculate the extinction curves in M33. We select the bright O-type and B-type supergiants in M33 from the Local Group Galaxies Survey (LGGS, Massey et al. 2016) as the extinction tracers following Paper I. Using the photometry available online, the spectral energy distribution (SED) for each tracer from UV to near-IR is constructed, the details of which are shown in Section 2. The method of forward modeling the SED to obtain the dust extinction law is described in Section 3. Section 4 presents the extinction curves derived in this work and the discussion. Finally, our conclusions are summarized in Section 5.

DATA AND SAMPLE
As in Paper I, we selected the isolated O-type and B-type supergiants from the LGGS catalog (Massey et al. 2016) as the extinction tracers in M33 because supergiants are usually free of circumstellar dust and relatively bright (Shao et al. 2018;Liu et al. 2019). The LGGS catalog contains 146,622 stars in M33, of which 130 and 471 are confirmed to be O-type and B-type stars, respectively (Massey et al. 2016). In this work, the isolated O-type and B-type supergiants from the LGGS catalog are also selected as extinction tracers to explore the extinction law in M33. 25 O-type supergiants and 318 B-type supergiants constitute the extinction sample for M33 in this work. Because of the limitation of ground-based telescopes, OB associations or binaries may be identified as single OB stars. The optical images obtained from the Hubble Space Telescope (HST) in the F 475W and V bands (F 547W , F 555W , F 569W ) are adopted to check the reliability of the extinction tracers.
There are 205 tracers in the extinction sample can be found in the HST/F475W image or the HST/V image, of which 200 stars seem to be single stars in the HST images. While 5 sources are suspect because they overlap with other celestial objects and cannot be distinguished in the HST images.
As a result, we suggest that 98% of the isolated supergiants from the LGGS catalog are reliable for calculating the dust extinction law in M33. The V − R/B − V diagram and B − V /V diagram for all LGGS sources and the supergiants in the extinction sample are plotted in Figure 1. We construct the observed SED for each tracer using the photometric data from the LGGS catalog (Massey et al. 2016) in the U, B, V, R, I bands, the United Kingdom Infrared Telescope (UKIRT, Irwin 2013) in the J, H, K bands 1 , the Panoramic Survey Telescope and Rapid Response System release 1 Survey (Pan-STARRS PS1, Chambers et al. 2016) in the g, r, i, z, y bands, the XMM-Newton Serendipitous ultraviolet source survey (XMM-SUSS, Page et al. 2019) in the U V W 2, U V M 2, U V W 1 bands and the swift ultraviolet and optical telescope (Swift/UVOT, Yershov 2015) in the U V W 2, U V M 2, U V W 1 bands, as mentioned in Paper I. The selection criteria for these catalogs are also the same as those in Paper I.
Instead of using the photometry from the Panchromatic Hubble Andromeda Treasury (PHAT) Survey (Williams et al. 2014) adopted in Paper I, we obtain the photometry in M33 from the Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (PHATTER, Williams et al. 2021). The PHATTER survey (Williams et al. 2021) presents panchromatic resolved stellar photometry for 22 million stars in the Local Group dwarf spiral Triangulum (M33), derived from HST observations with the Advanced Camera for Surveys in the optical bands (λ F 475W = 0.473 µm, λ F 814W = 0.798 µm) and the Wide Field Camera 3 in the near-UV (λ F 275W = 0.272 µm, λ F 336W = 0.336 µm) and near-IR bands (λ F 110W = 1.120 µm, λ F 160W = 1.528 µm). The survey covers ∼ 14 square kpc of the sky and extends to 3.5 kpc from the center of M33. The PHATTER catalog is the largest stellar catalog for M33. We check the GST ("good star") quality for each photometric point of each tracer and select the photometry with GST flag = '0', which means that the source passes the GST criteria in Williams et al. (2021). PHATTER photometry can thus be applied for 2 O-type supergiants and 36 B-type supergiants.
In addition, we adopt UV data from the Galaxy Evolution Explorer (GALEX, Bianchi et al. 2017) in this work. GALEX (Martin et al. 2005) performed the first sky-wide UV surveys with different coverage and depth (Morrissey et al. 2007;Bianchi 2009), yielding observations in the following two broad bands: far-UV (FUV, λ eff ≈ 1528Å) and near-UV (NUV, λ eff ≈ 2310Å) (Bianchi et al. 2017). Unfortunately, the M33 sky area is not completely covered by the observation of the GALEX.
In addition, it is estimated that most of the tracers in this work could be too faint for the GALEX with the typical depth of m F U V = 19.9 mag and m N U V = 20.8 mag to detect 2 . As a result, the UV data from the GALEX can be obtained for only a few tracers. We first check the artifact flag 2 Regardless of the dust extinction, the observed AB magnitudes in the F U V and N U V bands for a tracer with the median spectral type of B2 can be estimated with T eff = 18000 K, log(g) = 2.50 and d (distance) = 840 kpc. The derived values of m F U V = 19.84 mag and m N U V = 19.90 mag plus the effect of dust extinction could be greater than the typical depth of the GALEX. and the extraction flag for each photometric point of the tracers and eliminate spurious sources. We then eliminate the foreground photometry, which is too bright to fit the whole SED well. Finally, the GALEX data for only 1 O-type supergiant and 2 B-type supergiants are retained.
Above all, at most 27 bands of photometric data from UV to near-IR are obtained for each star.
We summarize the selection criteria and the number of photometric points adopted for each catalog in Table 1.

METHOD
For star-resolved galaxies, the pair method (Bless & Savage 1970) is extensively adopted to obtain the extinction law, which compares the spectrum of a reddened star with that of an unreddened (or slightly reddened) star with the same spectral type. In order to eliminate the influence of the limited unreddened standard stars and the mismatch error in the use of the pair method, Fitzpatrick & Massa (2005) proposed to use the stellar model atmospheres to derive the intrinsic SEDs rather than the unreddened standard stars.
Based on this "extinction without standards" technique, we first combine the intrinsic SED from the stellar model atmospheres extinguished by the model extinction curves to construct the model SEDs for the tracers, and then derive the extinction curves by fitting the model SEDs to the observed data.
Instead of the mathematical extinction models such as CCM (Cardelli et al. 1989), FM90 (Fitzpatrick & Massa 1990) and F04 (Fitzpatrick 1999(Fitzpatrick , 2004Fitzpatrick & Massa 2007) extinction laws that are widely used in many works, the classic silicate-graphite dust model is adopted to model the dust extinction law as Paper I, so that the dust properties can also be analyzed besides obtaining the extinction curve. In addition, the extinction curves derived from the dust model are more applicable in various interstellar environments than the parameterized extinction curves.
The detailed calculation process in this work can be found in Figure 3 of Paper I. The construction of the model SEDs for M33 is described in detail in Section 3.1, and Section 3.2 describes the fitting of the model SEDs to the observed data.

Model SEDs
Theoretically, the observed SED F obs λ of a reddened star can be expressed as follows: where F int λ is the intrinsic surface flux of the star at wavelength λ, θ ≡ (R/d) 2 is the angular radius of the star (where d is the distance and R is the stellar radius), and A λ is the absolute extinction/attenuation of the stellar flux by intervening dust at λ (Fitzpatrick & Massa 2005).
With 27 effective temperature values (27500 K ≤ T eff ≤ 40000 K with 2500 K steps for O-type supergiants, 10000 K ≤ T eff ≤ 30000 K with 1000 K steps for B-type supergiants), 7 surface gravity values (3.00 ≤ log g ≤ 3.50 with 0.25 dex steps for O-type supergiants, 2.25 ≤ log g ≤ 3.00 with 0.25 dex steps for B-type supergiants) and the solar value of the metallicity, a grid of intrinsic SEDs for M33 is constructed.
The model extinction curves in this work are also derived from the silicate-graphite dust model with the same exponential cutoff power-law grain size distribution proposed by Kim, Martin and Hendry (hereafter KMH, Kim et al. 1994) with a c fixed to 0.25 µm [dn/da ∼ a −α exp(−a/0.25)] for both components, as adopted in Paper I. Detailed information on the chemical abundances of the interstellar environment in M33 is still controversial, and it is difficult to quantify the chemical abundances of the dust in M33. However, some studies (e.g., Magrini et al. 2007a,b;Toribio San Cipriano et al. 2016;Ren et al. 2019Ren et al. , 2021a suggest that the abundances in M33 are close to the protosolar values (Asplund et al. 2009), which are adopted in some recent works (e.g., Neugent et al. 2017;Neugent 2021;Ren et al. 2021a). As a result, according to the previous works in other galaxies [e.g., MW ), M31 (Draine et al. 2014; Paper I), NGC 4722 (Gao et al. 2020), etc], we assume that the interstellar abundances in M33 are similar to the protosolar values (Asplund et al. 2009) and adopt a typical value of f cs = 0.3 for the mass ratio of graphite to silicate, which means the elements of Fe, Mg and Si are all in the solid phase and constrained in silicate dust, and the fraction of gas-phase carbon is 50% 3 . We then derive a grid of the model extinction curves with 56 values of α (0.5 ≤ α ≤ 6.0 with 0.1 dex steps) and 51 values of By combining the intrinsic SEDs and the model extinction curves, a large grid of monochromatic flux extinguished by dust can be derived as follows: To compare the model SEDs with the observed photometry, the model band flux can be calculated as follows: where B i (λ) is the bandpass response function for the ith band. The flux for each band is thus obtained from the response function and model SEDs, including the intrinsic SEDs and the model extinction curves.

Fitting Model SEDs to Observed Data
Since the tracers in M33 are paired with the stellar model atmospheres, the foreground MW dust extinction must be removed. We adopt an MW foreground extinction component of E(B −V ) ≈ 0.06 mag (Ruoyi & Haibo 2020), assuming an R V = 3.1 CCM dust, as a part of the fitting process.
The EMCEE fitting code (Foreman-Mackey et al. 2013) is used to fit the model SEDs to the observed data. It is a Markov-Chain Monte Carlo (MCMC) ensemble sampler and helps obtain the most suitable parameters and the corresponding confidence intervals. Gaussian likelihood is adopted, and flat priors are imposed on α, log(g) and A V in this work. A Gaussian prior is imposed on the effective temperature log(T eff ) based on the spectral type from the LGGS catalog (Massey et al. 2016) with one subclass in the spectral type considered the uncertainty (Clayton et al. 2015; Paper I). The calibration of the spectral type to log(T eff ) and log(g) refers to Cox (2000) and Conti et al. (2008).
With the model SEDs fit to the observed data, the fitting parameters for each tracer, including α in the dust model, extinction in the V band A V , the effective temperature log(T eff ) and the surface gravity log(g), can be derived, as well as the corresponding extinction curve, the average dust radius a, the color excess E(B −V ), and the total-to-selective extinction ratio R V . The average dust radius a is derived based on equation (8) in Nozawa (2016).
The results for one star in the extinction sample are plotted in Figure   2 as an example to show the comparison of the best-fitting model SED to the observed data with fitting parameters α, A V , the derived parameters R V , E(B − V ) and a marked, the corresponding extinction curve and the comparison of normalized model SED and normalized intrinsic SED.

Results Screening
As mentioned in Paper I, pcFlag is also introduced in this work to show the coverage of passbands adopted in the calculation for each tracers and the sample of the extinction tracers in this work is therefore divided into four subsamples (pcFlag = 'UVI ', pcFlag = 'UV ', pcFlag = 'VI ', pcFlag = 'V ') 4 . The numbers of the tracers in the four subsamples are 58, 20, 73 and 192, respectively. It is considered that the results for the sight lines with pcFlag = 'UVI ' are the most reliable, which are adopted to analyze the extinction curves in this work. The criteria used to select the reasonable results are similar to those in Paper I, as follows: III. The derived total-to-selective extinction ratio is in the range of R V = 1.5 − 7. In our extinction sample, there are also some reddened stars with good MCMC performance, yet the derived values of R V are unusually as high (R V > 7) as those in Paper I. The results may lead to acceptable models that fit the observations well, but the derived values of R V are unphysical because the values of R V are observationally in the range of 2 R V 6 (Mathis 1990;Welty & Fowler 1992;Fitzpatrick 1999;Draine 2011;Wang et al. 2017).
There are 58 tracers with pcFlag = 'UVI ' in the extinction sample, of which 39 tracers are selected for the further analysis based on the selection criteria mentioned above. It is anticipated that sufficient data covering UV to near-IR bands will bring more reliable results, because UV and near-IR data can effectively constrain the extinction curve and the intrinsic spectra. Although the results with pcFlag = 'UVI ' are recommended in this work, we also present the selected results for the tracers from the other three subsamples in the following section (see Table 2 and Table 3) as a reference.

The Extinction Curves in M33
The results for each selected tracer are partly listed in Table 2, of which the entirety is available in machine-readable form. The ID and spectral type for each tracer are extracted from the LGGS catalog (Massey et al. 2016). The fitting parameter [α, log(T eff ), log(g) and A V ] and the uncertainties are derived based on the 50%, 16% and 84% values of the parameter spaces generated from the EMCEE results. The corresponding values of a, E(B−V ) and R V derived in Section 3.2 are also listed in Table   2. The columns named "Total bands" and "pcFlag" in Table 2 present the number of passbands and the coverage of the passbands adopted in the calculation. The columns following the "pcFlag" column are the photometric information in each band, which are available in a machine-readable format. It should be noted that although extinction curves with smaller R V values usually have stronger 2175Å bumps and steeper far-UV rises, the CCM extinction law with only one parameter R V is not generally applied in external galaxies (Clayton et al. 2015); thus, R V cannot completely describe the extinction features on the extinction curves in external galaxies (Paper I). As a result, we prefer to adopt the parameter α in the dust size distribution function to describe the dust extinction and dust properties in M33 in this work. Table 3 summarizes the median values with the upper and lower limits of each parameter. As mentioned in Section 3.3, we divide our extinction sample into four subsamples based on the coverage of passbands adopted in the calculation (pcFlag). Lines 1 to 4 in Table 3 present the results of the four subsamples, and the fifth line shows the results of all the selected tracers in the extinction sample.
As illustrated in Paper I, the results for the sight lines with pcFlag = 'UVI' are the most reliable and are consequently adopted to describe the general extinction law in M33. Lines 6 to 9 in Table   3 show the influence of the lack of UV or near-IR data on the results (see Section 4.4 for details).
The last line in Table 3 is the results of fitting our model extinction curves derived directly from the silicate-graphite dust model to the MW extinction curve (Fitzpatrick et al. 2019, F19 hereafter) for comparison.
The extinction curves toward the sight lines with pcFlag = 'UVI' are plotted in Figure 3 with gray solid lines. The extinction curves in M33 cover a wide range of shapes, from flat extinction curves with large R V to steep curves with obvious 2175Å bumps, indicating the inhomogeneous interstellar environment and dust distribution in M33. The red solid line in Figure 3 shows the average extinction derived based on the median value of the fitting parameter α for the sight lines with pcFlag = 'UVI'.
The average extinction curve in M33 shows similarity to the extinction curve in the diffuse region of the MW and the average LMC extinction curve, but with a slightly weaker 2175Å bump and a slightly steeper rise in the UV bands. The average extinction law in M33 derived in this work can be applied to the general extinction correction in M33, and those toward individual sight lines can help with high-precision extinction correction (see Section 4.5 for details).
The dust size distributions toward the selected tracers with pcFlag = 'UVI' are plotted in Figure   4 with gray solid lines. As we know, grains with sizes that are comparable to the wavelength absorb and scatter light most effectively (2πa/λ ≈ 1, where a is the spherical radius of the grain, Li 2009).
Interstellar dust has long been considered to be "submicron-sized" (≈ 0.1 µm) since dust extinction was first confirmed by Trumpler (1930), because grain models that reproduce the observed extinction should have extinction in the visible bands (λ ≈ 0.55 µm) dominated by grains with a ≈ 0.1 µm (Draine 2011). However, it is now well recognized that the dust size actually spans a wide range from subnanometers to micrometers . In addition, the strong rise to λ ≈ 0.1 µm on the extinction curves requires a large abundance of grain with 2πa/λ 1; thus, interstellar dust must include a large population of grains with a 0.015 µm. Based on equation (8) in Nozawa (2016), we derive that the average dust size a toward individual sight lines with pcFlag = 'UVI' is in the range of 5.78-9.64 nm. The median value of the average dust size is a ≈ 7.54 nm, which is smaller than those of the MW (a MW ≈ 8.36 nm) and M31 (a M31 ≈ 8.41 nm). Table 2. An extracted list of the selected tracers in the sample with spectral type, fitting parameters (α, log T and information about observed photometry adopted in this work a .
LGGS ID LGGS SpT The average dust radius a is calculated from equation (8) in Nozawa (2016).
e This is the flag for the coverage of passbands adopted in the calculation. U = UV bands (here, this refers to the passbands bluer than the U band). V = Visual band (here, this refers to the passbands from the U to y bands). I = IR band (here, this refers to the passbands redder than the y band).  Note-a The superscript and the subscript in the table indicate the derived upper limit value and lower value of the derived parameters extracted from Table 2 for the selected tracers.
b The results of the sight lines with pcFlag = 'UVI' are adopted to derive the general extinction law in

M33.
c For tracers with near-IR (UV) data, the calculation is repeated without taking near-IR (UV) data into consideration, and reliable results are listed with * for comparison. e The Levenberg-Marquardt method is adopted to fit the model extinction curves to the F19 extinction curve with R V = 3.1. There is only one parameter (α) in our model extinction curves, and a grid ranging from 0.50 to 7.00 with a step of 0.01 is taken.

Comparison with the Local Galaxies
We compare the average extinction curve in M33 with those in the MW and other local group galaxies (SMC, LMC, M31) in Figure 3. The average LMC extinction law (Nandy & Morgan 1978; Fitzpatrick 1986; Gordon et al. 2003) resembles that in the MW, while most of the extinction curves in the SMC bar region (Prevot et al. 1984;Gordon et al. 2003) display a nearly linear rise with λ −1 and an absent 2175Å, similar to those in the starburst galaxies (Calzetti et al. 1994). The average extinction curve in M31 derived in Paper I (yellow dashed lines in Figure 3) shows similarity to that of the MW but rises less steeply in the far-UV bands. The average extinction curve in M33 is similar to that of M31 in shape but with a slightly larger slope.
The average dust extinction curve derived in this work is also compared with the attenuation curves derived in Gordon et al. (1999) and Hagen (2017) in Figure 3. As illustrated in Section 1, attenuation curves include both extinction and the assumed geometry of dust and stars, so attenuation is aimed at the effect of dust on an area instead of an individual sight line. Gordon et al. (1999) adopted radiative transfer modeling from UV to NIR of the M33 nucleus, which is an ideal interstellar environment of starburst, and found an MW-like attenuation curve with a strong 2175Å bump (solid cyan line in Figure 3). Hagen (2017) modeled the SEDs for 1170 large pixels in M33 from FUV to NIR and derived a steep median attenuation curve with a weaker 2175Å bump (solid green line in Figure 3).
The average extinction curve in M33 derived in this work presents a similar slope to Gordon et al. (1999) but with a weaker 2175Å bump as the median one in Hagen (2017). In this work, we map the derived A V of the selected tracers in Figure 5 and find that the median value of A V is ≈ 0.43 mag, which is slightly smaller than the median A V (≈ 0.53 mag) derived in (Hagen 2017) and larger than the mean amount of dust extinction (A V ≈ 0.25 mag) measured in Verley et al. (2009).
The discrepancy may be due to the different scales of dust and the different stellar models (Conroy 2013). In addition, we eliminate the results with E(B − V ) < 0.06 mag, as mentioned in Section 3.3, because slightly reddened stars may lead to larger errors. As a result, tracers with small A V values are excluded, increasing the median value of A V . (2022)

The 2175Å bump
The 2175Å bump, which is the broad excess in the extinction curve at a rest wavelength λ ≈ 2175 A, is the strongest signature of dust in the interstellar medium (Kashino et al. 2021). It has been considered as a unique probe of the nature of dust in galaxies since it was discovered by Stecher (1965). The 2175Å bump is obvious in the extinction curves toward the individual sight lines of the MW (e.g., Fitzpatrick & Massa 1986, 1990; F19), LMC (Nandy & Morgan 1978;Fitzpatrick 1986;Gordon et al. 2003) and M31 (Dong et al. 2014;Clayton et al. 2015; Paper I), while it was almost absent in the SMC (Prevot et al. 1984;Gordon et al. 2003). On galaxy scales, it was found that there is no significant 2175Å bump in the attenuation curves of nearby starburst galaxies (Calzetti et al. 1994;Gordon et al. 1997;Calzetti et al. 2000) and Lyman break galaxies at high redshifts (z > 2, Vijh et al. 2003). As a result, the attenuation curves with no 2175Å bump in Calzetti et al. (2000) are commonly adopted for both local and distant star-forming galaxies. However, the 2175Å bump has been detected and even measured for star-forming galaixes by many recent works (e.g. Noll et al. 2007Noll et al. , 2009Buat et al. 2011Buat et al. , 2012Scoville et al. 2015;Battisti et al. 2017;Salim et al. 2018;Battisti et al. 2020;Shivaei et al. 2020).
As to M33, although it was one of the star-forming galaxies in Calzetti et al. (1994) with no significant 2175Å bump in the attenuation curve, recent studies (Gordon et al. 1999;Hagen 2017) indicated that there exists a 2175Å bump in the attenuation curve. Since graphite is one of the possible candidates of the carriers of the 2175Å bump (Stecher & Donn 1965), we adopt the silicategraphite dust model in this work to derive the overall extinction curves from UV to near-IR toward individual sight lines in M33. In order to find out whether the 2175Å bump really exists in the extinction curves of M33, we also adopt the model extinction curves without a 2175Å bump derived from the silicate dust model (f cs = 0, no carbonaceous grains) to repeat the calculation for the tracers with ultraviolet data. By comparing the median values of χ 2 /d.o.f. 5 derived from both dust models for each tracer, it is found that the extinction curves derived from the silicate-graphite dust model can generally recover the observed SEDs better than the silicate dust model. We therefore suggest that there exists a 2175Å bump in the extinction curves of M33. The fine structure of the where N data is the number of the observed photometric points adopted in the calculation, N para is the number of adjustable parameters (see Section 3.1 for details), f observed is the observed flux of the photometric point, f model is the model flux of the photometric point and σ is the difference between logarithm of extreme and logarithm of f observed .
UV extinction curves can be analyzed and more comprehensive results can be expected, if the UV data is adequate in the future.

Influence of IR and UV Photometry
Photometric data that cover a wider range of passbands will constrain the observed SED better and bring more reliable results. Because of the observation limit, a number of tracers lack photometry in the U V W 2, U V M 2, U V W 1 and PHATTER bands. Meanwhile, UKIRT data are also not applied to all tracers in the extinction sample, as mentioned in Section 2. It is thus necessary to determine whether the lack of photometry in the UV and near-IR bands affects the derived extinction law.
Lines 6 and 8 in Table 3 summarize the results of the selected tracers with near-IR data and with UV data, respectively. We repeat the calculation for these two groups of tracers but ignore the near-IR data or UV data and list the number of tracers with reliable results as well as the derived results in the seventh and ninth lines of Table 3, respectively. As shown in Table 3, the number of tracers with reliable results is significantly reduced when UV data or near-IR data are not adopted in the calculation, indicating that photometric data in wider bands bring more reliable results.
On the other hand, we compare the reliable results derived with near-IR (UV) data ignored for the tracers with near-IR (UV) data [Tracers in Line 6 (8) of Table 3] and the results extracted from Table 2 for the same tracers in Figure 6. As Figure 6 shows, the lack of UV or near-IR data has little impact on A V and E(B − V ) in this work. However, the dust size parameter α and the average dust size a for most of the individual tracers are influenced by the coverage of the adopted passbands, indicating that UV and near-IR data are important to constrain the dust model.
To illustrate the reliability of the derived results for individual sight lines, as mentioned in Section 2, we introduce pcFlag in Table 2 to show the coverage of passbands adopted in the calculation.
The results of sight lines with pcFlag = 'UVI' are the most reliable, while those with pcFlag = 'V' are the least reliable. We anticipate that the results could be more comprehensive if the observed data in multiple bands are adequate, especially in UV bands, because UV data can provide a strong constraint on the extinction model. The coming 2 m-aperture Survey Space Telescope (also known as the China Space Station Telescope, CSST) will image approximately 17500 square degrees of the sky in the N U V , u, g, r, i, z and y bands (Zhan 2021) and will provide us with abundant data to explore the dust extinction law in M33 and other nearby star-resolved galaxies.

Prediction of Multiband Extinction
As mentioned in Section 4.1, the dust extinction curves derived in this work can help with the extinction correction in M33. Based on the average extinction curve of M33 in this work, multiband extinction values from UV to near-IR are predicted, which are shown in Table 4. High-precision extinction correction should refer to the extinction law of individual tracers derived in this work. ascension and declination listed in the first four columns of Table 5 for each tracer are obtained from the LGGS catalog (Massey et al. 2016). The α, a and E(B − V ) in Table 5 are extracted from Table   2. The column named pcFlag presents the passband coverage, as shown in Table 2. The following columns show the extinction values in multiple bands for each tracer. Although the results for some tracers are not affected by the lack of UV or near-IR data, we recommend multiband extinction toward individual sight lines with pcFlag = 'UVI'.
When applying the method and results in this work, there are three aspects that need to be noted. First, the extinction curves derived in this work are applicable from UV to near-IR (≈ 3 µm) bands. All dust models for the diffuse ISM predict that an extinction curve steeply declines with λ at 1 µm < λ < 7 µm and increases at λ > 7 µm because of the 9.7 µm silicate absorption feature (Mathis et al. 1977;Kim et al. 1994;Weingartner & Draine 2001;Li et al. 2015). However, many recent observations suggest that the extinction law in the mid-IR band (3 µm < λ < 8 µm) appears to be universally flat or gray in various interstellar environments (Lutz et al. 1996;Lutz 1999;Indebetouw et al. 2005;Flaherty et al. 2007;Gao et al. 2009;Nishiyama et al. 2009;Fritz et al. 2011;Wang et al. 2013). Although the µm-sized grain can be adopted to model the flat mid-IR extinction curve , it will make the model more complex and, thus, reduce the universality of the method. As a result, the derived extinction curves cannot be applied to mid-IR bands at present. repeating the calculation process but without taking the available near-IR (UV) data into consideration.
The deviation between the red (blue) dots and y = x presents the influence of the lack of near-IR (UV) data.
Moreover, the classic silicate-graphite dust model adopted in this work may not be applied to analyze the fine structure of the extinction curves in UV bands, although it can be adopted to derive the overall reliable extinction curves in M33 from UV to near-IR bands. As illustrated in Section 4.3, the 2175Å bump is known to be an important feature of the extinction curves in UV bands, which was first discovered by Stecher (1965). Since Stecher & Donn (1965) pointed out that small graphite particles would produce absorption very similar to this observed feature, some form of graphitic carbon has been an attractive candidate because the π → π * transition in graphite is responsible for the absorption feature at ∼ 2175Å (Draine 2003). However, this graphite hypothesis does not appear to explain the fact that the full width at half maxima (FWHM) of 2175Å varies with the interstellar environment while holding the central wavelength λ 0 nearly constant (Draine & Malhotra 1993). Currently, a polycyclic aromatic hydrocarbon (PAH) mixtures are carrier candidates for the 2175Å bump (Joblin et al. 1992;Li & Draine 2001;Xiang et al. 2011;Steglich et al. 2011;Mishra & Li 2015, 2017 because PAH molecules generally have strong π → π * absorption in the 2000 -2500 A region with variation in FWHM and small variation in λ 0 . As a result, we consider adding PAHs to the dust model in future work to analyze the fine structure of UV extinction curves and obtain a more detailed understanding of the dust properties in M33 and other nearby galaxies.
Finally, as shown in Figure 5, the size of the extinction sample adopted in this work is not adequate to cover the entire region of M33. Thus, it can only provide us with a low-resolution extinction map to help with a rough extinction correction for certain regions in M33. A major science project named CSST mentioned in Section 4.4 will provide us with larger extinction samples and adequate data in multiple bands. We can expect further exploration of the dust properties and extinction law in M33 and in other nearby star-resolved galaxies, as well as the development of higher-precision extinction corrections in the near future.
is the filter transmission function and V g(λ) is the Vega spectrum. b A λ is the average extinction in M33 based on the median value of A V . LGGS ID LGGS SpT Note-a This is an extracted table for the same tracers as listed in Table 2 with a portion of the multiband extinction values. The entire table is available in machine-readable form. b Flag for the coverage of passbands adopted in the calculation (the same as Table 2) 4.6. Application in Other Nearby Star-Resolved Galaxies The method adopted in this work and Paper I is extended to other nearby star-resolved galaxies.
The LGGS provides U BV RI plus the interference-image photometry of luminous stars in seven systems currently forming massive stars (IC 10, NGC 6822, WLM, Sextans A and B, Pegasus and Phoenix, Massey et al. 2007) in addition to the spiral galaxies M31 and M33 (Massey et al. 2006(Massey et al. , 2016. We can isolate O-type and B-type supergiants in NGC 6822 and WLM from the LGGS catalog (Massey et al. 2007), which are selected as the extinction tracers. The results for all the selected tracers are listed in Table 6, and the derived extinction curves in NGC 6822 and WLM are plotted in Figure 7   as well as the average extinction curve in M33 derived in this work (cyan dashed line). Table 6. A list of the selected tracers in NGC 6822 and WLM with spectral type, fitting parameters (α, log T eff , log g, A V ), derived a, E(B − V ),

R
V and information about the observed photometry adopted in this work a .

Galaxy
LGGS ID LGGS SpT These lines show the median values of the derived parameters for selected tracers in NGC 6822 and WLM.

SUMMARY
A sample of bright O-type and B-type supergiants from the LGGS catalog (Massey et al. 2016) are chosen as extinction tracers to derive the dust extinction curves in M33. This is the first study focused on the dust extinction curves toward individual sight lines in M33 rather than the dust attenuation curves in Gordon et al. (1999); Hagen (2017). The previous studies have been improved, and the main results of this work are as follows: 1. The extinction curves in M33 derived in this work cover a wide range of shapes, from curves with an obvious 2175Å bump (like the extinction curves with R V ≈ 2) to relatively flat curves with R V ≈ 6, implying the complexity of the interstellar environment and the inhomogeneous distribution of interstellar dust in M33. The derived parameter α in the dust size distribution ranges from ≈ 2.6 − 5.9, while the dust size ranges from ≈ 5.78 − 9.64 nm.
2. The average extinction curve in M33 (R V ≈ 3.39) is similar to the MW extinction curve with R V = 3.1 but with a slightly weaker 2175Å bump and a slightly steeper far-UV rise. The average dust size distribution in M33 is dn/da ∼ a −3.45 exp(−a/0.25), and the median value of the average dust size is a ≈ 7.54 nm, which is smaller than that of the MW (a MW ≈ 8.45 nm).
3. The derived A V in M33 is up to 2 mag with a median value of ≈ 0.43 mag, which is smaller than the median value (A V ≈ 0.53 mag) derived in Hagen (2017) and larger than the mean amount (A V ≈ 0.25 mag) measured by Verley et al. (2009). 4. The method adopted in this work and Paper I that combines the stellar model atmospheres and the dust models to calculate extinction curves and analyze dust properties toward individual sight lines is extended to the star-resolved galaxies NGC 6822 and WLM, but we can only derive the extinction curves toward a few individual sight lines. More observations are needed to gain a better understanding of the extinction law and dust properties in nearby star-resolved galaxies.