Slow dynamics of disordered zigzag chain molecules in layered LiVS2 under electron irradiation

Electronic instabilities in transition metal compounds often spontaneously form orbital molecules, which consist of orbital-coupled metal ions at low temperature. Recent local structural studies utilizing the pair distribution function revealed that preformed orbital molecules appear disordered even in the high-temperature paramagnetic phase. However, it is unclear whether preformed orbital molecules are dynamic or static. Here, we provide clear experimental evidence of the slow dynamics of disordered orbital molecules realized in the high-temperature paramagnetic phase of LiVS2, which exhibits vanadium trimerization upon cooling below 314 K. Unexpectedly, the preformed orbital molecules appear as a disordered zigzag chain that fluctuate in both time and space under electron irradiation. Our findings should advance studies on soft matter physics realized in an inorganic material due to disordered orbital molecules.


INTRODUCTION
Transition metal compounds with multiple electron degrees of freedom frequently form self-organized molecules in low temperature, called orbital molecules [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] . Examples include octamers in CuIr 2 S 4 1 and trimers in LiVO 2 and LiVS 2 5,6 . Recent structural studies clarified that the orbital molecules at low temperature can persist in the high-temperature paramagnetic phase as a disordered form in some substances [9][10][11][12][13][14][15][16] . Persistent orbital molecules contradict the general belief that the regular lattice should recover in the high-temperature paramagnetic phase, as speculated based on the average structure obtained from conventional diffraction experiments [1][2][3][4][5][6][7][8][9][10] . Because the long-range ordering of the arrangements of orbital molecules are restored upon cooling in these substances, it is expected that the persistent orbital molecules realize lattice dynamics. Kimber et al. predicted that a classical analog of the short-range Resonating Valence Bond state will be realized in the high-temperature phase of Li 2 RuO 3 as the existing dimer patterns resonate due to thermal fluctuations 13 . They explained that slower dynamics should appear compared to the characteristic timescale of their X-ray measurements, although these dynamics have yet to be experimentally clarified. By contrast, Browne et al. performed a quasi-elastic neutron scattering experiment up to 1100 K in GaV 2 O 4 with disordered trimer and tetramer pairs in anticipation of the appearance of fast dynamics. However, they concluded that clusters are well-defined and statically disordered even at 1100 K as they were unable to detect dynamic behaviors faster thañ 1 × 10 −11 s 12 .
Here, we present comprehensive structural studies of LiVS 2 , which exhibits a paramagnetic metal to nonmagnetic insulator transition at 314 K 5 . A recent structural study clarified that vanadium trimer molecules form in the low-temperature nonmagnetic phase 6 . The vanadium trimer molecules disappear above 314 K, and zigzag chain molecules consisting of dimers emerge with a finite correlation length. Cooling increases the correlation length, and sharp superstructure peaks grow in powder diffraction patterns below approximately 350 K.
Our annular dark-field scanning transmission electron microscope (ADF-STEM) experiment clearly shows that these existing zigzag chain molecules slowly fluctuate in both time and space on the order of seconds. The high-temperature phase with slow dynamics of some existing zigzag chain patterns is categorized as a plastic crystal phase, which is a mesophase between a crystal and liquid. The atoms retain the original position in the time average but the orientation of the zigzag chain fluctuates. Figure 1a displays the physical properties of LiVS 2 . It exhibits a metal to nonmagnetic insulator transition around 314 K. As a previous study has already clarified 5 , the high-temperature paramagnetic susceptibility increases upon heating, which is reminiscent of pseudogap behavior found in underdoped high-T c cuprates. Although the present data are consistent with those reported previously, the present entropy change of ΔS = 7.99 J/ mol K at 314 K is much larger than ΔS = 6.6 J/mol K in the previous study 5 . This difference may be due to the improved sample quality in this study. For checking the quality of the present sample, we performed the Inductively Coupled Plasma-Atomic Emission Spectrometry (ICP-AES) experiment, and confirmed that the Li content, x, is 0.97 (2) in Li x VS 2 . A subtle amount of off-stoichiometry does not have a significant effect on the physical properties, as mentioned above. Note that no discontinuous behaviors appear in these physical property data above 314 K for the discussion later. Figure 1b shows the synchrotron X-ray diffraction patterns as functions of temperature. Below 314 K, sharp superstructure peaks emerge around Q~3.2-3.3 Å −1 . Recently, the low-temperature crystal structure was solved using a Rietveld method by assuming a trigonal space group P31m. It was clarified that the vanadium trimers in a spin singlet state spontaneously form (Fig. 1c) 6 .

