Major elements of rock-forming minerals were analyzed using an electron microprobe analyzer (EPMA, JEOL-JXA8100) at the IGG-CAS. The operating conditions were a 15-kV accelerating voltage, a 20-nA beam current, and 3-μm spot diameter. The counting time was 20 s at peak and 10 s at the lower and upper background positions, respectively. All data were corrected online using a modified ZAF (atomic number, absorption, and fluorescence) correction procedure. The detection limits (1σ) were in the range of 0.008–0.02wt%. The precision of the major element analysis was better than 1.0%.

Zircon grains were separated from samples by using conventional heavy liquid and magnetic separation techniques. Approximately 200 zircons of each sample were handpicked out from under a binocular microscope and mounted in epoxy resin together with zircon U-Pb age reference materials. Then the epoxy mounts were polished to expose the interior of the crystals. All grains were photographed in reflected and transmitted light under petrographic microscope to avoid fissures and inclusions. Cathodoluminescence (CL) and back-scattered electron (BSE) images were obtained using field emission scanning electron microscope (Nova NanoSEM 450) equipped with Gatan MomoCL4 at the Institute of Geology and Geophysics, Chinese Academy of Sciences (IGG-CAS). Mineral inclusions in zircon were identified using Raman spectroscopy at the IGG-CAS. Homogeneous zircons without inclusions or fissures were chosen for U-Pb dating and trace element analyses. Zircon U-Pb ages were determined by using secondary ion mass spectrometry (SIMS) and laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS). Zircon trace elements were acquired in the same run of age dating by LA-ICP-MS.

Zircon U, Th, and Pb isotope analyses were performed with a Cameca IMS-1280HR SIMS at the IGG-CAS. Detailed instrumental parameters and analytical procedures are described by Li et al. (2009). The analytical spot is about 20 × 30 μm in size. Plešovice zircon standard (²⁰⁶Pb/²³⁸U age of 337.3 Ma; Sláma et al., 2008) was interspersed with unknown grains. Nonradiogenic ²⁰⁴Pb was used for common Pb correction of measured compositions, and the present-day crustal common Pb composition is used in the model of Stacey and Kramers (1975). An inhouse zircon standard Qinghu was alternately analyzed together with other unknown zircons for the purpose of monitoring external uncertainties of SIMS U-Pb zircon dating calibrated by the Plešovice standard. Excel and the add-in Isoplot 2.49 program (Ludwig, 2001) were used for data calculation. Uncertainties on individual analyses in Supplementary Table S3 are reported at 1σ level. The weighted mean U-Pb ages are quoted with 95% confidence interval. Seven measurements on Qinghu zircon yielded a Concordia age of 159.7 ± 1.8 Ma, which is identical within errors with the recommended value of 159.5 ± 0.2 Ma (Li et al., 2013).

U-Pb dating combined with in situ trace element analyses of zircons were carried out in a single run by Agilent 7500a ICP-MS instrument equipped with Geolas-193 UV laser ablation system at the State Key Laboratory of Continental Dynamics, Northwest University, China. Operating conditions and data processing are described by Liu et al. (2007) in detail. The spot diameter is 32 μm for two eclogites (16PK190 and 16PK196) and 44 μm for gneiss (16PK181) with a laser repetition rate of 6 Hz. Helium was used as the carrier gas. Laboratory standards (GJ-1, 91500, NIST 610) were interspersed with unknown grains. U-Th-Pb isotope ratios and trace element contents were calculated using the GLITTER 4.0 program (Macquarie University). Harvard zircon 91500 was used as external standard. To monitor the external uncertainties of U-Pb zircon dating calibrated against 91500, a second zircon standard GJ-1 was analyzed as an unknown together with other unknown zircons. Trace element concentrations were calibrated by using ²⁹Si as internal standard and NIST 610 as external standard. Data reduction was carried out using Isoplot/Excel version 2.49 (Ludwig, 2001). Fifteen measurements on GJ-1 and 32 measurements on 91500 yielded weighted mean ²⁰⁶Pb/²³⁸U ages of 602 ± 3.7 Ma (1σ) and 1,062.5 ± 2.8 Ma (1σ), respectively, which are in good agreement with the recommended value of 608 ± 0.4 Ma (Jackson et al., 2004) for GJ-1 and 1,065.4 ± 0.6 Ma (Wiedenbeck et al., 2004) for 91500.

