Light-induced hexatic state in a layered quantum material

The tunability of materials properties by light promises a wealth of future applications in energy conversion and information technology. Strongly correlated materials such as transition metal dichalcogenides offer optical control of electronic phases, charge ordering and interlayer correlations by photodoping. Here, we find the emergence of a transient hexatic state during the laser-induced transformation between two charge-density wave phases in a thin-film transition metal dichalcogenide, 1T-type tantalum disulfide (1T-TaS2). Introducing tilt-series ultrafast nanobeam electron diffraction, we reconstruct charge-density wave rocking curves at high momentum resolution. An intermittent suppression of three-dimensional structural correlations promotes a loss of in-plane translational order caused by a high density of unbound topological defects, characteristic of a hexatic intermediate. Our results demonstrate the merit of tomographic ultrafast structural probing in tracing coupled order parameters, heralding universal nanoscale access to laser-induced dimensionality control in functional heterostructures and devices.


C
ollective behaviour beyond the single-particle picture fosters a variety of fundamental physical phenomena, such as superconductivity 1 and ferromagnetism 2 , which involve the establishment of longrange order below a critical temperature.Due to weaker screening and reduced phase space volumes for scattering, collective excitations are particularly important in low-dimensional systems.Prominently, plasmon and spin waves constitute elementary excitations of the onedimensional Luttinger liquid 3 .
For ideal systems of reduced dimensionality 4 , the Mermin-Wagner-Hohenberg theorem prohibits long-range ordering at finite temperatures, owing to thermal fluctuations 5 .Strictly speaking, this precludes symmetrybreaking phase transitions.
Instead, the transition from a quasi-long-range ordered two-dimensional solid to a liquid phase is described within the framework of Kosterlitz-Thouless-Halperin-Nelson-Young (KTHNY) theory 6,7 .While the solid-liquid transition is always of first order in three-dimensional systems 8 , the theory predicts two successive second-order transitions in two dimensions via a hexatic phase.This intermediate is characterised by intact orientational but reduced translational order arising from the presence of unbound topological defects 7,9 .
Experimental realisations of this transition have been actively pursued with two-dimensional model systems of colloidal spheres 9 , particles physisorbed on crystalline substrates 8 , skyrmion lattices 10 and technologically relevant smectic liquid crystals 11 .Yet, quasi-twodimensional character is also readily found in layered materials 4 .Among those, transition metal dichalcogenides (TMDCs) stand out due to their highly anisotropic correlations 12 .This leads to the occurrence of directto-indirect band gap transitions 13 , layer-and dopingdependent superconductivity [14][15][16][17] , and the formation of charge-density waves (CDWs) 18 .Due to the delicate interplay between the various degrees of freedom, the properties of these systems, often referred to as quantum materials 19 , are highly tunable by external stimuli [20][21][22][23][24][25] .
The effective dimensionality of a CDW system is determined by the coupling between neighbouring layers.Hence, a phase transition that modifies interlayer correlations may lead to the occurrence of a dimensional crossover 50,51 .On ultrafast timescales, CDW stacking dynamics have recently been studied by means of ultrafast electron 39,44 and x-ray diffraction 43,52 .In particular, an optically induced breakdown of excitonic correlations was shown to quench out-of-plane order 44 , and careful tuning of the photoexcitation density was used to realise interlayer phase boundaries with two-dimensional characteristics 53 -raising the question if optical control may serve as a gateway to fundamentally low-dimensional phenomena such as the KTHNY transition.
In this article, we report on the observation of transient hexatic order in a nascent incommensurate CDW phase.We use ultrafast high-coherence nanobeam diffraction to track the three-dimensional phase ordering after femtosecond optical excitation.At early times and prior to the establishment of the equilibrium stacking sequence, we identify a quasi-two-dimensional state characterised by a pronounced anisotropic broadening of diffraction spots.We reproduce the experimentally observed behaviour in time-dependent Ginzburg-Landau simulations and attribute it to the presence of unbound topological defects, consistent with the predictions of KTHNY theory.
As the key experimental innovation of the present study, we harness these capabilities to conduct femtosecond electron diffraction combining high reciprocal-space resolution with a particularly narrow electron beam, in extended series of eight beam tilts.A concomitant transverse coherence length of up to 120 nm enables a precise measurement of in-plane spot profiles and allows to distinguish phases characterised by similar periodicities.Simultaneously, the nanometric probe beam guarantees diffraction from a sample region of sufficient homogeneity in terms of thickness and orientation, as required for the quantitative investigation of stacking dynamics.The fingerprint of such structural modifications is often encoded in lowintensity diffracted signals.We increase the sensitivity to these features by means of a sample design tailored to drive the transformation at an unprecedented laser repe-tition rate of 610 kHz, maximising the duty cycle of our measurement scheme.
In the experiments (Fig. 1), we excite a free-standing 70 nm thin film of the prototypical CDW material 1T -TaS 2 at room temperature, using ultrashort laser pulses (800 nm wavelength, 150 fs duration, between 0.3 mJ/cm 2 and 7.0 mJ/cm 2 fluence).As a function of a variable temporal delay ∆t, we capture the transient distribution of CDW spot intensities and profiles using high-coherence ultrashort electron pulses (120 keV beam energy, 500 fs duration, between 170 nm and 1.1 µm spot size, below 0.1 mrad convergence semi-angle, see Supplementary Information and Fig. S3 for fluence-dependent delay scans).
Incommensurate CDW phases in 1T -TaS 2 .We investigate the photoinduced transition between two incommensurate CDW superstructures in 1T -TaS 2 (grey and red regions in Fig. 1b).In thermal equilibrium and at temperatures above 353 K, the in-plane modulation wave vectors q i,IC (with i ∈ {1, 2}) of the CDW/PLD are aligned along the lattice vectors a i of the hexagonal host 66 .Due to the presence of a gap-less long-range phase fluctuation (or phason) mode, the CDW in this high-temperature incommensurate phase (in short IC phase) is effectively freefloating and only weakly coupled to neighbouring layers 67 .As such, the IC phase in 1T -TaS 2 is an ideal host for to- pological defects and, potentially, hexatic ordering 68 .
Below the phase transition temperature, the modulation remains incommensurate and the wave vectors q i,NC form an angle of ∼12°with the lattice vectors a i (nearly commensurate or NC phase).With significant contributions from higher harmonics, the realspace structure exhibits a network of discommensurations separated by domains of commensurate character in which the CDW/PLD is locally locked-in with the crystal structure 66,69 .Despite their different in-plane structures, both phases exhibit a commensurate, three-fold stacking periodicity along the out-of-plane direction ( a * 3 = 3 a 3 , see Fig. 1b) 66,69 .
The reciprocal lattice of the host material is given by basis vectors b i .The modulation wave vectors q i span the reciprocal CDW lattice (m 1 , m 2 ) around every individual main lattice point (h, k, l) (Fig. 2a).Hence, each lattice point k can be identified by 69 This notation absorbs the commensurate out-of-plane component of the CDW in a supercell, effectively leading to a shorter reciprocal lattice vector b * 3 = b 3 /3 and certain systematic absences in diffraction experiments.Specifically, whereas all main reflections (with m 1 , m 2 = 0) lie in planes l = 3n, first-order CDW spots with (m 1 , m 2 ) = (1, 0), (0, 1) and (−1, 1) (and those with opposite sign) are located in planes l = 3n ± 1 with non-vanishing k zcomponents.Thus, they only appear upon tilting the electron beam away from the [001] zone axis 39,69 .
A visualisation of the corresponding diffraction geometry is shown in Fig. 2b.Under tilted-beam conditions, electron diffractograms typically feature spots in more than one Laue zone (Fig. 2d).Main reflections are located in the zero-order Laue zone (ZOLZ) close to the unscattered beam and appear bright compared to the second-order CDW spots surrounding them.As a result of the CDW stacking sequence, first-order CDW reflections are found in the first-order Laue zone (FOLZ) and at larger wave vectors for moderate beam tilts (see also Supplementary Movies S1 and S2).3a shows a series of time-resolved diffraction patterns, averaged over several main reflections in the FOLZ.Before optical excitation (time-zero), we observe a triplet of sharp first-order spots of the initial NC phase.The additional low-intensity reflections in-between stem from higher CDW orders (> 2).After excitation, the NC spot intensity is largely suppressed within 1 ps.Simultaneously, nearby reflections with (m 1 , m 2 ) = (1, 0), (0, 1) and (−1, 1) emerge, which evidence the nascent IC phase and exhibit a significant anisotropic broadening 40,41,70 , as well as a few-picosecond increase in intensity 38 (cf.also ref. 39 ).Additionally, a slight increase of scattered intensity is detected at (m 1 , m 2 ) = (−1, 0), (0, −1) and (1, −1), i.e., opposite of the bright IC spots.Only after these earlystage dynamics, the spot shapes become isotropic, and solely the bright IC spots remain (see image at 100 ps).