Rietveld analysis
Upon heating above 314 K, the superstructure peaks associated with vanadium trimerization completely disappear in the powder diffraction data and additional superstructure peaks emerge (Fig.  1b, arrows), which has never been reported. The additional peaks gradually weaken upon further heating and eventually disappear at approximately 350 K. For the higher temperature data, refinement can be successfully performed by assuming a trigonal space group P3m1 with a regular triangular lattice, while the refinement can be successfully performed for powder diffraction data obtained at 320 K by assuming a monoclinic space group Pm with vanadium zigzag chain molecules (Fig. 1d). Strangely, the distinct superstructure peaks clearly show the monoclinic displacement of the consisting atoms, the monoclinic lattice strain from the triangular lattice is negligible. Both the temperature dependent lattice parameters and the details of the Rietveld refinement are summarized in the Supplementary Note 1. Figure  1e-g display the high-resolution X-ray diffraction patterns at 320 K with the indices. The emergences of the superstructure peaks with l ≠ 0 in the indices indicate that the zigzag chain appears as threedimensional ordering.

PDF analysis
Pair Distribution Function (PDF) analysis was performed to evaluate whether the preformed orbital molecules appear in a disordered form in the high-temperature paramagnetic phase. Figure 2a displays the reduced PDF data in the range of 2 ≤ r(Å) ≤ 10 obtained at 325 K, where the averaged Pm structure is realized, although the monoclinic lattice strain from the triangular lattice is negligibly small. At first, the fit assuming the average P31m model with the vanadium trimers is poor (R w = 18.6%) because the simulation generates an unrealized peak at r~3.0 Å, which corresponds to the inner-trimer V-V distance. This clearly indicates that disordered vanadium trimers are absent in the hightemperature paramagnetic phase. Secondly, the fit assuming the average P3m1 model with the regular vanadium lattice is also poor (R w = 14.6%). The peak at r~3.4 Å, which corresponds to adjacent V-V and S-S distances, is much broader than the peak at r~2.5 Å, which consists of the adjacent V-S distance. We cannot successfully fit both peaks by assuming the average P3m1 model even if we assume anomalous thermal parameters for V and S ions. If the average P3m1 model is realized, the peak at r~3.4 Å should be stronger and sharper. Therefore, we can also exclude the regular triangle from the candidates. The experiment confirmed that the monoclinic Pm crystal structure with a vanadium zigzag chain gives the best fit (R w = 5.37%), consisting with the Rietveld analysis results. Unexpectedly, sharply different from the powder diffraction patterns shown in Fig. 1b, where the superstructure peaks disappear at approximately 350 K, the peak profiles in the range of 2 ≤ r (Å) ≤ 10 are maintained without any significant changes upon heating at least up to 460 K (Fig. 2b). Hence, the average Pm model can still be fitted well to the reduced PDF patterns even at 460 K. The V-V distances obtained from the refined structures using Rietveld and the PDF analysis in the range of 2 ≤ r(Å) ≤ 10 exhibit significant differences (Fig. 1a, bottom). We can successfully explain it in terms of the correlation length of the Pm monoclinic domain, as described below.
The displacement of V obtained from Rietveld refinement is always smaller than that obtained from PDF analysis and decreases with increasing temperatures. Thus, the correlation length of the monoclinic domains with zigzag chain molecules has a finite value above 314 K and decreases upon heating. It should be noted that the superstructure peaks are as sharp as the fundamental peaks at 320 K, as seen in Fig. 1f, g. This means that the width of the superstructure peaks is limited by the instrumental width. From Scherrer's equation, we can roughly estimate the correlation length of monoclinic domains to be . Open circles are derived by analyzing the reduced PDF data in the range of 2 ≤ r(Å) ≤ 10, while the closed circles are from the Rietveld analysis of the powder diffraction data. Insets schematically depict the averaged structure obtained from Rietveld analysis in each temperature range. b Temperature dependences of powder X-ray diffraction data. Arrows indicate peaks that appear just above the trimer transition temperature of 314 K, which weaken upon further heating and disappear around 350 K. Detailed results of Rietveld analysis for each temperature regions were summarized in Supplementary Note 1. Inset shows the temperature dependences on superstructure peak intensity. c Schematic of vanadium trimers realized below 314 K. d Schematic of vanadium zigzag chain appearing with a finite correlation length in the high-temperature paramagnetic phase. e High resolution X-ray diffraction patterns at 320 K with indices of monoclinic phase. f, g Enlarged views of a fundamental 001 peak (f) and a superstructure 103 peak (g) with the correlation length estimated using the Scherrer equation. Instrumental resolution is not taken into account. These data are obtained from the high-resolution data.
longer than~2200 Å at 320 K. Although the superstructure peaks disappear above approximately 350 K in the powder diffraction data, the monoclinic domains persist up to at least 460 K (Fig. 2b, reduced PDF data). Because the superstructure peaks disappear in the continuous process where the correlation length shortens upon heating, anomalies do not appear in the physical property data around 350 K (Fig. 1a). As shown in Fig. 2c, the overall appearances of PDF data collected at 325 K and 370 K are quite similar between them even in the range of 90 ≤ r(Å) ≤ 100, and refinement is successfully performed by assuming the Pm crystal structure with a zigzag chain, indicating the long correlation length over 90 Å is maintained even at 370 K. Further discussion about the PDF analysis at high r regions is available at the Supplementary Note 1 of this article.