The ²⁰⁶Pb/²³⁸U weighted mean ages of ca. 284–280 Ma measured in the zircon cores from both eclogites represent the crystallization age of the magmatic protolith, which is demonstrated by their oscillatory zoning, high Th/U ratios, and steep REE patterns. This age is contemporaneous with the protolith formation of the Kaghan and Tso Morari UHP eclogites (Rajkumar, 2015; Rehman et al., 2016), which is ascribed to the eruption of the Panjal Traps in the north Indian continental margin (Shellnutt, 2018). Our results, coupled with the similar whole-rock chemistry and Nd isotopic compositions between the Stak eclogites and the Panjal Traps volcanic rocks (Kouketsu et al., 2017), indicate that the western Himalayan eclogites likely shared the same protoliths during the Permian magmatism of the Panjal Traps.

The rims of the eclogite zircons exhibit metamorphic characteristics and equilibrium texture with HP phases. The date of ca. 31 Ma most likely reflects eclogite-facies zircon recrystallization because, 1) the zircon rims have inclusions of garnet and rutile without plagioclase, 2) the disequilibrium texture of zircon with symplectitic diopside and plagioclase implies zircon inclusions presented in preexisting omphacite (Figure 4), and 3) the REE data indicate recrystallization in the presence of garnet and absence of plagioclase (Figure 6C). A similar metamorphic age (ca. 32 Ma) reported by Koukestsu et al. (2016) was interpreted as a recrystallization time after eclogite-facies metamorphism. However, the low Yb concentrations in those zircons only reflect equilibration with garnet, which could occur both under eclogite- and granulite-facies conditions. Moreover, metamorphic zircons from the country gneiss show REE characteristics and ages similar to those from the eclogites (Figures 6C,D). Therefore, the Stak massif was coherently buried beneath the Asian plate and underwent eclogite-facies metamorphism at ca. 31 Ma.

The HP metamorphic time of the Stak massif is at least 15 Ma later than the UHP metamorphic time (53–46 Ma) in the Kaghan and Tso Morari massifs (de Sigoyer et al., 2000; Kaneko et al., 2003; Donaldson et al., 2013). Therefore, the HP-UHP massifs in the western Himalaya were asynchronously buried to ∼2.5–3.8 GPa, indicating a continuous subduction/collision process spanning more than 15 Ma. After the Kaghan and Tso Morari massifs experienced UHP metamorphism (∼46 Ma) and rapid exhumation to the middle to lower crustal levels (∼44–40 Ma) (de Sigoyer et al., 2000; Parrish et al., 2006; Wilke et al., 2010), the Stak massif was buried to HP eclogite-facie conditions at ca. 31 Ma. Such a continuous burial process recorded by HP–UHP rocks is common in other collisional belts such as Dabie-Sulu (e.g., Zhang et al., 2009), Alps (e.g., Berger and Bousquet, 2008), and the Western Gneiss Region (e.g., Kylander-Clark et al., 2009). In the following section, we will focus on how this process proceeded through time.

