Conical intersection–regulated intermediates in bimolecular reactions: Insights from C(1D) + HD dynamics

See allHide authors and affiliations

Science Advances  26 Apr 2019:
Vol. 5, no. 4, eaaw0446
DOI: 10.1126/sciadv.aaw0446


The importance of conical intersections (CIs) in electronically nonadiabatic processes is well known, but their influence on adiabatic dynamics has been underestimated. Here, through combined experimental and theoretical studies, we show that CIs induce a barrier and regulate conversion from a precursor metastable intermediate (CI-R) to a deep well one. This results in bond-selective activation, influencing the adiabatic dynamics markedly in the C(1D) + HD reaction. Theory is validated by experiment; quantum dynamics calculations on highly accurate ab initio potential energy surfaces yield rate coefficients and product branching ratios in excellent agreement with the experiment. Quasi-classical trajectory calculations reveal that the CI-R intermediate leads to unusual reaction mechanisms (designated as C─H activation complex conversion and cyclic complex), which are responsible for large branching ratios. We also reveal that CI-R intermediates exist in other reactive systems, and the dynamical effects uncovered here may have general significance.


Conical intersections (CIs) are important in many chemical and biological systems (17), ranging from elementary reactions involving three to five atoms (4, 5) to biomolecules such as DNA (6) and proteins (7). In particular, CI-induced ultrafast nonadiabatic dynamics have attracted substantial interest over the past two decades (711). In contrast, the effects of CIs in adiabatic reaction dynamics have been mostly neglected when nonadiabatic effects are small. Here, we find that the CI between the bA1 and dA1 states in the C(1D) + HD reaction (Fig. 1A) induces a small barrier and a shallow precursor well (Fig. 1B) before the deep well on the lower potential energy surface (PES). These features have a profound effect on the adiabatic reaction dynamics. The shallow well supports a metastable-state complex corresponding to a cyclic CHD structure, whose stability is regulated by the CI-induced barrier, and thus, it is designated as a CI-regulated (CI-R) intermediate here. We further reveal that the CI-R intermediate triggers bond-selective activation through two unusual reaction mechanisms and has a strong influence on product branching ratios of the C(1D) + HD reaction. These findings are achieved through combined experimental and theoretical studies, including low-temperature experiments using a supersonic flow reactor and extensive dynamical calculations on highly accurate ab initio PESs. These results are important for acquiring a deeper understanding of complex-forming reactions and may provide us with a potential method to control chemical reactions by adjusting the location of CIs with lasers.

Fig. 1 Schematic reaction profile and C2v conical interaction.

(A) The bA1 potential profile along the reaction path for the C(1D) + HD reaction. MECP, minimum energy crossing point. (B) A two-dimensional view of the C2v CI between the bA1 and dA1 states. All indicated values and potentials are from our ab initio calculations, the details of which along with those about the g-h coordinates are given in the Supplementary Materials. The CI-R metastable-state intermediate is shown. The bond lengths are in bohr, and the angles are in degrees.

It is well known from the current understanding of chemical reactions that the transition state determines the rate coefficient of a direct-activated reaction, whereas the deep well intermediate complex dominates that of a complex-forming reaction, forming the basis of many approximate theories of chemical reaction rates. Here, we demonstrate a more complicated and interesting case, in which a precursor intermediate plays an important role and traditional approximate theories fail. This case involves multiple maxima and minima in the PES, presenting “the most severe challenge” (12) to PES construction and quantum dynamics (QD) calculations (13).

