Rock sample was crushed and washed, followed by conventional heavy-liquid and magnetic separation to separate zircon grains. Representative zircon grains were handpicked under a binocular microscope and then mounted in epoxy resin and polished to a smooth flat surface. Sample processing was conducted at the Langfang Regional Geological Survey (Hebei Province, China). Prior to U-Pb dating, internal factures of zircons were observed using transmitted, reflected, and cathodoluminescence images. These images were obtained at the Wuhan Sample Solution Analytical Technology Co., Ltd. (Hubei Province, China).

U-Pb dating and trace element analysis of zircon were simultaneously conducted by LA-ICP-MS at the Wuhan Sample Solution Analytical Technology Co., Ltd. Detailed operating conditions for the laser ablation system and the ICP-MS instrument and data reduction are the same as description by Zong et al. (2017). Laser sampling was performed using a GeolasPro laser ablation system that consists of a COMPexPro 102 ArF excimer laser (wavelength of 193 nm) and a MicroLas optical system. An Agilent 7700e ICP-MS instrument was used to acquire ion-signal intensities. Helium was applied as a carrier gas. Argon was used as the make-up gas and mixed with the carrier gas via a T-connector before entering the ICP. A “wire” signal smoothing device is included in this laser ablation system (Hu et al., 2015). The spot size and frequency of the laser were set to 32 µm and 8 Hz, respectively, in this study. Zircon 91,500 and glass NIST610 were used as external standards for U-Pb dating and trace element calibration, respectively. Each analysis incorporated a background acquisition of approximately 20–30 s followed by 50 s of data acquisition from the sample. An Excel-based software ICPMSDataCal was used to perform off-line selection and integration of background and analyzed signals, time-drift correction and quantitative calibration for trace element analysis and U-Pb dating (Liu et al., 2008, 2010). Concordia diagrams and weighted mean calculations were made using Isoplot/Ex_ver3 (Ludwig, 2003).

Experiments of in situ Hf isotope ratio analysis were conducted using a Neptune Plus MC-ICP-MS (Thermo Fisher Scientific, Germany) in combination with a Geolas HD excimer ArF laser ablation system (Coherent, Göttingen, Germany) that was hosted at the Wuhan Sample Solution Analytical Technology Co., Ltd. All the Hf analyses were done on the same or equivalent spots as those for U-Pb laser ablation analyses. A stationary laser ablation spot with a beam diameter of 32 μm was used for the analyses. Detailed operating conditions for the laser ablation system and the MC-ICP-MS instrument and analytical method were described in Hu et al. (2012). The major limitation to accurate in situ zircon Hf isotope determination by LA-MC-ICP-MS is the very large isobaric interference from ¹⁷⁶Yb and, to a much lesser extent ¹⁷⁶Lu on ¹⁷⁶Hf. It has been shown that the mass fractionation of Yb (β_Yb) isn’t constant over time and that the β_Yb that is obtained from the introduction of solutions is unsuitable for in situ zircon measurements (Woodhead et al., 2004). The under- or over-estimation of the β_Yb value would undoubtedly affect the accurate correction of ¹⁷⁶Yb and thus the determined ¹⁷⁶Hf/¹⁷⁷Hf ratio. We applied the directly obtained β_Yb value from the zircon sample itself in real-time in this study. The ¹⁷⁹Hf/¹⁷⁷Hf and ¹⁷³Yb/¹⁷¹Yb ratios were used to calculate the mass bias of Hf (β_Hf) and Yb (β_Yb), which were normalized to ¹⁷⁹Hf/¹⁷⁷Hf = 0.7325 and ¹⁷³Yb/¹⁷¹Yb = 1.132685 (Fisher et al., 2014) using an exponential correction for mass bias. Interference of ¹⁷⁶Yb on ¹⁷⁶Hf was corrected by measuring the interference-free ¹⁷³Yb isotope and using ¹⁷⁶Yb/¹⁷³Yb = 0.79639 (Fisher et al., 2014) to calculate ¹⁷⁶Yb/¹⁷⁷Hf. Similarly, the relatively minor interference of ¹⁷⁶Lu on ¹⁷⁶Hf was corrected by measuring the intensity of the interference-free ¹⁷⁵Lu isotope and using the recommended ¹⁷⁶Lu/¹⁷⁵Lu = 0.02656 (Blichert-Toft et al., 1997) to calculate ¹⁷⁶Lu/¹⁷⁷Hf. We used the mass bias of Yb (β_Yb) to calculate the mass fractionation of Lu because of their similar physicochemical properties. Off-line selection and integration of analyte signals, and mass bias calibrations were performed using ICPMSDataCal (Ludwig, 2003).

