The Chocolate Chip Cookie Model: Dust Geometry of Milky-Way like Disk Galaxies

We present a new two-component dust geometry model, the \textit{Chocolate Chip Cookie} model, where the clumpy nebular regions are embedded in a diffuse stellar/ISM disk, like chocolate chips in a cookie. By approximating the binomial distribution of the clumpy nebular regions with a continuous Gaussian distribution and omitting the dust scattering effect, our model solves the dust attenuation process for both the emission lines and stellar continua via analytical approaches. Our Chocolate Chip Cookie model successfully fits the inclination dependence of both the effective dust reddening of the stellar components derived from stellar population synthesis and that of the emission lines characterized by the Balmer decrement for a large sample of Milky-Way like disk galaxies selected from the main galaxy sample of the Sloan Digital Sky Survey (SDSS). Our model shows that the clumpy nebular disk is about 0.55 times thinner and 1.6 times larger than the stellar disk for MW-like galaxies, whereas each clumpy region has a typical optical depth $\tau_{\rm{cl,V}} \sim 0.5$ in $V$ band. After considering the aperture effect, our model prediction on the inclination dependence of dust attenuation is also consistent with observations. Not only that, in our model, the dust attenuation curve of the stellar population naturally depends on inclination and its median case is consistent with the classical Calzetti law. Since the modelling constraints are from the optical wavelengths, our model is unaffected by the optically thick dust component, which however could bias the model's prediction of the infrared emissions.


INTRODUCTION
Dust accounts for only a tiny fraction (∼ 0.1%) of the baryonic mass in star-forming galaxies (SFG) but plays a crucial role in many aspects of galaxy evolution and observational properties.Dust grains are formed in the ejected material of supernova (Ferrarotti & Gail 2006) and the stellar wind of low-mass stars at the end of their lives (Indebetouw et al. 2014;Dwek & Cherchneff 2011).They grow, coagulate, and are destructed in the interstellar medium (ISM).
Dust particles absorb and scatter ultraviolet (UV) and optical photons and re-emit the energy in the infrared (IR) wavelengths.The absorption and scattering of photons lead to a reduction of the emitted flux, known as attenuation, and also cause the reddening of the spectral energy distributions (SED) since the photons with shorter wavelengths are generally more attenuated than those with longer wavelengths.Because of the attenuation and reddening, dust attenuation is a critical effect to be quantified in deriving the physical properties, such as stellar mass and star formation rate (SFR) of galaxies (Kennicutt 1998;Popescu et al. 2011;Gadotti et al. 2010;Pastrav et al. 2013a,b).
In the optical spectral analysis, the reddening features of nebular emission lines and stellar continua are the two most commonly used probes of the dust attenuation of a galaxy.At low redshift, Balmer decrement is commonly used to calculate the reddening of the nebular emission, E g (B − V).For stellar continuum, the dust reddening E s (B − V) can be reliably deduced by using stellar population synthesis (SPS) models.
Many studies have shown that the nebular emission in galaxies is more attenuated than the stellar continua (Calzetti et al. 1994;Mayya & Prabhu 1996;Charlot & Fall 2000), which is usually explained by a two-component dust model (e.g.Calzetti et al. 2000;Charlot & Fall 2000).In this model, the dust distribution in the galaxy has two components: a diffuse ISM component and a clumpy birth-cloud component.Thus, the nebular emission originated from the star-forming regions suffers from the dust attenuation of both the envelope of birth clouds and the diffuse ISM, while the stellar emission, especially that of old stellar populations, most of which is radiated outside of the birth clouds, is extincted only by the diffuse ISM component.Based on the observations of a sample of local starburst galaxies, Calzetti (1997) proposed a typical ratio of the dust reddening for the stars to that for the nebular lines, which has been widely adopted in many later studies (Madau & Dickinson 2014;Peng et al. 2010;Daddi et al. 2007).However, many recent studies have shown that this stellarto-nebular dust attenuation ratio f varies systematically with the physical properties of galaxies.For example, Wild et al. (2011) find that f decreases with decreasing specific star formation rate (sSFR) for star-forming galaxies.Koyama et al. (2015) also show that more massive galaxies tend to have higher extra attenuation towards nebular regions, i.e., smaller f values (see also Zahid et al. 2017).These trends are consistent with the fact that the contribution of old stellar populations, which are less obscured than the young stellar populations that are more closely associated with nebular emission, increases in lower sSFR (or more massive) galaxies.
The dust attenuation of a galaxy is determined by the contents of its emission sources and dust particles as well as their geometry configuration.As a result, the dust attenuation features of disk galaxies have a strong dependence on their viewing angles.Many studies have examined the inclination dependence of the observed properties of disk galaxies (e.g.Shao et al. 2007;Maller et al. 2009;Masters et al. 2010;Yip et al. 2010;Li et al. 2021;Yuan et al. 2021), and reported that the stellar attenuation is proportional to the logarithm of the observed disk axis ratio (A λ ∝ −γ log(b/a)) with γ ∼ 1 (Yip et al. 2010;Chevallard et al. 2013;Battisti et al. 2017), while the nebular attenuation indicated by the Balmer decrements only shows a slight increase trend at low axis ratio.Moreover, Wild et al. (2011) shows that more inclined disks have smaller f values and flatter attenuation curves.On the other hand, these observational trends with disk inclination provide strong constraints on the geometry configuration of disk galaxies.For example, Yip et al. (2010) conclude that the inclination dependence of the stellar continuum attenuation features deviate significantly from the simple dust screen model, while the slab model (stars and dust mixed uniformly) and the sandwich model (a layer of dust+stars mixture sandwiched in-between two layers of stars) also cannot account for the observed features appropriately.
Many studies have used radiative transfer (RT) models to constrain the geometry parameters of the dust and stellar components for individual nearby disk galaxies(e.g., Xilouris et al. 1999;Misiriotis et al. 2001;Bianchi 2007;De Geyter et al. 2014).Popescu et al. (2000) proposed a multicomponent dust distribution model for edge-on disk galaxy NGC 891, which is composed of a bulge, a thick disk associated with diffuse dust and old stellar populations, and a thin disk associated with clumpy dust and newly-formed stars.Tuffs et al. (2004) (hereafter T04) extended the multicomponent of Popescu et al. (2000) and further refined the prescription for attenuation of clumpy dust component associated with star-forming regions in the thin disk.For the UV photons emitted from star forming regions, the T04 model assumes that they are either completely absorbed or are able to escape.The fraction of the UV photons being absorbed is then quantified by a free model parameter, clumpiness F. With this setting, the T04 model can further make a selfconsistent determination of the attenuation of the nebular lines.However, it is worth pointing out that the clumpiness factor F in T04 is not a geometry parameter, where the starforming regions are continuously distributed.
The T04 model has been widely used to study various dust attenuation properties of disk galaxies (Driver et al. 2008;Leslie et al. 2018b,a).For example, by properly setting few model parameters (e.g.F and the disk face-on dust opacity), the T04 model reproduce the observed inclination dependence of the broad band dust attenuation for both low and intermediate redshift galaxies successfully (Driver et al. 2007;Masters et al. 2010).However, a recent study of van der Giessen et al. (2022) (hereafter G22) finds that a completely optically thick dust component, as assumed in the T04 model, can not reproduce the observed inclination-dependence of the Balmer decrement.An extra component of optically thin dust within the birth clouds is required.In addition, it is worth noting that the attenuation properties derived from T04 model are integrated for global galaxy, while the observed Balmer decrement of SDSS is typically measured from fiber spectroscopy.This spatial mismatch will introduce extra bias into the model that tries to fit these two different attenuation features simultaneously (Chevallard et al. 2013).
In this study, we aim to build a new geometry model of disk galaxies to account for the inclination dependence of the nebular and continuum dust attenuation simultaneously.To avoid of the mismatch of the scales, we derive both of the nebular and continuum dust attenuation from the fiber spectroscopy.Correspondingly, our modeling of the dust attenuation effects will be mainly along the line of sight to galaxy center rather than make a global estimation.In particular, we will provide a refined modelling of the dust attenuation of the star-forming regions, where its self-extinction and clumpy distribution (as outlined by Charlot & Fall (2000)) are both quantified by model parameters.As a starting point of this new model, we will only probe the dust attenuation and reddening features in optical wavelengths, whereas the dust emission in IR bands will be left for subsequent studies.
The outline of this paper is as follows.In Section 2, we introduce the local star-forming galaxy sample and show the inclination dependence of two different dust reddening features.In Section 3, we model the inclination dependence features with two simple geometric models.In Section 4, we present the Chocolate Chip Cookie model, a novel two dust component model, and use it to fit the observed inclination dependence of two dust reddening features simultaneously.In Section 5, we discuss the fiber aperture effect in our sample galaxies and use it to explore the inclination dependence of the emission-line fluxes.In Section 6, we make further discussions on the best fits of our Chocolate Chip Cookie model.Finally, we summarize our main results in Section 7.

Milky-Way Like Galaxy Sample
Our galaxy sample is drawn from the Main Galaxy Sample (MGS) of the Sloan Digital Sky Survey Data Release 7 (SDSS DR7) 1 .We take the measurements of total stellar mass, fiber magnitude, Hα, Hβ, [OIII]λ5007 and [NII]λ6584 emission line fluxes from the MPA-JHU data release 2 (Kauffmann et al. 2003a;Tremonti et al. 2004;Brinchmann et al. 2004).To study the dust attenuation properties of both the stellar continua and nebular emission lines, we select starforming galaxies (SFGs) with emission line features and then use the classical BPT diagram (Baldwin, Phillips, & Terlevich 1981) to remove galaxies with signs of active galactic nuclei (AGNs).We adopt the criteria of Kauffmann et al. (2003b) and require signal-to-noise ratio (S/N) larger than 3 for the four emission lines used.This selection results in 260,856 SFGs.We take the stellar mass estimation from MPA-JHU and obtain the inclination angle θ of each SFG from Simard et al. (2011).According to Simard et al. (2011), θ is estimated from morphological fitting by applying a disk+bulge model and therefore gives a better approximation of the real inclination of disk galaxies than the simple axis ratio (b/a) parameter, especially for highly inclined disks.To avoid possible observational biases in the physical properties of our sample galaxies (Driver et al. 2007) , we further require sample galaxies at z < 0.1, where the PSF effects on the modelling of disk inclination could be minimized.To illustrate this effect, we show the cos θ distributions of all the SFGs in MGS and the selected SFGs at z < 0.1 as the dotted and solid histograms in the top panel of Figure 1 respectively .As can be seen, compared with all SFGs, SFGs at z < 0.1 show a much flatter cosθ distribution, indicating a better and unbiased measurement of disk inclination.In the bottom panel of Figure 1, we show the number density of the SFGs at z < 0.1 as function of their stellar mass and disk inclination.We see little systematic bias in the stellar mass of these SFGs at different inclinations.
In this study, our main goal is to construct a modeling framework for the geometry of different dust components in SFGs using a statistical approach.The Milky-Way like (MW-like) galaxies have stable and well-formed disks and the largest sample size in MGS, making them most suitable for this study.Therefore, we further select the SFGs with the stellar mass in the range 10 10.2 − 10 10.6 M .These selection criterics result our final MW-like galaxy sample of 33,273 galaxies.We leave the exploration of the geometry of disk galaxies with other stellar masses to an upcoming study.

