Cenozoic exhumation in the northeastern Tibetan Plateau provides insights into spatial-temporal patterns of crustal shortening, erosion, landscape evolution, and geodynamic drivers in the broad India-Eurasia collision system. The NW-SE trending West Qinling Belt has been a central debate as to when crustal shortening took place. Within the West Qinling Belt, a thick succession of Cretaceous sedimentary rocks has been deformed and exhumed along major basin-bounding thrust faults. We present new apatite (U-Th)/He ages from the hanging wall and footwall of this major thrust. Contrasting thermal histories show that rapid cooling commenced as early as ca. 45 Ma and continued for 15–20 Myr for the hanging wall, whereas the footwall experiences continuous cooling and slow exhumation since the late Mesozoic. We infer that accelerated exhumation was driven by thrusting in response to the northward growth of the Tibetan Plateau during the Eocene (ca. 45–35 Ma) based on regional sedimentological, structural, and thermochronological data.
Ongoing collision between India and Eurasia has exerted a profound influence on landscape development, fault activation, and seismicity in East Asia [1–3]. Continuing continental convergence may have reactivated pre-Cenozoic orogenic belts (e.g., Tian Shan-Altaids, Qilian Shan, and Qinling Belt) where preexisting faults of crustal scales exist [2, 4–8]. Studies have suggested that basin formation, mountain uplift, fault activity, and plateau growth may have protracted throughout the Cenozoic Indo-Asian collision, but details regarding crust shortening histories are still enigmatic [2, 4–10].
There are three schools of thoughts regarding the timing of deformation along the West Qinling Belt (WQB) in the northeastern Tibetan Plateau. First, oblique stepwise subduction of the Asian lithospheric mantle may have resulted in the growth and progressive development of the Tibetan Plateau toward the northeast from active plate boundaries . The onset of flexural basins (e.g., Linxia Basin) and provenance change in the Lanzhou Basin imply pulsed crust shortening and the growth of the WQB during the late Oligocene [11, 12]. Second, simultaneous contractional deformation in the northeastern plateau margin may have taken place in the early Cenozoic, accompanying the initial development of the Tibetan Plateau, implying a relatively constant bulk strain rate of the Asian lithosphere around the same time [8, 13]. Low-temperature thermochronologic data from hanging wall blocks of major thrust systems [14, 15], depositional history of major sedimentary basins [16–18], and vertical axis rotation of subblock/basins  indicate that the onset of mountain building took place along the northeastern plateau margin associated with the initial growth of the Tibetan Plateau. Third, the presence of unconformities between lower Cretaceous and Tertiary units along the WQB could be related to mountain building and significant uplift during the late Cretaceous, which may have no connection with the Cenozoic growth of the Tibetan Plateau [20, 21].
Whether significant contractional deformation took place in the WQB between late Cretaceous and Paleogene and whether such events were widespread of WQB remain critically unknown. When the Tibetan Plateau grew northward to the northeastern region remains enigmatic. For clear answers to this question, we used low-temperature thermochronology to examine the exhumation history and possible topographic evolution of the WQB by constraining the cooling of the hanging wall and footwall of a major thrust fault. Our new data, in combination with previously published geological data, are analyzed to advance the understanding of orogenic processes in the northeastern Tibetan Plateau.
2. Geological Setting
The Qinling Orogen, which trends E-W in central China, was formed during the middle-late Triassic collision between the North China Block and the South China Block [20, 22, 23]. The Qinling Belt is divided into the West Qinling Belt (WQB) and the East Qinling Belt, separated by the Foping Dome (Figure 1(a)) . The WQB may have experienced at least two major episodes of shortening during Mesozoic time [4, 20, 24–27]. The early episode has been attributed to the collision between the North China Block and the South China Block in the Mesozoic [22, 23, 28, 29], and the later stage is thought to be related to the northward growth of the Tibetan Plateau during the Cenozoic [2, 4, 12, 14, 16, 30]. Deep seismic reflection shows that the active deformation is dominated by upper crustal shortening, middle crustal shearing, and mantle-derived magmatism since the Oligocene (Figure 1(b)) [24, 31, 32].
The WQB is bounded by the West Qinling Fault (WQLF) and the Shangdan Suture Zone (SDSZ) to the north, the Mianlue Suture Zone (MLSZ) to the south, the Foping Dome to the east, and the Qilian Shan to the west-northwest (Figure 1(a)). The WQLF and SDSZ, which extend NW-SE for ~800 km, are characterized by early Paleozoic ophiolitic rocks that mark the Shangdan Ocean . Since the Miocene, the WQLF has been active as a left-lateral strike-slip fault. With an estimated slip rate during the Quaternary of 2–3 mm/yr , the WQLF is a major fault that connects the seismically active Qilian Shan in the west to the tectonically active Weihe Basin, which is located next to the relatively stable Longzhong Basin (Figure 1(b)). The north-dipping MLSZ lies north of an arcuate thrust that was active during the middle-late Triassic collision between the North China Block and the South China Block . Along the WQB, the elevation gradually decreases toward the east.
Four arcuate south-vergent thrust faults are found in the WQB (Figure 1(a)) . The Lintan-Tanchang Fault (LTF), which cuts Devonian sandstones and Carboniferous shales, is intruded by Triassic granitoids. Intensely folded lower to middle Triassic clastic rocks are found along the LTF and Guanggai Shan-Die Shan Fault (GDF). Rocks along the Diebu-Wudu Fault (DWF) include Silurian sandstones and siltstones, as well as minor Devonian and Carboniferous sandstones. The Kangxian-Wenxian-Maqu Fault (KMF) is a secondary fault of the SDSZ and reactivated by Mesozoic and Cenozoic intracontinental orogenesis.
In the WQB, a regional angular unconformity separates lower Cretaceous sediments of the Donghe Group from underlying pre-Jurassic strata (Figures 1(a) and 2(a)). Upper Cretaceous strata are absent. Lower Cretaceous rocks consist primarily of purple-red alluvial-fluvial-lacustrine deposits, dominated by conglomerate, sandstone, and mudstone, with intercalations of coal measures in the lower part. Paleogene strata (Guyuan Group), which are separated from the pre-Cenozoic rocks by an angular unconformity, crop out sparsely within intermontane basins in the WQB [21, 34]. The Guyuan Group is dominated by terrestrial red beds and massive matrix-supported conglomerate. Oligocene vertebrate fossils (the Longjiagou fauna) have been reported from the Guyuan Group [21, 34–36].
The Huicheng Basin (HCB in Figure 1) is a northeast-trending intermountain basin in the southeastern WQB, bounded by the KMF to the south and the Chenxian Fault (CXF) to the north. The Huicheng Basin consists of lower Cretaceous rocks (Donghe Group), whose ages are determined by spore-pollen assemblages [37, 38]. Rocks are dominated by purplish-red conglomerate, gray pebbly sandstone, grayish-green mudstone layers with shale interbeds, and coal seams. The succession is divided into the Tianjiaba (K1t), Zhoujiawan (K1z), and Huaya (K1h) Formations, whose ages were determined by fossil assemblages, lithofacies associations, and depositional contacts. The lowermost K1t consists of purplish-red massive conglomerate facies and well-sorted, coarse- to fine-grained glutenites, which have been interpreted to represent a transition from alluvial fan to sandy braided river lithofacies. The overlying K1z includes delta fan lithofacies, consisting of coarse sandstone interbedded with thin siltstone and mudstone. The uppermost K1h is dominated by a sequence of brown-red mudstone-sandstone and thickly bedded massive red conglomerate. This formation has been interpreted to represent a shallow lacustrine, sandy braided river and alluvial fan environment. The minimum total thickness of the Donghe Group is ca. 7000 m .
Early Cretaceous rocks in the Huicheng Basin were folded and thrust faulted (Figures 2(b) and 2(c)). The entire basin is an open synclinorium in the hanging wall of the Mianlue Suture Zone (Figure 3). Limbs of the synclinorium dip 40°–60° (Figures 2(c) and 3). The restored thickness of Cretaceous strata and the uncertainties concerning fault dips in depth conspire to give the minimum and maximum shortening estimates of ca. 10–20 km (Craddock et al.  and this study).
3. (U-Th)/He Thermochronology
3.1. Samples and Methods
Low-temperature thermochronological data derived from vertical elevation/depth profiles could provide constraints on the thermal and exhumation history [10, 14, 41, 42]. Apatite (U-Th)/He (AHe) low-temperature thermochronology has a partial retention zone (PRZ) from 55 to 80°C . Assuming an average regional geothermal gradient of 25°C/km , depths of the PRZ range from ca. 2 to 3 km. AHe data provide information on the cooling path through the PRZ, and rapid cooling might be interpreted as a result of uplift and following unroofing of rocks in response to tectonism. We collected seven samples from medium-grained lower Cretaceous sandstone (KB profile) across the Huicheng Basin for apatite (U-Th)/He analysis (Figure 3). These samples were collected along a roadcut that run perpendicular to the basin along a distance of >20 km, exposing a ~7 km section. The sampling strategy was to collect a full profile covering the largest paleodepth range of the hanging wall of the Huaqiao-Ganquan Fault (HGF) in order to use AHe ages to constrain the exhumation of the hanging wall of this thrust. We use the total thickness of Cretaceous strata to estimate the paleodepth of each sample. Additionally, five samples from a single Triassic pluton (TP profile) in the footwall of the HGF were dated with AHe. Five samples were collected from a ca. 1 km vertical transect (Figure 3). Based on our field observations, the two sampling profiles (TP and KB profiles) are not cut by brittle faults (Figure 3).
We obtained 60 single-grain ages from 12 samples. Typically, five single-grain analyses were used to calculate mean ages (Table DR1). Age uncertainties are reported as the standard error of replicate analyses for individual samples using an average standard deviation for each data subgroup . 14 samples with replicate ages that exceeded 30% standard deviation were rejected because we considered the mean age of these samples to be less well constrained than the rest of the dataset which had much higher reproducibility (Table DR1). Factors that cause age dispersion in the AHe system may include radiation damage, crystal size, and U and Th zoning [43, 47–49]. These effects become more pronounced for samples that cool slowly through the apatite partial retention zone [43, 50]. AHe age-eU or 4He plots show positive correlations for TP3, TP4, and KB5 (Figs DR1 and DR2), implying a role of radiation damage on helium diffusion in apatite . Most AHe age-eU or 4He plots do not show any correlation (Figs DR1 and DR2).
Apatite (U-Th)/He (AHe) ages from the TP profile range between 100 Ma and 45 Ma, which are younger than pluton crystallization ages (ca. 210 Ma) indicating complete resetting of the AHe system [20, 42]. AHe ages on the vertical profile increase with elevation with a distinct change in the apparent exhumation rate at ca. 1400 m (Figure 4(a) and Table DR1). Below 1400 m elevation, analyses define a relatively steep age/elevation gradient with an increase from 62 to 45 Ma (Figure 4(a) and Table DR1). Above 1400 m, analyses define a gentle age/elevation gradient with ages increasing from 66 to 100 Ma (Figure 4(a) and Table DR1). In the TP profile, age-elevation relationship defines a gentle gradient from ca. 100 to 50 Ma, which is consistent with slow exhumation.
Given the thickness (ca. 7 km) of the Donghe Group and the presence of coal seams in the lower part of K1t, we suggest that the KB profile could have been buried to completely anneal the apatite (U-Th)/He system and totally reset the AHe age (e.g., ). Seven samples from the KB profile with a narrow range of late Eocene age range from 40 to 25 Ma (Figure 4(b) and Table DR1), much younger than its stratigraphic ages. In the KB profile, AHe ages decrease monotonically with depth, ranging between ca. 40 Ma and 25 Ma (Figure 4(b) and Table DR1). The slope of the age-elevation relationship suggests an apparent exhumation rate 0.225 km/Myr during the Eocene.
The Bayesian transdimensional Markov Chain Monte Carlo method was applied to the TP and KB profiles, using the QTQt program (v. 5.4.0) . Data that were deemed inappropriate for inverse thermal modeling were excluded [47, 48], and helium diffusion in apatite is modeled using the diffusion parameters of Farley . Models represent 300,000 iterations, with the first 150,000 iterations used to stabilize or burn in the inversion process and the rest of the iterations used to form the posterior ensemble . Exploratory runs using larger numbers of iterations did not alter the inverse thermal model results. Geological constraints used for modeling include (1) a present-day temperature of °C; (2) a present geothermal gradient of °C/km ; (3) the initial temperature °C at Ma for the TP profile, according to published apatite fission track data [4, 25]; and (4) the initial temperature constraint at °C at Ma for the KB profile according to the regional geothermal gradient, maximum sedimentary thickness, and maximum burial age of the Huicheng Basin [37, 44, 52, 53].
Thermal history modeling for five samples from the TP profile indicates a protracted late Cretaceous cooling at a rate of ca. 0.5–1°C/Ma, followed by the early Cenozoic accelerated cooling at a rate of ca. 2°C/Ma (Figure 4(c)). Modeling for the KB profile results in a history of temperature offset between the uppermost and lowermost samples. The cooling history of the KB profile reveals episodic phases of cooling since late Cretaceous time, including a relatively prolonged slow cooling during late Cretaceous and Paleocene time (ca. 2–3°C/Ma), a distinct phase of relatively rapid cooling at ca. 9°C/Ma during the middle-late Eocene, and followed by a Neogene-Quaternary almost thermal quiescence with minor cooling (ca. 0.5–1°C/Ma). Modeled ages fit reasonably well with observed data (Figures 4(a) and 4(b)).
In summary, modeling that results from two profiles show episodic cooling. First, pre-Paleocene cooling appears to be slow, which is consistent with previous works in the WQB [4, 14, 25, 54, 55]. Second, the initiation of the accelerated exhumation is refed to as early as Eocene time (ca. 50 Ma), with a possible rapid cooling phase at middle-late Eocene time (ca. 45–35 Ma). Last, exhumation since the Oligocene seems to be less constrained. This cooling stage is ambiguous because it only involves relatively lower temperatures (<30°C), to which the AHe method is not sensitive .
5.1. Onset of the Huaqiao-Ganquan Thrust Fault
In a thrust faulting system, uplift and accelerated erosion of the hanging wall would result in enhanced exhumation and cooling, which could be reflected by low-temperature thermochronologic data [10, 14, 42, 48, 56, 57]. Contrasting cooling histories of the TP and KB profiles show the initial rapid cooling events as a response to the onset of rapid unroofing during the middle-late Eocene (ca. 45–35 Ma) thrust faulting on the Huaqiao-Ganquan Fault (HGF) following an extended period of slow cooling during late Mesozoic to Paleocene time.
The TP profile appears to have resided longer in the partial retention zone and has a different thermal history when compared to the KB profile. Considering that samples are adjacent to major thrust faults (KMF and MLSZ) where moderate hydrothermal fluid activity might have occurred (Figure 3), cooling ages may reflect the influence of local hydrothermal activities rather than possible regional cooling [26, 27]. This is consistent with a proposed slow cooling from the late Cretaceous to ca. 40–50 Ma for the Huicheng Basin and its adjacent region . Second, the closure temperature of AFT is higher than that of the apatite (U-Th)/He method. We suggest that the TP apatite (U-Th)/He age versus elevation profile marks the rapid Paleocene cooling that occurred in the hanging wall of the regional major thrust faults (KMF and MLSZ). This is consistent with the accelerated exhumation in the northwestern WQB during the Eocene .
Differential rates of cooling during ca. 45–35 Ma between the KB and TP profiles suggest that the former experienced more pronounced cooling. Assuming an average geothermal gradient of ca. °C/km in the WQB, the present value for the study area, the difference in cooling is equivalent to a ca. 2–6 km exhumation, which is similar to the upward offset and throw of less than 6 km of the KB profile relative to the TP profile across the HGF (Figure 3). Due to strike-slip faulting in the region, there is much uncertainty about amounts of shortening across the Huicheng Basin/or some region. We tentatively constrain a minimum magnitude of shortening through area balancing of early Cretaceous sequences, using the flexural slip folding model and line-length constant method. The inferred amount of shortening is ca. 20 km shortening. Considering an average magnitude of shortening of ca. 7–10 km across the Jungong Basin in the northwestern part of the WQB , we estimate that the Cretaceous rocks were subjected to ca. 10–20 km of north-south shortening in the WQB (Figure 3). The fault plane of the HGF dips at >30°, and the cutoff angle in the hanging wall strata is relatively low (Figure 2(b)). We suggest a fault-propagation fold developing in the hanging wall over a thrust ramp (Figure 4(e)) . In summary, combined AHe age-depth data suggest rapid cooling since ~45 Ma, corresponding to the removal of ca. 3–5 km of rocks (Figure 4(e)).
5.2. Timing of Deformation in the West Qinling Belt
During the late Cretaceous to Eocene (ca. 100–50 Ma), the TP profile is characterized by a slow exhumation (Figures 4(a) and 4(b)), consistent with a documented regional peneplanation period. First, the absence of the late Cretaceous and Paleocene sedimentary deposits along the Qinling Belt implies that this orogenic belt experienced erosion during a period of time lacking tectonism [22, 40]. Second, apatite fission track data suggest slow cooling from the late Cretaceous to ca. 55 Ma for the Taibai Shan, northwestern Qinling Belt . Third, apatite fission track data show slow cooling (ca. 2°C/m.y.) in the interior of the Daba Shan (southern Qinling Belt), with ca. 2 km of exhumation from the late Cretaceous to Paleocene (Li et al., 2014). Fourth, the Wudang Massif (eastern Qinling Belt) is argued to have experienced tectonic quiescence for at least 50 Myr (late Cretaceous to Paleocene) . From satellite imagery, we could identify a low-relief erosion surface that cuts across Mesozoic and older rocks. We note this planation surface on the basis of the topographic slope (ca. 10–15°) over contiguous areas greater than 25–30 km2 (Figure 1(b) in Clark et al. ). A low-relief geomorphic erosion surface across the WQLF is also identified and argued to have formed since the late Cretaceous . Collective evidence indicates that the geomorphology of the Qinling Belt was initially controlled by long-term weathering and denudation, characterized by peneplanation that might have formed low-amplitude, long-wavelength topography since the late Cretaceous.
Since the Eocene, mountain building events have developed in the Qinling Belt . Within the WQB, during the early Cenozoic, the initiation of localized thrust-induced rapid cooling  and 40Ar/39Ar age of fault gouge  lead to the proposal of an Eocene onset of contractional tectonism along the north portion of the WQB (Figure 5(a)). Similarly, rapid cooling since ~45 Ma is recorded along the southeastern edge of the WQB. A major piece of the evidence for the Eocene deformation comes from the Oligocene depositional history (Guyuan Group) in the Oligocene [21, 34–36]. Our data also argue for an Eocene exhumation along the southeastern WQB (Huicheng Basin). Meanwhile, similarities in Cretaceous deformation are observed among basins in the WQB, including the Jungong Basin in the northwestern WQB  and Tanchang Basin in the middle of the WQB . Widespread Cretaceous deposits due to contractional deformation of certain origin occurred within the WQB in the early Cenozoic (ca. 45–35 Ma).
To orogenic scales, the angular unconformity between lower Cretaceous and Eocene is ca. 45–35 Ma in the interior WQB (Figure 1); meanwhile, the disconformity between the lower Cretaceous and Eocene strata in the East Qinling Belt which indicate that Qinling Belt is characterized with segmentation and the geodynamic setting variations for the West and East segment (Figure 1(a)). These results indicate that the middle-late Eocene crustal shortening probably represents a phase of hinterland regional deformation that was reactivated by the growth of the Tibetan Plateau.
5.3. Implications for Deformation of the NE Tibetan Plateau
Several lines of evidence suggest that the northeastern Tibetan Plateau might experience crustal thickening at ca. 50–30 Ma. First, thermochronological data indicate that Taibai Shan (northeastern WQB)  and Ganjia Shan (northwestern WQB)  have undergone rapid exhumation since ca. 45–50 Ma, causing the development of topography above the WQLF. Second, a shift in provenance in the Lanzhou Basin from the eastern Qilian Shan to the WQB suggests that the deformation and topographic growth might have started in the late Eocene . Third, paleomagnetic data indicate that Eocene sedimentary successions in the northeastern Tibetan Plateau were subsequently subjected to clockwise rotation in response to the growth of the WQB (Figure 5(a)) [59, 61]. Fourth, the appearance of detritus shed from mountain ranges suggests the development of high topography. Therefore, the appearance of WQB and Kunlun Shan detritus from the Lanzhou-Xining Basin and Qaidam Basin at ca. 25–30 Ma (Honggou section) could be related to the WQB and Kunlun Shan tectonic uplift since the late Eocene [9, 12, 16]. Fifth, the Songpan-Ganzi Terrane , Longmen Shan , and Micang Shan (south of the WQB)  experienced accelerated exhumation during the Eocene, as indicated by cooling ages (Figure 5(a)). Sixth, 40Ar/39Ar geochronological data of granitoids from the East Kunlun Fault suggest that rapid growth of the Kunlun range started at ~30 Ma, in response to the rising topography of the Kunlun range [14, 63].
Spatial-temporal variations in the nature of orogenies in the WQB and surrounding areas may be linked to the growth of the Tibetan Plateau (Figure 5(a)). The ongoing continent-continent collision may have reactivated preexisting faults in central Asia. Continuous crustal thickening and uplift south of the WQLF at ca. 30 Ma may have resulted in northward thinning of sediments in the Linxia-Xunhua Basin (Figure 5(a)) [11, 41]. During the same period, the Xining-Lanzhou Basin records shifting sedimentation rates, paleocurrent directions, and detrital zircon age populations (Figure 5(a)) [12, 16]. These changes are attributed to the pulsed growth of the WQB at ca. 30 Ma. In contrast to the southeastern Longzhong Basin, the Sikouzi Basin started to develop during the Eocene in extensional half-grabens . Collective evidence suggests that during the Oligocene, the WQB and the area to the south might have already formed a high plateau, with the Longzhong Basin (sensu lato) being a foreland basin associated with this topography (Figure 5(a)).
Integrating our results with previous studies allows us to reconstruct the paleolandscape of the northeastern Tibetan Plateau. Thermochronological data from the Jishi Shan, Laji Shan , Liupan Shan , Min Shan , and Longmen Shan  further indicate that rapid exhumation and river incision  occurred during the Miocene (Figure 5(b)). We compile exhumation data in the northeastern Tibetan Plateau with information on the distribution and thickness of Miocene to Quaternary deposits. We exclude data on uplift and deposits since the Neogene [11, 12, 16, 33, 64, 67] and restore paleolandscape in the northeastern Tibetan Plateau. This reconstruction shows that the WQLF was the northeastern boundary of the Tibetan Plateau during the Oligocene (Figure 5(a)).
Thermochronological data from the northeastern Tibetan Plateau and the evidence from the foreland basin are consistent with the northward propagation of thrusts during the Paleogene. During the Miocene, the onset of eastward extrusion, in addition to the continued N-S contraction, may have led to the growth of N-S-striking mountains (e.g., the Jishi Shan)  and the development of pull-apart basins (e.g., the Wushan Basin) along the WQLF (Figure 5(b)) . Eastward crustal extrusion might reactivated long-lived lithospheric zones of weakness for evacuating low-volume asthenospheric melts and the formation of Miocene (ca. 20 Ma) mafic-ultramafic magmatism [24, 68, 69]. Shear wave velocity data show that low-velocity mantle lithosphere and asthenosphere exist at depths of 90–150 km beneath the central Tibetan Plateau and the WQB . Altogether, they indicate both northward and eastward lateral growth of the northeastern Tibetan Plateau (Figure 5(b)).
Low-temperature thermochronologic data from the West Qinling Belt help constrain thermal histories across a major thrust system in the northeastern Tibetan Plateau. In the Huicheng Basin, apatite (U-Th)/He ages from the hanging wall (Cretaceous sedimentary strata) and the footwall wall (Triassic pluton) of the Huaqiao-Ganquan thrust fault reveal an initial rapid cooling event that occurred at ca. 45–35 Ma. Contrasting thermal histories were also observed across the fault. Our observations were interpreted as results of thrusting in the Eocene time (ca. 45–35 Ma) with a possible decade-kilometer crustal shortening. The Eocene age thrusting is broadly consistent with published sediment record and fault gouge data that are thought to reflect deformation distributed across the entire West Qinling Belt. This result clarifies that the northeastern Tibetan Plateau experienced deformation and possible mountain uplift in the Eocene, which might reflect widespread crust shortening associated with the initial growth of the Tibetan Plateau.
This data in the manuscript is our new data and is tested in the School of Earth Sciences, the University of Melbourne, and Commonwealth Scientific and Industrial Research Organization, Australia. Details of the methodology and results are presented in the supplementary materials.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work is supported by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (2019QZKK0901), National Key Research and Development Program of China (2017YFC1500101), National Natural Science Foundation of China (41590861, 41872204), and International Program for PhD Candidates, Sun Yat-sen University. We thank Dr. Zhi-Gang Li, Dr. Chuang Sun, Dr. Yang Wang, and PhD student Gan Cheng for comments and encouragement.