Hao-Yu Yuan, Hou-Jun Lü, Ye Li, Bin-Bin Zhang, Hui Sun, Jared Rice, Jun Yang, and En-Wei Liang

1 Guangxi Key Laboratory for Relativistic Astrophysics, Department of Physics, Guangxi University, Nanning 530004, China; lhj@gxu.edu.cn

2 Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008, China

3 Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210093, China

4 Key Laboratory of Space Astronomy and Technology, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China

5 Department of Mathematics and Physical Science, Southwestern Adventist University, Keene, TX 76059, USA

Received 2022 February 22; revised 2022 May 9; accepted 2022 May 16; published 2022 June 14

Abstract The growing observed evidence shows that the long-and short-duration gamma-ray bursts(GRBs)originate from massive star core-collapse and the merger of compact stars, respectively. GRB 201221D is a short-duration GRB lasting ～0.1 s without extended emission at high redshift z=1.046. By analyzing data observed with the Swift/BAT and Fermi/GBM, we find that a cutoff power-law model can adequately fit the spectrum with a soft Ep = keV,and isotropic energyEγ,iso=1.36×10 51erg.In order to reveal the possible physical origin of GRB 201221D, we adopted multi-wavelength criteria (e.g., Amati relation, ε-parameter, amplitude parameter,local event rate density, luminosity function, and properties of the host galaxy), and find that most of the observations of GRB 201221D favor a compact star merger origin.Moreover,we find that α?is larger than2+β? in the prompt emission phase which suggests that the emission region is possibly undergoing acceleration during the prompt emission phase with a Poynting-flux-dominated jet.

Key words: (stars:) gamma-ray burst: individual (GRB 201221D) – stars: massive – acceleration of particles

1. Introduction

Observationally, supernovae (SNe) are associated with some long GRBs (or without short GRBs; Galama et al.1998; Hjorth et al. 2003; Soderberg et al. 2004; Campana et al. 2006) and the host galaxies of long (or short) GRBs are typically associated with irregular galaxies with intense (or little) star formation (Tanvir et al. 2005; Fruchter et al. 2006).These lines of observational evidence, as well as the joint detection of the gravitational-wave event GW170817 and the short GRB 170817A (Abbott et al. 2017; Goldstein et al.2017; Savchenko et al. 2017; Zhang et al. 2018b), encourage people to believe that long and short GRBs are likely related to the deaths of massive stars (Type II) and the merger of two compact stellar objects (Type I), respectively (Eichler et al. 1989; Woosley 1993; Zhang 2006). However, some apparently long-duration (or short-duration high-z) GRBs have been suggested to originate in compact stars mergers (or massive star core-collapse). Two counterexamples are GRB 060614 (Gehrels et al. 2006), and GRB 090426 (Levesque et al. 2010;Xin et al. 2011). In particular, Zhang et al.(2021)recently discovered another peculiarly short-duration GRB 200826A which seems to originate in massive star corecollapse.

In either the death of massive star or merger of compact stars, the catastrophic event leaves behind a hyper-accreting black hole or a rapidly rotating highly magnetized neutron star(called a magnetar), which serves as the central engine of a collimated outflow (or jet) with a relativistic speed toward Earth(Usov 1992;Dai&Lu 1998b;Zhang&Mészáros 2001;Zhang 2011; Lü & Zhang 2014; Kumar & Zhang 2015; Lü et al.2015;Chen et al.2017).One basic question is what is the composition of the relativistic jets? There are two models widely discussed in the literature (Lei et al. 2013). One is the fireball model with a matter-dominated outflow which dissipates its kinetic energy in internal shocks or external shocks to produce the observed GRB emission (Rees &Meszaros 1992, 1994; Meszaros et al. 1993; Kobayashi et al.1997). Within this model, the fireball has a rapid acceleration early on and can only reduce its kinetic energy at large radii from the central engine. The other one is Poynting-fluxdominated outflow. Within this scenario, the Poynting flux energy can be converted to kinetic energy(Drenkhahn&Spruit 2002; Komissarov et al. 2009), and then converted to particle energy and radiation via magnetic dissipation, such as reconnection, current instabilities, and internal collisions(Zhang & Yan 2011). In comparison to the fireball model, a Poynting-flux-dominated jet can undergo gradual acceleration in a large range of emission region (Gao & Zhang 2015; Uhm& Zhang 2015).