Balmer decrement of emission lines E g
To isolate the effect of attenuation curves, we take the Balmer decrement, which compares the ratio of the intensities of Hα and Hβ emission lines ( f Hα / f Hβ ) obs with its intrinsic value ( f Hα / f Hβ ) int , to represent the nebular attenuation, where the intrinsic Balmer decrement ( f Hα / f Hβ ) int is set to the typical Case B recombination value 2.86.Because the conversion from E g (Hα − Hβ) to classical E g (B − V) relies on the shape of the attenuation curve3 , we use E g (Hα − Hβ) instead of classical E g (B − V).For simplicity, we denote the nebular color excess E g (Hα − Hβ) as E g hereafter.We note that Equation 2 could result in negative color excess for some galaxies, which are likely to be caused by the measurement errors of the nebular emission line fluxes.For these galaxies, we set E g = 0.

Reddening of stellar population E s
We derive the dust attenuation of the stellar population by fitting the spectra of galaxies with SPS.In SPS, the observed spectrum of a galaxy is typically fit through where the component in bracket is the sum of the fraction of each single stellar population (SSP), and τ λ is the effective optical depth at wavelength λ that parameterizes the global dust attenuation effect.Also, by assuming an attenuation curve for the stellar population, k s (λ), the optical depth τ λ is parameterized by where E s (B − V) is the global dust reddening of the stellar continuum to be fit.We use E s (B − V) because it is less influenced by the attenuation curve shape than attenuation A V (see more discussions in Section 6.3).We use the full-spectrum stellar population fitting code STARLIGHT (Cid Fernandes et al. 2005) to fit the stellar continuum (i.e. with emission lines masked) of each galaxy and derive its E s (B − V) value.We adopt the BC03 stellar population (Bruzual & Charlot 2003) for SSP i (λ) and the standard Calzetti attenuation law (Calzetti et al. 2000) for k s (λ).Hereafter, we refer to the derived stellar color excess There is a known dust-age-metallicity degeneracy effect in SPS fitting.As shown by Li et al. (2020), among dustage-metallicity, the dust attenuation feature is the one that can be most accurately and unbiased recovered from SPS fitting.Moreover, since E s (B − V) is a strong function of disk inclination, if its measurement were biased by the dust-agemetallicity degeneracy effect, we would expect an inclination (dust attenuation) dependence of the average stellar population age and metallicity.We have checked this effect and find no inclination dependence (see appendix A for detail).

E g and E s
We show the joint distribution of E g and E s for our sample galaxies in the top panel of Figure 2, where we have removed 23 outliers with E g or E s values deviate more than 5σ from their median values.The final sample contains 33,250 galaxies.As can be seen, E g spans a wide range from 0 to ∼ 1 mag with a median value of ∼ 0.45 mag, while E s is mostly distributed in the range 0-0.4 mag and with a median value 0.17 mag.
In general, E g and E s show a monotonic correlation and with a mean ratio of E s /E g ∼ 0.4.Considering the conversion factor from E g to E g (B − V) (∼ 0.9 for the Calzetti attenuation curve), this mean ratio is consistent with the conical value f ∼ 0.44 of Calzetti et al. (2000).

Inclination Dependence of Dust Attenuation
We show the inclination dependence of E s and E g for our MW-like galaxies by color-coding each galaxy with their disk inclination θ in the top panel of Figure 2. As expected, there is a general trend that the heavily attenuated galaxies are mostly these galaxies with large inclination angle θ.However, at given θ, there are large scatters on both E g and E s , which reflects the variation of the intrinsic dust attenuation among different galaxies.
To better characterize the inclination dependence of these two types of dust reddening, we divide the sample galaxies into 90 θ bins and require that each θ bin includes the same number of galaxies.We show the median values (square dots connected with a solid line) of E s and E g for each θ bin in the bottom left and bottom right panels of Figure 2, where the 16 and 84 percentiles of the E g and E s distribution are indicated by the gray solid lines.Because of the large number of the galaxies in each θ bin (n ∼ 300), the uncertainty of the median values of E g and E s are both smaller than 0.01.As shown in Figure 2, at inclination θ < 75 • , both E g and E s increase monotonically with the disk inclination θ, while E g shows a sharper increase trend.For highly inclined disk galaxies (θ > 75 • ), the trend of E s becomes flat, whereas the trend of E g is reversed, i.e.E g decreases with increasing θ.
The monotonic increase trend of E s as a function of θ has been reported by many studies, and can be reasonably explained by either simple parametric models or dedicated RT models (e.g.Shao et al. 2007;Maller et al. 2009;Masters et al. 2010;Yip et al. 2010).Different trends of E g and E s as functions of disk inclination also have been found in observations (e.g.Yip et al. 2010;Chevallard et al. 2013;Battisti et al. 2017).In our study, because of the equal-size binning of θ and the large sample size, for the first time, we reveal a complicate inclination dependence of E g , especially at high θ values.

MODELLING OF INCLINATION DEPENDENCE OF E G AND E S : SIMPLE MODELS
From this section, we aim to build geometric models for disk galaxies so as to model the observed inclination dependence of E s and E g shown in the bottom-left and bottomright panels of Figure 2 .We start with the two most frequently used simple models, the uniform mixture model and the screen model, to model the observed E s (Section 3.1) and E g (Section 3.2), respectively.Before that, we introduce the common parts of the two different dust attenuation models.
We assume that our MW-like galaxies have the same geometry configuration so that the inclination dependence of attenuation being observed is only caused by the angle of view.In our modelling, we assume that the main components of disk galaxies follow a double exponential model in geometry, where R comp , h comp are exponential disk scale-length and scale-height, respectively.Here, depending on where the equation is used, the component could either be referred to emission source (e.g.stars) or dust component.Moreover, since the SDSS fiber spectra only cover the central part of the observed galaxies, we also consider only the attenuation of the line of sight through the center of the model galaxy.
For this line of sight, the integration of a double exponential component has inclination dependence For convenience, we set F comp = 1 when galaxies are viewed as edge-on (θ = 90 • ).
In addition, we make a convention for the subscripts of the symbols to be used.The order of the subscripts is 'wavelength', 'component' and 'specific'.The wavelength term is typically written as 'λ' .If not explicitly specified, the default wavelength is at V-band.The 'component' term includes 'stellar disk','nebular disk', and 'clumps', which are denoted by's', 'g ', 'cl' respectively.The 'specific' term includes'central', 'total', 'un-attenuated', denoted by '0', 'A', and '*', respectively.

uniform mixture Model for E s
On galactic scale, the dust attenuation of diffuse ISM dust on stellar emission can be simplified and parameterized by a uniform mixture model.We show a schematic map for such a simple mixture model in the top panel of Figure 3.In specific of a disk galaxy, we assume that the density profiles of stellar population and diffuse ISM dust follow the same exponential disk model (Equation 5), where R s and h s are the scale-length and scale-height of the stellar disk, ρ s,0 is the central effective density of dust particles or stellar population, depending on the usage of this Equation.
With this exponential disk modelling, for a disk galaxy with inclination θ, the total line of sight optical depth along the galactic center direction is where κ λ is the dust extinction coefficient (absorption crosssection) at given wavelength λ, and F s (θ) is the inclination  dependence term shown in Equation 6.We parameterize the wavelength dependence of dust extinction coefficient, namely the dust extinction curve κ(λ) with a simple powerlaw, which is normalized by V band extinction coefficient κ V at its effective wavelength 5500Å.We set the power-law index β to be 1.32, so that it has R V = A V E(B−V) = 3.1, a typical value for the diffuse ISM dust in the Milky-Way (Weingartner & Draine 2001;Fitzpatrick et al. 2019;Li et al. 2017).
Moreover, during the calculation of τ λ,s,A (θ), κ V is coupled with the central dust density parameter ρ s,0 .Therefore, for simplicity, we define a new parameter α s,0 ≡ ρ s,0 κ V , which directly represents the central dust absorption density.With this new parameter, we have For the uniform mixture model, the total emitted intensity I is linked to the total unextincted intensity I * through where τ λ,s,A is the total optical depth along the given line of sight.The dust reddening then follows ) .
(12) As shown by Equations 11 and 12, for the uniform mixture model, when the optical depth is thick (τ >> 1), there is a simple correlation of emitted intensity with τ (I ≈ I * /τ) and a nearly constant color excess The optical depth τ of the uniform mixture model (Equation 10) is composed of two terms.One is that describing the global amount of dust attenuation, which could be parameterized by the V band dust optical depth to the galactic center in edge-on case, τ s,A (90 • ) ≡ 2α s,0 R s .The other is the inclination dependence term F s (θ) (Equation 6), which has the only parameter h s /R s , i.e., the disk scale-height to scale-length ratio.Since the disk scale-length R s is included in both terms, for simplicity, we set R s = 1.Then, the two fitting model parameters become α s,0 and h s , which are in units of R −1 s and R s , respectively.We fit the observed inclination dependence of E s with these two free model parameters.The best fitting of the E s − θ relation is shown as the red solid line in the top panel of Figure 3 and the two best model parameters are listed in Table 1.
Figure 3 (top panel) shows that this uniform mixture model fits the observed E s − θ relation quite well.In our best fitting, the disk component has a height-to-radius ratio h s /R s = 0.27.The best estimate of α s,0 is 1.9, which means that the global V band optical depth along the galaxy center direction is varying from ∼ 1 (face-on) to ∼ 4 (edge-on), corresponding to the effective dust attenuation from A V ∼ 0.50 mag (face-on) to 1.41 mag (edge-on), which is broadly consistent with the observational values (0.5−1.5 mag, Bianchi 2007;De Geyter et al. 2014;Casasola et al. 2017).