The HP metamorphism of the Stak massif is broadly contemporaneous with that of quartz eclogites, HP granulites, and garnet amphibolites (mostly 40–25 Ma) in the central and eastern Himalayas (Figure 7A) (Corrie et al., 2010; Zhang et al., 2010; Rubatto et al., 2013; Wang et al., 2021). This indicates that the subducted Indian continental slab probably underwent a coherent process along the Indus–Yarlung Suture Zone since 40 Ma. However, UHP metamorphism spatially restricted to the western Himalaya (Kaghan and Tso Morari) was only recognized before ca. 46 Ma (Kaneko et al., 2003; Donaldson et al., 2013). The discrepancy of HP–UHP timing in the western and central Himalayas could be induced by continental subduction of a jagged Indian margin (Guillot et al., 2008), by different exhumation styles (residence times at crustal level) of eclogites (O’Brien, 2019; Rehman, 2019), or by changes in subduction dip angle (Figure 2). However, it is unlikely that an irregular northern margin of India would trigger the long discrepancy (>15 Ma) between UHP and HP peak metamorphic ages. The different P-T-t paths observed between UHP and HP eclogites were interpreted to result from the prolonged residence times in the central and eastern Himalayas but rapid exhumation of UHP rocks in the west (Rehman, 2019). However, the different exhumation styles are not the first-order mechanism responsible for the distinct peak metamorphic ages of HP–UHP rocks. The ∼31 Ma eclogite-facies metamorphic age in Stak is similar to those recorded in the central Himalayan HP eclogite (∼30 Ma; Wang et al., 2021) and the Zanskar HP granulite (31–28 Ma; Vance and Harris, 1999) and the Ayilari garnet amphibolite (∼32 Ma; Chen et al., 2021) in the western Himalaya (Figure 1B), indicating that these HP rocks were synchronously buried to shallow depths much later than the UHP rocks.

Notably, the old UHP eclogites and young HP rocks record different peak-pressure thermal gradients (temperature change with depth, referred to as metamorphic T/P) (Figure 8), implying a change in the subduction geometry of the Indian continental slab during 46–40 Ma. The subduction of oceanic lithosphere (Neo-Tethys) after the initial India–Asia collision (∼60 Ma; e.g., Parsons et al., 2020) would continuously release fluids or melts, generate convective corner flow in the mantle wedge, and, thus, induce magmatism at convergent plate margins (Figure 7B). The successive steep subduction of Indian continental crust driven by slab pull of the Neo-Tethys underwent UHP metamorphism (low T/P) at 53–46 Ma (Figures 2A and 9A). The steep subduction of the Indian continental slab is also supported by the presence of the UHP eclogite units in the western Himalaya now directly adjacent to the Indus–Tsangpo suture zone (O’Brien, 2019). However, after ∼40 Ma, the sole occurrence of HP metamorphism (high T/P) along the strike reflects an increase of geothermal gradient that requires a change in subduction dip angle (Figures 2B and 9B). These young HP rocks most likely resulted from low-angle underthrusting during continental collision (Guillot et al., 2008; Soret et al., 2021). Therefore, the change in subduction dip angle of the Indian continental slab during 46–40 Ma is critical for the transition from continental subduction to collision dynamics.