Similar to zircons, monazite separation and sample processing were conducted at the Langfang Regional Geological Survey. Transmitted, reflected, and BSE images, which were obtained at the Wuhan Sample Solution Analytical Technology Co., Ltd (Hubei Province, China), were used to determine internal features of monazite grains. Twenty-five eu- to subhedral monazite grains without obvious inclusions were selected for in situ LA-ICP-MS isotopic and compositional analyses. U-Pb dating of monazite was conducted by LA-ICP-MS at the Wuhan Sample Solution Analytical Technology Co., Ltd. Laser sampling was performed using a GeolasPro laser ablation system that consists of a COMPexPro 102 ArF excimer laser (wavelength of 193 nm) and a MicroLas optical system. An Agilent 7700e ICP-MS instrument was used to acquire ion-signal intensities. Helium was applied as a carrier gas. Argon was used as the make-up gas and mixed with the carrier gas via a T-connector before entering the ICP. A “wire” signal smoothing device is included in this laser ablation system, by which smooth signals are produced even at very low laser repetition rates down to 1 Hz (Hu et al., 2015). It is very useful for in-situ U-Pb dating of high-U mineral (Zong et al., 2015). The spot size and frequency of the laser were set to 16 µm and 2 Hz, respectively. The laser energy was set to 80 mJ. Monazite standard 44069 and glass NIST610 were used as external standards for U-Pb dating and trace element calibration, respectively. Each analysis incorporated a background acquisition of approximately 20–30 s followed by 50 s of data acquisition from the sample. An Excel-based software ICPMSDataCal was used to perform off-line selection and integration of background and analyzed signals, time-drift correction and quantitative calibration for trace element analysis and U-Pb dating (Liu et al., 2008, 2010). Concordia diagrams and weighted mean calculations were made using Isoplot/Ex_ver3 (Ludwig, 2003).

Major and trace element contents were analyzed at the Beijing Research Institute of Uranium Geology. The fresh rock samples were chipped and powdered to a mesh size of ∼200 using a tungsten carbide ball mill. The details of the analytical procedures are described by Gao et al. (2003). Major oxide analyses were conducted using a PANalytical Axios MAX X-ray Fluorescence Spectrometer with an analytical uncertainty of <5%. The trace element concentrations were analyzed using a Perkin-Elmer NexIon 300D ICPMS with an analytical precision of <1% for elements with concentrations >200 ppm and 1%–3% for elements with concentrations <200 ppm.

Whole-rock Sr-Nd-Pb isotopes were analyzed at the Beijing Research Institute of Uranium Geology. Approximately 200 mg of each sample powder was dissolved in an HF + HNO₃ acid mixture for 48 h in a Teflon beaker. All samples were prepared in duplicate. The digests were dried, dissolved in hydrochloric acid and heated in closed vials at 160 C for 1 h, and evaporated to dryness. The Sr and Nd were separated and purified by conventional cation-exchange techniques. The initial ⁸⁷Sr/⁸⁶Sr ratios and ε_Nd (t) values at the time of crystallization were calculated from the weighted mean zircon U-Pb age and the Rb, Sr, Sm and Nd contents. When calculating the ε_Nd (t) values, we assumed the model composition of a chondritic uniform reservoir at the estimated age. The T_DM1 and T_DM2 values are the estimated ages of extraction from the depleted mantle according to the one-stage and two-stage crustal pre-histories, respectively, as assumed by Depaolo (1988) and Depaolo et al. (1991). The precision of the calculated initial ⁸⁷Sr/⁸⁶Sr values is limited by the errors in the parent/daughter ratios calculated from the geochemical data. Nevertheless, in most cases, errors as high as ±10% in the Rb/Sr ratio introduce uncertainties below 0.0001 in the initial ⁸⁷Sr/⁸⁶Sr values (Dolgopolova et al., 2013), and a 10% error in the Nd isotopes results in an uncertainly of about 0.4–0.7 in the ε_Nd (t) values. The measured Pb isotopic ratios were corrected for the instrumental mass fractionation of 0.1 amu-1 by referencing to repeat analyses of the standard NBS-981.

