# Crystallization of polarons through charge and spin ordering transitions in 1T-TaS2

E. Bozin, M. Abeykoon, S. Conradson, Gianguido Baldinozzi, P. Sutar, D.

Mihailovic

► To cite this version:

E. Bozin, M. Abeykoon, S. Conradson, Gianguido Baldinozzi, P. Sutar, et al., Crystallization of polarons through charge and spin ordering transitions in 1T-TaS2. Nature Communications, 2023, 14 (1), pp.7055. 10.1038/s41467-023-42631-6. hal-04269727

HAL Id: hal-04269727

https://cnrs.hal.science/hal-04269727v1

Submitted on 3 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
1 Crystallization of polarons through charge and spin ordering transitions in 1T-TaS₂

2

3 E.S. Bozin¹, M. Abeykoon², S. Conradson³, G. Baldinozzi⁴, P. Sutar³ and D. Mihailovic³

4 ¹Condensed Matter Physics and Materials Science Division, Brookhaven National Laboratory, Upton, NY 11973,

5 USA

6 ²Photon Sciences Division, Brookhaven National Laboratory, Upton, NY 11973, USA

7 ³Dept. of Complex Matter, Jozef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia

8 ⁴Centralesupélec, CNRS, SPMS, Université Paris-Saclay, bât Eiffel, Gif-sur-Yvette, Île-de-France, 91190, FRANCE

9

10 ABSTRACT

11 The interaction of electrons with the lattice in metals can lead to reduction of their kinetic energy to the point

12 where they may form heavy, dressed quasiparticles – polarons. Unfortunately, polaronic lattice distortions are

13 difficult to distinguish from more conventional charge- and spin- ordering phenomena at low temperatures.

14 Here we present a study of local symmetry breaking of the lattice structure on the picosecond timescale in the

15 prototype layered dichalcogenide Mott insulator 1T-TaS₂ using X-ray pair-distribution function measurements.

16 We clearly identify symmetry-breaking polaronic lattice distortions at temperatures well above the ordered

17 phases, and record the evolution of broken symmetry states from 915 K to 15 K. The data imply that charge

18 ordering is driven by polaron crystallization into a Wigner crystal-like state, rather than Fermi surface nesting or

19 conventional electron-phonon coupling. At intermediate temperatures the local lattice distortions are found to

20 be consistent with a quantum spin liquid state.

21

22
23 INTRODUCTION

24 Quasi-two dimensional quantum-ordered materials display very diverse and often exotic spin<sup>1-3</sup> and charge
25 ordering behaviour<sup>4</sup>, including exciton condensation<sup>5</sup>, chiral spin liquids<sup>6</sup> and superconductivity<sup>7,8</sup>. They can also
26 display hidden phases which can be manipulated by pressure<sup>9,10</sup>, light<sup>11-15</sup> or doping<sup>16-18</sup>. An important aspect
27 that determines spin ordering behaviour is the competition between Coulomb repulsion and interlayer
28 interactions<sup>19</sup>. A common feature of these systems is that the strength of the electron phonon interaction (EPI)
29 and electron density control the electron kinetic energy. Increasing the EPI leads to bandwidth narrowing, and
30 lattice dressed electron quasiparticles - as originally introduced by Landau in 1933<sup>20</sup>. A further consequence of
31 electron localization is an enhanced Coulomb interaction that leads to correlations, ubiquitous mesoscopic
32 textures<sup>21,22</sup>, charge-density wave (CDW) formation<sup>23</sup> and superconductivity at high critical temperatures in the
33 fluctuation-dominated intermediate coupling regime of layered cuprate superconductors<sup>24</sup>. Recent model
34 calculations<sup>23</sup> applied to wide range of transition metal dichalcogenide (TMD) materials and comparisons with
35 scanning tunneling microscopy experiments highlighted the importance of the Coulomb interaction that lead to
36 ubiquitous charge ordering at magic filling fractions<sup>23</sup>. Specific features of the band structure may lead to further
37 effects, such as exciton condensation<sup>5</sup> and CDW ordering aided by Fermi surface nesting<sup>25</sup> on top of the
38 seemingly universal charge ordering<sup>23</sup>. Experimental evidence for polarons comes from many different
39 techniques in a very wide variety of materials<sup>26,27</sup>, but their role in forming quantum states, such as Wigner
40 crystals, superconductivity and quantum spin liquids is not experimentally clear, particularly on the transition
41 from the low to high density limit. Direct observation of symmetry-specific polaronic lattice distortions on the
42 low- to high- density limit are crucial for progress in understanding a wide class of quantum materials.

43 As a prototypical example, the layered dichalcogenide 1T-TaS₂ in equilibrium is of interest as a Mott insulator²⁸,
44 CDW system⁴, polaronic superconductor⁹ and quantum spin liquid². Additional recent interest in the material
45 comes from the discovery of hidden light- or charge injection- induced mesoscopic metastable phases ¹²,²⁹ that
46 are important both fundamentally and for applications³⁰⁻³². Structurally, it is a kinetically trapped 1T polymorph
47 of TaS₂, with octahedral stacking (Fig 1a) ⁴,³³⁻³⁵. It's deceptively simple electronic band structure in the
48 undistorted phase has a single Ta metal band crossing the Fermi level.¹²,³⁴ While nominally metallic at high
49 temperatures, it shows an incommensurate (IC) CDW (ICCDW) along the principal crystal axes already at high
50 temperatures (between 550 and 350 K), with an associated Kohn anomaly at a wavevector $q_{IC} \approx (0.283, 0,0)$
51 in units of the reciprocal lattice vector $a_0^{36,37}$ which implies a strong coupling between electrons and the lattice.
52 Below 350 K, the IC state breaks up into a patchy periodic network of commensurate islands separated by
53 discommensurations³⁸⁻⁴⁰. Upon further cooling, this nearly commensurate (NC) CDW (NCCDW) state becomes
54 fully commensurate (C) CDW (CCDW) after a first order transition at ~ 180 K, with the appearance of a Mott
55 gap in the charge excitation spectrum observed by single-particle tunneling⁴¹ and angle-resolved
56 photoemission²⁵,⁴². In contrast, gapless spin relaxation in the range 50-200 K is attributed to a possible quantum
57 spin liquid (QSL) phase, arising from the spin frustration on the triangular superlattice. While band structure
58 calculations³ suggest that the QSL phase may be suppressed by inter-layer coupling⁴³, this cannot explain the
59 marked difference of spin and charge excitation spectra². A further intriguing feature is a discontinuity in the
60 magnetic susceptibility and an anomalous peak in the spin-spin relaxation time concurrent with a crossover
61 around ~ 50 K of the spin-lattice relaxation time from a QSL-like, to a gapped T-dependence⁴⁴,⁴⁵. Near 50 K, a
62 resistivity upturn from nearly T-independent resistivity to strongly T-dependent variable-range hopping behavior
63 has been commonly reported, suggesting an onset of low-temperature charge localization⁴,⁴⁶. The common
64 picture of the 1T-TaS₂ lattice structure in the C state assumes simple distortions from the average lattice
65 structure (trigonal P\bar{3}m, Fig. 1a), whereby 12 Ta atoms are attracted symmetrically towards the central Ta atom
66 that carries an extra electronic charge, forming a star of David (SoD) pattern (Fig.1c). The SoDs are arranged in
67 a √13 × √13 unit cell superlattice structure (P\bar{3} symmetry⁴⁷). At such an electron density of 1/13, the
68 dimensionless ratio $r_s$ of the Coulomb energy $V$ to the kinetic (Fermi) energy $E_K$ $r_s = \frac{V}{E_K} = \frac{e^2m}{\hbar^2\sqrt{n}} \approx 70 \sim 100$,
69 where $e$ is the elementary charge, $n$ is the (2D) density, and $m \sim 3$ is the effective electron mass²⁵. Considering
that the Wigner crystal stability limit is $r_s = 31 \sim 38$ $^{23,48-50}$, 1T-TaS$_2$ is comfortably in the Wigner crystal regime, which means that electronic superlattice ordering on the basis of dominant Coulomb interactions and tell-tale lattice deformations around localized carriers may be anticipated. Room temperature extended X-ray absorption fine structure (EXAFS)$^{51}$, high-resolution transmission electron microscopy (HRTEM)$^{52-54}$ and X-ray structural measurements$^{44}$ indeed suggest the existence of symmetry-breaking atomic displacements, but studies over a wide range of ordering temperatures have not yet been performed.

Here we present the first systematic measurements of the local lattice structure that covers temperatures from 15 to 915 K. The data reveal local symmetry breaking from the established P$\overline{3}$m symmetry at all measured temperatures (Fig. 1a, inset), and unambiguous evidence for individual polarons in the dilute high-temperature limit, revising our current notions about the origin of charge and spin ordering, with clear implications for other materials of current interest.

## 81 RESULTS

The local structure of 1T-TaS$_2$ is investigated using X-ray pair distribution function (PDF) analysis$^{55,56}$, giving picosecond-timescale structural snapshots for $15 \text{ K} \le T \le 915 \text{ K}$. This encompasses the known electronically ordered states below 600 K portrayed by the resistivity vs. temperature plot shown in Fig. 1a$^9$, as well as irreversible polytype transformations above 600 K.

The anisotropic atomic displacement parameters, $U_{ij}$, of Ta atoms were retrieved from a fit to the established P$\overline{3}$m structure over a range of 20-60 Å (Fig. 1b). For $T < 600 \text{ K}$, the $U_{ij}$ reveal an anomalously large in-plane component, $U_{11}$, approximately a factor of 2 larger than the component along the stacking axis, $U_{33}$ (e.g. at 500 K $U_{11}$ is 0.025 Å$^2$ whereas $U_{33}$ is 0.01 Å$^2$). This is unusual for a layered system in which a larger out-of-plane component is expected to arise from interlayer disorder. At high $T$, in the metallic (M) phase, the large $U_{11}$ signifies an apparent intralayer 'disorder', unaccounted for by the P$\overline{3}$m structural model. At lower temperature, the P$\overline{3}$m model is inadequate at any temperature, revealed by the difference PDFs shown in Fig. 1c. The spatial extent and character of this symmetry breaking changes across the cascade of CDW transitions, tracking complex correlations in different electronic regimes. Intriguingly, as we describe below, nanoscale symmetry breaking is also evident in the polymorphs above 630 K (Fig. 2b,c).

As an aid to understanding, pair-specific contributions to the total PDF within the P$\overline{3}$m model are shown in Fig. 1e, where Ta-Ta pair contributions are dominant (see Methods). In the P$\overline{3}$m structure, the undistorted Ta sublattice features a single-valued nearest-neighbor Ta-Ta distance, seen as a sharp peak at ~3.4 Å, Fig. 1e. However, the data reveal a distinct additional shoulder, (e.g. Figs. 1e and 2b for 600 K), which is unexplained by the P$\overline{3}$m model (Figs. 1c and 2c (top)). At high temperature this distortion is most likely dynamic, as expected for symmetry-specific lattice fluctuations associated with carrier localization and polaron formation. Very specific distributions of Ta-Ta and Ta-S pairs in the 3.3-5.9 Å range are expected for SoD polaronic distortions$^{47}$. In the C phase with a single layer $\sqrt{13} \times \sqrt{13}$ supercell P$\overline{3}$ model (Methods and Supplementary Table 1) the SoD polaron structure is adequately described (Fig. 2f)$^{47,57}$. However, for all other electronic regimes it was necessary to add P$\overline{3}$m as a secondary phase to achieve good fits (see Fig. 3c, inset). The temperature dependence of the interatomic distances and bond angles, and deductions resulting from fits to the data are addressed sequentially below, starting from high temperatures, well above any known charge or spin ordered phase.

