Dual vortex charge order in a metastable state created by an ultrafast topological transition in 1T-TaS$_{2}$

Many body systems in complex materials undergoing non-equilibrium phase transitions may self-organise into ordered metastable emergent states with new and unexpected functionalities. Here, using large-area scanning tunneling microscopy we reveal an intricate chiral vortex structure and complex tiling of charged electron domains in the metastable metallic state in 1T-TaS$_2$ created by a non-equilibrium topological transition initiated by a single femtosecond optical pulse. A Moir\'e analysis shows that the interference of non-equilibrium nested Fermi surface (FS) electrons leads to the creation of charge vortices $\vec{\cal{D}}$ on a length scale of $\sim$70 nm. On a much smaller scale of $\sim5$ nm, domain configurational patterns appear, which show bound vortex-antivortex pairs, discommensurations, domain wall (DW) crossings and DW kinks, consistent with a rapidly quenched Berezinskii-Kosterlitz-Thouless (BKT) - transition. Revealing the detailed mechanism for the transition leads the way to design of long-range charge-ordered metastable states with intricate emergent properties under controlled non-equilibrium conditions.

Photoexcited metastable states in crystals have lifetimes which usually range from picoseconds to microseconds [1][2][3][4][5] , which makes it hard to investigate their structure in sufficient detail to reveal transient mesoscopic or microscopic order 8 . Consequently, the detailed mechanisms for photoinduced metastability have not been very well understood till now.
Even though it is of fundamental importance to prove the principle of the existence of LRO created of a transient emergent state, no direct experimental evidence has so far been presented in any system exhibiting a photoinduced phase transition. Uniquely, as a model system, 1T-TaS 2 (TDS) has a metastable state with a temperature-tunable lifetime, which is usefully long at low temperatures 6,9 for detailed experimental investigations. This opens the possibility of studying not only the fine structure, but also the origin and mechanism of the metastability with high resolution scanning tunnelling microscopy (STM). In this layered di-chalcogenide system (Fig. 1a) the competition of lattice strain, a FS instability and Coulomb interactions lead to a number of different phases and a multidimensional multi-parameter phase diagram, where strain, chemical pressure, doping and laser photoexcitation fluence play different roles [10][11][12] . In equilibrium, the high-temperature metallic state of TDS is unstable towards a three-pronged quasi-two dimensional FS nesting instability 13 causing a transition at T IC = 540 K from a single-band metal to a uniform three-directional incommensurate (IC) charge density wave (CDW) state. On cooling this transforms into a 'nearly-commensurate' (NC) state at T N C = 350 K as a result of incommensuration strain, giving almost commensurate patches separated by smooth discommensurations (DCs) 14 .
Eventually the state becomes fully commensurate (C) and insulating below 190 ∼ 220 K, as shown schematically 14,15 in Fig. 1b. The localized electrons in the C state have a remarkable star-of-David (SD) 'polaron lattice' structure, (Fig. 1c) where exactly one electron is localised on the central Ta of the SD and the 12 surrounding Ta atoms are symmetrically displaced towards it. Photoexcitation of the material by a single ultrafast optical pulse was shown to cause an insulator-to-metal transition (IMT) to a hidden metastable (H) state within a narrow range of photoexcitation fluences and pulse lengths. An indication of possible long range order (LRO) came from the narrowness of the collective amplitude mode of the CDW whose frequency shifts by 3% at the transition 6 . IMT switching caused by direct charge injection through electrodes was recently demonstrated in TDS [16][17][18][19][20][21][22] while switching by STM tip revealed that the uniform CCDW can be locally broken up into non-periodic domains. No LRO was observed however, in either case 18,19 . A hint of new states was found in ultrafast electron diffraction (UED) experiments 22 at significantly larger photoexcitation densities, but the k-space resolution was insufficient to resolve the ordering vector or the domain structure, so the fundamental question whether metastable LRO can form under non-equilibrium conditions remains unanswered. Here, for the first time we reveal LRO in an transient emergent state in a non-thermal transition using in-situ ultrafast optical switching in combination with a low-temperature high resolution STM. The studies reveal a remarkable dual vortex structure of polaron ordering in the emergent state, quite unlike any of the other states in this system, or in any other known material 23 .