Traditionally, it is assumed that the curvature effect6The curvature effect is due to the observer receiving progressively delayed emission from higher latitudes (Liang et al. 2006; Zhang et al. 2007).can be used to interpret the pulse decay(including both prompt emission and X-ray flare) if the emission region moves with a constant Lorentz factor. One can measure the decay index (α?) and the spectral index () during the pulse decay, and they satisfy a simple relationship7Throughout the paper, the notationfv (t )∝t- α?v -β?is adopted.in the Lab frame (Kumar &Panaitescu 2000).However,it is difficult to interpret the observed spectral lags of pulses with the curvature effect(Uhm&Zhang 2016b).Indeed,Uhm&Zhang(2016b)found that the dissipation of magnetic field energy in the shell via reconnection of magnetic field lines would result in a faster decrease than indicated by flux conservation if the magnetic field strength in the emitting region decreases with radius.This means that the emission region does not have a constant Lorentz factor but is accelerated, and the decay slopeα?should be steeper than2+β?(Uhm & Zhang 2015). Afterwards, evidence for rapid bulk acceleration was discovered in observations of both X-ray flares(Jia et al. 2016;Uhm&Zhang 2016a)and GRB prompt emission(Uhm&Zhang 2016b;Li&Zhang 2021).In fact,Jia et al.(2016)found that a large fraction of X-ray flares in GRBs are inconsistent with the predicted relation from the curvature effect.Li&Zhang(2021)invoked the same method to analyze the prompt emission lightcurves of single-pulse GRBs, and suggested that the emission region of at least some GRBs is undergoing acceleration during the prompt emission phase.

A bright short-duration GRB 201221D, triggered the Swift Burst Alert Telescope at 23:06:34 UT on 2020 December 21(BAT; Krimm et al. 2020) and located the source at R.A.=11h24m12sand=+ ° ′ ″decl.42 08 39 (J2000).This GRB was also detected by the Fermi Gamma-ray Burst Monitor (GBM;Hamburg et al. 2020) and Konus-Wind (Frederiks et al. 2020).Based on spectroscopy of the optical counterpart, it was measured at a redshift of 1.046 (de Ugarte Postigo et al. 2020)which is larger than 95% of short-duration GRBs without extended emission (Dichiara et al. 2020). In order to test the physical origin of the high-z short-duration GRB 201221D, we perform a comprehensive analysis of Fermi and Swift data on this burst as shown in Section 2.In Section 3,we compare some statistical relations of this burst with those of other long and short GRBs,and discuss its physical origin. In Section 4, we discuss the possible bulk acceleration in the prompt emission of this burst. A summary and conclusions are presented in Section 5.Throughout the paper, a concordance cosmology with parameters H0=71 km s?1Mpc?1, ΩM=0.30, and ΩΛ=0.70 is adopted.

2. Data Reduction and Analysis

2.1. Swift Data Reduction

GRB 201221D first triggered the Swift/BAT at 23:06:34 UT on 2020 December 21 (Page et al. 2020). The BAT data were processed using the HEASOFT package (v6.28). The light curves in different energy bands and spectra were extracted by running batbinevt (Sakamoto et al. 2008). The time bin size is fixed to 64 ms in this case due to the short duration, and the light curve shows a short-pulse with duration T90=0.15±0.04 s in the 15–350 keV (see Figure 1). The time-averaged spectrum from T0?0.06 to T0+0.17 s is best fit by a simple power-law model with spectral index 1.56±0.13 due to the narrow energy band.Moreover,we do not find any signature of extended emission even up to 100 s following the burst. The X-ray Telescope (XRT) began observing the field at 87 s after the BAT trigger, but the source is too faint to be detected with photon counting (Evans et al. 2020).

2.2. Fermi Data Reduction

We downloaded the corresponding Time-Tagged-Event data from the public data site for Fermi/GBM.8https://heasarc.gsfc.nasa.gov/FTP/fermi/data/gbm/daily/For more details on the light curve and spectra data reduction procedure see the discussion in Zhang et al. (2016). The light curves of the n8 and b1 detectors are shown in Figure 1,and the background is modeled via applying the “baseline” method (Zhang et al.2011) to a wide time interval before and after the signal and subtracting the GBM light curve. The lightcurves show a single-pulse emission with a duration of T90=0.13±0.01 in 50?300 keV. No significant signatures of precursor emission before the burst and extended emission(EE)after the burst are found in the GBM temporal analysis.