**Polymorphic regime** – Heating 1T-TaS$_2$ above 600 K results in irreversible polymorphic transformations altering average symmetry and stacking, $^{26,51}$where nominally 50% of octahedrally coordinated 1T layers convert into trigonal prismatically coordinated 1H layers $^{35,58}$ (Fig. 2a inset). This is accompanied by polymorph-specific stacking of the two layer-types $^{4,45}$ and changes of the local Ta environment in the prismatic layers. Two such transformations are seen in the PDF data (Fig. 2b), at 630 and 880 K (see Supplementary Figure 1 for more details). Remarkably, a shoulder-signal (~3.5 Å, Fig. 2b) associated with local lattice distortions (polarons) is
observable for 730 K < T < 915 K. The distorted fraction does not vanish in the polymorphs but instead drops to
half its value observed in the 1T regime (T < 630 K), with no significant reduction through the second
transformation (Fig. 2b). Notably, PDF data above 730 K conform to the 6R-TaS₂-like model (R3m symmetry)
featuring equally abundant alternating 1T and 1H layers (Fig. 2c, bottom)⁵⁹. The model features undistorted 1T
and 1H layers, but a clear distortion signal, albeit weaker, such as is seen at 600 K, persists in the fit differential,
implying the existence of polaron-specific deformations in the 1T layers well above 600K (Fig. 2c, top), with ~10
wt% distorted fraction (see below).

**High-temperature metallic regime** – the data at 600 K show multiple components at the position where the P3̄m model predicts a single Ta-Ta peak (Fig. 1c and red trace fit in Fig. 2c, R<sub>w</sub>=15%). The distortion spans ~5 Å with weaker signal extending up to ~12 Å, depicted by the difference PDFs. Such signatures, also seen in the IC phase, are weaker than in the NC and C regimes (Fig. 1c). However, their similarity and continuity imply a common origin and unambiguously reveals the presence of high temperature electron localization, detectable by charge-lattice coupling<sup>60,61</sup>. Despite similar P3̄m model misfits shown in Fig. 1c, explicit data differences, Fig. 1d, reveal subtle changes across the M-IC transition for r > 5 Å. The electron localization results in incomplete breakup of the P3̄m matrix, but without the formation of the full SoD motif. The light red profile in Fig. 2d and the intensity band in Fig. 1e at ~3.3 Å remain broad, consistent with incoherent distortions from weakly developed fluctuating charge correlations. In this diluted limit local distortions were approximated as P3̄ nanoscale inclusions, Fig. 2g, within the undistorted P3̄m matrix. This model (cyan traces in Fig. 2c (top), R<sub>w</sub>=8%) accounts for the misfit, resulting in ~20 wt % of P3̄ distortions involving heavily puckered hexagonal Ta discs surrounding central tantalum. (This is consistent with ~10 wt% of locally distorted environments in the 1T-layers in the polymorphic regime above 630 K discussed in the previous section).

**Incommensurate regime** – The growth of charge correlations and structural coherence of associated distortions
results in their signature extending over a longer length-scale (Fig. 1c). Stacks of experimental PDFs for 15-500
K (Fig. 2d,e) show that short length-scale peaks (r < 5 Å) sharpen gradually on cooling. However, the data still
does not show a clear SoD motif (Fig. 2d). The data at ~3.5 Å are broad and unresolved, and at ~3.75 Å are rather
featureless. Conversely, peaks at higher distances show the opposite temperature trend. In the light red stack,
Fig. 2e, the signal at ~23.5 Å is the sharpest just below the M-IC transition, portraying high symmetry, and
systematically broadens with intensity reduction as temperature decreases, evidencing growth of broken-
symmetry correlations. The change in slope of Te U₁₁(T) at 550 K (Fig. 1b) indicates an increased localized charge
density. Model-derived Ta-Ta distances for T > 350 K shown in Fig. 3a portray substantially distorted hexagonal
SoD cluster cores comprised of short Ta-Ta contacts (Fig. 3d,f). The distorted fraction increases to ~30 wt%
compared to the M state (Fig. 3c). The PDF peaks which are sensitive to interlayer stacking are very broad, as
shown in Fig. 4 a,c, implying inhomogeneous stacking. The distribution centroids of the closest interlayer
separation, Fig. 4a, are shifted to the left, well below 5.9 Å, indicating local compression of layers which are, on
average, closer together than they would be in the undistorted 1T structure. The observation of negative thermal
expansion jumps along the stacking axis is relevant in this context<sup>62</sup> (see Supplementary Figures 2 and 3 for
discussion).

**Nearly commensurate regime** – Significant local structure changes occur at the IC-NC transition. The distortion
fingerprint, Fig. 1c, grows dramatically and PDF intensities vividly redistribute below 5 Å (Figs. 1f and 2d). The
Ta-Ta peak at 3.25 Å, which is broad in the IC phase, splits into an incompletely resolved triplet of unequal
intensity peaks at ~3.15 Å, ~3.3 Å and ~3.8 Å, marking the formation of well-defined SoD clusters. The feature
at ~3.15 Å corresponds to the Ta₁₃ SoD intra-cluster configuration, whereas the other two features depict the
inter-cluster correlations (Fig. 3). Residual intensity at ~3.4 Å is associated with discommensuration regions of
approximately P3̄m character, separating hexagonal domains of well-ordered SoD clusters within the Ta
layers⁶³,⁶⁴, well observable by STM. Multiphase modeling with coexisting P3̄ and P3̄m phases confirms this
picture (Fig. 3c inset, see Methods). As SoD clusters emerge, the distortion magnitude increases (Fig. 3a,b),
consistent with increasingly resolved signal seen for 3 Å < r < 4 Å shown in Fig 2d. Importantly, the star vertices
161 exhibit sudden asymmetric contraction, accompanied by substantial puckering (Fig. 3d,f) and SoD twisting (Fig. 162 3e). The distorted P3̄ phase becomes the dominant fraction below 350 K, at ~60 wt% (Fig. 3c). This is consistent 163 with the patchy nature of the discommensurate NCCDW phase. Narrow thermal hysteresis, shown in Figs. 3 b- 164 e, reflects the 1st order character of the IC-NC transition. As SoD local order is established in the NC state, they 165 simultaneously form well defined bilayer correlations, Figs. 4 a,c: The SoDs in adjacent layers couple and align 166 along the c-axis. In contrast, inter-bilayer correlations are still fuzzy: different bilayers exhibit appreciable 167 stacking disorder, evident from inter-bilayer sensitive peak, Figs. 4 b,d, which probes off-axis correlations whose 168 distribution is still considerably broad in the NC state.

169 **Commensurate regime** – At the NC-C transition, the PDF intensity at ~3.4 Å diminishes, implying the gradual 170 disappearance of the ordered domain walls, while the intra-SoD-cluster peak sharpens (Fig. 2d), consistent with 171 increased SoD density and additional charge localization. Similar changes are observed in the far-r limit (Fig. 2e), 172 as a homogeneous crystal structure forms in the commensurate state. Closer inspection of the data unveils an 173 additional restructuring below ~50 K, embodied in subtle peak shifts and sharpening observed in the dark blue 174 traces in Figs. 2d and 4a. Modeling shows that the distribution of Ta-Ta distances depicting intra- and inter-SoD 175 Ta-Ta contacts changes character in the C phase, with tightly bound SoDs structurally further regularizing below 176 50 K (Fig. 3a,b). The intra-SoD angle θ reduces, approaching 120° at base temperature, Fig 3e. The intra-star 177 twist angle α increases, approaching 90° below 50 K (Fig. 3d), while the local inter-star angle β reaches 150° as 178 SoDs align and form fully commensurate order (Fig. 3f) and the undistorted phase diminishes, consistent with 179 vanishing domain walls (Fig 3c). Simultaneously, well defined inter-bilayer correlations develop (Fig 4 b,d). The 180 most dramatic changes are observed at ~13.1 Å where PDF intensity redistributes and sharpens, possibly 181 indicating an offset of pairs of SoDs, going from one bilayer to another, by one P3̄m lattice spacing along both 182 planar lattice vectors, Fig 4e. Notably, there are 6 choices of such offsets, consistent with 13x stacking period 183 and helical stacking arrangement suggested in some works⁴⁴. Further qualitative changes are seen below 50 K 184 in the PDF peak which includes the nearest neighbor intra-bilayer stacking peak just below 6 Å, whose intensity 185 distinctly and abruptly shifts to lower distance (Figs. 4a,c). In contrast, the inter-bilayer peak at ~13.1 Å sharpens, 186 but does not shift (Figs. 4b,d). This could imply enhanced intra-bilayer coupling, with neighboring bilayers 187 presumably adjusting to preserve the inter-bilayer spacings.³⁶

188 **Symmetry of the order parameter.** On the basis of a symmetry analysis of the observed PDF displacements we 189 can identify the symmetry of order parameters which transform as the $A_{2u}$ and $E_u$ irreducible representations 190 in the parent P3̄m structure space group. The displacements are described by two modes with wavevectors of 191 the form $q_1 = (a, b, 0)$ and $q_2 = (-b, a + b, 0)$ at 60° to $q_1$. In the C phase, $q_1 = \frac{3}{13}a^* + \frac{1}{13}b^*$, and $q_2 = 192 -\frac{1}{13}a^* + \frac{4}{13}b^*$, where $a^*$ and $b^*$ are the reciprocal lattice vectors of the undistorted lattice. It is possible to 193 decompose the atomic displacements of the Ta atoms at any temperature in terms of the independent 194 amplitudes of the $A_{2u}$ and $E_u$ modes described by $q_1$ and $q_1 + q_2$. The displacements at 15 K are shown in the 195 insert to Figure 3f. In the NC phase at room temperature⁴⁴, $a_{NC} = 0.2448(2)$ and $b_{NC} = 0.0681(2)$, which is close to 196 the C values ($a_c = 3/13$ and $b_c = 1/13$), while in the IC phase, $a_{IC} = 0.283(2)$, and $b = 0$, where $a_{IC}$ is close to 197 2/7=0.2857. Note that our PDF analysis does not reveal any information about long range order, and $q_1$ and $q_2$ 198 should not be interpreted to imply its existence. Rather, they describe the symmetry of the lattice distortions 199 associated with the twists and displacements of the Ta atoms in the distorted SoD polaron structures. (See 200 **Supplementary Fig. 4**).

201 **DISCUSSION**