Extensive experimental and theoretical efforts have been devoted to the C(1D) + H2 reaction and its isotopic variants. Experimentally, several kinetic and dynamical measurements (1419) have been reported. In particular, Fisher et al. (14) and Sato et al. (15) measured the CD/CH branching ratio of the C(1D) + HD reaction through pump-probe techniques. Moreover, intramolecular isotopic effects in A + HD–type insertion reactions under crossed beam or thermal conditions have been discussed by Liu (16). Recently, Hickson et al. (18, 19) measured low-temperature rate coefficients for the C(1D) + H2 (D2) reactions using the supersonic flow reactor technique. Theoretically, Bussery-Honvault et al. (20) were the first to develop a global ab initio PES [called Bussery-Honvault–Honvault-Launay (BHL)] for the aA1 state; the ab initio data were subsequently refitted with the reproducing kernel Hilbert space method (RKHS PES) (21). Later, Joseph and Varandas (22) developed a double many-body expansion aA1 PES. These surfaces have since been used in various dynamical calculations (1725). More recently, a highly accurate global ab initio PES (or ZMB-a) has been constructed by Bian’s group (26), which is unique in the accurate description of the regions around CIs and van der Waals interactions, and further dynamical calculations (27, 28) have been performed on this PES. For the bA1 state, a global ab initio PES [called Bussery-Honvault–Julien–Honvault-Launay (BJHL)] has been constructed by Bussery-Honvault et al. (29). This surface has been applied in the most recent studies (19), although the inclusion of contributions from the BJHL surface worsens the agreement between theory and experiment (17), indicating its probable deficiency. In contrast, the use of a bA1 PES newly constructed by Bian’s group has notably improved the agreement with crossed beam experiments (30).


Points of the C2v intersection between the bA1 and dA1 states are continuously connected, thus forming a seam line. The minimum energy crossing point (MECP) is determined by our ab initio calculations, with its geometry and energy indicated in Fig. 1A, and evidence that this intersection is conical is presented in Fig. 1B, fig. S1, and Supplementary Text. Figure 1 shows that the CI between the bA1 and dA1 states induces a barrier-like structure, which regulates the dissociation of a highly bent precursor intermediate supported by a shallow well (−14.37 kcal/mol relative to the reagent asymptote). The corresponding results of geometry optimization and frequency analysis are given in table S1. It is still unclear how this type of CI-R intermediate influences the reaction dynamics, and we will reveal its important role below.

The bA1 PES (hereafter called ZMB-b) constructed by Bian’s group provides an accurate description of the regions around the C2v CI (see fig. S2 and Supplementary Text), and an adjusted bA1 PES (hereafter called ZMB-bt) is also constructed for comparison (see below and fig. S2D). To fully characterize the role of the CI-R intermediate, we have performed extensive dynamical calculations on the ZMB-a, ZMB-b, and ZMB-bt PESs for the C(1D) + HD reaction. The C(1D) + HD isotopomer of the reactive system is chosen, since the CD/CH branching ratio is a sensitive probe of the PES topology and dynamics. These calculations are supported by our measurements (fig. S3) of the temperature-dependent rate coefficients and product branching ratios for both channels of the title reaction using a supersonic flow reactor over the temperature range of 50 to 296 K. Descriptions of the experimental and theoretical methodologies are provided in Materials and Methods and in the Supplementary Materials.

Figure 2A shows the experimental and theoretical rate coefficients as a function of temperature, where the exact time-dependent wave packet method including contributions from all angular momenta J (fig. S4) is used to perform QD calculations on the ZMB-a and ZMB-b PESs for the C(1D) + HD (v = 0, j = 0) reaction. Calculated initial state-specified rate coefficients are in excellent general agreement with the measured values, and both results present a slightly negative temperature dependence in the temperature range of 100 to 300 K. Although the calculated state-specified rate coefficient is slightly higher than the experimental one at room temperature, this can be attributed to the effects of higher HD initial rotational states, which will be discussed below.

Fig. 2 Comparison of experimental and theoretical results.

(A) The rate coefficients and (B) CD/CH product branching ratios as a function of temperature for the C(1D) + HD reaction. Accurate QD calculations are performed on the ZMB-a (aA1) and ZMB-b (bA1) PESs. The QD result on the BHL and BJHL surfaces is from (25), whereas the experimental result is from Fisher et al. (14).

Figure 2B shows the measured CD/CH product branching ratios alongside the calculated results on different PESs. The branching ratios derived by our QD calculations on the ZMB-a and ZMB-b surfaces are in excellent agreement with the experimental data at 50, 75, and 127 K, although a larger difference is observed at 298 K; a fact that could be due to the participation of higher HD(j) states at temperatures above 150 K. As the present QD calculations are too computationally expensive to take into account these higher HD rotational states, we performed quasi-classical trajectory (QCT) (31) calculations instead to elucidate the effects of j excitation on the calculated rate coefficients and branching ratios at 298 K. The QCT results show that the inclusion of higher HD rotational states lowers the branching ratio by 0.12, thereby reducing the difference between theory and experiment at this temperature. In addition, earlier experiments by Fisher et al. (14) yield a CD/CH ratio of 1.7 at room temperature, which is somewhat closer to our calculated value as shown in Fig. 2B. The QCT calculations also indicate that the calculated thermal rate coefficient at 298 K is slightly smaller when j excitation is included, bringing the theoretical results more into line with the experimental ones. The results from our QCT and QD calculations are compared in table S2, indicating that the two kinds of calculations are generally consistent with each other.