Figure 1. Swift/BAT and Fermi/GBM light curves of GRB 201221D in different energy bands with a 64 ms time bin.

We also extract both time-integrated and time-dependent spectral analyses of GRB 201221D between T0?0.03 and T0+0.1. This time interval is divided into four slices (see Table 1) based on brightness and the count statistical significance of the spectral fitting (Zhang et al. 2016, 2018a).The background spectra are extracted from the time intervals before and after the prompt emission phase and modeled with an empirical function (Zhang et al. 2011), and the spectral fitting is performed by using our automatic code “McSpecfit”in Zhang et al.(2018a).Several spectral models can be selected to test the spectral fitting of the burst,such as power-law (PL),cutoff power-law(CPL),Band function(Band),and Blackbody(BB).In order to test which model is the best fit of the data,weinvoke the Bayesian Information Criteria (BIC)9BIC is a criterion for model selection among a finite set of models,and it is defined asBIC=χ2+ k· ln(n ),where k is the number of model parameters,and n is the number of data points.The model with the lowest BIC is preferred.to judge the best model among different models. The comparison of the goodness of the fits for different models is shown in Table 2.We find that the CPL model is the best one for adequately describing the observed data. The CPL model fit of the timeintegrated spectrum is shown in Figure 2 for parameter constraints of the fit. It gives peak energyEp=113keV,and a lower energy spectral index ofΓ = -0.26.The bestfit parameters of the CPL fits are listed in Table 3. The CPL model is expressed as

Table 1 Properties of GRB 201221D

Table 2 BIC Values for Different Models We Adopted to Fit within Time-dependent Spectral Fitting

where Γ and N0are the photon index and the CPL spectral fitting normalization, respectively. To extract the time-dependent spectrum, we use a similar method to the one mentioned above. We find that the CPL model is also the best fit,and the fitting results are shown in Table 3. One can see that the tracking spectral evolution is observed during the burst (see Figure 3).

2.3. Host Galaxy and Other High-z Short GRBs

GRB 201221D was initially localized to R.A.=11h24m14 19, decl. = +42d08′35 5 with 3 9 uncertainty by Swift/XRT(Evans et al.2020).Later,with ther′ band afterglow from GTC/OSIRIS, it was localized to R.A.=11h24m14 09,decl.=+42d08′40 0 with 1″ uncertainty (Agüí Fernández et al.2021).A faint galaxy around it was identified as the host galaxy.Kilpatrick et al. (2020)analyzed stacked images of the Pan-STARRS data release (Flewelling & Alatalo 2016), and found the galaxy to have a g-band magnitude g=23.2±0.2 mag. This source was also observed by the Nordic Optical Telescope with r=23.1±0.3 mag (Malesani & Knudstrup 2020) and the Lowell Discovery Telescope with r=23.9 mag(Dichiara et al. 2020) which is consistent with emission identified in Kilpatrick et al.(2020).Recently,Agüí Fernández et al.(2021)observed the host and afterglow of GRB 201221D with GTC, and measured the redshift z=1.045±0.0008 which is consistent with the GCN report. In addition, it is detected in the Dark Energy Spectroscopic Instrument Legacy Survey (DESI/LS), and listed in the LS DR8 catalog. The DESI/LS g-band image with the localization from Swift/XRT(white dashed circle), the optical afterglow (green circle and dot), and the host galaxy (red cross) are presented in Figure 4.We analyzed the images and the host galaxy properties as follows.

Figure 2.The parameter constraints of the spectral fit with the CPL model for GRB 201221D.Histograms and contours in the corner plots show the likelihood map of constrained parameters by using our McSpecFit package. Red crosses are the best-fitting values, and pink, yellow, and green circles are the 1σ, 2σ, and 3σ uncertainties, respectively.

de Ugarte Postigo et al.(2020)obtained spectroscopy of the optical counterpart of GRB 201221D with the 10.4 m Telescope, and measured the redshift at z=1.046 based on a prominent emission feature. The high-z short GRBs play an important role in understanding the age of stellar progenitors,the cosmic chemical evolution, and formation channels of binary systems if we believe that short GRBs originate in the mergers of compact stars (Zheng & Ramirez-Ruiz 2007;Dominik et al. 2012; Anand et al. 2018). To date, more than 130 short GRBs have been detected by Swift/BAT,but less than 5% are found at z>1 (Dichiara et al. 2021).Dichiara et al. (2021) studied the properties of high-z short GRBs with z>1. They find that there are eight short GRBs with redshift z>1, with five short GRBs having EE and three short GRBs (GRBs 090426, 111117A, and 121226A) without EE. However, the short-duration GRB 090426 with z=2.609 seems to originate in a massive star core-collapse based on the properties of host galaxy and afterglow (Antonelli et al. 2009;Th?ne et al. 2011; Xin et al. 2011).