STEM experiment
The strong temperature dependence on correlation length expects us that the zigzag chain molecules have dynamic properties. To clarify whether the preformed zigzag chain molecules are dynamic or static, we observed a time series of high-resolution annular dark-field scanning transmission electron microscope (ADF-STEM) images of LiVS 2 at room temperature. In the thick part of the piece away from the edge, such as in the lower left portion of the inset of Fig. 3a, the diffraction pattern derived from the P31m phase having a trimer was clearly observed. The obtained diffraction pattern is similar with that previously reported 5 . On the other hand, unexpectedly, a diffraction pattern derived from the monoclinic Pm phase with a zigzag chain was observed despite the measurement at room temperature in the region near the edge of the sample. The ADF-STEM experimental results on this edge region will be presented below. Figure 3a (1024 × 1024 pixels) shows a high-resolution ADF-STEM image observed by [0001] incidence on the edge of the specimen. Image areas of approximately 130 Å × 130 Å were acquired with a pixel dwell time of 1 µs, which is a time interval of approximately 1 s. STEM images were taken continuously, and time series were acquired at intervals of almost 1 s. The bright and dark dots correspond to V and S atom columns projected along the [0001] orientation, respectively. The V and S atoms form a twodimensional triangular lattice. The Fourier transform pattern of the ADF-STEM image shows superlattice spots at positions ±1/2000, 0 ± 1/200, and ∓ 1=2 ± 1=200 for high temperature trigonal unit cell (Fig. 3b). This indicates the presence of three monoclinic domains, which are rotated by 120 degrees from each other (Fig.  3c). Figure 3d-l show a time series of Fourier-masked ADF-STEM images. The inverse FFT is performed using the superlattice spots of ±1/2000, 0 ± 1/200, and ∓ 1=2 ± 1=200, and the real space distribution of each domain is colored by red, green and blue. The intensity of the three types of lattice fringes is modulated in different ways to form a unique domain structure. Note that these images show the structures projected in the beam direction, where domain patterns appearing in several VS 2 layers are superimposed. Remarkably, the modulation and domain patterns change in every frame of the time series of STEM images, demonstrating that the zigzag chain molecules fluctuate on a timescale of seconds (The corresponding movie file is available). The slow dynamics on a timescale of seconds clearly indicates that the present state is caused by thermal fluctuations.
Our STEM experimental results raise the question of why the Pm phase, which should occur above 314 K, occurs in the edge region at room temperature. There are several possible scenarios. First, it is possible that the sample temperature rose due to electron beam irradiation and exceeded the phase transition temperature of 314 K. But this is probably not happening. When a carbon thin film is irradiated with a 5 nA beam with a diameter of 1 nm, the temperature rise is reported to be 1.4 K 17 . In this experiment, the beam diameter was about 0.1 nm and the electron beam current was 0.1 nA. This indicates that the sample has a much lower electron dose per unit time. Furthermore, this scenario cannot be explained that the P31m phase with vanadium trimers was also observed in the region away from the edge under the same experimental conditions. The second scenario concerns Li being unexpectedly missing at the edge of the sample. In the Li x VS 2 system, it has been reported that the trimer transition temperature decreases as the Li content x decreases from 1.0, and the Fig. 2 Reduced pair distribution function (PDF) analysis of LiVS 2 in the high-temperature paramagnetic phase. a Fitting (red) and residual (green) for reduced PDF data (blue) obtained at 325 K. Average P31m model with V trimers (top) generates unrealized peak around 3.0 Å, resulting in a poor fit. The P3m1 model with regular triangular lattice cannot be fitted well to the intensity data (middle). Best fit is obtained by assuming the Pm structure with a zigzag chain (bottom). Adjacent V-V distances can be separated into three types, which are drawn in blue, red, and green in the Pm structure. Blue is the shortest and forms a zigzag chain. Initial structural model is obtained from the Rietveld refinement of the powder diffraction experiment at the corresponding temperature. b Temperature dependences of PDF data in the hightemperature paramagnetic phase. Increasing the temperature up to 460 K does not cause a remarkable change in the PDF patterns, indicating that the zigzag chain model can be applied up to at least 460 K. c PDF analysis in the range of 90 ≤ r(Å) ≤ 100 by assuming the Pm structure both for 325 K and 370 K data. The reliable factors are R w = 6.04% (325 K) and 5.29% (370 K), respectively. trimer transition temperature falls below room temperature when x = 0.8. Note that our X-ray diffraction experimental results of the Li deficient samples are consistent with this tendency. The superlattice peaks survive at 300 K for x = 0.91(1), while they disappear at~250 K for x = 0.78(1) (Supplementary Note 1). When the transition temperature is significantly lowered from room temperature due to Li deficiency, the correlation length of the Pm phase at room temperature may be sharply shortened, and multiple domains may be superimposed in this ADF-STEM image at room temperature. Another possibility within this scenario is that the Li at the edges is strongly diffused at room temperature due to the thermal and/or electron irradiation effects. A third possible scenario is that the phase transition temperature can change with sample thickness, as is often the case with some transition metal dichalcogenides 18,19 . For example, in 1T-TaS 2 , the low temperature Commensurate Charge Density Wave (CCDW) phase is suppressed when the sample thickness is below the threshold of~40 nm 19 . Layer thicknesses at the edge of current samples are expected to be up to a few nanometers since it is not possible to perform ADF-STEM experiments on thicker samples. A fourth possible scenario is that the anomalous charge-fluctuating phase accompanied by lattice dynamics is induced by electron irradiation. A reference phenomenon has been observed in Ba 3 NaRu 2 O 9 , which contains the Ru 2 O 9 dimer, where the charge order melts under photoexcitation 20 . Further experimental studies should be required to clarify which scenario is realized in the future.

