Compact CubeSat Gamma-ray detector for GRID mission

Jia-Xing Wen ? Xu-Tao Zheng ? Jian-Dong Yu ? Yue-Peng Che ?Dong-Xin Yang ? Huai-Zhong Gao ? Yi-Fei Jin ? Xiang-Yun Long ?Yi-Hui Liu? Da-Cheng Xu? Yu-Chong Zhang ? Ming Zeng?Yang Tian ? Hua Feng1, ? Zhi Zeng? Ji-Rong Cang, ? Qiong Wu1, ?Zong-Qing Zhao ? Bin-Bin Zhang ? Peng An ? GRID collaboration

Abstract Gamma-Ray Integrated Detectors (GRID) mission is a student project designed to use multiple gammaray detectors carried by nanosatellites (CubeSats),forming a full-time all-sky gamma-ray detection network that monitors the transient gamma-ray sky in the multi-messenger astronomy era.A compact CubeSat gamma-ray detector,including its hardware and firmware,was designed and implemented for the mission.The detector employs four Gd2Al2Ga3O12 :Ce(GAGG:Ce)scintillators coupled with four silicon photomultiplier (SiPM) arrays to achieve a high gamma-ray detection efficiency between 10 keV and 2 MeV with low power and small dimensions.The first detector designed by the undergraduate student team onboard a commercial CubeSat was launched into a Sun-synchronous orbit on October 29,2018.The detector was in a normal observation state and accumulated data for approximately one month after on-orbit functional and performance tests,which were conducted in 2019.

Keywords Gamma-ray bursts · Scintillation detectors·SiPM·CubeSat

1 Introduction

Gamma-Ray Integrated Detectors (GRID) mission is a student project designed for the scientific purpose of monitoring the transient gamma-ray sky in the local universe,particularly,to accumulate a sample of gamma-ray bursts (GRBs) associated with gravitational waves (GWs).According to the estimation of the GW-GRB joint detection event rate,the maximum number of events that can be detected is more than a dozen per year.Therefore,GRID is designed to serve as full-time all-sky gamma-ray detection network,without Earth occultation or interruptions due to the South Atlantic Anomaly (SAA),with many compact and modularized gamma-ray detectors on a fleet of Cube-Sats in low Earth orbit.As a distributed system,GRID can localize detected GRBs via triangulation or flux modulation with simple scintillation detectors [1].

The scientific payloads of GRID comprise several modularized compact gamma-ray detectors.Many key technologies have been utilized to optimize the gamma-ray detection performance in the limited space of a CubeSat.In the first gamma-ray detector prototype of GRID,GAGG:Ce was used,and the package was optimized for high transmittance of low-energy X-rays down to 10 keV.To respect the power and space limitations of CubeSats,silicon photomultipliers (SiPMs) were utilized instead of traditional photomultiplier tubes (PMTs) owing to their attractive capabilities,such as their super miniature size,low weight,low power consumption,and insensitivity to magnetic fields [2].The strong dark noise of the SiPM array restricts the signal-to-noise ratio (SNR) at room temperature;hence,the current sensitive pre-amplifier was designed,modeled,and optimized to improve the SNR.Data acquisition (DAQ) electronics were designed based on an off-the-shelf ARM core microcontroller unit(MCU),which is sufficiently simple for undergraduate students.We note that numerous CubeSat-based missions similar to GRID mission have been proposed and are under development [3–8].The scientific payloads of BurstCube [4],CAMELOT [5],HERMES [7],and GRBAlpha (in-orbit demonstration for CAMELOT) [6] are scintillator-based detectors similar to those of GRID.BlackCAT [3] adopts silicon detectors and is sensitive to soft X-rays.LECX [8]is designed to employ four CdZnTe (CZT) detectors with high energy resolution.

The first detector prototype onboard a 6 U(30 cm×20 cm×10 cm) CubeSat developed by Spacety Co.Ltd,a commercial satellite company in China,was launched into a Sun-synchronous orbit [1] with an altitude of 500 km and inclination of 97.5°.It was in a normal observation state and accumulated data for approximately one month after its on-orbit functional and performance tests,which were conducted in 2019.Compared with the first detector,the second detector has fewer hardware modifications and improvements.An aluminized polyimide film was used to better insulate the sunlight,and the leakage current monitoring circuits of the SiPM arrays were modified to provide a wider measurement range in the second detector.The second detector was launched into a Sun-synchronous orbit with an altitude of 500 km and inclination of 97.3°C by the Long-March 6 rocket on November 6,2020.The second detector has accumulated data for more than 300 h of on-orbit observation.Multiple GRBs were observed.All the corresponding scientific data will be collected and published by the National Space Science Data Center (NSSDC) in the future.In this paper,we present the detailed design of the detector,electronics,and firmware.In addition,we discuss the energy resolution,low-energy X-ray detection performance,and high-rate performance of the detector.