Dynamics of the in-plane correlation length. Figure
The different spot profiles along azimuthal and radial directions (in relation to the nearby main reflection) is a most characteristic hallmark for the presence of a hexatic phase 7 .For a quantitative analysis, we recorded a second image series with a longer camera length, i.e., optimised reciprocal-space resolution (Fig. 3b).We fit the spot shape at every temporal delay (see insets) and extract the azimuthal and radial spot widths (indicated by brown and blue arrows, respectively).The results are shown in Fig. 4b (brown and blue data points).From a spot width ratio > 2 shortly after optical excitation (blue circles in Fig. 4c), the reflections assume an isotropic, yet broadened shape within 10 ps (temporal regime I).On longer timescales, the IC spot width finally approaches that of the NC peak measured before time-zero (temporal regime II; dashed grey line in Fig. 4b).Notably, this behaviour coincides with a monotonous growth of the IC wave vector, which initially is ∼2 % shorter than in the equilibrium phase at late delay times (black circles in Fig. 4c).
Reconstruction of the CDW rocking curve.The observed evolution of the in-plane spot profile suggests the presence of hexatic order during temporal regime I.In order to test this hypothesis, we probe the effective dimensionality of the system as induced by the optical excitation, analysing intensities of reflections with different out-of-plane momenta.
Under tilted-beam conditions, the lattice rod is densely sampled across the various spots present in the FOLZ of every single diffractogram (see visualisation in Fig. 2b and colour overlay in Fig. 2d for the experimental tilt angle of ∼2.4°).To further enhance sensitivity and fully cover reciprocal space, we perform delay scans with eight different beam tilts.Correlating the scattered intensities before time-zero with dynamical diffraction simulations, we determine the tilt angle for any of these scans with an accuracy better than 0.2°(see Supplementary Information and Fig. S1 for a more detailed description).
Sorting scattered intensities of spots with comparable CDW structure factors (Fig. 2c) by their associated outof-plane components k z , we reconstruct the CDW rocking curve from 11 individual IC reflections with a total of 72 values of k z throughout our tilt series at every stage of the dynamics (Fig. 3c and Supplementary Information).At early times, the intensity increases independently of k z , indicating a pronounced elongation of the reciprocal lattice rods along the b * 3 direction (see Fig. 3c at 3 ps and Supplementary Fig. S2).Subsequently, a well-defined spot profile emerges and narrows until approaching that of the NC phase measured before time-zero (red data points and dashed grey line in Fig. 4b).