202 The local structure data in 1T-TaS₂ reveals symmetry-specific local structural fingerprints of polaron formation 203 whose ordering evolves with temperature as shown schematically in Fig. 4f. Surprisingly, the polaron fingerprints 204 are already identifiable in the polymorphic 6R-TaS₂ regime (where individual 1T-monolayers appear between 1H 205 layers), in the metallic phase (> 600 K), then successively through the IC, NC and C ordering transitions, and
eventually in the non-QSL-like regime below 50 K, where – remarkably – SoD symmetry is restored in a layer-dimerized undistorted in-plane SoD structure. An important fundamental question arises in the context of conventional CDW viewpoint regarding the lattice distortions in the high-temperature M phase: Should the deformations be discussed in terms of a uniform reduced modulation amplitude in the IC-like phase, or in terms of dynamical polaron fluctuations? Figure 1d indicates the range of the structural correlations. The black trace shows a difference across the IC-M transition, while the red trace is a difference corresponding to two PDFs in the same (M) phase separated by the same $\Delta T'$. This comparison shows that dominant changes across IC-M occur at $r > 5$ Å, with smaller changes for $r < 5$ Å (Supplementary Fig. 5). This implies that at the M to IC transition both the number of distorted sites and the extent of IC spatial correlations grow. The M phase thus looks like a sparse polaron gas, whereas the IC phase is a sparse polaron crystal (incommensurate with the lattice). We can now consider the three possible scenarios for the IC melting transition based on the PDF data: In the first one, in the M phase, the long-range IC-CDW becomes dynamic, with amplitude and phase fluctuations on a timescale faster than the lattice can respond. In the second scenario, all relevant charges become fully itinerant on all length scales in the M phase without lattice distortions. In the third, polaronic, scenario the IC-CDW breaks up in the M regime but the lattice follows the localized charge fluctuations, which is equivalent to thermally activated polaron hopping. In all three scenarios, we would see a disappearance of the reflections at the IC modulation wave vector on going from IC to M in Bragg reflections, with the third scenario leaving some diffuse scattering in the M phase. In the PDF data, the first and second scenario would result in no characteristic polaronic distortions in the M phase. The fact that we *do* see them implies that only the polaron fluctuation scenario is relevant in 1T-TaS$_2$. The persistence of polarons in the polymorphic regime (Fig. 2 b,c) and the change from 20 wt% to 10 wt% on the transition from a uniform 1T to the 6R multilayer structure confirms this notion and implies that polarons exist only in individual 1T monolayers. These may be expected to appear also in free-standing monolayers or monolayers on substrates, provided that strain and interaction with the substrate do not interfere.

**Magnetic ordering.** It is interesting to consider the symmetry of the magnetic polaron structures that are compatible with the two possible in-plane magnetic space groups $P\bar{3} \#147.13$ and $P\bar{3}' \#147.15$. Two (ferro)magnetic $P\bar{3}m$ subgroups can be derived from the analysis of the symmetric irreducible representation at the $\Gamma$ point of the trigonal Brillouin zone of the 1T stack. For $P\bar{3}m$, inversion does not flip the magnetic moment, allowing a magnetic moment component along *c* for Ta1 and no constraints for the magnetic moments of the Ta2 and Ta3 atoms. In contrast, for the $P\bar{3}'m$ magnetic space group, inversion reverses time and thus flips the magnetic moment, the magnetic moments of the Ta atoms of the same kind are not flipped by the symmetry operators. By choosing different components for Ta2 and Ta3 magnetic moments it is possible to create different motifs maintaining the same symmetry. Figs. 4 g,h show the motifs compatible with the symmetry operations of the two magnetic groups for in-plane chiral and non-chiral arrangements of magnetic moments respectively. An in-plane QSL-like state is consistent with fluctuating in-plane moments, and *c*-axis antiparallel moments in the $P\bar{3}'$ magnetic group, consistent with dimerized stacking of CDW orders in the C state.

A remarkable feature of the data is that local in-plane mirror symmetries appear to be restored below 50 K, with $\alpha \approx 120^{\circ}$ and $\beta \approx 150^{\circ}$. The data below 50 K are strongly suggestive of dimerization along the *c* axis, consistent with inter-plane spin singlet formation and spin gap formation, that could explain the apparent cross-over from a ~$T^2$ QSL-like spin relaxation above 50 K to a gapped, exponential T dependence below this temperature$^2$. A plausible origin of the apparent symmetry restoration at low temperature comes from the disappearance of domain walls at low temperature. In-plane domain walls have been observed to persist well below the NC-C transition temperature, which topologically prevent interlayer stacking and spin dimerization. As long as domain walls exist, un-dimerized layers support an intra-layer QSL composed of spin structures shown in Figs. 4g and h, while regions with dimerized lattice layers are expected to exhibit inter-layer spin pairing$^3$. The restoration of SoD symmetry at lowest temperatures may thus be considered to be associated with the disappearance of DWs. A slightly different situation occurs in the metastable ‘hidden’ domain state$^{65}$, where a QSL is allowed to form within domains whenever inter-layer stacking is indirect.
The polaron crystallization paradigm<sup>66</sup> emerging on the basis of these data has the benefit of providing a common framework for the understanding the polaronic physics in both the equilibrium and metastable phases of 1T-TaS<sub>2</sub> and numerous related materials that display similar charge and spin ordering<sup>23</sup>. The local structure symmetry analysis is particularly relevant for revealing the origin of the still poorly understood spin polaron structure and reconciles the observation of QSL-like behavior with theoretical band structure modelling, which commonly suggests a spin-paired dimerized ground state in a perfect crystal<sup>3</sup>. Finally, we note that it was shown theoretically that Wigner-crystal ordering is not limited to systems with a half-filled electronic band, but may be also present at other fillings<sup>21,67</sup>, which suggest that the crystallization phenomenon reported here may be more general than previously thought.

## 263 METHODS

**Crystal growth and analysis.** The 1T-TaS<sub>2</sub> samples were grown using iodine vapor transport reaction, grown at 850 °C, and quenched to room temperature from the growth temperature. Single crystal XRD shows a pure single 1T phase is retained after the quench, with Ta:S composition determined by EDS to be 33:66 ± 1 at %.

**Synchrotron X-ray atomic pair distribution function (PDF) measurements:** Temperature dependent X-ray total scattering data were collected at 28-ID-1 beamline of the National Synchrotron Light Source II (NSLS II) at Brookhaven National Laboratory. Finely ground powders of 1T-TaS<sub>2</sub> were sealed in 1mm (outer diameter) polyimide capillary (referred to as experiment 1 or exp 1) and 1.5 mm (outer diameter) quartz capillary (experiment 2 or exp 2) within a glovebox under light vacuum. Measurements were carried out in capillary transmission geometry using a 2D PerkinElmer amorphous silicon area detector (2048 × 2048 pixels with 200 µm<sup>2</sup> pixel size) placed ~204 mm downstream of the sample. The setup utilized a monochromatic X-ray beam with 74.5 keV energy (λ = 0.1665 Å). Sample temperature control was achieved using a Cryo Industries of America cryostat (exp 1, 15 K ≤ T ≤ 500 K) and a FMB Oxford Hot Air Blower model GSB1300 (exp 2, 300 K ≤ T ≤ 915 K). Data for each experiment were collected in 5 K steps using 5 K/minute temperature ramp, 2 minutes thermalization, and 2 minutes data collection at each temperature. In exp1 data in temperature range 15 K ≤ T ≤ 300 K were collected on cooling, while these in temperature range 300 K ≤ T ≤ 500 K were collected on warming and cooling. In exp 2 data were collected on warming.

**PDF data processing and analysis:** Calibrations of the experimental geometry, momentum transfer range, and detector orientation were carried out by utilizing nickel standard measurements performed under the same conditions. Appropriate masking of the beam-stop shadow, inactive and outlier pixels, and subsequent azimuthal integration of the 2D images to obtain 1D diffraction patterns of intensity versus Q data were done using pyFAI software package<sup>68,69</sup> Standardized corrections to the data for experimental effects to obtain the reduced total scattering structure function, F(Q), and the subsequent sine Fourier transforms to obtain experimental PDFs, G(r), with Q<sub>max</sub>=25 Å<sup>-1</sup> (exp 1) and Q<sub>max</sub>=20 Å<sup>-1</sup> (exp 2) were carried out using the PDFgetX3 program within the xPDFsuite software package<sup>70</sup> The PDF analysis was carried out using the PDFgui<sup>71</sup> modeling platform. The PDF, derived from powder diffraction-based Bragg and diffuse scattering data collected over a broad range of momentum transfer, describes direction-averaged distribution of atom pairs in a material as a function of interatomic distance, r, and provides structural information across different length scales. Uncertainties in reported model-derived parameters were estimated from their maximum variations observed by changing either the r-range of modelling (Fig. 1b, Fig. 3a-f) or integration (Fig. 2b, inset) by 10%. A double of the maximal values of such uncertainties was adopted as the error bars shown. Due to large X-ray scattering contrast (Z(Ta)=73, Z(S)=16) the dominant contribution to total PDF (Fig. 1e) originates from Ta-Ta pairs (red trace), followed by the Ta-S pairs' contribution (blue trace). The S-S (green trace) contribution is about 20 times weaker than the Ta-Ta contribution and 4 times weaker than that of Ta-S pairs.

**Estimate of distorted fraction and structure modeling in polymorphic regime:** By using a calculated PDF based on undistorted P<span style="text-decoration: overline">3</span>m model as a reference, we form difference PDFs for all T>600 K data. Such differences are
then integrated over a 2.6-4 Å range (highlighted region on abscissa in Fig. 2b) and normalized to the 600 K value to quantify the relative change with temperature of the fraction of distorted sites across the polymorphic transformations, as shown in inset to Fig. 2b (for illustration see Supplementary Fig. 1). Above 730 K the short-range PDF data can be explained by a 6R-TaS₂-like model (R3m symmetry) featuring alternating 1T and 1H layers in equal abundance, Fig. 2c (bottom). Although the exact stacking and long-range character are likely more complex, as the model does not fully explain the data over longer length-scales, local analysis is rather robust. Importantly, the 6R-TaS₂ model features undistorted 1T & 1H layers, hence not accounting for local distortions responsible for nontrivial nearest neighbor Ta-Ta distance distribution.

**Local structure model with P<span style="text-decoration: overline">3</span> symmetry:** To map out temperature evolution of local interatomic distances, first data fits were carried out over 1 nm range and utilized a simple single S-Ta-S layer 1T-TaS₂ model conforming to P<span style="text-decoration: overline">3</span> symmetry. The model includes 19 independent structure parameters, constrained by symmetry to describe fractional coordinates of Ta (6 parameters) and S (13 parameters) atoms in the unit cell, two lattice parameters, two thermal parameters, and an overall scale factor. This explains the data features in the C state very well. The refinement range was established as a compromise between the model over-parametrization (data window too narrow) and the model deficiency (data window too wide) regimes. The 1 nm fitting range encompasses intra-bilayer correlations but excludes inter-bilayer correlations, which were subsequently explored in a model-independent manner as described below. (See Supplementary Table 1, Supplementary Figures 6 and 7, and Section 2 for additional model details.)

**Multiphase treatment:** The fits using only the P<span style="text-decoration: overline">3</span> phase in the NC regime were noticeably worse compared to these in the C state, with observable increase of fit residual, as shown in inset to Fig. 3c. This is consistent with the presence of discommensuration domain walls and lower SoD density in the NC state whose atomic structure is less distorted and not accounted for in our simple distorted P<span style="text-decoration: overline">3</span> model. We approximated domain wall contribution by adding an undistorted P<span style="text-decoration: overline">3</span>m structure as a secondary (minority) phase to account for local environment of the discommensuration regions. We utilized this approach to fit the PDF data in the IC and M regimes as well, where the distorted and undistorted phases swap roles and the distorted phase becomes a minority phase.