We also see from Fig. 2B that, at 127 K, the branching ratio from calculations on the BHL and BJHL PESs is close to that on the single ZMB-a surface and is lower than the measured value, whereas the agreement with the experiment achieved on the ZMB-a and ZMB-b surfaces is excellent below room temperature. Reaction over the ZMB-b surface yields very large CD/CH branching ratios (also see Table 1), a phenomenon that can be ascribed to the influence of the CI-R intermediate (Fig. 1A). Evidence is provided below describing the role of this intermediate species.

Table 1 Rate coefficients and CD/CH product branching ratios.

QCT and QD results on the ZMB-b and ZMB-bt surfaces for the C(1D) + HD (v = 0, j = 0) reaction are presented.

View this table:

To demonstrate the dynamical effect of the shallow well on the ZMB-b surface, we have constructed the ZMB-bt surface, which is identical to the ZMB-b one, except that the cyclic structure shallow well is removed. So the ZMB-bt surface omits the precursor well that supports the CI-R intermediate. We perform QCT and QD calculations on the ZMB-bt surface, and the rate coefficients and branching ratios obtained are summarized in Table 1, along with those obtained on the ZMB-b surface. The CD/CH branching ratios obtained from the ZMB-b surface are much larger than those from the ZMB-bt surface, although the rate coefficients obtained from both surfaces are very similar. It is clear that the CI-R shallow well plays the key role. Moreover, the CD/CH branching ratios from ZMB-b are also much larger than those from the BJHL surface, and this may be caused by the inadequate description of the CI-R barrier and shallow well on the BJHL surface (see the Supplementary Materials).

Table 1 also shows that the present QCT and QD results on the ZMB-b and ZMB-bt surfaces accord well, although the QCT method is more sensitive to the existence of the CI-R shallow well, leading to a more pronounced increase in the branching ratio. This justifies the following analysis with the QCT method, which is used to further account for the much larger branching ratios on ZMB-b than on other PESs below room temperature. Reactive classical trajectories on ZMB-b are compared with trajectories from the same initial conditions on other PESs. We find that the CI-R intermediate leads to unusual reaction mechanisms with obvious nonstatistical features.

We distinguish the newly found mechanisms into two types: Type I is designated as the C─H activation complex conversion mechanism, whereas type II is the cyclic complex mechanism (see the Supplementary Materials for details). At 0.01 eV collision energy, the type I mechanism dominates, accounting for 83% of all reactive trajectories. The main feature of the type I mechanism is that the C─H bond is selectively stretched or activated when the system leaves the shallow well region (before entering the deep well region). We attribute this selective activation of the C─H bond to two factors. First, the nascent cyclic complex has a large-amplitude C─H stretch vibrational mode (fig. S5A and table S1). Second, as shown in Fig. 1B, the top of the CI-induced barrier has a C2v symmetry, allowing the interaction potential to divert the system from the isosceles geometries. This effect helps to orient the C atom closer to the D atom before the system enters the deep well, promoted by the bending mode of the cyclic complex (fig. S5B).

