Sub-kpc scale gas density histogram of the Galactic molecular gas: a new statistical method to characterise galactic-scale gas structures

To understand physical properties of the interstellar medium (ISM) on various scales, we investigate it at parsec resolution on the kiloparsec scale. Here, we report on the sub-kpc scale Gas Density Histogram (GDH) of the Milky Way. The GDH is a density probability distribution function (PDF) of the gas volume density. Using this method, we are free from an identification of individual molecular clouds and their spatial structures. We use survey data of $^{12}$CO and $^{13}$CO ($J$=1-0) emission in the Galactic plane ($l = 10^{\circ}$-$50^{\circ}$) obtained as a part of the FOREST Unbiased Galactic plane Imaging survey with the Nobeyama 45m telescope (FUGIN). We make a GDH for every channel map of $2^{\circ} \times 2^{\circ}$ area, including the blank sky component, and without setting cloud boundaries. This is a different approach from previous works for molecular clouds. The GDH fits well to a single or double log-normal distribution, which we name the low-density log-normal (L-LN) and high-density log-normal (H-LN) components, respectively. The multi-log-normal components suggest that the L-LN and H-LN components originate from two different stages of structure formation in the ISM. Moreover, we find that both the volume ratios of H-LN components to total ($f_{\mathrm{H}}$) and the width of the L-LN along the gas density axis ($\sigma_{\rm{L}}$) show coherent structure in the Galactic-plane longitude-velocity diagram. It is possible that these GDH parameters are related to strong galactic shocks and other weak shocks in the Milky Way.


INTRODUCTION
Star formation is a sequential process composed of many stages from diffuse gas to a protostar, and each process takes place in different spatial scales from a few kpc in a galactic disk to 0.01 pc in a protostar core.Thanks to sub-arcsec submillimetre observations, many investigations into the final process have been made in recent decades.On the contrary, the first stages of the interstellar medium (ISM) and molecular cloud complex have been less studied.
In these early stages, there should be some physical connection between the ISM properties in kpc and sub-pc scales because some quantitative relations, such as the Kennicutt-Schmidt (KS) law, are found over these scales (Schmidt 1959;Kennicutt 1989).Another suggestion is a close relationship between the galaxy's star formation rate, efficiency, gas distribution, kinematics, and energetics (Handa et al. 1990;Gouliermis et al. 2010;Momose et al. 2010;Grasha et al. 2017;Fujimoto et al. 2019;Zhang et al. 2019;Villanueva et al. 2022;Maeda et al. 2023;Demachi et al. 2023).
★ E-mail:ren.matsusaka.jp@gmail.com† E-mail:handa@sci.kagoshima-u.ac.jp To address these early stages in star formation, statistics on the ISM density structure is a promising approach.The histogram of the gas density over an assigned area will give a big hint because it will provide information about the properties and evolution of the density fluctuation of the ISM, which is equal to the probability distribution function (PDF) of the gas density if the fluctuation is steady.Many observational investigations using a histogram of the column density, the so-called N-PDF, have been performed on the molecular cloud scale, or a few tens of pc (e.g.Stutz & Kainulainen 2015;Burkhart et al. 2015;Lin et al. 2017;Pouteau et al. 2023).Analogous to N-PDF studies, some investigations have used the brightness distribution index (BDI, Sawada et al. 2012aSawada et al. ,b, 2018)).The BDI is derived from the brightness distribution in a single velocity channel and calculated as the pixel fraction above an assigned threshold.However, although the BDI is the direct statistic index of observational data, it is too simplified to consider the density structure of the ISM.Therefore, it has been difficult to develop it to address the origin of the ISM density structure.Although with the N-PDF, it is necessary to convert the integrated intensity of a molecular line to the ISM density, it is a closer quantity to the simulation studies, as shown below.
Many of these studies have shown that the N-PDFs of a molecular cloud or a part of it have different shapes.For example, some N-PDFs show a log-normal (LN) component with a power-law tail (PL) in its high column density side (e.g.Schneider et al. 2013Schneider et al. , 2015aSchneider et al. ,b, 2016Schneider et al. , 2022;;Pokhrel et al. 2016;Ma et al. 2020;Wang et al. 2020;Jiao et al. 2022).The PL is often observed at the star-forming site, which is interpreted as a result of gravitational contraction (e.g.Klessen 2000;Ballesteros-Paredes et al. 2011;Chen et al. 2018;Khullar et al. 2021;Bieging & Kong 2022).However, some recent investigations show that the N-PDF comprises not one but multiple LN components and less correlation with star-forming activities (Murase et al. 2023).The other approach, (e.g.Leroy et al. 2017;Liu et al. 2021), modelled the molecular gas PDF in galaxies using large velocity gradient (LVG) radiative transfer calculations to estimate CO line ratios.They have established a framework linking global CO line ratios to the mean molecular hydrogen gas density and kinetic temperature.These studies highlight the importance of improved observational knowledge of the PDF shape.
The first observational studies on the ISM density structure in the galaxy were based on the near-IR extinction map (e.g.Lombardi et al. 2008;Kainulainen et al. 2009).However, this method can be affected by contamination from foreground and background emission in the galactic plane, leading to an artificially increased gas column density.This means that N-PDF studies derived from dust continuum observations are only effective for individual molecular clouds with no foreground or background emission.They may be less accurate in studying the ISM beyond the cloud scale.To address this regime, line emission data can separate the overlapping objects along the same line of sight.
Using the line data with distance estimation, we can assign the ISM location in a 3-dimensional voxel without contamination.The diffuse ISM under the accumulation process should be located in the outer envelope or the intercloud space.The diffuse component may be missed if we use the data after any cloud identification (e.g.clump find, CPROPS, Dendrograms Williams et al. 1994;Rosolowsky & Leroy 2006;Rosolowsky et al. 2008).To include this component for investigation on the early stage of ISM accumulation processes, we obtained a statistic of the gas density without any cloud identification.We called it the Gas Density Histogram (GDH), a fractional volume of every density diagram without morphological structure identification.This approach is a great advantage for investigating the early stage of ISM accumulation processes.The density distribution or the histogram of the gas density is the same as the PDF if all the processes are steady.In this paper, we use the GDH as a term of the PDF because the ISM density structure should be evolving and may not be steady.
One of the pioneering works on the 3-dimensional GDHs of the diffuse ISM was conducted by Berkhuĳsen & Fletcher (2008).They estimated the GDH of the atomic and ionised ISM toward pulsars and stars in the Milky Way.However, their sample was small, and their results indicated that the derived GDH was consistent with a single LN.Another pioneering work was carried out by Handa et al. (2012) using CO data along the galactic equator, i.e.,  = 0.They found that the GDH shows a single LN without any PL.
To reduce the Poisson noise for each density bin, more than 1000 pixels are required to make a GDH when the density structure is assumed to be uniform in an area.It means that an area for a GDH, i.e., the effective resolution of the GDH, is at least 30 times larger than the pixel resolution.Therefore, the GDH analysis requires data with large coverage and sharp pixel resolution.Our analysis contains sufficient data to cover more than 282 × 282 pixels by using the FUGIN dataset.Thanks to a sensitive non-biased survey with a large radio telescope, we can make GDHs for each section of the sky and investigate its variety, which may be connected to the galactic structure in the sub-kpc scale.In this paper, we present the GDH distribution over a large part of the galactic plane.
We organise this paper as follows.In Section 2, we show the specification of the data.In Section 3, we calculate the volume density using the FUGIN data.Section 4 shows our resultant GDHs in each area at an assigned velocity and the effects of beam resolution on the GDH shape.Section 5 discusses the differences in GDH shape at different places in the galactic plane using the multi-LN fitting to the resultant GDHs and the GDH parameter on the longitude-velocity diagram with our corresponding interpretations.