Results
We photoexcite freshly in-situ cleaved single crystals of TDS with focused 50 fs single pulses at 400 nm within a UHV STM chamber. A scanning electron microscope (SEM) is used to select a homogeneous region of the crystal without strain and for precise positioning the tip and the beam on the sample (Fig. 1e). The pulse fluence was adjusted to just above switching threshold at 0.9 mJ/cm 2 6 , taking care to avoid heating the lattice to the    In the C state, CDW wavevectors ±Q (i) C appear as the six brightest peaks in the FT in Fig. 3a. Less intense atomic peaks coincide exactly with the linear combination of the CDW vectors, a * = 3Q (1) C , defining the commensurability of the C CDW with the underlying lattice. The NC state has additional domain structure, giving rise to peak splitting. Fig. 3b shows one bright CDW superlattice peak surrounded by several less intense satellite peaks. The bright peak corresponds to fundamental CDW wavevectors Q N C , and in the NC state the order is limited to n = 1 14 . The interference of fundamental and the first harmonic wave results in "beatings" -large-scale periodic modulations of the distortion amplitude. The interference pattern has hexagonal domains with inverse period k sat separated by smooth domain walls. The CDW is locally commensurate with the lattice inside the domain, and its only distinction from the C-state is a decreased CDW amplitude near the domain walls (see SI).
In the photoexcited state, we first ascertain the relation between the atomic lattice peaks and the CDW lattice peaks from a FT of atomic resolution STM image of several domains 26 (Fig. 3d). From the position of the second order CDW peaks (for higher accuracy) we obtain In the H state the single satellite characteristic of the NC state is converted into a diffuse streak spanning a range of angles 1.7 • < Θ < 9.5 • (Fig. 3e,f) which reflects the distribution of the domain sizes and domain shapes. The cross-section AA along the streak (Fig. 3h) reveals that FT intensity is not random in k-space, but rather there are approximately equally separated peaks with different intensities, showing that some domain sizes are more favourable than others. If we zoom into aggregates of the small or large domains, individual peaks can be discerned (see SI for more detailed analysis). The fundamental vector Q (i) H stays the same, independent of the domain size, as a consequence of the LRO.
The existence of different domain sizes can be qualitatively described by extending the harmonic picture. The first peak position in the streak is at Θ ≈ 1.7 • and fits rather well the expected position of Q sat with n = 1, while peaks at larger Θ appear as higher harmonics.
Their positions can be fit by varying angle φ and length Q H /a * , which gives φ = 13.51 • and Q H /a * = 0.278 for n = 1 . . . 5 (see SI for fitting and real-space reconstruction). We note that these values are close to the first metastable state observed under optical excitation in time-resolved UED data 22 . There it was tentatively associated with the triclinic CDW appearing in the equilibrium phase diagram, but our higher resolution allows us to associate it with the emergent photoinduced state (see Supplementary Fig. 12 for comparison).
We continue by examining the charge displacements within domains and the associated DWs. To characterize the topological structure, we define a misfit displacement vector D as the charge displacement in the H state with respect to the C lattice (see Fig. 4a for definition).
Locally, the structure is composed of irregular domains with different D (Fig. 4b-d). There are 12 possible discrete displacement vectors and D = 0 (Fig. 4a), which can be separated by 6 different types of DWs. Commonly, three DWs meet at a Z 3 vertex resulting in pairs of vortices and antivortices (labeled Y andȲ ) such as shown in Fig. 4b 23 . For each individual Y defect the misfit displacement vector sum is nonzero, e.g. D 0 + D 2 + D 6 = 0, but for each pair it must be zero: D 0 + D 2 + D 6 + D 5 = 0. The latter condition sets the local topological rules for the misfit displacement vectors. In addition, a significant number of X wall crossings (Fig. 4c) and kinks (Ks) (Fig. 4d) are observed, higher energy non-trivial defects created by photoexcitation which are not expected in equilibrium hexagonal tiling ( Fig. 4b). A contour encircling an X defect results in a non-zero sum of the misfit vectors, D 10 + D 2 + D 5 = 0, which makes this a non-trivial topologically protected defect. The K defect corresponds to a change of domain wall type without a change of the misfit field. It contains one incomplete polaron hexagram, a nontrivial defect with zero Burgers vector (see Suppl. Fig. 18).
In addition to the local vortex structure, there exists a long wavelength structure which is even more intriguing (Fig. 5). Plotting the magnitude | D| and angle α of the vector field D in Fig. 5b and c respectively as a Moiré interference pattern between the C lattice and the H lattice (see Fig. 5a), we observe a clear hexagonal lattice of vortices with a period L H = 2π/(|Q H − Q C |) ∼ 70 nm (Fig. 5b-d). Repeated experiments show that the winding number 29 of the large-scale vortices is either n w = −1 or n w = +1, but never both simultaneously, implying that there is a macroscopically spontaneously broken mirror symmetry in the formation of the LRO D vector field. We also observe an accompanying weak, but unambiguous modulation of the domain size with a similar wavelength shown in