2 Detector structure

A schematic and photograph of the first detector fabricated for GRID mission are shown in Fig.1.The detector comprises four GAGG:Ce scintillators coupled with four SiPM arrays on one SiPM board.Each SiPM array comprises 4×4 J-60035-type SiPMs from SensL.The standard output pins of the SiPMs are connected in parallel on a front-end electronics (FEE) board.The following DAQ board provides the capabilities of four-channel signal digitization,power distribution,communication,and detector control. The detector has dimensions of 5 cm×9.4 cm×9.4 cm,which occupies half of a Cube-Sat with standard dimensions (units or ‘U’),and a power consumption of 5 V×0.6 A.The detector provides a novel electronic interface,supporting the serial peripheral interface (SPI),universal asynchronous receiver/transmitter(UART) protocols,and pulses per second (PPS) interface.The features of the first detector,which is a modularized CubeSat payload suitable for different nanosatellite platforms,are summarized in Table 1.

Table 1 Main features of the first compact gamma-ray detector of GRID

Fig.1 (Color online)3D model of a GRID(top)and photograph of a GRID (bottom)

3 Scintillator and FEE

3.1 Scintillator and package

We used GAGG:Ce scintillators manufactured by EPIC CRYSTAL Co.Ltd.Owing to the limited available size of GAGG:Ce scintillator,the detector comprises four GAGG:Ce scintillators to make full use of the area of a standard nanosatellite unit.A single GAGG:Ce scintillator has a surface area of 3.8 cm×3.8 cm and a thickness of 1 cm.GAGG:Ce has the advantages of a high density(6.5 g/cm3) and a high detection efficiency for gammarays of up to the MeV magnitude.The maximum emission spectrum for GAGG:Ce is approximately 530 nm,which is suitable for silicon-based photodetectors.Its high light yield (～30–70 ph/keV) and～100 ns decay time provide a reasonable SNR when the SiPM arrays are used.The reflection layer is an enhanced specular reflector (ESR),a 65 μm polymer with high reflectance (＞98%) manufactured by 3M.A series of aluminized polyester films,which are generally used as the outermost cover of satellites,is used as the light shield layer.Owing to the good mechanical characteristics and non-hygroscopic properties of GAGG:Ce [9],no further layers are required to ensure sufficient moisture proofing and mechanical reinforcement.The simplified scintillator package provides a high detection efficiency of low-energy X-rays at a low cost.

With an increase in the number of parallel connected SiPMs on the same channel,the output capacitance and dark count rate of the SiPM array increase,which reduce the SNR.Therefore,an SiPM array smaller than the bottom area of the scintillator is utilized,similar to many detectors where the scintillator is coupled with SiPMs.The ESR is cut into a particular shape and fully covers the top and side of the scintillator.Moreover,there is a 2.2 cm×2.2 cm square window at the bottom-center for the scintillation light collection,as shown in Fig.2.By measuring the 661.7 keV full-energy peak positions of the137Cs source using different packages,we found that the light collection efficiency of this type of package is 62%that of a full coupling of the scintillator bottom surface.The optical grease is used for the optical coupling between the scintillator and SiPM array in the first payload.It is replaced with a silicon sheet in the second payload because of the fixed shape and good light collection efficiency of the silicon sheet.

Fig.2 (Color online) ESR package (left),where the solid line is the outline,and the dashed line is the crease.Bottom view of a GAGG:Ce scintillator with the ESR package (right).The ESR is wrapped directly on the scintillator and taped on the side and bottom

3.2 SiPM

