X-ray liquidography decodes complex motions in azobenzene isomerization

Data collection

TRXL experiments were conducted at the XSS beamline of the Pohang Accelerator Laboratory X-ray free-electron laser (PAL-XFEL) facility. Key experimental parameters are summarized in Supplementary Table 1. A laser-pump and X-ray-probe scheme was used to initiate the photoisomerization reaction of azobenzene. Optical laser pulses at a wavelength of 400 nm, generated as the second harmonic of a femtosecond titanium:sapphire laser system, were used as the pump pulses. These pump pulses with an energy of 8.5 µJ were focused onto the sample with a spot size of 55 × 66 µm² (full-width at half-maximum (FWHM), horizontal × vertical), corresponding to a laser fluence of 1.5 mJ mm⁻². The repetition rate of the optical laser was 30 Hz. The X-ray beam, used as the probe, was focused onto the sample with a spot size of 20 µm × 20 µm (FWHM), delivering approximately 1 × 10¹² photons per pulse at a repetition rate of 60 Hz, with a centre energy of 15.01 keV. By utilizing the X-ray with twice the repetition rate of the laser, we alternately measured laser-on and laser-off scattering images. The optical laser and X-ray beams were overlapped at a 10° crossing angle at the sample position. The temporal resolution of the experiment, defined by the instrument response function (IRF), was determined by the temporal widths of optical laser and X-ray pulses, the velocity mismatch due to the capillary jet width and timing jitter. On the basis of our data analysis, the IRF value was found to be 169 ± 5 fs (FWHM; Supplementary Method 1).

The sample solution, 50-mM trans-AB (Sigma-Aldrich, 98%) dissolved in methanol, was flowed through a 100-µm capillary using a gear pump to create an open jet for the pump–probe measurement. The jet speed was approximately 3 m s⁻¹, sufficient to prevent multiple exposures of the sample to the pump laser pulses. Unlike typical TRXL experiments where the sample solution is recirculated, the sample solution was flowed only once and discarded after a single exposure to the laser and X-ray pulses, as the isomerization reaction to cis-AB is irreversible. In many solvents, cis-AB exhibits two- to three-times higher absorption than trans-AB at 400 nm (ref. ²¹). Thus, the generation of cis-AB under ambient light could potentially interfere with obtaining an unaltered signal from trans-AB. To address this, the preparation of trans-AB was performed in a darkroom, and the TRXL experiments were also conducted in the absence of ambient light. Furthermore, the absorption spectrum of trans-AB under ambient light showed no detectable changes within 5 h. For the trans-AB TRXL experiments, we used trans-AB solutions prepared within 5 h.

To obtain reference data representing the response of bulk solvent, a separate fs-TRXL measurement was conducted with neat methanol at a pump fluence of 5.3 mJ mm⁻². The neat methanol was heated via two-photon absorption of the 400-nm pump pulse, depositing energy equivalent to single-photon excitation at 200 nm. All chemicals, including the solvent, were sourced from commercial suppliers and used as received without further purification.

The X-ray scattering patterns were measured using a Rayonix MX225-HS detector operated in 6 × 6 binning mode, resulting in an effective pixel size of 234 µm and a frame size of 960 × 960 pixels. The sample-to-detector distance was set to 40.4 mm, providing a q range of 1.0 Å⁻¹ to 7.0 Å⁻¹. For the trans-AB solution, TRXL data were collected at 140 nominal time delays from −0.5 ps to 100 ps (uniform 0.05 ps steps to 2.0 ps, followed by logarithmic spacing of approximately 50 points per decade), with approximately 250 laser-on and 250 laser-off single-pulse images per delay per run. To achieve a sufficient signal-to-noise ratio (SNR) for data analysis, we acquired 18 runs for the target sample, trans-AB, and 4 runs for the neat solvent. In total, around 4,500 laser-on scattering images per time delay (with an equal number of laser-off images) were accumulated, resulting in fs-TRXL data with satisfactory SNRs.

The molar absorptivity of trans-AB at 400 nm was measured to be 292 M⁻¹ cm⁻¹ under the experimental conditions. To confirm that the photoisomerization is initiated by a single-photon process at our experimental fluence (1.5 mJ mm⁻²), we examined the pump-fluence dependence of the signal using femtosecond transient absorption experiments in which the pump fluence was varied up to 1.7 mJ mm⁻². The signal amplitude at 2 ps exhibited a linear relationship with pump fluence (not shown). This linearity confirms that multi-photon excitation or ionization processes are negligible under our experimental conditions.

Data processing

