999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

Can Turbulent, High-density Gas Form Stars in Molecular Clouds: A Case Study in Ophiuchus

2022-08-01 01:47:36SihanJiaoJingwenWuHaoRuanYuxinLinChaoWeiTsaiandLinjingFeng
Research in Astronomy and Astrophysics 2022年7期

Sihan Jiao, Jingwen Wu, Hao Ruan, Yuxin Lin, Chao-Wei Tsai, and Linjing Feng

1 National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China; sihanjiao@nao.cas.cn, jingwen@nao.cas.cn

2 University of Chinese Academy of Sciences, Beijing 100049, China

3 Centre for Astrochemical Studies, Max-Planck-Institut für Extraterrestrische Physik, Gie?enbachstra?e 1, D-85748 Garching, Germany

Received 2022 March 4; revised 2022 April 7; accepted 2022 April 12; published 2022 June 22

Abstract Star formation is governed by the interplay between gravity and turbulence in most of molecular clouds. Recent theoretical works assume that dense gas, whose column density is above a critical value in the column density probability distribution function(N-PDF),where gravity starts to overcome turbulence,becomes star-forming gas and will collapse to form stars. However, these high-density gases will include some very turbulent areas in the clouds. Will these dense but turbulent gases also form stars? We test this scenario in Ophiuchus molecular cloud using N-PDF analysis and find that at least in some regions, the turbulent, dense gas is not forming stars. We identified two isolated high-density structures in Ophiuchus,which are gravitationally unbound and show no sign of star formation. Their high densities may come from turbulence.

Key words: ISM: clouds – stars: formation – ISM: structure

1. Introduction

The interactions between gravity, turbulence, and magnetic fields determine the formation of the molecular cloud cores and how they collapse to form stars (Shu et al. 1991; McKee &Ostriker 2007).A fundamental problem of star formation is the very low star-forming efficiency,that the theoretical estimate of the star formation rate in the Milky Way is about 100 times greater than the rate inferred from observations, assuming the molecular clouds are gravitationally bound and will collapse in freefall(Zuckerman&Evans 1974;Zuckerman&Palmer 1974;Krumholz et al. 2014; Padoan et al. 2014). For example, the freefall time of a molecular cloud is

where ρ is the density of the cloud,G is the constant of gravity.Adopting a characteristic density of 100 cm?3for a molecular cloud,the freefall time is 3.34×106yr,corresponding to a star formation rate of 300 M☉yr?1. While based on the Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE)(Benjamin et al. 2003; Churchwell et al. 2009) survey,Chomiuk & Povich (2011) derived a star formation rate of 1.9±0.4 M☉yr?1by counting young stars,which is two orders of magnitude smaller than the freefall prediction. Turbulence and magnetic fields must have played important roles in counteracting gravity to cause this low star-forming efficiency in molecular clouds (Larson 1981; Brunt et al. 2009;Crutcher 2012; Lazarian et al. 2012).

Recently, Evans et al. (2021) pointed out that this disagreement could be partly relieved: only part of molecular clouds are gravitationally bound in the Milky Way. If only those gravitationally bound clouds are counted in a freefall collapse model,the predicted star formation rate over all clouds in a complete survey of the Galaxy is 46 M☉yr?1, a factor of 6.5 less than the value assuming all clouds are bound.But there is still more than one order of magnitude difference between the freefall model and real observations, leaving to turbulence and magnetic fields to account for.

The column density probability distribution function (NPDF) is a commonly used tool to quantify the distribution of gas in star-forming regions and to identify the gravitationally bound structure in the clouds. Observations and simulations have revealed that the high column density regime of the N-PDF normally traces self-gravitating dense molecular gas that will collapse to form stars,and often has a power-law form distribution (Kainulainen & Tan 2013; Lombardi et al. 2015;Schneider et al. 2015), while the lower column density part in the N-PDF traces turbulent, low-density gases, and is welldescribed by a lognormal form (Vazquez-Semadeni 1994;Kainulainen & Tan 2013). This gives a tentative approach to separate the turbulent-dominant gas and the gravity-dominant gas. Burkhart & Mocz (2019) identified the column density where the N-PDF turns from lognormal form to power-law form as a critical column density, above which gravity will dominate over turbulence, and argued that all higher column density gases in the N-PDF are star-forming gas, meaning that they will collapse under gravity and form stars.