U-Pb dating and trace elements results for zircons of the Liemai two-mica granite are shown in appendix (Supplementary Table S1). Cathodoluminescence (CL) images of representative zircon grains from the Liemai two-mica granite are shown in Figure 4A. Eu- to subhedral zircon crystals are long prismatic (80–200 μm), with aspect ratios of 2:1 to 5:1. Most zircon grains exhibit core-(mantle)-rim structure with gray subhedral to anhedral cores and dark rims (Figure 4A). The inherited cores normally show obvious oscillatory growth zoning (Figure 4A), suggesting magmatic origin (Wu and Zheng, 2004). In contrast, some rims display weak but visible oscillatory overgrowth zoning indicating magmatic origin, some others show no oscillatory overgrowth zoning and are homogeneous suggesting metamorphosed origin (Figure 4A; Wu and Zheng, 2004; Hoskin, 2005). Detailed description of zircon structure and analytical position are shown in appendix (Supplementary Table S1). Twenty-three analyses (with a concordance of >90%) yield ²⁰⁶Pb/²³⁸U ages ranging from 41.6 to 1721.5 Ma. Fifteen analyses on inhibited cores yield ²⁰⁶Pb/²³⁸U ages ranging from 130.0 to 1721.5 Ma. Eight analyses on dark rims or zircons without inherited cores yield ²⁰⁶Pb/²³⁸U ages ranging from 41.6 to 43.6 Ma with a weighted Eocene mean age of 43.0 ± 0.5 Ma (MSWD = 0.7) (Figure 4B).

Thirteen spots were also analyzed for Lu-Hf isotopes. The analytical results are listed in appendix (Supplementary Table S2). The analyzed zircon grains have ¹⁷⁶Lu/¹⁷⁷Hf ratios of 0.000018–0.003409 and ¹⁷⁶Hf/¹⁷⁷Hf ratios of 0.281963–0.282520. Except for one analysis (D5126-06-24), zircon crystals display ¹⁷⁶Lu/¹⁷⁷Hf ratios <0.002, indicating extremely low accumulations of the radioactive Hf isotope after crystallization (Wu et al., 2007). Seven analyses on juvenile zircon rims yield calculated ɛ_Hf (t) (t = 43 Ma) values of −14.3 to −8.0 (average of −11.1) with two stage Hf model ages (T_DM2) of 1,027–1,256 Ma.

The monazite U-Th-Pb dating and trace elements results for the Liemai two-mica granite are listed in appendix (Supplementary Table S3). The BSE images of representative monazites are showed in Figure 4C. All monazite grains are grey in BSE images, and most of them are homogeneous, although some grains show zoning (Figure 4C). The monazites in the Liemai two-mica granite have high contents of Th (37,794∼113,394 ppm) and U (1,663∼9,254 ppm), and relatively low Pb contents (83∼245 ppm). The ²⁰⁷Pb/²³⁵U ages of these monazites are scattered, ranging from 59.3 to 251.9 Ma. In contrast, the ²⁰⁶Pb/²³⁸U ages (42.2∼44.9 Ma) and ²⁰⁸Pb/²³²Th ages (43.1∼44.8 Ma with a weighted Eocene mean age of 43.7 ± 0.2 Ma, N=25, MSWD = 0.7; Figure 4D) are very concentrated. In addition, the ²⁰⁶Pb/²³⁸U ages and ²⁰⁸Pb/²³²Th ages are very concordant with each other (Figure 4D).