DATA
We used 12 CO and 13 CO (=1-0) datasets from the FUGIN project, which is an acronym of the FOREST Unbiased Galactic Plane imaging survey using the Nobeyama 45-m Telescope (Umemoto et al. 2017;Torii et al. 2019).In this study, we used the data in the region of 10 • ≤  ≤ 50 • , || ≤ 1.0 •1 .The angular and velocity resolutions of the published data are 20 ′′ and 0.65 km s −1 , respectively.To improve the sensitivity, we smoothed over 8 spectral channels, with resulting velocity resolution of ∼ 5.2 km s −1 and realised twodimensional spatial smoothing with a Gaussian function to obtain a resultant angular resolution of 1 ′ .The smoothed resolution is sufficient for the density structure of the sub-kpc scale (see Section 4.2).Fig. 1 shows the longitude-velocity diagram in 12 CO (=1-0) data for the entire area covered in this study.The yellow-shaded region is representative of the whole area presented in Section 4.

DATA ANALYSIS
In this paper, we present GDHs based on the volume density instead of the column density of the ISM.In this section, we show how to estimate the volume density and the volume of a voxel from the observed data.

H 2 Column density
The volume density is the column density divided by the length of the line-of-sight (LoS).As a first step, we estimate the column density of 13 CO,  ( 13 CO), from the 12 CO (=1-0) and 13 CO (=1-0) lines.We assume Local Thermodynamic Equilibrium (LTE) between these two transitions and the beam-filling factor () to be unity (Wilson et al. 2013), where  is the velocity in km s −1 ,  13 is the 13 CO opacity, and  ex is the excitation temperature.Then, the 13 CO optical depth ( 13 CO) at each voxel should be connected through the following equation: where  MB ( 13 CO) is the main beam temperature of 13 CO (=1-0) line, and  13 () is the equivalent brightness temperature of  in Rayleigh-Jeans approximation at the 13 CO (=1-0) frequency, and  CMB = 2.7 K is the cosmic microwave background temperature.The excitation temperature is primarily derived from the brightness temperature of 12 CO,  MB ( 12 CO) through the following equation, because of the large opacity where  12 () is the equivalent brightness temperature of  in Rayleigh-Jeans approximation at the 12 CO (=1-0) frequency.However,  ex should be warmer than 15 K, which is estimated by many investigations of the temperature of molecular gas (such as Planck Collaboration et al. 2011;Sokolov et al. 2017;Shimoikura et al. 2019;Murase et al. 2022).Therefore, we set  ex = 15 K if  ex < 15 K from Eq. (4).We have confirmed that the results shown in the next section are not affected if we use different threshold temperatures.
We used the abundance ratio of H 2 to 13 CO as [H 2 / 13 CO] = 3.8 × 10 5 (e.g.Pineda et al. 2008), although the different value gives a small shift for the resultant GDH along the density axis (see also Murase et al. 2023).