**Assignment and consideration of inter-layer PDF peaks:** The inter-layer PDF peaks considered in this study were assigned based on undistorted P<span style="text-decoration: overline">3</span>m model as reported for 300 K and its supercell variants. In this, interlayer spacings along the stacking direction are defined by multiples of c-axis lattice parameter, including ~5.9 Å which corresponds to intra-bilayer separation. The ~13.1 Å peak in undistorted P<span style="text-decoration: overline">3</span>m structure corresponds to a well-defined off-axis Ta-Ta contribution connecting two Ta atoms that are 2 layers apart. These Ta atoms belong to next nearest neighbor Ta-layers and are offset within planes away from the stacking axis. In the high temperature regime, if structural distortions are incoherent, this peak in the data would be smeared, as is observed. As order establishes, re-appearance (with respect to putative P<span style="text-decoration: overline">3</span>m structure that never materializes locally) of a sharp signal at ~13.1 Å implies re-emergence of structural coherence along this interatomic vector. The signal strength likely implies alignment of SoDs belonging to neighboring bilayers along this direction.

## DATA AVAILABILITY

All of the data supporting the conclusions are available within the article and the Supplementary Information.

Additional data are available from the corresponding authors upon request.

## REFERENCES
340 1. Law, K. T. & Lee, P. A. 1T-TaS2 as a quantum spin liquid. *Proc National Acad Sci* **114**, 6996–7000
341 (2017).

342 2. Klanjsek, M. et al. A high-temperature quantum spin liquid with polaron spins. *Nat Phys* **13**, 1130–
343 1134 (2017).

344 3. Darancet, P., Millis, A. J. & Marianetti, C. A. Three-dimensional metallic and two-dimensional insulating behavior in octahedral tantalum dichalcogenides. *Phys Rev B* **90**, 045134 (2014).

346 4. Wilson, J. A., Salvo, F. J. D. & Mahajan, S. Charge-density waves and superlattices in the metallic layered transition metal dichalcogenides. *Adv Phys* **24**, 117–201 (1975).

348 5. Kogar, A. et al. Signatures of exciton condensation in a transition metal dichalcogenide. *Science*
349 **358**, 1314–1317 (2017).

350 6. Persky, E. et al. Magnetic memory and spontaneous vortices in a van der Waals superconductor. *Nature* **607**, 692–696 (2022).
351

352 7. Hamill, A. et al. Two-fold symmetric superconductivity in few-layer NbSe2. *Nat Phys* **17**, 949–954
353 (2021).

354 8. Bednorz, J. G. & Müller, K. A. Possible high Tc superconductivity in the Ba–La–Cu–O system. *Zeitschrift Für Physik B Condens Matter* **64**, 189–193 (1986).

356 9. Sipos, B., Berger, H., Forro, L., Tutis, E. & Kusmartseva, A. F. From Mott state to superconductivity in 1T-TaS2. *Nat Mater* **7**, 960–965 (2008).

358 10. Dong, Q. et al. Structural phase transition and superconductivity hierarchy in 1T-TaS2 under pressure up to 100 GPa. *Npj Quantum Mater* **6**, 20 (2021).

360 11. Zhang, P. et al. All-optical switching of magnetization in atomically thin CrI3. *Nat Mater* **1–6**
361 (2022) doi:10.1038/s41563-022-01354-7.

362 12. Stojchevska, L. et al. Ultrafast Switching to a Stable Hidden Quantum State in an Electronic
363 Crystal. *Science* **344**, 177–180 (2014).

364 13. Ravnik, J. et al. A time-domain phase diagram of metastable states in a charge ordered quantum material. *Nat Commun* **12**, 2323 (2021).

366 14. Dąbrowski, M. et al. All-optical control of spin in a 2D van der Waals magnet. *Nat Commun* **13**,
367 5976 (2022).

368 15. Kogar, A. et al. Light-induced charge density wave in LaTe3. *Nat Phys* **52**, 1–6 (2019).

369 16. Tedstone, A. A., Lewis, D. J. & O’Brien, P. Synthesis, Properties, and Applications of Transition Metal-Doped Layered Transition Metal Dichalcogenides. *Chem Mater* **28**, 1965–1974 (2016).

371 17. Qiao, S. et al. Mottness Collapse in 1T-TaS2-xSex Transition-Metal Dichalcogenide: An Interplay
372 between Localized and Itinerant Orbitals. *Physical Review X* **7**, 041054 (2017).
373 18. Li, L. J. et al. Fe-doping-induced superconductivity in the charge-density-wave system 1T-TaS2.

374 *Epl-europhys Lett* **97**, 67005 (2012).

375 19. Wang, Y. D. et al. Band insulator to Mott insulator transition in 1T-TaS2. *Nat Commun* **11**, 1–7

376 (2020).

377 20. Landau, L. Electron motion in crystal lattices. *Phys. Z. Sowjetunion* **3**, 644–645 (1933).

378 21. Camjayi, A., Haule, K., Dobrosavljević, V. & Kotliar, G. Coulomb correlations and the Wigner–Mott

379 transition. *Nat Phys* **4**, 932–935 (2008).

380 22. Mertelj, T., Kabanov, V. V., Mena, J. M. & Mihailovic, D. Self-organization of charged particles on

381 a two-dimensional lattice subject to anisotropic Jahn-Teller-type interaction and three-dimensional

382 Coulomb repulsion. *Phys Rev B* **76**, 9 (2007).

383 23. Vodeb, J. et al. Configurational electronic states in layered transition metal dichalcogenides. *New

384 *J Phys* **21**, 083001 (2019).

385 24. Mihailovic, D. Why Does Maximum T c Occur at the Crossover From Weak to Strong Electron–

386 phonon Coupling in High-temperature Superconductors? doi:10.1007/s10948-022-06152-5.

387 25. Rossnagel, K. On the origin of charge-density waves in select layered transition-metal

388 dichalcogenides. *J Phys Condens Matter* **23**, 213001 (2011).

389 26. Franchini, C., Reticcioli, M., Setvin, M. & Diebold, U. Polarons in materials. *Nat Rev Mater* **6**, 560–

390 586 (2021).

391 27. Mihailovic, D. et al. Application of the polaron-transport theory to sigma (omega) in Ti2Ba2Ca1-

392 xGdxCu2O8, YBa2Cu3O7-d, and La2-xSrxCuO4. *Phys Rev B* **42**, 7989–7993 (1990).

393 28. Fazekas, P. & Tosatti, E. Charge Carrier Localization in Pure and Doped 1t-TaS2. *Physica B & C* **99**,

394 183–187 (1980).

395 29. Gerasimenko, Y. A. et al. Quantum jamming transition to a correlated electron glass in 1T-TaS 2.

396 *Nat Mater* **317**, 505–1083 (2019).

397 30. Mraz, A. et al. Charge Configuration Memory Devices: Energy Efficiency and Switching Speed.

398 *Nano Lett* **22**, 4814–4821 (2022).

399 31. Mihailovic, D. et al. Ultrafast non-thermal and thermal switching in charge configuration memory

400 devices based on 1T-TaS. *APL* **119**, 013106 (2021).

401 32. Liu, G. et al. A charge-density-wave oscillator based on an integrated tantalum disulfide–boron

402 nitride–graphene device operating at room temperature. *Nat Nanotechnol* **11**, 845–850 (2016).

403 33. Landuyt, J. V., Tendeloo, G. V. & Amelinckx, S. Electron diffraction study of inter- and

404 intrapolytypic phase transitions in transition metal dichalcogenides. *IV Cinematogeraphic*

405 Observations of the 1Tβ → 1Tγ-Phase Transition in TaS2. *Phys. Status Solidi A* **42**, 565–570 (1977).

406 34. Revelli, J. F. & Phillips, W. A. Studies of the system TaS2–xSex. *J Solid State Chem* **9**, 176–186

407 (1974).
408 35. Givens, F. L. & Fredericks, G. E. Thermal expansion of NbSe2 and TaS2. *J. Phys. Chem. Solids* **38**, 1363-1365 (1977).

410 36. Ziebeck, K. R. A., Dorner, B., Stirling, W. G. & Schollhorn, R. Kohn anomaly in the 1T 2phase of TaS
411 2. *J Phys F Metal Phys* **7**, 1139-1143 (1977).

412 37. Machida, Y. et al. Observation of Soft Phonon Modes in 1T-TaS2 by means of X-ray Thermal
413 Diffuse Scattering. *J Phys Soc Jpn* **73**, 3064-3069 (2004).

414 38. Nakanishi, K. & Shiba, H. Domain-Like Incommensurate Charge-Density-Wave States and
415 Collective Modes. *J Phys Soc Jpn* **45**, 1147-1156 (1978).

416 39. Nakanishi, K., Shiba, H., Takatera, H. & Yamada, Y. The nearly commensurate phase and effect of
417 harmonics on the successive phase transition in 1 T-TaS 2. *Journal of the Physical Society of Japan* **43**,
418 1509-1517 (1977).

419 40. McMillan, W. L. Theory of discommensurations and the commensurate-incommensurate charge-
420 density-wave phase transition. *Phys Rev B* **14**, 1496-1502 (1976).

421 41. Butler, C. J., Yoshida, M., Hanaguri, T. & Iwasa, Y. Mottness versus unit-cell doubling as the driver
422 of the insulating state in 1T-TaS2. *Nat Commun* **11**, 2477 (2020).

423 42. Hellmann, S. et al. Time-domain classification of charge-density-wave insulators. *Nat Commun* **3**,
424 1069-1069 (2012).

425 43. Ritschel, T. et al. Orbital textures and charge density waves in transition metal dichalcogenides.
426 *Nat Phys* **11**, 328-331 (2015).

427 44. Spijkerman, A., Boer, J. de, Meetsma, A. & Wiegers, G. X-ray crystal-structure refinement of the
428 nearly commensurate phase of 1T-TaS2 in (3+ 2)-dimensional superspace. *Physical Review B* **56**,
429 13757 (1997).

430 45. Jellinek, F. The system tantalum-sulfur. *J Less Common Metals* **4**, 9-15 (1962).

431 46. Svetin, D., Vaskivskyi, I., Brazovskii, S. & Mihailović, D. Three-dimensional resistivity and
432 switching between correlated electronic states in 1T-TaS2. *Sci Rep-uk* **7**, 46048 (2017).

433 47. Brouwer, R. & Jellinek, F. Superstructure of 1s-TaS2. *Mater Res Bull* **9**, 827-829 (1974).

434 48. Drummond, N. D. & Needs, R. J. Phase Diagram of the Low-Density Two-Dimensional
435 Homogeneous Electron Gas. *Phys Rev Lett* **102**, 126402 (2008).

436 49. Pankov, S. & Dobrosavljević, V. Self-doping instability of the Wigner-Mott insulator. *Phys Rev B*
437 **77**, 085104 (2008).

438 50. Radonjić, M. M., Tanasković, D., Dobrosavljević, V., Haule, K. & Kotliar, G. Wigner-Mott scaling of
439 transport near the two-dimensional metal-insulator transition. *Phys Rev B* **85**, 085133-7 (2012).

