Monte-Carlo radiative models of the diffuse Galactic light (DGL) in our Galaxy are calculated using the dust radiative transfer code MoCafe, which is three-dimensional and takes full account of multiple scattering. The code is recently updated to use a fast voxel traversal algorithm, which has dramatically increased the computing speed. The radiative transfer models are calculated with the generally accepted dust scale-height of 0.1 kpc. The stellar scale-heights are assumed to be 0.1 or 0.35 kpc, appropriate for far-ultraviolet (FUV) and optical wavelengths, respectively. The face-on optical depth, measured perpendicular to the Galactic plane, is also varied from 0.2 to 0.6, suitable to the optical to FUV wavelengths, respectively. We find that the DGL at high Galactic latitudes is mostly due to backward or large-angle scattering of starlight originating from the local stars within a radial distance of r < 0.5 kpc from the Earth. On the other hand, the DGL measured in the Galactic plane is mostly due to stars at a distance range that corresponds to an optical depth of ≈ 1 measured from the Earth. Therefore, the low-latitude DGL at the FUV wavelength band would be mostly caused by the stars located at a distance of r 0.5 kpc and the optical DGL near the Galactic plane mainly originates from stars within a distance range of 1 r 2 kpc. We also calculate the radiative transfer models in a clumpy two-phase medium. The clumpy two-phase models provide lower intensities at high Galactic latitudes compared to the uniform density models, because of the lower effective optical depth in clumpy media. However, no significant difference in the intensity at the Galactic plane is found. At wavelengths longer than the Lyman limit, dust radiative transfer processes play a crucial role in various astronomical systems, ranging from reflection nebulae with a single or a few stars (e.g., Witt 1977 ) to galaxies composed of stars, gas and dust (e.g., Bianchi 2008 ). Many different approaches have been developed to study the radiative transfer problems in dusty media (e.g., Baes & Dejonghe 2001 ). One example is to use an expansion in spherical harmonics to solve the radiative transfer equation ( Frannery et al. 1980 ; Roberge 1983 ). Monte-Carlo techniques are the most flexible and allow arbitrary geometries for stars and dust (e.g., Baes & Dejonghe 2001 ; Seon 2009c ). The main disadvantages of Monte-Carlo methods are that they are relatively slow and numerically demanding compared to other techniques. However, several clever Monte-Carlo radiative transfer algorithms have been developed and the exponential growth of computing speed has recently allowed to handle realistic three-dimensional geometries. Studies of the scattered light by interstellar dust are unique since they allow the absorption and scattering properties of dust grains to be separated. There are three different types of astrophysical objects suitable for studies of the dust scattering properties: reflection nebulae, dark clouds, and the diffuse Galactic light (DGL; e.g., Gordon 2004 ). The DGL results from coherent scattering of Galactic starlight by diffuse interstellar dust grains. There have been a variety of attempts to derive the scattering properties of the interstellar dust grains from observations of the DGL at different wavelengths (see Witt 1990 for a historical review). Brandt & Draine (2012) have investigated the spectrum of the DGL at optical wavelengths by using data from the Sloan Digital Sky Survey (SDSS; York et al. 2000 ). The DGL at far-ultraviolet (FUV) wavelengths has been observed through space missions. Witt et al. (1997) analyzed the DGL data at ~156 nm obtained with the Far Ultraviolet Space Telescope , flown as part of NASA’s ATLAS-1 space shuttle mission. Schiminovich et al. (2001) reported the results of a sounding rocket experiment, the Narrowband Ultraviolet Imaging Experiment for Wide-Field Surveys (NUVIEWS), designed to measure the DGL at ~174 nm. More recently, Seon et al. (2011a , b) presented the spectral properties of the DGL at FUV observed over ~80% of the sky with the Far-ultraviolet IMaging Spectrograph (FIMS; also known as the Spectroscopy of Plasma Evolution from Astrophysical Radiation , SPEAR) onboard of the first Korean science and technology satellite STSAT-1. The DGL observed with the Galaxy Evolution Explorer mission was presented in Murthy et al. (2010) , and Murthy (2014a , b) . Radiative transfer models for the DGL have been developed by Witt (1968) , van de Hulst & de Jong (1969) , Jura (1979) , Witt et al. (1997) , and Schiminovich et al. (2001) . In particular, van de Hulst & de Jong (1969) provided a numerical method to solve the radiative transfer equation for a plane-parallel geometry. However, they calculated the model for a slab geometry, in which both dust and stars are distributed uniformly over the same slab. Seon (2009b) has extended the iteration method of van de Hulst & de Jong (1969) to the case of exponentially distributed stars and dust. Witt et al. (1997) employed a Monte-Carlo model that utilized a realistic radiation field based on a measured star catalog and a cloud mass spectrum inferred from the observations of ISM. They calculated the scattered light by placing clouds randomly along each line of sight. A similar approach was also adopted in Schiminovich et al. (2001) . However, none of the previous studies has investigated in detail how far or close the stars mostly responsible to the DGL are located from the Earth. This is crucial when we try to make a realistic model of the DGL, because we need to determine to what extent the stellar and dust distributions should be appropriately modeled. It is also well-known that the ISM density distribution is not homogeneous, but highly structured or clumpy. Therefore, it may be also worthwhile to investigate the effect of a clumpy density structure on the DGL. In this paper, we use a fully three-dimensional MonteCarlo radiative transfer code to calculate the DGL and investigate the effect of inhomogeneity in the context of a two-phase clumpy medium. In Section 2, the MonteCarlo radiative transfer code is described. The detailed model results are given in Section 3. Section 4 summarizes the results. The radiative transfer models of the DGL are calculated using the three-dimensional Monte-Carlo simulation code MoCafe ( Mo nte- Ca rlo Radiative Trans fe r), which has been used in Lee et al. (2008) , Seon (2009a) , Jo et al. (2012) , Seon & Witt (2012 , 2013) , and Seon et al. (2014) . Detailed algorithms of the simulation code are described in Seon (2009c , 2010) . The basic MonteCarlo algorithms have been described in detail also by many authors ( Gordon et al. 2001 ; Baes et al. 2011 ; Steinacker et al. 2013 ). Here, we only recall the major techniques implemented in the code, such as the “first forced scattering” ( Cashwell & Everrett 1959 ) and the “peel-off” technique ( Yusef-Zadeh et al. 1984 ), and recent updates of the code. The interested reader is referred to these papers for more details. The concept of a weight ( w ) for every photon packet is used to increase the efficiency of a Monte-Carlo simulation. The first scattering of the photon is forced to insure that every photon contributes to the scattered light (the “first forced scattering” technique). The weight of the photon package is then decreased accordingly to correct for the escaped fraction of the photon packet. As in most radiative transfer codes, subsequent scattering events are not forced. While each photon experiences multiple scattering in random directions as it propagates through the dust medium, the code calculates which fraction of the photon would arrive at the location of the observer from each scattering (the “peeling-off” technique). The fraction is summed up for the prediction of the final output image. These two techniques greatly increase the signal-to-noise compared to the simple Monte-Carlo simulations, where only photon packages arriving at the observer or leaving the system are recorded for the output. Photon packets are isotropically emitted from a randomly selected location, according to a given stellar distribution. The observer is assumed to be located at the center of the Galactic plane and the all-sky intensity map of scattered light is constructed by binning the “peeled-off” photons using the Hierarchical Equal Area isoLatitude Pixelization (HEALPix) tessellation scheme ( Górski et al. 2005 ) with ~ 0.9° pixels (corresponding to the resolution parameter N _side = 64 and the number of pixels N _pix = 49,152). The all-sky intensity map is then median-filtered. The filtered map is rebinned to larger pixels with a ~ 1.8° resolution and a resolution parameter N _side = 32 to reduce Poisson noise, which is intrinsic in Monte-Carlo methods. Photon directions after each scattering event are randomly generated to follow the Henyey-Greenstein scattering phase function ( Henyey & Greenstein 1941 ) where θ is the scattering angle measured from the direction of the incident photon and g ≡ 〈cos θ 〉 is the phase-function asymmetry factor. The photon weight w is reduced by a factor of albedo a after each scattering event. Recently, MoCafe has been updated to employ a fast voxel traversal algorithm for ray tracing, developed for the computer graphics rendering by Amanatides & Woo (1987) . Tracking photon packets is the most timeconsuming part of Monte-Carlo radiative transfer simulations, and thus a fast ray-tracing algorithm is essential ( Kurosawa & Hillier 2001 ). The voxel traversal algorithm was indeed very fast, and dramatically increased the speed of radiative transfer calculations by a factor of three compared to the previous ray-traversal algorithm described in Seon (2009c) . Moreover, the algorithm provides a method to accurately track the paths of photon packets, which is crucial for simulations in highly inhomogeneous media, even with single precision floating-point numbers. For the radiative transfer model of the DGL, we need to specify the three-dimensional spatial distributions of dust and photon sources. In the present study, we examine two types of dust media: a smoothly varying, homogeneous medium, and clumpy two-phase medium. In the case of homogeneous dust distribution, we assume that interstellar dust is plane-parallel and exponentially distributed perpendicular to the Galactic plane. The dust density ρ ( z ) depends only on the vertical coordinate z and is given by where ρ ₀ and z _d are the dust density in the Galactic plane ( z = 0) and the scale-height of dust, respectively. The face-on optical depth, i.e., the total optical depth of a dust layer when our Galaxy is seen face-on, is then τ ₀ = 2 κρ ₀ z _d . Here, κ is the dust extinction coefficient. Photon sources are also assumed to be plane-parallel and distributed exponentially with a scale-height z_s . The stellar luminosity distribution is then given by where L ₀ = 1 is the stellar luminosity per square kpc at z = 0. The results presented here can be scaled up to an arbitrary value of the stellar luminosity per unit area. A two-phase clumpy medium consists of high-density clumps and a low-density inter-clump medium. The density structure of a clumpy medium can be characterized by the method of Witt & Gordon (1996) . In the algorithm, each voxel is assigned randomly to be either a high- or low-density phase. The statistical properties of a two-phase medium are determined by a volume filling factor, f , of high-density phase (clumps) and a density ratio, ρ_ℓ / ρ_h , of the low-density phase (inter-clump medium) to high-density clumps. The probability of any cell being in a high-density state or a low-density state is randomly determined, based on the filling factor f . The probability of each voxel being in high- or low-density phase is independent of the density phase of adjacent cells. Therefore, the medium is statistically homogeneous. In the present study, even in the cases of clumpy media, we assume that the “mean” dust density 〈 ρ 〉 obeys the exponential distribution in Equation (2). To ensure that the mean density of dust follows the exponential distribution as for the homogeneous model, we first assign the mean density distribution given by Equation (2) to each individual voxel. For a given mean density 〈 ρ 〉 at Galactic height z , the densities allocated to individual cells of high-density and low-density phases are determined by respectively. In this way, we can generate a clumpy, two-phase medium, which still statistically follows an exponential distribution. We considered the filling factors of f = 0.05, 0.1, and 0.5, that cover the extent of physical environments ranging from rare high-density cells in an extended low-density medium, till the case in which the high-density medium forms an interconnected structure, with voids filled by a low-density medium. The density ratios of ρ_ℓ / ρ_h = 0.01, 0.1, and 0.2 were considered. However, only the results for the extreme values of f = 0.05 and 0.5 and ρ_ℓ / ρ_h = 0.01 and 0.2 were presented (for clarity) in this paper. The models obtained with the intermediate values were found to have the properties easily inferred from other models. The physical size of each voxel bin was assumed to be 10 pc. However, since the density of a homogeneous plane-parallel model depends only on the Galactic height, a single bin along the x and y axes is adopted. In the present study, the dust scale-height was assumed to be 0.1 kpc and the physical size of the simulation box in the z direction ±1.005 kpc. The number of voxels for a homogeneous medium model is thus N_xN_yN_z = 1 × 1 × 201, regardless of the physical sizes of the grid along x and y . For the clumpy medium models, the number of bins along x and y was, however, increased up to 1800, proportional to the physical size of the simulation box along the axes. Thus, the number of voxels was up to N_xN_yN_z = (1800) ² × 201, depending on the physical size. The side lengths of the simulation box along x and y axes were always kept the same, e.g., N_x = N_y . The physical size of the simulation box in the x and y directions was increased from ±0.2 to 9 kpc, to investigate the effect of increasing the size of the simulation box. In the homogeneous models, we used 10 ¹⁰ photon packets. On the other hand, only 10 ⁹ photon packets were used for clumpy models because a huge number of voxels was needed in those models. We note that the computation time is linearly proportional to the bin number of box sides. Five random realizations were made for each set of parameters f and ρ_ℓ / ρ_h , and then the results were averaged to show the general properties for the clumpy medium, instead of a particular result specific to a single realization. We adopt the spatial distribution of the Milky Way dust layer found by Misiriotis et al. (2006) . They used the COBE/DIRBE maps and the COBE/FIRAS spectra to constrain the spatial distribution of dust, gas, and stars in our Galaxy. The dust scale-height is thus assumed to be z _d = 0.1 kpc. In their model, the faceon optical depth of our Galaxy in the V-band is about 0.2 at the location of the Sun. In the FUV wavelength band, the model gives a face-on optical depth of ~ 0.6. We, therefore, calculated the DGL models by varying the face-on optical depth from 0.2 to 0.6. The scaleheights of early- and late-type stars are z_d ~ 0.1 and ~ 0.3−0.4 kpc, respectively. We assume a stellar scaleheight of 0.1 or 0.35 kpc. We note that the former value is the same as that of the dust layer, while the latter is higher than that of the dust. In the models with z_s = 0.35 kpc, the starlight originating from above the dust layer would be forward-scattered and the resulting scattered light could be measured at high Galactic latitudes. However, this is not the case of the models with z_s = 0.1 kpc. Figure 1 shows all-sky intensity maps for a homogeneous medium model with z _s = 0.1 kpc and τ ₀ = 0.2. Note the Poisson noise signals in the left panel, mainly caused by stellar sources that happened to be located very close to the observer. This is a limitation of our algorithm, for which the locations of photon packets were uniformly selected within the simulation box. Even with a huge number of photon packets used in the present study, only a few of them were created to originate from very close points, producing strong signals at the points because of their proximity to the observer. However, these photon packets are supposed to represent the whole photons that should be spread all over the directions within the distance range instead of being concentrated only on a few points. The best way to avoid this random noise would be to generate more photon packets in a closer distance, but with a lower luminosity to compensate the increased number of photon packets at the distance range ( Seon 2010 ). In the present study, instead of applying this technique, a median-filter was applied to reduce the noise and the resulting map was rebinned to a lower resolution map, as described in the previous Section. The right panel of Figure 1 shows the median-filtered and rebinned map. The intensity profile versus the Galactic latitude was not significantly altered after the median-filtering. It is worthwhile to examine at what distance range the stars would give rise to most of the DGL intensity. In other words, the dependence of the scattered light intensity, measured on Earth, at the distance of the stellar source would be interesting. The DGL models were therefore calculated assuming that the stars are only confined within certain limits of the radial distance from the observer. The intensity profiles of the scattered light resulting from stars located at various distance ranges, as functions of Galactic latitude, are shown in the first and third columns of Figure 2 . The profiles were obtained by averaging the intensities at each Galactic latitude bin. In the figures, the stellar distance range was varied as follows: 0 < r < 0.2 kpc, 0.2 < r < 0.5 kpc, 0.5 < r < 1 kpc, 1 < r < 2 kpc, 2 < r < 3 kpc, 3 < r < 4 kpc, 4 < r < 5 kpc, 5 < r < 6 kpc, 6 < r < 7 kpc, 7 < r < 8 kpc, and 8 < r < 9 kpc. Here, r denotes the radial distance from the z axis in the cylindrical coordinate system. The intensity profiles normalized to the total intensity profile are also shown in the second and fourth columns. Colors of the curves were altered for clarity. The blue curve denotes the total intensity profile, obtained by adding the profiles calculated for the individual distance ranges. The left two panels (first and second columns) are the results for the models with the stellar scale-height of z_s = 0.1 kpc, while the right two panels (third and fourth columns) are for the models with z_s = 0.35 kpc. The face-on optical depth was assumed to be τ ₀ = 0.2 and 0.6 for the top and bottom panels, respectively. The most interesting result in Figure 2 is that the DGL at high-Galactic latitudes mainly originates from local stars close to the observer. More than half of the DGL at the Galactic pole (| b | = 90°) is due to local stars within a distance of ~ 0.5 kpc. The contribution from local stars to the DGL at high latitudes tends to be higher in models with smaller stellar scaleheight of z_s = 0.1 kpc than in cases with z_s = 0.35 kpc. The contribution from local stars within a distance of 0.2 kpc becomes dominant above the Galactic latitude of | b | ≈ 30°, when the stellar scale-height is z_s = 0.1 kpc. When z_s = 0.35 kpc, those stars start to dominate the DGL at a slightly higher latitude of | b | ≈ 50°. These results indicate that the DGL at highGalactic latitudes is mostly due to backward scattering or large-angle scattering of the local starlight. Since the stellar scale-height is smaller than or compatible with that of dust, the only chance that the scattered light can be detected at high Galactic latitudes is through backward or large-angle scattering. On the other hand, if the stars are present above the dust layer, forward scattering of the starlight originating from above the dust layer can also give rise to scattered light at high latitudes. Note also that the intensity profile originating from stars within 0.2 kpc is relatively flat, especially in models with z_s = 0.35 kpc, compared to profiles from distant stars. In addition, most of the DGL measured at the Galactic plane appears to come roughly from a distance range where stellar photons are scattered off only once by dust grains between the source location and the observer. This corresponds to a distance with an optical depth of τ ≈ 1 as measured from stars to the observer. Photons originating from a closer distance in the Galactic plane would have less chance to be scattered towards and measured at the observer. On the other hand, photons from a greater distance will be strongly attenuated. Therefore, photons from a distance corresponding to an optical depth of τ ≈ 1 are the most dominant source of DGL at the Galactic plane. The unit optical depth at the Galactic plane ( z = 0) corresponds to a distance of 1 kpc and 0.33 kpc for the face-on optical depth of τ ₀ = 0.2 and 0.6, respectively. This property can be confirmed in Figure 2 where it is apparent that stars in the distance range of 0.5 − 2 kpc are a predominant source of DGL models with τ ₀ = 0.2 at the Galactic plane, while stars in the distance range of 0.2 − 1 kpc dominate DGL models with τ ₀ = 0.6. It is also found from the model with τ ₀ = 0.2 and z_s = 0.1 in Figure 2 that all the stars within a distance of at least up to r = 5 ~ 6 kpc should be included in calculating the DGL model to achieve a result with an accuracy of ~ 2% or better at the Galactic plane. This can be further confirmed by comparing the DGL models that were sequentially calculated with increasing outer radial boundary for the stellar source. Figure 3 compares the results for a model with τ ₀ = 0.2 and z_s = 0.1 kpc, as the radial boundary of the stellar distribution is increased from r = 0.5 kpc to r = 9 kpc in various intervals. In the calculations, the inner radial boundary was always r = 0 kpc, unlike in Figure 2 . It is clear that, as the outer boundary increases, the resulting intensity gradually increases compared to the result of the previous step until the calculated DGL intensity starts to converge at ~ 5−6 kpc within a level, comparable to the Monte-Carlo noise. The radial boundary of the stellar distribution required for the convergence is clearly smaller for the models with a higher face-on optical depth (i.e., for those with τ ₀ = 0.6 as in Figure 2 ) because of the stronger attenuation at the higher optical depth. We also note, as can be seen in Figure 2 , that a model with a higher stellar scale-height requires a slightly larger outer boundary to obtain a similar accuracy. This is because the starlight emitted from above the dust layer has a higher chance to be scattered towards the observer with less attenuation compared to the model with a lower stellar scale-height. In the following models, the physical size of the simulation box was assumed to be 6 kpc so that all the necessary starlight is taken into account even in the worst case. In all of the above models, we have assumed an albedo of a = 0.6. We now investigate the dependence of the intensity profile on albedo by varying the albedo from 0.1 to 1.0 in steps of 0.1. Photon packets with a higher albedo are less diluted after each scattering event. In an extreme case, photon packets with an albedo of a = 1 would be never attenuated, but they are scattered forever unless they escape out of the system through the boundaries. Therefore, as the albedo increases, photons have more chances to be multiply scattered and the multiplicity of scattering would be the lowest at a = 0.1. Thus, the intensity profiles calculated for various albedos were normalized to the profile with a = 0.1, as shown in Figure 4 , where the multiplicity of the scattering is the lowest. It is clear from Figure 4 that the intensity level is not linearly proportional to the albedo. Namely, the normalized intensity profiles are not equally spaced, but the spacing increases with the albedo. This implies the increasing importance of multiple scattering with albedo. The multiple scattering process, especially at the Galactic plane, obviously becomes more significant as the optical depth τ ₀ increases. Therefore, the spacing between the normalized intensity profiles increases more rapidly in models with a higher τ ₀ compared to models with lower τ ₀ . This effect is more significant at the Galactic plane, where most of the dust is located. Note also in Figure 4 that the normalized intensity pro- files for models with a higher phase function asymmetry factor (i.e., g = 0.7 in Figure 7 ) is more strongly enhanced at the Galactic plane than the profiles of models with a lower asymmetry factor (i.e., g = 0.5). As the phase function asymmetry factor g increases, the scattering process becomes more strongly forward-directed. Therefore, photon packets with a higher g , traveling in the Galactic plane, will remain within the plane where most of dust resides, and experience more multiple scattering than photons packets with lower g . This is the main reason why the normalized intensity profiles for a higher g shows stronger peaks at the Galactic plane than those with a lower g . Figure 5 compares the intensity profiles for various models. In general, models with a higher g value show slightly higher peaks at the Galactic plane and much weaker intensity at high latitudes, compared to models with a lower g . This is because photon packets traveling in the Galactic plane, where most photons are produced, will remain in the Galactic plane and only a small fraction of those packets is scattered into high latitudes. We assumed a constant stellar luminosity per kpc ² along the Galactic plane and the total photons were spread out over the stellar scale-height. Therefore, models with a higher stellar scale-height (denoted by circles in Figure 5 ) appear to have a weaker and broader peak at the Galactic plane than those with a lower scale-height (denoted by solid lines). However, at high latitudes, the model with a higher scale-height results in a higher intensity because the starlight residing above the dust layer can be easily forward-scattered towards the observer. This effect is particularly noticeable in models with a higher asymmetric phase factor, i.e., models with g = 0.7 in the right panel of Figure 5 . On the one hand, a higher face-on optical depth would result in a stronger attenuation at the Galactic plane. On the other hand, more scattering would occur at high Galactic latitudes where the absorption effect is weak. Therefore, the intensity of the scattered light in the model with a higher optical depth would be lower at the Galactic plane, but higher at high latitudes, compared to the model with a lower optical depth. This can be confirmed by comparing the models denoted with different colors (black, red, and blue curves for τ ₀ = 0.2, 0.3, and 0.6, respectively) in Figure 5 . Figure 6 shows a few examples of all-sky intensity maps for the clumpy two-phase medium models, in which the filling factor of clumps f is 0.05 or 0.5 and the density ratio ρ_ℓ / ρ_h is 0.01 or 0.2. The stellar scale-height, the albedo, and the phase function asymmetry factor were assumed to be z_s = 0.1 kpc, a = 0.6, and g = 0.5, respectively. Clumpy structures are visible, especially, in models in which the filling factor of high density clumps f is smaller and the density ratio of clumps to inter-clump medium ρ_ℓ / ρ_h is lower. A higher f value yields a larger fraction of the lines of sight being covered by high density clumps. A higher ρ_ℓ / ρ_h value means no significant density contrast. Therefore, the higher f and ρ_ℓ / ρ_h values give a relatively smoother map, as shown in Figure 6 . By comparing Figure 6 with Figure 1 , it is noticeable that the scattered intensity at high-Galactic latitudes in clumpy density models is generally weaker than in uniform density models. It is a well known fact that the effective opacity of a clumpy dust distribution is less than that of a homogeneous distribution of the same dust mass ( Witt & Gordon 1996 ). The weakness of the scattered light in high-Galactic latitudes in clumpy models is due to a lower effective opacity in a clumpy dust medium. Both the filling factor f of the clumps and the density contrast between clumps and the inter-clump medium are important parameters to determine the effective optical depth. In the limit in which the density ratio ρ_ℓ / ρ_h is near unity, the effective optical depth is mainly determined by the opacity of the inter-clump medium ( Witt & Gordon 1996 ). On the other hand, in the limit of large density contrasts between the clumps and the inter-clump medium, the effective optical depth of the dusty medium is essentially determined by the clump filling factor f , as long as the optical depth per clump is substantial. The dependence of the intensity profile on the density contrast ρ_ℓ / ρ_h and the clump filling factor f is shown in Figure 7 . Five realizations of the clumpy density distribution were produced for each combination of the parameters f and ρ_ℓ / ρ_h , and then the intensity profiles calculated with these five realizations were averaged to yield an average intensity profile for each set of the parameters. In the figure, red and blue colors denote the clumpy models with f = 0.05 and 0.50, respectively. The models with ρ_ℓ / ρ_h = 0.01 and 0.02 are shown with lines and filled circles, respectively. The uniform density models are also shown with black lines, for comparison. A higher filling factor of clumps (blue color in Figure 7 ) yields higher intensities at high latitudes compared to the models with a lower filling factor (red color), because the high density clumps at high altitude give more scattering. A higher density ratio ρ_ℓ / ρ_h (filled circles in Figure 7 ) also gives higher intensities at high Galactic latitudes than models with a lower density ratio (lines). This is because a lower density contrast (a higher ρ_ℓ / ρ_h value) gives less departure from the uniform density model. In the Galactic plane, the DGL intensity for a clumpy medium is slightly suppressed, compared to that of a uniform medium. However, the suppression is not significant, regardless of the values of ρ_ℓ / ρ_h and f . If a clumpy system is not large enough that the photon packets starting from the system boundary are not fully extinguished after being traveled to the observer, then the intensity profile at the Galactic plane would be significantly suppressed. However, we assumed a large enough system so that the observer in a clumpy medium can detect the photons originating from a much farther distance than the observer in a uniform density medium. Namely, the photons from a far distance in a clumpy medium, which would have been completely extinguished in a uniform medium, can arrive to the observer. The contribution from a far distance compensates less scattering in the clumpy medium. Thus, no significant suppression was found in the Galactic plane. The dust radiative transfer model of the DGL was investigated using the Monte-Carlo radiative transfer code MoCafe. We found that the intensity of the scattered light at high Galactic latitudes are caused by a large angle scattering of relatively nearby stars. On the other hand, the DGL at the Galactic plane is mostly due to stars in a distance range that corresponds to an optical depth of ≈ 1 as measured from the observer. The DGL models in the clumpy density medium were found to provide lower intensity at high Galactic latitudes, compared to the uniform density medium. However, the Poisson random noise is unavoidable when the locations of photon packets are uniformly determined and the weights of all the photon packets are set equal. The most naive method to avoid noise would be to increase the number of photon packets. However, a more efficient algorithm can be elaborated to reduce noise. A promising algorithm is to generate more photon packets at a closer distance but with a lower luminosity (or a weight). In the present study, the stars were assumed to be distributed in a plane-parallel manner with an exponential functional form. It might be reasonable to describe distant stars as being distributed plane-parallel. However, the fact that the DGL at high-Galactic latitudes are predominantly caused by the relatively nearby stars implies that the DGL would strongly depend on the detailed three-dimensional distributions of local stars and dust. The local effect becomes even more important at a shorter wavelength band, where the optical depth is higher and thus the starlight from distant stars does not contribute to the DGL. Hence, it would be crucial, especially in short wavelengths such as FUV, to use a more realistic model for the spatial distributions of stars and dust. In addition, the Henyey-Greenstein function adopted in the present study may not be appropriate for the scattering phase function in FUV wavelengths. Draine (2003) found that the Henyey-Greenstein function underestimates the forward scattering intensity by a factor of ~ 2 for cos θ 0.98 ( θ 10°). Moreover, although he does not mention in the paper, his Figure 7 reveals that the backward scattering intensity is also underestimated in the Henyey-Greenstein function for cos θ −0.7 ( θ 134°). Since the DGL at high Galactic latitudes is largely caused by backward or largeangle scattering of the radiation field of nearby stars, it would be necessary to take into account of a correct phase function in calculating a model of the DGL. In summary, we may not only need to develop a detailed three-dimensional distribution of the local stars and dust, but also use a correct, numerical functional form of the scattering phase function, instead of using the widely adopted, approximate Henyey-Greenstein function. We plan to develop a more realistic, three-dimensional model of the stars and dust and compare the DGL model with observational data. We will also use the numerical phase functions that are calculated with the Mie theory, instead of the approximate Henyey-Greenstein function. It would be also interesting to adopt a more realistic, hierarchical density distribution than the simple two-phase medium, as in Seon & Witt (2013) .