Development of a wide-range and fast-response digitizing pulse signal acquisition and processing system for neutron flux monitoring on EAST

The neutron count rate fluctuation reaches six orders of magnitude between the ohmic plasma scenario and high power of auxiliary heating on an experimental advanced superconducting tokamak (EAST). The measurement result of neutron flux monitoring (NFM) is a significant feedback parameter related to the acquisition of radiation protection-related information and rapid fluctuations in neutron emission induced by plasma magnetohydrodynamic activity. Therefore, a wide range and high time resolution are required for the NFM system on EAST. To satisfy these requirements, a digital pulse signal acquisition and processing system with a wide dynamic range and fast response time was developed. The present study was conducted using a field-programmable gate array (FPGA) and peripheral component interconnect extension for instrument express (PXIe) platform. The digital dual measurement modes, which are composed of the pulse-counting mode and AC coupled square integral’s Campbelling mode, were designed to expand the measurement range of the signal acquisition and processing system. The time resolution of the signal acquisition and processing system was improved from 10 to 1 ms owing to utilizing high-speed analog-to-digital converters (ADCs), a high-speed PXIe communication with a direct memory access (DMA) mode, and online data preprocessing technology of FPGA. The signal acquisition and processing system was tested experimentally in the EAST radiation field. The test results showed that the time resolution of NFM was improved to 1 ms, and the dynamic range of the neutron counts rate was expanded to more than 106\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{6}$$\end{document} counts per second. The Campbelling mode was calibrated using a multipoint average linear fitting method; subsequently, the fitting coefficient reached 0.9911. Therefore, the newly developed pulse signal acquisition and processing system ensures that the NFM system meets the requirements of high-parameter experiments conducted on EAST more effectively.