440 51. Kochubey, D. I., Kim, T. K., Babenko, V. P. & Shaikhutdinov, S. K. Charge density waves in 1T-TaS2:
441 an EXAFS study. *Phys B Condens Matter* **252**, 15-20 (1998).
442 52. Ishiguro, T. & Sato, H. Electron microscopy of phase transformations in 1T-TaS2. *Phys Rev B* **44**,
443 2046–2060 (1991).

444 53. Ishiguro, T. & Sato, H. High-resolution electron microscopy of discommensuration in the nearly
445 commensurate phase on warming of 1 T- TaS2. *Phys Rev B* **52**, 759–765 (1995).

446 54. Ishiguro, Y. et al. Layer Number Dependence of Charge Density Wave Phase Transition Between
447 Nearly-Commensurate and Incommensurate Phases in 1T-TaS2. *J Phys Chem C* **124**, 27176–27184
448 (2020).

449 55. Billinge, S. J. L. & Levin, I. The Problem with Determining Atomic Structure at the Nanoscale.
450 *Science* **316**, 561–565 (2007).

451 56. Billinge, S. & Egami, T. *Underneath the Bragg Peaks*. (2003). doi:10.1016/s1470-1804(03)80011-1.

453 57. Brouwer, R., C, F. J. P. B. & 1980. The low-temperature superstructures of 1T-TaSe2 and 2H-TaSe2. *Phys B C* **99**, 51–55 (1980).

455 58. Tidman, J. P. & Frindt, R. F. Resistivity of thin TaS2 crystals. *Can J Phys* **54**, 2306–2309 (1976).

456 59. Achari, A. et al. Alternating Superconducting and Charge Density Wave Monolayers within Bulk
457 6R-TaS2. *Nano Lett* **22**, 6268–6275 (2022).

458 60. Koch, R. J. et al. Dual Orbital Degeneracy Lifting in a Strongly Correlated Electron System. *Phys
459 Rev Lett* **126**, 186402 (2021).

460 61. Bozin, E. S. et al. Local orbital degeneracy lifting as a precursor to an orbital-selective Peierls
461 transition. *Nat Commun* **10**, 3638 (2019).

462 62. Sezerman, O., Simpson, A. M. & Jericho, M. H. Thermal expansion of 1T-TaS2 and 2H-NbSe2. *Solid
463 State Commun* **36**, 737–740 (1980).

464 63. Ma, L. et al. A metallic mosaic phase and the origin of Mott-insulating state in 1T-TaS2. *Nat
465 Commun* **7**, 1–8 (2016).

466 64. Cho, D. et al. Nanoscale manipulation of the Mott insulating state coupled to charge order in 1T-TaS2. *Nat Commun* **7**, 10453 (2016).

468 65. Gerasimenko, Y. A., Karpov, P., Vaskivskyi, I., Brazovskii, S. & Mihailovic, D. Intertwined chiral
469 charge orders and topological stabilization of the light-induced state of a prototypical transition
470 metal dichalcogenide. *Npj Quantum Mater* **4**, 1–9 (2019).

471 66. Fratini, S. & Quémerais, P. Polaron Crystallization and Melting: Effects of the Long-Range
472 Coulomb Forces. *Mod Phys Lett B* **12**, 1003–1012 (1998).

473 67. Kapcia, K. J., Robaszkiewicz, S., Capone, M. & Amaricci, A. Doping-driven metal-insulator
474 transitions and charge orderings in the extended Hubbard model. *Phys. Rev. B* **95**, 125112 (2017).

475 68. Ashiotis, G. et al. The fast azimuthal integration Python library: pyFAI. *J Appl Crystallogr* **48**, 510–519 (2015).
477 69. Kieffer, J., Valls, V., Blanc, N. & Hennig, C. New tools for calibrating diffraction setups. *J Synchrotron Radiat* **27**, 558–566 (2020).

478

479 70. Juhás, P., Davis, T., Farrow, C. L. & Billinge, S. J. L. PDFgetX3: a rapid and highly automatable

480 program for processing powder diffraction data into total scattering pair distribution functions. *J Appl Crystallogr* **46**, 560–566 (2013).

481

482 71. Farrow, C. L. et al. PDFfit2 and PDFgui: computer programs for studying nanostructure in crystals.

483 *J Phys Condens Matter* **19**, 335219 (2007).

484

485
487 ACKNOWLEDGMENTS

488 We wish to acknowledge Igor Vaskivskyi, Jaka Vodeb and Viktor Kabanov for valuable discussions. Work at
489 Brookhaven National Laboratory is supported by the Office of Basic Energy Sciences, Materials Sciences, and
490 Engineering Division, U.S. Department of Energy (DOE) under Contract No. DE-SC0012704. This research used
491 beamline 28-ID-1 of the National Synchrotron Light Source II, a U.S. Department of Energy (DOE) Office of
492 Science User Facility operated for the DOE Office of Science by Brookhaven National Laboratory under Contract
493 No. DE-SC0012704. DM and PS wish to acknowledge funding from ARRS. The work at the Jozef Stefan Institute
494 was supported by the Slovenian Research Agency (P1-0040).

495 AUTHOR CONTRIBUTIONS

496 EB and MA designed and conducted the X-ray experiments, EB performed PDF analysis, GB and SC performed
497 group theoretical analysis. DM and EB conceptually proposed the experiments. EB and DM wrote the paper. PS
498 synthesized and characterized the crystals.

499 COMPETING INTERESTS

500 The authors declare no competing interests in this work.
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Electronic transitions & lattice symmetry breaking in 1T-TaS₂. (a) Temperature dependent resistivity, adopted from Sipos et al.⁹, highlighting electronic phases in 1T-TaS₂. The insert shows the undistorted P₃m lattice structure. (b) Temperature evolution of in-plane ADP of Ta from P₃m model fit to PDF data over 20-60 Å range. Red color indicates heating, blue color indicates cooling. Sloping black dotted line is a reference. Inset: Fit residual, R<sub>w</sub>, implicates IC-M and IC-NC transitions. (c) P₃m model fit to 600 K data. Green difference trace (offset for clarity) reveals short range deviations. Stack of difference curves for the same model against data at different temperatures (as indicated) is shown underneath. Red arrow indicates a length-scale coarsely corresponding to the SoD polaron and the in-plane nearest neighbor inter-star center-to-center distance (inset). (d) Comparison of difference PDF between experimental data that are 10 K apart within metallic phase (T > 550 K) and across the IC-M transition, implicating existence of local distortions in the M phase over at least a radius r. (e) Calculated total PDF and partial pair-contributions in the P₃m structure. (f) False color plot of PDF intensity versus temperature over a short r-range.</figcaption></figure>
516

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Temperature progression of atomic structure in 1T-TaS₂. (a) Stack of diffraction patterns revealing a polypomorphic transformation on warming at ~630 K. Inset: intensity collapse of the (012) P<span style="text-decoration: overline">3</span>m Bragg reflection during the restructuring from pure 1T to a hybrid 1T & 1H coordination. Further restructuring above 880 K results in additional Bragg reflections, marked by asterisk. (b) Corresponding PDFs over a narrow *r*-range. Data feature isosbestic points (circled) - nodes of temperature-invariant PDF intensity. Label 1T* indicates presence of local distortions in the 1T-type layers. Vertical blue arrow (1T*) marks polaronic distortion intensity. Inset: Change of fraction of distorted 1T* environments, relative to 600 K, estimated over the shaded region on the abscissa. (c) Polaron signature, seen in the M phase (top), persists in the restacked polymorph structure (bottom), and is associated with the 1T-type layers. Adding distortions, described in (g), as nanoscale inclusions results in improved fit (cyan trace). In (d) and (e) a stack of PDF data on different length scales is shown over the 15 K – 500 K range, color-coded according to the electronic phase they belong to. At high temperature longer length-scale peaks sharpen (indicated by vertical arrows) due to apparent average symmetry increase, despite expected broadening due to ADP increase at elevated temperature and in contrast to the shorter length-scale behavior. Features exhibiting negative thermal expansion (NTE) are indicated by block arrows. (f) At 15 K clear distortions associated with the star of David (SoD) are seen in the local structure, incompatible with the P<span style="text-decoration: overline">3</span>m symmetry (top) but explained well (bottom) by a simple P<span style="text-decoration: overline">3</span> broken symmetry model, depicted in (g) with focus on Ta atoms. Colored arrows show the Ta displacements. Out-of-plane degrees of freedom (exaggerated in the illustration for clarity) provide for puckering distortions.</figcaption></figure>
535

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Local structure quantification from distorted local SoD model. (a) Local Ta-Ta distances from a P$\bar{3}$ model fits to PDF data, obtained on cooling, over 1 nm length-scale in 15-500 K range. For reference, dashed horizontal gray line marks the undistorted P$\bar{3}$m value at 300 K. Intra- and inter-star distances are color coded as indicated in the legend. (b) Closeup of C phase behavior (left) and hysteretic response across the IC-NC transition (right). (c) In NC and IC phases secondary undistorted P$\bar{3}$m phase was necessary to reduce fit residual, R$_w$ (inset). Distorted phase fraction, displaying thermal hysteresis consistent with one seen in resistivity, is calculated from atomic weight-based phase content obtained from the fits. (d)-(f) show selected intra-SoD and inter-SoD bond angles, as indicated in the sketch. In the IC phase the SoD motif is heavily distorted, with local distortions resembling discs. Upon entering the NC phase, the SoD distortion becomes better differentiated, with stars exhibit twists, and with interatomic angles beginning to evolve towards ideal values. In the C phase twists start to diminish. Below ~50 K a sharp regularization of stars is observed, with intra-star Ta-Ta distances becoming equal, intra-star angles approaching ideal values, and no twisting. Colored arrows in (b)-(f) indicate data collection on warming (red) and on cooling (blue) cycles. Inset to (f): displacements of the Ta atoms for the $q_1$ (red) and $q_1+q_2$ (blue) modes from symmetry analysis at 15 K, compared to the parent P$\bar{3}$m phase (gray).</figcaption></figure>
550

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Model independent considerations of stacking correlations. (a) Temperature stack of experimental PDFs around ~5.9 Å peak containing intra-bilayer correlations of SoDs, indicated in (e) by blue and red double arrows. (b) Temperature stack of PDFs around ~13.1 Å peak depicting inter-bilayer correlations of SoDs, indicated in (e) by purple double arrows. The contribution assignment is made based on features in P$\bar{3}$m model, Fig. 1e. There, the 5.9 Å feature is the c-axis lattice repeat distance involving strong signal from the nearest neighbor interlayer Ta-Ta correlations, whereas the 13.1 Å feature is well isolated and dominated by the next nearest neighbor off-axis interlayer Ta-Ta correlations. In (c) and (d) temperature stacks of differential PDFs are shown for data in (a) and (b), respectively, with differentials calculated using 500 K reference. Vertical dashed lines mark positions of apparent maxima in the C phase. (e) Sketch of two bilayers highlighting correlations discussed in text. Inter-bilayer SoDs are separated by ~13.1 Å at low temperature, a peak forming on cooling from a broad distribution at high temperature. This separation results from a SoD-shift by approximately one P$\bar{3}$m lattice spacing in each P$\bar{3}$m lattice direction. Note 6 possible choices of relative positioning of nearest bilayers, indicated by the purple shaded hexagon. (f) Schematics of distortions across different electronic phases in 1T-TaS$_2$. (g, h) magnetic moments consistent with the distortions of the SoDs as determined by the PDF analysis for two different settings. The arrows at the center Ta1 atom pointing along the z axis are shown in green (only in P$\bar{3}$m).</figcaption></figure>

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567
Supplementary Information