The SiPM has numerous attractive features,such as small size,insensitivity to magnetic fields,low power consumption,and light weight,which are crucial in space mission applications [2],particularly nanosatellites.Therefore,four SiPM arrays are utilized as photoelectric converters.The SiPMs are the J-60035 type manufactured by SenSL and have a high photon detection efficiency(PDE) curve that matches reasonably with the emission spectrum of GAGG:Ce with a low dark count rate.Each SiPM chip has dimensions of 6.13 mm×6.13 mm×0.6 mm and consists of 22,292 single-photon avalanche diodes with a microcell fill factor of 75%.The bias voltage supplied to the SiPM is 28.5 V,which is considerably lower than that required by a PMT.To reduce the dark noise and output capacitance of the SiPM array and manage the costs,a 4×4 SiPM array is adopted with a 2.45 cm×2.45 cm area,which covers approximately onethird of the GAGG:Ce scintillator bottom surface area,as mentioned previously.

Four SiPM arrays are integrated on a single board designed and manufactured in-house (Fig.3).In addition,an SiPM through silicon via (TSV) package ensures that the SiPM fill factor of the printed circuit board (PCB)footprint is over 93%.Each array is poweredindependently,and the fast and standard outputs of every single SiPM chip are independently extracted to the FEE through a high-density connector QTE04003FDA manufactured by SamTec.Because the breakdown voltage of the SiPM changes with temperature,which affects the gain of the SiPM,there is a temperature-monitoring chip on the other side of the PCB for the correction of the SiPM gain.

Fig.3 (Color online) Top view(left)and bottom view(right)of the SiPM array board

3.3 FEE

A schematic of the FEE is shown in Fig.4.The standard output signals from one 4×4 SiPM array are directly shorted on the FEE and fed to a transimpedance amplifier(TIA) via an alternating current (AC) coupling,while the fast output signal pins are left floating.A 2 kΩ resistor connects the standard output of the array to the ground as the direct-current (DC) path of the standard output,which can restrain the current of the SiPM array to improve the system robustness.A standard high-speed amplifier OPA656 was adopted as the TIA amplifier,and the parameters of the TIA were optimized by detector modeling,as described in the next section.The TIA is followed by a low-pass filter circuit to reverse the signal and adjust its amplitude.The filter output is directed into two paths,which are connected to the trigger and peak hold circuits.The trigger circuit comprises a hysteresis comparator with an adjustable threshold and monostable pulse generator LTC6993,generating a high-level trigger signal for 2 μs without a retrigger.The high-bandwidth peak hold circuit comprises an operational transconductance amplifier(OTA) OPA615 and electronically controlled analog switch for discharge.Then,the four triggers and four peak hold signals are fed to the DAQ through a standard 2.54 mm connector.

Fig.4 Functional block diagram and output pulse shape of the FEE for one GRID channel

3.4 Detector optimization

For such a compact detector,considerable efforts have been invested to improve the SNR,such as the scintillator reflector design,SiPM array design,and high-speed amplifier application.However,we do not have flexibility in terms of the scintillator and SiPM,whereas the TIA parameters significantly affect the SNR and can be analyzed and optimized in detail.

The transient response of the GAGG:Ce scintillator can be expressed as a single exponential decay signal,and its normalized transfer function is

where τGAGG=100 ns is the GAGG:Ce decay time.

An accurate electrical model of SiPM is complex.However,neglecting the equivalent input resistance of FEE and quench capacitance can simplify the SiPM transient response to a single exponential decay signal,with a normalized transfer function as follows:

where τSiPM≈38 ns is the recovery time of the SenSL J-series SiPM [10,11].

With proper selection of feedback capacitance and resistance,the TIA can be treated as a second-order Butterworth filter with the following transfer function

where Ω=2πF,is the approximate–3 dB bandwidth of the TIA circuit,GBPis the gain bandwidth product of OPA656,RFis the feedback resistance,andCDis the output capacitance of the SiPM array.Here,Q=0.707 is the quality factor of the Butterworth filter.Thus,the output pulse waveform of the TIA can be expressed as

whereEis the incident photon energy,LYis the light yield of GAGG:Ce,CEis the scintillation light collection efficiency determined by the scintillator and its packaging,PDEis the PDE of the SiPM array,eis the elementary charge,andGis the gain of the SiPM.

The equivalent output noise voltage of the TIA contributed by the SiPM dark counts can be treated as a random pulse train.The standard deviation of the dark count noise can be obtained from Campbell’s theorem[12]

Considering that the output noise of the TIA is bandlimited,the equivalent output noise voltage contributed by the TIA can be estimated by a simple expression