Such a scenario defines the high-density gas in molecular clouds as star-forming gas. However, given the lognormal distribution of the turbulent gas,whose tail extents to the highdensity regime in N-PDF,are those high-density, but turbulent gases also star-forming gas? Will they eventually collapse to form stars?

Given the N-PDF have the advantage to recover the distribution of gas within some column density ranges in a column density map (for example, a dust emission map is a good approximation of a column density map for most of molecular clouds), we can test this scenario in some nearby molecular clouds to see whether such high-density yet turbulent regions can form stars.

In this paper, we use the N-PDF analysis to the Ophiuchus molecular cloud,a nearby star-forming region,to investigate their star formation properties at high-density regions. Ophiuchus is located at a distance of ~140 pc (Ortiz-León et al. 2018). It is made up of two dark clouds, L1688 and L1689 (Lynds 1962).L1688 contains more than 20 embedded protostars and two B stars, while L1689 contains five embedded protostellar systems(Enoch et al.2009;Dunham et al.2015).Because of its proximity and active star formation, the Ophiuchus molecular cloud is an ideal target to apply N-PDF analysis and to investigate the properties of high-density structures.

The rest of the paper is organized as follows: In Section 2,we summarize the continuum and spectral line observations of the cloud and the data reduction process. In Section 3, we calculate column density and generate the N-PDF based on dust continuum. The discussion of whether individual highdensity regions can form stars is presented in Section 4. Our main conclusions are given in Section 5.

2. Observations

2.1. Herschel Dust Continuum Emission

Dust emission is optically thin in most of the molecular clouds. We use Herschel4Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.data to generate the column density map of the Ophiuchus molecular cloud. We retrieved the Herschel Gould Belt survey (HGBS) data that were taken at 70/160 μm using the PACS instrument (Poglitsch et al. 2010)and at 250/350/500 μm using the SPIRE instrument (Griffin et al. 2010). The HGBS took a census in the nearby (0.5 kpc)molecular cloud complexes for an extensive imaging survey of the densest portions of the Gould Belt, down to a 5σ column sensitivityNH2~1021cm?2or AV~1 (André et al. 2010). We retrieved the level 2.5 processed, archival Herschel images(obsID: 1342227149, 1342205094). Since we are interested in extended structures, we adopt the extended emission products,which have been absolute zero-point corrected based on the images taken by the Planck space telescope.

2.2. The 12CO, 13CO and N2H+ Data

In order to investigate the dynamical information and virial parameters of the cloud, we also obtained some molecular line data for the Ophiuchus molecular cloud.12CO(1–0)and13CO(1–0) data were obtained from the COordinated Molecular Probe Line Extinction Thermal Emission Survey of Starforming Regions survey (COMPLETE survey) (Ridge et al.2006), which is carried out during 2002–2005 with the Five College Radio Astronomy Observatory (FCRAO) telescope at a resolution of 47”. Data were taken during a wide range of weather qualities, and system temperatures were generally between 500 and 1000 K at 115 GHz, between 200 and 600 K at 110 GHz. This OTF data cover 10 square-degree region toward the Ophiuchus molecular cloud with an effective velocity resolution of 0.07 km s?1,and a mean rms per channel of 1.0 K for12CO data and 0.3 K for13CO data.

The 2.5 square-degree N2H+map of Ophiuchus molecular cloud was observed with Delingha 13.7 m telescope. The system temperatures ranged from 160 to 200 K and the variation of the system temperature caused a±0.1 K km s?1variation of the rms values of the integrated intensity.The N2H+data have a velocity resolution of 0.2 km s?1and an integrated intensity rms of 0.8 K km s?1.The details of this N2H+observation can be found in Pan et al. ( 2017).