Discussion
The mechanism for the transition is quite unlike previously studied non-equilibrium second-order symmetry-breaking FS nesting transitions such as the one in TbTe 3 32 , and the trajectory of the C→H transition cannot be discussed in terms of conventional nonequilibrium symmetry-breaking with a Landau order parameter that vanishes at some critical time t → t c 32 . It requires consideration of the competing effects of FS hot-electron interference, Coulomb interactions and particularly the incommensurability strain which force the ordering of the system on different timescales. The state is created under very specific conditions of laser fluence 0.8 < F < 2.5 mJ/cm 2 , within a time and temperature window defined by the pulse length range 30 fs < τ p < 5 ps 6 . During this time, the lattice remains below the transition temperature to the NC state, while the peak electronic temperature is estimated ∼ 1700 K 6 , thermalizing with the lattice on a timescale of ∼ 7 picoseconds to 150 ∼ 185 K (see Fig. 5 in SI for lattice temperature estimates). Transient photodoping is an essential component: photoemission experiments show that short pulse photoexcitation causes melting and a rapid < 50 fs transient shift of the chemical potential above the switching threshold µ c 33,34 , which implies that FS nesting is modified while the system is out of equilibrium. The hot nested electrons interfere to condense into a transient incommensurate charge modulated structure with a wavevector Q * IC , reflecting the non-equilibrium FS. The accompanying change in low-energy electronic structure occurs on a timescale < 400 fs, as recently shown by coherent phonon spectroscopy 35 . Thereafter, as the incommensurate electronic CDW tries to adjust to the nearest commensurate lattice Q C , individual domains spontaneously form (Fig. 4e). However, local topological defects such as Ks, and Y vortices without associated nearby anti-defects 29,36-38 also appear, presumably as a consequence of the rapid quench through the transition. On a longer timescales 40,41 , a continuous tiling emerges with a LRO Q H , and the K and Y defects become bound to form X The validation of the concept of formation of emergent LRO through many-body interactions under non-equilibrium conditions, and the underlying topological transformation mechanism revealed by these experiments represents a large step towards understanding, and whence designing new emergent metastable states in complex materials. The present system, with its associated M-I transition may lead to advances in novel non-volatile all-electronic memory devices without the involvement of ions or magnetism, through controllable switching between electronic topologically ordered states.

Methods
The results presented here were measured in the custom-built Omicron Nanoprobe 4probe UHV LT STM with base temperature of 4.2 K and optical access for laser photoexcitation (Fig. 1e). Crystals of 1T-TaS 2 were synthesised by the iodine vapor transport method.
Samples were cleaved in situ at UHV conditions and slowly cooled down to 4.2 K.
Photoexcitation of 1T-TaS 2 single crystals was performed at 4.2 K (C-state) with a single 50 fs optical pulse at ∼ 400 nm (second-harmonic generation from 800 nm Ti-Sapphire laser) focused onto a 100 µm diameter spot within STM UHV chamber. This ensures highly spatially homogeneous excitation on the scale of the STM scans (∼ 200 nm). It is important that photoexcitation energy density is carefully adjusted to be slightly above the threshold value (∼ 0.9 mJ/cm 2 6 ) for switching to the hidden state. Significantly higher excitation energies result in heating of the lattice above the phase transition temperature T N C−C 6 , which we want to avoid.
The surface was characterised each time before photoexcitation to confirm a defect-free initial C state. The beam was aligned at low power to hit the apex of the STM tip. Then the tip was retracted and a high-power single pulse was applied (Fig. 1e). After approaching the tip, an area of the order of 500×500 nm 2 was checked for homogeneity.
Contrast adjustment in Fourier transforms is routinely used for analysis of STM data in       N C peak (red), surrounded by 3 neighbours (blue). The domain modulation vector is given by k sat . (see Supplementary Fig. 3 for more details).
Note that the FT amplitude is zero between the satellite peaks. with atomic lattice vectors, whereas in NC and C states they are shrunk and rotated. One can thus describe them as the following product: where a * , b * are reciprocal vectors of Ta triangular atomic lattice, γ = Q (i) /a * is star of David size and φ is the rotation angle between CDW and atomic lattices.
In CCDW state the fundamental Q (i) C is determined by the commensurability condition: or, equivalently,