2.4. Burst Energy

Figure 3.Zooming in to the GBM lightcurve of the pulse with photon and energy flux in time interval[T0 ?0.4,T0+0.5],FRED fit(red solid line),and power-law fit (red dashed line), respectively (Top panel). The evolution of photon index (middle panel) and peak energy (bottom panel) with CPL model is presented.

Table 3 Time-dependent Spectral Fitting Results of GRB 201221D with CPL Model

3. What is the Physical Origin of GRB 201221D?

The classification of GRBs remains an open question(Zhang 2011).Our purpose is to investigate the physical origin of GRB 201221D. In this section, we will discuss the origin of GRB 201221D by comparing the properties of GRB 201221D with those of other long- and short-duration GRBs, e.g., the Amati relation (Amati et al. 2002), luminosity function, properties of host galaxy, “tip of the iceberg” effect (Lü et al. 2014), and εclassification method (Lü et al. 2010).

3.1.Comparisons of the Empirical Relationships of GRB 201221D with other Type I/II GRBs

Figure 4. The DESI/LS g-band image with the localization of Swift/XRT (white dashed circle), optical afterglow (green circle and dot) and the host galaxy (red cross) of GRB 201221D.

3.1.1. Ep–Eγ,isoRelation at the outlier of 3σ uncertainty of fits for both Type I and Type II GRBs, but seems to be closer to the distribution of Type I GRBs in comparison to Type II GRBs. It is difficult to judge the progenitor based only on this empirical correlation.

3.1.2. Local Event Rate Density

The event rate density describes how many events happen per volume per unit of time.The observed event rate density is redshift-dependent,luminosity-dependent,and beaming-dependent.Most short-duration GRBs are believed to originate in the merger of compact stars, and the observed local event rate density has a large uncertainty. Sun et al. (2015) estimated the local event rate density of short GRBs with ～(0.5–3)Gpc?3yr?1above 1050erg s?1. We estimate the local event rate density of the source fromρ0=1(VmaxT), whereVmaxis the maximum volume from which the source can be detected weighted by the redshift evolution (see Equation (3) in Sun et al.2015)and T is the total exposure time of the telescope or survey. Assuming a flux limit of Fth=10?7erg cm?2s?1, and total operation time of 16 yr for Swift, the local event rate density of GRB 201221D is ～3.0×10-3Gpc-3yr-1for the peak bolometric luminosity of ～2.09×1052erg s?1in 1–104keV. The 1σ errors are derived from Gehrels (1986)based on one detection. This value is lower than that of other Type II GRBs within the same luminosity range.A comparison of the event rate density of GRB 201221D with other Type I/II GRBs is shown in Figure 6.Moreover,we do not consider the beaming factor effect in our calculations due to uncertainty of jet opening-angle in this case.12By considering the effect of beaming factor,one can roughly estimate local event rate densities which are 25 and 500 times greater than those ignoring the beaming factor,respectively when we adopt the typical beaming factors ～0.04 and ～0.002 for short- and long-duration GRBs cases (Frail et al. 2001; Fong et al. 2015). If this is the case, the local event rate density is nearly consistent with other Type I and II GRBs within the same luminosity range.

3.1.3. Luminosity Function

The luminosity function of high-luminosity (HL) longduration GRBs can be characterized as a broken power law function via a large enough sample of redshift measurements,but the luminosity function of low-luminosity (LL) longduration GRBs is not well constrained with a small sample size due to the detectors’ sensitivity limit (Liang et al. 2007;Virgili et al. 2009; Sun et al. 2015). In comparison with longduration GRBs, the luminosity function of short-duration GRBs is less well constrained with a small fraction of the data. Sun et al. (2015) found that the luminosity function of short-duration GRBs can be roughly fitted with a simple power law in the luminosity range of [7×1049, 1050] erg s?1by assuming that all short-duration GRBs have a compact star merger origin.GRB 201221D broadly follow the distributions of both the long and short GRB populations, as shown in the right panel of Figure 6. Moreover, Paul (2018) analyzes the luminosity function with a large sample of short-duration GRBs observed by CGRO/BATSE, Swift/BAT, and Fermi/GBM. They found that the luminosity function can be described with the exponential cutoff power law and broken power law models. However, the redshifts of most short GRBs in their samples are not measured, and they adopt pseudo-redshifts which are derived by the empirical relationship. Therefore, we still used the results from Sun et al. (2015).