Screen Model for E g
The nebular emission of star-forming (HII) regions is known to be distributed on a much thinner disk than the stellar component (e.g.Anderson et al. 2019).Considering that the SDSS fiber spectra only target the central regions of galaxies, these HII regions might be viewed as a pointlike source in galaxy center if they do not overlap each other along the line of sight.In this simplified model, the nebular emission lines emitted from the central star-forming regions are then only subject to the extinction of outer dust layers along the line of sight, which is consistent with a screen model.
For the dust screen model, the observed intensity of an object with intrinsic intensity I 0 is determined by the optical depth τ of the dust layer along the line of sight, I = I 0 e −τ .We show a schematic map for this simple dust screen model in the bottom panel of Figure 3.The line of sight dust consists of two parts.One is the dust component in diffuse ISM and the other is the dust shell of nebular region itself.For the former component, we take the best model estimate of Section 3.1, i.e. an exponential diffuse ISM dust disk layer with geometric configuration h s /R s = 0.27 and the central dust absorption density α s,0 = 1.9.Since the nebular emission region is assumed to be located at galaxy center, the inclination dependence of its line of sight optical depth from diffuse ISM dust is therefore half of the τ λ,s,A (θ) that parameterized by Equation 10.For the dust component of the nebular emission region itself, we assume that the dust layer is distributed in a spherical shell.In this case, the dust optical depth of the shell has no inclination dependence and therefore can be parameterized by a constant τ λ,cl .Putting these two dust components together, for disk galaxies with inclination θ, the total line of sight optical depth of the modelled nebular emission line region in galaxy center is We assume that the extinction curve of the dust particles in nebular regions is the same as that of diffuse ISM dust (Equation 9) and then parameterize τ λ,cl with the V band optical depth τ cl .In this simple dust screen model, the observed color excess of emission lines is (14) With above model, we fit the observed inclination dependence of E g shown in the bottom right panel of Figure 2 with the only parameter τ cl .We obtain the best fit with τ cl = 0.54, which is shown as the solid line in the bottom panel of Figure 3 and listed in Table 1.As can be seen, by including the dust component in nebular regions, the dust screen model of Equation 13 generally reproduces the observed E g − θ relation for low inclination disks (θ < 70 • ).However, for these high inclined disks (θ > 70 • ), the observed E g − θ relation shows a saturation effect, which can not be accounted by our simple dust screen model.This saturation effect resembles the behavior of a mixture model in optically thick case.That means, at very high inclinations, these nebular emission regions along the line of sight can no more be approximated by a point source in galaxy center, which suggests a more complex geometric configuration of these nebular line regions.* ρ g,0 and R cl are reduced as σ g,0 (σ g,0 = ρ g,0 πR 2 cl ).
As shown above in Section 3.2, a simple diffuse ISM dust component plus a simple spherical HII region can not properly account the observed E g −θ relation, especially for highly inclined disks.Moreover, the nebular dust component in principle will also have extinction effect on the stellar continua, which has not been taken into account by our simple mixture model presented in Section 3.1.In this section, we present a self-consistent two dust component model and use it to model the observed inclination dependence of the two different dust reddening features simultaneously.
As shown in Section 3.1, a uniform mixture model could generally reproduces the inclination dependence of the stellar reddening effect.Therefore, we still assume that the stellar emission comes mainly from a stellar disk , while the diffuse ISM dust component is uniformly mixed in.Indeed, there are observations show diffuse far infrared (FIR) emissions at high latitude of edge-on disk galaxies, whose scale-length and scale-height are linearly correlated with that of the stellar component (Mosenkov et al. 2022, and references therein).
The nebular emissions of SFGs are known to be mainly contributed by clumpy distributed HII regions.Because of heavy self-attenuation and short life-scales, the HII regions of local SFGs contribute few percent of stellar emissions.Observations show that the global structure of these HII regions (the youngest stellar population) also follows a disk geometry, but thinner and more extended than the general stellar populations (e.g.Bobylev & Bajkova 2021;Monteiro et al. 2021;Anderson et al. 2019).To approximate the geometry configuration of these clumpy HII regions, we assume that they are spherically symmetric and physically identical.The number density of these clumpy regions then follows an-other exponential disk distribution, which is co-planar and concentric with the stellar disk.We refer this geometry configuration of the clumpy HII regions as the clumpy nebular disk hereafter.Now, we have two different disk components for our model galaxy.One is the stellar disk composed of continuous stellar emission and its associated diffuse ISM dust, which has already been defined by Equation 7 in Section 3.1; the other is the clumpy nebular disk, composed of HII regions, each of which has an intrinsic optical depth τ cl for its nebular emission I g additionally.We refer to this two dust component model as the Chocolate Chip Cookie model (hereafter CCC model), where the dusty clumpy HII regions ("chocolate chip") are embedded in a stellar disk ("cookie").The schema of our newly proposed CCC model is shown in the top panel of Figure 4.
More specifically, we also parameterize the clumpy nebular disk as an exponential disk, where ρ g,0 now is the central number density of clumpy regions, R g , h g represent the nebular disk scale-length and scale-height respectively.However, unlike the continuous distribution of stellar disk, the distribution of nebular regions is clumpy, and the line of sight optical depth of clumpy regions is not a simple integration of Equation 15.Moreover, these clumpy distributed regions would also have dust extinction effects on the stellar emission.Therefore, we need further simplifications to model the dust attenuation effect of the clumpy regions on both of the nebular emission and stel-  lar emission, which we will discuss in the following subsections.

Modelling of the Clumpy Regions
We assume that each clumpy region is self-extincted by its dust shell with optical depth τ cl .When a clumpy region is in front of a emission source (e.g. a star or another clumpy region) along the line of sight, its mean (effective) optical depth is where R cl is the radius of a single clumpy region.We note that the optical depth τ is always a function of wavelength λ (Equation 9).Therefore, we henceforth omit the subscript λ for τ λ to avoid complexity of notation.
In our sample, the flux of a galaxy is collected from the SDSS fiber aperture with a radius of 1.5 arcsec4 , which corresponds to a radius of R a ∼ 2.2 kpc at redshift z ∼ 0.07 (the median redshift of our MW-like galaxy sample).We parameterize the number of clumpy regions in each SDSS aperture as N ∝ n g πR 2 a , where n g is the column number density of clumpy regions that equals to the integrate of ρ g along the line of sight.We note R cl is far smaller than R a and define a covering factor p as the ratio of the cross-sectional area of each clump to the aperture area, which also indicates the probability of a random emission source being covered by a foreground clumpy region along the line of sight.Because R cl << R a , the covering factor p << 1.We assume that there are N discrete clumpy regions randomly distributed within the aperture.Then, the number of foreground clumpy regions along a particular line of sight obeys a binomial distribution, Therefore, for an emitting source with intrinsic flux density I e , the observed flux density after extinction by these front clumpy regions is When τcl < 1 and p << 1, the discrete binomial distribution can be approximated by a continuous Gaussian distribution (N) with mean µ = N p and variance σ 2 = N p(1 − p) .Then, we obtain an approximation That is to say, the equivalent optical depth of foreground clumpy regions can be approximated by: Also, as p << 1 and if τcl < 1, the term pτ cl in Equation 21can be neglected and N p can be written as n g πR 2 cl .Thus, we finally get This final approximation shows that τcl is determined by the size R cl , the column number density n g and the mean optical depth τcl of an individual foreground clump.
During the derivation of Equation 19, we have assumed p << 1, N >> 1 and τcl < 1.In Appendix B, we provide a detailed comparison of the approximation of Equation 21to the numerical solutions of Equation 19 for different parameter sets, where we show that Equation 22 only has a difference from the numerical solutions at the level of ∼ 10% for τcl < 1.As we will show in Section 4.6.3,our best model estimation indeed have τcl < 1.
We remind that the continuous approximation of Equation 19 is only for solving Equation 20, not that the distribution of these star forming regions is itself continuous.Actually, the dust extinction effect of these clumpy regions is different from that of the continuously distributed dust.We show such a comparison in Appendix C.

Dust Extinction on Balmer Emission
With the geometric models of the clumpy nebular disk (Equation 15 and 22) and continuous ISM distribution (Equation 7), we model the overall dust extinction effect on a specific nebular region.For a given nebular cloud in a model galaxy along the line of sight to galaxy center, we parameterize its position with central distance l and disk inclination θ and get r = l sin θ, h = l cos θ, where l takes positive and negative values at the proximal and distal ends respectively.
For a given nebular region, its nebular emission is extincted by the dust of itself, foreground diffuse ISM dust, and other foreground clumpy regions.The total line of sight optical depth then is where τ cl , τ s and τcl are the optical depths of the clumpy region itself, the foreground diffuse ISM dust, and other foreground clumpy regions, respectively.For the double exponential diffuse ISM dust disk parameterized by h s and R s (Equation 7), the optical depth of the foreground ISM dust along the line of sight at (l, θ) is where τ s,A is the total optical depth along the central line of sight which follows Equation 10.For τcl , according to Equation 22, n g is the only parameter being function of (l, θ).With the global exponential disk distribution model of the clumpy regions (parameterized by h g , R g ), we have where n g,A (θ) is the total column number density of clumpy regions along the central line of sight:

Dust Extinction on Stellar Emission
Similar to Equation 23, the line of sight optical depth to a stellar emission region at (l, θ) is contributed by the foreground diffuse ISM dust and clumpy dust, Obviously, τ s (l, θ) and τcl (l, θ) follow Equations 24 and 22, respectively.

Off-axis Effect
In Sections 4.2 and 4.3, we have parameterized the foreground dust extinction on the stellar and nebular emissions at a specific region (l, θ) along the galaxy central line of sight.However, the SDSS fiber has an aperture of 1.5 arcsec instead of being an area-free point toward the center of each galaxy.Therefore, there is a certain fraction of photons not along the line of sight to galaxy center, which also enter the fiber.Because the SDSS fiber aperture is smaller than the radius of a typical disk galaxy in our sample, the on-axis (line of sight to galaxy center) assumption is a good approximation when galaxies are viewed in face-on cases.However, when galaxies are viewed edge-on, the scale-height of our sample galaxies becomes comparable or even smaller than the SDSS fiber aperture.In this case, the off-axis effect (line of sight not along the galaxy center), as we show next, will cause significant biases.
We assume that the scale-height of the nebular emission regions h g is smaller than that of the stellar disk h s and both scale-heights are smaller than the SDSS aperture size.In this case, as shown by the schema figure (top panel of Figure 4), along the line of sight to the galaxy center, the stellar emissions suffer from the dust extinction effects of both of the clumpy nebular dust and diffuse ISM dust.While for the stellar emissions off the disk plane, the dust extinction from nebular regions becomes negligible.
To approximate and correct for this off-axis effect, we assume that there are a fraction of photons ( f off ) emitted from the stellar disk being off-axis, which is defined as: These off-axis stellar emissions are not extincted by the nebular dust.On the other hand, the emissions from the nebular regions obviously do not have this off-axis effect.We consider a case where the disk inclination θ is relatively large but not completely edge-on.In this case, the effective disk height in the line of sight is expressed as As we will show later, the nebular disk has larger scale-length than that of the stellar disk.As a result, with the decreasing of disk inclination θ, the effective height of nebular disk (h g ) starts to approach that of the stellar disk (h s ).We write the specific inclination when h g = h s as θ crit .We assume that there is no off-axis effect on the stellar disk when θ < θ crit .Actually, when disk galaxy becomes face-on, both h s and h g will be larger than the SDSS fiber aperture, and thus there is no off-axis effect on both disks.In summary, because of the off-axis effect for the highly inclined disk galaxies, inside the SDSS fiber, a fraction of off-axis photons emitted from the thicker disk have not been extincted by the thinner component.In specific, this fraction is parameters by with above modeling, the dust extincted stellar intensity at given position (l, θ) is written as where τ s and τ s are given by Equations 24 and 27 respectively and I s represents the stellar emissivity in arbitrary unit.For nebular emission, there is no off-axis effect.Therefore, the dust extincted intensity is where τ g is given by Equation 23.