Crystallization of polarons through charge and spin ordering transitions
in 1T-TaS<sub>2</sub>

E.S. Bozin<fnref n="1" />, M. Abeykoon<fnref n="2" />, S. Conradson<fnref n="3" />, G. Baldinozzi<fnref n="4" />, P. Sutar<fnref n="3" /> and D. Mihailovic<fnref n="3" />

1. SUPPLEMENTARY FIGURES AND TABLES 2

Supplementary Figure 1. Polymorphic transformations in 1T-TaS<sub>2</sub> on warming in the M phase. 6

Supplementary Figure 2. Evolution of stacking correlations in 1T-TaS<sub>2</sub>. 3

Origin of c-axis negative thermal expansion (NTE) jumps 4

Supplementary Figure 3. Possible origin of the c-axis NTE jumps. 4

Supplementary Figure 4. Comparison of PDF data of 1T-TaS<sub>2</sub> in the low temperature regime. 5

Supplementary Figure 5. Local structure changes across M-IC transition in 1T-TaS<sub>2</sub>. 6

Supplementary Table 1. Details of the single layer 13 × 13 supercell P3 model. 7

Supplementary Figure 6. Features of the P3 structure fit to the 15 K PDF data of 1T-TaS<sub>2</sub> in the C state. 8

Supplementary Figure 7. Considerations of the fitting range and the applicability of the P3 approximant model. 9

2. ORDER PARAMETERS DESCRIBING THE SYMMETRY CHANGE LEADING TO THE COMMENSURATE CDW LOW TEMPERATURE PHASE OF 1T-TAS<sub>2</sub>. 11

Revisiting the 1T (single layer) C-CDW: atomic vibrations 11

The low T variants 12

Symmetry adapted coordinates of Ta displacements at q<sub>1</sub> and at q<sub>1</sub>+q<sub>2</sub>. 15

Cartesian components of the displacements of each independent Ta atom in the daughter phase 16

Supplementary table 2. Fractional displacement contributions originating from the amplitudes transported at q<sub>1</sub> and q<sub>1</sub>+q<sub>2</sub>. 17

<footnote n="1">¹Condensed Matter Physics and Materials Science Division, Brookhaven National Laboratory, Upton, NY 11973, USA</footnote>

<footnote n="2">²Photon Sciences Division, Brookhaven National Laboratory, Upton, NY 11973, USA</footnote>

<footnote n="3">³Dept. of Complex Matter, Jozef Stefan Institute, Jamova 39, SI-1000 Ljubljana, Slovenia</footnote>

<footnote n="4">⁴Centralesupélec, CNRS, SPMS, Université Paris-Saclay, bât Eiffel, Gif-sur-Yvette, Île-de-France, 91190, FRANCE</footnote>
# 1. SUPPLEMENTARY FIGURES AND TABLES

<figure><img src="image_6.png" /><figcaption>**Supplementary Figure 1. Polymorphic transformations in 1T-TaS₂ on warming in the M phase.**</figcaption></figure>

a Stack of diffraction patterns revealing two polymorphic transformations on warming, T > 600 K. Red stack data, T > 870 K, reveal additional transformation to that focused on in the main text. This additional restructuring likely involves a distorted mixture of 4Hb (P63/mmc) and 6R (R3m) polymorphs with significant stacking faults. Low resolution of our diffraction data did not allow unambiguous characterization. Tentative PDF assessment suggested a complex mixture of alternating 1T and 1H layers (inset). b Transformation temperatures charted by fit residuals of P3̄m and R3m models, as indicated in the panel. Inset: sketch of the free energy of TaS2, depicting metastable state (1T) and more stable configurations (1T+1H) achieved as temperature is elevated. c Stack of narrow *r*-range difference PDFs, ΔG, obtained by subtracting simulated PDF from a literature reference based on undistorted P3̄m at 300 K, from experimental PDFs (lower *r*-resolution, Qmax=20 Å⁻¹) collected at T > 600 K. The difference PDF is acquired for all M-phase data. To estimate the change in fraction across the polymorphic transformations, the difference PDF, ΔG, was integrated between 2.6 – 4 Å (isobestic points encompassing PDF feature of interest, between vertical arrows) for each temperature, and the final integral was normalized to the value obtained at 600 K. The uncertainties in the distorted fraction were estimated by extending the integration range over the grey window (nearest low-*r* zero crossings) and repeating the process. A double of the maximal value of such uncertainty was adopted as the error bars shown in the inset of Fig. 2b of the main text.
<figure><img src="image_7.png" /><figcaption>Supplementary Figure 2. Evolution of stacking correlations in 1T-TaS₂.</figcaption></figure>

a Idealized stacking of Ta-layers, with interlayer spacings as indicated, assigned based on $P\bar{3}m$ structure, and defined with respect to a reference layer. b Legend of color-coded temperatures of selected r-range PDF data shown in panels c-j. In these panels the r-ranges were selected to encompass features associated with the multiples of the $P\bar{3}m$ c-axis lattice parameter reported in literature for 300 K. While the PDF traces corresponding to the transitions to equilibrium electronic phases are colour coded, the other data in the stacks are shown in grey. The correlations up to 8c are shown. Profiles at 500 K are highlighted in red up to 4c range to emphasize the broadness of the high temperature signal which in part originates from interlayer disorder. Features with apparent negative thermal expansion characteristics are labelled as NTE.
<figure><img src="image_8.png" /><figcaption>Supplementary Figure 3. Possible origin of the c-axis NTE jumps.</figcaption></figure>

a Stack of PDF data around the nearest neighbour layer stacking distance at ~5.9 Å, as presented in Fig. 4a of the Main Text, and in Suppl. Fig. 3c. Dashed vertical line marks c-axis lattice parameter in undistorted P3̄m structure. Coloured traces mark 500 K (solid red line), 350 K (dashed red line), 200 K (dashed dark red line), 50 K (solid cyan line) and 15 K (solid blue line). b. In P3̄m structure layers are equidistant, as sketched. c. In distorted heavily disordered structure, such as that observed at higher temperatures, puckering distortions result in heavily distorted layers, resulting in interlayer stacking disorder and observed broad distribution of nearest neighbour interlayer distances, such as seen in a. In proposed NTE model this results in layer-shears and monoclinic compressions yielding shorter-than-expected average c-lattice parameter. d. As interlayer orders set in on cooling, the puff-pastry-like stacking disorder gets lifted, and the c-axis expands resulting in observed NTE jumps, as described in text.

Origin of c-axis negative thermal expansion (NTE) jumps

Dilatometry measurements in 1T-TaS2 [1] revealed anomalous NTE jumps in c-axis parameter, most significantly at the NC-C transition, but also at the IC-NC transition, origin of which has not been understood to date. Our data also reveal NTE behavior of certain PDF features, as marked in Figs. 2b,d,e and 4a of the Main Text, as well as in Supplementary Figure 3, and offer a simple explanation for this phenomenon. Based on these observations we propose a disorder-based rationale for NTE jumps. In putative undistorted P3̄m crystal structure 1T layers stacked along c-axis are equidistant, Suppl. Fig. 4b. However, as our PDF analysis demonstrates, 1T-TaS2 is locally distorted at all temperatures at which this metastable polymorph exists, including the highest temperature in the M phase, resulting in 1T*-TaS2. There, the distortions involve local puckering of Ta layers, such as schematically shown the top panel of Fig. 4f in the Main Text as well as in Supp. Fig. 4c. As a consequence, disordered layers exhibit a range of interlayer distances and concomitant sheer displacements due to Coulomb repulsion of polarons in neighbouring layers, evident from broad interlayer distance distributions that are, importantly, centered at nominal distance shorter than simple integer multiples of undistorted P3̄m c-axis lattice parameter (Main Text Fig. 4a,c at high temperature, and also Supplementary Fig. 3c and 4a) resulting in, on average, closer packing of the layers. Due to this, the experimentally measured reference c-axis lattice parameter at high temperature is anomalously short, resembling puff-pastry configuration as depicted in Fig. 4 f (top) and Supp. Fig. 4c. As polaron gas condenses into a
Wigner crystal, and the distortions become long range ordered and stacked, aided by Coulomb repulsion effects, the ordered layers space away, effectively removing monoclinic shears, resulting in increase of the average interlayer spacing (Main Text Fig. 4 a,c at low temperature), as sketched in Supp. Fig. 4d, leading to apparent NTE jumps observed in dilatometry.

<figure><img src="image_9.png" /><figcaption>Supplementary Figure 4. Comparison of PDF data of 1T-TaS₂ in the low temperature regime.</figcaption></figure>

a An overplot of experimental PDFs of 15 K (cyan trace) and 50 K (red trace), with corresponding differential (green trace) offset for clarity. The Ta-Ta contributions per P3m model used in this work are shown as blue ticks. b The r-region sensitive to intrastar correlations. Apart from sharpening, expected upon temperature reduction, centroid of the relevant peak also shifts to the left, consistent with SoD regularization revealed by modelling show in the Main Text. c The r-region sensitive to the nearest neighbour intra-bilayer SoD correlations. The peak also exhibits observable shifting, consistent with enhanced inter-bilayer binding. Vertical black line in b and c marks the centroid of feature of interest at 50 K. One possible origin of this behaviour, as discussed in the Main Text, is a formation of spin-singlet dimers associated with the stray spins, presumably from the $dz^2$ manifolds, in two adjacent SoDs within the bilayers. Notably, the amount of change in PDF signal seen in c is quite large. It is unlikely than only one Ta-Ta contact per SoD (1:13) contracts, implying that the mechanism is more complex than simple pairing of spins in $dz^2$ orbitals portrayed in d. Importantly, the features in the differential in a are well structured, stronger under some PDF peaks and weaker under others. Given that the differential is of two data sets, it is unlikely that they originate from artifacts such as termination ripples or noise in the data. Even if present, such effects would effectively cancel out in the subtraction, leaving only signature of subtle structural differences. We further note that as part of the experiment we carried out a counting statistics test to assess the signal to noise ratio. In this, data statistics was found to be satisfactory after counting for 10 seconds. The actual counting time was 120 seconds, an order of magnitude longer.
<figure><img src="image_10.png" /><figcaption>Supplementary Figure 5. Local structure changes across M-IC transition in 1T-TaS₂.</figcaption></figure>

a Comparison of difference PDF between data that are 10 K apart within metallic phase (T > 550 K) and across the IC-M transition, implicating existence of local distortions in the M phase over at least a radius r, as sketched in b. Note that these difficult-to-model distortions do not reflect star or David in the M phase, but heavily puckered hexagonal discs as indicated by modelling presented in the Main Text. Enhancement of difference signal on the length-scale r > 5 Å implies growth of structural correlations in the IC regime.
<figure><img src="image_11.png" /></figure>