Line-of-sight depth of a voxel and its volume density
For the line data, we use the kinetic distance to the ISM in the disk of the Milky Way (van de Hulst et al. 1954).To derive the kinetic distances, we assume the flat rotation curve and use the galactocentric distance to the Sun,  0 , and the rotation velocity of the local standard of the rest (LSR),  0 as  0 = 8.5 kpc and  0 = 220.0km s −1 (Kerr & Lynden-Bell 1986).
The line of sight velocity against LSR ( LSR ) of a molecular cloud (or interstellar gas) observed in the direction (,) rotating at the The rotation model provides exactly one distance  from the galactic centre for every longitude  and  LSR .There are two areas located at a galactocentric distance , with longitude  and latitude  = 0 • from the solar system, as shown by the grey shades.Note that these cannot be separated from the  LSR .galactocentric distance  is given by We used  = 0 because the observed range of the FUGIN dataset is only || ≤ 1.The heliocentric distance  of the cloud can be given as the solution of using  given by Eq. 5.An observed  −  −  voxel with a finite velocity width, , has a finite range of the galactocentric distance between  1 and  2 in the given line of sight.At the given galactic longitude, , the geometrical depth of a voxel, , at the near or far site is given by (see Fig. 2).Using the depth  we can estimate the volume density of the gas from the column density; We should note that each voxel is the sum of two separated volumes owing to the near/far ambiguity of the kinetic distance.Therefore, our data are the averages of two separate volumes.
The results are the almost same even when we use the values of  0 = 7.92 kpc and  0 = 238.9km s −1 (e.g.VERA Collaboration et al. 2020).If we use a different rotation curve, the shape of the resultant GDH is not changed.

RESULTS
This section presents the resultant GDH of each area along the galactic plane from FUGIN data, and the following sub sections present the artificial effects to the resultant GDH.

The gas density histogram
Previous studies on the N-PDF were performed for star-forming clouds and focused on its high-density end.However, the GDH is more powerful for diffuse and less dense gas because it does not require the identification of a "cloud" with a definite boundary.Therefore, we made the GDH of all (--) voxels in each geometrically square area.Fig. 3 shows 10 GDHs derived from channel maps of 18.0 • <  < 20.0 • , || ≤ 1 • with Δ = 5.2 km s −1 at different velocities as typical GDHs from the FUGIN data.The vertical lines are the sensitive limits estimated from 1 and 3 of the blank sky channel maps, respectively.Since the column density is primarily estimated from the 13 CO intensity and the difference between the lengths of LoS is negligible in the assigned (--) area, the sensitivity limit for the volume density can be estimated from the noise level of 13 CO.In Fig. 3, we can easily find that GDHs in the density range over the sensitivity limit change from place to place.This suggests that the ISM density structure differs from place to place.
We note that the gas density in the GDH is not a thermodynamical one because the volume filling factor should be much less than unity even for each sampling voxel.That is why we use  ⊙ pc −3 instead of H 2 cm −3 as the unit of the abscissa.The typical density is 10 −2  ⊙ pc −3 , which is equal to 2.0 × 10 −1 H 2 cm −3 , less than the CO critical density by three or four orders of magnitude 2 .Our estimation indicates the density averaged over a several-pc scale voxel under the low-volume filling factor.A low filling factor is consistent with lower brightness temperatures in the 12 CO (=1-0) line than the typical temperature of molecular gas or ∼ 15 K. Fig. 3 suggests that areas with different morphology of intense emission show a different GDH shape.The GDHs of areas with diffuse less-bright gas, such as those at  = 35.5 and 113.5 km s −1 show a sharp peak just above the sensitivity limit and a steep decline at the dense side.In contrast, the GDHs of areas with compact, bright features, such as those at 51.5 and 66.7 km s −1 , have another peak or excess component on the dense side.The simplest GDH is that at  = −11.3km s −1 ; it shows a close shape of Gaussian noise on the blank sky.At  = 4.3 and 129.1 km s −1 , the GDHs have real features over the Gaussian noise at the densest tail, and the channel maps there show a few real features.
We found that each GDH shows a straight line below the sensitivity limit.This is due to the Gaussian noise on the blank sky, which is quantitatively consistent with the noise estimation from the data.