Global Dust Attenuation
In the above sections, we have modelled the dust extinction effect on both nebular emission and stellar emission for a given point (l, θ) in model galaxy along the line of sight.Next, we model the global dust attenuation effect by integrating the dust extincted radiation along the line of sight while neglecting the dust scattering effect, For a given inclination θ of a disk galaxy, the unreddened intensity along the line of sight to galaxy center is In Equation 33and Equation 34, the subscript 'comp refers to 's for stellar disk and 'g for nebular disk.With the unattenuated and attenuated nebular emission line intensity, the dust attenuation then is easily defined as and so that the reddening is In above equations, we have assumed the dust extinction curves of both clumpy nebular and diffuse ISM dust follow the same simple power-law (Equation 9).

Reducing free parameters
In the CCC modeling, we have used 11 parameters in above equations, which are listed in Table 1.For modelling the stellar disk, we have used h s , R s , α s,0 , I s (Equation 10, 24, 27, and 31).For the nebular disk, we have used h g , R g , ρ g,0 , I g , τ cl , R cl , R a (Equation 22, 25, and 32), where we remind that τ cl is at the wavelength 5,500 Å and is directly related to the effective optical depth τcl through Equation 16.Since the dust attenuation describes the ratio of the emission with and without dust, the emissivity I s and I g are therefore unnecessary parameters.For simplicity, we set I s and I g to be 1.
Then, the dust attenuation value at a given inclination is determined by two sets of parameters: the geometric parameters:h g , R g , h s , R s , R a , R cl and dust-related parameters: α s,0 , ρ g,0 , τ cl .For these geometric parameters, their absolute values are actually degenerated with the central density parameter (e.g.Equation 10, 24, and 26).That is to say, the absolute sizes of our model galaxies can be varied arbitrarily, provided that their corresponding density parameters are adjusted accordingly.Therefore, following Section 3.1, we normalize these geometric parameter sets by setting the stellar scale-length R s = 1.Then, the rest scale parameters h g , R g , h s becomes equivalent scale-length and scale-height in unit of R s .The corresponding density parameter α s,0 and ρ g,0 become equivalent densities in unit of R −1 s .In addition, during the modeling of the dust extinction effect of the clumpy regions, the size R cl of the clumps is coupled with their column number density (Equation 22).Therefore, in our final modeling of global dust attenuation effect, R cl is coupled with the central number density ρ g,0 .
Considering this coupling effect, we define a new parameter σ g,0 ≡ ρ g,0 πR 2 cl , which represents the cross-section of the clumpy regions in galaxy center.
In summary, our model has 6 independent free parameters to be constrained.Among them, there are 3 equivalent scaleheight and scale-length parameters: h g , R g , and h s ; central absorption density of diffuse ISM dust α s,0 ; central crosssection σ g,0 of the clumpy regions and the optical depth of an individual clumpy region τ cl .

Parameter constraints
We use CCC model to fit the observed inclination dependence of E g (θ) and E s (θ) simultaneously, where a MCMC algorithm is used to make the best estimates of six free parameters.The probe of the model parameters is set in logarithmic space.Also, in order to obtain physically meaningful parameter values, we set the following constraints on the model parameters.
• Both nebular and stellar disk shall have a larger scalelength than scale-height: R s > h s , R g > h g ; • The scale-length of the two disks shall be comparable: 0.5 < R s /R g < 2; • The scale-height of nebular disk shall be smaller than that of stellar disk: h g < h s ; • The optical depth of SF regions shall follow the optically thin approximation: τ cl < 1;

Modelling results
After excluding invalid MCMC chains, we get excellent fittings to the observed E g − θ and E s − θ relations and strong constraints on all six model parameters.We show the corner plot of the MCMC fitting in the bottom panel of Figure 4 and list the fitting results in Table 1.As can be seen, the CCC model fits the E g − θ relation of these highly inclined disks (θ > 70 • ) very well, where the simple screen model has failed, as described in Section 3.2.
For the diffuse ISM component, we get the best fits that the height-to-length ratio h s /R s = 0.19 and the central absorption density of ISM dust α s,0 = 2.58.For CCC model, using Equation 10, we can easily estimate that the full optical depth of the diffuse ISM component along the galaxy center increases from τ s (0 • ) = 0.95 with face-on inclination to τ s (90 • ) = 5.16 with edge-on inclination.We note that the height-to-length ratio in our model results is slightly larger than that of reported in literature, which is around 0.1-0.15(Shao et al. 2007;Guthrie 1992).The specific reasons for this discrepancy will be discussed in sections 6.1 and 6.4.For the global distribution of nebular emission line regions, the CCC model presents a very thin clumpy nebular disk, which has a height-to-radius ratio,h g /R g = 0.066, and is much thinner than that of diffuse ISM disk.More specifically, in our modelling, the scale-length of the nebular disk is about 1.6 times larger than the diffuse ISM disk, while the scale-height is only about half of the diffuse ISM dust disk (h g /h s = 0.55).The decrease of the Balmer decrement at high inclination requires that the nebular clumpy disk is more extended than the diffuse ISM dust component, such that the outer disk clumps are less attenuated by ISM dust .
For individual clumpy regions, our model shows that the number of clump regions at galaxy center per stellar scalelength R s (σ g,0 ≡ ρ g,0 πR 2 cl ) is about 1.8.Therefore, we expect that there are only about 0.37 clumps along the central line of sight when our model galaxy is face-on ( Equation 26).This number increases to 5.61 when the model galaxy is edge-on, which means an overlapping effect of HII regions.That is right the reason of the failure of the simple screen model at high inclinations presented in Section 3.2, where no overlapping effect has been assumed.
For the optical depth of individual emission line region in the CCC model, we get the best estimate of τ cl = 0.5, which is very close to the result we have obtained in the simple screen model (Section 3.2).This result is not surprising.When the model galaxy is face-on, as we have discussed, there is essentially few overlapping effect of the HII regions along the line of sight, and the only extra dust extinction in addition to the ISM dust that affects an emission line region is only from the local dusty shell of this region itself.In this case, the CCC model degenerates to the dust screen model in Section 3.2.
Besides, in our modelling, the critical viewing angle where the effective scale-height of nebular disk equals to that of diffuse ISM disk happens at θ crit ≈ 83 • .That is to say, only for these very highly inclined disk, we have considered the off-axis effect for the nebular dust extinction observed by the SDSS fibers (Section 4.4).