Discussion
The two CDW modifications involved in the experiment are not commensurate with each other.As a result, there is no continuous global deformation of the crystal during the phase transition.Instead, photoexcitation gives rise to a high density of topological defects, causing the broadening of diffraction spots directly after time-zero.In the past, both CDW dislocations 40,42,43,70 and domain walls 41,43 have been invoked to describe the transiently disordered structure on a microscopic scale.While CDW states in 1T -TaS 2 are known to support both types of topological features 47 , the kinetics and interplay between defects in-and out-of-plane remain largely uncertain.
We interpret the evolution of the reconstructed diffraction spots in terms of a transient hexatic phase supported by a temporary suppression of interlayer correlations.Specifically, the observed anisotropic in-plane broadening of CDW reflections (Fig. 4d) indicates a pronounced loss of translational order.This feature is characteristic of KTHNY behaviour, relating the condensation of a twodimensional crystalline phase via a hexatic intermediate to the pair-wise annihilation of point-like defects 7 .In our measurements, the occurrence of such a mechanism is further substantiated by the temporal evolution of the IC wave vector (black circles in Fig. 4c) whose initial shortening has recently been attributed to a dislocation-induced self-doping of the electronic system 43 .The concurrent, highly uncorrelated out-of-plane structure is apparent from the elongation of reciprocal lattice rods (Fig. 3c and Fig. 4b), giving rise to the simultaneous visibility of all six IC reflections in the same Laue zone (Fig. 3a, 1 ps to 5 ps).
To further explore the notion that a gas of dislocationtype topological defects governs the initial stages of the IC phase formation kinetics, and to understand how threedimensional ordering affects this process, we implement a multilayer time-dependent Ginzburg-Landau simulation.The underlying free-energy functional is derived from Nakanishi's model 71 describing the most important CDW phases in 1T -TaS 2 and their out-of-plane configurations (see Methods for a more detailed description).We express the IC modulation in terms of three coupled order parameters φ j,l corresponding to the directions of the three in-plane modulation wave vectors q i,IC and the layer l in the atomic structure.The physical charge-density modulation (see Supplementary Fig. S4) is then given by ( This model is capable of closely reproducing the experimentally observed temporal evolution of the reciprocal lattice rod (solid lines in Fig. 4b).For the associated set of parameters, the early-stage dynamics are characterised by a high density of uncorrelated phase singularities (Fig. 4f).
In this state, local energy minimisation leads to a spatially anisotropic character of the phases of the individual order parameters (Fig. 4e), causing the typical anisotropic inplane spot shape shown in Fig. 3b.From here on, a pairing of singularities with opposite chirality in the phases of two of the three order parameters leads to the formation of physical CDW dislocations, and enables a local build-up of the PLD amplitude.The evolution of the average free energy associated with an individual phase vortex parallels that of a purely two-dimensional system, as seen for a reference simulation without interlayer coupling (red and grey curves in Fig. 4d, respectively).We therefore believe that the observed loss of stacking order at early times is a prerequisite for inducing the intrinsically two-dimensional hexatic intermediate in 1T -TaS 2 .It is intriguing to speculate whether such a state could in fact be realised as a thermodynamically stable phase during the metallic-to-IC transition in a 1T -TaS 2 mono-layer.
Beyond temporal delays of 10 ps and after the spot shape anisotropy has decayed, we observe the onset of three-dimensional behaviour.In contrast to the twodimensional reference system, the energy of individual vortices increases, and their resulting correlation in adjacent layers leads to the formation of 'flux lines' (black lines in Fig. 4d and g, see Methods for a more detailed description).In parallel, the establishment of long-range order is driven by the annihilation of singularities within each of the φ j,l , eventually reaching correlation lengths at late times that exceed the transverse electron beam coherence length in the experiment.
In conclusion, the identification of hexatic order and stacking dynamics in our study is enabled by ultrafast diffraction with a nanoprobe of exceptional collimation.These results tie in with recent observations of laser-induced transient phases 34 and dimensional crossovers 44,53 , exemplifying how optical control over interlayer correlations can be used to host low-dimensional states and transitions.Addressing a recurring experimental challenge in ultrafast and materials science, this approach promises to advance high-resolution non-equilibrium investigations in systems characterised by weak structural signatures, sub-micron sample sizes and considerable spatial heterogeneity 72 .As such, nanoscale structural analysis will guide the design of applications harvesting laser-induced functionality by material composition and tailored responses to external stimuli.