Theoretical considerations
After identifying the dynamics of disordered zigzag chain molecules appearing in the high-temperature paramagnetic phase of LiVS 2 , we considered the underlying physics generating such dynamics of disordered zigzag chain molecules. Figure 4a shows the band calculation result based on the parameters obtained from the Rietveld analysis of 360 K data by assuming the trigonal space group P3m1. The triply degenerate t 2g orbitals, which consist of xy, yz, and xz orbitals, are oriented toward the neighboring vanadium ions and are inherently separated into lower doubly degenerate e′ g orbitals and a higher nondegenerate a 1g orbital due to the trigonal splitting (Fig. 4d).
The corresponding first principles calculation result for a P3m1 structure with U = 4 eV suggests a Fermi surface with the possible common nesting vector of Q = a * /2 (Fig. 4e), which can cause multiple dimerizations of the lattice along the a-direction and (a + b)-direction. As shown in Fig. 4d, this leads to the zigzag chain pattern. Our band calculation result obtained using the parameters of the local Pm structure at 360 K clearly shows that the d yz + d xz component, which is due to the doubly degenerate e′ g bands in the high-temperature P3m1 phase, is separated into bonding and antibonding bands. Electrons fill the bonding bands (Fig. 4b). The Fermi surface survives due to the remnant d xy component even when the calculation assumes the Pm structure consisting of a metallic nature in high-temperature paramagnetic phase. We naively insist that the persisting Fermi surface in the Pm phase weakens the band energy stabilization associated with the zigzag chain formation as it may be related to the emergence of the observed slow dynamics. When vanadium trimers are formed in the lower temperature region, a large band gap consisting of the insulative behavior in an electrical resistivity experiment (Fig. 4c). We speculate that the transition at 314 K induces a large band energy at the expense of lattice energy.