3. Results

3.1. Deriving Column Density and Dust Temperature Based on SED Fitting

We performed single-component, modified blackbody spectral energy distribution (SED) fits to each pixel of input Herschel PACS/SPIRE images. Before performing any SED fitting,we smoothed all images to a common angular resolution of the largest telescope beam and all images were re-gridded to have the same pixel size. We weighted the data points by the measured noise level in the least-squares fits.we adopt the dust opacity per unit mass at 230 GHz of 0.09 cm2g?2(Ossenkopf& Henning 1994), and we assume a gas-to-dust mass ratio of 100. For the modified blackbody assumption, the flux density Sνat a certain observing frequency ν is given by

where the gas+dust column density N can be approximated by

Figure 1. Dust column density distributions of the Ophiuchus molecular and the N-PDf. (Left) Dust column density map based on Herschel observations. The red contour indicates the transitional column density when a log-normal distribution turns into a power-law form,corresponding to a column density 3×1021 cm?2.The orange contours mark the column density of 5×1021 cm?2, where the log-normal contribution drops to one percent of overall fitting at the density. The young embedded protostars (Class 0, I and Flatspectrum sources) are presented with magenta crosses (Dunham et al. 2015) while the star 22 Scorpius and the star HD 147 889 are shown with lime stars.The two dense and turbulent regions with only red contours but no orange contours(Region 1 and Region 2)are indicated with a cyan box and a pink box,respectively.(Right)Corresponding N-PDFs of Ophiuchus molecular cloud.The fitted lognormal and power-law distributions are indicated as a solid red line and a dashed green line. The red and orange vertical dotted lines indicate the column densities of the red and orange contours in column density maps, respectively.

3.2. The Column Density Distribution and Cloud Structures

For the Ophiuchus molecular clouds, we obtained the column density maps with an angular resolution of 36″,matched to the longest wavelength Herschel band. In Figure 1(left),we show the derived column density map of Ophiuchus.

To describe the N-PDF, we used the notation η, following the frequently used formalism from previous numerical works.The normalization of the probability function is given by

4. Discussion

4.1.Using N-PDF to Identify Substructures of the Clouds

The N-PDF method has been proved to be an effective way to separate self-gravity-dominated component and turbulent-dominated component. Burkhart & Mocz (2019) argued that for gas denser than the transitional column density in N-PDF,the gravity will overwhelm turbulence that the gas becomes star-forming gas,and would collapse to form stars. In Ophiuchus, this means that regions with column density greater than 3×1021cm?2should be gravitationally bound and will collapse to form stars. We can recover the pixels in the column density map whose column density is above this critical value and check the distribution of these high-density gases.

We plot the contours of column densityNH2=3×1021cm?2in red in Figure 1 (left). Within this contours should be starforming gas if Burkhart&Mocz(2019)’s claim is correct.In the N-PDF we also notice the overlapping region of lognormal and power-law distributions, that a large fraction of the gas within these density ranges is both dense and turbulent. To define and trace these dense, turbulent gas, we define a transitional region in N-PDF, staring from the transitional column densityNH2~3×1021cm?2, and ending at column density where the contribution from turbulence (lognormal form) drops to one percent to the overall contribution(power-law distribution),which corresponding to ~5×1021cm?2in Ophiuchus.We indicate the starting and ending column densities of these transitional regions in the N-PDF plot as vertical red and orange dotted lines,and plot another contour ofNH2=5×1021cm?2in orange in the column density map in Figure 1.Thus,we can see the distribution of these dense but turbulent gas in the column density map,by looking at the area between the red and orange contours.