Methods
Ultrafast nanobeam diffraction experiments.The Göttingen UTEM is based on a JEOL JEM-2100F transmission electron microscope modified to enable the investigation of ultrafast dynamics.Femtosecond laser pulses (515 nm wavelength after frequency doubling of the output of a 'Light Conversion PHAROS' femtosecond laser, 610 kHz repetition rate) are used to generate ultrashort electron pulses from the microscope's ZrO/W Schottky emitter via linear photoemission.A fraction of the laser output is converted to 800 nm wavelength by optical parametric amplification ('Light Conversion ORPHEUS-F') and incident on the sample at a variable temporal delay with respect to the electron pulses.Further technical details on the instrumentation are given in ref. 61 .Snapshots of the non-equilibrium dynamics are recorded on a direct electron detection camera ('Direct Electron DE-16') and processed by an electron counting algorithm.The diffractograms presented in this article have been integrated for 3 min (Fig. 3a), 7.5 min (Fig. 3b), and 5 min (Fig. 2d and Fig. 3c) per temporal delay, respectively.Specimen preparation.The investigated 70 nm thin film of 1T -TaS 2 has been obtained by ultramicrotomy.Details on the preparation process and a comprehensive characterisation of the specimen can be found in the Supplementary Information of ref. 35 .Time-dependent Ginzburg-Landau simulations.In the time-dependent Ginzburg-Landau simulations discussed in the main text, we describe the physical CDW modulation ρ l in layer l of the material as given in equation 2. This notation effectively strips the equilibrium in-plane CDW periodicity q j from the additional out-ofequilibrium modulation given by the three coupled order parameters φ i,l .Numerically integrating the equation of motion for the phenomenological free energy F of the sys-tem then yields the spatiotemporal dynamics shown in Fig. 4, starting from a stack of uncorrelated layers where both amplitude and phase of the φ j,l are randomised.The parameter Γ controls the overall timescale of the dynamics.For the free-energy functional, we choose an approach tailored to model the phase diagram of the different threedimensional CDW configurations in 1T -TaS 2 71 , i.e., + ae ig 2 φ * j,l φ j,l+2 . ( Therein, energy minimisation of the B and the C-term ensures a local equilibration of the CDW amplitude and the phasing-term D gives the relative phase relation between the order parameters such that well-defined CDW maxima emerge in a hexagonal arrangement (see Supplementary Fig. S4).The transition from a local lock-in of the CDW with the underlying main lattice in the C-phase to the incommensurate modulations found above a critical temperature T C is governed by a temperature-dependent competition between the commensurability energy E and the kinetic energy A. In an out-of-equilibrium scenario, a minimisation of the latter also determines the kinetics of topological defects.In order to account for the temperature-dependent relative orientation between the NC CDW wave vector and the main lattice periodicities 66 , the kinetic energy includes a softness towards a distortion of the order parameter along the azimuthal component, which is in conceptual agreement with the modelling of the KTHNY transition of a two-dimensional solid 7 .One finds 71 : with where ϕ j describes the angle between the wave vector q and q j , b j are the reciprocal lattice vectors of the undistorted structure and u, v and s are parameters.
Along the out-of-plane components, equation 4 perturbationally treats the coupling of the individual layers to their nearest neighbouring and next-nearest neighbouring layer via the parameters G and a, respectively, while the finite phase factors g 1 and g 2 ensure the establishment of the expected stacking periodicity 71 .
Starting from a randomised phase and amplitude pattern in all φ j,l , we find that a parameter set with B = 5 In order to extract spot widths from the simulation, we perform a three-dimensional Fourier-transformation of the real parts of the φ j at representative stages of the temporal dynamics to derive the simulated reciprocal lattice rod.Summing the squared modulus of the momentum distribution in planes corresponding to the outof-plane momentum | b * 3 | (considering the different rotations of the respective q j ) for all three φ j then gives the simulation analogue to the experimental diffraction spots depicted in Fig. 3b.Fitting a two-dimensional Lorentzian function to this in-plane spot-profile yields the corresponding spot widths γ sim along the azimuthal and radial directions.Along the out-of-plane directions, we assign the integrated intensity of the CDW diffraction spot in every reciprocal lattice plane to its out-of-plane component and fit a Gaussian function centred at the corresponding equilibrium stacking periodicities | b * 3 | to the resulting one-dimensional intensity distribution.The derived spot widths γ (FWHM) over the course of the dynamics are displayed in Fig. 4b.For adequate comparison between simulation and experimental results, we include the instrument resolution γ r in the form γ = γ 2 sim + γ r . 2 where γ r is the corresponding spot width measured at late times.Within this model, the phase formation is governed by the kinetics of point-like defects that appear as vortices in the φ j,l .For the 'flux lines' length shown in Fig. 4d, we assume vortices in neighbouring layers as correlated when their spatial separation along the in-plane coordinates is less than five CDW lattice vectors.At early times, this definition of the interlayer vortex correlation is additionally influenced by coincidental alignment stemming from the high initial density of defects within every layer.We correct for this influence by subtracting the temporal evol-ution of the 'flux line' length derived in the same manner for the case of uncoupled CDW layers.
The average energy of an individual vortex for both the coupled and uncorrelated stack of CDW layers is derived by integrating equation 4 at every step of the dynamics and additionally considering the ground-state energy of the equivalent fully equilibrated system.For the latter, we chose a uniform amplitude within the φ j,l and a relative phase shift of 2π/3 between neighbouring layers as the initial condition and let the local amplitude relax until the system reaches its energy minimum.