DISCUSSION
As discussed above, the short-range order of the zigzag chains appearing in the high temperature phase of LiVS 2 can be originating from the orbitally assisted CDW instability. C. Rovira and M.-H. Whangbo discuss that both zigzag chains and trimers are explained on the basis of the concept of both local chemical bonds and hidden Fermi surface nesting 21 : When t 2g orbital is triply degenerated and the proportionate band filling of (d xy , d yz , d zx ) = (2/3, 2/3, 2/3) is realized, the trimerization appears as a result of the Fermi surface nesting. On the other hand, when the disproportionate band-filling of (d xy , d yz , d zx ) = (1,1,0) is realized, the zigzag chain pattern appears as the ground state. Based on this argument, one may consider that the phase transition at 314 K in LiVS 2 originates from two energetically competing CDW states. However, the situation is not so simple. First, in the high temperature phase of LiVS 2 , the triply degeneracy t 2g orbital is not maintained, and the nesting instability that causes zigzag chain order is inherently present. In other words, the zigzag chain order can be formed without causing band disproportionation. Second, there is no band instability that stabilizes trimerization. This indicates that the phase transition at 314 K is accompanied by another ingredient which changes band structure, such as the band Jahn-Teller transition.
Both trimers and zigzag chains are composed of vanadium ions with two bonds, attributed by the V 3+ d 2 electronic state of LiVS 2 . Since the zigzag chain has a finite correlation length in the high temperature phase, vanadium with only one bond must be present at the end of the zigzag chain. This may cause pseudogaplike behavior of magnetic susceptibility temperature dependence. When bond formation occurs, spin singlets are formed inside the bond, which makes it essentially nonmagnetic. However, since the vanadium at the end of the chain can form only one bond, the extra electrons should contribute to the paramagnetism. As the correlation length of the zigzag chain becomes shorter upon heating, the number of vanadium at the edge increases, should be leading to the enhancement of Pauli paramagnetism upon heating, as shown in Fig. 1a. Previous studies clarified that the enhancement of paramnagnetic susceptibility continues up to at least 700 K, possibly indicating that the short-range order survives even 700 K.
Although our theoretical calculations indicate that the Fermi surface instability may be critical to form the zigzag chain patterns, the origin of the dynamics is unclear. Considering that the disordered zigzag chain consists of the orbital-coupled metal ions, the unconventional interaction between orbital and phonon should be an ingredient for generating the unusual lattice dynamics. However, the slow dynamics on a second timescale indicate that the fluctuation itself is not of an electronic origin. Phenomenologically, it is reminiscent of the physics of soft materials such as a polymer rotation. The high-temperature phase of LiVS 2 with slow dynamics of zigzag chain patterns does not belong to a crystal where both the center of gravity and the orientation of the atoms/molecules should be maintained. Considering that the atoms retain the original position in a time average but the zigzag chain orientation fluctuates over time, we categorize it as a plastic crystal, which is a mesophase 22,23 .
Conventionally, plastic crystals are realized in weakly interacting molecules or ions, where the consisting molecules/ions are thermally rotating at a fixed position. The dynamics of conventional plastic crystals are usually studied via NMR techniques. Considering that NMR covers the kHz-MHz order, we can estimate that the zigzag chain dynamics realized in the present LiVS 2 should be some order of magnitude slower than the rotation dynamics observed in conventional plastic crystals. We speculate that this is due to the strong interaction among neighboring atoms, which depend on the complex network structure of inorganic materials. Our findings should provide a new platform for investigating soft matter physics in inorganic materials, as well as expand the fundamental understanding of systems with preformed orbital molecules at high temperature.
In summary, we first observed the slow dynamics of disordered orbital molecules appearing in the high-temperature paramagnetic phase of LiVS 2 . The unconventional coupling between orbital and phonon should be an ingredient for generating the dynamics, surely impact the studies of conventional orbital physics, such as nematic state in iron selenides [24][25][26] . Motivated by this work, we expect that further explorations of dynamics targeting similar systems, such as Li 2 RuO 3 and AlV 2 O 4 , should be accelerated, leading to the research fields of disordered orbital molecules.