Figure 5. Top panel: Ep and Eiso correlation diagram. Black points and gray diamonds correspond to Type I and Type II GRBs,respectively.The red star is GRB 201221D,and other data are taken from Zhang et al.(2009).The best-fit Ep–Eiso correlations for both Type II(gray diamonds)and Type I(black points)GRBs are plotted (solid lines) with the 3σ boundary (dashed line) marked.Bottom panel: Ep distribution of GRB 201221D and other Type I/II GRBs observed by Fermi/GBM. The data of Ep values of other Type I/II GRBs are taken from Lu et al. (2017) and von Kienlin et al. (2020), respectively.

3.1.4. ε-parameter

3.1.5. Amplitude Parameter

Figure 6.Local event rate density(top)and normalized luminosity function(bottom)distribution for GRB 201221D(red),LL-LGRBs(blue),HL-LGRBs(olive),and short-duration GRBs (black) inferred Sun et al. (2015).

Figure 7. (a)1D and 2D distributions of GRBs samples in T90–ε space.The dotted line is ε=0.03.(b) f(feff)–T90 for type I and Type II GRBs which are from Lü et al. (2014). The vertical solid line is T90=2 s. Blue (black) solid circles represent the Type I(II) GRB candidates, and green triangles denote the nearby lowluminosity long GRBs. The red star is GRB 201221D.

3.1.6. Host Galaxy Properties—the Probability to Be a Type I GRB—log Odds

3.2. Merger of Compact Star Scenario

A leading progenitor for producing short-duration GRBs is binary neutron star mergers which were confirmed via the detection of the associated of gravitational-wave (GW)170817 and GRB 170817A (Abbott et al. 2017; Zhang et al. 2018b). Neutron stars likely receive a kick at birth from SNe explosions, and mergers of binary neutron stars are usually found to have a large offset from the galaxy center or even be outside of the host galaxy (Bloom et al.1999; Li & Zhang 2020). As a caution, Mandhai et al. (2021)point out that the observational offset is the minimum value due to projection effects, namely the real-offset is much larger. Most observations (amplitude parameter, properties of host galaxy, local event rate) support a compact star merger origin.

Moreover,black hole and white dwarf merger systems were also proposed as possible progenitors of short-duration GRBs(Fryer et al.1999),but these systems may not be favorable for launching a GRB jet with a small fraction of accreted mass(Narayan et al. 2001). MacFadyen et al. (2005) proposed that accretion-induced collapse of a neutron star can power a GRB,but this model can produce GRBs in non-star-forming galaxies with a small offset of GRB location. It is also inconsistent with the observations of the host galaxy of GRB 201221D.White dwarf-white dwarf mergers producing a short GRB were discussed by Lyutikov & Toonen (2017), but this model predicts the powering of long-lasting EE after the prompt emission, which is disfavored for GRB 201221D because there is no observed EE.

3.3. Massive Star Core-collapse Scenario

Figure 9.One-and two-dimensional projections of the posterior distributions of the Equation(2)free parameters by MCMC fit with pulse.Vertical dashed lines mark the the median and 1σ range. The contours are drawn at 68%, 95%, and 99% credible levels.

Figure 10. Similar with Figure 9, but adopting a power-law model to fit the decay segment of pulse.

Such a puzzle might be solved in three ways if we believe that GRB 201221D originated from massive star core collapse.One is that GRB 201221D is produced by the collapse of a supramassive NS (or magnetar) into a black hole (Rezzolla &Kumar 2015),and the supramassive NS is the initial remnant of massive star collapse due to the stiff equation of state of NS,but does not power the GRB during this process. If this is the case, it should leave clues in the X-ray band when it loses rotational energy via dipole radiation. Unfortunately, there is not enough X-ray observational data after the prompt emission of GRB 201221D. The second possibility is that the total central engine timescale(teng)is much longer than that of what we observed, but the majority of its time (tjet) is used in breaking the envelope out of the stellar surface,so the observed duration time of GRB (tGRB=teng?tjet) could be as short as 1 s(Bromberg et al.2012).However,this timescale is still one order of magnitude longer than the observed duration of GRB 201221D. The third possibility is that GRB 201221D is not intrinsically short but is actually long duration due to the tip of the iceberg effect (Lü et al. 2014). Within this scenario, the observed fluence is underestimated, which means that the intrinsic isotropic energy Eγ,isoshould be larger than the current value what we calculated in GRB 201221D.If this is the case,the location of GRB 201221D in both Ep–Eγ,isoand εparameter diagrams will move into the long GRB population.In that case, it is natural to explain the observed empirical correlations. However, the f-parameter of GRB 201221D should be consistent with that of other Type II GRBs, but the observed f-parameter of GRB 201221D does not support this hypothesis.