whereINandENare the input current and voltage noise of OPA656,respectively,4kTRFis the thermal noise of the feedback resistor,andF0≈1MHz is the band-limiting frequency of the TIA and low-pass filter system.Therefore,the SNR can be defined as

The SNR varies with the feedback resistanceRF,as shown in Fig.5.RFwas set to 500 Ω to satisfy both the SNR and appropriate gain.

Fig.5 SNR variation according to RF with an incident photon energy of 59.5 keV at 20 °C,SiPM operating voltage of 28.5 V,and scintillator light yield of 32 ph/keV.The maximum SNR value is 82 at RF=330 Ω

4 Data acquisition and processing

4.1 Data acquisition electronics

The DAQ board comprises the power regulator,analogue-digital converters (ADCs),MCU,embedded multimedia card (eMMC),and communication interface to process signals from the FEE,supply electricity for all analog and digital devices,supply and control bias voltage of SiPM arrays,process commands,format data for storage,and transmit data to the spacecraft,as shown in Fig.6.The power system regulates the +5 V input voltage to ±5 V for the analog devices,+2.5 V for the ADCs,and +3.3 V for the digital devices.An adjustable bias voltage of 0–40 V for the SiPM arrays can be generated through the SiPM bias voltage supply module,which can also monitor the bias voltage and leakage current.The DAQ core is a 32-bit ARM Cortex M0+MCU KEA128 with 16 kB static random-access memory (SRAM) and 128 kB flash manufactured by NXP,which is an automotive-level MCU optimized for cost-sensitive applications and focuses on exceptional electromagnetic compatibility (EMC) and electrostatic discharge (ESD) robustness.A 512 MB single-layer cell (SLC) eMMC stores the raw science and housekeeping data with high reliability.The eMMC can store approximately 12 h of data based on existing data format definitions and background count rates of 500 counts/s.Peak hold signals from the FEE are sampled by four individual ADCs with 1M sample rate and 16-bit precision.The four internal ADCs of the MCU are alternatives of each other and can be selected through a gating switch for redundancy.

As shown in Fig.6,the DAQ electronics provides an interface with spacecrafts,which comprises a data bus using differential SPI protocol with LVDS level for the raw science and housekeeping data transmission,a UART interface for firmware update,a PPS interface for time calibration,and some general-purpose input/output(GPIO)interface for MCU reset,data request,boot configuration,and burst trigger (denoted as Telemetry and RESET).

Fig.6 Functional block diagram of DAQ and its connection with FEE,SiPM Carrier,and payload on-board computer board (POBC).The red lines represent the power bus,green lines represent digital signals,and blue lines represent analog signals.The debug interface and TTL-LVDS converter are not shown

4.2 Flight firmware

The flight firmware operates on the MCU without an operating system and comprises two parts:the boot loader and application program,residing in the internal flash memory of the MCU.The firmware provides limited online data processing capacity because of the low MCU performance.However,the configuration of all hardware,response to the triggers,command reception and processing,data storage and transmission,control and monitoring of the detector,and the application program update are provided,which can satisfy all necessary on-orbit requirements.

After the MCU is powered on,it first runs the boot loader to check the Config&Trigger pin level,and updates the application program if high,or jumps to the application program if low.

The application program is interrupt driven.Interrupts are generated on the following events:

(1) FEE trigger;

(2) Timer interruption per second to record housekeeping data;

(3) PPS from the GPS module;

(4) Spacecraft commands.

When the FEE trigger signals interrupt the MCU,the peak hold signals are sampled by four individual ADCs or four internal ADCs of the MCU,and the internal clock,PPS count,and UCT time will be recorded.When the converter sampling is finished,the MCU discharges the peak holder.Four channels work in a single thread in the present configuration,so other channels are in dead time,while one channel is triggered.The peak values and time information are stored in the eMMC.Every 1 s,the MCU records the housekeeping data,as listed in Table 2.The internal clock is recorded,while the PPS triggers the MCU for accurate time reconstruction of incident photons.

Table 2 Housekeeping data package definition

Currently,the application program processes the following commands from the spacecraft.

(1) Test communication status.

(2) Set the bias voltage of each SiPM array and trigger threshold for each channel.

(3) Process all modules self-test.

(4) Update the Coordinated Universal Time (UTC).

(5) Read data from eMMC.