As we can see from Figure 1 (left), for most of the clouds,the orange contours are surrounded by the red contours,suggesting the gas within orange contours are very dense and really dominated by gravity,and turbulence may still influence some areas surrounding these very dense regions. The overall power-law shape of the high-density gas in N-PDF is a summation of the N-PDFs of all substructures inside the cloud complex, likely reflecting their overall contribution to the gravitational potential that triggers star formation. Detailed analysis of the anatomy of the N-PDF in Ophiuchus can be found in Chen et al. (2018).

However, we also find two gas structures that have only red contours but no orange contours inside, suggesting they are composed of only transitional gas that are denser than the transitional column density,but very turbulent that do not form even denser clumps inside. We name them Region 1 and Region 2 as indicated in the cyan and pink boxes in the column density map (Figure 1 (left)). These two gas structures are isolated from the main cloud of L1688.They will be classified as star-forming gas following the N-PDF classification method,but can they collapse to form stars?

4.2. Will Gas Structures in Region 1 and 2 Collapse to form Stars?

We generate the N-PDF plot for these two gas structures in Region 1 and Region 2, as shown in Figure 2. They are both located entirely in the transitional region between 3×1021cm?2and 5×1021cm?2, but they show a simple lognormal distribution that indicatives of turbulence contribution, no any higher-density power-law contribution presented in their N-PDFs.

We can further investigate whether these isolated substructures are gravitationally bound or unbound,to judge if they can collapse to form stars.The virial theorem(Mestel&Spitzer 1956)can be simply expressed as

without considering the surface kinetic term and magnetic fields. Here I is a quantity proportional to the inertia of the cloud, andI¨ determines whether the cloud is expanding or collapsing. For a cloud with a fixed mass M and velocity dispersion σν, the term

is the total kinetic energy in the cloud including both thermal and nonthermal (turbulence, rotation, infall, etc.) components.The gravitational termW of an ellipsoid cloud radius R is

Figure 2.The N-PDf of Region 1 and Region 2.(Left)Zoom-in on the right panel of Figure 1.The magenta histogram shows the N-PDF of Region 1.The transitional region (between 3×1021 cm?2 and 5×1021 cm?2) is shown as the slateblue area. (Right) Similar to the left panel but for Region 2.

the most likely value for a prolate ellipse, where1F is the Appell hypergeometric function of the first kind.

For Region 1 and Region 2,we need an optically thin tracer to illustrate the velocity dispersion of the cloud. Besides the12CO and13CO maps in the COMPLETE survey, the Ophiuchus molecular cloud has also been mapped with other denser gas tracers (e.g., C18O Yun et al. 2021, HCN Shimajiri et al. 2017, HCO+Shimajiri et al. 2017, and N2H+Pan et al.2017). However, Shimajiri et al. (2017) only focused on the center region of L1688 that did not cover Region 1 and Region 2. We see no C18O and N2H+detections toward Region 1 and Region 2 due to the very low column density (Pan et al. 2017;Yun et al.2021),which also suggests13CO could be treated as optically thin tracers in these two regions.Therefore we use the average spectrum of the13CO emission (see Figure 3) of Region 1 and Region 2 to obtain the velocity dispersion(σν)by fitting a Gaussian function. The effective radius R is yielded following

We can estimate the mass of Region 1 and Region 2 based on12CO and13CO data by assuming these two regions are under the conditions of the local thermodynamic equilibrium(LTE) and the molecules are uniformly excited (Qian et al.2012; Ma et al. 2020). Assuming that the12CO lines are optically thick,we calculate the excitation temperature with the equation:

where Tpeak(12CO) is the peak main-beam brightness temperature of12CO line (Li et al. 2018). Assuming that the13CO lines are optically thin and the excitation temperatures of the12CO and13CO are the same, we can estimate the column density of the13CO molecule with the equation:

whereN13COis the column density of13CO, Tmbis the main beam brightness temperature (Li et al. 2018). Assuming an abundance of[12CO]/[13CO]=69(Wilson&Rood 1994)and[H]/[CO] ~104(Solomon & Klemperer 1972), we obtain the mass of Region 1 and Region 2. Based on the mass estimated from CO data, the derived ratios ∣2T W∣of Region 1 and Region 2 are 37.7 and 50.7 respectively,which are greater than one by more than an order of magnitude as well. All the derived ∣2T W∣ratios of Region 1 and Region 2 indicate these two regions are unbound.