Spatial resolution effect
GDHs are basically free from the effect of image resolution of the data because they are made without identification of the spatial structure of the ISM.However, they are actually affected because compact peaks must be diluted with an observation beam.We call it the resolution effect.
The first significant difference is that the GDH made from the more smoothed image is the noisier.This is natural because of Poisson noise from the smaller number of voxels in the area.To obtain a reliable GDH, we need more than 10 3 pixels in a single GDH area, as mentioned in Section 1.A much better voxel resolution is required to obtain a reasonable resolution of the GDH, or the size of GDH area, by a factor of 30.
The second significant difference is the highest end of the GDH.In almost all cases, the ISM has a spatial structure with compact clumps or cores.When these clumps and cores are smaller than the observed beam size, the resultant GDH misses the highest end, and their volumes contribute to lower-density bins.The GDH near the highest end is important for discussing the final stages of molecular cores, where the self-gravity of the gas is the major process.However, the resolution of our data of a few pc is too large to discuss such a process.We, therefore, do not focus on the high-density end of the GDH (see Section 5.2).
Besides these differences, the three GDHs are consistent, and the resolution effect is negligible between the sensitivity limit and log (H 2 ) ≲ −0.8.The resolution effect shown above may severely affect the GDH near the high-density end, and it has large Poisson noises.
The voxel resolution also changes the sensitivity limit because the ordinate values of our GDHs primarily come from the 13 CO intensity.For the GDH shown in Fig. 7, the derived sensitivity limit from the rms noise level of the channel map is log (H 2 ) ≃ −2.6 ⊙ pc −3 , where the positive half of Gaussian noise with no real emission shows the peak in the log − log plot.In particular, at 61.5 km s −1 , there are two significant peaks in the 9 ′ and 3 ′ images.To discuss whether this is real, we should consider the random noise due to the observation system.For more details, see Section 5.1.
Fortunately, most of our GDHs show a peak or plateau in the resolution effect-free range.We, therefore, can discuss the shape of GDHs with sufficient accuracy.In this study, we discuss the GDH from the 1 ′ data and its shape above the sensitivity and Poisson noise limits, which shows the real density structure.

Velocity resolution effect
Although we convert the velocity difference to the distance, the cloud has internal turbulence.The kinematic distance is a nonlinear relation to the line-of-site velocity.Theoretical and observational studies suggest that velocity crowding affects the correlation analysis on the statistics of intensity fluctuations in Position-Position Velocity (e.g.The two vertical dotted and dashed lines correspond to 1 and 3, respectively.Along the ordinate axis, each plot has an error due to the number of pixels for each log  bin or the Poisson noise.However, it is not shown in this study as the effect is too small because the number of pixels in the GDH is sufficient.We show the line-of-sight velocity at the left-top corner of each panel.The full range of the grey scale is the same for all panels and is 0-3K.Log P (volume ratio) 7.8 km s 1 5.2 km s 1 2.6 km s 1 Figure 5. Velocity resolution effect of the GDHs.All GDHs in each panel are made from channel maps in 1 ′ angular resolution, but the velocity resolution is 2.6 km s −1 (magenta: dashed line), 5.2 km s −1 (black: solid line), and 7.8 km s −1 (orange: dash-dot line), respectively.The three panels correspond to the different line-of-sight velocities, as in Fig. 4. Lazarian & Pogosyan 2000, 2004;Yuen et al. 2021).We, therefore, evaluate the velocity resolution effect on the resultant GDH.
As shown in Fig. 5, the global shape of each GDH is very similar.Therefore, there is no velocity crowding effect on the GDHs covered in this study.This looks strange because a molecular cloud has a line width of several km s −1 (e.g.Simon et al. 2001;Colombo et al. 2015;Rosolowsky et al. 2021;Fujita et al. 2023).However, it can be understood because the GDH resolution, i.e., the cell size of the GDH is 2 • × 2 • , or 300 pc × 300 pc at 8.5 kpc.This resolution is much larger than the typical cloud size, or 10 pc.Although each cloud has velocity width and random velocity against the pure galactic rotation, statistics of many clouds in the single GDH cell should be less affected by the velocity resolution.
This also means that the GDH is less affected by the geometrical depth and velocity crowding.This is supported by the no systematic difference of GDHs near the terminal velocity (see Fig. 10).