3Q
(1) The above commensurability condition allows to calculate the ideal length of CDW vector, C /a * ≈ 0.277, and its rotation w.r.t. atomic lattice, φ ≈ 13.9 • . Real space images of IC and C state can be readily built as the following sum: where R is in-plane real-space vector. The sum term gives CDW modulation, the last two -atomic modulation. Suppl. Fig. 1a-d shows the real and reciprocal space models of IC and C states. In the commensurate case each star of David is centered at a certain Ta atom position, resulting in the same atomic structure of all CDW units across the sample (Suppl. Fig. 1c). This feature can be checked experimentally in atomic resolution STM images and is indeed observed in C, NC and H states (main text, Fig. 2b,d,g). In contrast, violation of commensurability will result in irregular atomic configuration for each star of David, as illustrated in Suppl. Fig. 1a. In NCCDW state the periodic domain structure appears, a model example of which is shown in Suppl. Fig. 2a N C basis of CDW. It can be modelled with simple approach used to interpret the early X-ray results 2 , which considers domain structure as a result of the interference of two triple CDWs. One of them is given by the three fundamental NCCDW vectors, whereas the other has smaller amplitude and is given by three additional vectors chosen in such way, that their linear combination with fundamental ones will be equal to some vector of reciprocal atomic lattice. This gives the new commensurability condition: Eq. 2 The corresponding vectors are shown in Figure 2b,c. This equation gives only one possible commensurability condition for reciprocal lattice vector G = a * , but it can be extended to the case of so-called higher harmonics G = na * , n = 1, 2, . . ., as shown in the main text.
It is worth noting, that the domain wall is always two stars of David thick, independent of Q (i) N C , as long as (Eq. 2) is fulfilled. This model also gives slowly varying modulation of height, peaking at the domain centre (cf. Suppl. Fig. 2a). Real space image is then given by: Full details of the nearly commensurate state are captured with Nakanishi-Shiba theory 1 , which considers NC state as a product of CCDW and periodic domain structure overlayed on top of it: C are fundamental CCDW wavevectors and is the modulation responsible for the domain structure.
Modulation periods are given by: Eq. 3 Here the last term is the discommensuration vector: N C being the fundamental NCCDW wave vector (see Suppl. Fig. 2d). q (i) therefore represent the deviation from commensurability. The correct domain periodicity should result in the commensurate CDW inside the domains and thus can be calculated using commensurability condition similar to that for fundamental vectors (Eq. 1): Eq. 5 The above equations are illustrated in Figure 2d.
These two approaches are equivalent and connected by the relation (Suppl. Fig. 2c): sat . Eq. 6 Finally, additional satellites can appear from linear combinations of Q sat where i = j. This is equivalent to several of (l, m, n) coefficients being non-zero in NS model (Eq. 3). For example, the satellite −k (3) domain will be seen near Q Eq. 7 Below we will refer to these satellites as −k   . Fig. 4b). The results are shown in Suppl. Fig. 4a, where the first three panels correspond to David stars PDF and the last one is atomic. In the latter, peaks are marked by the respective atomic configurations, i.e. 2 + 1 is the sum 2 a + b. In the case where just one satellite, k      Ta atoms (c) used for PDF calculations.