Figure 3. Observed 13CO (black) and N2H+ (blue) lines toward Region 1 and Region 2.

4.3. Star Formation Activities in Substructures

The active star-forming regions are often well correlated with young embedded protostars. The catalog of protostar in the Ophiuchus molecular cloud is quoted from Dunham et al.(2015),which generates a full catalog of Young Stellar Objects(YSOs)by using Spitzer in nearby 18 molecular clouds.In the left panel of Figure 1, we present the young embedded protostars (Class 0, I and the Flat-spectrum sources) distribution in the Ophiuchus cloud. It is clear that most of the young embedded protostars are located in the orange contours,confirming they are experiencing on-going star formation, that these high-density,purely gravity-dominated gases are directly related to star-forming activities. On the contrary, we see no young embedded protostar within the red contours in Region 1 and Region 2, showing these areas are not active in star formation.

We also show the distribution of all young stellar objects(YSOs), including protostars, Class II and Class III sources in Figure 4. The Class II and Class III YSOs are at a later evolutionary stage than the protostars(Dunham et al.2015)and show a more dispersed spatial distribution than the protostars.Assuming a Class II duration of 2 Myr(Dunham et al.2014),it is possible that the Class II and Class III YSOs have left their birth places and streamed into the diffuse region. In Figure 4 some Class II and Class III sources have shown up in the transitional regions between the red and orange contours, as well as some move outside the clouds.Only one Class II YSO is found within the red contour in Region 1, but close to the edge. No YSOs are found in Region 2. The YSO distribution agrees with our judgement that the gas structures in Region 1 and Region 2 are not undergoing star-forming activities.

4.4. Uncertainties in Using N-PDF to Quantify Cloud Properties

Now we have a full consideration of all the possible uncertainties while we use the N-PDF method to quantify the properties of Region 1 and Region 2 to see if any of these uncertainties may significantly affect our major conclusions:Region 1 and Region 2 are dense but turbulent regions; using N-PDF with N greater than a transitional column density to define star-forming gas over a large area could be problematic;some dense but turbulent regions may not be forming stars.

These uncertainties may come from the measurement in creating the N-PDF or the N-PDF itself. For the measurement uncertainties, the first part comes from the noise level of the measurements. As seen from Section 2.1, 1σ sensitivity for column density isNH2~2×1020cm?2, which is lower than the transitional density of ~3×1021cm?2by a factor of 15. The other uncertainty comes from the limited field of view covered by the N-PDF statistics.For example,Chen et al.(2018)pointed out that the limited field of view of the target region (L1689 in Ophiuchus) may affect N-PDF in the column density range of 0.5–3×1021cm?2. Our work has analyzed a much bigger area,including L1688 and L1689, thus sampling the N-PDF more completely, especially toward the low-density regime. Following the same method used in Chen et al. (2018), we found that the limited field of view will lead to uncertainties of the N-PDF at N<2.0×1021cm?2. Hence the uncertainty introduced by the limited field of view will not affect our discussed transitional column density range of 3–5×1021cm?2in Region 1 and Region 2. Therefore, the measurement uncertainties will not change our conclusions that Region 1 and Region 2 are dense and turbulent regions.

Figure 4.Similar to the left panel of Figure 1 but for all YSOs(protostars,Class II and III sources).While the young embedded protostars(Class 0,I and Flatspectrum sources) are presented with magenta crosses, Class II and III sources are presented with cyan crosses (Dunham et al. 2015).