The effect of the beam filling factors and 13 CO abundance
As shown in Section 3, we assume that the beam-filling factors of both CO lines are unity to derive our GDHs.In the case that the beam-filling factors are less than unity, the factors of  13 and  12 reduce both  MB ( 13 CO) and  MB ( 12 CO) from Eqs. 2 and 4, respectively.However, the beam-filling factors have little effect on the GDH shapes, even if they are less than unity; they only cause some shift along the horizontal axis.This can be easily understood when the line intensities are much stronger than the CMB and  () ≃ , i.e.,  ex ≫ 5.5 K.In this case,  ex is proportional to 1/ 12 and  13 averaged over the beam increases proportionally to  13 / 12 .It means that  for 13 CO, estimated with  12 =  13 = 1, is also applicable for other cases if  12 =  13 .Therefore, the resultant GDH is shifted along the horizontal axis and does not change its shape.To confirm the non-linear effects, we made additional GDHs using three models with different values of  12 and  13 (Table 1).These GDHs are consistent with the expectation shown above (Fig. 6).Model A and B are different sets of excitation temperature assumptions; A and B are the estimates  ex from the optically thick 12 CO and  ex = 15 K constant, respectively.Models A1, A2, and A3 are the different beam-filling factors.
13 CO abundance has a similar effect on the GDH because it changes the conversion factor linearly from  ( 13 CO) to  (H 2 ), which only produces a shift of the GDH along the horizontal axis in the log scale.The 13 C abundance only shifts by a factor of 2 in the inner galaxy (e.g.Milam et al. 2005;Giannetti et al. 2014;Jiao et al. 2021) without the galactic centre (FUGIN dataset).Even if the abundance is different for each GDH, their change is the shift of each GDH, but not a change in shape.

Subtraction of random noise contribution
The most striking feature of all GDHs is a straight line in the least dense region, although it is below the sensitivity limit.This comes, of course, from the random noise (see also Ossenkopf-Okada et al. 2016), as mentioned in the previous section.In a log-log plot, Gaussian noise appears as a straight line with the power index of unity below the peak, which corresponds to its rms value.In addition, the effect of the noise decays rapidly beyond the peak (see Fig. 7).We can easily estimate its contribution to the observed GDHs with a single free parameter, i.e., the volume fraction of blank sky.
As the first step, we set a GDH component of blank sky with a Gaussian noise given by where  is the rms noise level.When we plot on the log-log plane, the horizontal axis should be described by  = log   norm , i.e.,  =  norm exp(), where  norm is the normalisation factor to set  = 0.
The vertical axis should be described by  = log ().Therefore, the curve of the Gaussian is described by which is shown by a green line in each panel of Fig. 7.We can fit by the noisy blank sky because it should occupy a large portion of the low-density region.Therefore, we can remove it from the observed GDH.However, we cannot apply it for the noise with very weak emission below the noise level, because we do not know the actual brightness distribution below the noise level.Beyond the noise level, the actual emission must be dominant.Therefore, we will discuss the profile beyond the noise level (Fig. 7: higher than the dashed line, Fig. 9: non-shaded area).
Using this estimation, we subtract its contribution and determine whether there are components of GDH close to the sensitivity limit.It can give the ISM density structure in an outer region beyond the completeness limit defined by previous N-PDF works on a molecular cloud (e.g.Schneider et al. 2015b;Alves et al. 2017;Ma et al. 2020;Spilker et al. 2021).Although it is very important to investigate the density structure in the least dense envelope of ISM, it is beyond the scope of this paper.