Supplementary note 2: fluence dependence of amplitude mode temperature
The peak lattice temperature is determined from the frequency of the amplitude mode (AM) (Suppl. Fig. 5a), whose temperature dependence is determined from independent low-fluence measurements (Suppl. Fig. 5b). A plot of T AM vs. fluence (Suppl. Fig. 5c) shows that at threshold, the temperature of the AM reaches 150 ± 10 K. The AM is the mode which is most strongly coupled to the electrons, so its temperature is inevitably by far the highest of all the lattice modes. In order to get high value of satellite spread angle Θ and small deviation of the CDW rotation angle φ from its C-state value within the NS model, we use k (i) domain harmonics (see also Suppl. Fig. 16d). The fitting procedure of the real data with NS model is given below.
In the first step we estimate the required number of harmonics. To this end, we find the smallest angular distance between the fundamental, Q H , peak and the satellite peaks around, which in this case is Θ 1 ≈ 1.7 • . This value appears the same, independent of the fundamental peak (i.e. i = 2, 3) It gives a rough estimate of k (1) domain length and the number of groups as Θ max /Θ 1 = 5. Indeed, one can see from the cross-section that the peaks along k (1) domain direction can be separated into five groups (Suppl. Fig. 8).
In the second step, we determine the averaged peak positions for each group. Peaks along the k (i) domain direction in the experimental data (Suppl. Fig. 9a,b) are additionally split in the radial direction by 2-3 pixels. This gives the characteristic scale of the additional modulation equal to 1/3 of the scan size or ∼ 70 nm. The only periodicity that can be found on such scale in real space image corresponds to the distance between groups of small domains, separated by groups of larger domains and vice versa. This kind of modulation cannot be described by NS model. We collapse them into 'average' single peak using the simple centroid procedure described below.
To this end we estimate seed positions, that should be located at multiples of k  domain . The result of the fit is shown in Suppl. Fig. 9b. The comparison of the STM data with the real-space simulation of the individual harmonics is presented in the Suppl. Fig. 9c-f and shows nice correspondence of the domain sizes.
The emergence of higher harmonics underlines the differences between hidden and equilibrium states. The distinction that takes place in the reciprocal space can be illustrated on the three-dimensional phase diagram, which shows simultaneously the fundamental CDW vector position, (φ, Q (i) /a * ), and corresponding harmonics position, Θ (Suppl. Fig. 10).
It appears, that with temperature and photodoping 1T- TaS  The height profile (Suppl. Fig. 11c) confirms that David star distortion is the same independent of the distance to the domain wall (x = 7 nm). Next, we apply band pass filter to the Fourier transform of the image: the band is a ring including only the first CDW Brillouin zone and all the Nakanishi-Shiba harmonics in it. Suppl. Fig. 11b shows the same two domains as before but in the bandpassed inverse Fourier transform. The respective height profile (Suppl. Fig. 11c) clearly shows, that David star distortion becomes smaller towards the domain wall. Thus, the harmonics observed in the H state are not related to the domain flatnees. This result also shows that the latter comes from the increased intensity in the higher Brillouin zones of CDW. Switching from CCDW to metastable metallic hidden state can be done also with electrical pulse 9 applied between two contacts on top of 1T-TaS 2 crystal. Alternatively, the pulse can be applied between the crystal and STM tip 4,10,11 , and the resulting state can be imaged in situ. Below we perform the same analysis of electrically switched state as the one presented in the main text for the optically switched state.
Here the best switching results were obtained for small tip-sample separation with a low pulse voltage (V p ∼ 3.5V ). Suppl. Fig. 14a shows that the homogeneous CCDW state To illustrate the real space details of the domain structure, we have performed the mapping of David stars misfit vectors. The results are shown in Suppl. Fig. 15. Misfit vector field shows no distinct structure. Its amplitude changes randomly in space (Suppl. Fig. 15b), and its direction jumps between very different, sometimes opposite, values (Suppl. Fig. 15b). Such a behavior is qualitatively different from the Moiré-like vortex structure observed for the long-range ordered H or NC states (see main text): gradual rotation of the direction (dark-blue to dark-red) and gradual increase of the amplitude from 0 (green) to 1 (cyan) and then 1 + 1 (blue). This allows us to conclude that the H state obtained by electrical switching does not have long range order.
It should be noted, that domain walls in this state also show X and K crossings in addition to the standard Y −Ȳ ones characteristic of the pure hexagonal structure. This resemblance could suggest similar mechanism of transition. Since NS theory describes well the domain structure of NCCDW state, we can use it to model various outcomes that could be found for supercooled NCCDW state. Here we consider two relevant scenarios: (i) "frozen" domains with single Q N C and (ii) sum of "frozen" domains corresponding to different Q N C .
The first scenario is equivalent to the trajectory, where the system is heated to some temperature above C-NC transition and then is quenched to low temperatures. The model FT corresponding to this scenario is shown in Suppl. Fig. 16a, where the parameters are chosen to obtain large Θ angle (length of the streak) similar to that in experiment. This corresponds to the temperatures as high as T = 300 K, since at lower temperatures Θ is known to be smaller. In FT picture we should observe single fundamental Q This situation is discussed below.
The second scenario corresponds to the trajectory, where high-temperature state is quenched at intermediate rates and different Q N C are present in real-space picture. Another option is that single-Q rapidly quenched state has relaxed at low temperatures, as described above. In both trajectories the FT of such image will result in a number of Q  Fig. 16c). This scenario has two distinct features in the FT picture. First, fundamental peak will be smeared and satellites will no longer be equidistant (due to random choice of local fundamental vectors). Second, satellites will not lie on the same line with fundamental peak. Both these features are in contrast to experimentally observed FT picture. For comparison the model CDW state with 5 harmonics used to fit experimental data for the hidden state is shown in Suppl. Fig. 16d   domain . All of them correspond to either of three lines, which originate from fundamental peak (compare to panel (c)).