The change in subduction dip angle can be induced by the continuous subduction of buoyant continental lithosphere (Parsons et al., 2021) or by the breakoff of the oceanic (Neo-Tethyan) lithosphere (Davies and von Blanckenburg, 1995; Magni et al., 2017). Breakoff of the Neo-Tethyan lithosphere would lead to fast exhumation of UHP rocks (Kohn and Parkinson, 2002), uplifting of the buoyant continental lithosphere, and shifting to gently dipping subduction or underthrusting (Davies and von Blanckenburg, 1995; Chemenda et al., 2000). In addition, slab breakoff mostly occurs at greater depths (∼130–240 km) than the base of the overriding lithosphere and would not trigger significant magmatism at convergent margins (e.g., Freeburn et al., 2017), which is consistent with the weak magmatism within the Asian plate at 46–40 Ma (Figure 7B). There are many lines of evidence supporting the breakoff of the Neo-Tethyan slab around 46–40 Ma after the initial India–Asia collision. 1) According to paleomagnetic reconstructions and tomographical observations, Negredo et al. (2007) estimated the breakoff of oceanic slab occurring at 48–44 Ma. 2) The rapid exhumation of the Kaghan UHP eclogite from ∼100 to ∼35 km during 46–44 Ma (Parrish et al., 2006) also agrees with the breakoff model (Kohn and Parkinson, 2002). 3) Based on tomography, the Asian tectonics reconstructions and the Indian plate kinematics, Replumaz et al. (2014) proposed that the major breakoff between India and the Tethys ocean occurred at ∼45 Ma. However, Parsons et al. (2021) proposed a different timing for slab breakoff based on interpretations of the same tomographic anomalies. 4) The occurrence of 42–40 Ma intraplate-type mafic dykes in eastern Tibet supports the breakoff of subducting Neo-Tethyan slab from the Indian continental slab during the middle Eocene (Xu et al., 2008). 5) The 45 Ma oceanic island basalt-type gabbro in southern Tibet was used to constrain the breakoff time of the Neo-Tethyan slab (Ji et al., 2016). Numerical studies suggest that the timing of slab breakoff after the onset of continental collision varies from 10 to 25 Ma, which is largely affected by the strength and age of the subducting oceanic slab (e.g., van Hunen and Allen, 2011). Our proposed slab breakoff timing (46–40 Ma) for the old Neo-Tethyan slab after the onset of continental collision (∼60 Ma) broadly fits this modeled result. Therefore, we suggest that the different T/P and ages between the Himalayan UHP and HP rocks were induced by the breakoff of the Neo-Tethyan slab and subsequent transition from steep continental subduction to low-angle underthrusting (collision) at 46–40 Ma (Figure 9B).

The synchronous HP metamorphism induced by laterally large-scale low-angle underthrusting is also supported by palaeomagnetic and magmatic data. Since the Late Cretaceous, the India–Asia convergence rates in the eastern and western Himalayas show high consistency (Figure 7C), decreasing from 140–160 to 80–100 mm/year at 52–50 Ma and then decreasing to 40–60 mm/year at ca. 45 Ma (van Hinsbergen et al., 2011). The first deceleration is attributed to the initial India–Asia collision (van Hinsbergen et al., 2011), whereas the second most likely resulted from the loss of oceanic slab pull after breakoff (e.g., Bercovici et al., 2015; Ji et al., 2016) or the India–Asia second collision (45–40 Ma) after an earlier collision (∼60 Ma) of either Indian continent + intraoceanic arc or greater India microcontinent + Asia (e.g., Parsons et al., 2020, 2021). A low convergence rate benefits sufficient heat exchange between the subducting plate and the overlying mantle, thus inducing HP metamorphism with high T/P (Guillot et al., 2008). In addition, the Himalayan magmatic activities from east to west were systematically complementary to the HP metamorphic events after 40 Ma. The Kohistan–Ladakh–Gangdese arc system and south Karakoram exhibit a coherent magmatic lull at 40–25 Ma (Figure 7B). The Stak eclogite in the western Himalaya and most of the other HP rocks in the central and eastern Himalaya also formed at this stage (Figures 7A and 9C). The low-angle underthrusting after slab breakoff would have driven the asthenosphere beneath the Asian plate away and then shielded the active continental margin from convective heat, which would have led to the magmatic lull in the convergent margin (Chung et al., 2005; Ji et al., 2016; Ji et al., 2020). The driving force for the continuous convergence and underthrusting after the breakoff of the Neo-Tethyan slab could be subduction of Australian oceanic lithosphere (e.g., Parsons et al., 2021) or subduction of Indian continental lithosphere (e.g., Capitanio et al., 2010). In addition, the continuous low-angle underthrusting during continental collision at 40–25 Ma potentially enhanced crustal thickening (Soret et al., 2021) and elevated the thermal structure of orogenic wedge due to radiogenic heat production in the thickened crust (Berger et al., 2011), contributing to the high- or ultrahigh-temperature overprints in HP eclogites during exhumation (Wang et al., 2021). In this regard, the change in subduction dip angle of continental slab would have influenced both magmatic activities and peak geothermal gradients of HP–UHP rocks. To summarize, a lateral (>2,500 km) change in the subduction geometry of the Indian continental slab during the middle Eocene would be broadly consistent with the igneous and metamorphic records.