Supplementary Figs. 1 and 2 detail the entire data-processing workflow, from raw images to the final structural refinement, explicitly showing each separation and correction step. The spatial overlap between the pump and probe pulses was established by adjusting the optical laser pointing to maximize the transient signal amplitude. The temporal overlap between the two pulses was ensured before data collection using an yttrium aluminium garnet (YAG) crystal, by monitoring the X-ray-induced change in its optical transmittance to identify the coincidence of the two pulses. No dedicated timing tool was used in this experiment. This YAG measurement was used to determine only the nominal time zero of the recorded delays; the precise time zero was refined afterwards, during data processing, in the two stages described later. Before all subsequent processing steps, each recorded diffraction image was corrected by subtracting a dark image (background) acquired under X-ray beam-off conditions, followed by corrections for polarization and solid angle effects. The polarization correction accounts for the linearly polarized nature of XFEL X-rays, which causes azimuthally anisotropic scattering intensity not related to sample anisotropy. Specifically, the measured intensity (I_measured) was corrected for polarization as I_polar = I_measured × 1/[1 − sin²(2θ)cos²(φ)], where φ is the azimuthal angle relative to the X-ray polarization direction and 2θ is the scattering angle. The solid angle correction compensates for the reduced solid angle subtended by detector pixels at higher scattering angles in a flat detector geometry. This compensation requires an intensity correction with the formula: I_corrected = I_polar × 1/cos³(2θ). These corrections were applied using standard formulas implemented in the pyFAI package⁴³.

Each corrected two-dimensional scattering image was then processed to obtain one-dimensional scattering curves as a function of q, the magnitude of the momentum transfer vector (q = 4πsinθ/λ, where λ is the X-ray wavelength). For each single-shot image, we obtained either (1) azimuthal scattering curves, S_azi(q, t)^single, through azimuthal averaging using the ‘integrate1d’ function from pyFAI, or (2) separated isotropic and anisotropic scattering curves, S_iso(q, t)^single and S_aniso(q, t)^single, via the anisotropic decomposition method^{36,40,44,45,46,47} using the ‘integrate2d’ function from pyFAI. The isotropic component contains information on changes in interatomic distances independent of molecular orientation, whereas the anisotropic component reflects orientational dynamics of the molecules. The detailed procedure for the anisotropic decomposition method is provided in the literature^{36,40,44,45,46,47}. This two-dimensional-to-one-dimensional conversion was applied identically to every detector image, regardless of whether it was recorded laser-on or laser-off; each image yields one isotropic and one anisotropic curve, and the laser-on/off distinction enters only at the subsequent difference and pairing step (below).

Before calculating difference scattering curves, each scattering curve was normalized by its area under the curve (AUC). Specifically, S_azi(q, t)^single and S_iso(q, t)^single were normalized by dividing by the AUC within the normalization range of q = 4–6 Å⁻¹, while S_aniso(q, t)^single was normalized by the AUC of the corresponding S_iso(q, t). For each laser-on (pumped) measurement, we then calculated difference scattering curves by subtracting an appropriately paired laser-off (unpumped) signal. Rather than using simple time-sequential pairing, we used an AUC-based pairing method to minimize systematic errors. Although normalization is applied to account for shot-to-shot fluctuations in beam intensity and sample thickness, detector nonlinearity (that is, imperfect linearity between incident photon count and recorded signal intensity) can introduce systematic artefacts when subtracting laser-on and laser-off images with substantially different signal amplitudes. To suppress these artefacts, each laser-on curve was subtracted by its matched laser-off curve, the one with the most similar AUC value, thereby minimizing the influence of detector response nonlinearity on the resulting difference signal. This procedure yielded single-shot difference scattering signals, ΔS_[x](q, t)^single ([x] = azi, iso or aniso), for each time delay t.

To improve the SNR, multiple single-shot difference curves measured at the same time delay were averaged in two stages: (1) within each run and (2) across runs. Before averaging, an outlier rejection procedure was applied to ensure data quality. For each single-shot measurement, we calculated the AUC of the absolute difference signal (∫|ΔS(q, t)^single|dq) over the q range of 1.0–7.0 Å⁻¹. For the azimuthal data, the AUC was computed from ΔS_azi(q, t)^single. For the isotropic and anisotropic data, where each shot yields a paired isotropic and anisotropic curve, the AUC was computed from the isotropic difference, ΔS_iso(q, t)^single; shots identified as outliers were then excluded from both the isotropic and their paired anisotropic averages, with the anisotropic curves not filtered separately. Single-shot curves with AUC values exceeding an empirically determined threshold were identified as outliers and excluded from subsequent averaging. This procedure typically retained 95–99% of the collected data, removing only extreme outliers arising from beam instability or sample flow irregularities. The outlier-filtered difference curves were first averaged within each run to obtain the run-averaged ΔS_[x](q, t)^run and its associated standard error of the mean, σ_[x](q, t)^run. The standard error of the mean was calculated as:

Here σ_[x](q, t)^single is the standard deviation of the single-shot difference curves within that run for the specific time delay and q value, and N_pair is the number of difference scattering curves within that run that survived the outlier rejection.

The run-averaged curves from the 18 runs of the trans-AB TRXL data, all acquired on the same set of 140 nominal time delays (−0.5 ps to 100 ps), were then averaged across runs into a single dataset (across-run-averaged). First, the nominal time delays were referenced to the experimental time zero by subtracting a single, run-independent preliminary offset, t_0,pre, which was determined from the across-run-averaged anisotropic signal before time-pooling, giving the preliminary pump–probe delays. The 18 runs differed slightly in the overall amplitude of their difference signals, reflecting small run-to-run variations in the experimental conditions, such as the sample (liquid jet) thickness and the position and intensity of the X-ray and laser pulses. To place them on a common amplitude scale, each run was first characterized by a single amplitude: the integrated magnitude of its isotropic difference signal, ∫|qΔS_iso|dq over q = 1–7 Å⁻¹ averaged over the 90–100 ps time delays. The run with the largest amplitude was taken as the reference. Each run-averaged curve was then multiplied by a scaling factor equal to the reference amplitude divided by that run’s own amplitude, so that run-to-run amplitude fluctuations were removed. The corresponding uncertainty, σ_[x](q, t)^run, was scaled by the same factor. The t_0,pre-corrected, scaled curves were then binned onto a common analysis grid of 137 time delays (spanning −0.35 ps to 100 ps: uniform 0.05 ps spacing up to 2.0 ps, then approximately logarithmic spacing of approximately 50 points per decade). Specifically, each run’s t_0,pre-corrected delay was assigned to the nearest delay of the analysis time grid, and the contributions from all runs at each grid delay were then averaged with weights proportional to the N_pair of the run-averaged curve (an image-count-weighted average). This across-run averaging onto the analysis time grid is distinct from the pooling along the time axis applied later (below), and gave the across-run-averaged difference scattering curves ΔS_iso(q, t)^raw, ΔS_aniso(q, t)^raw and ΔS_azi(q, t)^raw (Extended Data Fig. 1a,c,e); their uncertainties σ_[x](q, t)^raw were propagated from the σ_[x](q, t)^run according to standard error propagation rules.

For quantitative analysis, these across-run-averaged curves were placed on an absolute scale, normalized to the scattering intensity of a single solvent molecule, following standard methods in the TRXL community^{36,38,40,45,46,48,49,50,51,52,53,54,55,56,57,58,59,60,61}. The total scattering curve of the solvent was calculated from the pair distribution functions g(r) obtained from molecular dynamics simulations, giving the scattering intensity in electron units (e.u.), and normalized by the number of solvent molecules in the molecular dynamics simulation box to yield the intensity in e.u. per solvent molecule. A scaling factor was determined by comparing the normalized experimental total scattering curve with this simulated curve. This use of the simulated scattering curve of the pure solvent to calibrate the experimental solution data is justified because, at the low solute concentration (50 mM), the experimental total scattering intensity of the solution arises almost exclusively from the bulk solvent. This single scaling factor was applied identically to all three experimental difference scattering curves, ΔS_iso(q, t)^raw, ΔS_aniso(q, t)^raw and ΔS_azi(q, t)^raw, and to their uncertainties σ_[x](q, t)^raw, expressing the signals in e.u. per solvent molecule. This normalization process enables the extraction of quantitative information regarding the reaction dynamics, including the concentration of solute molecules participating in the isomerization pathways.

The PEPC step then removes the solvent contribution from the signal, propagating the raw-data uncertainty σ_[x](q, t)^raw into the solvent-independent data, ΔS_[x]^⊥(q, t_k), whose per-point standard error is σ_[x]^⊥(q, t_k). For the subsequent quantitative analysis, these solvent-independent data at adjacent time delays were pooled along the time axis to further improve the SNR. Rather than an unweighted two-point moving average, the difference scattering curves at each pair of neighbouring time delays (t_k and t_k+1) were averaged with weights proportional to the number of images contributing to each delay (an image-count-weighted mean of the signals). Each pooled time delay was assigned the corresponding image-count-weighted mean of the two time delays. This time-pooling reduced the 137 original time delays on the analysis grid (spanning from −0.35 ps to 100 ps) to 136 pooled time delays spanning about −0.33 ps to 98 ps. The same image-count time-pooling was also applied to the pre-PEPC anisotropic data, ΔS_aniso(q, t)^raw, which was then used to determine the refined time zero and the IRF, as described below. The isotropic curves in Figs. 2 and 3 are the time-pooled dataset, whereas the curves in Extended Data Figs. 1 and 2 are shown before time-pooling. Throughout this paper, difference curves are presented as qΔS(q, t) for visualization in the figures, to emphasize high-q features associated with solute structural changes. This q weighting was applied only for display and not used in quantitative analysis.