<table><thead><tr><th>Tantalum</th><th>x</th><th>y</th><th>z</th></tr></thead><tbody><tr><td>1</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2</td><td>p1</td><td>p2</td><td>p3</td></tr><tr><td>3</td><td>1-p2</td><td>p1-p2</td><td>p3</td></tr><tr><td>4</td><td>1-p1+p2</td><td>1-p1</td><td>p3</td></tr><tr><td>5</td><td>1-p1</td><td>1-p2</td><td>-p3</td></tr><tr><td>6</td><td>p2</td><td>1-p1+p2</td><td>-p3</td></tr><tr><td>7</td><td>p1-p2</td><td>p1</td><td>-p3</td></tr><tr><td>8</td><td>p4</td><td>p5</td><td>p6</td></tr><tr><td>9</td><td>1-p5</td><td>p4-p5</td><td>p6</td></tr><tr><td>10</td><td>1-p4+p5</td><td>1-p4</td><td>p6</td></tr><tr><td>11</td><td>1-p4</td><td>1-p5</td><td>-p6</td></tr><tr><td>12</td><td>p5</td><td>1-p4+p5</td><td>-p6</td></tr><tr><td>13</td><td>p4-p5</td><td>p4</td><td>-p6</td></tr></tbody></table>

<table><caption>Supplementary Table 1. Details of the single layer $\sqrt{13} \times \sqrt{13}$ supercell $P\overline{3}$ model.</caption><thead><tr><th>Sulfur</th><th>x</th><th>y</th><th>z</th><th>Sulfur</th><th>x</th><th>y</th><th>z</th></tr></thead><tbody><tr><td>14</td><td>p7</td><td>p8</td><td>p9</td><td>27</td><td>1-p14</td><td>p13-p14</td><td>p15</td></tr><tr><td>15</td><td>1-p8</td><td>1+p7-p8</td><td>p9</td><td>28</td><td>1-p13+p14</td><td>1-p13</td><td>p15</td></tr><tr><td>16</td><td>-p7+p8</td><td>1-p7</td><td>p9</td><td>29</td><td>1-p13</td><td>1-p14</td><td>1-p15</td></tr><tr><td>17</td><td>1-p7</td><td>1-p8</td><td>1-p9</td><td>30</td><td>p14</td><td>1-p13+p14</td><td>1-p15</td></tr><tr><td>18</td><td>p8</td><td>-p7+p8</td><td>1-p9</td><td>31</td><td>p13-p14</td><td>p13</td><td>1-p15</td></tr><tr><td>19</td><td>1+p7-p8</td><td>p7</td><td>1-p9</td><td>32</td><td>p16</td><td>p17</td><td>p18</td></tr><tr><td>20</td><td>p10</td><td>p11</td><td>p12</td><td>33</td><td>1-p17</td><td>1+p16-p17</td><td>p18</td></tr><tr><td>21</td><td>1-p10</td><td>p10-p11</td><td>p12</td><td>34</td><td>-p16+p17</td><td>1-p16</td><td>p18</td></tr><tr><td>22</td><td>1-p10+p11</td><td>1-p10</td><td>p12</td><td>35</td><td>1-p16</td><td>1-p17</td><td>1-p18</td></tr><tr><td>23</td><td>1-p10</td><td>1-p11</td><td>1-p12</td><td>36</td><td>p17</td><td>-p16+p17</td><td>1-p18</td></tr><tr><td>24</td><td>p11</td><td>1-p10+p11</td><td>1-p12</td><td>37</td><td>1+p16-p17</td><td>p16</td><td>1-p18</td></tr><tr><td>25</td><td>p10-p11</td><td>p10</td><td>1-p12</td><td>38</td><td>-</td><td>-</td><td>p19</td></tr><tr><td>26</td><td>p13</td><td>p14</td><td>p15</td><td>39</td><td>-</td><td>-</td><td>1-p19</td></tr></tbody></table>

The sketch of the Ta plane, top left, expanded to show four SoDs, with Ta atoms in the model enumerated. Fractional coordinates of Ta atoms are described by 6 independent positional parameters, with symmetry constraints as tabulated in top right. Fractional coordinates of S atoms are described by 13 independent positional parameters, that were symmetry constrained as tabulated (bottom). The model does not account for the c-axis stacking correlations beyond the nearest neighbour layers and does not explain PDF data beyond the length-scale of ~1 nm. Since the fitting range is shorter than the unit cell size, lattice parameters (a=b, c) provide the distance metrics rather than describing translational symmetry. The parameters are labelled as $p_i$, where i=1-19 represent the index of each symmetry-independent parameter.
<figure><img src="image_12.png" /><figcaption>**Supplementary Figure 6.** Features of the P3/2 structure fit to the 15 K PDF data of 1T-TaS2 in the C state.</figcaption></figure>

The illustration focuses on two SoD motifs. Tantalum atoms are shown as larger gray spheres, whereas sulphur atoms are depicted as smaller coloured spheres, where colouring reflects their location with respect to the central Ta in SoDs. Moving away from the central Ta the TaS6 octahedra are progressively more distorted, as portrayed by the colour bar which could be considered as coarse indicator of strain associated with the SoD formation. Although our local model is relatively simple approximant that does not describe the long-range ordering, it does incorporate the degrees of freedom necessary to account for the modulation of sulphur sublattice, as shown in the bottom most panel where the sulphur atoms are represented by very small spheres to emphasize the effect seen in the fit. The solid lines are the guides to the eye and the modulation they depict is out of proportion. Such modulation is expected in the long range ordered structure and is enabled in our approximant model focused on the evolution of local Ta-Ta correlations.
<figure><img src="image_13.png" /><figcaption>Supplementary Figure 7. Considerations of the fitting range and the applicability of the $P\bar{3}$ approximant model.</figcaption></figure>

In order to establish suitable range for the local structure study, the model was tested against the 15 K dataset. In this, 19 structural parameters, described in Suppl. Fig. 6, were varied, while the lattice parameters, thermal parameters, and the overall scale were kept fixed. The fitting test was conducted over a variable range of data by incrementing the upper limit of refinement, $r_{MAX}$, by 1 Å, from 5 Å to 15 Å, as shown in the figure. In this process we monitored two quantities: the fit quality captured by the fit residual, $R_w$, and the number of parameter-pairs whose correlation is larger than $|0.8|$ as reported by the fitting program. In the former, large value of $R_w$ indicates inadequacy of the model used. In the latter, parameters exhibiting correlation $>|0.8|$ are typically considered to be unreliable since they are not independent when too narrow a data range is considered. Their report by the fitting platform indicates an over-parametrization of the model or, conversely, insufficient information in the data to constrain the model. In our test, fits of exceptionally good quality were obtained for $r_{MAX} \le 8$ Å, as shown by the RW(T) dependence (red symbols in the bottom panel, right ordinate) and illustrated in the top left panel. However, fits over these ranges result in heavily over-parametrized fits, as indicated by the temperature dependence of the number of correlated parameter pairs (blue symbols in the bottom panel, left ordinate). Increasing the fitting range results in parameters becoming uncorrelated, as sufficient amount of data is added to the fitting. While this initially
helps, with a moderate jump in R<sub>w</sub>, for r<sub>MAX</sub> ≥ 11 Å the fits begin to deteriorate, as illustrated in the upper right panel.
The reason behind this is that wider range of PDF data provides access to significant contributions of the NNN
interlayer correlations which are beyond the scope of our present model. We identified r<sub>MAX</sub> = 10 Å as a sweet spot
offering a reasonable compromise between the over-parametrization (model too complex for the data range) and
underfitting (model inadequacy). The result is a reliable determination of the local structure, which is the focus of
the present work.
## 2. ORDER PARAMETERS DESCRIBING THE SYMMETRY CHANGE LEADING TO THE COMMENSURATE CDW LOW TEMPERATURE PHASE OF 1T-TaS₂.

The analysis of the symmetry of the high temperature structure of 1T-TaS₂ can help understanding the structural changes of the low temperature daughter structures in terms of the symmetry adapted coordinates belonging to the irreducible representations of the parent space-group symmetry. The physical order parameters describing the change (in this case, the atomic displacements) can be classified as basis functions of the irreducible representations of the parent symmetry group. In the current case, multiple order parameters contribute to the final observed distortion: the individual degrees of freedom provided by an instance of these multi-dimensional order parameters are physically meaningful because they provide a natural description of the distortions of the parent structure.

Lowering the temperature, some of the symmetry elements of the symmetry group of the parent phase are lost (e.g the mirror planes), but the surviving symmetries that describe the daughter phase are a subgroup of the parent phase that is uniquely identified by its combination with the new superlattice basis (and eventual new supercell origin choice). Using a frozen-phonon picture, the structural parameters describing the low temperature daughter structures can be organized as a superposition of polarization vectors containing the cartesian components of the displacements of each one of the independent atoms of the daughter phase.

### Revisiting the 1T (single layer) C-CDW: atomic vibrations

In the following, we totally neglect any c-stack rule and consider the symmetry changes occurring in a single layer of the TaS₂ structure. We start from the standard description of the prototype phase that has space group symmetry P-3m1 (#164, $D_{3d}^3$, $a_p=b_p\approx3.34$ Å, $c_p\approx5.894$ Å, $\gamma=120^\circ$). In this structure, there are two independent atoms with positions described by the coordinates summarized in the table below:

<table><thead><tr><th>atom</th><th>site</th><th>x</th><th>y</th><th>z</th><th>sym</th></tr></thead><tbody><tr><td>TaO</td><td>1a</td><td>0</td><td>0</td><td>0</td><td>D3d</td></tr><tr><td>SO</td><td>2d</td><td>1/3</td><td>2/3</td><td>0.22074</td><td>C3v</td></tr></tbody></table>
The mechanical representation Γ at the Brillouin zone centre of this prototype structure is reducible and it consists of 2 acoustic and 4 optic mode frequencies as shown by the decomposition into irreducible representations.

$$ \Gamma = (A_{2u} \oplus E_u) \oplus (A_{1g} \oplus A_{2u} \oplus E_g \oplus E_u) $$

In particular, the acoustic phonons of this prototype structure belong to the ($A_{2u} \oplus E_u$) irreducible representations and involve the displacements of the Ta atoms. The $A_{2u}$ irrep describes a vibration of the Ta atom along z, while the $E_u$ symmetry adapted displacements describe degenerate vibrations in the xy plane.

## The low T variants

All the low temperature phases have axes rotated when compared to the original frame of the high temperature prototype structure. The lattice changes produced by rotation around the z-axis (which has $A_{2g}$ representation) can maintain the -3 symmetry operations, but they systematically remove the mirrors and 2-fold axes of the parent group. The low temperature phases also involve superlattice reflections and therefore require the definition of a superlattice cell in real space. The superlattice reflections characteristic of the low temperature variants (IC, NC, and C-CDW phases) involve a condensation of Bragg intensity at some q-points in the high temperature Brillouin zone.

In the simplest description of the low temperature variants (incommensurate or commensurate), two additional independent vectors of the high temperature Brillouin zone are required. The low temperature phase (C-CDW) is obtained by the following matrix transformation of the prototype lattice

$$ (a \ b \ c) = (a_P \ b_P \ c_P) \begin{pmatrix} 4 & -1 & 0 \\ 1 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix} $$