The whole-rock major and trace element data for the Liemai two-mica granite are listed in appendix (Supplementary Table S4). Samples are characterized by high contents of SiO₂ (66.88–67.48 wt.%), Al₂O₃ (16.24–16.44 wt.%), Na₂O (3.72–3.83 wt.%) and K₂O (3.03–3.27 wt.%), and relatively low contents of CaO (2.42–2.52 wt%), total iron (TFeO = 2.06–2.11 wt.% TiO₂ (0.31–0.33 wt.%), MnO (0.03–0.04 wt.%), MgO (1.24–1.29 wt.%), and P₂O₅ (0.13–0.14 wt.%) with relatively high Na₂O/K₂O mass ratios (1.14–1.25) and Mg# [Mg# = MgO×100/(MgO+FeO^T) (molar ratio)] values (51.75–52.20) (Supplementary Table S4). The total alkali contents (K₂O+Na₂O) range from 6.82 to 6.99 wt%, indicative of sub-alkaline series in the total alkalis vs. silica (TAS) diagram (Figure 5A; Middlemost, 1994). The A/CNK values [A/CNK = Al₂O₃/(CaO + Na₂O+ K₂O) (molar ratio)] vary from 1.15 to 1.17, whereas the A/NK values [A/NK = Al₂O₃/(Na₂O+ K₂O) (molar ratio)] range from 1.69 to 1.72, indicating peraluminous features (Figure 5B; Maniar and Piccoli, 1989). On the K₂O vs. SiO₂ diagram, the samples show high-K calc-alkaline characteristics (Figure 5C; Peccerillo and Taylor, 1976). On the A-F-C diagram, the Liemai two-mica granite samples plot in the peraluminous granitoid-type field (Figure 5D; Chappell and White, 1992).

The Liemai two-mica granite is characterized by LREE (light rare Earth element) enrichment and HREE (heavy rare Earth element) depletion with high LREE/HREE ratios (16.0–16.7; Supplementary Table S4), suggesting significant differentiation of LREE from HREE. This is in accordance with the right-leaning chondrite-normalized REE patterns (Figure 6A). The Liemai two-mica granite displays an indistinctive negative Eu anomaly (Eu/Eu*=0.7∼0.8; Supplementary Table S4; Figure 6A). The rock is enriched in LILE (large ion lithophile element) such as Rb and K and depleted in HFSE (high field strength element) such as Nb, Ta, Zr, Hf, and Ti (Figure 6B). Importantly, the Liemai two-mica granite has low concentrations of Y (6.9∼7.4 ppm) and Yb (∼0.5 ppm) and relatively high concentrations of Sr (269∼278 ppm) and La (32.3∼36.0 ppm) with high Sr/Y (36.8∼39.2) and (La/Yb)_N (44.5∼50.9) ratios (subscript N represents chondrite-normalized value based on Sun and McDonough, 1989), showing adakitic features (Figures 6C, D; Martin, 1986; Defant and Drummond, 1990).

The whole-rock Sr-Nd-Pb isotopic compositions of the Liemai two-mica granite are listed in appendix (Supplementary Table S5). The samples have low (⁸⁷Sr/⁸⁶Sr)_i ratios (0.718414–0.718423) and negative ε_Nd (t) values (–14.7 to –13.9), with T_DM2 ranging from 1,803 to 1,852 Ma. The (²⁰⁶Pb/²⁰⁴Pb)_i, (²⁰⁷Pb/²⁰⁴Pb)_i, and (²⁰⁸Pb/²⁰⁴Pb)_i values are 18.781–18.791, 15.687–15.695 and 39.214–39.236, respectively. All initial values are calculated for 43 Ma.

As mentioned above, the juvenile zircon rims with oscillatory zoning in the Liemai two-mica granite are magmatic while those without oscillatory zoning are of metamorphosed origin. Interestingly, these rims have undistinguishable ²⁰⁶Pb/²³⁸U ages (Supplementary Table S1). This indicates that metamorphism and partial melting of the source rock for the Liemai two-mica granite were possibly coeval. This phenomenon has been recorded in the Yardoi dome (Zeng et al., 2011). Therefore, the weighted ²⁰⁶Pb/²³⁸U age of juvenile zircon rims (43.0 ± 0.5 Ma) can be interpreted as the timing of partial melting of the source rock and the formation age of the Liemai two-mica granite.

Monazite incorporates a significant amount of ²³⁰Th into its crystal structure during crystallization, which will lead to generation of “excess” ²⁰⁶Pb and result in disequilibrium in the ²³⁸U→²⁰⁶Pb decay. In juvenile monazite (Cenozoic), the “excess” ²⁰⁶Pb will lead to apparent ²⁰⁶Pb/²³⁸U ages that are older than measured ²⁰⁷Pb/²³⁵U and ²⁰⁸Pb/²³²Th ages. Besides, due to their young age, such monazites normally contain low contributions of radiogenic ²⁰⁷Pb, which will result in ²⁰⁷Pb/²³⁵U ages that are less precise than ²⁰⁶Pb⁄²³⁸U and ²⁰⁸Pb/²³²Th ages. In contrast, owing to high contents of ²³²Th (wt.% levels) and therefore significant ²⁰⁸Pb in monazite, ²⁰⁸Pb/²³²Th ages can be taken as the best estimate of the monazite crystallization age (Cottle et al., 2015). Therefore, the weighted ²⁰⁸Pb/²³²Th age (43.7 ± 0.2 Ma) are interpreted as crystallization age of monazites in the Liemai two-mica granite. Owing to relatively low closure temperature and high sensitivity, inherited monazites are very rare even in low-temperature granites (Parrish, 1990). Thus, the weighted ²⁰⁸Pb/²³²Th age (43.7 ± .2 Ma) of monazites in the Liemai two-mica granite can be interpreted as formation age of the rock.

In summary, the weighted ²⁰⁶Pb/²³⁸U age of juvenile zircons are consistent with the weighted ²⁰⁸Pb/²³²Th age of monazites within errors, indicating that the Liemai two-mica granite was generated at ∼43.5 Ma, coeval to two-mica granites in Yardoi, Dala, Quedang (Figure 2).

A pronounced magmatic flare-up with intensive input of mantle materials along the southern Gangdese at ca. 52–51 Ma (Zhu et al., 2015), together with a sudden drop of the India-Asia convergence rate (Patriat and Achache, 1984; van Hinsbergen et al., 2011), marked the slab breakoff of the subducting Neo-Tethyan oceanic lithosphere. In addition, the crustal thickness beneath the Gangdese reached a thickness>50 km at the beginning of the Eocene (Zhu et al., 2017; Tang et al., 2020). These observations signify that the Neo-Tethys Ocean has been completely closed prior to the Eocene and the Gangdese already had an abnormally thick crust.

Ultra-high pressure metamorphism (∼55 Ma) and high amphibolite facies to granulite facies metamorphism (∼43–47 Ma) in the Tethyan Himalaya suggest that intensive collision and crustal thickening occurred following closure of the Neo-Tethyan Ocean (Zeng et al., 2011; Gao et al., 2012). This is followed by intensive anatexis of the amphibolitic lower crust to generate Early Eocene adakitic magmas in the Yardoi area (Zeng et al., 2011, 2015; Hou et al., 2012; Dai et al., 2020a). It is worth mentioning that the garnet amphibolites, as the major source rocks for the adakitic rocks, show negligible Eu anomalies and differentiation of REE with low Sr/Y values (< 5) (Figure 6; Hou et al., 2012). In contrast, the garnet amphibolite-derived adakitic rocks are obviously depleted in MREE, HREE and Y, in despite of their comparable LREE and Sr contents and Eu anomalies to the garnet amphibolites (Figure 6). Therefore, MREE-enriched amphibole and HREE-enriched garnet should serve as residual minerals while the plagioclase should serve as melting phase during partial melting to produce melts with remarkable depletion of MREE and HREE and high Sr/Y ratios. Therefore, generation of these Early Eocene adakitic granites in the Yardoi area are indicative of a thickened crust.

In the past, the Sr/Y ratio was used to qualitatively estimate the paleo-crustal thickness. Recently, Chiaradia (2015) found that the Sr/Y ratios of magmatic rocks show positive correlation with crustal thickness in young arcs. Subsequently, Chapman et al. (2015) and Profeta et al. (2015) reconstructed the global and regional correlations between the whole-rock Sr/Y and (La/Yb)_N ratios of intermediate-felsic arc magmatic rocks and arc crustal thickness, providing quantitative constraints on the paleo-crustal thickness in ancient arcs. More recently, Hu et al. (2017) successfully extended this method to continental collisional zones. The Eocene granites in the Yardoi area were formed posterior India-Eurasia collision and slab breakoff of the Noe-Tethys oceanic lithosphere thus can be used to estimate the paleo-crustal thickness using the equations from Hu et al. (2017). Interestingly, granite sills (e.g., Liemai two-mica granite) or small intrusions (e.g., Quedang two-mica granite) show relatively lower Sr/Y ratios than big intrusions in Yardoi and Dala (Figure 14A). This is probably because small volume of magmas cooled faster than big ones, thus preserved relatively primitive Sr/Y values. In contrast, magmas with large volume could remain a relative long period at crystal mush state, which could result in elevated Sr/Y ratios due to extraction of low-Sr/Y melts. Therefore, although the extraction of low-Sr/Y melts (LMGs) could have increased the Sr/Y ratios of the crystal mush (HMGs) to some degree (Gao et al., 2021), the small intrusions (e.g., Liemai and Queang two-mica granites) with relatively primitive Sr/Y values can be preferably used to estimate the paleo-crustal thickness. This is supported by the close sample plotting sites of these granites to that of the assumed initial melt in the Sr vs. Ba diagram (Figure 13). The calculated CT_Sr/Y values (crustal thickness based on whole-rock Sr/Y ratios) are broadly range from 50 to 90 km (Figure 14B) using the equation from Hu et al. (2017). Thus, the crustal thickness beneath the Himalaya was inferred to ∼50 km during the Early Eocene.

Therefore, both sides of the IYTS experienced notable crustal thickening and gained an abnormally thick crust (∼50 km) in the Early Eocene due to the India-Asia collision. Importantly, ascent of hot asthenospheric materials though slab windows caused by slab breakoff of the Neo-Tethys oceanic lithosphere triggered generation of mafic magmas (Ji et al., 2016). Mafic magmas that underplated to the crust-mantle boundary caused partial melting of the thickened lower crust to form adakitic rocks while those intruded in the middle-upper crust formed mafic rocks (Figure 15).