The experimental time zero was accounted for in two stages. The preliminary offset above (t_0,pre) served only to reference the nominal delays to time zero before merging. The refined time zero (t_0,ref) and the IRF were then determined, after merging and time-pooling, by fitting the time-pooled, pre-PEPC anisotropic signal. Near time zero, the solvent produces a coherent response, the optical Kerr effect, a sharp feature that marks the arrival of the pump pulse and therefore provides an accurate temporal reference for the time zero and the IRF. Because PEPC removes the solvent contribution, it would also remove this feature. The fit therefore uses the anisotropic data before PEPC, where the feature is preserved. The refined time zero is consistent with t_0,pre to within about 15 fs (well below the IRF). These t_0,ref and IRF were then held fixed as parameters in the kinetic analysis of the time-pooled, PEPC-treated isotropic (Fig. 3) and anisotropic (Extended Data Fig. 3) signals.

From these isotropic, time-pooled, solvent-independent data, ΔS_iso^⊥(q, t_k), the kinetic analysis involving SVD and kinetics-constrained analysis was performed first, followed by structural refinement. Fitting the RSVs obtained from SVD yielded time constants used to set up the kinetic model. The time-dependent concentration profiles of the species were calculated from the kinetic model and then used by the kinetics-constrained analysis to extract the SADSs for the structural refinement.

The standard error of the RSVs enters the analysis at two stages: (1) the fit of the RSVs and (2) the subsequent estimation of the uncertainties of the resulting time constants. This standard error is denoted σ_n(t_k) for the nth RSV at time delay t_k. The two stages depend on σ_n(t_k) in different ways. The fit is sensitive only to the relative magnitudes of σ_n(t_k) across the time points, as these relative magnitudes determine the relative weights of the data points and hence the location of the χ² minimum, that is, the best-fit time constants. A uniform rescaling of σ_n(t_k) multiplies the entire χ² by a constant and therefore leaves the location of its minimum, and hence the fitted time constants, unchanged. The uncertainties of the time constants, by contrast, are obtained by displacing each time constant from its best-fit value until the χ² increases by a fixed amount, and the required displacement is taken as the uncertainty. Such a displacement worsens the agreement between the RSVs and the model, so that the residuals, and with them the χ², increase. Because the χ² weights the residuals by 1/σ_n(t_k)², the rate of this increase is governed by the absolute magnitude of σ_n(t_k): a smaller σ_n(t_k) produces a faster increase, and hence a smaller displacement to the threshold. The resulting uncertainty therefore depends on the absolute magnitude of σ_n(t_k), and not merely on its relative magnitudes. We therefore used two estimates of σ_n(t_k): an analytic estimate σ_n^ana(t_k) for weighting the fit, and a bootstrap estimate σ_n^boot(t_k) for estimating the uncertainties.

The model fitted to each RSV is a sum of exponential functions convolved with the IRF, and each RSV value is weighted by 1/[σ_n^ana(t_k)]². The analytic estimate σ_n^ana(t_k) was obtained by propagating the per-point error of the data, σ_[x]^⊥(q, t_k), through the time-pooling and the SVD. It reproduces the correct relative magnitudes of σ_n(t_k), and therefore yields the correct time constants, but it underestimates the absolute magnitude of σ_n(t_k), because it does not fully account for the correlation introduced between neighbouring time points by the two-point time-pooling (Supplementary Method 1).

The bootstrap estimate σ_n^boot(t_k), by contrast, does capture this correlation. It was obtained by generating repeated synthetic datasets, each formed by adding random noise to the data at the level of the per-point error σ_[x]^⊥(q, t_k), and re-running the time-pooling and the SVD on each (Supplementary Method 1). Using σ_n^boot(t_k) as the standard error, the uncertainties of the time constants were then obtained by the profile-likelihood analysis, with each time constant displaced until the χ² increased by Δχ² = 1 and the remaining parameters re-optimized at each displacement (Supplementary Methods 1 and 4). The standard error of the SADSs, denoted σ_n(q_j) for the nth SADS at the jth q point q_j, was obtained by an analogous bootstrap (Supplementary Methods 2 and 3). In the structural refinement, the structural parameters were determined by minimizing a χ² in which the difference in the SADSs at each q_j was divided by its standard error σ_n(q_j) and squared. Their 1σ uncertainties were obtained by the same profile-likelihood criterion (Δχ² = 1). The σ_n(q_j) are shown as the error bars on the SADSs (Fig. 3f), and the profile-likelihood limits as the error bars on the refined structural parameters (Fig. 4a).

Obtaining solvent-independent TRXL data