Multi log-normal fitting
Although some GDHs can be fitted by a single LN3 , others require excess components on the single LNs (Fig. 9).Most previous studies on molecular cloud's N-PDFs used one or two PL components to explain the excess.However, we propose that the second LN component naturally fit the excess because of the following two reasons.First, the double LN model expresses the observed GDH better than the LN plus PL (see Fig. 8).We evaluated the wellness of the fitting using the adjusted coefficient of determination (adj 2 ) defined as follows: where  is the number of free parameters,  is the number of data points in the fit, y  is the volume fraction of the th data point in the fit, y ,model is the th data point by each model fit, and ȳ is the average of y.The adj 2 closer to unity means that they can perfectly explain the observed GDH.This adj 2 method is often used when the best number of parameters is unknown and it should be found by actual data (e.g.Shao et al. 2010;Mann et al. 2014).Many GDHs fitted by single or double LN models have values of adj 2 above 0.95 and are much better than a single LN plus one or two PLs.We also found that the GDH has no systematic deviation from double LN fittings.
Although a small number of GDHs can be fitted slightly better with the LN + PL model.We choose the double LN model for all GDH to simplify the interpretation, which will give a uniform model for the whole galaxy (Fig. 8).
Second, many previous works interpreted that these PL components are made by self-gravity.However, as Murase et al. (2023) pointed out, the self-gravity effect on GDH appears in 1 pc scale (Khullar et al. 2021), which is much smaller than the pixel resolution in this work.
We, therefore, fitted our GDHs with single or double LN distributions, following the approach of Gazol & Kim (2013); Brunt (2015); Schneider et al. (2022); Murase et al. (2023).Therefore, the fitting function we used is the following: where  L and  H are the mean volume density of the high-and lowdensity LN distributions, respectively.We call these two components as L-LN and H-LN for the lower and higher density LN components, respectively.
We also introduce a GDH parameter of the H-LN fraction.We calculated the volume ratios of the H-LN components to the total for each GDH defined by which indicates the fraction of the H-LN component in the total gas.Some GDHs can be well-fitted without the second LN component.In this case, we call it a GDH with no H-LN component or  H () = 0.
We use an adj 2 to perform model selection.More than 90% of all areas determined the number of LN components by this statistical method.However, the statistical method is not sufficient (e.g.Ma et al. 2022;Schneider et al. 2022).We, therefore, also check the residuals of the best fit of the two models for each area.The residuals on the log −P scale are plotted below the GDHs, as shown in the panel attached to Fig. 9.If the residuals are not distributed around zero, or if there is a large residual around the peak density that is the most highly weighted (see Fig. 9,  = 66.7 km s −1 ), we consider the fitting to be not good.If the selected model based on adj 2 unsatisfies Log P (volume ratio) Log P (volume ratio) the residual criterion, the GDH model is selected by eye.However, these situations occur only in a few percent of cases.
Some GDHs throughout the galaxy suggest the presence of third LN-like components, for example, at  = 51.1 km s −1 .However, these examples are limited, so we only consider up to the second components.

Interpretation of double LN components
Using numerical simulation, Vazquez-Semadeni (1994) showed that a GDH under a nearly pressureless supersonic flow without selfgravity is a single LN.They interpreted that the "random walk" modulation of the gas density amplification makes Gaussian distribution along the logarithmic scale of the gas density, i.e., LN.In this case, the peak and the width of one LN component are the ini- tial density of the gas and the evolutional stage of the random walk process or the modification factor at a single modification process, respectively.
We found that nearly all GDHs of 2 • × 2 • areas can be expressed by double LNs (Fig. 9).This means that the ISM in the single 2 • × 2 • area is a mixture of two different density structures.However, we cannot separate them spatially.We can decompose the set of voxels into components based on the fitting and determine the fraction of each component at an assigned gas density.However, we cannot assign any voxel to a member of a given component.A method based on morphology with the help of the GDH analysis is required to separate them spatially in a sophisticated manner.
Although we cannot separate each component spatially, our results suggest that H-LN and L-LN correspond to compact and diffuse gas, respectively (a similar model is also suggested by Gazol & Kim (2013)).We found that a single LN area shows only a few intense clumps embedded in the diffuse emission in its channel map, and a double LN area shows some intense clumps (see Fig. 3).
The GDH comprises only LN, and no PL component is natural.It is not inconsistent with previous investigations, which showed the GDH or N-PDF of molecular clouds.We note that the resolution of our GDH is 2 • or typically 300 pc, and any structure less than 20 ′′ or 0.8 pc can affect an observed GDH negligibly because this resolution is larger than the scale where gravitational collapse is the dominant process in the density evolution (e.g.Chen et al. 2018;Khullar et al. 2021), as Murase et al. (2023) pointed out.
As we shown in Section 3, our estimation of volume density is an average of two parts, i.e., near and far sides.However, this produces only small effects on the GDH because of the following two evidences.One is that we found no common feature in GDHs near the terminal velocity, where the kinematic distance is given as the double root.The prominent features of the GDH parameters shown in Fig. 10 are not limited to the  −  area along the terminal velocity but continue to the other.The other evidence is that the GDHs are less affected by resolution effects both in velocity and space (see Section 4).Therefore, the different source beam couplings do not severely affect the GDHs.