Reconstruction and fitting of CDW rocking curves
In our experiments, we reconstruct the out-of-plane shape of the CDW reciprocal lattice rod by assigning the respective k z -component to the measured intensity of every CDW diffraction spot in each individual image of the image series.In addition to a modulation by k z , the intensity of a CDW spot of order n at a scattering vector k depends on the structure factor V k , which can be expressed in terms of a Bessel function of first kind J n (x): For the longitudinal distortion with comparably small amplitude A i found here ( A i || q i , where q i is the CDW wave vector) 2 and considering high-energy electrons, the diffractograms mainly contain information about the modulation of the tantalum sublattice 3 .Additionally, CDW spots, to a good approximation, can be sorted by their respective scattered intensity into three groups such that the CDW rocking curve can be reconstructed from multiple CDW spots within a single image (Fig. S1d-f).For the CDW rocking curves displayed in the main text, we only consider the brightest reflections where the angle ϕ between k and q i (see Fig. S1b) is smaller than 45°(see Fig. S1d for the resulting NC rocking curve at ∆t < 0 ps).In comparison, CDW reflections with ϕ > 45°are less intense and more heavily affected by dynamical scattering effects (Fig. S1e-f).The reconstruction of the CDW rocking curves is based on the k z -dependent temporal evolution of scattered CDW intensities shown in Fig. S2.Following optical excitation, the NC CDW amplitude is suppressed within 1 ps irrespective of the out-of-plane component.In contrast, IC diffraction spots situated at k zcomponents deviating from the FOLZ display an initial intensity overshoot, whereas those spots probed directly in the FOLZ experience a monotonous increase of scattered intensity over the entire temporal regime (Fig. S2a).
The length of the reciprocal lattice rod at every temporal delay is determined by fitting a Gaussianshaped profile to the rocking curves inferred from these delay curves.Every rocking curve consists of intensities scattered into 11 individual IC reflections probed at varying k z in our tilt series, resulting in a total of 72 out-of-plane momenta considered for the fitting of the rod * Corresponding author: claus.ropers@mpinat.mpg.deshape.The majority of these CDW reflections lie at outof-plane wave vector components smaller than the equilibrium stacking periodicity of the CDW.In order to enhance the fitting weight of the brightest spots close to the centre of the rocking curve, we linearly interpolate the measured CDW intensity to gain a homogeneously spaced distribution along k z .The fit results with the corresponding error bars are displayed in Fig. 3c and Fig. 4b in the main text.

Dynamical electron diffraction simulations
Within a single diffractogram, scattered CDW intensities depend on the excitation error, i.e., the mismatch to an exactly fulfilled Laue condition, that is determined by the relative tilt angle between specimen and electron beam.In order to determine this tilt angle from measured intensities in our tilt-series, we perform dynamical diffraction simulations based on the scattering matrix approach.Therein, the incident electron wave function inside the crystal is divided into individual beams n scattered by a respective crystal periodicity k n 4 .The beam amplitude Ψ n as a function of the penetration depth z can then be described as Equation S2 is the Darwin-Howie-Whelan equation for elastic scattering of high-energy electrons.It can be expressed in terms of a scattering matrix S N 5 and then integrated numerically in close analogy to the treatment of N coupled harmonic oscillators.The diagonal elements of this matrix account for the scattering geometry in the form of the excitation error s n , describing the deviation from an exactly fulfilled Laue condition.The transfer of scattering amplitude between beams n and n is mediated by the off-diagonal elements ξ n,n representing the effective scattering length, thereby incorporating the influence of the specimen.Specifically, 1 ξ k ∝ λ|V k | where λ is the electron wavelength and V k the Fourier component of the scattering crystal potential, i.e., the structure factor associated with the wave vector k.
Our simulations are based on the reported crystal structure of the NC phase 6 where we approximate the incommensurate lattice distortion by a supercell comprised of 45 × 45 × 3 unit cells of the undistorted structure (see ref. 7 for a justification of this approximation), while the outof-plane periodicity of the NC phase is given by the σ h2 stacking described in detail in ref. 8 .The structure factor of the material is calculated with the atomic scattering arXiv:2207.05571v2[cond-mat.mtrl-sci]21 Jul 2022 factors for tantalum and sulphur listed in table I 3 , yielding scattering lengths of, e.g., ξ = 108 nm for a scattering vector relating a first-order main lattice reflection to the direct beam, and ξ = 1545 nm for a first-order NC spot relative to its associated main lattice spot.
In our simulations, we consider 61 main lattice reflections situated in the ZOLZ, in addition to a total of 786 CDW spots of first and second-order in the FOLZ and in the ZOLZ, respectively.Thereof, only those wave vectors yielding a scattering strength of more than 2 % of the maximum found at q = 0 are taken into account.Furthermore, we include the effects of finite temperatures on scattered intensities in the form of the reported Debye-Waller factors of tantalum and sulphur at roomtemperature 9 .A numerical integration of equation S2 with varying sample thickness and electron beam incidence then allows us to fit the overall specimen tilt following the approach described in the subsequent section.
A resulting simulated diffractogram for an exemplary beam incidence from our measurements is depicted in Fig. S1b.We account for a spatially varying sample orientation and the finite electron beam convergence angle by averaging over beam tilts within a range of 1.5 mrad centred around the fitted beam tilt.Since electron diffractograms, especially those recorded close to zone-axis orientation, always feature multiple diffraction spots of significant intensity, the average 70 nm thickness of our sample results in non-negligible effects of dynamical scattering on the measured intensities, leading to pronounced deviations from kinematical scattering theory.In particular, diffraction spots of the undistorted lattice situated in the vicinity of the Laue circle in the ZOLZ display strong intensity oscillations during a (virtual) beam tilt series, and, additionally, are highly susceptible to small variations of the tilt direction (dotted rocking curves in Fig. S1c).In contrast, first-order CDW reflections found at larger wave vectors are less affected by multiple scattering events (Fig. S1g).For these spots, the simulations yield only minor intensity modulations along k z and a small shift of the rod centre.