APERTURE EFFECT AND DUST ATTENUATION
As we have shown in Figure 2 and discussed in Section 3.2, when a galaxy disk is highly inclined, there is a saturation effect on the observed dust reddening of the nebular emission lines, which is caused by the fact that the inner emissions would be nearly completely extincted by the outer dust layers along the highly inclined line of sight.That is to say, in observation, we can not distinguish whether such a saturation effect has occurred or not from color excess values alone.On the other hand, if we further consider the extincted flux, we would easily find how much of the emission flux has been extincted by dust and then distinguish the saturation effect of color excess.We show the inclination dependence of the r band absolute magnitude and Hα flux inside fiber aperture of our sample galaxies in Figure 5.At given inclination, the median absolute fiber magnitude and Hα flux are shown by the solid dots, which can be viewed as the inclination dependence of the observed fiber magnitude and Hα flux for a typical MWlike SFG at typical redshift z ∼ 0.07.As can be seen, both of the fiber magnitude and Hα flux do not change when the disk inclination is low (θ <∼ 50 • ) and then increases/decreases monotonically with increasing disk inclination for the fiber magnitude and Hα flux respectively.Overall, the fiber magnitude dims about 1 mag and the Hα flux drops about 0.5 dex from the face-on view to edge-on view.
We use the best CCC model to calculate the attenuated nebular emission line flux densities and r band stellar surface brightness along the central line of sight for different θ val-ues.The output (attenuated) r band stellar surface brightness and Hα flux density are plotted as the dashed curves in the top and middle panels of Figure 6.The emission line flux and surface brightness have been both normalized to zero at θ = 0 • .When dust attenuation effect is considered, the attenuated central flux density from the CCC model prediction generally shows a plateau at low inclination (θ <∼ 50 • ), which is consistent with observations.While for high inclinations (θ >∼ 50 • ), the CCC model predicts a increasing central flux density with disk inclination, which apparently is contrary to the decreased trend of the aperture flux seen by the observational data.
The inconsistency between the increasing trend of the central flux densities and the decreasing trend of the aperture fluxes can be explained by the aperture effect.In other words, what we observe is the total flux inside the fiber aperture, while our model prediction is the flux density along the central line of sight.To estimate the observed flux correctly, we need to take the aperture effect into account.
For our sample galaxies, the physical aperture size of SDSS fibers (R a ∼ 2.2 kpc) brings important differences between the observed flux inside a fiber and the central line of sight flux density in our model.On one hand, the projected nebular emission line (or stellar continuum) flux inside a fiber is not a constant, but has a profile depending on both the disk inclination and dust geometry in a complicated way.On the other hand, when the disk is highly inclined, the projected disk height (especially that of clumpy nebular disk) will not fill all of the fiber aperture.We show the fiber aperture effect for galaxies with different inclinations (three cases) with schematic icons in the bottom panel of Figure 6.
Detailed analytical modeling of the projected emission line and stellar continuum fluxes inside fiber aperture as a function of disk inclination is highly complex and beyond the scope of this study.Here, we make approximations based on simplified assumptions.For brevity, we illustrate below the process of estimating the emission line flux as function of different inclinations.For the stellar emission (r band fiber magnitude), the estimation process can be referred exactly to that of the emission lines.
We assume that the observed nebular emission line flux is proportional to the product of the central surface brightness I cen and disk projected area in aperture S , With CCC model, we have derived the central flux density I cen,comp as a function of disk inclination (dashed line in middle panel of Figure 6).Next, we discuss the values of S at different inclinations.For a MW-like disk galaxy with stellar mass 10 10.4 M , the scale-length of stellar disk is about 2 kpc (Shen et al. 2003).Adopting this value as the scale-length of stellar disk the observed median fiber magnitude in θ bins (dots), the attenuated central surface brightness predicted by the CCC model (dashed line), the aperture corrected fiber magnitude after applying aperture corrections for the model galaxy with R a /h g = 10(solid line).Middle panel: the observed median Hα flux in θ bins (dots), the attenuated central nebular line fluxes predicted by the CCC model (dashed line), the aperture corrected nebular line fluxes after applying aperture corrections for three model galaxies with different physical sizes (solid lines, different colors represent different h g values).Bottom panel: the projection of a double-exponential disk into a fiber with radius R a (2.2 kpc) at different inclinations.The three solid lines show the projected area S as function of disk inclinations for three model galaxies, which are all set to unit value at θ = 0 • .The three schematic icons show the projection of a model galaxy with scale-length 3.2 kpc and scale-height 0.2 kpc into the fiber aperture (2.2 kpc) for face-on (case I), edge-on (case II) and critical (case III) inclinations respectively.R s in CCC model, we obtain that h s , R g , h g of the disk of our model galaxy are 0.37 kpc, 3.17 kpc and 0.21 kpc respectively.When the model galaxy is face-on, the emission line disk scale-length is larger than the fiber aperture, the projected disk will fill all of the fiber aperture and thus S = πR 2 a (case I in Figure 6).When the model galaxy is edge-on, the projected area is determined by the scale-height of the disk.In this case, we approximate the projection area as S ≈ πR a h g (case II in Figure 6).For a disk with an arbitrary inclination θ, following Equation 29 of Section 4.4, its equivalent projection scale-height can be written as h comp (θ) = h comp sin θ + R comp cos θ as the approximation of the projected height of disk, we can consistently write S as: where 'Min' takes the minimum value between R a and h comp (θ).With this approximation, there is a critical inclination angle θ crit where h comp (θ crit ) = R a (case III in Figure 6).When the model galaxy is observed from face-on to critical inclination θ crit (θ < θ crit ), we have h comp ≥ R a , and S equals to the covering area of fiber, πR 2 a .When θ > θ crit , S = πR a h comp .
In the middle panel of Figure 6, we plot the model predicted Hα fluxes inside fibers as functions of disk inclinations for three different chosen R a /h g values (5,10,15) using solid lines of different colors.Since R a has a physical size ∼ 2.2 kpc, these three different R a /h g ratios imply three different physical sizes of the model galaxies.Because the best estimate of the CCC model have R g /h g = 1.6, for these three R a /h g values (5,10,15), we have θ crit ∼ 75 • , 55 • , 25 • respectively.For each R a /h g value, we then normalize the model-predicted mean Hα fluxes for inclinations between 0 and θ crit to unit value.When R a /h g = 5, θ crit ∼ 75 • , because the central flux density I cen still increases in this inclination range (larger than 0.1 dex), our model does not predict a very flat plateau in this inclination range.Moreover, from θ crit to edge-on, the drop of the model predicted Hα flux is also not as large as that being observed.On the other hand, when R a /h g = 15, because then R g < R a , our model predicts a continuous decrease of the effective projection area (Equation 37) and resulted Hα flux, which looks also not be in good consistence with observations.Finally, we see that R a /h g = 10 provides a fairly good prediction on the global behavior of the observed nebular emission line flux as a function of disk inclination.In fact, we also have calculated the sum squared residual values of these three different model lines to the observational values in each θ bin, and find that the line of R a /h g = 10 is indeed the best.Moreover, adopting R a /h g = 10 and replacing the parameters with those of the stellar disk, we can make a similar aperture effect correction for the observed r band fiber magnitude using Equation 36,37.The result is shown as the solid line in the top panel of Figure 6.Again, we see that, after correction of aperture effect, our CCC model makes an excellent prediction on the inclination dependence of dust attenuation effect for r band stellar continuum.It is worth emphasizing that the constraints on the parameters of our CCC model are from the the reddening features of the emission lines and the stellar continuum, without using their attenuation features.In the above, we show that the CCC model predicted relations between the attenuation features and disk inclination are also very consistent with the observations once the fiber aperture effect is probably accounted.This result further illustrates the internal self-consistency of our CCC model in predicting the dust attenuation and reddening features of local disk galaxies.
The excellent consistence of R a /h g = 10 model with observations provides an interesting constraint on the physical size of our model galaxy.In Section 4, because of the degeneration of model parameters, we can only obtain the relative geometry parameters (all in unit of R s ) for our model galaxy.Here, because of the average aperture size R a ∼ 2.2 kpc as we have discussed, we naturally obtain h g = 0.22 kpc and then get estimates of all the model parameters in physical units.We list them in the last column of Table 1.With the physical sizes of all these geometric parameters, we can explore further the physical implications of the CCC model.

DISCUSSION
In Section 4, we have obtained the best estimates of the geometric parameters of the CCC model in units of R s .Besides, we also get a constraint that each clumpy region has optical depth of ∼ 0.50 in V band.In Section 5, we find that, by using the aperture effect, we get estimates of the model parameters in physical units: R s ∼ 2.1 kpc, R g ∼ 3.33 kpc, h s ∼ 0.41 kpc, h g ∼ 0.22 kpc, α s,0 ∼ 1.22 kpc −1 , σ g,0 ∼ 0.84 kpc −1 for modeled MW-like disk galaxies.In this section, we compare our model estimates with other observational results of local massive disk galaxies and make further discussions.

Dust Geometry
We first compare the geometry parameter of our model galaxy with observational or modelling results of the Milky-Way and other nearby disk galaxies.
For the overall structural parameters of the dust component of the Milky-Way, many studies have reached consistent conclusions that dust is thinner and more extend than stars (e.g.Drimmel & Spergel 2001;Misiriotis et al. 2006;Li et al. 2018).For the extra-galactic galaxies, most of the studies on the dust geometry in optical wavelengths also have assumed only one global dust component and obtained similar conclusions.For example, Xilouris et al. (1999)(hereafter X99) investigated the surface brightness profiles of five nearby edge-on galaxies in B, V, and I bands, Bianchi (2007) analyzed another seven nearby edge-on galaxies in V and K bands.By applying radiative transfer analysis on these galaxies, they find that the radial scale-length of dust is about 1.4 times larger than that of stars, while its vertical scale-height is about half of the stellar disk.More recently, De Geyter et al. ( 2014) studied the SDSS g, r, i, and z band images of 12 edge-on spiral galaxies selected from the CALIFA survey and found similar conclusions.Besides edge-on galaxies, Casasola et al. ( 2017) investigated 18 face-on spiral galaxies in DustPedia from UV to sub-millimeter bands and found that the dust scale-length is about 1.6 times of the stellar one.For our model galaxy, the clumpy nebular disk has a larger scale-length and smaller scale-height than the stellar disk, while the diffuse ISM dust component has been assumed to follow the same geometry as the stellar disk.Qualitatively, merging the dust in the clumpy regions with that in the diffuse ISM will give a global dust component that is larger in scalelength and smaller in scale-height than the stellar component, which is consistent with other studies.However, considering the non-linearity of the dust attenuation effect, the geometric parameters of different dust models are not quantitatively comparable.
To quantify the global dust attenuation effect of the two dust components in our CCC model, we reconstruct the projected image of our model galaxy from the edge-on view using the best model parameters.More specifically, the image projection process is an integral of the dust attenuated stellar emissions along the line of sight, which is given by where ρ s and τ denote the density of stellar emission and dust optical depth at given position (X, Y, l z ) respectively, and (X, Y, l z ) is a Cartesian coordinate system with the centre of the model galaxy as the origin.(X, Y) forms the projection plane and X, Y is the major and minor axis along the disk plane,separately, while l z is the axis perpendicular to the projection plane.This Cartesian coordinate system can be easily mapped to the cylindrical coordinate system (r, h) that has been defined in Equation 5through As can be seen from the Equation 38, the resulted image is not only a function of geometric parameters (e.g.scalelength and scale-height, discussed in detail below), but also is a function of the normalization parameter of dust optical depth (or density).In addition, depending on the model assumptions, ρ s (X, Y, l z ) and τ(X, Y, l z ) can both be combined by multiple components.For example, ρ s of our CCC model is only an exponential disk characterized by R s , h s , while τ is combined of two disks, the diffuse dust (c.f.Equation 8) and the clumpy dust (c.f.Equation 22).
With above equations, the resulted V-band edge-on model galaxy image of our best fit CCC modelling is shown in the top panel of Figure 7.To do the projection, we have set the apparent magnitude of our model galaxy without any dust attenuation to be 17.5 mag and its stellar disk scale-length R s = 1.5 arcsec (typical values of our sample galaxy).In this reconstructed image, a dust lane structure along the mid-dle plane of the model galaxy is clearly seen, which is originated from the extra obscuring effects of the clumpy dust component.To further quantify the global properties of the projected image, we plot the effective optical depth τ eff5 and the surface brightness profiles along the galactic mid-plane in the bottom two panels of Figure 7 as dotted lines.With these two quantitative profiles, we make more detailed discussions on the dust geometry of our best fit CCC model and further compare it with other dust attenuation models of disk galaxies.

Equivalent single dust component model
Since many early studies of the dust geometry of disk galaxies considered only one dust component, we are interested in testing whether the two dust components in our CCC model can be equated with one dust component.To test this idea, we consider a simple model of a continuously distributed double exponential dust disk mixed with a stellar emission disk which is assumed to be the same as that of the CCC model.To predict the observed surface brightness profile for this single dust component model, as that shown in the bottom panel of Figure 7, three model parameters are needed: the dust-to-stellar scale-length ratio R d /R s , the dust-to-stellar scale-height ratio h d /h s and a normalization parameter representing the dust density.For this dust density parameter, follow convention, we take the V-band the central optical depth of the model galaxy in face-on view τ f .By adjusting these three model parameters, we find that R d /R s ∼ 1.1, h d /h s ∼ 0.7 and τ f ∼ 0.4 provide almost the same surface brightness profile as that of the CCC model prediction, which is shown as the red dashed lines in the bottom two panels of Figure 7.For this equivalent single dust component, the dust scale-length is about 10% larger than the stellar scale-length.This result is in good agreement with that of Muñoz-Mateos et al. (2009) for nearby disk galaxies, where the global dust scale-length is obtained from the modelling of the infrared emissions.
From Figure 7, we conclude that our two-component CCC model can indeed be equated by a single dust component, which is thinner and larger than the stellar component, as expected.We compare this equivalent dust component quantitatively with other studies in the next subsection.On the other hand, this result also implies that the projected image alone may not be sufficient to recover the detailed structure of the dust component.To reveal the three dimensional dust geometry of galaxies, a comprehensive study on the different dust attenuation properties is needed.

Compare with other dust attenuation models
In this sub-section, we compare the dust geometry of our model galaxy with the result of other model galaxies.There are two widely used dust geometric models of disk galaxies in literature, one is the single dust component model of X99 and the other is the two dust component model of T04.
To have a consistence comparison with the CCC model, we set the stellar disk scale-length R s = 1.5 arcsec and the total un-attenuated stellar emission to 17.5 mag for both of the model galaxies in X99 and T04.