There are also uncertainties introduced by using the N-PDF methodology. One of the main caveats of utilizing N-PDF is the cloud geometry and the projection effect to convert volume density to column density (Stutz & Gould 2016). The main issue would be to avoid there are two or more cloud structures along the line of sight that could affect the integrated column density. For the possible source confusion problem, we can investigate molecular line mapping observations to resolve them by checking any additional components in velocity space.We checked13CO spectra,13CO moment one map,and column density map of Region 1 and Region 2,confirming there are no apparent overlapping components along the line of sight in Region 1 and Region 2.

Different turbulence driving mechanisms may also introduce uncertainties in the N-PDF method. Based on numerical models and observations, Burkhart & Lazarian (2012) shows that the width of the lognormal (ση) can be expressed as a function of the sonic Mach number (sM) and the forcing parameter,

Figure 5. Second-order velocity structure function of Region 1 and Region 2. The grayshaded region is the spatial-resolution limit which corresponds to 47″. The power-law fitting result is indicated by the black dashed line, while the lime lines represent the the classic Larson relation (γ=0.76; Larson 1981).

where b is the forcing parameter, varying from ≈1/3 (purely solenoidal forcing) to 1 (purely compressive forcing), and A is the scaling constant from volume density to column density(Federrath et al. 2008; Burkhart et al. 2017; Pan et al. 2019).Assuming that the gas and dust are in thermal equilibrium, we can use the dust temperature from Herschel data and the velocity dispersion(σv)from13CO data to derivesM by using where csis the sonic sound speed. The derived value ofsM is≈7 for the entire Ophiuchus cloud and ≈9 for Region 1 and Region 2. If Region 1 and Region 2 have similar forcing parameters and scaling factors,we should expect comparableσηfor Region 1, Region 2, and the entire Ophiuchus cloud.However,as seen in Figure 2,the measuredσηfor Region 1 and Region 2 is only about half of the measuredσηfor the entire Ophiuchus cloud. It may suggest the parameter b is quite different for Region 1 and 2, and the Ophiuchus cloud,implying Region 1 and Region 2 may have different origins of turbulence from the rest of the cloud.However,when we cover the whole area to make an overall N-PDF, Region 1 and Region 2 will also be counted as dense gas and may be judged as star-forming gas,following Burkhart&Mocz(2019),which remains to be problematic. So this does not change our major argument.

Turbulence in molecular clouds is usually non-isotropic and will appear differently in different directions (Carroll et al.2010; Vázquez-Semadeni et al. 2019). This will lead to nonisotropic clustering of the molecular cloud, which also affects the accuracy of the N-PDF measurements, and cause a nonisotropic velocity field.Following previous studies(e.g.,Heyer et al. 2008), we used the second-order velocity structure function, S2, to quantify the velocity anisotropy. The secondorder velocity structure function is a two-point correlation function that quantifies the mean velocity difference:

where l is the spatial lag between two positions,x and x+l,and γ is the power-law index. We fitted ellipses to Region 1 and Region 2 by using the Dendrogram program (Rosolowsky et al. 2008) and then calculated the second-order velocity structure function along the major axes(S2major)and minor axes(S2minor). The amplitude and index forS2majorandS2minorare similar in Region 1 and Region 2(see Figure 5),indicating the turbulence is close isotropic. The resulting γ is ≈0.36 for Region 1 and 0.40 for Region 2, which is similar to previous reported values of molecular clouds (Heyer & Brunt 2007;Roman-Duval et al. 2011; Hacar et al. 2016). Thus the nonisotropic turbulence does not have significant influence on the line width we obtained.

5. Conclusions

Recent theoretical works (e.g., Burkhart & Mocz 2019)claim that dense gas whose column density is greater than a transitional column density in N-PDF becomes star-forming gas and will collapse to form stars.We test this scenario in the nearby star-forming cloud Ophiuchus, using N-PDF analysis based on Herschel observations. The high-density, turbulent gases are distributed around dense clumps within molecular clouds, or form small isolated structures outside big clouds.