Introduction
The neutron yield reflects the internal nuclear reaction rate of the experimental fusion device. The neutron flux monitoring (NFM) measurement result is a significant feedback parameter related to the acquisition of radiation protection-related information and rapid fluctuations of neutron emission induced by plasma magnetohydrodynamic activity [1][2][3]. Therefore, NFM is crucial in the physical research of high-temperature plasma and safety control for the experimental advanced superconducting tokamak (EAST). When high-parameter and long-pulse experiments are carried out on EAST, neutral beam injection (NBI) auxiliary heating systems are added [3][4][5].
During these experiments, fusion neutrons are not only obtained from the reaction of main plasmas but also from the reaction of injected neutral beam ions with main plasmas and neutral beam ions. These reactions are the main source of fusion neutrons [5][6][7][8]. During the entire experiment, the neutron counts rate can exceed 10 6 counts per second [3,5]; therefore, a wide measurement range is required by the NFM system utilized on EAST. However, traditional reactors generally use different types of detectors to extend their range through a range relay [9,10]. Because different types of detectors have different reaction cross sections for neutrons of different energies, the detection efficiencies for neutrons of the same energy are different as well, which are not conducive to the conversion and continuity of the measurement data.
In addition, an NFM system with a high time resolution is required for acquiring information on rapid fluctuations [5,[11][12][13]. However, in the existing NFM systems [5], the analog-to-digital converter (ADC) sampling rate of the signal acquisition and processing system is 250 million samples per second (MSPS), core processing module adopts the EP3C16Q240C8N (Cyclone III) field-programmable gate array (FPGA) chip produced by Altera, and communication is based on the PXI bus with programmed input-output (PIO) mode. These factors result in a slow data acquisition and communication speed, which limits the time resolution. The time resolution of the existing NFM system is only 10 ms [5], which cannot adequately observe the rapid fluctuations of neutron emission induced by plasma magnetohydrodynamic activity [2]. Moreover, the time resolutions of NFM systems based on the 235 U fission ionization chamber on JT-60U and DIII-D are 5 ms and 10 ms, respectively [13,14]. The time response of an NFM on a large helical device (LHD) is upgraded to 0.5 ms [15,16]. The NFM on joint European torus (JET) and international thermonuclear experimental reactor (ITER) are preparing for upgrading to meet the requirement of 1 ms time resolution [2,17].
Therefore, the NFM on EAST must be upgraded. In this study, a digital pulse signal acquisition and processing system with a wide dynamic range and fast time response was developed based on FPGA. In the proposed system, digital dual measurement modes that are composed of the pulse-counting mode and AC coupled square integral's Campbelling mode were designed. The Campbelling mode was calibrated with a multipoint average linear fitting method using the overlap region data of the dual modes. This significantly expanded the measurement range of the NFM system and maintained the continuity of the measurement data. Moreover, this approach adopted highspeed ADCs to improve the sampling rate from 250 to 500 MSPS. The core processing module adopted a high-performance xc7k325tffg900 (Kintex-7) FPGA chip manufactured by Xilinx. The old communication based on a PXI bus with a PIO mode was replaced by a PXIe bus communication with a direct memory access (DMA) mode. Accordingly, the time resolution of the NFM system was improved from 10 to 1 ms. The pulse signal acquisition and processing system was tested experimentally in the EAST radiation field. Moreover, the test results were analyzed and discussed. The measurement results of the dual modes were compared with the variation trend of the injection power of neutral beam injection (P NBI ), and the effective measurement range of the dual modes was discussed based on the experimental data for different injection powers of NBI.

Design and algorithm
The NFM system used on EAST consists of 235 U fission ionization chambers (FCs), charge-sensitive preamplifiers, signal acquisition and processing system, and upper computer software. An FC765 fission ionization chamber manufactured by Centronic company was utilized, which has a neutron sensitivity of 1.0 cps/nv and an output pulse width of 100 ns. The low-noise preamplifier (CAP60) matched to the FC produced by Cook electronics company has a rapid rise time (10 ns), typical output pulse width of approximately 100 ns, and output impedance of 50 X. When the input current pulse is 1 lA, the output-voltage amplitude of the preamplifier is 50 mV. In addition, because the FC output signal amplitude is small and needs to be transmitted over a long distance, the preamplifier in the radiation-shielding box was positioned close to the FC to effectively amplify the output signal. Both the FC and preamplifier were installed on top of the EAST host. To avoid radiation damage [18], a signal acquisition and processing system was installed in a special neutron diagnostic laboratory with a radiation-shielding wall.
The workflow of the NFM system is as follows: First, the neutron signals are properly amplified by a preamplifier and then transmitted to the signal acquisition and processing system for processing. Later, the data processed by the signal acquisition and processing system are transmitted to the upper computer for storage through the PXIe bus on the back panel of the chassis [19,20]. Thus, the signal acquisition and processing system is an important carrier for the NFM system function realization.

Hardware design of a signal acquisition and processing system
A signal acquisition and processing system was designed based on the PXIe bus architecture, as shown in Fig. 1a. It adopts the structure of mother-daughter board that are connected through printed circuit board (PCB) connectors manufactured by TE Connectivity. The board was designed with four channels, only one of which channel was used for testing during the experiment.
The daughterboard shown in Fig. 1a, mainly includes an analog signal processing module and waveform digitization module. The analog signal processing module is composed of a SMA interface, AC coupling, impedance matching, overvoltage protection, amplifier, single-ended to differential, and anti-aliasing filtering. The output signals of the preamplifier are transmitted to the daughterboard through the SMA interface. The signals are then sequentially processed by AC coupling, impedance matching, and overvoltage protection. AC coupling can isolate the DC component of the signal, help compress the amplitudes of the piled-up signals, and suppress the drift of the baseline. Impedance matching uses a 50X resistor that matches the characteristic impedance of the transmission line in parallel with the terminal load. Thus, the reflection coefficient of the terminal is 0, and the pulse waveform distortion caused by signal reflections at both ends of the transmission line is avoided. To avoid damaging the device owing to excessive input voltage, the overvoltage protection circuit uses a pair of clamping diodes to limit the signal amplitude to AE5 V: A amplifier (LMH6703) [21] is used to adjust the signal amplitude within -1-? 1 V to maximize the use of the ADC's dynamic range and improve the signal-to-noise ratio of the ADC. The singleended differential circuit converts a single-ended analog signal coming from the amplifier (LMH6703) [21] into a differential signal, which enters the waveform digitization module after anti-aliasing filtering and is converted into a digital signal through sampling, quantization, and coding.
The waveform digitization module includes high-speed ADCs and a clock jitter cleaner. Owing to the narrow pulse width and fast rising edge of the neutron signal described above, as well as the Campbelling integration accuracy requirement, a single-channel pipeline ADC (ISLA214P50) [22] with a sampling rate of up to 500 MPS, resolution of 14 bits, and high conversion rate is selected. The ADC can conduct sampling at a sufficient rate in each pulse signal for algorithm analysis. Furthermore, the ADC provides a 250 MHz associated clock. The input double data rate (IDDR) primitive of the FPGA is used for double-edge (rising and falling edges) sampling. This method reduces the FPGA's processing time; thus, more data can be processed and a high time resolution can be achieved. To achieve the desired effective number of bits, a clock jitter eliminator is used to reduce the jitter of the ADC sampling clock. Here, the chip (LMK04828B) [23] of TI is selected, which can generate an ultralow jitter clock to meet the requirement of the ADC sampling clock.
The motherboard, which is shown in Fig. 1b, mainly includes a core processing module, external trigger module, and bus interface. The core processing module mainly processes the ADC output digital signals. Considering factors such as speed grades, pins, resources, number of high-speed transceiver channels, and price, the xc7k325tffg900 [24] FPGA chip (Kintex-7) manufactured by Xilinx was selected. Important functions such as the pulse count algorithm of the pulse-counting mode, square integral algorithm of the Campbelling mode, and data buffering were all realized in the FPGA. The external trigger module receives and processes the external trigger signal taken from the EAST central control system and initiates the signal acquisition and processing system. The bus interface adopts a PXIe framework. It is designed based on the PCIe IP core of the Xilinx 7 series FPGA integrated block. Considering the 8 b/10 b encoding overhead, the theoretical unidirectional bandwidth can reach 32 Gb/s, which is sufficient to meet the current application of the NFM system and requirements of future upgrades. The GTX high-speed transceiver, PCIe hard core, and related IP cores are integrated into the FPGA. It greatly facilitates the hardware and FPGA programming and improves the integration of the board. In this study, the communication between the signal acquisition and processing system and upper computer is realized via a highspeed PXIe protocol in DMA mode. The DMA mode of the PXIe communication adopts a high-speed data packet transmission and is controlled through interruption.

FPGA design
A logical design block diagram of the FPGA used in the signal acquisition and processing system is presented in Fig. 2. It mainly includes the trigger control logic module, clock management module, jitter eliminator configuration module, ADC configuration and data preprocessing module, pulse-counting mode, Campbelling mode, first input first output (FIFO) module, and PXIe communication module with a DMA mode.
The trigger control logic module realizes the hardware trigger function of the NFM system. If the signal acquisition and processing system receives the hardware trigger signal from the central control system, the NFM system starts to operate. The clock management module realizes a low-skew clock distribution based on high-speed buffers and routing, and supports frequency synthesis and phase shift. It uses the 100 MHZ differential clock provided by the chassis backplane as the backplane triggering clock signal, which can be applied to the synchronous triggering mechanism of multiboard cards in the future. The clock jitter eliminator configuration module mainly utilizes the configuration of the clock jitter eliminator constructed via a serial peripheral interface (SPI) communication, such that the clock jitter eliminator outputs low-jitter differential 500 MHZ clocks as the sampling clock frequency of the ADC. The ADC configuration and data preprocessing module mainly complete the SPI communication and register mechanisms for the ADC. It also converts high-speed raw data into preprocessed data that can be recognized by the next module according to the defined rules. Then, the preprocessed data are divided into two parts. One part enters the pulse-counting mode to directly count the neutron pulses, while the other enters the Campbelling mode to obtain a Campbelling square integral value (hereafter called the Campbelling value). The two modes work at the same time. Their measurement results are acquired every 1 ms and buffered in the FIFO module. Finally, the buffered data are transmitted at high speed to the upper computer for storage by the PXIe communication module with a DMA mode.

Count algorithm
In case of low neutron counts rate, the pulse-counting mode can be adopted by counting the pulses in the timesharing scale, as the FC outputs approximately discrete pulse signals. However, there are many gamma rays and X-rays present in the environment where the FC of the NFM system is located. Furthermore, alpha particles emitted from spontaneously decaying 235 U can produce many small irregular pulses in the FC. These will result in deviations in the counting. To decrease these deviations, it is necessary to eliminate the interference of non-neutron signals. Because most neutron and non-neutron particles are in two non-overlapping amplitude ranges, a suitable pulse amplitude discrimination threshold (high  Fig. 2 Logical design block diagram of FPGA threshold in this study) is set to decrease the influence of interfering particles. Furthermore, various electronic noises with small amplitudes and wide spectrum ranges are superimposed on the pulse waveform, which can easily cause false counting measurements near the discrimination threshold. A low threshold can be set to reduce the influence of electronic noise. Therefore, a dual-threshold discrimination method is adopted to decrease counting error.
To better distinguish the two piled-up neutron pulses, the low threshold is generally set slightly lower than the high threshold. Only when an input pulse signal completely passes the high and low thresholds, the neutron pulse count is increased by one. Thus, amplitude discrimination can eliminate most of the interference and noise; subsequently, the pulse-counting mode can achieve good measurement results at low counts rate.

Campbell algorithm
However, in case of high neutron counts rate, many pulse signals accumulate, which can result in a large deviation in the pulse-counting mode. Fortunately, the Campbelling mode is suitable for high neutron counts rate. An FPGA is adopted to calculate the Campbelling value of the sampling data obtained from the ADC. The calculation is based on the Campbell theorem [10,[25][26][27][28].
The output signal voltage of the FC can be approximated using the impulse current Qd t ð Þ, where Q represents the total charge of each current pulse output by the actual detector and d t ð Þ is the unit step function. Because the number of pulse signals generated by the FC in unit time obeys the Poisson distribution principle, the average contribution of the pulse impact within s $ s þ ds of the output voltage at time t is defined as follows: where k represents the average count rate of the pulse signal and h(t) represents the impulse response of the electronic system. The root-mean-square deviation of dv caused by the randomness of the pulse number k is stated as: Because small ds intervals are independent of each other, the mean squared error produced by the count fluctuation is given as: When t \ 0 and h(t) = 0, Eqs. (1) and (3) integrate s to obtain Eqs. (4) and (5), respectively: The charge quantity Q of each current pulse output by the actual detector suffers from a fluctuation. The fluctuation is related to the energy spectrum of the incident particles and background and statistical fluctuation produced by the detector. Generally, the average value of charge Q is used instead of Q. The integral expressions of the Campbell theorem are as follows: v ¼ kQ According to Eqs. (6) and (7), the Campbelling mode utilizes two calculation methods. One method is the direct integration of the signal, which is called the current mode. The other method, which reflects neutron flux information by measuring the mean squares error of the voltage signal, is called the mean square voltage (MSV) mode. The Campbelling values of the two methods are proportional to the average count rate k of the pulse signal. In other words, their measured values are proportional to the neutron flux density. Because the first method is more affected by baseline drift, and the second method has a stronger suppression effect on the non-neutron background with a lower amplitude, this study adopts the second calculation method.
In addition, the integral of MSV can be as depicted in Eq. (8): In this study, AC coupling is used to isolate the DC offset of the signal. Thus, v ¼ 0, and Eq. (8) can be converted into Eq. (9) as follows: Because v 2 is proportional to the count rate k, the value of MSV can be used to replace the count rate k. In addition, the DC offset of the signal generated by the piled-up pulse, background, and baseline of the preamplifier is eliminated in this method. Therefore, the method has a high n=c rejection ratio, which can help compress the signal amplitude and improve the dynamic range upper limit of the Campbelling mode.

Experimental results and discussion
To test the performance of the signal acquisition and processing system, an experiment was conducted on the EAST. The data from shot No.102884 were selected to qualitatively analyze the dual modes of the signal acquisition and processing system, which is depicted in Fig. 3. In the shot, the basic parameters are plasma current (I p * 0.4 MA), central electron density (n e * 4:8 Â 10 19 m À3 ), injection power of the low hybrid wave (P LHW * 2 MW), and injection power of NBI (P NBI * 1.6 MW) at the flattop. The dual modes obtained measurement results every 1 ms, indicating that the time resolution of the NFM system reached 1 ms.
In addition, when fusion neutrons are generated because of the added NBI auxiliary heating system, the change trend in the dual-mode test results should be consistent with the change trend of P NBI . Figure 3 demonstrates that NBI starts to be injected at 2.021 s, and the P NBI continues to remain constant in the subsequent 8.019 s. The neutron counts rate and Campbelling values measured in the dual modes start to increase at 2.021 s and both last approximately 8.019 s. Finally, at 10.04 s, when NBI injection stops, the dual-mode measurement results change accordingly. The duration of the neutron pulse signal measured by the NFM system is consistent with that of the NBI. Therefore, we can conclude that the change trends of the neutron counts rate and Campbelling values in the dual modes are consistent with the change trend of P NBI . Based on the abovementioned experimental results, we confirmed that the dual modes of the signal acquisition and processing system are reliable and suitable for neutron measurement and can meet the requirement of attaining a 1 ms time resolution.

Multipoint average linear calibration for Campbelling mode
When the neutron counts rate is low, the background noise occupies a larger proportion of the Campbelling values, which results in larger deviations in the Campbelling values. The counts rate measured by the pulsecounting mode can better reflect the neutron counts rate. When the neutron counts rate is in the middle region [ (approximately 40-500 counts per millisecond (cpms) [5,25]], the results measured by the pulse-counting and Campbelling modes can accurately reflect the neutron counts rate. Thus, this region is called the overlapping region in the dual mode [29,30]. When the neutron counts rate is high, the pulse signals will accumulate significantly because of the pulse signals with a certain pulse width. The neutron counts rate measured in the pulse-counting mode is lower than realistic counts rate. Fortunately, the Campbelling mode directly calculates the squared integral of the voltage signal, thus, its results will not be affected by the piled-up pulse and can accurately reflect the neutron counts rate. However, the Campbelling value is proportional to the neutron count rate value. Therefore, it is necessary to convert the Campbelling value into the neutron count rate via proper calibration.
As mentioned previously, there is an overlapping region in the dual modes. Thus, the Campbelling mode can be calibrated using the neutron count rate of the pulsecounting mode. The relationship between these two modes in the overlapping region is defined as follows: Here, k is mainly determined by the neutron pulse waveform, amplitude, and number of signal sampling points, and b reflects the contribution of background noise, which mainly depends on the magnitude of the electronic noise and various non-neutron background rays. Count and Campbell represent the neutron count rate of pulse-counting mode and Campbelling value of Campbelling mode, respectively. To accurately obtain the parameters k and b, this study adopted a multipoint average linear fitting method. First, all the data with neutron counts rate ranging from 40 to 500 measured in the pulse-counting mode were extracted. Then, to suppress the random error in the Campbelling mode, the Campbelling values under the same count rate are extracted and averaged. Finally, the Campbelling average values and corresponding counts rate are linearly fitted using Eq. (10). In this study, shot 104218 (where P NBI is approximately 870 kW) was selected to calibrate the Campbelling mode. The results are presented in Fig. 4 and Eq.  Fig. 5. In the overlap region, when neutron counts rate fluctuates, the change trends are consistent between the counts rate converted by the Campbelling values and those directly measured in pulsecounting mode. These values can overlap significantly. These results indicate that the results calibrated by Campbelling mode can be used to calibrate different shot data with different neutron count rate fluctuations while attaining significant reliability.

Comprehensive test
In this study, to obtain the effective measurement range of the dual modes, the experimental data of different P NBI were analyzed. The dual-mode measurement results of shot No.101955 with a high level of injection power are presented in Fig. 6a. The basic parameters of the shot are plasma current (I p * 0.4 MA), central electron density (n e * 7:2 Â 10 19 m À3 ), and injection power of NBI (P NBI * 3.5 MW) at the flattop. The red line represents the neutron counts rate measured in the pulse-counting mode, while the blue line represents the neutron counts rate converted by the Campbelling values through calibration. The neutron counts rate measured in the Campbelling mode was significantly higher than those in the pulse-counting mode. The reason for this is as follows: the neutron pulse signals are piled up in the case of high neutron counts rate, which results in smaller measured values in the pulse-counting mode; however, the piled-up pulse signals have almost no effect on the Campbelling mode. In addition, because the pulse-counting and Campbelling modes perform well in the case of low and high neutron counts rate, respectively, they can be used as realistic counts rate for evaluating the deviations of another mode. In this study, the shot data are further analyzed. When cpms is lower than 500 (including the overlap region), the counts rate measured in the pulse-counting mode can be used as realistic counts rate and those measured by the Campbelling mode can be used as measured counts rate. When cpms is higher than 500, the opposite is true, which is demonstrated in Fig. 6b. In the overlap region, the counts rate converted by the Campbelling values through calibration is approximately equal to the realistic counts rate, which again indicates that the dual modes attain good working performance, and the calibration method and results of the Campbelling mode are sufficiently reliable. In the high-count-rate region, the deviation increases with increasing count rate and is approximately 19.8% as the count per millisecond reaches 1500. These results prove that the pulse-counting mode is unsuitable for high neutron counts rate.
The data for shot No. 104216 are selected to analyze the measurement effect of the Campbelling mode with low counts rate, as shown in Fig. 6c, d. The basic parameters of the shot are the plasma current (I p * 0.4 MA), loop voltage (V loop * 0.1 V), and injection power of NBI (P NBI * 0.9 MW) at the flattop. When the neutron counts rate is lower than those measured in the overlap region, the counts rate in the Campbelling mode is higher than the realistic counts rate (counts rate measured in pulse-counting mode). This is because the counts rate measured in the Comprehensive tests with different injection powers of NBI. a The comparison between the converted counts rate in Campbelling mode and neutron counts rate directly measured in pulse-counting mode inside the high neutron count rate region; b The converted counts rate in Campbelling mode is used as the realistic counts rate for evaluating the deviations of the neutron counts rate directly measured in pulse-counting mode inside the high neutron count rate region; c The comparison between the converted counts rate in Campbelling mode and neutron counts rate directly measured in pulse-counting mode inside the low neutron count rate region; d The neutron counts rate directly measured in pulsecounting mode is used as the realistic counts rate for evaluating the deviations of the converted count counts in Campbelling mode inside the low neutron count rate region Campbelling mode is interfered by non-neutron background such as noise, which will increase the neutron counts rate measured in the Campbelling mode and lead to inaccuracies in the measurement results.
In summary, when the neutron counts rate is lower than those measured in the overlap region, the results measured via the pulse-counting mode can be selected as realistic neutron counts rate. In the overlap region, the results measured by pulse-counting mode or Campbelling mode can be selected. However, it is best to select the results measured by the pulse-counting mode in the overlap region with low neutron counts rate, and select the results measured by the Campbelling mode in the overlap region with high counts rate to avoid greater deviations. When the neutron counts rate is higher than those measured in the overlap region, the results measured by the Campbelling mode can be utilized. Therefore, combining the dual modes expands the measurement range of the NFM system to 10 6 cps.

Conclusion
In this study, a digital pulse signal acquisition and processing system with a wide dynamic range and fast response time was developed for EAST. Dual modes were used in the signal acquisition and processing system to expand the measurement dynamic range of the neutron counts rate to more than 10 6 cps. The time resolution of the signal acquisition and processing system was improved from 10 to 1 ms. Moreover, a multipoint linear calibration for the Campbelling mode was conducted. The fitting coefficient reached 0.9911, and the calibration results of the Campbelling mode demonstrated good reliability. Thus, the newly developed pulse signal acquisition and processing system can ensure that the NFM system meets the requirements of high-parameter experiments. The fluctuation of neutron counts rate can reach six orders of magnitude in a field experiment involving EAST. However, in the future, the fluctuation is expected to be as high as eight orders of magnitude in the D-T operation phase. The system should be further verified for a higher range expansion of the NFM system with higher neutron counts rate.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Author contributions All authors contributed to the study conception and design. Material preparation, data collection, and analysis were performed by Li Yang, Hong-Rui Cao, Jin-Long Zhao, Zi-Han Zhang, Qiang Li, and Guo-Bin Wu. The first draft of the manuscript was written by Li Yang and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.