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

    Detection of Quasi-periodic Oscillations in SGR 150228213

    2023-09-03 01:37:00RunChaoChenCanMinDengXiangGaoWangZiMinZhouXingYangDaBinLinQiWangandEnWeiLiang

    Run-Chao Chen,Can-Min Deng,Xiang-Gao Wang,Zi-Min Zhou,Xing Yang,Da-Bin Lin,Qi Wang,and En-Wei Liang

    Guangxi Key Laboratory for Relativistic Astrophysics,School of Physical Science and Technology,Guangxi University,Nanning 530004,China dengcm@gxu.edu.cn,wangxg@gxu.edu.cn

    Abstract The detection of quasi-periodic oscillations(QPOs)in magnetar giant flares(GFs)has brought a new perspective to studies of the mechanism of magnetar bursts.Due to the scarcity of GFs,searching for QPOs in magnetar short bursts is reasonable.Here we report the detection of a narrow QPO at approximately 110 Hz and a wide QPO at approximately 60 Hz in the short magnetar burst SGR 150228213,with a confidence level of 3.35σ.This burst was initially attributed to 4U 0142+61 by Fermi/GBM on location,but we have not detected such QPOs in other bursts from this magnetar.We also found that there was a repeating fast radio burst associated with SGR 150228213 on location.Finally,we discuss the possible origins of SGR 150228213.

    Key words: methods: statistical– (stars:) pulsars: individual (4U 0142+61)– X-rays: bursts

    1.Introduction

    Magnetars are a class of young neutron stars that have the strongest magnetic fields in the universe found so far.They have typical magnetic fieldsB~1014G,spin periodP~2–12 s,and spin-down rate~ 10?13–10?11s s?1(Turolla et al.2015).These isolated neutron stars emit a wide array of electromagnetic radiation in radio,optical,X-ray,and gamma-ray bands by the decay of their enormous internal magnetic fields,which also yields the name “magnetar”(Duncan &Thompson 1992;Kaspi &Beloborodov 2017).Magnetars can be divided into soft gamma repeaters (SGRs)and anomalous X-ray pulsars (AXPs) judging from burst activities and other aspects.

    Bursts from magnetars can be divided into three categories:short bursts are the most common type and have typical duration ~0.1 s and peak luminosities ~1039–1041erg s?1;intermediate flares are rare events that usually last 1–40 s with peak luminosities ~1041–1043erg s?1;giant flares (GFs) are the most violent and unusual activities in magnetars and have an extremely bright hard peak lasting 0.1–0.2 s with a luminosity of 1044–1047erg s?1,which is usually followed by a long pulsating tail lasting a few hundred seconds and modulated by the magnetar spin period (Turolla et al.2015).Only four such GFs have been confirmed: GRB 790 305 from SGR 0526-66 (Mazets et al.1979;Cline et al.1980),GRB 980827 from SGR 1900+14 (Feroci et al.1999;Hurley et al.1999;Mazets et al.1999),GRB 041227 from SGR 1806-20(Cameron et al.2005;Gaensler et al.2005;Hurley et al.2005;Palmer et al.2005),and GRB 200415A (Yang et al.2020;Zhang et al.2020;Roberts et al.2021;Svinkin et al.2021).

    The associated events of SGR 1935+2154 and FRB 200428 on 2020 April 28 (Lin et al.2020b;Bochenek et al.2020;CHIME/FRB Collaboration et al.2020;Mereghetti et al.2020;Li et al.2021;Ridnaia et al.2021)had established that at least some fast radio bursts (FRBs) are produced during magnetar bursts (Lyubarsky 2014,2021;Katz 2016;Yang &Zhang 2018;Yu et al.2021),but the mechanism behind these phenomena is unclear.Starquakes have been invoked to explain the occurrence of hard X-ray bursts and FRBs from magnetars (Thompson &Duncan 1995;Wang et al.2018a).This kind of crustal oscillation would leave imprints in the form of quasi-periodic oscillations(QPOs)in the temporal profiles of magnetar bursts (Huppenkothen et al.2014b;Miller et al.2019).

    QPOs have been found during the pulsating tails and the main peak of magnetar GFs (Barat et al.1983;Israel et al.2005;Strohmayer &Watts 2005,2006;Castro-Tirado et al.2021),and have also been found in some short bursts from SGRs (Huppenkothen et al.2014a,2014c;Li et al.2022).These investigations have opened up the possibility of studying magnetars using asteroseismology(Huppenkothen et al.2013).At present,due to the scarcity of GFs,searching for QPOs from short bursts is reasonable,although the duration of short bursts would limit the minimum frequency for such a search(Huppenkothen et al.2013).In this paper we conduct a comprehensive analysis of SGR 150228213 and report the detection of a(quasi-)periodic signal in this burst.The structure of this paper is as follows.In Section 2 we describe the Bayesian framework for searching for (quasi-)periodic signals in the observed periodogram of magnetar bursts and estimating the significance.Section 3 is the periodogram analysis of SGR 150228213.We describe how to select samples and choose the appropriate time interval to conduct such analysis.We also discuss the results of (quasi-)periodic research in this section.In Section 4 we discuss possible origins of SGR 150228213 and Section 5 is a summary of this work.

    2.Methods for Periodogram Analysis

    2.1.Generate the Periodogram

    The observed periodogram analyzed in this work is based on the fast Fourier transform(FFT)of the light-curve data from the selected time interval.Powers in observed periodogram correspond to the squared Fourier transform of the data,and we make use of the stingray Python package (Huppenkothen et al.2019) to perform this conversion to get the Leahynormalized periodograms.

    A periodogram generated from a pure noise process can be seen as the conversion of a stochastic time series.It is well known that the periodogram of any stochastic time series of lengthN,denotedIj=I(fj) at Fourier frequencyfj=j/NΔT(withj=1,...,N/2),is exponentially distributed about the true spectral densitySj=S(fj):

    (Groth 1975;Leahy et al.1983;Timmer&Koenig 1995).Thus we sampled the exponential distribution corresponding to the model power to generate (see Vaughan 2010) the simulated periodograms in this work.

    2.2.Model the Periodogram

    There are two alternative approaches to modeling the periodogram: one relies on models of the original light curve to generate the periodogram,the other uses the models of the observed periodogram directly.Modeling the original light curve is based on an accurate understanding of the burst mechanism,otherwise artificial model selection would have an immeasurable impact on potential QPO detection.Owing to the unknown emission mechanism of magnetar bursts,we chose to model the observed periodogram generated from the original light curve to search for (quasi-)periodic signals in magnetar bursts.

    While modeling the observed periodogram,we made a simple but conservative assumption that all broadband powers in the periodogram are supplied by a noise process without QPO,which is the combination of red noise at low frequencies and white noise at high frequencies (Huppenkothen et al.2013).Based on this assumption,searching for(quasi-)periodic signals through periodogram research can be followed by the Bayesian approach developed by Vaughan (2010);such a method provides a statistically rigorous framework to test whether additional model components (such as Lorentzian QPOs) are required by the data.And as was stated in Castro-Tirado et al.(2021),such an assumption will cause weak signals at low frequencies to be buried in the higher variance of the broadband noise but would yield a very low false-positive detection rate in return.

    A theoretical pure red-noise profile follows a broken powerlaw model,but in many cases the break frequency is relatively low,and the red-noise profile would be fitted better by the power-law model (Belli 1992;Lazzati 2002).Therefore,we need to select the preferred noise model of the observed periodogram from these two nested models below.We defined the PL model as a red-noise power-law function plus a whitenoise (Poisson noise) constant as

    where ν is the frequency,P(ν)is the power,Ais the amplitude,α is the power-law index,andCis the constant representing the white-noise level.The BPL model is the combination of a broken power law and a white-noise constant,which is described as

    whereNis the normalization value,νbis the break frequency,and β is the power-law index after νb.

    As for the model fitting,we obtained the optimum model parameter set from the maximum a posteriori estimates,which could be computed by minimizing the maximum likelihood estimation (MLE) function (Vaughan 2010;Huppenkothen et al.2013)

    wherep(I∣θ,H)=is the joint likelihood function,Ijis the individual power in the observed periodogram,andSjis the power in the noise model for a parameter set θ.

    To select the preferred noise model,we make use of the likelihood ratio test (LRT).Thenull hypothesisis that the periodogram can be described by a simple model,PL(H0);then we estimated whether theH0model could be replaced by a more complex model,BPL (the alternative hypothesis,H1)through the LRT statistic

    We can generatensets of simulated periodograms by sampling the posterior distribution ofH0model parameters.We then compute the correspondingTLRTby fitting each fake periodogram with bothH0andH1models.The preferred noise model can be judged from the tail area probability (p-value) of the observedTLRTin the distribution of the simulatedTLRT.It is necessary to emphasize that this test cannot be seen as direct evidence in favor of theH1model (usually the more complex one),but only strictly as evidence against theH0model(Huppenkothen et al.2014c).

    2.3.Search for (Quasi-)periodic Signals

    After the selection of the preferred noise model,we use it to search for periodic signals or QPO candidates.We computed residuals of the observed power to noise model power in the logarithmic periodogram from the selected noise model with the optimum parameter set.Such a residual is equivalent toIj/Sj,for which we can use theTRstatistic to estimate the chance probability of the candidates.TRis the maximum ratio of observed to model power given by

    In this step,we generatednsets of simulated periodograms by sampling the posterior distribution of the selected noise model parameters;from each periodogram we could obtain the newTR.These statistics would be distributed as χ2and we can get thep-value ofTRby computing the tail area probability(Vaughan 2010;Huppenkothen et al.2013).

    The search for QPOs in the observed periodogram is similar to the selection of a noise model.In this step,the null hypothesis,H0,becomes that the periodogram can be well described by the selected noise model,and the alternative hypothesis,H1,model is the superposition of theH0noise model and one or several Lorentz lines to account for QPOs(Castro-Tirado et al.2021).A Lorentz line is described by(Arnaud 1996)

    whereKis the normalization factor,σ is the FWHM of the line,and νpis the centroid frequency of QPO.We can take thepvalue of theH1model obtained by LRT statistics as the significance of such QPO based on the establishment ofH0.

    3.Periodogram Analysis for SGR 150228213

    3.1.Sample Selection

    Fermi/GBM is an all-sky monitor for any burst event and covers the energy range from 8 keV to 40 MeV(Meegan et al.2009),which is suitable for the detection of short bursts from magnetars.After years of accumulation,we have collected information on 524 bursts that were classified as SGR by machine from the official website of Fermi1https://fermi.gsfc.nasa.gov;177 of them are certified from known sources and the other 347 bursts are certified from unknown sources.

    Studies of magnetar bursts based on the observation data from Fermi/GBM usually target specific magnetars for batch analysis,and especially for those active SGRs,e.g.,SGR 1935+2154 (Lin et al.2020a)and SGR 1550-5418 (Huppenkothen et al.2014c).In this work,we focused on those bursts that were certified from unknown sources and preferred those associated with known magnetars or FRBs (sources),which are more likely to originate from magnetars.Therefore,we compared the location information for the 347 bursts from unknown sources with 30 magnetars from the McGill Magnetar Catalog2https://www.physics.mcgill.ca/pulsar/magnetar/main.html(Olausen &Kaspi 2014) and 626 FRBs from the CHIME/FRB catalog3https://www.chime-frb.ca(CHIME/FRB Collaboration et al.2021),FRBCAT4https://www.frbcat.org(Petroff et al.2016),and the Transient Name Server.5https://www.wis-tns.orgAfter the comparison,except for the samples associated with SGR 1935+2154,we found only one burst,SGR 150228213,that is related to a known magnetar,AXP 4U 0142+61,and a repeating FRB source,FRB 180916,on location.

    However,4U 0142+61 is not associated with FRB 180916,but the periodogram analysis for SGR 150228213 has revealed a possible periodic or quasi-periodic signal.The later content of this section describes our periodogram analysis for SGR 150228213 and estimation of the significance of related results;two different origins will be discussed in Section 4.

    3.2.Light curve Analysis

    Considering that short bursts from magnetars usually have short duration and a soft energy spectrum,we combined the time-tagged event data files from all triggered NaI detectors(n4,n8) and rebinned the data to 2 ms time resolution to analyze the light curve in the energy range 8?100 keV.

    We useT90to describe the main part of this burst,which is the time interval within which the accumulated counts of the burst increase from 5%to 95%of the total counts(Kouveliotou et al.1993).Since the estimation ofT90will be affected by background level,and SGR 150228213 was triggered during the active phase of 4U 0142+61,which had caused the background to fluctuate greatly,we selected a relatively long time interval near the burst to estimate the average background level in order to neutralize the effects of some potentially weak bursts;this is the time intervals from ?25 s to–1 s and 1?25 s relative to the trigger timeT0.The light curve from ?0.2 to 0.2 s is shown in Figure 1,in which we also plot the accumulated net counts corresponding to the light curve.TheT90we computed is ~98 ms.

    Figure 1.Light curve of SGR 150228213.(a)The black solid line represents the light curve obtained from the combination of events data from detectors n4 and n8 in 8?100 keV.The red solid line shows the background level.(b)The black points represent the variation in accumulated counts,and the red solid lines show the 0%and 100%levels of the total accumulated counts.The two regions marked by the blue and yellow vertical dashed lines are the T90 and Bayesian block time(Tbb,2)intervals of the light curve;the time interval Tbb,1 within green vertical dashed lines is the total burst duration we select to conduct the periodogram research.

    In addition,to depict the local characteristics of the burst and select a suitable interval to conduct the periodogram research,we adopted the Bayesian blocks algorithm described in Scargle et al.(2013)to analyze the light-curve data from ?25 s to 25 s relative toT0in the same time resolution.The accumulated net counts of each“block”are also drawn in Figure 1 in the form of a ladder graph.For the light curve in Figure 1 that includes the total duration to calculateT90of the burst,there is a long“block” covers the 95% and 100% points of the total accumulated net counts.According to the description in Yang et al.(2021),the Bayesian block duration timeTbbfor bursts from magnetar SGR 1935+2156 has a power-law trend with.Following this correlation,Tbb,2~126 ms in Figure 1 is the suitable Bayesian block duration of SGR 150228213.However,since the intervalTbb,2does not contain the main part of the burst (T90),we treat the intervalTbb,1(?80 ms to 174 ms relative toT0) as the total duration of SGR 150228213 for periodogram analysis;this contains the main part of the burst and the part that could not be distinguished as burst or background.

    3.3.Periodogram Analysis

    According to temporal analysis of SGR 150228213,we select two different time segments to compute the observed periodograms: the interval from ?80 ms to 174 ms relative toT0denotes the total duration of the burst based on Bayesian blocks,and the interval from ?44.8 ms to 72.8 ms relative toT0is the burst duration based onT90,refer to Huppenkothen et al.(2013).We combined the event data from all detectors in 8–100 keV and rebinned the light-curve data to 0.2 ms time resolution(corresponding to a Nyquist frequency of 2500 Hz).

    Referring to Huppenkothen et al.(2014c),the specific process for noise model selection and the LRT statistic is as follows.

    1.We made use of the Python package emcee (Foreman-Mackey et al.2013) to perform a suit of Markov Chain Monte Carlo simulations (MCMCs) and sampled the posterior predictive distribution of theH0model (PL)with 50 MCMC ensemble walkers and 1000 samples for each walker (containing 20% of samples in the burn-in phase for each walker).

    2.We simulated 1000 sets of periodograms from the MCMC sample of the PL model and fit each periodogram with PL and BPL models to compute the distribution ofTLRTfor those fake periodograms.

    3.If thep-value of rejecting the PL model (H0) from the observed periodogram falls below 0.05,we selected the BPL model as the preferred noise model.Otherwise,we preserve the PL model as the preferred one.After selection of the noise model,we found that the preferred noise model of both segments is PL.We then use the PL model with optimum parameter set to calculate a boundary frequency between the red-noise-dominated part and the white-noisedominated part through ν=(A/C)1/α.Divided by this boundary,we can computeTRin each part on the observed periodogram and obtain the corresponding (quasi-)periodic candidates.Using the MCMC sample of the PL model,we simulated 1000 sets of periodograms to compute the distribution ofTRin each part on the fake periodograms and then estimate the correspondingp-value of each candidate.

    Results for noise model selection and periodic research in different time segments are presented in Table 1.It can be seen that there might be a possible periodic signal or QPO candidate at ~110 Hz for each observed periodogram,which is located within the red-noise-dominated part.The signal with minimump(TR) appears in the time interval from ?80 ms to 174 ms.

    Table 1Preferred Noise Model and Potential Periodicities in SGR 150228213

    Considering that these candidates could be a narrow QPO signal at ~110 Hz,we add one Lorentz line to the PL model as a newH1model to fit the observed periodograms in each time segment.The frequency of each QPO candidate is set as the initial value of the centroid frequency of the Lorentz line,and the width of this QPO was limited within a very narrow range(less than three times the minimum frequency in each periodogram).As can be seen from Table 2,the centroid frequency of this narrow QPO we suspect in all segments is still at about 110 Hz.We then drew 5000 sets of simulated periodograms from the MCMC sample of the PL model (newH0for QPO research) to compute the distribution of LRT statistics from PL and PL+QPO models;such a QPO with the lowestp-value(the tail area fraction of)of ~0.0008 exists in the interval from ?80 ms to 174 ms,which is also consistent with the result above.

    Table 2Parameter Posteriors and Chance Probabilities for Potential QPOs in SGR 150228213

    Figure 2 is the periodogram of the observed data from ?80 ms to 174 ms relative toT0and the corresponding models with the optimum parameter sets in each step of the periodogram analysis.We noticed that there still exists a potential wide QPO signal at about 60 Hz.However,such a signal is not significant enough for the periodogram from ?44.8 ms to 72.8 ms relative toT0(Figure 3).In order to find this potential wide QPO,we continued to use the QPO model with two Lorentz lines as a newH1model to fit the observed periodograms.In this case,we no longer restrict the width parameter for the wide QPO component but still set its initial centroid frequency at 110 Hz.We still use 5000 sets of simulated periodograms generated from the MCMC sample of the PL model to compute the distribution of LRT statistics and estimate the correspondingp-value of the newH1model with a wide QPO and a narrow QPO.The fitting results corresponding to each time segment are presented in Table 2.We can see that the frequency of the narrow QPO is still ~110 Hz in the time intervals from ?80 to 174 ms and from ?44.8 to 72.8 ms,and the wide QPO component is located at approximately 60 Hz.The result with the lowestp-value of ~0.0004 still exists in the interval from?80 to 174 ms.

    Figure 2.The observed periodogram from the time interval from ?80 ms to 174 ms relative to T0.The left panel presents the diagram for the QPO model based on the assumption that the periodic signal is a potential QPO signal.The right panel presents the diagram for the model of a wide QPO at 57.54 Hz and a narrow QPO at 112.20 Hz.

    3.4.Duration of the QPOs

    To depict the variation of QPOs we discovered in the burst light curve,we employed the Lomb–Scargle method (Lomb 1976;Scargle 1982) to analyze the detrended light curve of SGR 150228213 in 8–100 keV.Considering the weak signalto-noise ratio in some untriggered detectors,we only analyzed data from n4,n8,the combination of n4+n8,and the combination of all NaI detectors.A time window of length 0.1 s was used to produce Lomb–Scargle periodograms,which were combined into a spectrogram with time step of 0.2 ms.The corresponding diagram is presented in Figure 4.The analysis for data combined from all NaI detectors is consistent with the detection of QPOs: the wide QPO at about 60 Hz appeared from about ?0.06 to 0.05 s,and the narrow QPO at about 110 Hz appeared from about ?0.05 to 0.05 s.Such QPOs are also visible in the results of n8 and the combination of n4+n8,and we can see that the most significant result exists for the single detector n8.In addition,the result for n4 presented a continuous power excess or no support for the existence of QPOs.

    Figure 4.Lomb–Scargle periodogram analysis of SGR 150228213 between ?0.2 and 0.2 s.Different panels denote analysis of the detrended light curve from the combination of n4+n8 (top left),the combination of all NaI detectors (top right),detector n4 (bottom left),and detector n8 (bottom right).The energy range is 8–100 keV and the time resolution is 0.2 ms.

    From the results in Table 2,we can see that the significance is much lower in the shorter time interval centered on this burst,while we would usually expect the opposite behavior if the QPOs were a real property of the burst.However,as can be seen from Figure 4,the most significant QPOs appeared at about ?0.05 s,which may cause the lower significance in the shorter time interval.In addition,the centroid frequencies of these two QPOs seem to have a relation of an integral multiple,which indicates that the high frequency of 110 Hz (or 120 Hz)might be the second harmonic of the 55 Hz (or 60 Hz)fundamental.

    3.5.Gaussian Process Analysis

    Gaussian processes(GPs)have been employed for searching QPOs in transient astrophysical events in recent years(Hübner et al.2022;Xiao et al.2022).The GP models QPOs as a stochastic process on top of a deterministic shape,and the latter can be understood as a mean model describing the overall trend of the burst light curve.Since the QPOs at about 60 and 110 Hz lie in the red-noise-dominated part and their confidence level was insufficient based on the noise model,we can use GPs to verify whether such QPOs are generated from the red-noise process in the time domain.

    Following the procedure described in Hübner et al.(2022),we defined the kernel function describing a QPO as

    where τ is a time constant,ais the amplitude of the oscillation,fis its frequency,andcis the inverse of the decay time of the QPO.The kernel function describing the red noise is defined as

    As the function for the mean model,since the physical mechanism of SGR 150228213 is unknown,we adopted three phenomenological models that can describe the trend of light curves for gamma-ray bursts (GRBs) or flares,i.e.,skewed Gaussians,skewed exponentials,and FRED models (Norris et al.1996;Huppenkothen et al.2015;Hübner et al.2022).The significance of the QPO can be described by the Bayes factorBFqpo,defined as

    wherekqpo+rn=kqpo(τ)+krn(τ) is the kernel function describing the QPO and the red-noise process with differentc,Z(d|kqpo+rn,μ) andZ(d|krn,μ) are the respective evidences in the QPO+red-noise and red-noise models,μ is the parameter of the mean function,anddis the data.

    In this section,we applied GPs to the light-curve data from the combination of all detectors in the time interval from ?80 to 174 ms relative toT0with 1 ms time resolution,and we made use of the publicly available code6https://github.com/MoritzThomasHuebner/QPOEstimationof GPs released by Hübner et al.(2022) to obtain the results.According to the results,the QPO is disfavored under the mean models of one FRED(BFlnqpo=?1.8),two FREDs (BFlnqpo=?0.97),and one skewed Gaussians (BFlnqpo=?0.78).And the QPO is favored under one skewed exponential (BFlnqpo=0.21),two skewed exponentials (BFlnqpo=3.02),and two skewed Gaussians(BFlnqpo=1.49).The light curve under different mean models is presented in Figure 5,and the frequency posterior is presented in Figure 6.We found that the results of the analysis based on different mean models may (or may not) be favorable to the existence of QPOs,and the skewed exponentials performed better than other models for the burst profile if the Bayes factor is used for comparisons of mean models.

    Figure 5.Gaussian process analysis of SGR 150228213 between ?80 and 174 ms.Different panels denote results using the kqpo+rn kernel and different mean models,which contain one FRED (top left),one skewed exponential (top middle),one skewed Gaussian (top right),two FREDs (bottom left),two skewed exponentials(bottom middle),and two skewed Gaussians(bottom right).In each panel,black error bars denote the total light curve with 1 ms time resolution after zero correction,the dark green line is the mean function from the maximum likelihood sample,light green lines denote 10 other samples from the posterior,and the orange line is the prediction based on the maximum likelihood sample and the 1σ confidence band.The energy range is 8–100 keV and the time resolution is 1 ms.

    In addition,the QPO frequency posterior for SGR 150228213 is constrained in all models,and the results are consistent with the detection of QPOs through frequency domain analysis.As we concluded in Section 3.4,the QPO at about 110 Hz may be a second harmonic of the 55 Hz fundamental,and such a conjecture seems also supported by the frequency distributions in Figure 6.However,since the significance of such QPOs varies under different mean models,we reserve the results of frequency domain analysis as final judgment.We can see the potential of GPs for detecting QPOs in magnetar bursts;after all,the significance based on frequency domain analysis is usually recommended under the premise of infinitely long time series.

    4.Discussion on Possible Origins of SGR 150228213

    4.1.SGR 150228213 as a Magnetar Burst from 4U 0142+61

    In the trigger report for SGR 150228213,Fermi/GBM attributed this burst to the activity of 4U 0142+61,and the location of SGR 150228213 is close to this magnetar (Roberts 2015).In addition,Swift has detected a series of hard X-ray bursts from 4U 0142+61 ~800 s before the trigger time of SGR 150228213(Barthelmy et al.2015);these bursts have also been detected by Fermi/GBM.

    4U 0142+61 is a prominent emitter in hard X-rays,optical,and infrared(Hulleman et al.2004;den Hartog et al.2008).It is the only magnetar with a debris disk but it is still debated whether it is an active gaseous one or a passive dust disk(Wang et al.2006;Ertan et al.2007).

    For typical magnetar bursts,it is not clear whether burst spectra are predominantly thermal or nonthermal (Lin et al.2011;van der Horst et al.2012).Table 3 shows the spectral fitting results for short bursts from 4U 0142+61 detected by Fermi/GBM in 2015 collected from G??ü? et al.(2017).Here we select the “preferred” model for each burst following the Bayesian information criterion (BIC,Schwarz 1978),the numerical value of which is calculated from

    Table 3Comparison of Spectral Parameters of SGR 150228213 with Bursts from 4U 0142+61a

    wherekis the number of parameters in the model,d.o.f.is the data points used in fitting,andL is the maximum likelihood.When we compare the BIC of two models,if ΔBIC <6,we consider there is no significant preference between them,but if ΔBIC>6,we prefer the model with smaller BIC (Jeffreys 1939;Mukherjee et al.1998).

    The“preferred”model parameters for each burst are marked in bold in Table 3.According to the results from fitting of the energy spectrum,SGR 150228213 is not significantly different from other bursts from 4U 0142+61.Moreover these short bursts from 4U 0142+61 detected in 2015 mostly have a harder energy spectrum than “regular” bursts from short magnetars(which usually haveEpbelow 50 keV in COMPT model fitting),which may indicate different physical origins of these bursts.Unfortunately,we did not find any QPOs similar to SGR 150228213 in other bursts from 4U 0142+61,which may be because the burst intensities are too low to provide sufficient significance for the potential QPOs.

    Combined with the relationship for the active phase of 4U 0142+61 and location,4U 0142+61 is undoubtedly the most likely origin of SGR 150228213.If this QPO signal is not a false detection,it would be the first observation of QPOs in bursts from AXPs.Considering the special feature of 4U 0142+61 itself,it may bring us a new perspective on the burst mechanism of this magnetar.

    4.2.The Relation between SGR 150228213 and FRB 180916

    Since SGR 150228213 is associated with FRB 180916 on location,we try to discuss the different origins of SGR 150228213 from a more interesting perspective.

    FRB 180916 is an active repeating FRB source with a period of ~16.35±0.15 days and a phase window of 5 days(CHIME/FRB Collaboration et al.2020).It was localized to a star-forming region in a nearby massive spiral galaxy at redshift z ~0.0337±0.0002 (Marcote et al.2020).If this connection exists,SGR 150228213 may be a short GRB event generated from a newborn magnetar,which can also explain the highly active features of FRB 180916.

    4.2.1.Spectral Analysis and the Amati Relation

    If we treat SGR 150228213 as a possible short GRB,we can use the Amati relation (Amati 2006) to check whether it is correlated with the trend of short GRBs based on analysis of its energy spectrum.In this case,we used the COMPT model and the multicolor blackbody (mBB) model to fit the energy spectrum of SGR 150228213 between 8 keV and 40 MeV,and check which model fits the burst better to compute its fluence.We extract the source spectra and background spectra,and generate the instrumental response matrix from the detectors n4,n8,and b0.All spectra are fitted using Xspec (Arnaud 1996).We use the maximum likelihood for Poisson data with Gaussian background to estimate the best-fit parameters and choose the optimum model parameters through the MCMCs.

    The COMPT model is defined as

    whereKis the normalization factor,α is the photon index,andEpis the peak energy in the νFνspectrum.The mBB model we used corresponds to diskpbb in Xspec,and it is defined as(Iyyani &Sharma 2021)

    whereKis the normalization factor,ζ is the power-law index of the radial dependence of temperature (T(r)∝r?ζ),Tpis the peak temperature in keV,andTminis the minimum temperature of the underlying blackbodies and is considered to be well below the energy range of the observed data.

    The spectrum of SGR 150228213 and the model fitting results are presented in Figure 7.According to these results,the nonthermal origin of SGR 150228213 is supported further and we can use the COMPT model fitting results to computeEγ,isoof SGR 150228213 as ~1.25×1048erg based on the redshift of FRB 180916.

    According to the Amati relation,the correlation between isotropic bolometric emission energy(Eγ,iso)and the rest-frame peak energy (Ep,z) could be written as

    whereCis around 0.8–1 andmis around 0.4–0.6.This relation is initially found in long GRBs with known redshifts,but similar relations for short GRBs have also been found in later works (Ghirlanda et al.2009;Zhang et al.2009).

    Figure 8 is theEp,z–Eγ,isodiagram of short GRBs.The position of SGR 150228213 atz=0.0337 is within the 1σ to 2σ error region of the distribution of short GRBs,and the“best” redshift range for SGR 150228213 corresponding to short GRBs isz=?+0.150.0970.35.

    Figure 8.SGR 150228213 in the Ep,z–Eγ,iso correlation diagram of short GRBs.The blue solid line denotes the relation for short GRBs,and the blue and gray dashed lines denote the 1σ and 2σ regions.The orange dashed line denotes the position of SGR 150228213 if its redshift were taken from 0.0337 to 3.Red diamonds denote the position of SGR 150228213 at z=0.0337(the redshift of FRB 180916),z=0.15(the“best”position for SGR 150228213 in current correlation for short GRBs),z=0.053,and z=0.5.Other data on short GRBs are taken from Zhang et al.(2009) and Wang et al.(2018b).The best correlation of short GRBs is taken as log (E p/erg)=(3 .24 ±0.07) +(0 .54 ± 0.04) log (E γ,iso/ 10 52erg) (Zhang et al.2018).

    4.2.2.Chance Probability

    Apart from the possibility of verifying SGR 150228213 as a short GRB from the Amati relation,we need to estimate the chance probability of the association between FRB 180916 and SGR 150228213.However,the calculation may suffer from some uncertainties.Nevertheless,we simply assume that SGR 150228213 is a candidate for a short GRB associated with FRB 180916.Following the methods in Wang et al.(2020),the chance probability of the association may be calculated as

    where λ=ρSis the number of FRBs in the regionS(≈[4 1,252.96 (1 ? cosδR)]/2).The surface number density of our FRB samples is ρ ≈626/41,252.96 ≈0.015 deg–2.For the centering angular distance from FRB 180916 to SGR 150228213 δR~ 0°.4975,one gets the chance probability~1.16%.7However,if one takes δR ~5°.45 as the radius of the position error circle(with 90% confidence) of SGR 150228213,then one gets a much higher chance probability ~24.69%.It can be seen that the chance probability of ~1%is relatively slight,which implies the possibility of association,but it is not significant.Therefore,combined with the physical analysis in the previous section,we leave open the possibility of a true association between SGR 150228213 and FRB 180916.

    5.Summary

    We applied a Bayesian framework to the observed periodogram of SGR 150228213 based on the assumption that all broadbandpower in a periodogram comes from the noise process without QPOs.We detected a narrow QPO at 112.20 Hz with a width of 5.64 Hz and a wide QPO at 57.54 Hz with a width of 26.48 Hz in SGR 150228213,with a significance level of 0.0004 (corresponding to a confidence level ?3.35σ).

    We have also discussed the possible origins of SGR 150228213,and consider that it most likely comes from the known magnetar 4U 0142+61.If it indeed comes from 4U 0142+61,this would be the first detection of QPOs in bursts from AXPs,which may lead to new insights into the physical mechanisms of magnetar bursts.However,we still do not rule out the possibility that it is a short GRB associated with FRB 180916.

    Acknowledgments

    We acknowledge the use of the Fermi archive public data.This work is supported by the National Natural Science Foundation of China (grant Nos.12203013,12273005,and U1938201),China Manned Spaced Project(CMS-CSST-2021-B11),and the Guangxi Science Foundation (grant Nos.AD22035171 and 2023GXNSFBA026030).

    免费高清视频大片| 伊人久久大香线蕉亚洲五| 婷婷亚洲欧美| 99riav亚洲国产免费| 最好的美女福利视频网| 日本免费一区二区三区高清不卡| 99久久久亚洲精品蜜臀av| 老熟妇乱子伦视频在线观看| 夜夜夜夜夜久久久久| 制服丝袜大香蕉在线| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲 国产 在线| 日韩 欧美 亚洲 中文字幕| avwww免费| 亚洲国产精品sss在线观看| 亚洲美女视频黄频| 欧美成人性av电影在线观看| 中出人妻视频一区二区| 白带黄色成豆腐渣| 草草在线视频免费看| 亚洲欧美一区二区三区黑人| 精品99又大又爽又粗少妇毛片 | 精品免费久久久久久久清纯| 又爽又黄无遮挡网站| 欧美高清成人免费视频www| 一进一出抽搐gif免费好疼| 91麻豆精品激情在线观看国产| 欧美黑人巨大hd| 夜夜夜夜夜久久久久| 99久久99久久久精品蜜桃| 嫩草影视91久久| 国产成人福利小说| 国产成人一区二区三区免费视频网站| 中文字幕av在线有码专区| 国产高潮美女av| 久久欧美精品欧美久久欧美| 国产黄a三级三级三级人| 日本免费a在线| 黄片大片在线免费观看| 国产 一区 欧美 日韩| 国产精品1区2区在线观看.| 一a级毛片在线观看| 51午夜福利影视在线观看| 欧美3d第一页| 一个人看视频在线观看www免费 | 亚洲国产精品成人综合色| 国产97色在线日韩免费| 国产伦一二天堂av在线观看| av国产免费在线观看| 国产伦人伦偷精品视频| 成人亚洲精品av一区二区| 韩国av一区二区三区四区| 神马国产精品三级电影在线观看| 黄频高清免费视频| 欧美中文综合在线视频| 精品国产亚洲在线| 九色成人免费人妻av| 后天国语完整版免费观看| 亚洲av成人精品一区久久| 亚洲熟女毛片儿| 一进一出好大好爽视频| 午夜免费激情av| 最新美女视频免费是黄的| 小说图片视频综合网站| 十八禁网站免费在线| 国产亚洲精品久久久com| 国产 一区 欧美 日韩| 亚洲电影在线观看av| 岛国视频午夜一区免费看| 一本久久中文字幕| 精品久久久久久久人妻蜜臀av| 大型黄色视频在线免费观看| 免费大片18禁| 国产激情久久老熟女| 老熟妇乱子伦视频在线观看| 国产精品,欧美在线| 男插女下体视频免费在线播放| 国产 一区 欧美 日韩| 亚洲电影在线观看av| 亚洲精品色激情综合| 午夜精品在线福利| 日韩欧美在线乱码| 亚洲成人中文字幕在线播放| www.www免费av| 久久午夜综合久久蜜桃| 三级国产精品欧美在线观看 | 特级一级黄色大片| 久久99热这里只有精品18| 欧美成人一区二区免费高清观看 | 综合色av麻豆| 色av中文字幕| 老司机在亚洲福利影院| 波多野结衣高清作品| 成人av在线播放网站| 窝窝影院91人妻| 一本久久中文字幕| 两个人看的免费小视频| 国产爱豆传媒在线观看| 一夜夜www| 精品久久久久久成人av| 久久精品91蜜桃| 亚洲精品美女久久av网站| 18禁黄网站禁片午夜丰满| 国产亚洲精品av在线| 无遮挡黄片免费观看| 日日干狠狠操夜夜爽| 国产激情久久老熟女| 九九在线视频观看精品| 亚洲精品美女久久久久99蜜臀| 18禁裸乳无遮挡免费网站照片| 真实男女啪啪啪动态图| 欧美又色又爽又黄视频| 天天躁日日操中文字幕| 日韩大尺度精品在线看网址| 日本与韩国留学比较| 最近最新中文字幕大全免费视频| 精品一区二区三区av网在线观看| 精品电影一区二区在线| 香蕉丝袜av| 欧美一级毛片孕妇| 精品一区二区三区视频在线观看免费| 国产精品自产拍在线观看55亚洲| 法律面前人人平等表现在哪些方面| 又紧又爽又黄一区二区| 综合色av麻豆| www.熟女人妻精品国产| 日韩精品青青久久久久久| 国内毛片毛片毛片毛片毛片| 成人三级做爰电影| 欧美日韩黄片免| 久久久久久国产a免费观看| 两人在一起打扑克的视频| 亚洲最大成人中文| 亚洲,欧美精品.| 日韩有码中文字幕| 日本熟妇午夜| 日本精品一区二区三区蜜桃| 老司机午夜十八禁免费视频| 精品不卡国产一区二区三区| or卡值多少钱| 黄色丝袜av网址大全| 亚洲中文字幕日韩| 久久中文字幕一级| 最新美女视频免费是黄的| 精品电影一区二区在线| 久久亚洲真实| 真人做人爱边吃奶动态| 国内揄拍国产精品人妻在线| 熟女人妻精品中文字幕| 国产精品1区2区在线观看.| 亚洲国产欧洲综合997久久,| 99久国产av精品| 在线观看免费视频日本深夜| 99国产综合亚洲精品| 亚洲avbb在线观看| 人妻夜夜爽99麻豆av| 国产精品精品国产色婷婷| 人妻久久中文字幕网| 久久香蕉精品热| 亚洲精品国产精品久久久不卡| 国产精品亚洲av一区麻豆| 黄色日韩在线| a级毛片a级免费在线| 国产v大片淫在线免费观看| 国产精品亚洲一级av第二区| 久久中文字幕人妻熟女| 黑人欧美特级aaaaaa片| 欧美av亚洲av综合av国产av| 午夜精品久久久久久毛片777| 国产午夜精品论理片| 国产高清激情床上av| 真人一进一出gif抽搐免费| 久久久久久久久中文| 18禁黄网站禁片免费观看直播| 18禁裸乳无遮挡免费网站照片| 精品久久久久久久人妻蜜臀av| www日本黄色视频网| 美女cb高潮喷水在线观看 | 老鸭窝网址在线观看| 精品一区二区三区av网在线观看| 99re在线观看精品视频| tocl精华| 亚洲自偷自拍图片 自拍| 久久天堂一区二区三区四区| 老司机午夜十八禁免费视频| 国产一区二区三区视频了| 极品教师在线免费播放| 一本精品99久久精品77| 国产精品久久久av美女十八| 成人性生交大片免费视频hd| 手机成人av网站| 亚洲av成人av| 国产精品久久电影中文字幕| 日韩欧美国产一区二区入口| 午夜福利视频1000在线观看| 婷婷亚洲欧美| av在线天堂中文字幕| 国产精品永久免费网站| 无遮挡黄片免费观看| 嫁个100分男人电影在线观看| 欧美日本视频| 亚洲精品美女久久久久99蜜臀| 狂野欧美激情性xxxx| 黑人巨大精品欧美一区二区mp4| 亚洲av日韩精品久久久久久密| 熟女少妇亚洲综合色aaa.| 国产一区二区三区视频了| 制服丝袜大香蕉在线| 亚洲人成网站高清观看| www日本黄色视频网| 麻豆国产av国片精品| 免费看光身美女| 久久天堂一区二区三区四区| 亚洲欧美精品综合久久99| 1000部很黄的大片| 性色avwww在线观看| netflix在线观看网站| 欧美午夜高清在线| 狂野欧美激情性xxxx| av在线蜜桃| 两个人视频免费观看高清| 高潮久久久久久久久久久不卡| 五月伊人婷婷丁香| 很黄的视频免费| 中文亚洲av片在线观看爽| 天天添夜夜摸| 看免费av毛片| 18禁国产床啪视频网站| 成熟少妇高潮喷水视频| 国产男靠女视频免费网站| 午夜福利视频1000在线观看| 亚洲av五月六月丁香网| 综合色av麻豆| 啦啦啦韩国在线观看视频| 国产一区在线观看成人免费| 可以在线观看的亚洲视频| 成年版毛片免费区| 51午夜福利影视在线观看| 免费观看精品视频网站| 最新中文字幕久久久久 | 999久久久精品免费观看国产| 一夜夜www| 成人欧美大片| 欧美黑人欧美精品刺激| 黄色丝袜av网址大全| 午夜影院日韩av| 国产淫片久久久久久久久 | 日本三级黄在线观看| 波多野结衣高清无吗| 丝袜人妻中文字幕| 啦啦啦观看免费观看视频高清| 国产一区二区在线观看日韩 | 精品99又大又爽又粗少妇毛片 | www.熟女人妻精品国产| 国产又色又爽无遮挡免费看| 成人鲁丝片一二三区免费| 悠悠久久av| 色精品久久人妻99蜜桃| 久久久久国产精品人妻aⅴ院| 国产午夜精品久久久久久| 欧美国产日韩亚洲一区| 国产一区二区在线观看日韩 | 国产久久久一区二区三区| 日韩欧美国产一区二区入口| 操出白浆在线播放| 黄色 视频免费看| 亚洲精华国产精华精| 久久精品亚洲精品国产色婷小说| 久久精品综合一区二区三区| 最新美女视频免费是黄的| 国产一区二区在线观看日韩 | svipshipincom国产片| 久久人人精品亚洲av| 男人和女人高潮做爰伦理| 岛国在线观看网站| 国产成人福利小说| 中文字幕久久专区| 国产精品,欧美在线| av片东京热男人的天堂| 国产乱人视频| 午夜福利视频1000在线观看| 在线观看免费视频日本深夜| 18美女黄网站色大片免费观看| 无限看片的www在线观看| 国产一区二区在线av高清观看| 高清在线国产一区| 成人三级黄色视频| 久久精品国产99精品国产亚洲性色| 一区二区三区高清视频在线| 午夜精品在线福利| 精品一区二区三区视频在线观看免费| 俺也久久电影网| 亚洲欧美日韩东京热| 国产精品爽爽va在线观看网站| 欧美日韩综合久久久久久 | 免费一级毛片在线播放高清视频| АⅤ资源中文在线天堂| 日韩大尺度精品在线看网址| 十八禁网站免费在线| www.熟女人妻精品国产| 黄色女人牲交| 国产精品久久电影中文字幕| 欧美在线黄色| 欧美一级毛片孕妇| 18禁黄网站禁片免费观看直播| 久久久久性生活片| 精品电影一区二区在线| 国产精品美女特级片免费视频播放器 | 婷婷六月久久综合丁香| 欧美成人免费av一区二区三区| av中文乱码字幕在线| 精品一区二区三区视频在线 | 国产免费av片在线观看野外av| 欧美中文日本在线观看视频| 国产精品久久久久久亚洲av鲁大| 久久这里只有精品19| 日本黄大片高清| 免费看十八禁软件| 精品福利观看| a级毛片在线看网站| 超碰成人久久| 性色av乱码一区二区三区2| netflix在线观看网站| 久久久久久国产a免费观看| 怎么达到女性高潮| 国产欧美日韩精品亚洲av| 两个人看的免费小视频| 亚洲精品色激情综合| 中文字幕最新亚洲高清| 国产91精品成人一区二区三区| 99久国产av精品| 精品久久久久久,| 亚洲成a人片在线一区二区| 男人舔女人的私密视频| 男人的好看免费观看在线视频| 免费一级毛片在线播放高清视频| 99精品欧美一区二区三区四区| 法律面前人人平等表现在哪些方面| 999精品在线视频| 亚洲av熟女| 搞女人的毛片| 91在线精品国自产拍蜜月 | 精华霜和精华液先用哪个| 色视频www国产| 国产精品影院久久| 国产亚洲欧美98| 婷婷精品国产亚洲av在线| 熟女人妻精品中文字幕| 亚洲av成人不卡在线观看播放网| 很黄的视频免费| 此物有八面人人有两片| 亚洲成a人片在线一区二区| 此物有八面人人有两片| e午夜精品久久久久久久| 国产精品久久久久久亚洲av鲁大| 桃色一区二区三区在线观看| 女人高潮潮喷娇喘18禁视频| 欧美中文日本在线观看视频| 国产不卡一卡二| 欧美成人免费av一区二区三区| 97人妻精品一区二区三区麻豆| e午夜精品久久久久久久| 亚洲精品一区av在线观看| 国产日本99.免费观看| 男人舔女人的私密视频| xxxwww97欧美| 亚洲成人精品中文字幕电影| 99久国产av精品| 国产一区二区在线观看日韩 | 国内精品久久久久精免费| 国产精品一区二区三区四区久久| 国产精品电影一区二区三区| a级毛片a级免费在线| 午夜日韩欧美国产| 黄色视频,在线免费观看| 国产熟女xx| 极品教师在线免费播放| 国产精品乱码一区二三区的特点| 欧美在线一区亚洲| 在线a可以看的网站| 色精品久久人妻99蜜桃| 欧美三级亚洲精品| 欧美乱妇无乱码| 身体一侧抽搐| 久99久视频精品免费| 亚洲专区中文字幕在线| 夜夜夜夜夜久久久久| 村上凉子中文字幕在线| 久久草成人影院| 精华霜和精华液先用哪个| 色av中文字幕| 午夜免费观看网址| 中文字幕久久专区| 国产精品爽爽va在线观看网站| netflix在线观看网站| 亚洲性夜色夜夜综合| 国产午夜精品论理片| 搡老妇女老女人老熟妇| 国产精品98久久久久久宅男小说| 欧美成人一区二区免费高清观看 | 老司机午夜十八禁免费视频| 搡老岳熟女国产| 欧美日韩亚洲国产一区二区在线观看| 一进一出好大好爽视频| 一级毛片女人18水好多| 久久久精品欧美日韩精品| 一本综合久久免费| 国产高清有码在线观看视频| 一级毛片高清免费大全| 日本精品一区二区三区蜜桃| 成年女人毛片免费观看观看9| 变态另类成人亚洲欧美熟女| 极品教师在线免费播放| 看免费av毛片| 在线观看一区二区三区| 国产视频内射| 亚洲国产高清在线一区二区三| 亚洲av熟女| 国产亚洲精品综合一区在线观看| 日本精品一区二区三区蜜桃| 欧洲精品卡2卡3卡4卡5卡区| 久久人人精品亚洲av| 操出白浆在线播放| 搡老岳熟女国产| 毛片女人毛片| 国产精品一区二区精品视频观看| 美女高潮喷水抽搐中文字幕| 亚洲,欧美精品.| 欧美黑人欧美精品刺激| 夜夜躁狠狠躁天天躁| 久久这里只有精品中国| 久久精品国产99精品国产亚洲性色| www.自偷自拍.com| 亚洲欧美日韩高清专用| 一级黄色大片毛片| 色播亚洲综合网| 欧美精品啪啪一区二区三区| 国产精品爽爽va在线观看网站| 脱女人内裤的视频| 国产精品免费一区二区三区在线| 搡老妇女老女人老熟妇| 狂野欧美激情性xxxx| 久久精品人妻少妇| 国产亚洲av嫩草精品影院| 老熟妇乱子伦视频在线观看| 精品国产乱码久久久久久男人| 国产探花在线观看一区二区| 国产91精品成人一区二区三区| 国产久久久一区二区三区| 亚洲国产精品sss在线观看| 久久久久性生活片| 曰老女人黄片| 亚洲九九香蕉| 欧美高清成人免费视频www| 变态另类丝袜制服| 国产aⅴ精品一区二区三区波| 最新中文字幕久久久久 | 国产探花在线观看一区二区| 99热只有精品国产| 在线播放国产精品三级| avwww免费| 成熟少妇高潮喷水视频| 99re在线观看精品视频| 欧美黄色淫秽网站| 国产精品久久久av美女十八| 黄色丝袜av网址大全| 国产av不卡久久| 国产伦一二天堂av在线观看| 18美女黄网站色大片免费观看| 12—13女人毛片做爰片一| 国产亚洲欧美98| 久久精品aⅴ一区二区三区四区| 91老司机精品| 床上黄色一级片| 国产伦一二天堂av在线观看| www.www免费av| 日韩免费av在线播放| 老司机午夜福利在线观看视频| 一进一出好大好爽视频| 一级毛片高清免费大全| 亚洲黑人精品在线| 国内精品美女久久久久久| 在线看三级毛片| 在线观看一区二区三区| 男女那种视频在线观看| 最近视频中文字幕2019在线8| 亚洲av五月六月丁香网| 1024手机看黄色片| 少妇人妻一区二区三区视频| 国产人伦9x9x在线观看| 久久久久久大精品| 国语自产精品视频在线第100页| 毛片女人毛片| 色在线成人网| 黑人欧美特级aaaaaa片| 免费看美女性在线毛片视频| 午夜两性在线视频| 国产视频内射| 18禁国产床啪视频网站| 国产亚洲精品综合一区在线观看| 最近最新免费中文字幕在线| 精品电影一区二区在线| 熟女少妇亚洲综合色aaa.| 亚洲国产欧洲综合997久久,| 亚洲自偷自拍图片 自拍| 18禁裸乳无遮挡免费网站照片| 97碰自拍视频| svipshipincom国产片| 国产午夜精品久久久久久| 欧美精品啪啪一区二区三区| 欧美色视频一区免费| 在线视频色国产色| 国产成人av激情在线播放| 日韩欧美免费精品| 精品国产亚洲在线| 成人鲁丝片一二三区免费| 宅男免费午夜| 欧美乱码精品一区二区三区| 久久久久国内视频| 国产精品影院久久| 伦理电影免费视频| 综合色av麻豆| 又紧又爽又黄一区二区| 日本撒尿小便嘘嘘汇集6| 巨乳人妻的诱惑在线观看| 淫妇啪啪啪对白视频| 嫩草影院精品99| 免费av不卡在线播放| 国产成人一区二区三区免费视频网站| 国产主播在线观看一区二区| 看片在线看免费视频| cao死你这个sao货| 国内少妇人妻偷人精品xxx网站 | 亚洲熟妇熟女久久| 色老头精品视频在线观看| 精品一区二区三区视频在线 | 一个人观看的视频www高清免费观看 | 精品乱码久久久久久99久播| 好看av亚洲va欧美ⅴa在| 男插女下体视频免费在线播放| 日本与韩国留学比较| 观看免费一级毛片| 国产黄色小视频在线观看| 国产一级毛片七仙女欲春2| 久久精品国产99精品国产亚洲性色| 亚洲18禁久久av| 欧美乱色亚洲激情| 后天国语完整版免费观看| 巨乳人妻的诱惑在线观看| 日本一本二区三区精品| 亚洲精品中文字幕一二三四区| 亚洲成人久久性| 99热这里只有精品一区 | 精品一区二区三区av网在线观看| 18禁美女被吸乳视频| 亚洲成人久久性| 国产成人啪精品午夜网站| 淫秽高清视频在线观看| 黄色成人免费大全| 日本撒尿小便嘘嘘汇集6| 亚洲,欧美精品.| 老司机深夜福利视频在线观看| 亚洲成av人片免费观看| 搡老岳熟女国产| 99热精品在线国产| a级毛片在线看网站| 精品人妻1区二区| 国产人伦9x9x在线观看| 在线观看美女被高潮喷水网站 | xxx96com| 18禁美女被吸乳视频| 久久久久久九九精品二区国产| 99热这里只有精品一区 | 精品久久久久久久久久免费视频| 日日干狠狠操夜夜爽| 18禁观看日本| svipshipincom国产片| 亚洲性夜色夜夜综合| 午夜福利在线观看免费完整高清在 | 亚洲狠狠婷婷综合久久图片| 日本黄色视频三级网站网址| 色综合欧美亚洲国产小说| 久久久水蜜桃国产精品网| 国产精品98久久久久久宅男小说| e午夜精品久久久久久久| 国产高清视频在线观看网站| 亚洲av五月六月丁香网| 免费在线观看亚洲国产| 中文字幕人妻丝袜一区二区| 国产毛片a区久久久久| 欧美日韩福利视频一区二区| 久久久成人免费电影| 熟妇人妻久久中文字幕3abv| 99riav亚洲国产免费| 精华霜和精华液先用哪个| 一区二区三区国产精品乱码| 美女cb高潮喷水在线观看 | 精品一区二区三区av网在线观看| 欧美在线一区亚洲| 人人妻,人人澡人人爽秒播| 天堂网av新在线| 最近最新免费中文字幕在线| 村上凉子中文字幕在线| 啦啦啦韩国在线观看视频| 18美女黄网站色大片免费观看| 最好的美女福利视频网| 高潮久久久久久久久久久不卡| 久久亚洲真实| 2021天堂中文幕一二区在线观| 欧美最黄视频在线播放免费| 日本 欧美在线| 免费看a级黄色片| 99精品在免费线老司机午夜| 性色avwww在线观看| 美女大奶头视频|