For dense, turbulent gas within big clouds, they probably will still contribute to the gravitation potential of the dense clumps they are surrounding, they also host some YSOs, but their contribution to the overall star formation rate is complex.We identified two isolated gas structures (i.e., Region 1 and Region 2)with column density ranging between 3×1021cm?2and 5×1021cm?2, their gas are high-density, but turbulent.The N-PDFs of these two regions show only lognormal distribution that is indicative of turbulence contribution,and we cannot exclude the possibility that they have different turbulence origin from the main cloud. A virial analysis suggests that these two regions are gravitational unbound, and we see no sign of active star formation in these two regions.We conclude that these two gas structures,although having column densities greater than the transitional column density,are likely not star-forming gas and are not forming star.

Acknowledgments

We thank Dr. Hauyu Baobab Liu for helpful discussions on the method design and science explanation. This work is supported by the National Natural Science Foundation of China Grant Nos. 11988101 and 12041302, and the National Key R&D Program of China No. 2017YFA0402600. This research has made use of data from the Herschel Gould Belt survey(HGBS) project (http://gouldbelt-herschel.cea.fr). S.J. thanks Ningyu Tang for helpful discussions on the12CO,13CO data,and Yan Gong for helpful discussions on the second-order velocity structure function analysis. The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAF-IFSI Rome and INAF-Arcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC).


登錄APP查看全文

主站蜘蛛池模板: 亚洲无码高清视频在线观看| 男人天堂伊人网| 中日韩欧亚无码视频| 成人福利在线视频免费观看| 高潮爽到爆的喷水女主播视频 | 极品国产在线| 尤物视频一区| 亚洲人成影视在线观看| 久久国产精品波多野结衣| 好紧好深好大乳无码中文字幕| 亚洲欧洲AV一区二区三区| 中文毛片无遮挡播放免费| 爽爽影院十八禁在线观看| 91网站国产| 国产幂在线无码精品| 黄色在线不卡| 亚洲成A人V欧美综合天堂| 尤物精品国产福利网站| 成人福利在线免费观看| 亚洲黄色网站视频| 中国特黄美女一级视频| 精品自窥自偷在线看| 成人精品区| 国产交换配偶在线视频| 91麻豆久久久| 丁香婷婷激情综合激情| 色偷偷综合网| 亚洲精品爱草草视频在线| 欧美国产日产一区二区| 超级碰免费视频91| 国产日韩av在线播放| 亚洲精品爱草草视频在线| 一级毛片免费观看久| 亚洲一级无毛片无码在线免费视频| 黄色网在线| 国产美女在线观看| 人妻精品全国免费视频| 国产尹人香蕉综合在线电影 | av一区二区无码在线| 色婷婷在线播放| 中文字幕资源站| 久久久精品国产SM调教网站| 欧美在线观看不卡| 久久婷婷五月综合色一区二区| 最新亚洲人成无码网站欣赏网 | 91精品小视频| 亚洲AⅤ无码国产精品| 中文字幕在线看| 伊人久久精品无码麻豆精品| 国产黄色免费看| 色综合热无码热国产| 亚洲国语自产一区第二页| 天堂网亚洲综合在线| 亚洲码在线中文在线观看| 中文字幕久久精品波多野结| 亚洲精品在线观看91| 99热这里都是国产精品| 亚洲精品日产AⅤ| 亚洲a免费| 露脸真实国语乱在线观看| 99精品伊人久久久大香线蕉| 国产在线自乱拍播放| 伊人久久久久久久久久| 国产在线小视频| 亚洲日本一本dvd高清| 91福利在线观看视频| 国产激情国语对白普通话| 国内a级毛片| 色噜噜在线观看| 久久国产精品麻豆系列| 亚洲高清在线天堂精品| 亚洲最大看欧美片网站地址| 19国产精品麻豆免费观看| 免费国产高清精品一区在线| 国产对白刺激真实精品91| 五月婷婷综合网| 韩日无码在线不卡| 最新国产你懂的在线网址| 国产男女免费视频| 亚洲欧美日韩视频一区| 性色一区| 亚洲精品视频在线观看视频|