• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Physical Properties of 29 sdB+dM Eclipsing Binaries in Zwicky Transient Facility

    2022-05-24 14:21:26MinDaiXiaoDianChenKunWangYangPingLuoShuWangandLiCaiDeng

    Min Dai ,Xiao-Dian Chen,2,3 ,Kun Wang ,Yang-Ping Luo ,Shu Wang,2 ,and Li-Cai Deng,2,3

    1School of Physics and Astronomy,China West Normal University,Nanchong 637002,China;chenxiaodian@nao.cas.cn

    2 CAS Key Laboratory of Optical Astronomy,National Astronomical Observatories,Chinese Academy of Sciences,Beijing 100101,China

    3 School of Astronomy and Space Science,University of the Chinese Academy of Sciences,Beijing 101408,China

    Abstract The development of large-scale time-domain surveys provides an opportunity to study the physical properties as well as the evolutionary scenario of B-type subdwarfs(sdBs)and M-type dwarfs(dMs).Here,we obtained 33 sdB+dM eclipsing binaries based on the Zwicky Transient Facility (ZTF) light curves and Gaia Early Data Release 3(EDR3)parallaxes.By using the PHOEBE code for light curve analysis,we obtain probability distributions for parameters of 29 sdB+dMs.R1,R2 and i are well determined,and the average uncertainty of mass ratio q is 0.08.Our parameters are in good agreement with previous works if a typical mass of sdB is assumed.Based on parameters of 29 sdB+dMs,we find that both the mass ratio q and the companion’s radius R2 decrease with the shortening of the orbital period.For the three sdB+dMs with orbital periods less than 0.075 days,their companions are all brown dwarfs.The masses and radii of the companions satisfy the mass–radius relation for low-mass stars and brown dwarfs.Companions with radii between 0.12 R⊙and 0.15 R⊙seem to be missing in the observations.As more short-period sdB+dM eclipsing binaries are discovered and classified in the future with ZTF and Gaia,we will have more information to constrain the evolutionary ending of sdB+dMs.

    Key words:(stars:) binaries:eclipsing–stars:evolution–(stars:) subdwarfs–(stars:) brown dwarfs

    1.Introduction

    Hot subdwarf stars are located between the main sequence and the white dwarf (WD) sequence in the Hertzsprung-Russell diagram (Heber 2016).They are classified as B-type (sdB,Teff<40,000 K)and O-type(sdO,Teff>40,000 K)according to their different spectral types (Hirsch et al.2008).sdB stars are located in the blue tail of the horizontal branch(HB),which is also known as the extreme horizontal branch(EHB;Heber 1986).sdB has a helium-burning core and an extremely thin residual hydrogen envelope (mass <0.01 M⊙).Its progenitor is likely an evolved star that loses most of its hydrogen envelope when it ascends the red giant branch.sdB will not enter the asymptotic giant branch phase,instead,it will become a WD directly.The average lifetime on the EHB is on the order of 108yr (Dorman et al.1993).Maxted et al.(2001)suggested the binary ratio of sdB is~60%.About half of these binary systems are close binaries with periods of a few days (Allard et al.1994;Maxted et al.2001),while the others are wider binaries with periods of several years (Stark &Wade 2003;Chen et al.2013;Vos et al.2018).

    The evolution of sdB is still controversial (Geier &Heber 2012),although the binary evolution scenario is favored.Han et al.(2002,2003)proposed three binary evolution channels for sdBs (common envelope evolution,stable Roche-lobe overflow and the merger of two He WDs).Subsequently,a number of works tried to test these theories through observations(Luo et al.2019,2020;Kupfer et al.2020;Kramer et al.2020;Pelisoli et al.2020).For single sdB stars,the formation has always puzzled us.About 40%of sdB stars have been found to be single.Although it has been argued that all sdB stars evolved from binary stars(Pelisoli et al.2020),several facts suggest that the merger channel of two He WDs cannot satisfactorily explain single sdB stars(Burdge et al.2020;Van Grootel et al.2021).First,the number of low-mass WD binaries is small (Ratzloff et al.2019b).Second,the predicted broad mass distribution of single sdB stars from the merger channel seems to be inconsistent with observations(Fontaine et al.2012).Third,observations show that most single sdB stars have a very slow rotation,which counters the fast rotation predicted by the merger channel (Geier &Heber 2012;Charpinet et al.2018).Recently,a new channel for the formation of single sdB stars has been proposed (Meng et al.2020,2021).

    When an sdB+dM binary system is eclipsed,it exhibits an HW Vir-type light curve (LC).The number of such stars is relatively small,and since they possess a characteristic LC,they provide a direct way to measure the physical properties of components by modeling the LC and radial velocity curve,which helps to test the theory of formation and evolution of such systems.For a more detailed review on HW Vir-type stars,see Heber (2016).Such stars have been studied continuously and extensively over the last few decades,and their number is increasing (see e.g.,Kilkenny et al.1978;Menzies &Marang 1986;Drechsel et al.2001;?stensen et al.2007;For et al 2010;Barlow et al.2013;Schaffenroth et al.2014;Kupfer et al.2015;Almeida et al.2017;Ratzloff et al.2019a;Koen 2019;Ratzloff et al.2020;Sahoo et al.2020a;Baran et al.2021).The orbital periods of HW Vir-type stars are easily obtained from LCs,while the physical parameters such as mass and radius of sdB and the companion stars can be determined based on radial velocity measurements and LC analysis.Based on these parameters,it is possible to constrain the evolution of HW Vir-type stars and the mass–radius relations of two components.The mass–radius relation of brown dwarfs was studied by Baraffe et al.(2003) using the COND model from Chabrier et al.(2000).The mass–radius relation of low-mass stars was studied by Knigge et al.(2011) using cataclysmic variables(CVs).HW Vir-type stars are suitable objects for updating the physical properties of low-mass stars and brown dwarfs.

    There are currently about 200 HW Vir-type stars in the database of AAVSO International Variable Star Index4https://www.aavso.org/vsx/index.php(Watson et al.2006,VSX),which are mainly collected from the periodic variable catalog of the Zwicky Transient Facility (Chen et al.2020,ZTF),the fourth phase of the Optical Gravitational Lensing Experiment(Udalski et al.2015,OGLE)and the catalog of the Asteroid Terrestrial-impact Last Alert System (Heinze et al.2018,ATLAS).Schaffenroth et al.(2019)studied a sample of HW Vir-type stars in OGLE and ATLAS.ZTF is a northernsky time-domain survey with limiting magnitudes of 20.8 mag in g band and 20.6 mag in r band (Masci et al.2019).By monitoring the northern sky with several hundred detections for each object over a 3 yr period,ZTF will discover a large number of HW Vir-type stars.The presence of more short-period(P <0.2 days) HW Vir-type stars will help constrain the evolutionary timescale and ending of such systems.

    In this paper,we analyze the physical and orbital parameters of 29 sdB+dM eclipsing binaries using LCs from ZTF Data Release 5(DR5).By fixing an sdB’s mass and temperature,we obtain the orbital inclination,mass ratio,radius of sdB and dM,and the dM temperature.We also determine the mass–radius relation of dM and investigate how an sdB+dM evolves with the shortening of the orbital period.In Section 2,we describe the selection of sdB+dM eclipsing binaries and their LCs.We explain how to use PHOEBE to obtain orbital parameter solutions from LCs,and the reliability of our determined parameters in Section 3.In Section 4,we present our results and make a comparison with previous works.We also discuss the physical properties and the evolutionary stage of our sdB+dM eclipsing binaries in this section.Finally,we conclude this work in Section 5.

    2.Sample and Data Selection

    We cross-matched the HW Vir-type stars in the VSX database with the ZTF catalog of periodic variable stars and obtained a sample of 31 HW Vir-type stars.Two other HW Vir-type stars not included in the VSX were identified by eye and added to the sample.About two-thirds of this sample have periods smaller than 0.2 days.We compared them to~900 other short-period variables from the ZTF Catalog of Periodic Variable Stars through the Gaia-band intrinsic color versus absolute magnitude diagram (CMD).This helps determine the type of primary component for HW Vir-type stars.In Figure 1,the red dots represent the 33 HW Vir-type stars,and the gray dots are the ZTF short-period variables with P <0.2 days.We adopted the extinction values of Green’s 3D extinction map (Green et al.2019) and converted them to AGand E(GBp-GRp) using the updated extinction law (Wang &Chen 2019).The absolute magnitude and intrinsic color are estimated by MG=G-AG-5respectively.? is the corrected Gaia Early Data Release 3(EDR3) parallax in mas (Lindegren et al.2021).The corrected EDR3 parallax was found to reduce the systematic bias to <10 μas (Ren et al.2021).

    Figure 1.Gaia-band absolute magnitude vs.intrinsic color diagram.The red dots represent 33 HW Vir-type stars,and the error bars of four brighter HW Vir-type stars are added.The gray dots signify the locations of ZTF shortperiod variables with P <0.2 days and distance uncertainty less than 50%.

    From Figure 1,we found that HW Vir-type stars concentrate as a clump with an absolute magnitude around MG=5 mag.The clump is bluer than the main sequence variables and brighter than the cataclysmic variables.This means our HW Vir-type stars are sdB+dM eclipsing binaries rather than WD+dM eclipsing binaries.To further confirm that primary components in our sample are all sdB stars,we adopt the method in Gentile Fusillo et al.(2015),where the authors calculated the reduced proper motion defined as H=G +5 logμ+5.G and μ are G band magnitude and proper motion respectively.We found that H values of our sample are all less than 14.5,which demonstrates that they are sdBs rather thanWDs (Schaffenroth et al.2019).The distribution of HW Virtype stars is consistent with constraints of-1.0 <MG<7.0 and(Geier et al.2019;Geier 2020),but more concentrated in MG.Four sdB+dMs are brighter than the others,and these deviations are due to the large uncertainties propagated from Gaia parallaxes.We display their error bars in Figure 1.The 1σ uncertainties of the absolute magnitude for other sdB+dMs are less than 0.35 mag.We also note that in the CMD,there are dozens of gray dots distributed in the location of our sdB+dM eclipsing binaries.By examining their LCs,we found that they are non-eclipsing sdB+dM binaries,and their brightness variations are due to reflection.Reflection occurs when the hot component illuminates part of the spherical surface of the cooler component when the temperature difference between the two components is large.The size of the illuminated surface varies with the projection angle,which leads to an eventual sinusoidal-like LC without eclipse.This suggests that many sdB binaries orbit at such a small inclination that we cannot see the eclipse.The HW Vir-type star is a subtype of detached eclipsing binaries (EA-type) that have similar LC shapes.The specific distribution in the intrinsic color versus absolute magnitude diagram helps to separate HW Vir-type stars from other main sequence eclipsing binaries.With these properties,we expect to find hundreds of HW Vir-type stars in future ZTF-based variable star searches.Moreover,the CMD is useful to separate sdB+dM eclipsing binaries and WD+dM eclipsing binaries when Gaia parallaxes are accurate enough.

    Figure 2.MCMC results for ZTF J204046.56+340702.8.In the diagram,i,R1,R2 and q present a Gaussian-like distribution.

    Figure 3.MCMC results for ZTF J224547.36+490824.7.In the diagram,i,R1 and R2 present a Gaussian-like distribution,but q does not obey a Gaussian distribution.

    Figure 4.LC analysis of ZTF J204046.56+340702.8.The top and bottom panels are the LC fitting diagram and residual diagram respectively.In the top panel,the green and yellow lines are modeled LCs,and the blue and black dots are photometric measurements in g and r bands,respectively.The error bars represent the uncertainties of the photometry.In the bottom panel,the blue and black dots are the residuals between the modeled and observed LCs respectively.The green and yellow lines indicate the average of the residuals.

    We downloaded the g,r band LCs of 33 sdB+dM eclipsing binaries from the ZTF DR5 database.To ensure the quality of LCs,we selected only good photometric data with catflag <10’.ZTF DR5 contains photometric data obtained during a twoand-a-half year survey,with more than 200 detections per band for most objects.Given that typical photometric uncertainties are around 0.015 mag,we believe that a reliable analysis can be performed on the ZTF LCs.LCs were folded with periods from Chen et al.(2020),and the ephemeris of the primary minimum T0was estimated by the lc_geometry package in PHOEBE 2.3(Conroy et al.2020).We also performed a careful visual inspection of all LCs to make sure that the period and T0were not significantly biased.We excluded photometric outliers(less than 10) for each LC according to the fit line of the Gaussian processes.Magnitude was converted into normalized flux by assuming the maximum flux in g,r bands is equal to 1.The primary eclipses of four sdB+dMs are fainter than the detection limit of ZTF (r >20.5 mag).We excluded them and performed LC analysis on the remaining 29 sdB+dMs.

    3.Light Curve Analysis

    sdB+dM eclipsing binaries are the best objects to study the physical properties of sdBs and dMs.However,only the radial velocity curve of the primary component is available,so the mass ratio cannot be directly determined from radial velocity curves.One way to solve this problem is to assume a typical mass of sdB,e.g.,M1=0.47 M⊙.Based on this assumption,parameters of two components can be constrained well by analysis of the single-line radial velocity curve and LCs.The shortcoming of this method is that M1cannot be studied.The other way is to constrain the mass ratio q based onloggand fix it in the subsequent analysis(Drechsel et al.2001;Heber et al.2004;Vu?kovi? et al.2007;For et al 2010;Schaffenroth et al.2014;Almeida et al.2017;Barlow et al.2013;Schaffenroth et al.2021).In this method,the accuracy of q depends on the extent to which it is constrained by LCs.When q does not converge in the iterations,the uncertainty is relatively larger.In this work,we used the PHOEBE 2.3 code to study the probability distribution of parameters for 29 sdB+dMs and to evaluate how well q is constrained by LCs.

    With some exceptions,most sdB stars have a mass around 0.47 M⊙,which is also known as the typical mass (Han et al.2002,2003;Sahoo et al.2020b).Most sdB stars have a temperature in the range of 27,000–31,000 K (Geier et al.2012).Three of 29 sdB+dMs were studied previously with spectra and an sdB’s temperature is in this range.The three sdB+dMs are ZTF J162256.66+473051.1 with a temperature T1,eff=29,000±600 K(Schaffenroth et al.2014),ZTF J153349.44+375927.8 with a temperature T1,eff=29,230±125 K (For et al.2010) and ZTF J223421.49+245657.1 with a temperature T1,eff=28,500±500 K(Almeida et al.2017) or T1,eff=28,370±80 K (?stensen et al.2007).The companion star of an sdB+dM is an M dwarf or a brown dwarf,and its temperature is around 3000 K (Menzies &Marang 1986;Schaffenroth et al.2014;Vu?kovi? et al.2014).

    Our idea is to assume that the primary components of sdB+dMs are all canonical sdB stars with a mass of 0.47±0.03 M⊙,a radius of 0.175±0.025 R⊙and temperature about 30,000 K(Van Grootel et al.2021).We fixed the sdB mass M1to 0.47 M⊙and the sdB temperature T1to 30,000 K in our LC analysis.Small deviations in T1will directly affect the poorly constrained T2(typically with 30% uncertainty),but have little effect on R2and other parameters.With these assumptions,we can still study the properties of the companion star and the evolution stage of the system in terms of statistics.We set the orbital eccentricity e to 0 since all of our sdB+dM eclipsing binaries have a very short period (less than 0.5 days).We set the linear limb-darkening coefficient of sdB to 0.25 and 0.2 in the g band and r band,respectively,while that for the companion is 1.0 in both the g band and r band (Claret &Bloemen 2011).We set the gravity-darkening factor to 1.0 and 0.32 for sdB and the companion,respectively(Claret&Bloemen 2011).The reflection factor of sdB and the companion was set to 1.0.The range of the initial orbital inclination i is 60-90 degrees,and the other parameters are unrestricted.All preset parameters are listed in Table 1.

    Table 1 Preset MCMC Parameters in PHOEBE 2.3.

    We used the Markov chain Monte Carlo (MCMC) sampler based on EMCEE (Foreman-Mackey et al.2013) to determine parameters of 29 sdB+dM eclipsing binaries.For sampling,we relied on 50 walkers and 2000 iterations,and each run cost 30 hours on a 4-core CPU.We also tested 5000 iterations for five sdB+dMs and 20,000 iterations for one sdB+dM and found the results were almost the same as those based on 2000 iterations.In the first run,the mass ratios qiare not well constrained for most of the objects,as it is difficult to obtain the true q from many local optimal mass ratios.In contrast,the radii R1,i,R2,iand inclination iiare constrained well in two-band LC analysis.From the R1,i-q and R2,i-q distributions,R1,iand R2,ionly increase by 20%when q increases from 0 to 1,so we decided to use R2,ito establish a prior distribution of q.According to the mass–radius relationship of dM,the ratio of mass to radius is slightly less than 1 in solar units.Given that q and R2are usually overestimated in the first run(q of the sdB+dM is usually less than 0.3),we assumed a prior distribution of q as [0,R2,i/0.47].We again performed MCMC estimation and obtained reliable parameter distributions for most of the sdB+dMs.Eleven sdB+dM eclipsing binaries that have temperatures T2?3000 K are not consistent with their masses and radii.For these objects,we added another prior T2<4000 K and performed the MCMC estimation again to obtain the final distributions of parameters.

    Figure 2 displays the parameter distributions of ZTF J204046.56+340702.8 as an example of one class.We can see that i,q,R1and R2obey Gaussian-like distributions and the correlations between these parameters are small.Parameters of the other three sdB+dM eclipsing binaries are similarly distributed(see Figures A1–A3 in Appendix A).Figure 3 depicts the parameter distributions of ZTF J224547.36+490824.7.It differs from Figure 2 in that q does not satisfy a Gaussian distribution.For this class,q is uniformly distributed,or has a wide tail,or is biased toward one side in the prior range(see Figures B1–B24 in Appendix B).Considering that the prior range is narrow,i.e.,qiis roughly distributed between 0 and 0.2–0.5,the deviation between the maximum likelihood and the true q is not large.This deviation is already included in the error of q(mean uncertainty is 0.08).In addition,the error in R2hardly increases with the error in q due to the weak correlation between R2and q.Based on multi-band LC analysis,the determined masses and radii of the companions in sdB+dM eclipsing binaries are reliable for studying their statistical properties.Nevertheless,in most cases,uncertainty in q or M2is larger.

    Figure 5.Similar to Figure 4,but for ZTF J224547.36+490824.7.

    In both Figures 2 and 3,T2does not exhibit a Gaussian distribution,which means that T2is difficult to be constrained in the LC analysis.The typical uncertainty of 1500 K also indicates that T2is rather uncertain.The upper panels of Figures 4 and 5 feature a comparison of the modeled and observed LCs of ZTF J204046.56+340702.8 and ZTF J224547.36+490824.7,while the lower panels show the residual diagrams.Diagrams for other sdB+dM eclipsing binaries are supplemented in Appendix A and B.

    Figure 6.The relationship between the mean values and the residuals when q is fixed to the maximum error and minimum value.The top panel is for q.The leftbottom,middle-bottom and right-bottom panels are similar to the top but for i,r1 and r2,respectively.

    Figure 7.MCMC results for ZTFJ011338.79+431154.9 when fixed q is set to minimum error.

    Figure 8.Similar to Figure 7,but fixed q is set to maximum error.

    4.Results and Discussions

    In this section,we present a table including parameters of 29 sdB+dM eclipsing binaries and make a comparison with previous works.We also discuss how to use our result to constrain the evolutionary stage of sdB+dM eclipsing binaries and the physical properties of companion dM stars.

    4.1.Results of 29 sdB+dM Eclipsing Binaries

    The results of our sample,containing 29 sdB+dMs,were included in Table 2.It contains the ZTF ID,position(J2000 R.A.and decl.),period,physical parameters and the corresponding uncertainties.The physical parameters include i,q,R1,R2,T2,L1and L2.The bolometric luminosities of the sdB and companion are calculated using the PHOEBE code,converted from the passband luminosities.The uncertainty of luminosity is simply estimated based on the propagation uncertainties of the radius and temperature through the equationThe mean radius and luminosity of sdB are R1=0.188±0.018 M⊙and L1=21±4 L⊙respectively,which are consistent with the typical radius and luminosity of sdB (Heber 2016).

    4.2.External Errors of Orbital Parameters

    Among the obtained orbital parameters,the mass ratio q has the largest error.To check whether the accuracy of the other parameters is affected by q,we fix q to the maximum and minimum valuessee Section 4.1and Table 2) and performed the PHOEBE analysis again,respectively.Figures 7 and 8 display the MCMC results when q is fixed to the minimum and maximum values using ZTFJ011338.79+431154.9 as an example.The residuals of i,r1and r2were estimated by the difference between the new values and the mean values in Table 2.In Figure 6,the top panel shows the relationship between the mean q and the residual values of q,while the left-bottom,middle-bottom and right-bottom panels are similar,but for i,r1and r2respectively.The average percentage error of q is 32.8%,but only 3.1% and 1.4% for r1and r2respectively.These imply that for the sdB+dM system,the radii are reliable even if q has a large error.We also note that i,r1and r2all exhibit Gaussianlike distributions for a fixed q in Figures 7 and 8.

    Table 3 ZTF J184847.05+115720.3 as an Example to Test the Effect of Different Initial Values

    The selection of initial values in PHOEBE might lead to a bias in the solution of the orbital parameters.Our default initial values of q,r1and r2are 0.3,0.175 and 0.15 respectively,where q obeys a uniform distribution and r1and r2obey Gaussian distributions with a standard deviation of 0.05.To investigate the effect of the initial values on the final results,we also adopted unreasonable initial values of q=0.1,r1=0.05 and r2=0.05.From Table 3,we can see that the choice of the initial value has little effect on our results.

    4.3.Comparison with Previous Work

    Three of our sdB+dM eclipsing binaries were studied in previous works,which used both the single-line radial velocity curve and LCs to determine parameters.We compared our results with theirs to validate our method.As shown in Table 4,parameters of ZTF J162256.66+473051.1 (PG 1621+476) and ZTF J153349.44+375927.8 (FBS 1531+381) agree well with Schaffenroth et al.(2014) and For et al (2010),respectively.For et al(2010)yielded a relatively small mass for the sdB star,but it does not conflict with the typical mass of an sdB if the error is taken into account.We also confirmed that using a temperature of 30,000 K or temperatures from the literature hardly affect our parameter determination.

    As displayed in Table 5,the case of ZTF J223421.49+245657.1 is more complicated.?stensen et al.(2007)suggested that the primary component is an sdB star with a mass of 0.47 M⊙or 0.499 M⊙,while Almeida et al.(2017) suggested that it is a low-mass WD with a mass of 0.19 M⊙or 0.288 M⊙.Our q agrees with these two works while R1and R2do not.It is worth noting that the radii of ?stensen et al.(2007) were given in units of orbital separation not solar radius.Our radius ratio is consistent with Almeida et al.(2017),which means that the reason for the different parameters is due to the difference in the adoption of M1.We prefer the primary component to be a typical sdB when considering its location on the CMD.

    Table 4 Physical Parameters of ZTF J162256.66+473051.1 and ZTF J153349.44+375927.8 Compared to Previous Works.

    Table 5 Physical Parameters of ZTF J223421.49+245657.1 Compared to Previous Works

    In summary,when the primary star is a typical sdB,parameters i,q,R1,R2and M2obtained from only LCs are consistent with those obtained from both LCs and the singleline radial velocity curve.With the single-line radial velocity curve,q and M2can be determined with better precision,especially when the error of the mass ratio based on LCs alone is large.According to the discussion in Section 3 and this section,the mass estimates of sdB stars based on the single-line radial velocity curve have a large uncertainty if q is not constrained well by LCs.

    4.4.Evolutionary Stage of sdB+dM

    Figure 9 depicts the evolutionary diagram of our sdB+dM eclipsing binaries.In the top panel of Figure 9,the mass ratios are randomly distributed in q=0.1-0.6 when the orbital periods are longer than 0.1 days.However,we only found low mass-ratio(q~0.2) sdB+dM eclipsing binaries when P <0.08 days.In the bottom panel of Figure 6,the decrease of R2with the orbital period is more significant than the decrease of q,given the much smaller uncertainty of R2.According to Figure 6,we ascertained that in sdB+dM,the companion stars are found to be less massive as the orbital period shortens.This implies that sdB+dM with more massive companions merge earlier in the orbital decay process.When P <0.075 days,the companions of all three sdB+dMs in this work are brown dwarfs.Geier(2015)suggested that the orbital period P <~0.2 d and the minimum mass of companion stars M2<~0.06 M⊙is a critical region where no companion stars would survive in the common envelope phase.Instead they may merge with the sdB or be evaporated(Soker 1998).In the top panel of Figure 9,the blue dashed lines are the period and mass limits of the companion stars,and the absence of stars in this region suggests that our results support Geier’s arguments.More sdB+dM eclipsing binaries with periods shorter than 1 hour will be detected in future data releases of ZTF,and we will gain more knowledge about the evolutionary end of sdB+dMs.

    4.5.Physical Properties of Companion Stars

    Figure 10 shows the mass–radius relation of the companion stars.The blue,yellow and green solid lines are the theoretical mass–radius relation for brown dwarfs with ages of 1 Gyr,5 Gyr and 10 Gyr (Baraffe et al.2003).The black solid line is mass–radius relation for M dwarfs in CVs from Knigge et al.(2011),and the black dot–dashed line represents the theoretical mass–radius relation for single M-type stars from Baraffe et al.(1998).Most companion stars in our sample agree with the mass–radius relations.The three exceptions are due to underestimation of the mass ratio error or not being a genuine sdB+dM that hosts a 0.47 M⊙sdB.A gap was found between R2=0.12 R⊙and R2=0.15 R⊙and this gap is likely the boundary of low-mass stars and brown dwarfs.Our comparison reveals that the dM mass distribution in our sample is similar to that of Kupfer et al.(2015).Besides,the absolute magnitude and intrinsic color of five companion stars with R2<~0.12 R⊙are uncharacteristic.These suggest that the gap is not the result of a selection effect.In the future,based on hundreds of sdB+dM eclipsing binaries,it will be possible to update the mass–radius relation for low-mass stars and brown dwarfs.This will help to study the minimum mass of lowmass stars,and the maximum mass of brown dwarfs.

    Figure 9.Evolutionary diagram of 29 sdB+dMs.The distribution of mass ratio with period is in the top panel,while the distribution of the companion’s radius R2 with period in the bottom panel.Our sdB+dMs are shown as red dots.The blue dashed lines signify the period and mass limits of companion stars.

    Figure 10.Diagram of the mass–radius relation for the companion stars.The blue,orange and green solid lines signify theoretical mass–radius relations for brown dwarfs with ages of 1 Gyr,5 Gyr and 10 Gyr respectively.The black dotted line is theoretical mass–radius relation for M dwarfs in CVs.The black dot–dashed line represents the theoretical mass–radius relation for single M-type stars.Our 29 sdB+dMs are shown as red dots.

    5.Conclusion

    In this work,we have selected a sample of 33 HW Vir-type stars in ZTF DR5.Based on Gaia EDR3 parallax and extinction correction,we found that these HW Vir-type stars are concentrated in a clump in the intrinsic CMD.Their locations in the CMD imply that they are sdB+dM eclipsing binaries.By fixing the mass and temperature of sdB to M1=0.47 M⊙and T1=30,000 K respectively and setting the prior of q and T2,we analyzed LCs using the PHOEBE 2.3 code to obtain the probability distributions of parameters q,R1,R2,i and T2.We obtained LC solutions for 29 sdB+dM eclipsing binaries with full primary eclipse detection.R1,R and i obey Gaussian-like distributions and have little correlation with other parameters,which means that they are well constrained by the LC analysis.q does not show a Gaussian distribution in most cases,and the mean uncertainty is 0.08.

    Our parameters for three sdB+dMs are comparable with previous works that used both LCs and the single-line radial velocity curve if the mass of sdB is M1=0.47 M⊙.This means that parameters of sdB+dMs determined from LCs are suitable for statistical analysis.Based on 29 sdB+dMs,we found that both q and R2decrease with the decrease of the orbital period.sdB+dMs with larger mass companions are likely to merge early during the shortening of the orbit.It is worth mentioning that companions of all three sdB+dMs are brown dwarfs when the orbital period is less than 0.075 days.The masses and radii of the companion stars are consistent with the mass–radius relation for low-mass stars and brown dwarfs.We found a gap between R=0.12 R⊙and R=0.15 R⊙which can be explained as the boundary between low-mass stars and brown dwarfs.

    sdB+dM eclipsing binaries are important objects to study the nature of sdBs and dMs,and their evolutionary endings are very interesting.With deeper photometry and more detections,the ZTF project will detect and classify hundreds of shortperiod sdB+dM eclipsing binaries.Before the availability of large-scale and deep spectroscopic surveys,the statistical properties of sdBs and dMs can be obtained from the LC analysis of large samples.

    Acknowledgments

    We thank the anonymous referee for the useful comments.This work is supported by Sichuan Science and Technology Program (grant No.2020YFSY0034) and National Natural Science Foundation of China (NSFC) through the projects 12003022,12173047,11903045,12003046,and U1731111.This work is also supported by the Major Science and Technology Project of Qinghai Province 2019-ZJ-A10.

    This work has made use of PHOEBE software for the analysis of light curves.PHOEBE is funded in part by the National Science Foundation (NSF #1517474,#1909109) and the National Aeronautics and Space Administration (NASA 17-ADAP17-68).The PHOEBE project website is http://phoebeproject.org/.This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC,https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions,in particular the institutions participating in the Gaia Multilateral Agreement.This publication is based on observations obtained with the Samuel Oschin 48 inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project.Z.T.F.is supported by the National Science Foundation under grant AST-1440341 and a collaboration including Caltech,IPAC,the Weizmann Institute for Science,the Oskar Klein Center at Stockholm University,the University of Maryland,the University of Washington,Deutsches Elektronen-Synchrotron and Humboldt University,Los Alamos National Laboratories,the TANGO Consortium of Taiwan,the University of Wisconsin at Milwaukee,and Lawrence Berkeley National Laboratories.Operations are conducted by COO,IPAC,and UW.

    Appendix A.Figures for sdB+dMs with Gaussian-like Distributed mass Ratio

    Figure A1.Diagram of MCMC results and LC fitting for ZTF J072905.44-183703.4.

    Figure A2.Diagram of MCMC results and LC fitting for ZTF J184042.41+070321.9.

    Figure A3.Diagram of MCMC results and LC fitting for ZTF J192055.46+041619.5.

    Appendix B.Figures for sdB+dMs with Non-Gaussian Distributed mass Ratio

    Figure B1.Diagram of MCMC results and LC fitting for ZTF J050114.39+424741.3.

    Figure B2.Diagram of MCMC results and LC fitting for ZTF J054744.02+304732.2.

    Figure B3.Diagram of MCMC results and LC fitting for ZTF J153349.44+375927.8.

    Figure B4.Diagram of MCMC results and LC fitting for ZTF J162256.66+473051.1.

    Figure B5.Diagram of MCMC results and LC fitting for ZTF J011339.09+225739.0.

    Figure B6.Diagram of MCMC results and LC fitting for ZTF J183431.88+061056.7.

    Figure B7.Diagram of MCMC results and LC fitting for ZTF J183522.73+064247.1.

    Figure B8.Diagram of MCMC results and LC fitting for ZTF J184847.05+115720.3.

    Figure B9.Diagram of MCMC results and LC fitting for ZTF J185207.60+144547.1.

    Figure B10.Diagram of MCMC results and LC fitting for ZTF J014600.90+581420.4.

    Figure B11.Diagram of MCMC results and LC fitting for ZTF J190705.22+323216.9.

    Figure B12.Diagram of MCMC results and LC fitting for ZTF J192240.88+262415.5.

    Figure B13.Diagram of MCMC results and LC fitting for ZTF J192513.66+253025.6.

    Figure B14.Diagram of MCMC results and LC fitting for ZTF J193555.33+123754.8.

    Figure B15.Diagram of MCMC results and LC fitting for ZTF J193604.87+371017.2.

    Figure B16.Diagram of MCMC results and LC fitting for ZTF J193737.06+092638.7.

    Figure B17.Diagram of MCMC results and LC fitting for ZTF J195403.63+355700.6.

    Figure B18.Diagram of MCMC results and LC fitting for ZTF J195908.44+365041.1.

    Figure B19.Diagram of MCMC results and LC fitting for ZTF J200411.60+141150.2.

    Figure B20.Diagram of MCMC results and LC fitting for ZTF J203535.01+354405.0.

    Figure B21.Diagram of MCMC results and LC fitting for ZTF J204638.16+514735.5.

    Figure B22.Diagram of MCMC results and LC fitting for ZTF J210401.41+343636.3.

    Figure B23.Diagram of MCMC results and LC fitting for ZTF J221339.18+445155.8.

    ORCID iDs

    亚洲精品美女久久av网站| 人人澡人人妻人| 亚洲美女黄色视频免费看| 久久鲁丝午夜福利片| 老熟女久久久| 国产精品九九99| 国产日韩欧美视频二区| 晚上一个人看的免费电影| 国语对白做爰xxxⅹ性视频网站| 男女边吃奶边做爰视频| 两人在一起打扑克的视频| 国产精品国产三级专区第一集| 国产有黄有色有爽视频| 99re6热这里在线精品视频| 乱人伦中国视频| 又紧又爽又黄一区二区| 少妇猛男粗大的猛烈进出视频| 激情五月婷婷亚洲| 亚洲国产av影院在线观看| 最新的欧美精品一区二区| 最近中文字幕2019免费版| 久久久欧美国产精品| 国产免费又黄又爽又色| av在线app专区| 亚洲五月色婷婷综合| 亚洲自偷自拍图片 自拍| 国产av国产精品国产| 大码成人一级视频| 日本午夜av视频| 大香蕉久久成人网| 成人免费观看视频高清| 日韩熟女老妇一区二区性免费视频| 考比视频在线观看| 日本欧美视频一区| 精品人妻1区二区| 少妇 在线观看| 咕卡用的链子| 精品第一国产精品| 菩萨蛮人人尽说江南好唐韦庄| 精品福利永久在线观看| 国产一区二区三区综合在线观看| 在线观看人妻少妇| 精品国产超薄肉色丝袜足j| 一区二区三区乱码不卡18| 中文字幕最新亚洲高清| 777米奇影视久久| 观看av在线不卡| av网站免费在线观看视频| 五月开心婷婷网| 久久毛片免费看一区二区三区| 女人高潮潮喷娇喘18禁视频| 久久人人爽人人片av| 欧美人与性动交α欧美软件| 国精品久久久久久国模美| 2018国产大陆天天弄谢| 国产极品粉嫩免费观看在线| 1024视频免费在线观看| 大陆偷拍与自拍| 亚洲男人天堂网一区| 国产淫语在线视频| 女人被躁到高潮嗷嗷叫费观| 精品福利观看| 老司机影院成人| 首页视频小说图片口味搜索 | 这个男人来自地球电影免费观看| 亚洲人成77777在线视频| 18禁观看日本| 五月天丁香电影| 国产亚洲午夜精品一区二区久久| 国产精品亚洲av一区麻豆| 日本色播在线视频| 麻豆国产av国片精品| 欧美成人午夜精品| 99国产精品99久久久久| 欧美精品av麻豆av| 男女边吃奶边做爰视频| 桃花免费在线播放| 美女主播在线视频| 免费看十八禁软件| 久久人妻福利社区极品人妻图片 | 欧美国产精品一级二级三级| 亚洲av男天堂| 99热国产这里只有精品6| 精品久久久久久电影网| 亚洲欧洲国产日韩| 久久精品亚洲av国产电影网| 久久免费观看电影| 国产精品成人在线| 日韩av在线免费看完整版不卡| 免费人妻精品一区二区三区视频| 国产黄频视频在线观看| 亚洲精品自拍成人| 免费看十八禁软件| 日本色播在线视频| 一区二区av电影网| 两人在一起打扑克的视频| 性少妇av在线| 精品一区二区三区四区五区乱码 | 亚洲国产av新网站| 久久久精品国产亚洲av高清涩受| 国产一区亚洲一区在线观看| 国产高清videossex| 色视频在线一区二区三区| 99国产精品99久久久久| av国产久精品久网站免费入址| 欧美日韩国产mv在线观看视频| 校园人妻丝袜中文字幕| 久久午夜综合久久蜜桃| av国产精品久久久久影院| 久久影院123| svipshipincom国产片| 成年女人毛片免费观看观看9 | 久久鲁丝午夜福利片| 99re6热这里在线精品视频| 熟女av电影| 国产一区二区三区综合在线观看| 看十八女毛片水多多多| 亚洲欧美精品自产自拍| 肉色欧美久久久久久久蜜桃| 少妇裸体淫交视频免费看高清 | 美女中出高潮动态图| 久久国产精品男人的天堂亚洲| 国产一区有黄有色的免费视频| 天天躁狠狠躁夜夜躁狠狠躁| 午夜免费成人在线视频| 黑人巨大精品欧美一区二区蜜桃| 精品人妻在线不人妻| 日韩 欧美 亚洲 中文字幕| 精品一品国产午夜福利视频| 热99久久久久精品小说推荐| 国产不卡av网站在线观看| 亚洲精品一二三| 国产精品.久久久| 成在线人永久免费视频| 成在线人永久免费视频| 国产av精品麻豆| 大香蕉久久网| 大香蕉久久网| 精品国产一区二区三区久久久樱花| 少妇裸体淫交视频免费看高清 | 久久国产精品男人的天堂亚洲| 在线观看免费视频网站a站| 午夜免费观看性视频| 大香蕉久久网| 亚洲成色77777| 国产成人精品在线电影| 成在线人永久免费视频| 香蕉丝袜av| 51午夜福利影视在线观看| cao死你这个sao货| 电影成人av| 国产国语露脸激情在线看| 又粗又硬又长又爽又黄的视频| 国产精品一区二区在线不卡| 国产1区2区3区精品| 午夜91福利影院| 我要看黄色一级片免费的| 麻豆乱淫一区二区| 国产成人欧美| 一个人免费看片子| 99久久人妻综合| 欧美人与性动交α欧美软件| 免费在线观看视频国产中文字幕亚洲 | 久久狼人影院| 美女中出高潮动态图| 五月开心婷婷网| 色精品久久人妻99蜜桃| 精品一区二区三区四区五区乱码 | 日日夜夜操网爽| 男男h啪啪无遮挡| 我要看黄色一级片免费的| 婷婷色av中文字幕| 999精品在线视频| 黄色 视频免费看| 国产在线一区二区三区精| www.av在线官网国产| 久久国产精品男人的天堂亚洲| 国产精品久久久av美女十八| 99久久综合免费| 久久精品久久久久久噜噜老黄| 国产成人影院久久av| 亚洲人成77777在线视频| 亚洲精品美女久久av网站| 精品人妻在线不人妻| 美女午夜性视频免费| 亚洲专区中文字幕在线| 亚洲伊人色综图| 亚洲国产毛片av蜜桃av| 久热这里只有精品99| 国产一卡二卡三卡精品| 久久精品国产亚洲av涩爱| 少妇的丰满在线观看| 久久精品成人免费网站| 日本色播在线视频| 十分钟在线观看高清视频www| 午夜激情久久久久久久| 欧美另类一区| 国产99久久九九免费精品| 丝袜美足系列| 午夜激情av网站| 一级黄片播放器| 国产爽快片一区二区三区| tube8黄色片| 亚洲,欧美精品.| 精品人妻熟女毛片av久久网站| 捣出白浆h1v1| 国产av一区二区精品久久| 国产一卡二卡三卡精品| 国产在线视频一区二区| 51午夜福利影视在线观看| bbb黄色大片| 亚洲免费av在线视频| 国产精品免费大片| av又黄又爽大尺度在线免费看| 成年女人毛片免费观看观看9 | 精品福利永久在线观看| 久久性视频一级片| 久久精品久久久久久噜噜老黄| 99久久精品国产亚洲精品| 免费久久久久久久精品成人欧美视频| 丰满人妻熟妇乱又伦精品不卡| 久久午夜综合久久蜜桃| 91精品国产国语对白视频| 久久久久网色| 啦啦啦啦在线视频资源| 人人妻,人人澡人人爽秒播 | 五月开心婷婷网| 人体艺术视频欧美日本| 久久久久久久久久久久大奶| 久久99精品国语久久久| 国产欧美日韩综合在线一区二区| 中文精品一卡2卡3卡4更新| 丰满迷人的少妇在线观看| 欧美精品av麻豆av| 青青草视频在线视频观看| 咕卡用的链子| 午夜激情av网站| 99久久人妻综合| 在线观看一区二区三区激情| 成人亚洲欧美一区二区av| 黄色片一级片一级黄色片| 国产精品欧美亚洲77777| 91精品国产国语对白视频| 欧美 亚洲 国产 日韩一| 日韩一区二区三区影片| 日本五十路高清| 精品一区二区三卡| 成年人免费黄色播放视频| 男女下面插进去视频免费观看| 一区二区三区精品91| 99热国产这里只有精品6| 亚洲国产中文字幕在线视频| 久久 成人 亚洲| 亚洲欧美一区二区三区国产| 中文字幕制服av| 亚洲人成网站在线观看播放| 国产三级黄色录像| 日韩 欧美 亚洲 中文字幕| 50天的宝宝边吃奶边哭怎么回事| 免费在线观看影片大全网站 | 欧美人与性动交α欧美软件| 久久女婷五月综合色啪小说| 少妇精品久久久久久久| 丝袜喷水一区| 高清视频免费观看一区二区| 久久精品人人爽人人爽视色| 最新的欧美精品一区二区| 一本色道久久久久久精品综合| 久久综合国产亚洲精品| www日本在线高清视频| 日韩制服丝袜自拍偷拍| 国产欧美日韩一区二区三 | 97人妻天天添夜夜摸| av线在线观看网站| 深夜精品福利| 亚洲av国产av综合av卡| 色综合欧美亚洲国产小说| 精品卡一卡二卡四卡免费| 一区二区三区激情视频| 美女午夜性视频免费| 国产精品av久久久久免费| 别揉我奶头~嗯~啊~动态视频 | 亚洲精品国产色婷婷电影| 老司机深夜福利视频在线观看 | 晚上一个人看的免费电影| 亚洲成人国产一区在线观看 | 老司机在亚洲福利影院| 亚洲黑人精品在线| 视频区图区小说| 久久精品久久精品一区二区三区| 97人妻天天添夜夜摸| 午夜福利乱码中文字幕| 亚洲一码二码三码区别大吗| 一级毛片 在线播放| 欧美日韩国产mv在线观看视频| 久久久久久人人人人人| kizo精华| 麻豆av在线久日| 在线观看国产h片| 99热全是精品| 欧美老熟妇乱子伦牲交| 欧美亚洲 丝袜 人妻 在线| 午夜久久久在线观看| 丝袜脚勾引网站| 一区二区日韩欧美中文字幕| 母亲3免费完整高清在线观看| 日韩av免费高清视频| 我要看黄色一级片免费的| 人人妻,人人澡人人爽秒播 | 波多野结衣一区麻豆| xxxhd国产人妻xxx| 伊人久久大香线蕉亚洲五| 日韩熟女老妇一区二区性免费视频| 精品国产乱码久久久久久小说| 啦啦啦中文免费视频观看日本| 亚洲精品国产av蜜桃| 19禁男女啪啪无遮挡网站| 欧美日韩一级在线毛片| 精品卡一卡二卡四卡免费| av在线老鸭窝| 亚洲国产精品一区二区三区在线| 亚洲精品国产一区二区精华液| 777久久人妻少妇嫩草av网站| 免费少妇av软件| 亚洲黑人精品在线| 精品一区在线观看国产| 性色av一级| 日韩,欧美,国产一区二区三区| 中文欧美无线码| 高清不卡的av网站| 国产亚洲精品第一综合不卡| 成在线人永久免费视频| 久久久久国产一级毛片高清牌| 老司机在亚洲福利影院| 国产精品国产av在线观看| 一区二区三区精品91| 亚洲色图综合在线观看| 涩涩av久久男人的天堂| 国产成人精品无人区| 伦理电影免费视频| 亚洲精品国产区一区二| 看十八女毛片水多多多| 性少妇av在线| 中文字幕精品免费在线观看视频| 国产日韩欧美在线精品| 亚洲人成电影观看| 色视频在线一区二区三区| 丝袜在线中文字幕| 亚洲av电影在线进入| 国产成人免费观看mmmm| 中文乱码字字幕精品一区二区三区| 午夜免费男女啪啪视频观看| 波多野结衣一区麻豆| 一区二区三区乱码不卡18| 美女午夜性视频免费| 电影成人av| 成人18禁高潮啪啪吃奶动态图| 日日爽夜夜爽网站| 国产在线一区二区三区精| 久久人妻熟女aⅴ| 亚洲欧洲精品一区二区精品久久久| 精品国产超薄肉色丝袜足j| 午夜福利免费观看在线| 国产一区亚洲一区在线观看| 久久国产精品影院| 啦啦啦 在线观看视频| 少妇人妻久久综合中文| 99国产精品一区二区三区| 亚洲人成网站在线观看播放| 日本一区二区免费在线视频| 无限看片的www在线观看| av视频免费观看在线观看| 亚洲熟女精品中文字幕| 王馨瑶露胸无遮挡在线观看| 免费黄频网站在线观看国产| 婷婷色麻豆天堂久久| 亚洲成人手机| 亚洲精品久久成人aⅴ小说| 国产精品一国产av| av电影中文网址| 亚洲精品一二三| 欧美日韩福利视频一区二区| 欧美乱码精品一区二区三区| av欧美777| 99精品久久久久人妻精品| 国产成人a∨麻豆精品| 日本欧美视频一区| 亚洲人成电影观看| 91精品国产国语对白视频| 搡老岳熟女国产| 青草久久国产| 精品少妇内射三级| 日韩电影二区| 日本av免费视频播放| 国产成人免费无遮挡视频| 久久国产精品男人的天堂亚洲| 亚洲精品日韩在线中文字幕| 啦啦啦视频在线资源免费观看| 乱人伦中国视频| 欧美精品人与动牲交sv欧美| 青青草视频在线视频观看| 国产爽快片一区二区三区| 欧美精品一区二区大全| 亚洲国产精品一区三区| 国产爽快片一区二区三区| 在线看a的网站| 我的亚洲天堂| 无遮挡黄片免费观看| 电影成人av| 美女福利国产在线| 亚洲精品久久成人aⅴ小说| 日本一区二区免费在线视频| 美女国产高潮福利片在线看| 啦啦啦 在线观看视频| 国产男人的电影天堂91| 黄色 视频免费看| 黄频高清免费视频| 久久影院123| 国产野战对白在线观看| 超色免费av| 亚洲中文日韩欧美视频| 日韩视频在线欧美| 日本av免费视频播放| 国产精品久久久久成人av| 亚洲成人国产一区在线观看 | 欧美成人精品欧美一级黄| 亚洲欧美一区二区三区黑人| 久久ye,这里只有精品| 韩国精品一区二区三区| 天天操日日干夜夜撸| 国产成人一区二区三区免费视频网站 | 久久性视频一级片| 一区二区三区乱码不卡18| 亚洲精品一二三| 国产97色在线日韩免费| 建设人人有责人人尽责人人享有的| 国产av精品麻豆| 国产高清国产精品国产三级| www.自偷自拍.com| 老鸭窝网址在线观看| 国产免费一区二区三区四区乱码| 久久久精品免费免费高清| 操出白浆在线播放| 秋霞在线观看毛片| 看十八女毛片水多多多| √禁漫天堂资源中文www| 国产在线一区二区三区精| 国产精品国产三级专区第一集| 一区二区三区乱码不卡18| 十八禁网站网址无遮挡| 一边亲一边摸免费视频| 自线自在国产av| 一级,二级,三级黄色视频| 午夜福利视频精品| 成人午夜精彩视频在线观看| 国产麻豆69| 黄网站色视频无遮挡免费观看| 国产熟女欧美一区二区| 欧美久久黑人一区二区| 国产亚洲av片在线观看秒播厂| 日韩制服丝袜自拍偷拍| 日本欧美视频一区| 男人操女人黄网站| 97人妻天天添夜夜摸| 叶爱在线成人免费视频播放| 看免费av毛片| 国产又色又爽无遮挡免| 成人国产av品久久久| 色94色欧美一区二区| 国产高清videossex| 亚洲成人免费av在线播放| 青草久久国产| 老司机亚洲免费影院| 女人久久www免费人成看片| 亚洲国产中文字幕在线视频| 97在线人人人人妻| 视频区图区小说| 又紧又爽又黄一区二区| 免费女性裸体啪啪无遮挡网站| 高清不卡的av网站| 尾随美女入室| 国产视频首页在线观看| 一区二区日韩欧美中文字幕| 免费在线观看黄色视频的| 亚洲一区中文字幕在线| 久久精品人人爽人人爽视色| 久久久久精品国产欧美久久久 | 亚洲欧洲精品一区二区精品久久久| 亚洲少妇的诱惑av| 亚洲美女黄色视频免费看| 热99国产精品久久久久久7| 深夜精品福利| 亚洲欧洲日产国产| videos熟女内射| 中文欧美无线码| 2021少妇久久久久久久久久久| 国产欧美日韩精品亚洲av| 成人三级做爰电影| av在线app专区| 婷婷色综合大香蕉| 天堂俺去俺来也www色官网| 欧美在线一区亚洲| 亚洲自偷自拍图片 自拍| 九草在线视频观看| 国产片内射在线| 欧美黑人精品巨大| 女人精品久久久久毛片| 老司机靠b影院| 丝袜人妻中文字幕| 午夜激情av网站| 亚洲一码二码三码区别大吗| 欧美日韩国产mv在线观看视频| 99国产精品免费福利视频| 一级毛片电影观看| 欧美97在线视频| 亚洲国产av新网站| 久久精品久久久久久久性| 亚洲av欧美aⅴ国产| 国语对白做爰xxxⅹ性视频网站| 国产免费现黄频在线看| 久久国产亚洲av麻豆专区| 久久久欧美国产精品| 欧美中文综合在线视频| 国产成人系列免费观看| 日韩一卡2卡3卡4卡2021年| 悠悠久久av| 国产成人影院久久av| 九色亚洲精品在线播放| 美女脱内裤让男人舔精品视频| 黄色毛片三级朝国网站| 国产精品二区激情视频| 亚洲情色 制服丝袜| kizo精华| 91麻豆av在线| 国产日韩欧美在线精品| 黄色片一级片一级黄色片| 欧美人与性动交α欧美软件| 啦啦啦 在线观看视频| 中国国产av一级| 91字幕亚洲| 久久精品成人免费网站| 一二三四在线观看免费中文在| 久久久亚洲精品成人影院| 精品福利观看| 最近最新中文字幕大全免费视频 | 国产精品九九99| 啦啦啦视频在线资源免费观看| av有码第一页| 啦啦啦在线免费观看视频4| 欧美人与善性xxx| 精品一区二区三区av网在线观看 | 老司机午夜十八禁免费视频| 国产爽快片一区二区三区| 婷婷色综合大香蕉| 老司机靠b影院| 在线观看免费午夜福利视频| 久久久久久免费高清国产稀缺| 只有这里有精品99| 超碰97精品在线观看| 欧美成狂野欧美在线观看| 女人被躁到高潮嗷嗷叫费观| 婷婷色综合大香蕉| 成人国产一区最新在线观看 | 男人爽女人下面视频在线观看| 久久人妻福利社区极品人妻图片 | 亚洲九九香蕉| 婷婷色综合www| 亚洲欧美一区二区三区黑人| 欧美日韩精品网址| 婷婷色麻豆天堂久久| 啦啦啦啦在线视频资源| 午夜av观看不卡| 黄色 视频免费看| 免费在线观看黄色视频的| 国产深夜福利视频在线观看| 波多野结衣一区麻豆| 亚洲精品第二区| 成年人黄色毛片网站| 每晚都被弄得嗷嗷叫到高潮| 操出白浆在线播放| 丝瓜视频免费看黄片| 免费不卡黄色视频| 少妇精品久久久久久久| 久久久久国产一级毛片高清牌| 精品国产一区二区三区四区第35| 日本猛色少妇xxxxx猛交久久| 80岁老熟妇乱子伦牲交| 亚洲国产精品一区三区| www.av在线官网国产| 午夜91福利影院| 亚洲国产精品一区二区三区在线| 手机成人av网站| 婷婷色麻豆天堂久久| 国产精品久久久久久精品古装| 国产成人av教育| 男女午夜视频在线观看| xxx大片免费视频| 考比视频在线观看| 青青草视频在线视频观看| 美女高潮到喷水免费观看| 人妻人人澡人人爽人人| 99国产精品一区二区蜜桃av | 在线观看人妻少妇| 大码成人一级视频| 日韩欧美一区视频在线观看| svipshipincom国产片| 免费在线观看完整版高清| 校园人妻丝袜中文字幕| 久久女婷五月综合色啪小说| 欧美性长视频在线观看| 美女大奶头黄色视频| 久久精品久久久久久久性| 国产国语露脸激情在线看| 97在线人人人人妻| 成人三级做爰电影|