Distribution of the GDH parameters on 𝑙 − 𝑣 plane
GDHs can be expressed by some parameters because almost all observed GDHs can be well-fitted by single or double LN components, as shown in the previous section.We choose two independent parameters from five: the fraction of the H-LN component  H and the half-width of the L-LN component  L .Fig. 10 shows their distribution and the 12 CO integrated intensity on the  −  plane.We found three features, named ridges A, B, and C, which connect from (, ) ≃ (10 • , 25 km s −1 ) to (25 • , 100 km s −1 ), from (, ) ≃ (10 • , 25 km s −1 ) to (30 • , 80 km s −1 ), and from (10 • , 10 km s −1 ) to (50 • , 60 km s −1 ), respectively.These ridges correspond to the spiral arms features traced in the CO intensity distribution (Fig. 10a).Ridges A, B, and C, are the Norma, Scutum, and Sagittarius arms, respectively (e.g.Sanders et al. 1985;Dame et al. 2001).
As it is well known, Hii regions are associated with spiral arms, and they are also associated with ridges A, B, and C (e.g.Georgelin & Georgelin 1976;Downes et al. 1980;Anderson et al. 2012).We note that there is no direct relation between  H and  L .Actually  H and  L show the different distribution between 20 • and 40 • along  ≃ 10 km s −1 , which is the Local Arm.The Local Arm is a clear feature in the 12 CO intensity from (, ) ≈ (20.0 • , 10 km s −1 ) to (40 • , 10 km s −1 ) (Fig. 10a; e.g.Reid et al. 2016Reid et al. , 2019)).The Local Arm can also be traced in  L , although no feature is found in  H .This suggests that the Local Arm is a physically different entity from the other arms.Although Xu et al. (2013) suggest that the Local Arm has similar kinematic properties as those found for the prominent spiral arms, it shows less active star formation (see also Nakanishi & Sofue 2016;Swiggum et al. 2022).This suggests that the broad  L is not made by star formation activity but by a kinematical effect in the sub-kpc scale.