Fitting of specimen tilt angles
In any TEM, the electron beam angle of incidence can be set by beam deflectors, the resulting shift of the diffractogram yields a reliable measure of the relative beam tilt between two settings.The absolute beam tilt, however, is encoded in diffracted intensities and requires fitting routines based on diffraction simulations 10 .
Using the simulation approach outlined in the previous section, we calculated more than 32 000 diffractograms with variable angles of incidence, while averaging over the thickness distribution of our specimen (see supplementary material of ref. 11 ) and a relative angle of incidence between specimen and electron beam 1.5 mrad, considering both the electron beam convergence angle as well as a spatial variations of the sample orientation.The simulated intensity T i,j of a diffraction spot j for every incidence i is then compared to the experimental intensity P j measured before time-zero.
The strong impact of dynamical scattering leads to ambiguous fitting results for the individual T i .To overcome this limitation, we additionally make use of the known relative beam tilts ∆ in our tilt series, consisting of delay scans recorded under eight different angles of incidence.Summing over the intensities of all main lattice and firstorder CDW reflections, this approach yields the image correlation index 10 where Q i = 1 indicates a perfect agreement between experiments and simulation.
The calculated Q i are depicted in Fig. S1h.We derive a maximum correlation index of Q i = 0.89.The corresponding beam tilt results in the simulated intensities of main lattice spots shown in Fig. S1c and yields rocking curves for first-order CDW reflections that are centred around the wave vector of the equilibrium CDW stacking sequence (Fig. S1d-f).

Fluence-dependent ultrafast nanobeam diffraction
In our experiments, we employ a sample design (described in more detail in ref. 11 ) that is tailored to minimise the cumulative heating of our 1T -TaS 2 specimen.This enables laser pump/electron probe experiments at exceptionally high laser repetition rates (up to 610 kHz) over a wide range of pump fluences (here: up to 7.0 mJ/cm 2 ).The thusly enhanced duty cycle in our measurements allows us to further limit the diameter of the collimated electron beam by apertures down to 170 nm in the specimen plane.The resulting delay curves of CDW intensities scattered into first and second-order CDW spots are shown in Fig. S3a and b, respectively.Following the optical excitation, intensity scattered into both first and second-order NC CDW spots is suppressed within 1 ps.For higher pump fluences, IC CDW spots emerge on a similar timescale.
Movie S1 Electron beam tilt-series around the ZOLZ.Tilting the electron beam around the ZOLZ modulates diffracted intensities.For higher beam tilts, first-order CDW reflections appear at larger wave vectors.
Movie S2 Electron beam tilt-series around the FOLZ.Higher beam tilts shift the probing of first-order CDW reflections to smaller wave vectors. .h Fitted quality factors of the template matching described in the text as a function of the k x , k y -components of the normalised incident wave vector k 0 .The black circle indicates the fitting results, the tilt direction for the rocking curves in c and g is given by the solid black arrow.A similar virtual tilt series only shifted by 0.29°along the direction given by the dotted arrow yields the dotted rocking curves displayed in c.The scattered intensity at the IC spot position before time-zero is subtracted for both types of curves as a measure of the inelastic background in the images.The NC curves are additionally normalised with respect to the measured intensity before time-zero and show similar behaviour as a function of temporal delay (cf.Fig. 3b).For the IC curves, we normalise to the measured intensity after 300 ps.During the phase formation, IC spots probed at a k z -component closer to the ZOLZ display an initial overshoot (blue curve).In contrast, reflections in the FOLZ are characterised by a continuous intensity increase (red curve), leading to the sharpening of the reciprocal lattice rod along the out-of-plane direction described in the main text and Fig. 3. b Average intensity of first-order NC (blue), first-order IC (red) and second-order NC spots (green).Curve normalisation as described in a.  S2.We observe an initial suppression within 500 fs that scales with the pump fluence.Above a threshold of around 2 mJ/cm 2 , the measured NC intensity after time-zero equals the inelastically scattered background.For such higher fluences, we measure an increase of IC spot intensity within 5 ps with a temporal characteristic independent of the pump fluence.b Delay scans of second-order NC intensity for varying pump fluence with the sample tilted into the zone-axis.In contrast to a, the electron beam has a larger diameter (1.5 µm), adding a spatiotemporal component to the dynamics following the initial suppression.For intermediate fluences, we witness a partial recovery of NC intensity after around 5 ps, followed by a second suppression at later delays.We attribute this behaviour to a locally varying excitation density leading to the formation and subsequent growth of IC phase domains on longer timescales as described in ref. 11 .