In this reference cell, there are 13 replicas of the $TaS_2$ motif. The atomic positions of the prototype structure can be described in this new metric that has symmetry P-3 (#147, a=b≈12.04254 Å, c≈5.894 Å, γ=120°) by a new set of independent atomic positions summarized in the table below:
<table><thead><tr><th>atom site</th><th>x</th><th>y</th><th>z</th></tr></thead><tbody><tr><td>Ta1</td><td>1a</td><td>0.</td><td>0.</td><td>0.</td></tr><tr><td>Ta2</td><td>6g</td><td>10/13</td><td>1/13</td><td>0.00000</td></tr><tr><td>Ta3</td><td>6g</td><td>11/13</td><td>5/13</td><td>0.00000</td></tr><tr><td>S1</td><td>6g</td><td>-1/3+6/13</td><td>1/3-2/13</td><td>0.22074</td></tr><tr><td>S2</td><td>6g</td><td>2/3+3/13</td><td>1/3-1/13</td><td>0.22074</td></tr><tr><td>S3</td><td>6g</td><td>2/3+4/13</td><td>1/3+3/13</td><td>0.22074</td></tr><tr><td>S4</td><td>6g</td><td>-1/3+7/13</td><td>1/3+2/13</td><td>0.22074</td></tr><tr><td>S5</td><td>2d</td><td>2/3</td><td>1/3</td><td>0.22074</td></tr></tbody></table>

This phase corresponds to the lock-in of the components of two vectors $q_1$ and $q_2$; these reciprocal space vectors have generic components (ab0) and form an angle of 60° to maintain the symmetry of the threefold axis. Their components have the form $q_1 = (a, b, 0)$ and $q_2 = (C_{-3+}q_1) = (-b, a+b, 0)$. At the lock-in, described by the lattice change transformation matrix, $4q_1-q_2=a^*$ and $q_1+3q_2=b^*$ so that $q_1=3/13 a^*+1/13 b^*$ and $q_2=-1/13 a^*+4/13 b^*$.

Since there is no component normal to the layer, the $q_1$ vector defines a star of 6 vectors (including $q_2$) related by the symmetry operations of the parent space. The same considerations actually apply also to the symmetry analysis of the incommensurate variant studied by Spijkerman et al. In that case, the components a and b of the vectors are incommensurate but not far from their lock-in values:

$$a_{IC}=0.2448(2) -> a_{C-CDW}=0.2308=3/13$$

$$b_{IC}=0.0681(2) -> b_{C-CDW}=0.0769=1/13$$

As we can see, at the lock-in, the $a_{IC}$ component has decreased and $b_{IC}$ has increased, so the $q_1$ and $q_2$ vectors orientations relative to the parent phase have rotated while the temperature was lowered. A
commensurate approximant of the incommensurate phase with P-3 symmetry could be constructed
eventually taking a=1/4=0.25 and b=1/16=0.0625, or some other convenient rational approximation of
those values.

To calculate the frozen phonons displacements in the C-CDW phase, we have to consider two stars of
vectors. The first star is the one already discussed, originating from q₁ (q₂ obviously belongs to this star by
construction) and represented by blue arrows in Supplementary Figure 9. The second star of vectors is
the one originating from q₁+q₂ with components (2/13 5/13 0) that indexes the 2nd order superlattice
reflections of this phase, so that all the 12 extra spots of the reciprocal lattice of the C-CDW phase can be
indexed. Each one of these two stars has in principle 12 branches and so the dimension of the Order
Parameter for each one of the two stars. The 12 branches of the star of these vectors in P-3m1 have
components E(a,b,0), i(-a, -b, 0), -3⁺(-b,a+b,0), 3⁺(b, -a-b, 0), -3⁻(a+b,-a,0), 3⁻(-a-b, a, 0), m₁₀₀(-a,a+b,0),
2₁₀₀(a, -a-b, 0), m₀₁₀(a+b,-b,0), 2₀₁₀(-a-b, b, 0), 2₁₁₀(b,a,0), and m₁₁₀(-b, -a, 0).

<figure><img src="image_14.png" /><figcaption>Supplementary Figure 9. Representation of the two stars indexing the superlattice reflections of the C-CDW phase.</figcaption></figure>

Elements of the q₁ star are displayed as blue vectors, elements of the q₁+q₂ star are shown as green vectors. The shaded area is the reciprocal cell of the prototype phase.

For the q₁ star we can simply take the spots indexed by the elements (E, -3⁺ and -3⁻) with explicit
components (3/13 1/13 0) (-1/13 4/13 0) (4/13 -3/13 0): therefore, the OP describing the superpositions
of the atomic displacements has components (a0a0a0000000). The same choice can be made for the q₁+q₂
star, or explicitly (2/13 5/13 0) (-5/13 7/13 0) (7/13 -2/13 0). The OP describing the transformation from
the parent phase to the C-CDW lock-in phase has overall symmetry components
(a0a0a0000000;b0b0b0000000). It is then possible to decompose the atomic displacements of the Ta
atoms at any temperature in terms of the amplitudes of the E<sub>u</sub> and A<sub>2u</sub> modes transported at q<sub>1</sub> and at
q<sub>1</sub>+q<sub>2</sub>.

Symmetry adapted coordinates of Ta displacements at q₁ and at q₁+q₂.

These displacements of the Ta atoms are particularly relevant for the description of the dimerisation of
Ta distances that is a characteristic signature of the C-CDW phase.

Eu (q₁) 1st component – phases and relative displacement amplitudes along x,y for Ta2 and Ta3 atoms

Ta1 0.0000 0.0000 0.0000

Ta2 -0.9708 -1.0000 0.0000

Ta3 -0.2939 -0.3028 0.0000

Eu (q₁) 2<sup>nd</sup> component – phases and relative displacement amplitudes along x,y for Ta2 and Ta3 atoms

Ta1 0.0000 0.0000 0.0000

Ta2 -0.5222 0.4778 0.0000

Ta3 -0.1581 0.1447 0.0000

Eu (q₁+q₂) 1st component – phases and relative displacement amplitudes along x,y for Ta2 and Ta3 atoms

Ta1 0.0000 0.0000 0.0000

Ta2 -0.0088 0.2939 0.0000

Ta3 0.0292 -0.9708 0.0000

Eu (q₁+q₂) 2<sup>nd</sup> component – phases and relative displacement amplitudes along x,y for Ta2 and Ta3 atoms

Ta1 0.0000 0.0000 0.0000

Ta2 -0.3028 -0.1581 0.0000

Ta3 1.0000 0.5222 0.0000

A2u (q₁) – phases and relative displacement amplitudes along z for Ta2 and Ta3 atoms

Ta1 0.0000 0.0000 0.0000

Ta2 0.0000 0.0000 0.3028

Ta3 0.0000 0.0000 1.0000
A2u ($q_1+q_2$) – phases and relative displacement amplitudes along z for Ta2 and Ta3 atoms

Ta1 0.0000 0.0000 0.0000

Ta2 0.0000 0.0000 -1.0000

Ta3 0.0000 0.0000 0.3028

This analysis can be also extended to S atoms, but:

- more components are involved, reducing the value of the decomposition procedure

- optic phonons do not seem particularly relevant for this phase transition mechanism

## Cartesian components of the displacements of each independent Ta atom in the daughter phase

The supplementary table 2 summarizes the fractional displacement contributions originating from the amplitudes transported at $q_1$ and $q_1+q_2$. These displacements add to the undistorted positions of the Ta2 and Ta3 atoms inherited from the high temperature parent phase. The displacements of the Ta1 atom at the origin of the lattice are fixed by the symmetry element of the daughter phase, so they are not reported in the table. Some trends can be remarked:

- The out-of-plane z component decreases with decreasing T, resulting in puckering of the Ta atoms

- The Eu $q_1$ component is more important for Ta2 atoms whereas the Eu ($q_1+q_2$) component is the dominant displacement component for Ta3 atoms.

<table><thead><tr><th></th><th>undistorted</th><th>q1</th><th>175K</th><th>150K</th><th>125K</th><th>100K</th><th>50K</th><th>15K</th></tr></thead><tbody><tr><td>Ta2_x</td><td>0.07692</td><td>Eu Ta2-cl</td><td>-0.00957</td><td>-0.00957</td><td>-0.00898</td><td>-0.00903</td><td>-0.00838</td><td>-0.00779</td></tr><tr><td>Ta2_y</td><td>0.30769</td><td>Eu Ta2-dy</td><td>-0.01577</td><td>-0.01576</td><td>-0.01536</td><td>-0.01541</td><td>-0.01548</td><td>-0.01529</td></tr><tr><td>Ta2_z</td><td>0</td><td>A2u Ta2-dz</td><td>-0.00219</td><td>-0.00210</td><td>-0.00151</td><td>-0.00162</td><td>-0.00162</td><td>-0.00169</td></tr><tr><td>Ta3_x</td><td>0.15385</td><td>Eu Ta3-cl</td><td>-0.00188</td><td>-0.00188</td><td>-0.00193</td><td>-0.00193</td><td>-0.00215</td><td>-0.00227</td></tr><tr><td>Ta3_y</td><td>0.61538</td><td>Eu Ta3-dy</td><td>0.00291</td><td>0.00290</td><td>0.00272</td><td>0.00273</td><td>0.00254</td><td>0.00236</td></tr><tr><td>Ta3_z</td><td>0</td><td>A2u Ta3-dz</td><td>-0.00727</td><td>-0.00695</td><td>-0.00498</td><td>-0.00535</td><td>-0.00536</td><td>-0.00560</td></tr><tr><td></td><td></td><td>q1+q2</td><td>175K</td><td>150K</td><td>125K</td><td>100K</td><td>50K</td><td>15K</td></tr><tr><td></td><td></td><td>Eu Ta2-cl</td><td>0.00311</td><td>0.00320</td><td>0.00341</td><td>0.00359</td><td>0.00392</td><td>0.00410</td></tr></tbody></table>
<table><tr><td>Eu</td><td>Ta2DY</td><td>0.00405</td><td>0.00404</td><td>0.00443</td><td>0.00445</td><td>0.00446</td><td>0.00446</td></tr><tr><td>A2u</td><td>Ta2_DZ</td><td>0.00894</td><td>0.00877</td><td>0.00797</td><td>0.00755</td><td>0.00722</td><td>0.00731</td></tr><tr><td>Eu</td><td>Ta3DX</td><td>-0.00313</td><td>-0.00278</td><td>-0.00337</td><td>-0.00285</td><td>-0.00178</td><td>-0.00118</td></tr><tr><td>Eu</td><td>Ta3DY</td><td>0.01027</td><td>0.01056</td><td>0.01125</td><td>0.01185</td><td>0.01295</td><td>0.01355</td></tr><tr><td>A2u</td><td>Ta3_DZ</td><td>-0.00271</td><td>-0.00265</td><td>-0.00241</td><td>-0.00229</td><td>-0.00218</td><td>-0.00221</td></tr></table>

Supplementary table 2. Fractional displacement contributions originating from the amplitudes transported at $q_1$ and $q_1+q_2$.

[1] O. Sezerman, A. M. Simpson, and M. H. Jericho, Solid State Commun 36, 737 (1980).