The type I mechanism can be further divided into two submechanisms, C─H preactivation and reactivation ones, and their difference is that reactivation trajectories undergo at least one C─H bond reactivation by traveling back to the shallow well from the deep well before dissociation. The ratio between the former and latter mechanisms is about one-half at 0.01 eV. To illustrate the type I mechanism, snapshots of a typical reactive trajectory undergoing a simple C─H preactivation complex conversion are shown in Fig. 3. When the reagent C atom approaches the HD molecule, the reagents initially form an intermediate species bound by the shallow well, followed by a rapid increase in the H─D distance. Since the lighter H atom moves much faster than the heavier D atom, the newly formed C─H bond is notably longer than the C─D bond, leading to the formation of one highly activated (excited) bond (C─H) and one inactivated bond (C─D) before the trajectory enters the deep well region. Here, C─H bond activation is retained as the intermediate bends and inverts before dissociating into CD + H. Note that Fig. 3 shows a simple direct mechanism, accounting for only 8% of the C─H preactivation complex conversion mechanism. While many other reactive trajectories (92%) remain longer in the deep well before dissociating (fig. S6A), our analysis indicates that C─H activation is retained for more than 40% of these indirect trajectories. The C─H preactivation mechanism is beneficial for the formation of CD + H products and a correspondingly large CD/CH branching ratio. Actually, in the type I mechanism, the C─H reactivation submechanism dominates, and it is common for long-lived trajectories on ZMB-b to travel back and forth between the shallow and deep well regions (fig. S6B). Approximately half of the C─H reactivation trajectories preserve C─H bond activation before dissociation, which promotes CD + H formation.

Fig. 3 Revealing the C─H bond activation process.

(A) Snapshots of a typical reactive trajectory and (B) the corresponding internuclear distance and potential energy variations as a function of time. The trajectory is undergoing the type I mechanism and is calculated at 0.01 eV collision energy for the C(1D) + HD (v = 0, j = 0) reaction on the ZMB-b surface.

The type II mechanism also favors CD + H formation. Its dominant feature is that reactive trajectories experience the shallow well without entering the deep well region during the course of reaction. The type II mechanism can be further divided into two submechanisms: direct and indirect ones (see the Supplementary Materials). Snapshots of a typical reactive trajectory, which belongs to the cyclic complex direct mechanism, are shown in fig. S6D. Here, the C─H bond is stretched and activated in the cyclic complex configuration region, yielding CD + H products. At 0.01 eV, about 40% of the reactive trajectories follow this direct mechanism, with the indirect mechanism (fig. S6C) that also preferentially forms CD + H accounting for the rest. Reactive trajectories occurring through the type II mechanism lead to a large CD/CH branching ratio of 3.92. In summary, by examining all reactive trajectories at low collision energies, we find that the dynamical effect of the CI-R intermediate elaborated above is quite general, resulting in very large CD/CH branching ratios on the ZMB-b surface. At higher temperatures, trajectories with higher HD rotational excitation are involved, which are less strongly influenced by the CI-induced well. Consequently, this preference is weakened for a thermal population of initial HD(j) states at temperatures in excess of 150 K.

Although there are five electronic states correlating to the reagents C(1D) + HD, our calculations and analysis indicate that the adiabatic dynamics on the lowest two (aA1 and bA1) is dominant at low collision energies. The three higher excited states can only participate in the reaction via nonadiabatic pathways through CIs and have been neglected in our dynamics calculations. Two types of mechanisms contribute to reaction on the ZMB-b surface. These mechanisms are influenced by the CI-R intermediate, which has large-amplitude C─H stretch and bending vibrational modes (fig. S5 and table S1), diverting the system from the isosceles geometries required for the C2v bA1-dA1 CI. Consequently, the reaction system should quickly pass through the C2v CI regions with negligible nonadiabatic transition probabilities. In addition, the MECP on the linear CI seam between the aA1 and cA1 (or bA1 and dA1) states is located 12.4 kcal/mol above the reagent asymptote, below which the cA1, dA1, and eA1 states are not accessible at linear geometries. Furthermore, recent work (30) indicates that the inclusion of Renner-Teller effects is not expected to have a large influence on the results.

Large-scale QD calculations on the ZMB-a and ZMB-b surfaces are also performed to predict detailed scattering properties encoded in differential cross sections (DCSs). Some nonstatistical behavior is observed in the calculated DCSs. Typical three-dimensional product angular distribution plots are presented in fig. S7. More sideways distributions including small peaks are observed, which may be caused by the types I and II mechanisms mentioned above. Meanwhile, interesting oscillation patterns are revealed in fig. S7 (C and D). These arise from quantum interferences between partial waves according to our analysis. This phenomenon has recently been observed experimentally in the H + D2 (32) and H + HD (33) reactions. The oscillation phenomenon in DCSs predicted in fig. S7 is prominent enough to be observed experimentally. However, the types I and II mechanisms do not cause clear bias in the forward or backward directions, and approximately forward-backward symmetric product angular distributions are still obtained.