TRXL signals generally consist of three main components: (1) the signal arising from changes in the solute itself (solute term, ΔS_solute(q)), (2) the signal resulting from changes in the solvent–solute interatomic distances (cage term, ΔS_cage(q)), and (3) the signal from the response of the bulk solvent (solvent term, ΔS_solvent(q)), which includes contributions from changes in the temperature and density of the bulk solvent induced by heat dissipation from the photoexcited solute molecules, as well as the rearrangement of solvent molecules due to the optical Kerr effect⁴⁶. The sum of the solute and cage terms is referred to as solute-related terms (ΔS_{solute-related}(q) = ΔS_solute(q) + ΔS_cage(q)), and the kinetics of ΔS_{solute-related}(q) reflect the behaviour of the solute molecules. In contrast, ΔS_solvent(q) may exhibit different kinetics associated with the hydrodynamic response of bulk solvent and optical Kerr effect of the bulk solvent⁴⁶. To accurately analyse the kinetics of the solute molecules, it is important to eliminate the influence of the solvent term from the TRXL signals.

As shown in Fig. 2, more than 90% of the TRXL signal of azobenzene is attributed to the solvent signal. This large contribution from the solvent makes it crucial to accurately account for and subtract this component to extract the solute-related kinetics. To address this, the PEPC was used⁴¹. The core concept of this approach is to project the data into a transformed space where the contribution of ΔS_solvent(q) becomes zero, thereby effectively screening out the influence of ΔS_solvent(q). We denote the data projected into this space using the superscript ⊥ symbol. In this transformed space, the solvent contribution is mathematically expressed as ΔS^⊥_solvent(q) = 0. As a result, the PEPC-treated data simplifies to ΔS^⊥(q) = ΔS^⊥_solute(q) + ΔS^⊥_cage(q) + ΔS^⊥_solvent(q) = ΔS^⊥_solute(q) + ΔS^⊥_cage(q) = ΔS^⊥_{solute-related}(q), leaving only the solute and cage contributions. By projecting the data into a space orthogonal to the solvent terms, this transformation isolates the solute-related terms, enabling the accurate extraction of solute-specific kinetics and structural information, as detailed in previous work^41,60,62. For the PEPC treatment, we used heating signals and artefact components obtained from a separate experiment with a neat solvent as known reference components, allowing for the removal of their contributions from the data. By analysing the resulting PEPC-treated data, we were able to isolate the kinetics of the solute molecules. The solvent contributions were effectively screened out from the original TRXL data (Extended Data Fig. 1) to obtain the PEPC-treated data presented in Extended Data Fig. 2, which we refer to as the solvent-independent data. The term ‘solvent-independent data’ is equivalent to the term ‘solvent-contribution-free data’ used in a previous study⁶²; both describe TRXL data that has been projected into a subspace orthogonal to the solvent response. To obtain the reference components representing the solvent response for this PEPC procedure, we performed a separate fs-TRXL measurement of neat methanol under two-photon heating conditions. The LSVs extracted directly from the SVD of this neat solvent data were used as the reference components without further modification. These components sufficiently capture the signal variations corresponding to the solvent response in the azobenzene solution, enabling effective removal of solvent contributions through orthogonal projection. Both ΔS_iso(q, t)^raw and ΔS_aniso(q, t)^raw from the neat solvent dataset were independently subjected to SVD. The reference components were selected by individually inspecting each LSV and RSV pair against two criteria: (1) autocorrelation coefficient of the RSVs ≥ 0.5, ensuring temporal consistency, and (2) singular values well above the noise plateau, ensuring the components represent genuine signal variations rather than noise. Through this evaluation, 7 LSVs (LSVs 1–6 and 9) were selected for ΔS_iso(q, t)^raw, and 4 LSVs (LSVs 1–4) were selected for ΔS_aniso(q, t)^raw (Supplementary Figs. 4 and 5). LSV₄ for ΔS_iso(q, t)^raw was included despite a lower overall autocorrelation because it exhibited a dominant contribution near time zero, capturing the initial ultrafast response. By separately projecting the ΔS_iso(q, t)^raw and ΔS_aniso(q, t)^raw of the trans-AB solution data onto the subspaces orthogonal to their respective reference components, PEPC effectively isolates the signal components strictly independent of the solvent response, thereby retaining only the signal components associated with solute structural changes that are orthogonal to the solvent response. The resulting solvent-independent data (ΔS^⊥(q, t)) are presented in Extended Data Fig. 2. During this PEPC procedure, the standard errors of the PEPC-treated data, σ^⟂(q, t), were calculated by propagating the standard error of the mean of the original TRXL data (σ(q, t)^raw) and the uncertainties associated with the reference solvent components. Specifically, as PEPC is a linear projection operation, the covariance matrix of the projected data was computed using standard error propagation rules to derive σ^⟂(q, t), ensuring that the uncertainties in both the raw experimental data and the solvent basis set were rigorously accounted for. We note that the oscillations observed in the low-q region of some SADS curves are a mathematical consequence of the PEPC procedure, which projects the data into a space orthogonal to the solvent response, and therefore should not be interpreted on their own as direct structural fingerprints of the solute. A flow chart illustrating the data analysis pipeline is shown in Supplementary Fig. 2.