4. Possible Bulk Acceleration in Prompt Emission?

Uhm & Zhang (2016b) proposed that the relationship betweenα?and2+β?is not consistent with the predicted of curvature effect when the emission region itself is undergoing acceleration or deceleration in the prompt emission phase. In this section, we describe how to test whether the prompt emission of GRB 201221D is undergoing acceleration or deceleration to diagnose the composition of the jet.

As the light curve of GRB 201221D shows a single-peaked structure,we adopt a fast-rising and exponential decay(FRED)empirical function to fit in order to describe the shape of a pulse(Kocevski et al. 2003). It reads as,

Figure 11.α? –β?relation of the curvature effect during the initial decay phase of prompt emission. The solid line is = 2+β?.

5. Conclusions

Combined with the above analysis, the origin of GRB 201221D is favored to be a binary neutron star merger. This model can be used to interpret some observations of this case,except for the ε-parameter and the inconclusiveness of the Amati relation.If this is the case,the remnant of such a merger should leave either a magnetar or black hole which produces a short-duration GRB with a Poynting-flux-dominated jet.It also naturally explains undergoing bulk acceleration in the prompt emission phase. Moreover, another possibility is that GRB 201221D originates in a black hole–neutron star (BH–NS)merger at high redshift. Mandhai et al. (2021) pointed out that the probabilities of BH–NS systems for producing GRBs peak at higher redshifts (e.g., z>1), although production of shortduration GRBs is rarer with BH–NS mergers (Shibata &Taniguchi 2011).

On the other hand, the high-z short-duration GRBs without extended emission are very rare,with only five GRBs(060801,090426, 111117A, 121226A, and 190627A) observed so far.Within these five GRBs,GRB 111117A has been discussed as having a high probability of originating in a compact binary merger(Dichiara et al.2021).Moreover,Wiggins et al.(2018)estimated the rate for short GRBs peaks in the redshift range z=0.6–1 by studying the population synthesis, and this is similar to observations of short GRBs with a known redshift distribution(Fong et al.2015).Also,we expect more and more short GRBs with high-redshift to be observed in the future,and they can be used to study the star formation delay time(Wanderman&Piran 2015).If the binary neutron star mergers indeed occurred in high-z region of early universe, there are important implications for understanding binary stellar evolution, heavy element nucleosynthesis, and chemical evolution,or the possibility of the high-z short-duration GRBs belonging to a different population of bursts (e.g., NS–BH or WD–WD merger systems). However, owing to the relative lack of expected mergers in the high-z region, to understand those implications remains an open question. We hope that there are more neutron star mergers at high redshifts than that expected in the future. The third generation of GW detectors, such as Cosmic Explorer(Reitze et al.2019)or the Einstein Telescope(Punturo et al. 2010), may play a crucial role in understanding the mergers of these objects out to redshift z ～2–3.

Acknowledgments

We acknowledge the use of the public data from the Swift and Fermi data archive.This work is supported by the National Natural Science Foundation of China (grant Nos. 11922301 and 12133003), the Guangxi Science Foundation (grant Nos.2017GXNSFFA198008 and AD17129006). LHJ thanks support by the Program of Bagui Young Scholars Program, and special funding for Guangxi distinguished professors (Bagui Yingcai & Bagui Xuezhe). B.B.Z. acknowledges support by the National Key Research and Development Programs of China (2018YFA0404204), the National Natural Science Foundation of China (grant Nos. 11833003 and U2038105),and the Program for Innovative Talents, Entrepreneur in Jiangsu. H.S. acknowledges support by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant Nos. XDA15310300, XDA15052100 and XDB23040000).Y.L. acknowledges support by the National Natural Science Foundation of China (grant Nos. 12041306 and 12103089),and the Natural Science Foundation of Jiangsu Province(grant No. BK20211000).