We also perform extensive ab initio calculations on other similar systems, finding a series of CI-R intermediates. The results of some of our geometry optimization and frequency analyses are given in table S1. In particular, a cyclic CS2 intermediate and the associated CI have been located, with the calculated vibrational frequencies of this cyclic species in excellent agreement with the experimental measurements of Bahou and coworkers (34). In addition, the calculated vibrational frequencies of the CI-R intermediate for PbO2 are also in very good agreement with the experiment, and its corresponding CI is shown to exist in a similar way to the CHD case.

Here, we introduce a general definition of the CI-R intermediate in a bimolecular complex-forming reaction as a special metastable-state intermediate whose conversion to a deep well complex is regulated by a CI-induced barrier. In the title reaction, we find that it has the effect of bond-selective activation, markedly increasing the CD/CH product branching ratios, which are validated by the excellent agreement between theory and experiment. We reveal that the CI-R intermediate leads to two types of reaction mechanisms, having clear nonstatistical features. We have also shown that the CI-R intermediate exists in other systems, and its dynamical effects may have general significance. It is our hope that the CI-R intermediate may become a useful concept in the understanding of unusual phenomena in various complex-forming reactions and in the laser control of chemical reactions.


Experimental measurement

The kinetics of C(1D) + HD reaction were studied experimentally using a continuous supersonic flow reactor over the temperature range of 50 to 296 K using three Laval nozzles based on argon as the carrier gas. Excess concentrations of HD were introduced into the main flow upstream of the Laval nozzle along with CBr4 vapor, which was carried into the reactor by passing a small argon flow over solid CBr4. C(1D) atoms were generated by the pulsed multiphoton dissociation of CBr4 at 266 nm. Previous work (35) has shown that, while ground electronic state C(3P) atoms are the major photolysis product, C(1D) atoms are also created at the 10 to 15% level with respect to C(3P). The C(3P) + HD → CD(CH) + H(D) reaction itself is endothermic by 22.8 kcal/mol, so that this process is negligible in the present work.

The title reaction was investigated by following D atom production. D(2S) atoms were detected through on-resonance vacuum ultraviolet (VUV) laser-induced fluorescence (LIF) using the 1s 2S → 2p 2P0 transition at 121.534 nm. VUV LIF signals from product D(2S) atoms were recorded as a function of time between photolysis and probe lasers. The evolution of the D atom signal (ID) is described by the expressionID=A{exp(kL(D)t)exp(kt)}(1)k′ = k[HD] + kL(C), where k is the second-order rate coefficient for the C(1D) + HD reaction, [HD] is the HD concentration, and kL(C) is the first-order loss of C(1D) by other processes. kL(D) is the first-order rate coefficient for D(2S) loss, which is dominated by diffusion.

Two types of experiments were performed. The first type of experiment allowed us to determine the temperature-dependent rate coefficients. Here, D(2S) kinetic profiles, such as those shown in fig. S3A, were recorded at different HD concentrations, yielding a range of k′ values. Plots of the individual k′ values versus [HD] allowed the second-order rate coefficient to be determined at a given temperature through linear least squares fits to the data. Typical second-order plots for data recorded at 296 and 50 K are shown in fig. S3B.

The second type of experiment allowed us to derive absolute yields for D(2S) atom production as a function of temperature. Here, the time-dependent relative D atom signal intensity produced by the C(1D) + HD reaction was compared to the one produced by a reference process, namely, the C(1D) + D2 reaction, which is considered to have a 100% D atom yield. A pair of these D atom formation curves recorded at 50 K are shown in fig. S3C. More details of experimental measurements can be found in table S3 and Supplementary Methods.

Ab initio calculations and PESs

Full geometry optimizations were run to locate the MECP of the C2v CI between the bA1 and dA1 states and the potential minimum corresponding to the CI-R intermediate for the C(1D) + HD system through ab initio calculations. The internally contracted multireference configuration interaction (36) and five state–averaged complete active space self-consistent field methods (37) were used. Five reference states were used to generate the internally contracted pairs. The active space consists of six electrons distributed among seven orbitals. The basis set used was the Dunning’s correlation consistent quadruple-zeta basis augmented by diffuse functions (aug-cc-pVQZ) (38). The corresponding ab initio calculations for other similar systems (CS2, SiO2, etc.) were performed at lower levels. All ab initio calculations were performed with the MOLPRO package (39).