Supplementary Method 2: mapping misfit vectors in the H state
David stars positions in the H state can be determined with high accuracy using image processing methods. Each David star appears as a well-defined peak in STM images. The accuracy of finding the peak could be compromised mostly by the noise. The noise spatial scale can be divided into two categories: (i) slowly varying background on the scale of several DS and (ii) spikes smaller than DS. Type-1 noise makes detection of DS inside domain walls hard. This problem is solved by using gradient-based algorithms, e.g. Laplacian of Gaussian.
Type-2 noise can be overcome using the known separation of neighboring DS. Still algorithm performance can be improved by preprocessing the image, based on the known global and local structure.
Both noise sources can be most substantially reduced in the bandpassed inverse Fourier transform of an experimental STM image. To this end, the band is selected to include only the fundamental and harmonic peaks. It is important that such procedure preserves the whole domain wall structure (see Suppl. Fig. 17a,b for comparison). Another approach is to build cross-correlated image using the simulated David star shape and template matching algorithm. We have checked that different preprocessing methods result only in quantitative differences of several pixels in peak position, whereas qualitative behavior is the same. Below we use the maximum filter based on the dilation algorithm 13 applied to bandpassed inverse FT image (Suppl. Fig. 17c).
To map the misfit vectors of David stars in the H state with respect to their original positions in the C state, we have to know the latter. Here we note, that C order is locally present in the single domain. Therefore, we create the simulated C lattice which exactly coincides with that in the largest domain we can find in the experimental image (Suppl. Fig. 17d). In this way two arrays containing H and C state centers are obtained. Misfit vectors D are then calculated by finding the nearest C neighbor for each H state David star and measuring length and angle of the shift between them (Suppl. Fig. 17d,e). Angle is then corrected by the difference between atomic axis and experimental x axis. Each H state David star is then given two discrete indices: one shows the shift on (0, 1, 1 + 1) atomic scale and another -the angle on the 12-category scale (see Suppl. Fig. 17f,g). To assign the shift index we assume that detected position can have the error up to ±2 pixels. In the next step angle index is assigned based on angle and the known shift index. Here we assume that angle can have error up to ±30 degrees. STM fine X-Y piezoscanners were calibrated at low temperatures within ten percents before measurements. Later each picture was corrected with FFT peaks to match the triangular lattice of stars of David, i.e. using 2 fundamental CDW vectors. Correction for scanner calibration, drift and X-Y crosstalk was done either by calculating the affine transform between the observed and ideal lattice. In all the cases 1.174nm was taken for CDW period, but this choice affects neither ratio between CDW and atomic periods nor the angles. The ratios were checked for consistency with positions of atomic peaks in FT whenever possible. The ratios were checked for all -H, C, NC -states.