Structure analysis

We refined the structures of I₁, I₂ and cis-AB by minimizing χ_ν², which quantifies the discrepancy between the experimental and simulated SADSs. As described in Supplementary Method 5, simulated SADSs, SADS_n′(q)s, were calculated based on the structures of four species: trans-AB, I₁, I₂ and cis-AB. For trans-AB, the structure derived from quantum chemical calculations (CASPT2) was used without further refinement. For I₁ and I₂, structural refinement was performed using a motion-based approach, considering four key structural motions as degrees of freedom: N–N rotation (ROT), inversion (INV), N–N bond-length change (NND) and C–N torsion (CNT). ROT describes the out-of-plane rotation of one phenyl ring about the N–N bond axis. This motion is quantified by the CNNC dihedral angle (∠CNNC), which defines the twisting of the molecule around the azo linkage. To model this, one phenyl group was rotated by a specified angle relative to the N–N axis, thereby introducing a controllable change in torsional strain. In INV, the CNN angles of both phenyl rings (∠CNN₁ and ∠CNN₂) are allowed to vary independently, with each ∠CNN treated as a separate fitting parameter during refinement. The concerted inversion of azobenzene, which has been extensively proposed in theoretical studies (Extended Data Table 2), can be described as a special case of INV, in which both ∠CNNs increase simultaneously and to the same extent, accompanied by a widening separation of the two rings. The NND motion is quantified simply by the N–N bond length (r_NN). The CNT motion describes changes in the ∠CCNN dihedral angles, parameterized by two ∠CCNN angles. This internal coordinate change can be realized either by torsion of the phenyl rings about their respective C–N bond axes or by concerted displacement of the azo group (bicycle-pedal motion); both scenarios produce identical ∠CCNN changes and are therefore indistinguishable in the internal coordinate description. Each motion acts exclusively on a distinct structural degree of freedom (∠CNNC, ∠CNN, r_NN and ∠CCNN, respectively) and preserves the geometric reference frames of the others. Although CNT has received relatively little attention in previous studies of azobenzene isomerization, it was included to capture ∠CCNN changes that may accompany the structural rearrangement. Detailed definitions of these motions are provided in Supplementary Method 14. The cis-AB structure was refined separately by a direct two-parameter adjustment starting from the CASPT2 reference (see below).

For I₁ and I₂, we used SSPM, in which all possible combinations of the four motions are systematically evaluated and ranked by the BIC to identify the most parsimonious model consistent with the data (Supplementary Method 12). The intermediate structures I₁ and I₂ were constructed by applying the four elementary structural motions (ROT, INV, NND and CNT) defined above to their respective precursors (CASPT2-calculated trans-AB for I₁, and the TRXL-refined I₁ for I₂). The structures of the intermediates generated by applying these four atomic motions were used to calculate the solute-only and cage scattering signals. The structural parameters corresponding to these four motions were optimized through χ_ν² minimization to ensure the best agreement between the simulated SADS′(q) and experimental SADS(q). The I₁ structure was generated by applying motions to the CASPT2-calculated trans-AB geometry, whereas the I₂ structure was generated by applying motions to the TRXL-refined I₁ geometry. For cis-AB, rather than performing an SSPM scan, a direct two-parameter refinement (∠CNNC and ∠CCNN) was carried out starting from the CASPT2-calculated cis-AB geometry with Tikhonov regularization towards the CASPT2 reference (Supplementary Method 12.7). In all cases, the structural parameters and scaling coefficients were optimized simultaneously through nonlinear least-squares minimization (Supplementary Methods 3 and 12).

The structural refinement was carried out iteratively to achieve self-consistency among the solute atomic charges, the molecular dynamics simulations and the TRXL-refined geometry. In the first iteration, natural population analysis (NPA) atomic charges were computed for the CASPT2-calculated structures and used to parameterize GROMACS molecular dynamics simulations, from which the cage terms were obtained and applied in the structural refinement. The NPA charges were computed at different levels of theory according to the electronic state of the species: B3LYP-D3BJ/def2-TZVP for the ground-state species (trans-AB and cis-AB), and state-specific CASSCF(10,8)/aug-cc-pVDZ corrected by strongly contracted NEVPT2 (SC-NEVPT2) for the S₁ intermediates (I₁ and I₂); in all cases, solvation was described using an implicit solvent model, the solvation model based on density. In each subsequent iteration, the NPA atomic charges were recomputed from the geometry refined in the preceding iteration, updated cage terms were generated and the structure was refined again. This cycle was repeated until the structural parameters showed negligible change between consecutive iterations, confirming convergence. Full details are provided in Supplementary Methods 8 and 9.