Double LN components and the spiral arms
As shown in Section 5.4, we found that the distribution of the CO line intensity is similar to those of the high  H and wide  L (see Fig. 10b,c).What do they suggest?
Based on the interpretation by Vazquez-Semadeni (1994), an LN component is made by a "random walk" process in the ISM.A typical "random walk" process consists of multiple shock passages produced by random supersonic turbulence.This shock must be on a smaller scale than the size of a GDH area.This is the natural interpretation of a single LN area.However, many areas show double LNs.Double LNs cannot be produced by the same mechanism as that of a single LN.Another process to make additional LNs with different peak densities is to give a different initial or mean density in the same area.What is this process?
The gravitational collapse is inadequate.Numerical simulation suggests it happens on much smaller scale (e.g.Chen et al. 2018;Khullar et al. 2021).Although star formation feedback may compress the ISM, it is also inadequate.The feedback can affect the surrounding ISM only in a 1 pc order (Rumble et al. 2021;Murase et al. 2022, and references therein).Furthermore, Egusa et al. (2018) suggest that there is no significant effect of supernovae (SNe) on the kpc scale GDH.However, more detailed studies of nearby galaxies using pc scale resolution observations and simulations are needed to discuss the effect of the SNe.
A shock on a similar or larger scale than a GDH area is promising.The passage of the galactic shock (Fujimoto 1966(Fujimoto , 1968;;Williams et al. 1994) is a single shot and can produce coherent compression, which produces a shift of an LN component.Therefore, it can produce H-LN from L-LN efficiently.In this case, high  H at the major spiral arm can be interpreted naturally (e.g.Wilson & Scoville 1991;Tosaki et al. 2007).
A wider  L is produced by more steps of density modification or by a larger change at each modification step.We cannot distinguish them only from the GDH shape.In the former case,  L shows the evolutional stage of L-LN (Ward et al. 2014).The width should be wider along the gas flow.Unfortunately, we cannot see it in galactic data due to the near-far ambiguity.GDH analysis of nearby galaxies using sensitive and high spatial resolution data will provide the answer.
In the latter case, the density modification is activated on a comparable or larger scale than a GDH area.The galactic shock may activate the small-scale shock via the cascading process.However, this mechanism must work in a short time scale because a wide  L has no systematic shift from the spiral arms traced as intensity ridges, which should be associated with the galactic shock.To clar-ify this, time evolution studies using large-scale and high-resolution ISM simulations with small CO clumps are required.
Fig. 11 shows a relation between  H and  L .There are correlations but large variance between  H and  L .This implies that the pointto-point correlations are poor.However, as can be observed in Fig. 10, both  H and  L are well associated with the spiral arms on the  −  plane.We guess this is due to a slightly different way of the association.As one important result, we found that the Local Arm is different from the major spiral arms.It shows wide  L but low  H . On the contrary, the part of 30 • ≲  ≲ 35 • of ridge C shows high  H but narrow  L .They show no direct relation between high  H and wide  L .This gap is beyond our study.To clarify it, we need a more detailed investigation and it should be difficult to address this issue only with the Milky Way data.We plan to address this issue using a nearby face-on galaxy, and it will be presented in the forthcoming papers.

CONCLUSIONS
We used all of the inner galaxy side data from FUGIN to investigate the characteristics of the gas density distribution in the Milky Way.Our conclusions are summarised as follows: (i) We present the volume density histograms (GDHs) of 2 • × 2 • areas of the inner disk of the Milky Way Galaxy using FUGIN data with the kinematic distance.
(ii) All GDHs above the sensitivity limit can be well-fitted with single or double LN distributions.
(iii) The two GDH parameters,  H and  L , show three or four coherent enhancements on the  −  plane.The pattern well traces the prominent spiral arms.However, the local spiral arm is traced by the wide  L but not by the high  H .
(iv) The GDH features traced by the spiral arms are not due to the star formation feedback, such as Hii regions, because the star formation feedback can produce an effect on its GDH only in a smaller-scale, and it has little effect on our GDH.
(v) These results provide information about what happens in the spiral arms.The interstellar shock waves, both supersonic and subsonic, will affect the structural evolution of the ISM.The large-scale strong shock will produce a peak density shift of an LN component in the GDH, and the small-scale weak shock or subsonic pressure modification will result in a wider density dispersion of the lowdensity LN component.This means that the GDH is a powerful tool to investigate the density evolution of the ISM, although we need finer simulations.This work suggests that the GDH analysis is useful to address the ISM response to large-scale dynamics induced by galactic potential.We can confirm the relation between the GDH parameters and spiral arms based on this GDH analysis using high-fidelity CO images of nearby galaxies observed with the Atacama Large Millimeter/submillimeter Array (e.g.Leroy et al. 2021;Koda et al. 2023;Muraoka et al. 2023).This should open the door to a comprehensive picture of the evolution of the ISM in a galaxy.
Longitude-velocity diagram in 12 CO (=1-0) data obtained from smoothed FUGIN datasets.The data were integrated with the galactic latitude over || ≤ 1.0 • .The units on the intensity scale are K degrees.The vertical feature near  = 31 • is an artefact from the on-the-fly mapping.The yellow shaded region shows the  = 18.0 • -20.0 • region.

Figure 3 .
Figure 3. GDHs and intensity distribution on the sky of the corresponding areas, centred at (, ) = (19.0•, 0 • ) at different velocities.The total number of pixels is 282 × 282 (=79542).The two vertical dotted and dashed lines correspond to 1 and 3, respectively.Along the ordinate axis, each plot has an error due to the number of pixels for each log  bin or the Poisson noise.However, it is not shown in this study as the effect is too small because the number of pixels in the GDH is sufficient.We show the line-of-sight velocity at the left-top corner of each panel.The full range of the grey scale is the same for all panels and is 0-3K.

Figure 4 .
Figure 4. Resolution effect of the GDHs.The three GDHs in each panel are made from channel maps at 1 ′ resolution (black: solid line), 3 ′ resolution (red: dashed line), and 9 ′ resolution (blue: dash-dot line).These three panels correspond to different line-of-sight velocities of 19.9 km s −1 , 35.5 km s −1 and 61.5 km s −1 , respectively.The vertical dot-dashed line corresponds to 1.

Figure 6 .
Figure6.The four GDHs in each panel are derived using different assumptions of  ex and beam filling factors.The angular resolution of all GDHs is the same as that of 1 ′ .Model A1 is the solid (black, this work) line, Model A2 is the blue (dotted) line, Model A3 is the dashed (magenta) line, and Model B is the dash-dotted (green) line.

Figure 7 .Figure 8 .
Figure7.Same as Fig.4, but after subtraction of the random noise contribution shown in thick green lines.The GDH before subtraction is shown in a solid grey line.The noise model is based on blank sky with the random noise estimated from the observations; the fitting parameter is only along the ordinate axis.The filled circle mark is subtracted from the random noise effects.The 13 CO rms level is 0.07 K in (, )=(19.0• , 0 • ) areas.

Figure 9 .
Figure 9. GDHs after removing the Gaussian noise component and multi LN fitting at  = 18.0 • -20.0 • ,  = ±1 • .The vertical dashed lines correspond to 1.The L-LN and H-LN components are marked by blue and red dashed lines, respectively.Light grey dots show GDHs outside the discussion because of the low reliability of the noise removal procedure or poor pixels per bin.All GDHs were fitted using the least-squares method over 1 and sufficient pixel density per bin (N>10 1.5 ) for the random noise subtracted data.The log −P scale residuals of the single or double LN fitting are shown in an attached panel just beneath each GDH.In the residual plot, black circles show double LN residuals, and open circles show single LN residuals.In the lower right of the attached panel, we show adj 2 for each of the models used.

Figure 10 .
Figure 10.(a)  −  diagram of the FUGIN 12 CO =1-0 data (same as Fig. 1).The black line shows the curve of the terminal velocities.The coloured lines show the loci of the galactic arms constructed (Reid et al. 2016; Torii et al. 2019) and the Hii regions (red points: Anderson et al. 2009).(b)  −  distribution of the  H averaged over Δ × Δ = 2 • × 2 • .The 2 • × 2 • area where the H-LN cannot be defined (i.e., only diffuse gas, see Section 3) is shown as  H =0. (c) The same as (b), but the colour is  L .

Figure 11 .
Figure 11.Correlation between  H and  L .The correlation coefficient is 0.80.The mean value and standard deviation of  L in each 0.1 steps of  H are over-plotted as white-open circles and bars, respectively.It is noted that the case  H () = 0 is not plotted.