Fig. 1
Fig. 1 High-coherence ultrafast nanobeam diffraction.a Schematic of the Ultrafast Transmission Electron Microscope (UTEM).A specimen (inset in b) is excited using ultrashort laser pulses (red).After a variable temporal delay ∆t, ultrashort electron pulses (green) are emitted from a tip-shaped photocathode via linear photoemission (blue) and capture the specimen's transient state on a direct electron detection camera (inset in c).b Visualisation of the 1T -TaS 2 specimen located within the objective lens of the UTEM instrument.The optical excitation drives a transformation of the layered CDW lattice between two incommensurate CDW phases (represented by grey and red spheres), probed by a nanometer-sized electron beam of variable tilt (double arrow).c Schematic diffraction patterns for two temporal delays.CDW spots (red) are arranged around a main lattice peak (grey).The temporal evolution of CDW spot intensities and profiles (radial as well as azimuthal) encodes both in-and out-of-plane CDW dynamics on ultrashort timescales.

Fig. 2
Fig. 2 Reciprocal lattice of 1T -TaS 2 , diffraction geometry and structure factor.a Reciprocal lattice of the NC (open grey) and IC phase (red) of 1T -TaS 2 (first-order spots only).Due to the CDW stacking periodicity, the CDW spots are located in the first-order Laue zone (FOLZ, l = 3n ± 1 in Eq. 1) of the CDW lattice, while the related host lattice reflections (solid grey) lie in the zero-order Laue zone (ZOLZ).b Ewald construction visualising the diffraction geometry.The intensity scattered into individual first-order CDW spots is governed by proximity to the Laue condition, i.e., where the Ewald sphere (green) intersects the CDW reciprocal lattice rods (red).Accessing the FOLZ requires an electron beam slightly tilted away from the [001] zone axis.The low curvature of the Ewald sphere for 120 keV electron energy results in a dense sampling of the rod shape.c A second important contribution to the CDW spots intensity is given by the CDW structure factor depending on the angle between scattering vector and CDW wave vector (see Supplementary Information).d Temporally averaged experimental diffractogram.The high electron beam coherence leads to a clear separation of NC and IC spots.The coloured overlay indicates the local height k z of the Ewald sphere above the ZOLZ in units of b * 3 .While first-order spots of both phases are the dominant feature in the FOLZ and contain information on the stacking periodicity, the low-intensity second-order NC spots appearing between bright host reflections in the ZOLZ are a sensitive measure of the CDW amplitude.

Fig. 3
Fig. 3 Temporal evolution of the in-plane and out-of-plane CDW spot shapes.a Section of reciprocal space in the FOLZ at representative temporal delays ∆t.A triplet of first-order CDW diffraction spots surrounds the associated structural reflection, satellites with lower intensities are higher-order spots of the NC phase attributed to neighbouring reflections.b Average over all visible first-order CDW diffraction spots close to the FOLZ.The equilibrium IC spot position is given by dashed red circles.Brown and blue arrows indicate the evaluation axes for the azimuthal and radial spot widths.Insets: Two-dimensional pseudo-Voigt fit of the spot shape.c Out-of-plane spot profiles at representative temporal delays ∆t reconstructed from in-plane spot intensities measured at different k z (red circles).The scaled equilibrium spot profile of the NC phase before time-zero is displayed for reference (grey circles).Solid lines represent Gaussian fits of the rocking curves.The spot widths of the NC phase before time-zero and of the IC phase at 400 ps are indicated by arrows and amount to 0.23/nm (FWHM), respectively.

Fig. 4
Fig. 4 Dynamics of the IC spot profile and simulation results.a Schematic temporal evolution of the reciprocal lattice rod shape.b Azimuthal (brown) and radial (blue) IC spot width (FWHM) extracted from the in-plane data shown in Fig. 3b (see brown and blue arrows at 10 ps delay) and out-of-plane IC spot widths (red) derived from the reconstructed rocking curves shown in Fig. 3c.All curves approach the NC equilibrium spot widths at late delay times (dashed grey lines).Solid lines represent the corresponding temporal evolution of the spot widths in the simulation corrected for the instrument resolution (see Methods).c Left axis: Temporal evolution of the ratio between azimuthal and radial spot width (blue).Right axis: Dynamics of the CDW wave number | q IC | extracted from the fits in Fig. 3b (black).d Left axis: Simulated energy per phase vortex for a stack of correlated (red) and uncorrelated (grey) layers, respectively.Right axis: Average flux line length corrected for the influence of coincidental alignment present at early times (see Methods).e Exemplary simulated diffractograms in a section of reciprocal space as in Fig. 3a.In accordance with the experimental results, we find six azimuthally broadened diffraction spots in the FOLZ at early times, while only three CDW reflections indicative of the three-dimensional CDW are present at late times.f Simulated phase modulation of φ 3,l within an individual layer.g Three-dimensional CDW phase modulation of φ 3,l in 4 of the 51 simulated layers.The stacking sequence favours a local phase shift of 2π/3 between adjacent layers and the alignment of phase vortices into 'flux lines' (marked as black circles and diamonds for opposite chirality, respectively).The images display the same region as in f.