SADS_n′(q)s calculated from the refined structures of I₁, I₂ and cis-AB demonstrated excellent agreement with the corresponding experimental SADS_n(q)s (Fig. 3f). For cis-AB, the refinement confirmed that the solution-phase structure of cis-AB deviates only mildly from the CASPT2-calculated structure, showing a relatively small improvement in χ_ν² from 1.20 (CASPT2-calculated) to 0.97 (TRXL-refined), corresponding to a 19% reduction. This indicates both the sufficient accuracy of the CASPT2 structures of trans-AB and cis-AB, and the high quality of the TRXL data. In contrast, the refinement substantially improved the fits for I₁ and I₂, reducing the χ_ν² values from 0.46 to 0.24 for SADS₁ (47% reduction) and from 0.79 to 0.33 for SADS₂ (58% reduction). This indicates that the CASPT2-calculated structures of I₁ and I₂ deviate from their TRXL-refined counterparts to a greater extent than in the case of cis-AB. The structures of I₁, I₂ and cis-AB obtained from (1) CASPT2 calculations and (2) structural refinement using fs-TRXL data are compared in Extended Data Table 3, and a detailed comparison of these intermediate structures is provided in Supplementary Note 6.

It is worth noting that some of the χ_ν² values are near or well below unity. This is a well-known characteristic of X-ray scattering data analysis, arising from the fact that the effective number of independent information channels is fundamentally smaller than the number of measured q points. According to the Shannon sampling theorem, the independent information content of a scattering curve is governed by the Shannon number, N_s ≈ (q_max − q_min) × D_max/π, where q_max and q_min are the upper and lower limits of the measured q range, respectively, and D_max is the maximum interatomic distance contributing to the scattering signal. For a solute molecule such as azobenzene (D_max ≈ 12–15 Å), N_s is estimated to be approximately 25 for the q range of 1–7 Å⁻¹ used in this study. However, in solution-phase scattering experiments, the signal also includes bulk solvent response and cage contributions, which extend to distances much larger than the solute dimensions (>30–60 Å). This effectively increases the independent information content to N_s ≈ 60–120. Because the distance range over which these solvent contributions affect the scattering signal is not precisely defined, the effective D_max, and thus the true Shannon number, cannot be uniquely determined for solution scattering data. As adjacent q points carry redundant structural information beyond the Shannon limit, defining the degrees of freedom from the number of q points (n_q − d_f, with the small number of fit parameters d_f « n_q; Supplementary Method 3), rather than from the far smaller number of independent information channels (N_eff ≈ 110), results in systematically small χ_ν² values. Therefore, the observed χ_ν² values below unity do not indicate overfitting but rather reflect this inherent property of solution scattering data. A rigorous determination of N_eff using a BIC-based polynomial analysis is described in Supplementary Method 13. Moreover, one important caveat is that directly comparing the absolute χ_ν² values obtained for different SADSs (for example, comparing the χ_ν² of SADS₁ with that of SADS₃) can be misleading. This stems from the fact that the χ_ν² statistic is normalized by the uncertainties (σ_n(q_j)) of the extracted SADS, and in linear least-squares analysis, including the fitting of RSVs from SVD, the uncertainty of a signal component S_k is approximately inversely proportional to the norm of its concentration profile C_k(t):

This relationship implies that a species with a higher concentration and longer duration (that is, a larger ΣC²) will therefore have a much smaller uncertainty \({\sigma }_{{S}_{k}}\) and correspondingly amplified χ_ν² values. In our dataset, the cis-AB species persists up to the latest measured time delay without decay, giving it the largest integrated concentration profile among all species. As a result, the σ associated with SADS₃ (cis-AB) is systematically much smaller than those for the transient intermediates SADS₁ and SADS₂. Consequently, χ_ν² values for SADS₃ (or fits dominated by it) are mathematically amplified and tend to appear larger than those of SADS₁ and SADS₂, even for strictly comparable fit qualities. Therefore, comparisons of χ_ν² values among candidate structural models for the same species should be the primary criterion for model assessment, rather than cross-species absolute comparisons.

Systematic scan for parsimonious motions

