Configurational Electronic States in Layered Metallic Dichalcogenides

Mesoscopic irregularly ordered and even amorphous self-assembled electronic structures were recently reported in two-dimensional metallic dichalcogenides (TMDs), created and manipulated with short light pulses or by charge injection. Apart from promising new all-electronic memory devices, such states are of great fundamental importance, since such aperiodic states cannot be described in terms of conventional charge-density-wave (CDW) physics. In this paper we address the problem of metastable mesoscopic configurational charge ordering in TMDs with a sparsely filled charged lattice gas model in which electrons are subject only to screened Coulomb repulsion. The model correctly predicts commensurate CDW states corresponding to different TMDs at magic filling fractions $f_m=1/3,1/4,1/9,1/13,1/16$. Doping away from $f_m$ results either in multiple near-degenerate configurational states, or an amorphous state at the correct density observed by scanning tunnelling microscopy. Quantum fluctuations between degenerate states predict a quantum charge liquid at low temperatures, revealing a new generalized viewpoint on both regular, irregular and amorphous charge ordering in transition metal dichalcogenides.


Introduction
whether the state is an ICCDW or a CCDW in 2D dichalcogenides 10 . More recently, arguments have also been presented for a significant role of the EPC in 1T-TaS2 42,43 . Thus, in spite of the fact that the problem is known since the 1970s, there are many open issues that have not been resolved even in 1T-TaS2, let alone universally in the TMD materials.
Recent scanning tunneling microscopy (STM) experiments have shown that doping introduces densely packed domains and sharp domain walls (DWs) that are not easily understood using the established CDW physics discussed above. Such textures also seem to appear under various external perturbations such as photodoping 5,38,44 , charge injection [2][3][4][5]45 , non-isovalent 46 and isovalent 47,48 transition metal substitution, chemical doping by chalcogen ion substitution 49 or even intercalation of transition metal ions 50 . Apart from DW textures, a few examples of amorphous states have been reported as a consequence of photodoping or charge injection 51 as well as alkali or N2H4 vapor deposition 52,53 and possibly Nb intercalation 47 .
The states created through non-equilibrium methods are of particular interest for ultrafast memory devices 5,54-61 , but the mesoscopic metastability associated with these phenomena are not easily explained with conventional CDW theory.
Superconducting 2H-NbSe2 exhibits a domain wall structure even in its pristine state with no lock-in transition to a CCDW as in other 2H compounds 62 . Superconductivity thus seemingly coincides with the appearance of a textured state not only in 2H-NbSe2, but also in 1T-TiSe2 15 (and 1T-TaS2 11 , another generic feature which conventional theories do not address.
Thus, in spite of intense current interest, it is clear that none of the above-mentioned models which discuss k-space interactions in a homogeneous periodic lattice can describe the intricate irregular topologically non-trivial mesoscopic features that universally appear in doped and metastable photoinduced or charge-injected structures [2][3][4]38,44 . Such discommensurations were discussed theoretically to some extent by McMillan 23 based on Landau theory in 2H-TaSe2, and Nakanishi and Shiba 63,64 for 1T-TaS2, but only as regular structures.
In this paper we approach the problem of describing irregular charge configurational states from the viewpoint of strong correlations between electrons involved in CDW formation in the spirit of the Mott state introduced by Fazekas and Tossatti 37,65,66 .
Our approach has its origins in CLG modelling of strongly correlated polarons on a square lattice as applied originally to oxides [67][68][69][70] and extended by Brazovskii 71 to investigate the CCDW state of 1T-TaS2. DW formation resulting from doping was suggested by Karpov and Brazovskii 45,71 . Here we examine the phase diagram of a sparse charged lattice gas (CLG) of electrons on a triangular lattice at different filling (i.e. doping) levels using extensive Monte Carlo CLG simulations.
The experimental evidence for local lattice distortions associated with the CDW implies that the electrons are effectively coupled to the lattice, which is usually discussed in terms of polarons, justified in part by the large difference between static and high-frequency dielectric constant that implies substantial screening (as discussed later).
The charged lattice gas (CLG) interpretation of a physical system is typically a very crude oversimplification, yet its intuitive nature is very appealing. In our case it consists of a system of charged repulsively interacting electrons whose motion is slowed down by interaction with the lattice.
Despite its simplicity, it has been utilized successfully as a model for a variety of seemingly unrelated physical systems. Superconducting vortices either in Josephson junction arrays 72,73 or superconducting layers in the presence of periodic pinning [74][75][76][77][78] have been studied extensively by Monte Carlo simulations and analytical approaches [79][80][81][82][83][84][85][86][87][88][89] as well as molecular dynamics simulations 90,91 . Adsorbed atoms on a crystalline surface [92][93][94][95][96][97] , colloidal systems on templates 98 , organic conductors 99,100 , as well as magnetic bubble arrays 101 also exhibit CLG characteristics. Recently, the configuration of a system of polarons in 1T-TaS2 at very low temperatures has been investigated using the CLG model. Remarkably, this model describes well not only the equilibrium state 71 but also two different photoinduced states 45,51 . This success may be considered as justification for attempting to generalize the CLG model for understanding the phase diagram and properties of polaron gases.
Remarkably, we find that the low temperature CCDW states in all the aforementioned systems correspond to electronic crystals, where the filling of the system is one of the many possible magic fillings . We find that there exists an infinite number of polaron arrangements at non-magic fillings, where the energy gap between the analytically predicted ground state and the first excited long range ordered state can be made arbitrarily small in the thermodynamic limit. Therefore, there exists a whole spectrum of energetically low lying states and this is where an amorphous state is very close in energy to the analytically predicted ground state. Monte Carlo simulations extrapolated into the thermodynamic limit show for a finite amount of cases that this is indeed the case. Moreover, we prove that a triangular lattice cannot exist at ≠ due to the frustration. We conclude that a translationary invariant state with hexagonal symmetry at ≠ indeed needs to exhibit a specific domain wall structure. This is consistent with previously developed theories concerning domain walls 102,103 .
From these findings we can draw parallels with a quantum spin liquid. On a triangular lattice a spin liquid is typically understood as a superposition of an infinite number of valence bond states, which all break hexagonal symmetry and have the same energy. The spin liquid ground state then entangles all the states preserving hexagonal symmetry while lowering the energy at the same time 104 . We find a similar situation with the polaron charge configurations in our calculations in the sense that the low lying energy configurations entangle into a quantum charge liquid. The MC simulations also confirm a clear absence of a first order phase transition in a wide temperature range in all cases where ≠ , which supports the idea that there are two quantum critical points present at = 0 around each magic filling. This is shown on the phase diagram of the CLG on a triangular lattice presented in this paper.
Our conjecture also has important consequences regarding the current paradigm of glass formation 105 .
The paradigm states that upon cooling a system down to low temperature, a glass is formed due to the system getting stuck in a local minimum of free energy. There is an overwhelming probability in favor of this phenomenon due to an extensive number of local minima present in the energy landscape of the system. One of the main assumptions in the glass formation paradigm is that the global minimum of the system is the crystalline state, which is challenged by our calculations.

The Model
The model we employed in this paper is based on the CLG Hamiltonian

Analytical Considerations
From the structure of the Hamiltonian in Eq. (1) it is obvious that the polarons will always tend to arrange themselves in a triangular lattice with the lattice constant = /� , as the most closely packed structure in 2D. However, we need to take into account that polarons can only reside on a triangular underlying atomic lattice (AL). The problem is trivial in the case where the polaron lattice (PL) coincides with the AL. In this case, every point of the PL can be mapped onto a point on the AL.
This condition will be dubbed as the hexagonality condition and considered in the following paragraph.
Imagine a 2D plane on which the AL resides. Let the origin of the coordinate system used to describe the points of the plane lie on some AL point and let an arbitrary point on the plane be parametrized by polar coordinates ( , ). Also, let an arbitrary AL point be parametrized by polar coordinates ( , ), where runs through the infinite set of AL points. In order to find all the possible triangular PLs we have to find all the possible triangular lattice unit cells, which map onto the AL. Therefore, we need to find two primitive lattice vectors of a triangular lattice unit cell such that both correspond to some two points on the AL. The hexagonal symmetry of the problem helps because we only need to find one primitive vector. The second one is immediately determined as a /3 rotation of the first and will The fillings which satisfy this condition are dubbed as the magic fillings. It is easy to see that 1/ = 13 corresponds to = 3 and = 1 and 1/ = 7 corresponds to = 2 and = 1 for example.
One interesting fact can be deduced from the mirror symmetry of the problem with respect to the =

Monte Carlo Simulations
We performed MC simulations (for details on the parallel tempering method we employed see SI) with periodic boundary conditions and analytical considerations in order to find the ground state of the system at fixed fillings ranging from 1/2 to 1/21 with the increment of 1 in the denominator. The screening radius ranged from to 100 but here we focus on just one value of  There is also a significant amount of intensity present in between the peaks.
We also performed MC simulations for different values of . The overall trend is that the system becomes increasingly harder to equilibrate at larger values of . These calculations also reveal the presence of phase coexistence and ordering of domain walls, as well as higher energy amorphous states (see SI for examples).
In Fig. 4 we plot the phase diagram of the CLG model on the basis of MC simulations together with the effective (normalized) gap Δ/ . Clearly, the transition temperature increases significantly at magic fillings. This is also where the phase transition is of the first order, and not a glassy transition as it is for ≠ . By glassy transition we mean that MC simulations indicate the presence of a weak continuous transition at a certain temperature as is depicted in Supplementary Figures 8 through 11, even though it is obvious from the low temperature configuration snapshots that no ordering is present (see SI). We interpret this transition as the onset of glassy dynamics, where correlations slow down the movement of polarons and eventually freeze the system. Therefore, the transition temperature may be regarded as the interaction scale at which the entropy from all the different low energy long range ordered polaron configurations starts to dominate.
We can speculate that in the quantum case of spinless fermions there are two quantum critical points in the vicinity of each magic filling. Karpov et al. 45 have shown that in the case of 1/ = 13 the crystal lattice is still present in the system at small negative or positive dopings with the addition of domain walls. If we define the order parameter of the system as the FT peaks corresponding to the 1/13 lattice, then it is obvious that the order parameter does not disappear immediately under doping.
It is very natural to assume that this is the case for all magic fillings. However, at the same time our MC simulations clearly show that if we dope the system enough, the first order phase transition becomes a glassy transition (see SI). This implies that the order parameter disappears in this case.
Likely, there is a quantum critical point somewhere in between. This is the case for positive and negative doping but the two critical points are not necessarily symmetric around each . Also, it is possible that at low enough fillings the two critical points merge and even overlap and therefore disappear, as the magic fillings get closer and closer to one another.
In order to avoid the scaling of the gap Δ with we overlayed the plot of the relative gap Δ/ on top of the phase diagram. The peaks in the plot agree well with the peaks in the phase diagram, which explains the emergence of a first order transition at . There is a well pronounced free energy minimum at , which is not true in the frustrated case. The cases 1/4, 1/9 and 1/16 stand out due to the fact that a square lattice is energetically very close to the triangular lattice. This is probably true for

Exact Diagonalization of a Charged Lattice Gas With Quantum Tunneling
In order to substantiate our claim regarding the quantum charge liquid, we carried out exact diagonalization calculations using the Lanczos method in order to find the ground state. The code used in the calculations was from the ALPS initiative, version v2.3b1 107-110 . We were interested in how does the introduction of quantum effects such as the hopping of electrons alter the charge order compared to the classical case. The model we used is the model of spinless fermions with a fixed filling and both nearest and next-nearest neighbor interaction where ( † ) is the annihilation (creation) fermion operator on the -th triangular lattice site, = † , is the hopping parameter, 1 the interaction between nearest neighbors, 2 between nextnearest neighbors and is the number of polarons in the system. We set the system size to 5 × 5 and = 5 in order to explore the frustrated case of = 1/5 and also set 1 = 10, 2 = 5 and varied .
Supplementary Figure 13 shows the density correlation function 〈 〉 for the classical case = 0 and two quantum cases = 0.5 and = 10.
The correlation function clearly shows that for small values of , the interaction driven charge ordering dominates in the system due to the fact that the correlation functions in the = 0 and = 0.5 cases are practically the same. This is not the case for larger values of , where a metallic phase is present as one would expect for a system where kinetic energy of electrons dominates. These calculations by themselves do not prove the existence of a QCL because system size limitations limit the accuracy.
Indeed, the density correlation function of the small 5 × 5 system is substantially different from that from classical Monte Carlo simulations. However, they do confirm that the presence of a large number of nearly degenerate configurations in the frustrated cases most likely plays a significant role.

Experimental Realizations of Predicted Charge Orders
We have collected all of the experimentally studied systems for which our model is applicable in  Table 1 A collection of all the experimentally studied systems for which our model is applicable.
Every system in left column for which we were able to extract the filling within the scope of our model is associated with its and corresponding theoretical phase in the middle and right column, respectively. In the six examples in which this was not possible, we specify the corresponding theoretical phase and only the closest in the case of a domain state and the approximate in the amorphous case with a ∼.
The first case (Fig. 5a,e,i) of photoexcited 1T-TaS2 shows the pristine state ( = = 1/13), the lightly photoexcited "hidden" state 116 and moderately photoexcited "amorphous" state 51 . With light photoexcitation, the domain state is reached, while for moderate photoexcitation, an amorphous state is reached. The measured values of are 1/13, 1/12.6 and 1/11 respectively. For ≠ , the density is determined by counting the polarons in the respective images using standard procedures described in 51 . In the amorphous state, falls in between two magic numbers ( = 1/9 and 1/12).
The final example that we highlight here is doping by alkali or N2H4 vapor deposition on the surface of These examples were chosen because they display diverse doping mechanisms, but show generic behavior with doping as predicted by the phase diagram in Fig. 4.
There are also two other compounds worth mentioning, which could not be presented either in Fig. 5 or in the list of polaronic crystals. The case of intercalated 1T-NbxTaS2 47 exhibits a 1/13 polaronic crystal at = 0, a domain state at = 0.04 and = 0.07 and a possible amorphous state at = 0.1.
However, we could not reliably determine the system's filling and therefore we excluded this case from the analysis in Fig. 5. It should also be noted, that the role of structural disorder is not clear in this case. The second case is 2H-NbSe2, which even in its pristine form does not exhibit a perfectly commensurate CDW. However, the work of Chatterjee et al. 62  often not well understood. Resistivity relaxation measurements suggest that quantum configurational reordering processes may play a role in transport 41 .
The microscopic effect of doping in all the different compounds is not yet well understood and is not universal. Substitution seems to cause the formation of more polarons in the system. Intercalation also seems to cause the formation of more polarons, except for the 1T-Cu0.08TiSe2 case. The role of fluence in photodoping has also not yet been fully explored, however the work of 1,51 suggest that larger fluences lead to more polarons being formed. We also note that the isovalent Se for S substitution in 1T-TaS2-xSex effectively introduces doping though strain induced by the size mismatch, but the end members (1T-TaS2 and 1T-TaSe2) both have = 1/13 17 .
The apparent outlier of all the compounds is 1T-TiSe2, where exciton condensation is the proposed mechanism for CCDW formation 25 . However, this approach on its own cannot explain the presence of the observed DW structure which emerges under pressure 15 or doping with transition metal ions 48 , nor does it take into account Coulomb interactions. We propose that 1T-TiSe2 might be a good example of cooperation between strong polaron correlations enhanced by exciton crystallisation that favours the CCDW with = 1/4.
The model presented in this paper is consistent with the observations in the work of Dai et al. 111 on Fe doped 2H-TaS2, 2H-TaSe2 and 2H-NbSe2. The change of the superlattice type in experiments can be explained as the change of the magic filling of the present polaronic crystal. They also report on distortions in the superlattice, which is consistent with the domain wall state as well as the amorphous state in between magic fillings.
One interesting fact is the absence of the 1/7 superlattice predicted by the model, which has not been observed anywhere else. Doping also seems to increase the importance of the role of correlations in the system due to the onset of commensurability from incommensurability in 2H-TaSe2 and 2H-NbSe2. The induced frustration in all the non-magic filling cases causes the system to break hexagonal symmetry and an amorphous state becomes very close in energy to the analytical ground state. The resulting hyperuniform jammed packing is the densest possible maximally correlated state in 2D that is not periodic. We interpret the fact that an amorphous state is so close in energy to the analytical ground state to be a consequence of the presence of a large amount of entropy even at very low temperatures. There seem to be a large number of states available to the system at low temperatures and they are mixed into an amorphous state even if the temperature is only slightly above zero. We conjecture that in every case where the hexagonality condition (see section Analytical Considerations) is not satisfied, the gap between the ground state and the first excited state can be made arbitrarily small.
From the fact that there are many broken-symmetry states present at low energies, we can speculate that for the quantum case of the CLG model (spinless fermions), where one introduces even a very small amount of nearest neighbor hopping ( ≪ 0 ) and fixes the system's filling, the ground state will entangle the low lying energy states (| 〉) into a quantum charge liquid (| 〉 = ∑ | 〉, where 〈 | 〉 = 1) analogue of a quantum spin liquid with no long range order. The system would therefore have an amorphous ground state which breaks the current paradigm of glass formation by which a crystalline state is the eventual ground state of a glassy system. It is important to note that needs to be small enough so that correlations are still dominant in the system but large enough to retain entanglement between the low-lying states. At larger values of ≠ a so called "pinball liquid" state was proposed to be the ground state 100,117 . However, due to limitations in the system size of quantum calculations, we believe that true long range order is still elusive in the two abovementioned studies. A much larger system size would be necessary in order to achieve conclusive results. An amorphous ground state has already been proposed within the context of frustration between two different length scales of the interaction 118,119 . In our case frustration emerges from the competition between the interaction-driven polaron lattice and the underlying atomic lattice on which the polarons reside (see Analytical Considerations). The QCL in both the domain state or the glassy regime would in the theoretically ideal case exhibit a homogeneous real space structure, assuming quantum coherence across the whole system. However, in actual experiments impurities would force the system to choose a specific configuration as observed in 38,51 . Therefore, only states local to the energy of the chosen state would be susceptible to entanglement. An interesting implication which follows is that any kind of observed polaron tunneling or even tunneling between different configurational states could be explained by the QCL concept.
Finally, let us estimate the macroscopic tunneling rate between the near-degenerate configurational states. The barrier tunnelling rate depends on the barrier height as ∼ (− √ ), where is a constant. The energy of the system is proportional to the number of particles involved , and therefore in the thermodynamic limit as → ∞, the tunneling rate → 0. However, there may be processes in which a finite number of particles are involved, so is finite. If the configurations do not differ significantly, the tunneling rate may be significant, and observable.
Reconfiguration processes should show a cross-over from tunneling to thermally-activated classical behavior above some threshold temperature. Experimentally this should be observable as a cross-over from temperature-independent tunneling to thermally activated behavior at some characteristic temperature. Such processes may be observed when a small number of polarons are involved in a slight reconfiguration of domains, for example.
In the amorphous state jamming may inhibit the relaxation process 51 , and the system is exceptionally stable, in agreement with observations.

Conclusions
Remarkably, the generic physical picture which emerges appears to predict the CCDWs wavelengths at:   111 , that reveal a transition from one magic filling to another under substantial doping show that it is possible to move between magic fractions along the doping axis within a single compound.
Finally, given that superconductivity appears in between and coincides with the appearance of domain textures, this work suggest that new superconductors may be searched for in other areas of the phase diagram that were not previously investigated for superconductivity. This may also reveal some valuable insight in answering the question how the appearance of textures is related to superconductivity -a question that has wider and persistent significance in cuprates and other textured high-temperature superconductors.    1T-Ti0.07Ta0.93Se2 ( ≈ 1/12.6) and 1T-TaSeS ( ≈ 12.6) respectively. The last case i represents a glassy state in photodoped 1T-TaS2 (fluence > 3.5 / 2 , ≈ 1/11).