X99 model -
The dust attenuation model of X99 is the same as the single dust component model we discussed in Section 6.1.1,with the only difference being the inclusion of an additional bulge component in its stellar emission.
For the specific model parameters of X99, we adopt the typical values of their sample galaxies: B/T ∼ 0.2, stellar scale-height to scale-length ratio h s /R s ∼ 0.1, dust to stellar scale-height ratio h d /h s ∼ 0.5, dust to stellar scale-length ratio R d /R s ∼ 1.4,V-band face-on optical depth τ f ∼ 0.5.We then make its V-band edge-on view image and show it in the second top panel of Figure 7.The corresponding effective optical depth and surface brightness profiles along the galactic plane are shown as the solid curves in the bottom two panels.
As can be seen from the projected images, the X99 model shows a more pronounced dust lane than CCC model.The effective dust optical depth along the galactic plane of the X99 model is systematically larger than that of the CCC model (∆τ eff ∼ 1).However, for the surface brightness profiles, it is interesting to see that our CCC model prediction is quite close to that of the typical galaxy in X99 (∆µ < 0.5mag • arcsec −2 at all radii).
Before further discussion, it is worth reminding that the geometry of the stellar emission of the typical X99 model galaxies is different from that of the statistical MW-like galaxy in our CCC model.First of all, the CCC model galaxy does not include a bulge component whereas the X99 model galaxy has B/T ∼ 0.2.Second, the stellar scale-height to scale-length ratio of the CCC model galaxy (h s /R s ∼ 0.2) is larger than that of the X99 model galaxy (h s /R s ∼ 0.1).The reasons for these differences are twofold.On the one hand, the higher h s /R s of the CCC model is partly a compensation for the absence of the bulge component in its model assumption (see more discussions in Section 6.4.1).On the other hand, the h s /R s of CCC model is fitted from a statistical sample of MW-like galaxies, while X99 model only fits 7 nearby edge-on galaxies.That is to say, the sample galaxies in these two models may not be comparable.

Lu et al.
The comparable µ profiles of two models are combined results of different effects.First, due to the higher h s of the CCC model, the un-attenuated surface brightness profile along the galactic plane of the CCC model is fainter.On the other hand, the scale-height of the equivalent single dust component of the CCC model galaxy (h d ∼ 0.14R s , Section 6.1.1)is about twice of that of the X99 model (h d ∼ 0.07R s ), which makes effective optical depth along the galactic plane to be significantly larger in the X99 model.Because of this significant large τ eff , the stellar emission of the X99 model galaxy is more attenuated, so that it has a comparable and even fainter µ profile along the galactic plane.These complementary effects indicate the degeneracy between the stellar emission and the dust component in the modeling of the dust attenuation of galaxies.In other words, in the dust attenuation model, it is better not to predetermine the geometry of the stellar emission component, otherwise the geometric properties of the dust component obtained by modelling could be biased.

T04 model -
The T04 study provides a framework of dust attenuation process with multi-components, where the stellar composition is composed of a bulge, a thick disk of old stars, and a thin disk of young stars, and the dust composition also includes both a thick and thin disk respectively.The thick dust disk, which represents the continuously distributed ISM dust, is thinner and more extended than the old stellar disk, while the thin dust disk has the same geometry as the young stars and represents the dust associated with new-born stars.Most of the geometric parameters of the T04 model use the values of the nearby galaxy NGC 891 obtained from radiative transfer modelling (see T04 for detail).The only three free parameters excepted inclination remained in T04 are: the face-on optical depth τ f , bulge-to-total ratio B/T, and clumpiness F, where F is defined as the total fraction of UV light being locally absorbed by dust in the thin disk.By giving different parameter settings for these three parameters, the T04 model is capable of describing the attenuation of galaxies at different inclinations and and therefore is widely used to study the attenuation-inclination dependence of disk galaxies(van der Giessen et al. 2022;Driver et al. 2007;Masters et al. 2010).
We also make the edge-on projection image for MW-like galaxies using the T04 model.Specifically, we first take the basic geometric parameters from Table 1 of T04.For the free parameters, we take the values from Table 3 of G22, where the T04 model parameters for the MW-like galaxies in SDSS have also been constrained from the stellar attenuation-inclination relation: B/T ∼ 0.21, τ V,f = 3.05, and F = 0.34 .Here, we convert the B band face-on dust optical depth τ B,f to V band by dividing a factor of 1.32 (for R V = 3.1 extinction curve).The clumpiness factor F is not related to the dust attenuation of old stellar population, and therefore does not plays a role in the projection.The edgeon view image projected from this T04-G22 model galaxy is shown at the top third panel of Figure 7, whereas its τ eff and µ profiles are shown as the two dashed lines in the bottom two panels respectively.
As can be seen, the dust lane of this T04-G22 model galaxy is extremely prominent.As a result, its µ profile is significantly fainter than both of the CCC and X99 model galaxies.Moreover, in this T04-G22 model galaxy, the µ profile is almost a constant out to 4 times R s , which means that the galactic plane is optically thick even to its very out region.This very optically thick surface brightness profile along the galactic plane (or the very prominent dust lane) of the T04-G22 model galaxy is a combined result of the relatively large dust content (τ V,f ∼ 3) and the very small height-to-length ratio (0.016) of the thin dust disk preset in T04.A detailed discussion of the origin of this atypical profile of the T04-G2 model galaxy is beyond the scope of this study.However, we will discuss more about the dust attenuation on emission lines of the T04 model in Section 6.4.2).

Optical Depth
In the CCC model, besides the geometric parameters, we also have obtained constraints on the optical depth of two dust disks.
For the diffuse dust disk, we get an estimate of the central absorption coefficient α s,0 ∼ 1.22 kpc −1 .However, this central parameter α s,0 in our model is used more as a normalization parameter to describe the optical depth of the diffuse dust component at different regions rather than has an unambiguous physical implication of its own.The reason is that our exponential disk model is too simplified to describe the central region of a real galaxy that contains other complicate physical components, e.g., bulge, nuclear star cluster, and active galactic nuclei etc., which have not been taken into account.On the other hand, it is worth using α s,0 to estimate the optical depth of typical regions of our the model galaxy.For example, at "solar neighborhood" (the location with a distance of ∼8.3 kpc from the model galaxy center on its galactic plane), our model predicts an absorption coefficient of 0.02 kpc −1 .Then, by integration of the absorption coefficient along the line of sight to high galactic latitude regions, we obtain a line of sight optical depth ∼ 0.01 for diffuse ISM dust, which is in good consistence with the SFD map data of our MW (Schlegel et al. 1998).
For the clumpy nebular disk, we get estimates of the optical depth τ cl ∼ 0.5 for each clump and the central absorption cross-section of clumpy regions σ g,0 ∼ 0.84 kpc −1 .We remind that σ g,0 is the product of the central number density of the clumpy regions and projection area of each clump, ρ g,0 πR 2 cl .For σ g,0 , similar to the argument for the diffuse ISM dust, our estimate also can not be directly compared with the observations.However, we may use it to further probe the overall properties of the clumpy HII regions.Before that, we need an estimate of the size of the individal clumpy region, R cl .The sizes of Galactic HII regions detected by their middle infrared (MIR) emissions are about ∼ 10 pc (Anderson 2014), while the sizes of giant HII regions detected in the nearby disk galaxies are shown in the range 10 ∼ 100 pc in the optical wavelength (Gutiérrez & Beckman 2008).The different sizes of HII regions detected in Galactic and extra-galactic disks may reflect the clumpy nature of the HII regions and possible selection effect 6 .Theoretical studies on the Strömgren sphere also show that the size varies greatly with different H atom density and in different luminosity class of OB stars, ranging from several pc to nearly a hundred pc (Gutiérrez & Beckman 2008).
Here, for simplicity, we assume R cl ∼ 30 pc and then obtain ρ g,0 ∼ 300 kpc −3 .With this number, we can further obtain an estimate of the total number of clumps inside the clumpy disk, N = 4πR 2 g h g ρ g,0 ∼ 9.2 × 10 3 , where we have adopted R g = 3.33 kpc, h g = 0.22 kpc in physical units.If we assume that each clump has 10 3 M newly formed stars (a HII region or open cluster usually contains ∼ 10 3 stars, initial mass function(IMF) also shows that only about 0.1% stars are OB stars), we obtain a total of ∼ 10 7 M new-born stars in our model galaxy.Considering the fact that the typical age of a HII region is about 10 Myr, and assuming that all new-born stars are formed inside the HII regions, the star formation rate of our model galaxy is then S FR ∼ 1M /yr, which is in excellent agreement with that of the local mainsequence star-forming galaxies (Brinchmann et al. 2004).In above discussion, we have used approximations for both R cl and the number of newborn stars in each HII region.However, we would like to emphasize that our final estimate of the star formation rate does not significantly depend on the specific values of these two parameters.The reason is that these two parameters are physically correlated.If we assume a smaller R cl , the number of newborn stars in each HII region should also be smaller.
For the optical depth of clumpy regions, our model constraint, τ cl ∼ 0.5, is close to the typical value of the Galactic HII regions (Sun et al. 2021), but at the lower limit of the observed HII regions in extra-galactic galaxies.(Gutiérrez & Beckman 2008).Considering the clumpy nature of HII regions and the selection and resolution effects of observations, the HII regions observed in extra-galactic galaxies are more likely biased to giant HII regions or HII region groups.Indeed, as we have discussed, the sizes of HII regions reported in extra-galactic studies are systematically larger than that of Galactic ones.
With the CCC model, we may further have an estimate of the total amount of dust in each dust component.The total amount of dust in the diffuse ISM is where κ V is the mass absorption coefficient and in unit of kpc 2 /M .For clumpy HII regions, the total dust amount is Combining these two estimates, we conclude that, for MWlike galaxies, the total amount of dust in clumpy regions is less than but comparable to that of the diffuse ISM dust.
Here, it is worth mentioning that we have not considered the fully optically thick HII regions in above discussion (see more discussions in Section 6.4.2).