Fig. S1 Fitting
Fig. S1 Fitting of specimen tilt angle and reconstructed rocking curves.a Exemplary electron diffractogram recorded with one of the overall eight different specimen tilt angles used to reconstruct CDW rocking curves.The displayed image is a temporal average over all stages of the phase formation kinetics as can be recognised from the presence of spots corresponding to both CDW phases.Direct electron detection enables simultaneous recording of scattered intensities spanning five orders of magnitude (see intensity scale).b Simulated electron diffractogram for the NC phase of 1T -TaS 2 for the fitted tilt angle in a. Coloured circles mark the spots belonging to the measured and simulated rocking curves in c, d and g, respectively.Scattered CDW intensity is additionally modulated by the CDW structure factor which predominantly depends on the angle ϕ (cf.equation S1).c Measured (black circles) and simulated (solid lines) rocking curves of the main lattice reflections indicated in b.The dotted curves describe the rocking curve for an additional sample tilt of 0.29°along the direction indicated by the dotted arrow in h, exemplifying the pronounced influence of dynamical scattering on bright intensity reflections situated in the ZOLZ.Measured data is rescaled by a global factor for all curves to match the simulated intensities.d-f Measured NC CDW rocking curves for the fitted specimen tilt angle and sorted by the CDW structure factor (cf. equation S1).The bright-intensity CDW spots in d are used to extract the out-of-plane shape of the reciprocal lattice rod discussed in the main text.Especially CDW spots with smaller structure factor (e and f) are subject to dynamical scattering.g Simulated CDW rocking curves for the spots highlighted in b.Overall, dynamical scattering has a smaller influence on first order CDW intensities than on main lattice reflections in the ZOLZ.Still, it causes an additional broadening of the spot profile and a shift of the k z -component with highest CDW intensity (dotted curves).h Fitted quality factors of the template matching described in the text as a function of the k x , k y -components of the normalised incident wave vector k 0 .The black circle indicates the fitting results, the tilt direction for the rocking curves in c and g is given by the solid black arrow.A similar virtual tilt series only shifted by 0.29°along the direction given by the dotted arrow yields the dotted rocking curves displayed in c.

Fig. S2 CDW
Fig. S2 CDW diffracted intensity for reconstruction of rocking curves.a Intensity of first-order NC (left) and IC (right) diffraction spots as a function of temporal delay for the data set shown in Fig.3in the main text.The curves are additionally sorted and binned by the k z -component of the respective reflections.The scattered intensity at the IC spot position before time-zero is subtracted for both types of curves as a measure of the inelastic background in the images.The NC curves are additionally normalised with respect to the measured intensity before time-zero and show similar behaviour as a function of temporal delay (cf.Fig.3b).For the IC curves, we normalise to the measured intensity after 300 ps.During the phase formation, IC spots probed at a k z -component closer to the ZOLZ display an initial overshoot (blue curve).In contrast, reflections in the FOLZ are characterised by a continuous intensity increase (red curve), leading to the sharpening of the reciprocal lattice rod along the out-of-plane direction described in the main text and Fig.3.b Average intensity of first-order NC (blue), first-order IC (red) and second-order NC spots (green).Curve normalisation as described in a.

Fig. S3
Fig. S3 Ultrafast nanobeam diffraction of CDW dynamics.a Intensity of first-order NC (top) and IC (bottom) diffraction spots as a function of temporal delay for varying pump fluence.The data was recorded under tilted-sample conditions with an electron beam size of 170 nm and a repetition rate of 610 kHz.The curves are normalised as in Fig.S2.We observe an initial suppression within 500 fs that scales with the pump fluence.Above a threshold of around 2 mJ/cm 2 , the measured NC intensity after time-zero equals the inelastically scattered background.For such higher fluences, we measure an increase of IC spot intensity within 5 ps with a temporal characteristic independent of the pump fluence.b Delay scans of second-order NC intensity for varying pump fluence with the sample tilted into the zone-axis.In contrast to a, the electron beam has a larger diameter (1.5 µm), adding a spatiotemporal component to the dynamics following the initial suppression.For intermediate fluences, we witness a partial recovery of NC intensity after around 5 ps, followed by a second suppression at later delays.We attribute this behaviour to a locally varying excitation density leading to the formation and subsequent growth of IC phase domains on longer timescales as described in ref.11  .

Fig. S4 CDW
Fig. S4 CDW formation in time-dependent Ginzburg-Landau simulations.CDW pattern ρ l (see main text) at variable temporal delay derived from the time-dependent Ginzburg-Landau simulations.Orange and black circles highlight several disclinations with seven-fold and five-fold coordination, respectively (2 ps and 5 ps temporal delay).The late stage dynamics is governed by the annihilation of dislocations (black markers at 10 ps and 30 ps temporal delay).

Table I Parameters used in the diffraction simulation.
The atomic scattering factors f of tantalum and sulphur to calculate the structure factor for scattering by a wave vector with magnitude k are derived following f = i a i exp −b i k 2 .The coefficients apply to k < 6/ Å.The structure factor is additionally modulated by the Debye-Waller factor exp −2W k 2 , the specified values for W apply to a broad temperature range around ≈ 300 K.