Sample growth and preparation
Powder samples of LiVS 2 were prepared by a soft-chemical method and a subsequent solid-state reaction. Li-deficient Li ∼0.75 VS 2 was synthesized initially by a reaction with an appropriate amount of Li 2 S, V, and S in an Arfilled quartz tube at 700°C for 3 days. Note that the obtained precursors include small single crystalline samples with the order of μm. The products were put in 0.2 M n-BuLi hexane solution for 2 days to attain the maximum Li content 27 . The Li content, x, was confirmed to be 0.97 (2) in Li x VS2 from the Inductively Coupled Plasma-Atomic Emission Spectrometry (ICP-AES) using Hitachi SPCTRO ARCOS, equipped at the Institute of Solid State Physics (ISSP), Japan. The samples were crashed well to make powder sample to perform the X-ray diffraction experiments. In a STEM experiment, the sample of the same batch was diffused in anhydrous b Orbital-decomposed PDOS calculated for the local Pm structure. c Orbital-decomposed PDOS calculated for the trimer structure. d Orbitalassisted multiple dimerizations consisting of zigzag chain patterns. e Fermi surfaces calculated using the GGA+U scheme at U = 4 eV.
hexane and roughly ground for about 30 s. The supernatant was placed on a sample holder and used for measurement.

Physical property measurements
Differential scanning calorimetry (DSC) was conducted using a DSC 204 F1 Phoenix (Netzsch). The magnetic susceptibility was measured by a SQUID magnetometer (Quantum Design). The electrical resistivity was measured by the four-probe method. The powder samples were sintered at a low temperature of 300°C because the inserted Li ions were partially deintercalated at higher temperatures. Experiments were performed in an Ar atmosphere.
To investigate the average structure, synchrotron powder X-ray diffraction experiments were performed with the incident X-ray energy of E = 19 keV in BL5S2 beamline at Aichi Synchrotron, Japan, and with the incident X-ray energy of E = 30 keV in BL02B2 beamline at SPring-8, Japan. While Fig. 1b was constructed by the data obtained in BL5S2, the data of Fig. 1e, f was collected in BL02B2 beamline. A two-dimensional semiconductor detector, PILATUS 100 K, was used for high-resolution measurements and high-speed data collection in the experiments performed in BL5S2 beamline, while the six MYTHEN detectors were used in BL02B2 beamline. RIETAN-FP software 28 was employed for the Rietveld analysis, while VESTA was used for graphical purposes 29 . High-energy synchrotron x-ray diffraction experiments with E = 61 keV were performed at BL04B2 at SPring-8, Japan, to study the local structure from PDF analysis. Hybrid point detectors of Ge and CdTe were employed for data collection. The reduced PDF G(r) was obtained via the conventional Fourier transform of the collected data 30 . The PDFgui package 31 was used to analyze G(r). Capillaries containing samples and empty capillaries were measured under similar optical conditions, and the measurement results of the empty capillaries were subtracted from the data as the background. Polarization correction, Compton correction, and absorption correction were applied to obtain the reduced PDF data. For making reduced PDF data, the diffraction data with the Q range of 0.2 ≤ Q (Å −1 ) ≤ 25.7 with the ΔQ = 0.011 were used. Delta 1 and delta 2 options were not utilized in PDFGui software.

Scanning transmission electron microscopy
A time series of high-resolution scanning transmission electron microscope (STEM) images were obtained with a spherical aberration-corrected electron microscope equipped with a cold-type field emission gun, JEOL JEM-ARM 200 F, operating at 200 kV. The convergent semi-angle of the incident beam was set to 20 mrad. The inner and outer semi-angles of the annular detector were set to 65 mrad and 265 mrad, respectively. A series of the STEM images with 1024 × 1024 pixels were acquired with a pixel dwell time of 1 µs. Although the specimen did not show irradiation damage in the STEM observation time series, the electron beam might increase the temperature in the illuminated areas, inducing the thermal fluctuation observed in the study.
The fast Fourier transform (FFT) patterns and Fourier-masked images were obtained by computer software, Gatan Digital Micrograph. The diameters of the masks used to select the superlattice reflections were set to 1 nm −1 .

Computational details
We employed the WIEN2k code 32 based on the full-potential linearized augmented plane-wave method for our density functional theory (DFT) calculations. The calculated results obtained in the generalized gradient approximation (GGA) for electron correlations were presented with the exchange-correlation potential of ref. 33 . To improve the description of electron correlations in V 3d orbitals, we used the rotationally invariant version of the GGA+U method with the double-counting correction in the fully localized limit 34,35 . The main text presents the results obtained at U = 4 eV. In the self-consistent calculations, we used 16 × 16 × 7, 10 × 10 × 18, and 13 × 13 × 11 k-points in the irreducible wedges of the Brillouin zones for the P3m1, Pm, and P31m phases, respectively. Muffin-tin radii (R MT ) of 2.15 (Li), 2.35 (V), and 2.03 (S) Bohr were used assuming a plane-wave cutoff of K max = 7.00/R MT . We used VESTA 29 and XCrySDen 36 for graphical purposes.