Attenuation Curves
In the CCC model, we have used a power law extinction curve (Equation 9) with R V = 3.1 for both the diffuse ISM dust and clumpy dust.However, during the SPS fitting to derive the E(B − V) from the stellar continua of our sample galaxies, we used the Calzetti attenuation curve with R V = 4.05 (Calzetti et al. 2000).We plot these two different extinction (attenuation) curves on the top panel of Figure 8, which are significantly different from each other at short wavelengths (λ < 5000Å).The difference is mainly due to the differences of their definitions, which we discuss in detail in following.
The extinction curve is only determined by the physical and chemical properties of the dust particles, while the attenuation curve, or the effective extinction curve, is further correlated with the geometrical distributions of both the dust particles and radiation sources (e.g., Calzetti 1997;Calzetti et al. 2000;Witt et al. 1992;Witt & Gordon 1996, 2000).Since the CCC model gives a full description on the geometry of the dust particles and radiation sources, we can easily derive the shape of the attenuation curves for our model galaxy at different inclinations with the assumed power-law extinction curve (R V = 3.1).In our modelling, we used E g (Hα − Hβ) to represent the dust reddening of the nebular emission (Equation 5), which is independent of the extinction curve being assumed.Therefore, we only need to consider the shape of the attenuation curves for stellar continua.
We use the CCC model and take the best estimates of the model parameters (Table 1) to calculate the effective dust attenuation for the stellar continua of our model galaxy at dif-  9) used in the CCC model (red curve) and the Calzetti attenuation curve (blue curve).Bottom panel: the attenuation curves derived from the CCC model for the edge-on case (dashed, θ = 90 • , R V = 7.0), face-on case (dotted, θ = 0 • , R V = 3.7) and median case (solid,θ = 60 • , R V = 4.1), the blue curve also shows the Calzetti attenuation curve for comparison.ferent wavelengths and then derive the shape of the attenuation curve.We show the resulted attenuation curves for three representative inclination angles, θ = 0 • , 60 • , and 90 • , as the red dotted, solid and dashed lines in the bottom panel of Figure 8, which have R V = 3.7, 7.0, and 4.1, respectively.We see the CCC model naturally predicts an increasing of R V with increasing disk inclination, which has been reported in observational studies (e.g.Battisti et al. 2017).The increasing of R V of the attenuation curve with disk inclination (optical depth) is caused by the saturation effect when the emission sources are mixed with the dust (see Equations 11 and 12).Also, for the median disk inclination (θ ∼ 60 • ), the CCC model predicts an attenuation curve with R V = 4.1, which is in excellent agreement with the classical Calzetti attenuation curve (Calzetti et al. 2000) with R V = 4.05.Moreover, the ranges of the attenuation curves R V ∈ (3.7, 7.0) predicted from the CCC model is also matched with the observational results of nearby galaxies (Calzetti 1997).
In Section 2.2.2, we used the Calzetti attenuation curve with an constant R V = 4.05 in the SPS fitting.To make the picture fully self-consistent, we should apply different attenuation curves for galaxies with different inclinations.How-ever, introducing the inclination effect comprehensively into the SPS fitting process is beyond the scope of this work.Nevertheless, we expect that this inconsistency does not have a significant impact on our conclusions.For example, the derived median color excess E s (B − V) of a sample of nearby galaxies only changes from 0.15 to 0.16 when the attenuation curve is changed from R V = 4.05 to 4.88 (Calzetti 1997;Calzetti et al. 2000).

Caveats
Our CCC model has not only presented excellent fits to the complicate inclination dependence of both the nebular and continuum reddening features(Figure 4), but also give consistent predictions on their dust attenuation effects (Figure 6).However, our model is also subject to limited observational constraints and uncertainties in the model assumptions.We discuss these caveats in the CCC model below.

Bulge component
For simplicity, our CCC model has not considered the bulge component of disk galaxies.A reasonable assumption is that the bulge component does not have star formation, so there is no cold gas and dust associated.In this case, the dust attenuation of the nebular lines E g would be independent of whether or not including a bulge component in our model.For stellar emission, if we assume that the bulge component is concentrated and spherical, then the combination of a spherical bulge and a thin disk is essentially equivalent to a slightly thicker disk.This is precisely the reason for the relatively high stellar disk scale-height to length ratio, h s /R s ∼ 0.19，obtained by our CCC model.A more refined geometric model with a bulge component introduces more free parameters, and therefore requires more observational constraints, e.g. the combination of multi-wavelength images, and/or sub-samples of galaxies with similar bulgeto-disk ratios.

dense clouds
In our CCC model, the HII regions have been assumed to be identical.The best fit of the optical depth, τ cl ∼ 0.5, can be considered as a statistical average of different HII regions.However, in real galaxies, the optical depths of star forming regions are related to their evolution phase (McKee & Ostriker 2007).The star forming regions at their early stage, which are still embedded in molecular clouds, could be extremely dense and optically thick and make no contribution to the observed optical nebular emission lines at all.Not only that, these optically thick star forming regions would also block the stellar emission behind them and thus cause a significant fraction of dark area when galaxies are viewed from edge-on.
Our modelling is based on the observational constraints from the Balmer decrement of the optically thin HII re-gions and the reddening of stellar continuum in optical wavelengths.Therefore, these optically thick star forming regions will not directly bias our model results.Moreover, our results in Figure 6 also show that the CCC model predicted dust attenuation effects are in good consistence with observations once the aperture effects have been probably accounted.That is to say, at least inside the SDSS fiber aperture, there will not be many of these completely optically thick HII regions that could result in significant areas being completely obscured.However, we cannot exclude that, when the galaxies are completely edge-on, these dark areas will have a nonnegligible effect on the projected image.Indeed, the dust optical depth along the galactic plane predicted by our CCC model is slightly smaller than that of the X99 model, which is obtained from RT modelling and should be unaffected by these dark clouds.(Section 6.1).
The optically thick star forming regions have been properly accounted in T04 model with a parameter F (volume fraction of optically thick HII regions), which has a typical value of 0.2, but varying significantly in different galaxies.As the setup in the T04 model, these optically thick star forming regions would have a significant impact on the UV radiation from very young stellar objects and the corresponding transferred infrared emission.For future studies, these optically thick components need to be taken into account, especially when there are observational constraints from either the UV or IR.

Resolution of fiber spectroscopy
In this study, we mainly model the dust attenuation effects along the line of sight to galactic center since our observational constraints are mainly from the fiber spectroscopy of sample galaxies.As we have discussed in Section 4.4 and 5, the fiber aperture, although much smaller than typical sample galaxy size and can approximate the galaxy central region well for face-on galaxies, which also brings significant biases on interpreting the observed dust attenuation features of edge-on galaxies.
Here, we present another bias effect from the resolution of fiber aperture that has not been discussed yet.For edgeon galaxies, because the galaxy scale-height is much smaller than the fiber spectroscopy, the observed dust attenuation is not simply along the galactic plane but rather an integral along the vertical direction of the disk.As shown by the projected images of edge-on galaxies, the optical depth decreases significantly from the galactic plane to high latitudes, which thus makes the observed dust attenuation be systematically smaller than that completely along the galactic plane.In our modeling, for edge-on galaxies, we have only counted the dust optical depth and calculated the dust attenuation effect on the galactic plane.That is to say, for edge-on galaxies, our model fitting values are biased towards lower dust attenu-ation, which would bias the clumpy dust disk to be relatively thick.Indeed, as we have shown in Figure 7, the effective optical depth along the galactic plane predicted from our CCC model is smaller than that of the X99 model.However, the observational constraints of the CCC model are the reddening features of all inclinations (Figure 2, 90 θ bins), we therefore do not expect such a bias from edge-on galaxies only could significantly change our model results.

IR properties
In this study, our model has only investigated the dust absorption process without considering its emission.To better discuss and constrain the dust properties of galaxies, a radiative transfer model that takes into account IR emissions of dust is required, which however is beyond the scope of current work.Here, based on the framework and basic results of our CCC model, we give a brief outlook on the infrared emissions of our model galaxy.
In our CCC model, we have assumed that the dust in both clumps and ISM have the same properties (extinction curves), thus the interstellar radiation field (ISRF) plays the key role in dust temperature.Considering that the clumpy HII regions are radiated by central young stellar objects, we expect that the clumpy dust will absorb more UV photons and therefore constitute the warm dust that could be traced by MIR emissions.Indeed, in observation, the warm dust is spatially correlated with the molecular gas H 2 and star forming regions (Stevens et al. 2005;Hippelein et al. 2003).It is also shown that the scale-length of the MIR disk is similar to that of Hα (Vogler et al. 2005).
On the other hand, our CCC model shows that the clumpy dust is more extended than that of diffuse dust.If we assume that the clumpy dust in the HII region has averagely higher temperature than that of ISM dust, the disk will be more extended in the MIR than in the FIR, which is on the contrary to the observations that the disk scale-length of the IR emission increases with wavelength (Hippelein et al. 2003).One of the reasons for this contradiction is that our model does not contain any optically thick clumpy components, as discussed in Section 6.4.2, which has no impact on the model fitting since all our modelling constraints are in the optical wavelengths.These optically thick star-forming regions (e.g., molecular clouds) are one of the main contributors of cold dust (Planck Collaboration et al. 2011).If we take the IR properties of these optically thick star-forming regions into account and assume that they have the same geometrical distribution as the HII regions in CCC model, we would naturally get a more extended cold dust component.Another reason is that a uniform and low dust temperature assumption for the diffuse ISM dust is oversimplified.The temperature of the diffuse ISM dust is positively correlated with the intensity of ISRF, which decreases from the galactic center to the outer regions.
Therefore, the spatial distribution of the lower temperature component of diffuse ISM dust will also be biased to a larger scale-length.To fully quantify the IR emissions of our CCC model, we need to model the radiative transfer of both the clumpy dust (optically thin HII regions and optically thick molecular clouds) and the diffuse ISM dust in detail, which is beyond the scope of this work and requires further investigation in future works.