To identify the minimal set of structural motions sufficient to describe each intermediate, we performed an SSPM. This procedure systematically evaluates all 15 possible combinatorial subsets of the four basis motions (4 single, 6 pairs, 4 triples, 1 quadruple) by fitting the simulated SADS from each candidate structure against the experimental SADS. Details are given in Supplementary Method 12. Model selection was performed using a χ²-based BIC, defined as BIC(χ²) = χ² + (k + 4) × ln(N_eff), where k is the number of structural motion parameters and the four additional parameters account for the mean-square displacement parameter σ², which accounts for thermal broadening of interatomic distances, the scaling coefficient, and the two baseline coefficients a_n and b_n, and N_eff is the effective number of independent data points determined by polynomial-order optimization (Supplementary Method 13). The model with the lowest BIC was selected as the optimal description; when two models were statistically indistinguishable (ΔBIC < 2), the more parsimonious model was preferred following the guideline of ref. ⁶³. The structural parameters, symmetry constraints and BIC values for all candidate models are tabulated in Supplementary Table 2.

Basis-motion analysis

To identify the motions governing structural transitions between pairs of three-dimensional molecular structures, we developed BMA. Given an initial structure A and a target structure B, BMA evaluates how closely B can be reproduced by applying restricted combinations of the four basis motions to A. The similarity is quantified by the Kabsch-aligned root-mean-square deviation (RMSD) of atomic coordinates. All 15 combinatorial models were systematically evaluated with multi-phase global optimization, and an RMSD-based BIC was used for model selection (Supplementary Method 14). The RMSD and ΔBIC(RMSD) values from BMA analysis for TRXL-refined structures and CASPT2-calculated structures are shown in Extended Data Fig. 7a–c and 7d–f, respectively, and summarized in Supplementary Tables 3 and 4, respectively.

Motional contribution analysis

To quantify the importance of each motion within an identified multi-motion optimal model, we developed a MoCA framework. The contribution of each motion was evaluated using a leave-one-out approach, which measures the increase in the optimized RMSD (or χ²) when a specific motion is excluded from the model while all remaining motions are re-optimized. This provides a robust measure of the functional necessity of each motion for reproducing the target structure, independent of the order in which the motions are applied (Supplementary Method 15). The MoCA results for the best models of the trans-AB → I₁, I₁ → I₂, and I₂ → cis-AB transitions are presented in Supplementary Table 5.

Ab initio calculations

We performed multi-reference quantum chemistry calculations to determine the molecular structures of key species involved in the trans-to-cis photoisomerization of trans-AB. Before geometry optimizations, the underlying energetic landscape of the S₀ and S₁ states was explored by constructing two-dimensional potential energy surfaces obtained from rigid scans along the CNNC dihedral angle and the symmetric CNN bending coordinate, while keeping all remaining internal coordinates fixed (Supplementary Fig. 7). These coordinates provide a physically transparent description of the photoisomerization pathway and allow straightforward identification of the trans minimum and the two S₁ local minima, as well as the approximate location of the cis minimum. These include the two ground-state structures, trans– and cis-AB, as well as the two intermediates in the S₁(nπ*) state, denoted as I₁ and I₂. Although previous studies^7,22 reported two minima on the S₁ state surface, the optimized geometries of these minima have not been thoroughly discussed. To address this gap, we carried out full geometry optimizations for the four key structures directly at the extended multistate CASPT2 (XMS-CASPT2) level, using the BAGEL package (version 1.2.0)⁶⁴, which provides fully analytic XMS-CASPT2 nuclear gradients^65,66; these make unconstrained geometry optimizations at this level of theory feasible for a molecule of this size.

The active space for these calculations included ten electrons distributed across eight orbitals, comprising two π and two π* orbitals from the phenyl rings, π and π* orbitals from the azo group, and the lone-pair orbitals on the nitrogen atoms. The CASSCF reference was optimized using a state-averaged CASSCF approach including the two lowest singlet states (S₀ and S₁(nπ*)). Previous studies have shown that substantial CASPT2 corrections are necessary to accurately describe the low-lying excited-state manifold of azobenzene, particularly for states such as S₁(nπ*) and S₂(ππ*). To incorporate these essential corrections, we used the extended multistate version of the perturbative method (XMS-CASPT2), which effectively accounts for interactions among the CASSCF states. An imaginary shift of 0.2 a.u., a standard value for CASPT2 excited-state calculations, was used to prevent numerical instabilities caused by intruder states (accidental near-degeneracies) in the perturbative treatment⁶⁷. Furthermore, symmetry constraints were disabled throughout the calculations, including geometry optimizations, to avoid potential bias. All calculations were performed in the gas phase using the aug-cc-pVDZ basis set with the cc-pVDZ-JKFIT auxiliary basis for density fitting. Solvent effects were not included in these geometry optimizations, but were incorporated in the atomic partial charge calculations for the cage-term calculation (Supplementary Methods 8 and 9). Supplementary Note 5 compares the CASPT2-calculated structures with previously reported theoretical structures.

Ethics statement

This study did not involve human participants, animals or cell lines, and no ethical approval was required.