(6) Read the housekeeping data.

(7) Set the ADC chosen switch.

(8) Erase the eMMC.

In daily operation,the instruction sequence will be sent to the spacecraft from the Earth station.Then,the spacecraft adjusts its attitude,powers on the detector,and controls the detector to set the SiPM bias voltage starting the observation,and powers off the detector after a specificobservation time.The data stored in the detector eMMC will be read to the POBC’s eMMC and downloaded at the appropriate time.

4.3 Data format

The DAQ produces two types of data packets:the raw science and housekeeping data.All the raw data are stored in the eMMC as a series of 512-byte packages.The definitions of housekeeping and raw science data are summarized in Tables 2 and 3,respectively.In a 512-byte raw science data package,the first trigger event occupies 3–25 bytes,and the other 43 trigger events occupy 26–499 bytes in the same format as shown,from 26 to 36 bytes.

Table 3 Raw science data package definition

In a 512-byte housekeeping data package,there are seven housekeeping datasets occupying bytes 3–492 in the same format as shown in bytes 3–72.‘‘UTC’’,‘‘PPS count’’,and ‘‘Internal clock’’ are indicated when the housekeeping data package is recorded.‘‘PPS to UTC’’indicates the corresponding PPS count last time UTC is received,and ‘‘Internal clock to PPS’’ denotes the corresponding internal clock last time PPS is received.From these data,we can precisely correspond the MCU internal clock to the UTC time,which is significant for GRB triangulation.

5 Performance

The detector performance can be investigated through experimental calibration and simulation [14].The performance of the GRID detector,including the energy-channel relations at different temperatures and biases,space uniformity,energy resolution,effective areas,angular responses,and detector noise,was calibrated experimentally in detail on the ground,and the effective areas and angular responses were investigated using a Monte Carlo simulation.The detector was irradiated with collimated radioactive sources in the laboratory (from 14 keV to 1.4 MeV) for calibration.Because of the low number of radioactive source emission lines in the low energy range,calibration measurements from 10 to 160 keV were performed with X-ray radiometry at the National Institute of Metrology of China.The detailed calibration and simulation results will be described in a future study.Here,we summarize the key features of low-energy performance,energy resolution,and high-rate performance,as well as the temperature dependence of the GRID detector.

5.1 Low-energy performance

Based on the SNR expression derived previously,if the minimum SNR required is 6,which means that the peak value of the signal amplitude is six times the standard deviation of noise,the lower limit for the low-energy detection can be derived theoretically.The light yield of the GAGG:Ce is calibrated experimentally in-house.The GAGG:Ce and a LaBr3:Ce scintillator with a known light yield are irradiated by a241Am source,respectively,and the same PMT and electronics are used for the scintillation light readout.Considering the quantum efficiencies of the PMT cathode,the light yield of GAGG:Ce is estimated to be 16 ph/keV by comparing the 59.5 keV peak positions.Under normal observation conditions,where the temperature is 20°C and SiPM operating voltage is 28.5 V with a scintillation light collection efficiency of 62% and PDE of the SiPM array,defined by the PDE of a single SiPM multiplied by the array’s fill factor,of approximately 27.9%,the lower limit of the detector is 13 keV.The light yield of the crystal is lower than that reported in a previous study [9].The provider of the GAGG:Ce crystal has improved the crystal growth process.Therefore,a crystal with a higher light yield will be used in the next detectors,with the lower limit of the GRID expected to extend to 10 keV.In addition,with the development of low dark count rates [15] and high PDE SiPMs [16],the low-energy performance of these types of detectors can be further improved.

The trigger threshold was set to 20 mV (approximately 13 keV) and a spectrum of241Am with 59.5 and 13.5 keV X-rays was measured,as shown in Fig.7.The performances of the four channels are not exactly the same owing to the difference in the scintillator light yield and light collection efficiency.However,the peak of the 59.5 keV X-ray can be observed clearly.The peak at 13.5 keV and dark count noise are mixed together near the threshold.

Fig.7 Spectrum of 241Am with a trigger threshold of 20 mV for four channels

5.2 Energy resolution