The two highly accurate ab initio global PESs for the aA1 and bA1 states of the title system constructed by Bian’s group were applied; in addition, an adjusted global PES for the bA1 state, which omits the CI-induced barrier and well, was constructed. The aA1 PES, called ZMB-a, was constructed in our previous work (26), in which the main scheme of PES construction was described. The bA1 PES constructed by Bian’s group was applied in our recent dynamical studies (30) and, for convenience, is referred to as the ZMB-b surface in this work, about which further analysis and more details are given in Supplementary Text and fig. S2. The adjusted bA1 PES is called the ZMB-bt surface, which is the same as the ZMB-b surface in various regions, except that it does not have the C2v shallow well (one typical contour plot is given in fig. S2D). It is clear that the ZMB-bt surface leaves out the precursor well that supports the CI-R intermediate complex.

QD and QCT calculations

The QD calculations including all Coriolis couplings were performed using the real wave packet approach as implemented in the DIFFREALWAVE code (40). Numerical details are given in Supplementary Methods, and numerous test calculations were performed to check the convergence. Partial wave contributions up to J = 50 and J = 45 were included in the calculations for the aA1 and bA1 states, yielding converged integral cross sections and DCSs for collision energies up to 0.40 and 0.45 eV, respectively. About 120,000 iteration steps were implemented for each partial wave calculation, which were reduced at high J values.

The QCT calculations were performed using a modified version of the VENUS96 code (41) customized to incorporate the PESs. The QCT method is described in the literature (31), and some computational details are given here. All trajectories were initiated with the reagents separated by 15 bohr. An integration time step of 0.01 fs was chosen, which guarantees an energy conservation of better than 1 in 106. Batches of 10,000 trajectories were calculated for each set of initial conditions with the reagent HD set in different rovibrational states, including those calculated at the temperatures of 50, 75, 127, and 298 K, with the impact parameter randomly sampled between 0 and 14 bohr and with the relative translational energy randomly selected to mimic the Boltzmann distribution at the given temperature.


Supplementary material for this article is available at

Supplementary Methods

Supplementary Text

Fig. S1. Sign changes of the wave function main configuration coefficient for the bA1 state.

Fig. S2. Contour plots of the ZMB-b and ZMB-bt surfaces.

Fig. S3. Experimental LIF profiles and second-order rate coefficient plots.

Fig. S4. Total reaction probabilities for the C(1D) + HD (v = 0, j = 0) reaction.

Fig. S5. Diagrams of the vibrational modes of the CI-R intermediate in the C(1D) + HD reactive system.

Fig. S6. Internuclear distances and potential energy as a function of time for typical trajectories.

Fig. S7. Three-dimensional product flux plots in the center of mass velocity space for the C(1D) + HD (v = 0, j = 0) reaction.

Table S1. The properties of CI-R intermediates and the CIs associated with them.

Table S2. Comparison between QD and QCT results on the ZMB-a and ZMB-b surfaces.

Table S3. Measured rate coefficients and CD/CH branching ratios for the C(1D) + HD reaction.

References (4251)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was supported by the National Natural Science Foundation of China (nos. 21773251, 21473218, and 91741106). J.C. acknowledges the Youth Innovation Promotion Association CAS (no. 2018045). K.M.H. and D.N.-R. acknowledge support from the French program “Physique et Chimie du Milieu Interstellaire” (PCMI) of the CNRS/INSU with the INC/INP cofunded by the CEA and CNES and funding from the “Program National de Planetologie” (PNP) of the CNRS/INSU. W.B. and J.C. acknowledge Beijing National Laboratory for Molecular Sciences and the CAS Key Research Project of Frontier Science. Numerous computations were carried out at National Supercomputing Centers in Shenzhen and Tianjin. Author contributions: W.B. and K.M.H. conceived and supervised the research. K.M.H. and D.N.-R. performed the experiment and analyzed the data. Y.W. and J.C. carried out the QD and QCT calculations and analyzed the data. H.M., C.Z., and Y.W. performed the ab initio calculations. C.Z., H.M., J.C., and W.B. did the work associated with the PESs. W.B., K.M.H., and J.C. wrote the paper with contributions from all co-authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article