CONCLUSION
In this study, we have measured the dust reddening features from the fiber spectra of a sample of 33,065 MW-like disk galaxies in the SDSS, where the stellar reddening is determined using the full-spectrum SPS code STARLIGHT (Cid Fernandes et al. 2005), and the nebular emission line reddening is measured from the Balmer decrement.We explore the variation of these two different dust attenuation tracers as a function of disk inclination and then use it as a constraint to build geometric models for the dust attenuation of MW-like disk galaxies.
We find that the stellar attenuation shows a monotonic increase with disk inclination while the nebular attenuation does not.For highly inclined disks (θ > 75 • ), although the fluxes of emission lines (e.g.Hα) continue to decrease with increasing of disk inclination, the nebular attenuation (reddening) shows a saturation effect: E g (B − V) max ∼ 0.6.
For the model part, which is also the focus of this study, we find that a single uniform mixture model can generally reproduce the observed inclination dependence of the stellar attenuation, while a single screen model can only partly reconstruct the inclination dependence of the nebular attenuation for low inclination disks.Based on the results of these two simple models, we construct a new two-component dust geometry model, the Chocolate Chips Cookie (CCC) model.In the CCC model, the clumpy nebular regions are embedded in a diffuse stellar/ISM disk, like chocolate chips in cookies.By adopting a clumpy approximation and considering the off-axis effect of highly inclined disks, the CCC model successfully reproduces the observed inclination dependence of the stellar and nebular reddening simultaneously.Moreover, after proper accounting for the fiber aperture effect, the CCC model prediction on the inclination dependence of the dust attenuation effects of both Hα flux and r band magnitude are also in good consistence with observations.In addition to the observational properties from fiber spectroscopy, the global edge-on photometric properties predicted from our CCC model are also broadly consistent with the RT studies for nearby disk galaxies.
The best estimates of the geometric parameters of our model galaxies are as follows: the stellar disk scale-height to scale-length ratio h s /R s ∼ 0.19, nebular disk scale-height to scale-length ratio h g /R g ∼ 0.06, the ratio between the nebular disk scale-height to diffuse disk h g /h s ∼ 0.56, and the ratio between the nebular disk scale-length to the diffuse disk R g /R s ∼ 1.6.Moreover, we obtain model constraints on the optical depth of two dust disks.For the diffuse dust component, the CCC model predicts a line of sight dust reddening to high galactic latitude at "solar neighborhood" E(B − V) ∼ 0.02.This result indicates that statistically, our MW-like disk galaxies selected from the SDSS according to their stellar mass have similar geometry and dust properties to the MW as well.For the clumpy regions, we conclude that there are about 10 4 clumpy regions in our model galaxy if taking 30 pc as the clump size, and each region has an optical depth τ cl ∼ 0.50 in V band.These parameter estimates give a self-consistent inference that these clumpy regions can properly represent the HII regions in SFGs.
The CCC model framework also has some limitations.First, in order to reduce the number of free parameters, we have not considered the bulge component, which makes the scale-height of the stellar emission R s in the CCC model is an effective parameter, containing the bulge contribution.Secondly, our model has not considered the optically thick star forming regions, which may prevent the direct application of our model to radiative transfer studies involving UV and/or IR data.Finally, our CCC model is developed to study the dust attenuation process along a single line of sight, which poses difficulties in detailed comparison of our model results with the finite resolution observational data in a straightforward way.
The CCC model has various applications.For example, our model can be further extended and applied to the disk galaxies with other stellar masses.With such modeling, the structural parameters of the dust components of the galaxies with different masses can be compared and explored, enabling us a better understanding of the formation and evolution of the dust components in galaxies with different masses.Moreover, by assuming a physical prescription of the properties of the dust particles (e.g.Draine & Li 2007), our model can be applied on more detailed dust extinction and emission processes in different wave-lengths, which then can be compared with the results from numerical simulations using Monte-Carlo processes (e.g.Camps & Baes 2015).The CCC model is an analytical model based on several reasonable approximations.It simplifies the calculation of radiative transfer process when considering the star-dust geometry.Comparing with the common Monte-Carlo method for calculating the detailed radiative transfer process, the CCC model is much faster, and thus can be applied for large samples of galaxies conveniently.which is then approximated with Equation 21.
In specific, we have approximated the Binomial distribution B(N, p) with a Gaussian distribution N(N p, N p(1 − p)) under the assumptions of N >> 1 and p << 1.However, due to the extinction term exp(−kτ cl ) in Equation B1, a very large optical depth may lead to deviations from this approximation.In this appendix, we make detailed comparisons between the numerical calculation of the binomial distribution extinction (Equation B1) and our Gaussian model approximation(Equation 21).
For our numerical calculations, we examine the deviation between Equation 21and Equation B1 in terms of 1 − τcl /τ Bio .There are 3 parameters in Equation B1: number of clumps N, covering factor p and effective optical depth τcl of individual clump.We show the deviation as functions of different sets of N, p, τcl in three different panels of Figure B1.In top, middle and bottom panel, we set n = 10, 000, p = 10 −4 , τcl = 0.67 (the best estimates of the CCC model, τcl = 4/3τ cl ) and let the other two parameters being free respectively.The ranges of these free parameters are set as follows: 0 to 0.01 for p, 100 to 100,000 for N and 0 to 1 for τcl .
As can be seen, inside the parameter range, the deviation is almost independent of N and p, but significantly dependent on τcl .For the best estimate of the CCC model, τcl ∼ 0.67, the deviation is about 8%.However, we note that when τcl ∼ 1.0, the deviation of the approximation Equation 21 can be up to 20%.Moreover, this bias is a negative bias, i.e. it underestimates the effective optical depth of the clumps.

C. DUST ATTENUATION OF CLUMPY CLOUDS: DIFFERENCE FROM CONTINUOUS DISTRIBUTION
We present a quantitative comparison between the dust attenuation of two different models: the clumpy distributed dust and continuously distributed dust.In specific, we consider two models that both have the same amount of dust inside a given column V = S * l (much larger than the size of dust cloud itself), but one is continuously distributed and the other is clumpy distributed.We illustrate the clumpy distribution and the continuous distribution in Figure C1.For the clumpy dust distribution, if there are N dust clouds uniformly distributed inside V, the total amount of dust is where R cl is the radius of the dust cloud, τ cl is the optical depth of a emission source at the center of the cloud, and κ is the dust extinction coefficient.In this case, as we have shown in Section 4.1 (Equation 22), the line of sight dust optical depths for a region outside cloud and a region inside cloud can be approximated by As can be seen, there are two major differences between these two models.On the one hand, for the clumpy distributed dust, the dust attenuation of the inside cloud regions and outside cloud regions are different with a term of τ cl , while the continuous dust does not have this effect.On the other hand, for the diffuse ISM region, the clumpy distributed dust has a systematically smaller effective dust extinction than the continuously distributed dust, with a factor of (2 − τ)/2, which is also known as 'mega-grain' effect (Városi & Dwek 1999).
Combining these two effects, the attenuation of the emission lines by the clumpy dust can be either greater or smaller than the case of continuous dust, depending on the column number density of the clumpy clouds along the line of sight.For example, if we take τ cl = 0.5 and assume n g πR 2 cl = 0.5, corresponding to the face-on view of a typical disk galaxy, the resulted τl,cl,in is larger than τl,cont with a factor of 0.38.If we take n g πR 2 cl = 4 (approximating edge-on view), τl,cl,in , on the contrary, is smaller than τl,cont with a factor of 0.5.

Figure 1 .
Figure 1.(a) The histograms of disk inclinations.The dotted line shows all the SFGs in MGS, while the sold line represents SFGs at z < 0.1.Both histograms are normalized to unit area.(b)The density map of the stellar mass and disk inclination of the SFGs at z < 0.1, where the black dots represent median stellar mass at different inclinations.

Figure 2 .
Figure 2. The distributions of nebular color excess E g and stellar color excess E s as functions of disk inclination θ.Panel (a): Joint distribution map of E g and E s of our sample galaxies, where θ is indicated by the color.The black line shows the regression line with forced zero intercept, E g = 0.4E s .Panel (b): E s as function θ, where the median of E s distribution at each θ bin are shown as square dots connected with a solid line and the 16 and 84 percentiles are indicated by two solid curves.Panel (c): E g as function θ, the same structure as Panel (b).

Figure 3 .
Figure 3. Modelling of the inclination dependence of E s and E g with simple models.Panel (a): uniform mixture model for E s .The small dots show the observed median E s − θ relation, while the solid line shows the best model fitting presented in Section 3.1.Panel (b): screen model for E g .The small dots show the observed median E g − θ relation, while the solid line shows the best model fitting presented in Section 3.2.The schematic diagram of each model and the values of the best model parameters are inserted in each panel.

Figure 4 .
Figure 4. Modeling of the E g − θ and E s − θ relations with the Chocolate Chip Cookie (CCC) model.Panel (a) : schematic map of the CCC model, where R g , h g , R s , h s are the scale-length and scale-height of the clumpy nebular disk and diffuse ISM disk respectively and τ cl is the optical depth of each individual clumps.The dashed arrow line shows the off-axis effect of a line of sight, where the stellar emission has not been extincted by the nebular disk (see Section 4.4 for detail).Panel (b): the corner plot of the MCMC fitting for the six free model parameters (Table1) of the CCC model, where the best fittings of E g − θ and E s − θ relations are inserted at the top right corner of this panel.

Figure 5 .
Figure 5.The absolute r band fiber magnitude(top panel) and H α flux(bottom panel) as function of disk inclination θ for our MWlike SFGs, where the number density of the sample galaxies in the parameter space are color coded.The median values at different inclinations are shown as solid dots in each panel.

Figure 6 .
Figure6.The inclination dependent Hα flux and r band fiber magnitude explained by the CCC model and aperture effect.Top panel: the observed median fiber magnitude in θ bins (dots), the attenuated central surface brightness predicted by the CCC model (dashed line), the aperture corrected fiber magnitude after applying aperture corrections for the model galaxy with R a /h g = 10(solid line).Middle panel: the observed median Hα flux in θ bins (dots), the attenuated central nebular line fluxes predicted by the CCC model (dashed line), the aperture corrected nebular line fluxes after applying aperture corrections for three model galaxies with different physical sizes (solid lines, different colors represent different h g values).Bottom panel: the projection of a double-exponential disk into a fiber with radius R a (2.2 kpc) at different inclinations.The three solid lines show the projected area S as function of disk inclinations for three model galaxies, which are all set to unit value at θ = 0 • .The three schematic icons show the projection of a model galaxy with scale-length 3.2 kpc and scale-height 0.2 kpc into the fiber aperture (2.2 kpc) for face-on (case I), edge-on (case II) and critical (case III) inclinations respectively.

Figure 7 .
Figure7.The edge-on projection images, effective optical depth and surface brightness profiles of the CCC, X99 and T04 model galaxies.The top three image panels: the edge-on projection images of the CCC,X99 and T04 models.The middle plot panel: the effective optical depth profiles along the galactic plane of the three models.The bottom panel: the surface brightness profiles along the galactic plane of the three models.

Figure 8 .
Figure 8. Extinction and attenuation curves of the CCC model and the Calzetti law represented by A(λ)A V .Top panel: the extinction curve (Equation9) used in the CCC model (red curve) and the Calzetti attenuation curve (blue curve).Bottom panel: the attenuation curves derived from the CCC model for the edge-on case (dashed, θ = 90 • , R V = 7.0), face-on case (dotted, θ = 0 • , R V = 3.7) and median case (solid,θ = 60 • , R V = 4.1), the blue curve also shows the Calzetti attenuation curve for comparison.
Figure B1.The consistence of the model predicted optical depth of clumpy regions between the Gaussian approximation and numerical calculation of binomial distributions for different parameter sets, which are color-coded by 1 − τcl /τ Bio in each panel.Panel (a): τcl = 0.67, p ranges from 0 to 0.01 and n ranges from 100 to 100,000; Panel (b): p = 10 −4 , n ranges from 10 to 100,000 and τcl ranges from 0 to 1; Panel (c): n = 10, 000, p ranges from 0 to 0.01 and τcl ranges from 0 to 1.

Figure C1 .
Figure C1.Cartoon of clumpy dust distribution(left) and continuous dust distribution(right) for the same amount of dust in the same volume V = S * l.
n g = N/S is the column number density of the dust clouds along the line of sight.For the continuous dust distribution, the optical depth along the same line of sight can be expressed as τl,cont = M tot S κ = n g πR 2 cl τcl .(C4)