The energy resolution results were calibrated by radioactive sources in the laboratory,which handled the detector and electronics noise,as well as the statistical fluctuation and energy nonlinearity.The sources used for calibration with their emitted photon energies are listed in Table 4.As shown in Fig.8,the energy resolutionis approximately proportional toand is approximately 9%at 662 keV.The poor fitting in the low-energy region is due to the nonlinearity of GAGG:Ce,particularly the distinct inconsistency at 81 keV,which is caused by the X-ray absorption edge of GAGG:Ce at approximately 70 keV[9].

Fig.8 (Color online) Dependence of the detector energy resolution of four channels,and fitting curves and equations in terms of y=axb

Table 4 Radioactive sources used for calibration

5.3 High-rate performance

Two effects typically impair the performance of scintillation detectors at high photon rates:dead time and pulse pile-up.The pulse pile-up occurs when the count rate is so high that the pulses from successive events overlap in the FEE,which causes distortions in the measured spectrum that are difficult to characterize.These types of distortions are generally treated as systematic errors in the determination of the gamma-ray spectrum.Owing to the high bandwidth of the TIA and filter,the pulse of a single event lasts less than 1 μs from generating a trigger to recover to baseline,which incurs little distortion in the measured GRB spectrum according to the relevant discussion about the Fermi GBM detector [13].However,there is another type of pulse pile-up occurring at the peak-hold circuit;thatis,a small signal will be overridden by a larger signal.This problem can be solved by a limit switch that restricts the peak-hold circuit input voltage to the ground level when the channel is triggered.

The nominal detector dead time is approximately 50 μs per event,which mainly consists of the MCU–ADC communication and MCU–eMMC communication,as shown in Fig.9.However,in the last firmware version,which is used in the second GRID detector,the dead time is optimized to 15 μs with the same hardware design.

Fig.9 Contribution of dead time in GRID detector.When the MCU is triggered by the FEE trigger output signal,a conversion input(CNV) is generated after 6 s to initiate the ADC conversions.The ADC conversion takes 4 μs,which is followed by a 26 μs MCU–ADC SPI communication time.Then,the discharge signal is generated after 2 μs and lasts for 12 μs until the ADC value is wrapped in the eMMC

5.4 Temperature dependence

As the breakdown voltage of an SiPM varies with temperature,the gain of SiPM arrays is significantly affected by the temperature.The relationships between the gain,temperature,and bias voltages should be calibrated[17].Figure 10 shows the calibration results for the temperature dependence of the GRID detector channel 0.The gain decreases with increasing temperature,and other channels show the same changing rule as channel 0.In a GRID detector,the bias voltage does not vary with temperature to stabilize the gain.The temperature of SiPM arrays will be recorded for offline correction of the energychannel relation of the detector.

Fig.10 (Color online) Calibration results for the temperature dependence of the GRID detector channel 0 using 241Am source under bias voltages between 27.0–29.0 V

6 Conclusion

GRID mission is a student project with a dedicated and straightforward scientific goal:to detect and locate GRBs produced by neutron star mergers jointly with groundbased GW detectors in the local universe.In this paper,the detailed design of a GRID detector,electronics,and firmware,as well as the energy resolution and low-energy and high-rate performances of the detector is introduced.Further calibration of the detector,including the angular response,detection efficiency,and temperature response,using both simulations and experiments,will be reported in detail in a future study.

GRID mission was initially proposed and developed by students,with a considerable contribution from undergraduate students and shall continue to operate as a student project in future.The current GRID collaboration involves more than 20 institutions and continues to grow.The purpose of GRID mission is twofold.In addition to its scientific goals,we hope to attract excellent students from different disciplines into astrophysics and train them to organize and participate in a multi-disciplinary collaboration,while learning how to build a real science project that covers hardware,data,and science.

In conclusion,GRID mission is a scientific collaboration that accepts students and scientists worldwide.Members can launch their own detectors,share the data and produce scientific results under certain agreements.The detailed hardware and firmware design materials described in this paper will be part of a standard design package to be delivered and shared within the GRID collaboration community.

Author contributionsAll authors contributed to the study conception and design.System design,material preparation,data collection and analysis were performed by Jia-Xing Wen,Xu-Tao Zheng,Jian-Dong Yu,Yue-Peng Che,Dong-Xin Yang,Huai-Zhong Gao,Yi-Fei Jin,Xiang-Yun Long,Ji-Rong Cang and Ming Zeng.The first draft of the manuscript was written by Jia-Xing Wen,and all authors commented on previous versions of the manuscript.All authors read and approved the final manuscript.