*Article* **A Database for Crystalline Organic Conductors and Superconductors**

**Owen Ganter 1,\*,†, Kevin Feeny 1, Morgan Brooke-deBock 1, Stephen M. Winter <sup>2</sup> and Charles C. Agosta 1,\***


**Abstract:** We present a prototype database for quasi two-dimensional crystalline organic conductors and superconductors based on molecules related to bis(ethylenedithio)tetrathiafulvalene (BEDT-TTF, ET). The database includes crystal structures, calculated electronic structures, and experimentally measured properties such as the superconducting transition temperature and critical magnetic fields. We obtained crystal structures from the Cambridge Structural Database and created a crystal structure analysis algorithm to identify cation molecules and execute tight binding electronic structure calculations. We used manual data entry to encode experimentally measured properties reported in publications. Crystalline organic conductors and superconductors exhibit a wide variety of electronic ground states, particularly those with correlations. We hope that this database will ultimately lead to a better understanding of the fundamental mechanisms of such states.

**Keywords:** organic superconductor; superconductivity; charge density wave; spin density wave; spin liquid; FFLO state; materials database; data science

#### **1. Introduction**

Machine searchable databases that contain structural properties of related materials, calculated electronic structure, and measured electromagnetic properties, are providing a new way to design advanced functional materials. In addition, consolidating structural and functional information may lead to a better understanding of the microscopic mechanisms of correlated electron materials. Herein, we detail the launch of a new database of crystalline organic materials, many of which are conducting or superconducting, with the goal of motivating data-centered research to enhance the understanding of lower dimensional correlated electron materials. The database can be accessed through a website at osd.clarku.edu.

The crystalline organic materials (COM) are well suited to create this type of database first of all because they are interesting experimentally. Partially driven by their low dimensionality, this class of materials exhibits a variety of competing electronic behaviors [1–5] including metallic conductivity, Mott insulators [6–8] , antiferromagnetic states [9–11], and superconductivity [12,13]. Other forms of long range charge order have also been observed, such as charge density waves (CDWs) [14–16] and spin density waves (SDWs) [17,18]. More exotic long range order, such as the quantum hall effect [18] and some of the first believable evidence for field induced inhomogeneous superconductivity (the FFLO state) [19–23] were also found in COM. There has also been discussion about the existence of tilted Dirac points [24] and spin liquids [25–27] in these organic salts.

In addition to their rich correlated electron behavior, COMs are easy to access theoretically because they form regular stoichiometric crystals based on a few common cation

**Citation:** Ganter, O; Feeny, K; Brooke-deBock, M; Winter, S.M.; Agosta, C.C. A Database for Crystalline Organic Conductors and Superconductors. *Crystals* **2022**, *12*, 919. https://doi.org/10.3390/ cryst12070919

Academic Editor: Andrej Pustogow

Received: 10 May 2022 Accepted: 21 June 2022 Published: 28 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Department of Physics, Clark University, 950 Main Street, Worcester, MA 01610, USA; kfeeny@clarku.edu (K.F.); mbrookedebock@clarku.edu (M.B.-d.)

<sup>2</sup> Department of Physics and Center for Functional Materials, Wake Forest University, 1834 Wake Forest Road, Winston-Salem, NC 27109, USA; winters@wfu.edu

molecules held together with various anion complexes. It is the electron deficient cation layer that contains holes, which enables itinerant electron behavior within the layer. It is evident that the geometric arrangement of the cation molecules is a principal factor that determines the electronic ground state of the system. This so-called packing of the cation molecules can be altered either by placing different molecules into the anion layer, or by applying external pressure. By adjusting the physical parameters of the cation layer, competing correlated electron states are selectively enabled and a complex electronic phase diagram can be traversed [28–30].

There has been much theoretical progress towards understanding how the various degrees of freedom of these materials lead to the bulk electronic states that are observed [31–33]. From a theoretical perspective, the materials offer a unique window into the physics of correlated electrons, because (i) COMs span the full range of states from interaction-induced insulators to superconductors and metals, (ii) the electronic structures are relatively simple, with only one orbital per molecule typically being relevant, (iii) nearly all the materials are stoichiometric and exhibit a high degree of crystalline order, so that simple models may closely approximate experiments, and (iv) COMs are built from common molecular entities, so that variations in properties may be directly related to structural variations across vast numbers of compounds.

The theoretical study of these materials via the tight binding model, or density functional theory (DFT), is leading to a better understanding of the fundamental physics behind correlated electron systems and of quantum materials in general. Various packing sub-families are known, in which particular structural degrees of freedom, and their relationship with the underlying hopping integrals, are key inputs to predict the behavior of the system [34–37]. Highlights of the strong exchange between experiment and theory include, e.g., quantitative agreement between results of high level calculations (dynamical mean field theory) and measurements of Mott critical scaling [7] and correlation-driven crossovers in the optical response [38].

This new database was built to be used both as a way to find representative materials for targeted experiments investigating particular correlated electron states, or the proximity of competing states, and as a research tool for discovering structure–function relationships [39,40]. An aspirational use of the database would be to predict and design new materials with targeted electronic properties. To serve these functions, the database contains experimentally measured properties, crystal structures, and calculated electronic structures of quasi one and two-dimensional crystalline organic conductors and superconductors. We hope that the enhanced accessibility of information that our database provides will serve the scientific community and lead to new discoveries.

#### **2. Website and Database**

When arriving at the home page of the website, the user is presented with a list of available materials and a window in which one may specify a desired type of packing, cation molecule, and anion molecule in order to search for any matches in the database. Leaving a field blank will act as a wildcard. The user can then click on a material to navigate to the home page for that material, which shows its crystal parameters followed by an interactive set of graphics starting with rotatable views of the crystal structure, the calculated electronic structure, and the cation morphology. These views are followed by available measurements of the material properties. Individual pages for each measurement enable the user to view associated information in greater detail. When observing crystal structures, they can be viewed and filtered to include all of the atoms, only the cations, only the anions, or rectangles to represent the cations. Electronic structure diagrams can also be customized by the user to show different k-paths or results from various types of calculations.

In order to identify relevant information to include in the database, the websites of journals were automatically searched for key phrases such as "BEDT-TTF", "organic conductor", and "organic superconductor". The papers resulting from those searches were then recorded as potentially containing measurement information relevant to the database. Each paper of interest was parsed to determine its relevance, and encode any reported measurements into the database. Candidates for crystal structures were obtained from the Cambridge Structural Database (CSD) by performing a substructure search on the set of known cations. An algorithm was then used to analyze each crystal structure and determine if it was a lower dimensional charge transfer salt of interest. Relevant crystal structures were then added to the database, and their electronic structures were computed automatically. For a smaller selection of materials DFT (WIEN2k) was also used to calculate the electronic structure to compare to the tight binding results. The organization is depicted in Figure 1. The database currently contains 110 materials, 184 crystal structures obtained from the CSD and 440 measured properties. A link is provided for each material to the corresponding CIF entry on the CSD website for the crystal structure information, and a link is provided to the paper where each measurement was found. Below we describe in more detail how the electronic calculations were made, and the methods for encoding the measurements into the database.

**Figure 1.** Diagram showing the structure and flow of information into the database.

#### **3. Crystal Structure Analysis Algorithm**

To automatically identify relevant crystal structures from the Cambridge Structural Database, we created an algorithm to assess relevance, and to perform some preliminary diagnostics. The algorithm reads a CIF file, which contains a list of the atomic coordinates within the unit cell, and the lattice information. The distance between each pair of atoms is calculated, and if that distance is below the bonding threshold distance for the given atomic species, it is assumed that a bond exists between them. In this manner, the molecules within the unit cell are identified (see Figure 2a,b). The structure of each molecule is then compared to a list of predefined structures of cation molecules such as BEDT-TTF. If the structures match, it is then known that the molecule is a cation molecule of interest. At this point it can be determined whether the material is quasi two-dimensional or not, and if so, what crystal axis is perpendicular to the layers. This is achieved by examining the overlap between cation and non-cation molecules for each Cartesian axis. If along a certain axis there are only overlapping cation molecules with no non-cation molecules, we make the assumption that this axis is in the plane of a conducting layer.

With the orientations and identities of the cations and anions automatically identified, this algorithm can automatically suggest the chemical formula and packing type of a given crystal structure. This is useful in grouping multiple structures for the same compound; however, we did not use this detected chemical formula as the compound label in the database. The name of a crystal structure is entered as denoted in the original paper to ensure that any special chemical naming conventions used by the authors are preserved in the database.

In order to display the packing geometry of the cation molecules to the user, the best fit plane of each cation is computed, and a minimum area rectangle algorithm is applied to generalize a cation molecule as a rectangle in three dimensional space. This type of generalization is common in cartoon diagrams of these materials. The resulting unit cell of cation-rectangles, shown in Figure 2c–f, is useful to inspect the packing geometry of the material. Given that the overlap of the cations is the major determinant of the electronic structure of the crystalline organics, seeing the cations as blocks to visually show the morphology of the crystal symmetry of the cations is instructive. The relative distances and angles of the rectangles can then be calculated to quantitatively analyze the packing. We also generated two-dimensional diagrams, which more simply showed the geometric orientation of the cation molecules within the conducting layers based on the angles calculated before.

**Figure 2.** Views of various stages of the crystal structure analysis algorithm. The individual images show the initial unit cell consisting of atoms within a lattice (**a**), detected molecules (**b**), detected cation molecules shown as rectangles (**c**), detected anisotropy and layers (**d**), a new unit cell consisting of only the cation layer (**e**), and a two-dimensional depiction of the cation layer (**f**). In this example, the material is *β*"-(ET)4[(H3O)Cr(C2O4)3]2[(H3O)2]5H2O, as reported in [41].

#### **4. High-Throughput Electronic Structure Calculations**

High-throughput electronic structure calculations using density functional theory (DFT) have found applications in several materials databases; however, the large unit cells of these materials and low crystal symmetries make full-scale DFT calculations with planewave basis sets computationally expensive. This is especially problematic for cases with open-shell anions, which feature localized unpaired electrons. Unless the local correlations in the anion layer are treated explicitly (via DFT + U), anion bands may appear erroneously near the Fermi energy, yielding incorrect Fermi surfaces. More importantly however, many of our crystal structures have missing or disordered atoms, especially in the anion layer. This is particularly prevalent in anion layers that contain solvents, which are often disordered across unit cells. It is therefore necessary to reduce the computational expense and focus exclusively on the cation layers.

To calculate the electronic structure of every crystal structure in our database, we construct a two-dimensional tight binding model [42] for highest occupied molecular orbitals in the layer of cation molecules using a series of local DFT calculations. Solving the tight-binding (TB) model produces the band structure. To carry out the calculation, we used the crystal structure analysis algorithm previously described to identify all symmetrically equivalent molecules and pairs of molecules with the cation layers. We then used the method employed in [43] to estimate tight-binding hopping integrals using quantum chemistry packages (in this case ORCA [44]). Results are shown in Figure 3, for the example of *α*-(ET)2KHg(SCN)4. The method is based on calculations on pairs of molecules in which the local crystal environment is otherwise ignored, which significantly reduces computational expense. For this purpose, we used basis sets including 3-21G, 6-31G, 6-311G, and def2-SVP in conjunction with the B3LYP hybrid density functional. Localized Wannier molecular orbitals (MOs) are constructed for each molecule via maximizing the overlap with the corresponding orbital of the isolated molecules. The procedure is as follows:

1. **Obtain Isolated MOs:** For each molecular pair (labeled *i*, *j*), a calculation is first performed on the isolated molecules. From this, the MO coefficients (in the basis of Gaussian atomic orbitals) for each molecule are obtained as Φ<sup>0</sup> *<sup>i</sup>* and <sup>Φ</sup><sup>0</sup> *<sup>j</sup>* . These are combined as:

$$
\Phi\_0 = \begin{pmatrix}
\Phi\_i^0 & 0 \\
0 & \Phi\_j
\end{pmatrix} \tag{1}
$$

2. **Construct Wannier Functions:** For each molecular pair, a calculation is then performed in the geometry corresponding to the crystal structure. This produces the diagonal MO energies **E**, the overlap matrix **S**, and the MO coefficients Φ. In ORCA, **S** is output in the atomic orbital basis. It is first rotated into the basis of the isolated MOs:

$$
\dot{\mathbf{S}} = \Phi\_0 \mathbf{S} \,\Phi\_0^\dagger \tag{2}
$$

In this geometry, the basis of isolated MOs are no longer orthonormal. Thus, the local Wannier functions are constructed via symmetric orthornormalization, Φ¯ <sup>0</sup> = **S**¯ <sup>−</sup>1/2Φ0.

3. **Rotate Fock Matrix:** The diagonal orbital energies are then rotated into the abovedefined localized MOs:

$$\mathbf{F} = \Phi\_0 \,\,\Phi^{-1} \,\,\mathbf{E} \, (\Phi^\dagger)^{-1} \,\, \Phi\_0^\dagger \,\tag{3}$$

The resulting Fock matrix has the structure:

$$\mathbf{F} = \begin{pmatrix} \mathbf{F}\_{ii} & \mathbf{F}\_{ij} \\ \mathbf{F}\_{ji} & \mathbf{F}\_{jj} \end{pmatrix} \tag{4}$$

The on-site terms **F***ii* and **F***jj* now contain both the diagonal Wannier orbital energies, and small off-diagonal "crystal field" contributions. It is advantageous to remove the latter terms via unitary transformation:

$$\mathbf{F} = \begin{pmatrix} \mathbf{U}\_i & 0 \\ 0 & \mathbf{U}\_j \end{pmatrix} \mathbf{F} \left( \begin{array}{cc} \mathbf{U}\_i^\dagger & 0 \\ 0 & \mathbf{U}\_j^\dagger \end{array} \right) \tag{5}$$

where **F**¯*ii* = **U***<sup>i</sup>* **F***ii* **U**† *<sup>i</sup>* and **<sup>F</sup>**¯*jj* <sup>=</sup> **<sup>U</sup>***<sup>j</sup>* **<sup>F</sup>***jj* **<sup>U</sup>**† *<sup>j</sup>* are diagonal. The intersite hoppings can then be read from **F**¯*ij* = **U***<sup>i</sup>* **F***ij* **U**† *j* .

We note, because this latter unitary transformation is different for every molecular pair, the hopping integrals obtained for different pairs represent slightly different definitions of the local Wannier functions. Nonetheless, this approximation is no more severe than the pairwise construction inherent to the method. Although this approach neglects the anion layer, the results agree well with full-scale calculations performed with Wien2k (at the GGA level) and experimental electronic structure as well (see Figure 3c); however, the former approach is much faster. With Wien2k, for example, a full DFT calculation using GGA functionals can take several days to complete (with 100 processors), while construction of a TB model with ORCA calculations takes approximately five minutes per compound (with 10 processors), even when using more expensive hybrid functionals. Such a speed-up is desirable when making high-throughput calculations for each crystal structure entry in the database.

Because the pairwise calculations are made separately, we had to adjust the signs of the resulting charge transfer integrals such that the phase of the molecular orbital on each symmetrically equivalent molecule was the same. We used the centroids of the cation molecules for the positions of the sites, disregarding the out of layer component. The filling of each model was deduced from the stoichiometry and charge of the cation molecules. Solving the tight binding Hamiltonian at each point in k-space produced the energy eigenvalues that constitute the electronic structure. In this manner, the band structure, Fermi surface, and density of states are automatically computed and may be viewed on the website. Users can interact with these data directly by selecting the k-path to use and which basis sets to display. In addition, the computed hopping integrals, and their locations in the unit cell are presented to the user in order to serve as a basis for further theoretical modeling.

**Figure 3.** Views of various stages of the tight-binding electronic structure calculation. The individual images show the quantum chemistry calculation of inter-molecular charge transfer integrals (**a**), site-hopping integrals model (**b**), and electronic structure (**c**). In this example the material is *α*-(ET)2KHg(SCN)4, as reported in [39]. Note that the WIEN2k band structure in red is close to the TB band structures, especially near the Fermi level.

#### **5. Measurements**

	- **–** Metal insulator.
	- **–** Superconductivity, *Tc*, *Hc*1, *Hc*2, *HP*, (where *HP* is the Pauli paramagnetic limit).
	- **–** Charge density wave.
	- **–** Spin density wave.
	- **–** Magnetic ordering.

As a function of:


Each measurement entry in the database consists of three blocks of information: a block specifying the state of the system being measured, a block specifying the value of the measurement, including the error bars if available, and a block specifying the method by which the measurement was made. Our goal behind the data entry is to create a digital copy of measurement information in the precise manner in which it was specified by the authors who made the measurement. We implemented support for as many measurements as we could for this process. For example, the ability to specify a numeric value using an exact decimal, a range between two decimals, or an average decimal with a plus or minus value. Although laborious, we found that manual data entry was the most reliable way to extract measurement information from papers. We used a web application on the website for this purpose, shown in Figure 4.

Once we had a sufficient number of measurements, there were many decisions that needed to be made about how they were presented. For a deep understanding of a single material, it is necessary to see a detailed view of the actual measurements labeled by the method of measurement; for example, the superconducting transition temperature, *Tc*, found by resistance or specific heat, and the point on the transition curve, e.g., onset or midpoint, used to locate the transition, with citations for each measurement. For that reason, the details of the measurement method and the error bars, if given, are stored in the database. For that reason, it is also necessary to present an average value for a material when a number of materials are being compared to each other. We made the decision to discount some of the grossly outlying measurements in cases where we thought the data were not convincing. These rules for curating the data are constantly being reconsidered to present the most useful data to the community; however, the full collection of measurements will always be available so that a user of the database can analyze the published measurements with their own algorithms.


**Figure 4.** View of the web application used for data entry. In this example, the *Tc* of *β*-(TMTSF)2PF6 under 12 kbar of applied pressure is specified to be 0.9 Kelvin as reported in [12].

#### **6. Discussion and Future Outlook**

We present this database of crystalline organic conductors and superconductors as an evolving tool that will be continuously updated with new materials, features, and metrics. The goal of gathering all relevant calculation and measurement information into one central location is an arduous one, but the progress that we have made so far illustrates that it is possible. In further developing this database, we have two main goals.

Our first main goal is to populate the database with as many entries as possible. The limiting step in this process for inclusion of experimental measurements is data entry. Manual data entry of measurements from scientific articles is the only method with a

high enough degree of reliability to be useful in our database. We have used members of our laboratory to perform data entry, and have trained undergraduate students as well. Currently, our database includes only a small fraction of all the relevant data that exist. In order to increase the number of measurements in our database, we will need to increase the size of our data entry team. We are considering crowdsourcing the process so that verified database users from across the world can also contribute. In addition, verified users will have the option of submitting CIF files for automatic calculation of electronic structure and tight-binding parameters. We would also like to populate the database with results from explicitly correlated theoretical methods suitable for high-throughput applications (such as density matrix embedding theory [45]). Presently, a tight-binding electronic structure is provided for each crystal structure in the database. For some of these crystal structures, missing or corrupt atoms in the anion layer prevent the use of a full LAPW DFT calculation; however, we eventually plan to include full LAPW DFT calculations for as many of the crystal structures in the database as possible. Any persons interested in becoming involved with the project can click on the orange button at the bottom of the home page to request an account. An option is also available to provide anonymous feedback.

Our second main goal regarding development of the database is to implement new features. We plan to add an interactive web interface by which users can analyze the contents of the database as a whole by correlating various calculated and measured properties. We also would like to improve the search feature of the website so that more detailed searches can be performed. There are many different avenues by which our existing work can be further developed. We are trying to create as many tools as possible to perform simple visualization and analysis of data online, such as the feature shown in Figure 5. A number of parameters can be extracted from the band structure calculations, such as the density of states at the Fermi level, and the cross sectional area of the Fermi surface. We are working on finding robust universal algorithms to calculate these and other representative values. Given this collection of the unit cell parameters, extracted electronic values, and measurement parameters, any set of data can be graphed against any other set of data, and scatter plots can be created including markers labeled with the material names. It is also possible to combine parameters with common arithmetic operations to create additional metrics. We will continue to enhance the user interface to create a more versatile and expansive analysis interface.

Following the invention and widespread availability of computers, an increasing trend towards digitization in science has taken place. Scientific databases have emerged in practically every area of study because they enable the analysis of many pieces of information, and the distribution of that information to individuals around the world. The field of data science has also grown to develop new ways of analyzing the large amount of data available. Many databases for materials science currently exist, particularly those that focus on electronic structure, crystal structure, and other measured properties [46]. Our inspiration to create this database of organic conductors and superconductors was drawn in part by the success of other databases containing density functional theory electronic structure calculations for many crystal structures [47–50]. Our goal is to gather as many different types and pieces of information related to quasi two-dimensional organic conductors and superconductors as possible. We foresee this database as having a number of different applications. Primarily, it will serve as a useful reference tool for the scientific community of organic conductors and superconductors. Database users can easily find and view crystal structures, electronic structures, and other measured properties for the materials that they are interested in. We hope to cultivate a community of scientists from across the world who are interested in using the database, contributing data to the database, and requesting new features for the database. We also look forward to analyzing the data that are stored in the database. Many techniques in the field of data science are appropriate for this purpose. In particular, certain types of data mining and machine learning have proven useful in the analysis of other materials databases [51–54]. We hope that the identification of trends between various parameters will ultimately lead to a better

understanding of the fundamental mechanisms of correlated electron systems in quasi two-dimensional organic conductors and superconductors.

**Figure 5.** View of the 3D in-browser crystal structure analysis tool. The section on the right provides interactive features to the users. In this example the material is *β*"-(ET)4[(H3O)Cr(C2O4)3]2[(H3O)2]5H2O, as reported in [41].

**Author Contributions:** Conceptualization, O.G. and C.C.A.; Methodology, O.G., K.F., M.B.-d. and S.M.W.; Software, O.G.; Validation, K.F. and M.B.-d.; Investigation, O.G.; Resources, C.C.A. and S.M.W.; Data Curation, K.F., M.B.-d. and O.G.; Writing—Original Draft Preparation, O.G., C.C.A. and S.M.W.; Supervision, C.C.A.; Project Administration, C.C.A.; Funding Acquisition, C.C.A. and S.M.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by NSF grant number DMR-1905950. S.M.W. and O.G. also acknowledge support from a pilot grant from the Center for Functional Materials, Wake Forest University.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available in the described database, accessible via osd.clarku.edu.

**Acknowledgments:** The authors acknowledge discussions with Ben Powell. The authors thank the following people for their work in manually adding measurement entries to the database: Alireza Alipour, Abdulai Gassama, Ahad Ali Khan, Brett Laramee, Gwynnevieve Ramsey, Jade Consalvi, Meherab Hossain, and Raju Ghimire. Computations were performed on the Clark University High Performance Computing Cluster and the Wake Forest University DEAC Cluster, a centrally managed resource with support provided in part by Wake Forest University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **How to Recognize the Universal Aspects of Mott Criticality?**

**Yuting Tan 1,\*, Vladimir Dobrosavljevi´c 1,\* and Louk Rademaker 2,\***


**Abstract:** In this paper we critically discuss several examples of two-dimensional electronic systems displaying interaction-driven metal-insulator transitions of the Mott (or Wigner–Mott) type, including dilute two-dimension electron gases (2DEG) in semiconductors, Mott organic materials, as well as the recently discovered transition-metal dichalcogenide (TMD) moiré bilayers. Remarkably similar behavior is found in all these systems, which is starting to paint a robust picture of Mott criticality. Most notable, on the metallic side a resistivity maximum is observed whose temperature scale vanishes at the transition. We compare the available experimental data on these systems to three existing theoretical scenarios: spinon theory, Dynamical Mean Field Theory (DMFT) and percolation theory. We show that the DMFT and percolation pictures for Mott criticality can be distinguished by studying the origins of the resistivity maxima using an analysis of the dielectric response.

**Keywords:** Mott transition; quantum criticality; resistivity maxima; dielectric response; dilute 2DEGs; Mott organics; twisted transition-metal dichalcogenide bilayers; dynamical mean field theory; percolation theory; spinon theory

#### **1. Introduction**

The physics of strongly correlated matter has many faces. Still, for a majority of systems the underlying theme is the role of "*Mottness*" [1]. It is clear that if one aspect of strong correlations should be understood first, it should be the fundamental nature of the Mott metal-insulator transition [2]. Its simplest reincarnation is the transition induced by tuning the bandwidth at half-filling, a setup that produced rather spectacular advances in recent years. Several systems were identified as nearly-ideal realizations of this paradigm, allowing systematic study using a wide arsenal of experimental probes.

In this article we present an overview of three classes of two-dimensional experimental systems that exhibit bandwidth-controlled Mott criticality: dilute two-dimensional electron gases in semiconductors, "Mott-organic" compounds, and transition-metal dichalcogenide moiré systems. Thereby we aim to present the experimental facts as "bland" as possible, in Section 3, without favoring one or the other theoretical explanation. The remarkable similarities between these model systems suggests a robust universality, including characteristic behavior such as the appearance of resistivity maxima.

Possible explanations of two-dimensional Mott criticality follow in the section thereafter (Section 4), where the experimental distinguishable features of each theory takes the forefront. This is followed by a separate discussion of the largely-overlooked utility of dielectric spectroscopy in Section 5, in not only identifying phase segregation and spatial inhomogeneity, but also in revealing the thermal destruction of coherent quasiparticles associated with Landau's Fermi liquid theory.

However, first we need to address the demarcation of our topic. What makes the metal-to-insulator transition in these systems stand out from 'traditional' metal-to-insulator transitions [3,4]?

**Citation:** Tan, Y.; Dobrosavljevi´c, V.; Rademaker, L. How to Recognize the Universal Aspects of Mott Criticality? *Crystals* **2022**, *12*, 932. https://doi.org/10.3390/ cryst12070932

Academic Editor: Dmitri Donetski

Received: 4 May 2022 Accepted: 3 June 2022 Published: 30 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **2. In Search of Mott Criticality**

Condensed matter physics, or recently for sales purposes re-branded as "quantum matter physics", is the study of electric and magnetic properties of materials that surround us. The grand question that trumps all others is: how to understand which materials conduct electricity and which are insulating? Traditionally, a conducting material is called a *metal*—not to be confused with the chemical, metallurgical or astronomical meaning of that word. Our main question (metal or insulator?) has not only tremendous technological applications (in fact, all modern electronic technology depends on our ability to rapidly switch materials between metallic and insulating behavior), but also requires a thorough understanding of the problem of emergent behavior of many interacting quantum-mechanical electrons and ions.

Only at zero temperature does there exist a sharp difference between insulators and metals [4]. There are three distinct possibilities: zero conductivity *σ*(*T* = 0) = 0 means insulating; zero resistivity *ρ*(*T* = 0) = 0 means a superconductor; and anything in between is a metal *σ*(*T* = 0) = 1/*ρ*(*T* = 0) = 0. At any nonzero temperature, an insulator typically has activated behavior *<sup>ρ</sup>*(*T*) <sup>∼</sup> *<sup>e</sup>*Δ/*<sup>T</sup>* whereas the standard Fermi liquid theory of a metal predicts a temperature-squared increase of the resistivity *ρ*(*T*) = *ρ*<sup>0</sup> + *AT*2. It has therefore become common-place to use the *derivative* of the resistivity *dρ*/*dT* as a measure of whether something is conducting (*dρ*/*dT* > 0) or insulating (*dρ*/*dT* < 0)—but this is highly misleading! As we will show later in Section 3, close to a Mott metal-insulator transition we often find non-monotonic behavior of the resistivity as a function of temperature, making the 'derivative' criterion useless. Even worse, there exist cases where the resistivity has *dρ*/*dT* < 0 but at zero temperature it does not diverge, signalling that this is not a true insulator (see e.g., [5,6]). Another example is the case of Mooij correlations [7,8], where the temperature-derivative of the resistivity in a metal can become negative. Consequently, since only at zero temperature the insulator/metal distinction is well-defined, we must stick with that definition. Regardless of the slope, a material is a metal if its resistivity *does not diverge* as *T* → 0.

Many materials can be understood within the framework of *band theory* and its extensions such as Fermi liquid, Boltzmann transport, and density functional theory. This framework provides a very simple answer to the metal-or-insulator question: if the Fermi level lies in the middle of a band gap, the system is insulating; otherwise, the system behaves as a metal. This concept has the important consequence that for a crystalline material with (up to some weak disorder) well-defined unit cells, insulators can only appear when there is an even number of electrons per unit cell. Consequently, within the band theory picture, there exist only three possible routes to induce a metal-to-insulator transition: by changing the electronic density; via spontaneous symmetry breaking; or via band overlap when the filling is even. An example of the first is doping a semiconductor, which is the metal-insulator transition we induce on a daily basis inside transistors. An example of the second is the transition into antiferromagnetic ordering: when the system is at half-filling of a band (meaning one electron per unit cell), after antiferromagnetic unit cell doubling there are two electrons per unit cell, and the system can become a band insulator. The third case can be realized by for example straining a system such that the band gap changes from positive to negative.

There are, however, two main exceptions to the paradigm of band theory. On the one hand, disorder can become so large as to prevent the motion of the charge carriers—this is known as Anderson localization [9]. On the other hand, the presence of very strong electron-electron interactions can force the electrons to become "stuck" like in a traffic jam—this is known as *Mott insulation* [2]. The standard model of Mott insulation is the Hubbard model with a tight-binding Hamiltonian:

$$H = -t\sum\_{\langle ij\rangle\sigma} c\_{i\sigma}^{\dagger} c\_{j\sigma} + lI \sum\_{i} n\_{i\uparrow} n\_{j\downarrow\prime} \tag{1}$$

where *t* is the nearest-neighbor hopping on some lattice and *U* is the onsite repulsion. When *U* = 0, the system is a metal when half-filled. When *U t*, it becomes energetically favorable to occupy each site with exactly one electron rather than to fill bands up to the Fermi level. The resulting Mott state can therefore not be described by band theory!

Mott insulators have been observed in a wide variety of materials, most famously transition-metal oxides, including high *Tc* superconducting cup rates [3]. Observing a clear transition from a standard Fermi liquid metal to a Mott insulator, however, is quite elusive. This transition can be induced either by changing the electronic density ("filling-controlled") or by changing the ratio *U*/*t* ("bandwidth-controlled"). The filling-controlled Mott transition [10] notoriously leads to a whole zoo of different instabilities, pseudogaps, and strange metal behavior, and is typically masked at low temperatures by superconductivity. The bandwidth-controlled Mott transition is, in contrast, often masked by (antiferromagnetic) spin order that hides any Mottness behind the veil of unit cell doubling.

This might, at first, suggest that *Mott criticality* is something unattainable. By "criticality" we mean that approaching the Mott transition we find vanishing energy scales, and that the resistivity curves display scaling behavior. There are, however, two clever tricks to realize Mott criticality. The first trick is *dimensionality*: a transition that is strongly first-order in *d* = 3 dimensions often becomes continuous or weakly first-order in *d* = 2 dimensions. The most striking example of this is, of course, the solidification of 3He. The second, and perhaps even more important trick is *frustration*: if the lattice structure is highly frustrated (with competing magnetic interactions) one can avoid [11] antiferromagnetic ordering altogether—revealing the true Mott transition.

In this review we, therefore, focus on three classes of systems that are indeed (quasi) two-dimensional as well as frustrated: Wigner crystals in extremely dilute two-dimensional electron gases; layered Mott organic compounds; and the more recent addition of transitionmetal dichalcogenide (TMD) moiré bilayers. Indeed, as we will show in Section 3, these systems all seem to exhibit remarkably similar distinct features, including clear signatures of critical resistivity scaling. Because these systems all have a fixed electron density per unit cell of *n* = 1 (at least in the insulating limit), the observed transitions are plausibly within the universality class of bandwidth-tuned Mott transitions.

A brief side-note is in order: we briefly mentioned superconductivity and disorderinduced insulators. These phases can also have a continuous transition between them, the so-called superconductor-to-insulator transition [12,13]. This, however, is an interesting topic that falls outside the scope of this review. Similarly, we also will not consider disorderdriven metal-insulator transitions [14,15], since this regime typically does not include any Mottness. More general but also somewhat older reviews of metal-insulator criticality can be found in Refs. [3,4,16–18].

#### **3. Experiments**

Given that experimental results should always be leading, the aim of this section is to introduce three material systems that are likely exhibiting a bandwidth-tuned Mott metalinsulator transition: dilute 2DEGs, organics, and moiré systems. To support the clarity of interpretation, we will stress *experimental* similarities between these systems without much room for theoretical guesswork—that is the next section's realm.

While each system has a different tuning parameter (density, pressure, or field), the *electrical resistivity* through the transition is the key observable, see Figure 1. Its behavior reveals how the transport gap Δ decreases when we approach the transition from the insulating side; as well as how the resistivity behaves on the metallic side, where Fermi liquid behavior *ρ*(*T*) = *ρ*<sup>0</sup> + *AT*<sup>2</sup> is typically seen at *T* < *TFL* with an enhanced effective mass *m*∗. Remarkably, in all systems one also observes distinct *resistivity maxima* at *T* ∼ *Tmax* > *TFL*, signalling the breakdown of coherent transport. Crossover to the quantum critical regime is described by an additional temperature scale *To*, which is extracted from the scaling collapse of the resistivity curves as shown in Figure 2.

**Figure 1.** The key observable revealing a metal-insulator transition is the resistivity. Here we show *ρ* vs. *T* resistivity curves as a function the tuning parameter, for representative examples of the three material systems considered: (**a**) 2DEG in Si-MOSFET tuned by electronic density (reprinted with permission from Ref. [19] Copyright 2019 American Physical Society), (**b**) Mott organic material *κ*-(BEDT-TTF)2Cu2(CN)3 tuned by pressure [20], and (**c**) TMD moiré bilayer MoTe2/WSe2 tuned by displacement field (Data imported from [21]). In all cases, one observes distinct resistivity maxima on the metallic side, at a temperature *Tmax* that decreases towards the transition.

**Figure 2.** Critical scaling has been observed in all three experimental systems, when the resistivity is plotted versus *T*/*To* where *To* is the characteristic crossover (quantum critical) energy scale. Note that in all cases a strong "mirror" symmetry [22,23] exists between the insulating (upper) and metallic (lower) scaling branch. (**a**) In a dilute 2DEG, scaling of the bare resistivity *ρ*(*T*) was achieved by simply rescaling *T* with *To* ∼ |*δ*| 1.6 (Adapted with permission from Ref. [24] Copyright 1995 American Physical Society); (**b**) In organic compounds, the normalized resistivity *ρ*˜ is obtained by normalizing the resistivity by the critical resistivity along the Widom line. This leads to excellent scaling collapse with *To* ∼ |*δ*| 0.60±0.01 (Adapted with permission from Ref. [25] Copyright 2015 Springer Nature); (**c**) A similar approach was followed in TMD moiré bilayer MoTe2/WSe2, with similar *To* ∼ |*δ*| 0.70±0.05 (Data imported from [21]).

A practical summary of the experimental results is presented in Table 1.

**Table 1.** A summary of available experimental results for the three classes of systems considered. The sources (references) are given in the text below. Question-marks indicate the lack of reliable data. Fermi liquid (*T*2) transport behavior has not been documented in 2DEG systems, in contrast to strong evidence for it in Mott organics and TMD moiré bilayers. Note that the characteristic energy scales Δ, (*m*∗)−1, *TFL*, *Tmax*, as well as *To* display similar continuous decrease towards the transition in all three systems, consistent with general expectations for quantum criticality. One should keep in mind that the error bars on the estimated exponent could be substantial, since the results typically depend strongly on the utilized fitting range.


#### *3.1. Dilute 2DEG in Semiconductors*

In dilute two-dimensional electron gases (2DEG) [26], the electron density can be quantified by the dimensionless parameter *rs* = 1/ <sup>√</sup>*πnaB* where *<sup>n</sup>* is the electron density and *aB* the Bohr radius. The ratio of interaction energy versus kinetic energy scales as *rs*, and therefore at large enough *rs* (of the order *rs* ∼ 40 in 2D) the electrons will spontaneously crystallize into a Wigner solid. In a two-dimensional Wigner crystal, the electrons form a triangular lattice with exactly one electron per unit cell—essentially forming a frustrated Mott insulator. When the electron density *n* is varied, the size of the unit cell changes accordingly so that the Wigner crystal always remains fixed at one electron per unit cell. The transition from an insulating Wigner crystal to a metal can therefore be plausibly viewed as a bandwidth-tuned Mott transition. Note that this is counter-intuitive: after all, one tunes the electron density! However, what matters is the electron density *counted per unit cell* and that remains constant. This idea suggests [27–29] that the melting of a Wigner solid by increasing density should be viewed as a Wigner–Mott transition, possibly bearing many similarities to Mott transitions in narrow-band crystalline solids such as Mott organics or transition-metal oxides. If this viewpoint is correct, then the resulting metal above the transition should display resemble other strongly correlated Fermi liquids, a notion that is starting to gain acceptance on the base of recent experiments [30–32].

Experimentally, high-quality 2DEGs can be realized in metal-oxide-semiconductor fieldeffect devices (MOSFETs) in various semiconductors [24,30,32,33]. Through electrostatic gating the electronic density can be elegantly tuned, typically in the range of *<sup>n</sup>* ∼ 1010–1012 cm<sup>−</sup>2. The peak electron mobility in ultra-clean samples can be as high as 10<sup>4</sup> cm2/V s [19], which implies that down to very low temperatures the transport properties are dominated by electron-electron interactions (like Wigner crystallization) rather than extrinsic disorder effects. Lower-mobility devices have also been extensively studied (for a review see Chapter 5 of Ref. [18]), displaying different types of metal-insulator transitions displaying electron glass dynamics [34], which we will not discuss here.

Indeed, tuning the electronic density leads to insulating transport below a critical density, typically around *nc* ∼ <sup>10</sup><sup>11</sup> cm−<sup>2</sup> [24,33], see Figure 1a. Activated behavior is often observed close to the transition [34,35], with the activation energy Δ ∼ |*n* − *nc*|. Further on the insulating side disorder effects may become important, where Efros-Shklovskii hopping

(and other effects of disorder) is also observed [24], but only at the lowest temperatures. On the metallic side, a pronounced resistivity drop (often by a factor of 10 or more) is observed [24,30] below the temperature *Tmax* ∼ |*n* − *nc*| which decreases as the transition is approached. Characteristic scaling of the resistivity maxima has been reported in several systems [30–32], see Figure 3, which has been interpreted as evidence for strong correlation effects. However, the expected *T*<sup>2</sup> dependence of the resistivity has not been observed, despite the reported effective mass enhancement (*m*∗)−<sup>1</sup> ∼ |*<sup>n</sup>* − *nc*| [36], characteristic of correlated Fermi liquids. Quantum critical scaling collapse of the resistivity curves has also been demonstrated [24] around the critical density, albeit excluding the lowest temperatures data, as shown in Figure 2a. This is achieved by rescaling *T* by a crossover scale *To* ∼ |*n* − *nc*| *<sup>ν</sup>z*, with *<sup>ν</sup><sup>z</sup>* <sup>≈</sup> 1.6. The resulting scaling function reveals surprising "mirror symmetry" [22], which was phenomenological interpreted [23] as evidence "strongcoupling quantum criticality". Similar systems to these 2DEGs include the observation of a Wigner crystal in low-density doped monolayer WSe2 [37], where more detailed experiments still need to be performed.

**Figure 3.** Characteristic scaling of the resistivity maxima has been reported in several 2DEG electron systems in semiconductors: (**a**) Si-MOSFETs (adapted with permission from Ref. [29] Copyright 2012 American Physical Society); (**b**) p-GaAs/AlGaAs quantum wells (adapted with permission from Ref. [29] Copyright 2012 American Physical Society); (**c**) SiGe/Si/SiGe quantum wells (adapted with permisssion from Ref. [30] Copyright 2020 American Physical Society); (**d**) few layered-*MoS*<sup>2</sup> (Adapted with permission from Ref. [32] Copyright 2020 American Physical Society). All data collapse to the same (theoretical) scaling function [29] obtained from the Hubbard model at half-filling, in the vicinity of the Mott point.

#### *3.2. Organic Compounds*

An organic compound [38] refers to a crystalline system where each unit cell contains an entire *molecule*, rather than just loosely bound ions. A particularly interesting class of organic compounds is based on the molecule bis-(ethylendithio)-tetrathiafulvalen (BEDT-TTF or ET), which can be fabricated with other ions into quasi-two-dimensional layered systems. Compounds based on BEDT-TTF exhibit a spectrum of interesting quantum matter phenomena, ranging from superconductivity [39] to electron glass behavior [40].

Our interest goes out especially to *κ*-(ET)2Cu[N(CN)2]Cl and *κ*-(ET)2Cu2(CN)3, where the molecules are organized in triangular lattice layers [38]. These materials are strongly correlated, and indeed, despite being half-filled they are insulating at ambient pressures. Due to the geometric frustration of the triangular lattice [11], no magnetic order has been observed in *κ*-Cu2(CN)3 and antiferromagnetic order only at relatively low temperatures (*T* < *TN* ≈ 20*K*) in *κ*-Cl. The absence of magnetic order is the strongest indication that *κ*-Cu2(CN)3 might realize a spin liquid ground state [41,42].

Upon applying pressure, a zero-temperature *first-order* phase transition brings the system into a paramagnetic metallic phase at *pc* = 122 MPa (*κ*-Cu2(CN)3) or *pc* = 24.8 MPa (*κ*-Cl), see Figure 1a [20,25]. The first order phase boundary ends in a critical point at *Tc* = 20 K or *Tc* = 38 K, respectively. It is important to emphasize that these temperatures are very small compared to the electronic energy scales. The Hubbard repulsion *U* and bandwidth *W* are both on the order of a fraction of eV [42], which implies *Tc U*, *W*. As such, even though the observed Mott criticality appears at nonzero temperatures, much of the observed phenomena above *Tc* can be described *as if* the system resides in the vicinity of a quantum critical point [25].

Above *Tc*, a crossover pressure *Pc*(*T*) can traced [25] where the measured resistivity exhibits an inflection point, see Figure 4. This defines the *"quantum Widom line"* (QWL) [43] by analogy to the standard liquid-gas crossover. Defining the critical resistivity *ρc*(*T*) to be the resistivity along the QWL, all resistivity curves collapse onto each other when plotted as *ρ*(*P*, *T*)/*ρc*(*T*) vs. *T*/*T*0(*P*), as shown in Figure 2b. Here the scale *To*(*P*) reflects a critical energy scale that vanishes at the critical pressure, *To* ∼ |*P* − *Pc*| 0.6 for Cu and *To* ∼ |*<sup>p</sup>* − *pc*| 0.5 for *κ*-Cl [25]. On the insulating side of the transition, the resistivity is approximately activated *ρ* ∼ exp(Δ/*T*) [42]. On the metallic side, it follows [44] the standard Fermi liquid behavior at low temperatures *ρ*(*T*) = *ρ*<sup>0</sup> + *AT*2, up to a temperature scale *TFL*, see Figure 5a. This destruction of the Fermi liquid seems to correspond to the appearance of a maximum in the resistivity [20].

**Figure 4.** Finite temperature phase diagram of the Mott organic materials. (**a**) first-order phase transition line, as observed in *κ*-Cu2(CN)3 (adapted with permission from Ref. [25] Copyright 2015 Springer Nature) at *T* < *Tc* ∼ 20*K*, displaying "Pomeranchuk" behavior [45], by "sloping" towards the metallic phase. The corresponding "Quantum Widom Line" [43] arises at *T* > *Tc*, which is identified as the center of the quantum critical region [46] with resistivity scaling [25]. (**b**) Phase diagram [44] for *κ*-Cu2(CN)3 over a broader *T*-range, displaying the convergence of the quantum Widom line (QWL) on the insulating side, and the "Brinkman-Rice" line (*TBR* = *Tmax*, which intersect at the critical end-point *T* = *Tc*. The Fermi-Liquid line *TFL* < *TBR* is also shown. (**c**) The universal phase diagram for a series of spin-liquid Mott organics compounds was established [42] by rescaling the temperature *T* and the interaction strength *U* by the respective electronic bandwidth *W*. The parameters *W* and *U* were independently measured [42] for each material using optical conductivity.

In addition to transport measurements, and in contrast to other systems we consider, Mott organics have also been carefully investigated using optical probes. This allowed to directly identify [42] the quantum Widom line, which is back-bending towards the insulating side at higher temperatures following the closing of the Mott gap. In addition, the "Brinkman–Rice" line traced by *Tmax* was identified as marking the thermal destruction of Landau quasiparticles, as seen by the vanishing of the Drude peak in the optical conductivity [44], see Figure 6.

Finally, the controversy about the presence or absence of the low-*T* phase coexistence region has been resolved in Mott organics, by using dielectric spectroscopy [20]. Its precise location on the phase diagram has been identified by the observation [20] of colossal dielectric response, as a smoking gun for percolative phase coexistence. In addition, the same technique was able to demonstrate the coincidence of the resistivity maxima in the (uniform) metallic phase, with the thermal destruction of Landau quasiparticles. This is seen as a dramatic drop and a change [20] of sign of the dielectric function at *T* < *TBR* = *Tmax*. These experimental results are shown in Figure 6, and discussed in more detail in Section 5.

**Figure 5.** Fermi liquid behaviour at low temperatures, for (**a**) Mott organic material *κ* − [(BEDT-TTF)1−*x*(BEDT-STF)*x*]2Cu2(CN)3 [44] and (**b**) MoTe2/WSe2 moiré bilayers (Data imported from [21]). Clear *ρ* = *ρ*<sup>0</sup> + *AT*<sup>2</sup> behavior is observed in both cases, up to a temperature scale *TFL* that seems to decrease linearly towards the metal-insulator transition. The resistivity curves can be collapsed by plotting *ρ*(*E*, *T*)/*ρc*(*T*) vs. *T*/*T*<sup>0</sup> where *T*<sup>0</sup> ∼ |*E* − *Ec*| 0.70±0.05, see Figure 2c. Note that this crossover scale seems to follow both the gap size on the insulator, as well as the destruction of the Fermi liquid on the metallic side.

**Figure 6.** Transport behavior vs. dielectric response across the phase diagram of *κ*-Cu2(CN)3 [20]. (**a**) DC transport shows only very gradual change across the BR line (resistivity maxima), and cannot one see any clear indication of the phase coexistence region. (**b**) In dramatic contrast, the lowfrequency dielectric function <sup>1</sup> assumes small positive values in the Mott insulator (pale pink), and large negative values in the quasiparticle regime (deep blue); we clearly see the boundaries of these regimes tracing the QWL and the BR line (following *Tmax*), as observed in transport. Remarkably, "resilient" quasiparticles [47] persist past the Fermi Liquid line, at *TFL* < *T* < *TBR* = *Tmax*, where bad metal behavior [48] (metallic transport above the Mott-Ioffe-Regel [49] limit is observed). At low temperature, the Mott point is buried below the phase coexistence dome, which is vividly visualized through colossal dielectric response (<sup>1</sup> <sup>∼</sup> <sup>10</sup>3–104).

#### *3.3. Moiré Materials*

The most recent addition to the field of strongly correlated systems are *moiré materials*. These are bilayer structures made of Van der Waals materials such as graphene and transitionmetal dichalcogenides (TMDs). A lattice mismatch or relative twist angle between the layers causes a large-scale geometric "moiré" pattern. This larger unit cell (typically in the range of 5–10 nm) drastically reduces the effective electron kinetic energy such that the bandwidth is on the order of *W* ∼ 10 meV. As a result, the systems become strongly correlated, with *U*/*W* ∼ 10 or larger. Using electrostatic gating one can tune the electronic density, typically in the range of a few electrons or hole per moiré unit cell (corresponding to *<sup>n</sup>* ∼ 1012 cm<sup>−</sup>2).

While the most famous of correlated moiré materials is without a doubt twisted bilayer graphene, more convincing evidence for Mott correlations has so far only been observed in TMD bilayers. Here, we focus on one particular system: the heterobilayer MoTe2/WSe2 [21]. The lattice mismatch between MoTe2 and WSe2 gives rise to a moiré period of *aM* ∼ 5 nm. At half-filling of the first valence band, an insulating phase appears which can be tuned into a metal by applying a vertical displacement field *E*. This flat valence band can be described by a spin-orbit coupled triangular lattice Hubbard model, where the displacement field *E* tunes the bandwidth [50].

The temperature-dependent resistivity across the transition is shown in Figure 1c. On the insulating side, the system has well-defined activated behavior of the resistivity with a gap Δ continuously vanishes as the critical displacement-field value is approached, Δ ∼ |*E* − *Ec*| 0.60±0.05. At the critical point, the resistivity is claimed to follow a powerlaw, *<sup>ρ</sup><sup>c</sup>* ∼ *<sup>T</sup>*<sup>−</sup>1.2, although the reliability of the low-*<sup>T</sup>* data may be questionable. On the metallic side, the low-*T* resistivity follows the Fermi liquid law *ρ*(*T*) = *ρ*<sup>0</sup> + *AT*<sup>2</sup> where the quadratic prefactor diverges *A* ∼ |*E* − *Ec*| <sup>−</sup>2.8±0.2. This clear Fermi liquid does not persist up to all temperatures, instead, a resistivity maximum appears at *Tmax* ∼ |*E* − *Ec*| 0.70±0.05, see Figure 5b. No magnetic order has been observed, which might be due to the geometric frustration of the triangular moiré lattice structure. Remarkably, this experiments uses the recently-developed *excitonic sensor* [21], which allows the measurements of the spin susceptibility across the transition. This reveals Curie-law behavior over a broad temperature range, thus demonstrating the presence of localized magnetic moments, as expected for a Mott system.

Finally, among other moiré systems it is worth mentioning the twisted hetero-bilayer WSe2 [51]. Though it was claimed to exhibit some sort of Mott-related criticality, it is not certain whether a true insulating phase has indeed been observed given that the activated transport does not continue down to the lowest measured temperatures. In addition, insulating behavior seems to disappear above a relatively low of the order of *T*<sup>∗</sup> ∼ 5–10 K, which is similar to twisted bilayer graphene, but much smaller than the estimated bandwidth of the order of *W* ∼ 100–200 K. Furthermore, no clear resistivity maxima on the metallic side have been observed. Whether or not tbWSe2 can be classified as a true Mott insulator is therefore quite controversial. Alternatively, the observed behavior could be a result of some sort of magnetic order, which may arise close to half filling even in weakly-coupled systems.

#### *3.4. Universal Criticality*

The resistivity curves of the three systems, as shown Figures 1 and 2, show remarkable universality, reflected in the fact that all curves can be collapsed by scaling with *T*<sup>0</sup> ∼ |*δ*| *zν* where *δ* is the tuning parameter and *zν* the critical exponent. It is important to realize, however, that the precise scaling procedure applied was not identical in the three cases, and the resulting critical exponent also somewhat depend on the system.

So what is *different* between these systems? Let us first focus on the energy scale. The typical bandwidth *W* ranges from ∼100 s meV in Mott organics, to ∼10 s meV in the moiré systems, to 0.1–1 meV in the dilute 2DEGs. A finite temperature critical point is only observed in the organics, though at about *Tc* ∼ 1%*W*—which leaves open the possibility that a finite *T* critical endpoint exists in the other two systems. Indeed, most experiments

are performed (so far) above the Kelvin range in moiré materials, and above 100 mK in 2DEG, which makes it hardly possible to reliably explore the *T*-range below few percent of the bandwidth.

Secondly, in order to achieve quantum critical scaling in Mott organics, one needs to first identify a Widom line as a demarcation of the finite-temperature crossover from insulator to metal. A similar analysis was carried out for moiré materials, although the obtained Widom line displayed no apparent "curving" as a function of temperature. This is manifestly not performed in the dilute 2DEGs. It might be interesting to see whether better collapse can be achieved through such a method.

Thirdly, with a bit of good-will, the critical exponents in organics and Moiré systems are in the same ballpark; whereas the critical exponent in the dilute 2DEGs with *zν* = 1.60 is significantly larger. It is also important to realize that in 2DEGs there has not been a clear observation of a Fermi liquid regime—unlike in organics and Moiré systems, see Figure 5. These ways in which 2DEGs stand out might be related to the fact that there is no underlying (Wigner) lattice on the metallic side, which could point to a perhaps nontrivial role of significant charge density fluctuations on the metallic side, an effect not present in lattice Mott systems.

#### **4. Competing Theoretical Pictures**

As we mentioned in Section 2, the observation of Mott criticality and scaling opens big questions on the theoretical front. Currently, there exist *three* main different physical pictures to address these issues.

A true Mott transition should not be hidden by some period doubling symmetry breaking. The Lieb-Schultz-Mattis theorem states that in absence of spin order, the ground state of Mott insulator must be a spin liquid [52]. This leads directly to the first theoretical picture: Mott criticality can only occur if the Mott phase is a spin liquid, where inter-site spin correlations play an important role. The theory of Senthil [53] chooses this path, by introducing an explicit spinon theory of the Mott spin liquid.

Alternatively, one focuses on local electronic processes only, ignoring inter-site spin correlations. Then the Mott transition at low temperature becomes first-order; however, it is only *weakly* first order. A first order transition line always ends at a critical point *Tc*, and as long as *Tc* is sufficiently low compared to any experimental scale, one still finds criticality and scaling. This is the picture emerging from Dynamical Mean Field Theory (DMFT) [54], a strong-coupling self-consistent approach to calculate the local electronic self-energy.

The third picture again accepts the first-order nature of a Mott transition, but this time embraces it. A first-order transition is always accompanied by a region where both phases coexist. Minor disorder or self-generated pattern formation [55,56] can smear this phase coexistence region into a continuous-looking transition exhibiting nontrivial electron dynamics. This is the 'percolation theory' picture of Mott criticality.

The goal of this review paper is to put the main theoretical predictions next to the experimental findings. As such, we will not dive into the pros and cons of each theoretical picture. A summary of the main theoretical predictions is provided in Table 2.


**Table 2.** A summary of predictions from competing theoretical pictures. The expected transition type differs between the three pictures, with observable differences in the behavior of the mass enhancement *m*∗, the Kadowaki–Woods ratio *A*/(*m*∗)2, the destruction of the Fermi liquid at *TFL*, and the appearance of a resistivity maxima at *Tmax*. Details are provided in the text below.

#### *4.1. Spin Liquid Picture of the Mott Point*

A popular approach to describe a spin liquid state is through *spin-charge separation*. In Ref. [53], the electron is split into a charge-0 spin-1/2 fermionic *spinon f* and a charge-*e* spin-0 bosonic *chargon b*. The Mott transition, in this picture, amounts to the condensation of the chargon field, whose critical behavior falls within the 3D *XY* universality class. The Fermi liquid corresponds to the condensed phase of the chargon, whereas the Mott insulator corresponds to a gapped phase of the charged boson. The splitting of the electron leads to redundant degrees of freedom described by an emergent gauge field. Fluctuations of this gauge field lead to a logarithmic enhancement of the quasiparticle effective mass,

$$m^\* \sim \ln \frac{1}{|\mathcal{g} - \mathcal{g}\_c|},\tag{2}$$

where *g* is the tuning parameter and *gc* is the critical value. However, as in any theory with a non-local electronic self-energy, the quasiparticle residue *Z* is *not* simply proportional to the inverse effective mass; instead *Z* ∼ |*g* − *gc*| *<sup>β</sup>*/ ln <sup>1</sup> |*g*−*gc* | . Furthermore, approaching the Mott transition from the metallic side the spin susceptibility *χ* remains constant whereas the compressibility *κ* vanishes. Physically, these effects result from important inter-site spin correlations, where a gapless spin liquid can be viewed as a certain superposition of spin singlets formed by pairs of spins in the Mott insulating state. As a result, there emerges a finite gap *δ* to charge excitations, while the rearrangement of singlets leads to characteristic gapless spin excitations with fermionic quasiparticles. This picture is a specific realization of the famous RVB picture of Baskaran and Anderson [57], first proposed in the context of high-*Tc* superconductors.

Another significant consequence of describing the Mott transition as chargon condensation, is that the *T* = 0 conductivity is not continuous. The electron resistivity will display a *universal* jump from a (disorder)-dependent constant value *ρ* = *ρ*<sup>0</sup> in the Fermi liquid; to *ρ* = *ρ*<sup>0</sup> + *Rh <sup>e</sup>*<sup>2</sup> (with *<sup>R</sup>* of order one) at the critical point; to *<sup>ρ</sup>* = <sup>∞</sup> in the Mott insulator. On the metallic side, the Fermi liquid is predicted to break down above *TFL* ∼ |*g* − *gc*| <sup>2</sup>*<sup>ν</sup>* and give rise to a *marginal* Fermi liquid state, which in turn survives up to *TMFL* ∼ |*g* − *gc*| *ν*. In both cases, *ν* = 0.67 is the 3D *XY* correlation length exponent. On the insulating side, the boson condensation picture implies that the charge gap vanishes as Δ ∼ |*g* − *gc*| *<sup>ν</sup>*. The spinons, however, remain gapless and form a spinon Fermi surface, with low-temperature specific heat scaling as *<sup>C</sup>* ∼ *<sup>T</sup>*2/3.

Note that the original work in Ref. [53] does not directly provide a detailed description for finite temperature dependence of the resistivity, and thus no explicit prediction for a possible deviation from the Kadowaki-Woods (KW) law (*A*/(*m*∗)<sup>2</sup> ≈ constant) [58]. On the other hand, Ref. [58] presents arguments that the physical requirement for the validity of the KW law is the locality of the electronic self-energy (as in DMFT theory) [54], a condition which is not obeyed by the RVB-type spin-liquid theories such as the Senthil's spinon picture.

It should be stressed that the spinon theory makes one sharp prediction about finite temperature transport in the critical regime. Namely, the *critical resistivity curve* is predicted to assume a universal power-law form *ρc*(*T*) ∼ 1/*T* in *d* = 3 [59], but remain *T*-independent in *d* = 2, see Figure 7. This therefore leads to distinct resistivity maxima in 3D, but *not* in 2D [53], where monotonic *T*-dependence should be found on both sides of the transition, albeit with opposite slope. Physically, this difference reflects the proposed importance of "infrared" (IR, long distance) effects due to gauge fields, which should have strong (spatial) dimensionality dependence. Concerning quantum critical scaling, it is interesting that this theory proposes the emergence of *two* crossover temperature scales, both of which vanish at the transition.

**Figure 7.** Predictions of the spinon theory (reprinted with permission from Ref. [59] Copyright 2009 American Physical Society). (**a**) The phase diagram features a quantum critical point at *T* = 0, and two distinct finite-*T* crossover scales *T*∗ (above which the system is quantum critical) and *T*∗∗ (below which the system is either a metal or a gapless spin liquid). (**b**) and (**c**) Resistivity and conductivity along the lines A, B and C in the phase diagram in (**a**). Critical resistivity is predicted to diverge as *ρc*(*T*) ∼ 1/*t* in *d* = 3, leading to resistivity maxima on the metallic side (coductivity minima). In contrast, the same theory predicts finite critical resistivity *ρc*(*T*) ∼ *ρ*<sup>∗</sup> in *d* = 2 [53], hence monotonic behavior on both sides of the transition and no resistivity maxima.

We finally mention that a similar spin-charge separation theory has been very recently proposed to also describe the Wigner–Mott transition in TMD bilayers, where a possible role of charge fluctuations has also been discussed for the metallic side [60,61].

#### *4.2. Dynamical Mean Field Theory Picture of the Mott Point*

Our second theory of interest is *Dynamical Mean Field Theory* (DMFT), which explicitly ignores all nonlocal (spin or charge) spatial correlations, and therefore aims to self-consistently calculate the local electronic self-energy Σ(*ω*) [54]. Physically, its real part describes the modifications of the electronic spectra, while its imaginary part encodes the frequency and temperature dependence of the electron-electron scattering rate. In this way, this theory is not limited to low-temperature excitations only, but is able to capture strong inelastic scattering at high temperatures, and therefore describe both the (coherent) Fermi liquid regime, and also the incoherent high-temperature transport, for example the famed bad metal behavior [48,62] above the MIR limit [49].

While there exist some limiting cases where an analytic solution is possible, it is mainly a numerical approach at finite temperature. DMFT is exact in the limit of large coordination, which physically corresponds to maximal magnetic frustration. Therefore, in the simplest implementation, DMFT describes Mott physics in absence of any magnetic order, nor does it include any (inter-site) spin liquid correlations. We should mention that extensions of DMFT have recently been proposed [63] that include spinon effects, based on an alternative (matrix M, N) rotor representation. This theory, which includes some dynamical effects even at the saddle-point level, suggest that coherent spinon excitations are very fragile to charge fluctuations emerging upon the closing of the Mott gap, suppressing the spin liquid correlations not only on the metallic side, but also within the critical region. We will not further discuss these most sophisticated approaches here, but will limit our attention to the predictions of the simplest single-site DMFT theory.

When applied to the single-band Hubbard model on a frustrated lattice (such as the triangular lattice), DMFT predicts that on the metallic side the quasiparticle mass diverges linearly *m*<sup>∗</sup> ∼ |*U* − *Uc*2| <sup>−</sup><sup>1</sup> at a critical value *Uc*<sup>2</sup> (similar to the prediction of the Brinkman–Rice (BR) theory of the Mott transition [64]). The quasiparticle weight *Z* is inversely proportional to (*m*∗). Similarly, other features of the Fermi liquid such as the Kadowaki-Woods law *<sup>A</sup>* ∼ (*m*∗)<sup>2</sup> are upheld. This Fermi liquid behavior persists up to a temperature *TFL* that vanishes linearly when approaching *Uc*2. Interesting, at *Tmax* ∼ *TFL* the resistivity exhibits a maximum [29]. On the insulating side, there exist no well-defined quasiparticles as the self-energy diverges, Σ(*ω*) ∼ 1/*ω*. The electronic spectrum is split into an upper and lower Hubbard band, separated by a gap that remains nonzero at *Uc*2. The insulating state becomes unstable at a lower value of the interaction *Uc*<sup>1</sup> <sup>&</sup>lt; *Uc*2, where the gap closes *<sup>δ</sup>* ∼ |*<sup>U</sup>* <sup>−</sup> *Uc*1|, *<sup>ν</sup>z*, *<sup>ν</sup><sup>z</sup>* <sup>≈</sup> 0.8 [65]. As a result, there emerges a low-*T* first-order metal-insulator transition, and an associated phase coexistence region at *Uc*<sup>1</sup> < *U* < *Uc*<sup>2</sup> [54]. These main predictions are summarized in Figure 8.

At nonzero temperature, the first-order transition line ends at a critical point at a temperature *Tc* ≈ 0.015*W*, significantly smaller than the bare bandwidth *W*. At temperatures *T Tc* the results can be viewed as effectively quantum critical [46,65]. This quantum critical regime is centered around the so-called quantum Widom line (QWL) [43], which physically represents a finite-temperature instability trajectory of the insulating phase, as shown in Figure 8a. It extends the first-order line past *T* = *Tc*, and can experimentally be detected from an inflection point analysis [25] of the resistivity curves. In this regime, the resistivity satisfies the scaling law *ρ*(*T*, *δU*) = *ρc*(*T*)*f*(*T*/*T*0(*δU*)) with a crossover temperature scale *To* ∼ |*δU*| *<sup>ν</sup><sup>z</sup>* where *<sup>ν</sup><sup>z</sup>* <sup>≈</sup> 0.6, see Figure 8c. The crossover scale *To* is a property of the quantum critical regime and should not be confused with the low-temperature scaling of the Fermi liquid temperature *TFL*. DMFT therefore predicts two different regimes of scaling: the quantum critical regime at *T Tc*, and the metal regime dominated by scaling in *TFL*.

**Figure 8.** Predictions of DMFT theory. (**a**) Phase diagram featuring a phase coexistence region at *T* < *Tc*, and a Quantum Critical region centered around the Quantum Widom Line (QWL) (adapted with permission from Ref. [46] Copyright 2011 American Physical Society). (**b**) Resistivity (normalized by the Mott-Ioffe-Regel (MIR) limit) as a function of temperature *T* across the transition. Note the pronounced resistivity maxima on the metallic side (adapted with permission from Ref. [66] Copyright 2010 American Physical Society. (**c**) scaling collapse of the resistivity curves, displaying pronounced "mirror symmetry" of the two branches (adapted with permission from Ref. [46] Copyright 2011 American Physical Society).

#### *4.3. Percolative Phase Coexistence Picture*

In the early theories of both Mott insulators [2] and Wigner crystals, the transition from insulator to metal was often assumed to be robustly first order, at least at sufficiently low temperature. However, even the presence of weak disorder or medium-ranged interactions will create an "emulsion" (microscopic phase coexistence) with Mott/Wigner insulating "islands" in between metallic "rivers" as proposed by Spivak and Kivelson in the context of 2DEG systems in semiconductors [55,67]. If so, then tuning bandwidth and/or temperature should produce a continuous variation of the metallic fraction *x*. As long as it exceeds the percolation threshold (*x* > *xc*), the system is conducting. At *x* < *xc* the metallic domains no longer connect across the system, and conduction stops, at least at *T* = 0. Critical behavior now arises because we dealing with a classical percolation transition.

In this picture, the *T* = 0 metal-insulator transition may occur without actually closing the insulating gap Δ. Similarly, at the percolation threshold *x* = *xc* the metallic Fermi liquid at low *T* is still stable, and consequently there is no strict divergence of the effective mass *m*∗, nor the vanishing of *TFL* at the percolation threshold. We note, however, that such a (classical) percolation picture should apply only if the characteristic domain size is sufficiently larger than the characteristic correlation (or dephasing) length, therefore strictly speaking not at the lowest temperatures. However, finite-temperature variation of the transport properties should be adequately captured, as abundantly documented [68] in other systems featuring microscopic phase separation, such as Colossal Magneto-Resistance (CMR) manganites, for example.

There is, however, an interesting and nontrivial feature of the percolation picture pertaining to finite temperatures. Because of its localized spins, the entropy of the Mott or Wigner–Mott insulating phase should be higher than that of the of the metal (Fermi liquid). As a result, when raising the temperature in the metallic regime, the insulating volume fraction will *increase*: a manifestation of the Pomeranchuk effect [45]. As a result, the resistivity should increase up to some *T* = *Tmax*, above which the metallic domains no longer connect, and resistivity will decrease again [56], leading to resistivity maxima. This qualitative picture has been advocated in the work of Spivak and Kivelson, but no concrete prediction of the precise temperature dependence of the resistivity has been made, nor how the corresponding family of curves should scale as the transition is approached.

Interestingly, the same physical picture should in fact apply not only to Wigner–Mott transitions in semiconductors, but also the to conventional Mott transition, provided that there exists a well-defined metal-insulator phase coexistence region around the Mott point. Indeed, recent work on Mott organics [20] revealed precisely such a phase coexistence region, albeit only at very low temperatures of the order of at most few percent of the (bare) Fermi energy. Here, careful theoretical modeling [20] firmly established the validity of the percolation picture, but only within a well-defined phase coexistence region. In contrast, in all the systems studied (2DEG, Mott organics, moiré), the pronounced resistivity maxima persist even much further onto the metallic side, where *Tmax* can reach a substantial fraction of *TF*, where phase coexistence is very unlikely. Furthermore, recent experimental work on Mott organics by Kanoda and collaborators demonstrated [69,70] the extreme fragility of such a phase coexistence region to disorder, as generally expected in 2D systems [71]. Nevertheless, it is extremely useful to have an independent experimental method to distinguish the phase coexistence region (where percolative effects are likely) from the regimes where a more uniform electron fluid/solid resides. The possibility to do so was spectacularly demonstrated in the context of Mott organics. In the next section we discuss how the dielectric response can tell which mechanism (quasiparticle destruction or percolation) is at play in a given regime.

We briefly mention that percolation effects have been also discussed in the context of spinon theory in a recent paper [72], which does require however significant disorder. On the other hand, the Spivak-Kivelson theory does not require disorder as the micro-emulsion of insulators and metals can be self-generated. This seems to be more in line with the experiments of Section 3: at least the 2DEGs [19] and the Mott organics [38] are displaying Mott criticality in the cleanest samples possible. It is therefore very plausible that most universal features observed in all critical Mott systems are not the result of disorder, but are instead the inherent manifestations of strong correlation physics.

#### **5. Interpreting Resistivity Maxima**

As we have seen from our brief theory overview above, several scenarios were proposed, with sometimes similar predictions for characteristic features seen in experiments. A notable example is the clear emergence of the resistivity maxima on the metallic side, at a temperature *T* = *Tmax* - *TFL*, which is seen to decrease towards the transition. What is its physical content? The three theoretical pictures propose very different physical perspectives on what goes on here.

As we mentioned in Section 4.1, spinon theory [53] predicts the presence of resistivity maxima only in *d* = 3, but not in *d* = 2. However, robust resistivity maxima are clearly seen all the material systems of Section 3. An understanding of the resistivity maxima must therefore come from either the DMFT perspective on Mott physics, or from the percolative scenario. Both mechanisms provide reasonable albeit very different routes to explain the resistivity maxima. How should one distinguish them and thus identify the precise mechanism at play in a given system? Luckily, important clues were provided by recent experiments on Mott organics [20]. Here one finds two distinct regimes, both featuring similar resistivity maxima, but with very different dielectric response. One such regime is corresponds to the (spatially inhomogeneous) metal-insulator phase coexistence region, where colossal enhancement of dielectric response has been found. The other regime was found further on the metallic side, where a dramatic drop and a change of sign the dielectric constant signaled thermal destruction of coherent quasiparticles due to strong correlation effects.

In the following we show how general scaling arguments can be used within each of the two proposed scenarios, to demonstrate the general robustness of these trends, thus providing a new window of what precisely goes on near the Mott point.

#### *5.1. Resistivity Maxima from Thermally Destroying Coherent Quasiparticles*

Both experiments and theory provide evidence that a strongly correlated Fermi liquid forms on the metallic side of the Mott point, with a characteristic "Brinkman–Rice" (BR) energy scale *TBR* ∼ 1/*m*∗, which decreases towards the transition, thus characterizing the heavy quasiparticles. Inelastic electron-electron scattering increases with temperature, eventually leading to the thermal destruction of the quasiparticles around *T* ∼ *TBR*, and the associated modification of both the single particle (ARPES) spectra and the optical conductivity. At higher temperatures, transport assumes incoherent character, which can no longer be understood in terms of the quasiparticle picture or Fermi Liquid ideas alone. Its precise form generally depends on band filling and the correlation strength, but more precise predictions require a specific microscopic model and a theoretical picture.

Concrete and quantitative results, in this regime, were given by DMFT theory, which provided first insight into the origin of the resistivity maxima in certain Mott organic materials at half-filling, as well as in certain oxides. Subsequent DMFT studies stressed that the characteristic temperature scale for the resistivity maxima indeed tracks the BR scale of the quasiparticles (*Tmax* = *TBR*), while preserving the functional form of the resistivity curves across this coherence-incoherence crossover. This revealed the scaling behavior of the resistivity curves in the correlated metallic regime, with a universal scaling function of *T*/*Tmax*. The predicted scaling behavior has been confirmed by a number of experiments on various systems [29,30,32], displaying even quantitative agreement with the theoretical scaling function, with no adjustable parameters.

Further optical and dielectric studies [20] in Mott organics also confirmed the predicted destruction of the Drude peak around *TBR*, again signaling the thermal destruction of quasiparticles. They established that it dramatically affects not only DC transport, but also the dielectric response, which in this metallic regime is seen to display a dramatic drop from moderate positive values at *T* > *Tdrop* ∼ *Tmax* to very large but negative values at *T* < *Tdrop* ∼ *Tmax*. These studies, combining experiments and DMFT theory, have firmly established that the dielectric response can be used to directly reveal the thermal destruction of quasiparticles around the BR temperature.

In the following, we extend the systematic studies of Ref. [29], to stress that within DMFT both DC transport and the dielectric response display the characteristic crossover behavior across *TBR*, and the associated scaling behavior upon approaching the Mott point. To do this we calculate the dielectric function <sup>1</sup> as a function of temperature and interaction *U*, using the same setup as in our recent work [20]. For simplicity, we focus on a simple semi-circular band model at half filling, and carry out DMFT calculations using the standard CTQMC impurity solver with the the Maximum Entropy method for analytical continuation to the real axis. Just as in Ref. [29], once we get the single particle self energy from our DMFT equations, we calculate the real part of optical conductivity *σ*1(*ω*) from the standard Kubo formula, and the imaginary part of optical conductivity *σ*2(*ω*) using the Kramers-Kroning tranform; the (complex) dielectric function is then obtained via [73]

$$
\epsilon(\omega) = 1 + 4\pi i \frac{\sigma(\omega)}{\omega}.\tag{3}
$$

The results for the single-particle density of states and the optical conductivity are shown in Figure 9, and results for the DC transport and the low-frequency dielectric response are displayed in Figure 10, for the parameter range corresponding to the correlated metallic phase (*U Uc*1). Here panels (a) and (b) reproduce the results of Ref. [29], showing the characteristic scaling behavior of the resistivity maxima near the Mott point. The analogous behavior for the dielectric function <sup>1</sup> is shown in panels (c) and (d), firmly establishing that the observed crossover behavior assumes a universal scaling form in the correlated metallic regime. The notion that the thermal destruction of quasiparticles lies behind both phenomena is seen even more clearly in Figure 11, where we show how the scale *Tmax* for the resistivity maxima, and the scale *Tdrop* of the dielectric response, both scale with the quasiparticles weight *Z* = *m*/*m*∗, as the transition is approached.

These results establish a way to experimentally recognize the thermal destruction of quasiparticles as a dominant mechanism behind the resistivity maxima within a correlated but uniform metallic phase. Since the correlation processes captured by DMFT are essentially *local* (i.e. "ultraviolet, UV"), these effects should not display significant dependency on spatial dimensional. Indeed, experiments have shown that similar resistivity maxima are seen within correlated metallic phases both in 2D and in 3D systems.

**Figure 9.** (**a**) DMFT results for the evolution of the single-particle Density of States (DOS) for several values of the temperature (reprinted with permission from Ref. [66] Copyright 2010 American Physical Society), as well as (**b**) that of the optical conductivity, in the strongly correlated metallic regime. Different colors correspond to the four distinctive transport regimes (inset in (**b**)). DOS features a distinct quasiparticle peak at low temperatures, which is thermally destroyed at temperature *Tmax* <sup>=</sup> *TBR* <sup>∼</sup> (*m*∗)−1, where the resistivity (inset of right panel) reaches a maximum. The optical conductivity displays the corresponding suppression of the low-frequency Drude peak around the same temperature.

**Figure 10.** (**a**) DC resistivity as a function of temperature for several interaction strengths. (**b**) Scaled resistivity curves. (**c**) Real part of dielectric function <sup>1</sup> at *ω*/*D* = 0.01, as a function of temperature for several interaction strengths. (**d**) Scaled dielectric function curves. Results are obtained for a half-filled Hubbard model solved within DMFT.

**Figure 11.** (**a**) *Tdrop* as a function of *Tmax*. (**b**) *Tmax* as a function of *Z*. (**c**) *Tdrop* as a function of *Z*.

Note, however, that DMFT predicts very different behavior closer to the Mott point, specifically within the phase coexistence region at *T* < *Tc* and *Uc*<sup>1</sup> < *U* < *Uc*2. Here, just as around any first-order phase transition line, we expect hysteresis phenomena and inhomogeneous phase separation, where metallic and insulating domains coexist on a nano-scale. As stressed in the seminal work by Spivak and Kivelson [55], thermal effects can modify the relative volume fraction of the two coexisting phases, producing under appropriate conditions the characteristic resistivity maxima. In recent work motivated by experiments, a microscopic "hybrid-DMFT" approach was developed [20] to quantitatively describe this regime in the context of Mott organics, resulting in spectacular agreement with experiments. In the following section, however, we wish to stress that the *qualitative* aspects of this regime display a number of universal scaling features, which can be precisely understood from the perspective of percolation theory.

#### *5.2. Percolation Scenario Due to Phase Coexistence*

To focus on the universal scaling aspects of percolative processes within the metalinsulator phase coexistence region, we follow the seminal ideas of Efros and Shklovskii [74], and set up a two-component random resistor network model, with characteristic lowfrequency form for the (complex) conductivity for each component:

$$\begin{aligned} \sigma\_I &= \sigma\_I^o \exp(-\Delta/T) - i \mathbb{C} \omega\_\prime \\ \sigma\_M &= \frac{\sigma\_M^o}{1 - i \omega \tau} . \end{aligned} \tag{4}$$

Here we assumed activated DC transport for the insulating component, with capacitance *C* and a standard Drude form for the conducting component, with finite DC conductivity *σ<sup>o</sup> <sup>M</sup>*. To leading order near the percolation point, we ignore the *T*-dependence of *σ<sup>o</sup> <sup>M</sup>*, *<sup>σ</sup>*<sup>0</sup> *<sup>I</sup>* , and *C*, since the dominant effects come from the variation of the respective volume fractions, and the activated form of insulating transport. The temperature is expressed in the units of the activation energy Δ, which is also taken to be a constant. The corresponding expressions for the (complex) dielectric functions of the two components are given by:

$$\begin{aligned} \varepsilon\_I &= 1 + 4\pi \text{C} + \frac{4\pi i}{\omega} \sigma\_I^o \exp(-\Delta/T), \\ \varepsilon\_M &= 1 - 4\pi \tau \sigma\_M^o + \frac{4\pi i}{\omega} \sigma\_M^o. \end{aligned} \tag{5}$$

Here we ignored the capacitance of the metallic domains, which can be neglected if *τσ*<sup>0</sup> *<sup>M</sup>*/*C* 1.

Such a random resistor network model is appropriate for any percolating two component metal-insulator system. To describe formation of the resistivity maxima, an additional physical condition has to be met, as emphasized by Spivak and Kivelson in the context of Wigner–Mott transitions, but which is in fact valid for any Mott-like system in general. As we mentioned before, this "Pomeranchuk effect" [45] requires that the first-order line (and the entire phase coexistence region) be "tilted" towards the metal, so that the higher-entropy phase emerges at higher temperatures. To schematically represent such a situation we assume that, within the phase coexistence region, the volume fraction *x* of the metallic component *decreases* with temperature. As an illustration, we take the following simple model:

$$\mathbf{x}(\mathbf{x}\_o, T) = \mathbf{x}\_c + \frac{1}{2} \tanh[\frac{\mathbf{x}\_o - \mathbf{x}^\*(T)}{w(T)}],\tag{6}$$

where *xc* represents the percolation threshold, *w*(*T*) = *a*(*Tc* − *T*)/*Tc* defines the width of the coexistence region, and *x*∗(*T*) = *xc* + *b*(*T*/*Tc*), as illustrated in Figure 12a, for *a* = 0.4 and *b* = 1. In this model, the parameter *xo* controls the metallic volume fraction at *T* = 0, which decreases at *T* > 0, and reaches the percolation threshold *x* = *xc* at *T* = *T*∗(*xo*) = *Tc*(*xo* − *xc*)/*b*. Physically, the the DC resistivity will first increase with *T* as the metallic volume fraction decreases. Past percolation threshold, however, the metallic domains no longer connect. Transport then assumes insulating (activated) form, resulting in subsequent resistivity decrease at *T* > *T*∗, and the emergence of resistivity maxima around *T* ∼ *T*∗. Similarly, the dielectric constant <sup>1</sup> grows (diverges) as the percolation threshold is approached from the insulating side, due to the formation of large metallic clusters with increased polarizability. On the metallic side, however, it displays a rapid decrease, dropping to large negative values within the metallic phase. As a result, dielectric response displays colossal enhancement around the percolation threshold, a phenomenon that can be viewed as a smoking gun for percolative charge dynamics.

**Figure 12.** (**a**) The red line is *x* = *x*(*T*∗). For *T* larger than the blue dashed line, *x* = 0. We calculate the percolation results along the grey dashed lines. (**b**) *R*/*R<sup>o</sup> <sup>M</sup>* as a function of *T* for different *T*∗. (**c**) Scaled resistivity curves.

To illustrate these ideas, we use the Effective Medium Approximation (EMA) for percolation, which solves the following nonlinear equations for the complex dielectric function:

$$\text{tr}(\frac{\epsilon\_M - \epsilon}{\epsilon\_M + (z/2 - 1)\epsilon}) + (1 - x)(\frac{\epsilon\_I - \epsilon}{\epsilon\_I + (z/2 - 1)\epsilon}) = 0,\tag{7}$$

and for illustration selected *σ<sup>o</sup> M*/*σ<sup>o</sup> <sup>I</sup>* = 100, *τσ<sup>o</sup> <sup>M</sup>* = 1000, *C* = 1, *Tc*/Δ = 0.4, and *z* = 4 corresponding to 2D transport (*xc* = 0.5). Precisely the anticipated behavior is observed from numerically solving the EMA equation for the corresponding DC resistivity *R* = *σ*−1, as shown in Figure 12b. Here, we select several values of *xo*, corresponding to *x* > *xc* (low temperature metallic regime), and plot the resistivity as a function of temperature (following dashed lines in Figure 12a. We observe distinct resistivity maxima around the temperature *T*∗(*xo*) corresponding to the percolation threshold. Note how the maxima become sharper and sharper as *T*∗ is reduced, corresponding to the exponential (activated) decrease of the "field" *h* ∼ exp{−Δ/*T*∗}. The expected behavior is also seen in dielectric response, as shown in Figure 13a, where we observe sharp maxima at *T* ∼ *T*∗. Here again we see the increased "rounding" of these maxima at higher *T*∗, corresponding to larger *h*(*T*∗). This behavior can be seen even more clearly in Figure 13b, where <sup>1</sup> is plotted as a function of the reduced concentration *τ* = (*x*(*T*) − *xc*)/*xc*, which vanishes at *T* = *T*∗.

**Figure 13.** (**a**) The dielectric constant <sup>1</sup> as a function of *T* for different *T*∗. (**b**) <sup>1</sup> as a function of *τ* for different *T*∗. (**c**) Scaled dielectric function curves.

These qualitative trend can be even more rigorously described within the scaling theory for percolation, where the DC conductivity as well as the *ω* = 0 dielectric constant are known to satisfy the following scaling relations:

$$
\sigma\_1(\tau, h) = \sigma\_M^o h^s F\_\sigma(\tau/h^m); \quad \epsilon\_1(\tau, h) = h^{s-1} F\_\sigma(\tau/h^m), \tag{8}
$$

where *τ* = (*x*(*T*) − *xc*)/*xc* measures the distance to the percolation threshold, and *h* = *σI*/*σ<sup>M</sup>* plays a role of the "symmetry breaking field", which leads to the rounding of the transition. The critical exponents *s* and *m*, as well as the crossover scaling

functions *F<sup>σ</sup>* and *F* are universal quantities within percolation theory. To illustrate this scaling behavior within EMA, we collapse the family of resistivity curves by plotting *Rh<sup>s</sup>* as a function of *τ*/*hm*, as shown in Figure 12c, and 1*h*1−*<sup>s</sup>* as a function of s a function of *τ*/*hm*, as shown in Figure 13c. Note how a perfect scaling collapse is observed here, but only around the peak of the dielectric response, i.e. only close to the percolation threshold. Such behavior is, in fact, not surprising, since we expect scaling phenomena to arise only within a given critical region, and not further away from the critical point.

We should stress again that all our qualitative results are rigorously valid within general percolation theory for our two-component phase coexistence model, and EMA was simply used as an illustration. EMA correctly captures the general crossover phenomena associated with percolation, but only introduces approximate values for the critical exponents *sEMA* = 0.5 and *mEMA* = 0.5, which are otherwise know even more accurately from numerical simulations. These details, however, are not of direct relevance for our purposes. What is important is the result that, within our "Pomeranchuk" model for phase coexistence, the percolation scenario predicts distinct resistivity maxima but also striking colossal dielectric anomalies, at the same temperature scale of *T* = *T*∗ which decreases towards the MIT. This behavior is in distinct contrast to the behavior we found from the DMFT picture of a correlated but uniform metallic phase, which also leads to resistivity maxima, but very different behavior of the dielectric response. This observation, which was quantitatively validated in recent experiments on Mott organics [20], thus reveals a distinct criterion to settle the long-lasting controversies between the origin of the resistivity maxima in different systems.

#### **6. Conclusions**

In this paper we discussed three different classes of physical systems which all display very similar phenomenology expected for Mott-like metal-insulator transitions. We stressed that most qualitative features are clearly seen in all these examples, including the continuous decrease of the characteristic energy scales *TFL*, *TBR* = *Tmax*, *To*, Δ towards the transition, the phenomenon of quantum critical scaling seen in transport, as well as the emergence of distinct resistivity maxima on the metallic side. These observations, which are starting to portray a robust and consistent phenomenology of Mott criticality, is putting serious constraints on theory. We discussed which of these features seem compatible with various proposed theoretical pictures of the Mott point, and which ones do not.

In the final section of this paper we also presented new theoretical results, which open the possibility to precisely determine, from experiments, which mechanism dominates in which regime. We argue that the *dielectric response* offers unique insights, which so far have not been appreciated enough, as a powerful tool to distinguish between different phase coexistence and the thermal destruction of quasiparticles.

A class of issues we did not discuss in any detail in this paper is the (explicit) role of disorder around the Mott point. Given the fact that new classes of ultra-clean material are starting to emerge, with even more pronounced salient features of Mott criticality, it is becoming possible to plausibly minimize the role of disorder on experimental grounds. On the other hand, new experimental efforts are starting [69,70] to emerge in the opposite direction: to systematically add and to control the level of disorder, for example by highenergy X-ray irradiation. These fascinating research directions are guaranteed to open entirely new chapters in the study of metal-insulator quantum criticality. This will require theorists to rekindle the efforts to understand the interplay of strong correlation with disorder [75], and perhaps to develop new ideas in the process.

**Author Contributions:** Y.T. carried out theoretical calculations and performed the analyses. V.D. and L.R. designed the project. All authors discussed the data, interpreted the results, and wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the NSF Grant No. 1822258, and the National High Magnetic Field Laboratory through the NSF Cooperative Agreement No. 1644779 and the State of

Florida. LR acknowledges support by the Swiss National Science Foundation via Ambizione grant PZ00P2\_174208.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are provided: https://github.com/yutingtanphysics/ (accessed on 3 May 2022) How-to-Recognize-the-Universal-Aspects-of-Mott-Criticality.

**Acknowledgments:** We thank Pak Ki Henry Tsang for providing technical support in doing computer cluster calculations.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **Neural Network Solver for Small Quantum Clusters**

**Nicholas Walker 1, Samuel Kellar 1, Yi Zhang 1,2, Ka-Ming Tam 1,3,\* and Juana Moreno 1,3,\***


**Abstract:** Machine learning approaches have recently been applied to the study of various problems in physics. Most of these studies are focused on interpreting the data generated by conventional numerical methods or the data on an existing experimental database. An interesting question is whether it is possible to use a machine learning approach, in particular a neural network, for solving the many-body problem. In this paper, we present a neural network solver for the single impurity Anderson model, the paradigm of an interacting quantum problem in small clusters. We demonstrate that the neural-network-based solver provides quantitative accurate results for the spectral function as compared to the exact diagonalization method. This opens the possibility of utilizing the neural network approach as an impurity solver for other many-body numerical approaches, such as the dynamical mean field theory.

**Keywords:** metal insulator transition; anderson localization; random disorder; typical medium theory; dynamical mean field theory; coherent potential approximation; neural network; quantum impurity solver; Anderson impurity

#### **1. Introduction**

A single quantum impurity is the simplest quantum many-body problem for which interaction plays a crucial role [1,2]. It was designed as a model to describe diluted magnetic impurities in a non-magnetic metallic host. In the 1960s, it was showed that the perturbation series in coupling strength diverges even with an infinitesimal anti-ferromagnetic coupling value [1]. This early unexpected result motivated the development of innovative approaches to model the strongly correlated systems [1–3].

While the physics of a single impurity problem has been rather well studied, interest in the quantum impurity problem was revived during the 1990s. This was partly due to the interest in mapping lattice models onto impurity models [4–8], since at infinite dimensions, the lattice models are equivalent to single impurity models in a mean-field represented by the density of states of the host. This approximated mapping is known as the dynamical mean-field theory. It has been further generalized to cluster impurity models which include some of the effects present in finite dimensional systems [9–11].

These mappings provide a systematic tractable approximation for the lattice models and have become a major approach in the field of strongly correlated systems [11]. Combined with density functional theory, they provide one of the best available methods for the study of the properties of materials in which strong interactions are important [12].

The mean-field where the single impurity is embedded, i.e., the density of the bath, can be rather complicated. Therefore, there is in general no analytic method for an accurate solution. Many different methods for solving the effective impurity problem have been proposed. They can be broadly divided in two categories: semi-analytic and numeric.

Between the semi-analytic methods, the most widely used one is the iterative perturbation theory. It interpolates the self-energy at both the weak and strong coupling limits and incorporates some exact constraints, such as the Luttinger theorem [13]. Another

**Citation:** Walker, N.; Kellar, S.; Zhang, Y.; Tam, K.-M.; Moreno, J. Neural Network Solver for Small Quantum Clusters. *Crystals* **2022**, *12*, 1269. https://doi.org/10.3390/ cryst12091269

Academic Editor: Andrej Pustogow

Received: 2 August 2022 Accepted: 29 August 2022 Published: 6 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

widely used semi-analytic method is the local moment approximation. It considers the perturbation on top of the strong coupling limit represented by the unrestricted Hatree Fock solution [14].

Numerical methods can also be divided into two main classes, diagonalization-based and Monte-Carlo-based. Diagonalization methods usually require discretizing the host by a finite number of so-called bath sites. The Hamiltonian which includes the bath sites and the impurity site are diagonalized exactly [15]. Another digonalization-based method is the numerical renormalization group in which the bath sites are mapped onto a one-dimensional chain of sites. The hopping amplitude decreases rapidly down the chain. The model is then diagonalized iteratively as more sites are included [16]. Density matrix renormalization group and coupled cluster theory have also been used as impurity solvers [17–20].

Hirsh and Fye were the first to propose the quantum Monte Carlo method for solving impurity problems [21]. The time axis of the simulations is divided up using the Trotter– Suzuki approximation. The interaction in each time segment is handled by the Hubbard– Stratonovich approximation [22]. The Monte Carlo method is then used to sample the Hubbard–Stratonovich fields. Continuous time quantum Monte Carlo methods have seen a lot of development over the last decade. They directly sample the partition function without using the Trotter–Suzuki approximation [23], similar to the Stochastic Series Expansion in the simulation of quantum spin models. Notably, the continuous time quantum Monte Carlo using the expansion with respect to the strong coupling limit has been proposed and complicated coupling functions beyond simple Hubbard local density-density coupling terms can now be studied [24].

On the other hand, the past few years have seen tremendous development of machine learning (ML), both algorithms and their implementation [25–27]. Many of the ML approaches in physics are designed to detect phase transitions or accelerate Monte Carlo simulations. It is a tantalizing proposal to utilize ML approaches to build a solver for quantum systems.

To build a quantum solver based on the ML approach, we need to identify the feature vector (input data) and the label (output data) for the problem. Then a large pool of data must be generated to train the model, specifically a neural network model. The Anderson impurity problem is a good test bed for the validity of such a solver. We note that similar ideas have been explored using machine learning approaches [28]. This paper is focused on using the kernel polynomial expansion and supervised ML, specifically a neural network, as the building blocks for a quantum impurity solver.

While it is relatively inexpensive to solve a single impurity problem using the above methods in modern computing facilities, current interest in the effects of disorder warrants the new requirement of solving a large set of single or few impurities problems in order to calculate the random disorder average [29–35]. Our hope is that a fast neural-networkbased numerical solver in real frequency can expand the range of applicability of the recently developed typical medium theory to interacting strongly correlated systems, such as the Anderson–Hubbard model [36–39].

This paper is organized as follows. In the next section, we map the continuous Green function into a finite cluster as has been performed in many dynamical mean field theory calculations. In Section 3, we discuss the expansion of the spectral function in terms of Chebyshev polynomials. In Section 4, we discuss how to use the results from Sections 2 and 3 as the feature vectors and labels of the neural network. In Section 5, we present the spectral function calculated by the neural network approach. We conclude and discuss future work in the last section.

#### **2. Representing the Host by a Finite Number of Bath Sites**

We first identify the input and the output data of a single impurity Anderson model. The input data includes the bare density of states, the chemical potential and the Hubbard interaction on the impurity site. For a system in the thermodynamic limit, the density

of states is represented by a continuous function. Since inputting a continuous function to the neural network presents a problem, we describe the continuous bath by a finite number of poles as it is performed within the exact diagonalization approach [15,40–43]. We approximate the host Green function by a cluster of bath sites,

$$\mathcal{G}\_0(i\omega\_n) \approx \mathcal{G}\_0^{cl}(i\omega\_n). \tag{1}$$

In the exact diagonalization method, the continuum bath is discretized and represented by a finite number of so-called bath sites, see Figure 1. Assuming that there are *Nb* bath sites, each bath site is characterized by a local energy (*i*) and a hopping (*ti*) term to the impurity site. Two additional variables, one the local Hubbard interaction (*U*) and the other the chemical potential (*f*), are required to describe the impurity site. Therefore, there is a total of 2 + 2*Nb* variables representing the impurity problem.

The full Hamiltonian in the discretized form (Figure 1) is given as

$$\begin{aligned} H &= \sum\_{i,\sigma} t\_i (c\_{i,\sigma}^\dagger c\_{0,\sigma} + H.c.) + \sum\_{i,\sigma} \epsilon\_i c\_{i,\sigma}^\dagger c\_{i,\sigma} + \\ \mathcal{U}(c\_{0,\uparrow}^\dagger c\_{0,\uparrow} - 1/2)(c\_{0,\downarrow}^\dagger c\_{0,\downarrow} - 1/2) - \epsilon\_f \sum\_{\sigma} c\_{0,\sigma}^\dagger c\_{0,\sigma} \end{aligned} \tag{2}$$

*c*† *<sup>i</sup>*,*<sup>σ</sup>* and *ci*,*<sup>σ</sup>* are the creation and annihilation operators for site *i* with spin *σ*, respectively. The impurity site is denoted as the 0-th site. The sum of the bath sites goes from 1 to *Nb*, and the spin *σ* = ±1/2.

The host Green function represented in such a finite cluster can be written exactly as,

$$\mathcal{G}\_0^{cl}(i\omega\_\hbar) = (i\omega\_\hbar + \varepsilon\_f - \sum\_{k=1}^{N\_b} \frac{t\_k t\_k^\*}{i\omega\_\hbar - \varepsilon\_k})^{-1}.\tag{3}$$

Many different prescriptions for the parameterization of the host Green function have been investigated in detail for exact diagonalization solvers [44]. Conceptually, practical applications of the numerical renormalization group method also require the approximated mapping of the problem onto a finite cluster chain. Unlike the exact diagonalization method, the cluster chain can be rather large; therefore much higher accuracy can be attained in general.

These approaches do not mimic the continuum in the time dimension as it is done by continuous time Quantum Monte Carlo methods, and the mapping onto the finite cluster may represent a nuisance. Nonetheless, this is a necessity for any diagonalization-based method. In our case, the mapping presents an opportunity to naturally adapt the method to a machine learning approach in which a finite discretized set of variables is required.

Under the above approximation, the finite set of variables, {*ti*}, {*i*}, *U*, *<sup>f</sup>* , can be treated as the input feature vector for the machine learning algorithm. The next question is what is the desired output or label. We will focus on the spectral function in this study. For this purpose, the next step is to represent the spectral function with a finite number of variables instead of a continuous function. The kernel polynomial method fulfills this goal [45].

#### **3. Expanding the Impurity Green Function by Chebyshev Polynomials**

In this section, we briefly discuss the kernel polynomial method for the calculation of the spectral function of a quantum interacting model. Once the host parameters, the impurity interaction and chemical potential are fixed, the ground state of the cluster is obtained using a diagonalization method for sparse matrices. We use the Lanczos approach in the present study [46,47]. Once the ground state is found, the spectral function can be calculated by applying the resolvent operator, 1/(*ω* − *H*), to the ground state. A popular method is the continuous fraction expansion [47]. The challenge is that the continuous fraction tends to be under-damped and produce spurious peaks [46,47]. A more recent method is to use an orthogonal polynomial expansion. We will argue that for the application of the ML method, the polynomial expansion method tends to produce better results as we will explain later [45].

The zero temperature single particle retarded Green function corresponding to a generic many-body Hamiltonian is defined as

$$G(\omega) = \langle GS|c \frac{1}{\omega + i0^{+} - H} c^{\dagger}|GS\rangle,\tag{4}$$

where |*GS* is the ground state of *H*, and *c*† and *c* are the creation and annihilation operators, respectively [45,48]. The spectral function is given as *A*(*ω*) = −(1/*π*)*Im*(*G*(*ω*)).

Consider the Chebyshev polynomials of the first kind defined as *Tn*(*x*) = *cos*(*n arccos*(*x*)). Two important properties of the Chebyshev polynomial are their orthogonality and their recurrence relation. The product of two Chebyshev polynomials integrated over *x* = [−1, 1] and weighted by the function *wn* = <sup>1</sup> *π* √ <sup>1</sup>−*x*<sup>2</sup> is given as

$$\int d\mathbf{x} w\_n(\mathbf{x}) T\_n(\mathbf{x}) T\_m(\mathbf{x}) = \frac{1 + \delta\_{n,0}}{2} \delta\_{n,m}.\tag{5}$$

The recurrence relation is given as

$$T\_n(\mathbf{x}) = 2\mathbf{x}T\_{n-1}(\mathbf{x}) - T\_{n-2}(\mathbf{x}).\tag{6}$$

The Chebyshev polynomials expansion method is based on the fact that the set of Chebyshev polynomials form an orthonormal basis as defined in Equation (5). Thus a function, *f*(*x*), defined within the range of *x* = [−1, 1] can be expanded as

$$f(\mathbf{x}) = \mu\_0 + 2\sum\_{n=1}^{\infty} \mu\_n(\mathbf{x}) T\_n(\mathbf{x}),\tag{7}$$

and the expansion coefficient can be obtained by the inner product of the function *f*(*x*) and the Chebyshev polynomials as follow

$$
\mu\_{\text{fl}} = \int\_{-1}^{1} d\mathbf{x} f(\mathbf{x}) T\_{\text{fl}}(\mathbf{x}) w\_{\text{fl}}(\mathbf{x}). \tag{8}
$$

Practical calculations involve truncation at a finite order. The truncation is found to be problematic when the function to expand is not smooth. In our problem, the function is the spectral function of a finite size cluster, which is a linear combination of a set of delta functions. For this reason, a direct application of the above formula will not provide a smooth function. This is analogue to the Gibbs oscillations in the Fourier expansion. The remedy is to introduce a damping factor (kernel) in each coefficient of the expansion [45,49–52]. Here we use the Jackson kernel given as

$$f(\mathbf{x}) \approx \sum\_{n=0}^{N} g\_n \mu\_n(\mathbf{x}) T\_n(\mathbf{x}), \text{where} \tag{9}$$

$$\mathbf{g}\_{\mathbb{H}} = \frac{(N - n + 1)\cos(\frac{\pi n}{N + 1}) + \sin(\frac{\pi n}{N + 1})\cot(\frac{\pi}{N + 1})}{N + 1}. \tag{10}$$

We refer the choice of the damping factor to the review by Weiße et al. [45]. We list the steps for calculating the coefficients as follows.


$$
\mu\_n = \preccurlyeq \alpha\_0 | \alpha\_n > ,\tag{11}
$$

where |*α*<sup>0</sup> = *<sup>c</sup>*†|*GS* > and |*α<sup>n</sup>* = *Tn*(*H*)|*α*<sup>0</sup> . With the |*α*<sup>0</sup> and the |*α*<sup>1</sup> = *H*|*α*<sup>0</sup> ready, all the higher order coefficients can be obtained via the recurrence relation

$$|a\_n\rangle = 2H|a\_{n-1}\rangle - |a\_{n-2}\rangle.\tag{12}$$

#### 5. The spectral function is obtained by feeding the coefficients into Equation (7).

All the coefficients can be obtained by repeated use of Equation (12) which involves matrix vector multiplication. The matrix for an interacting system is usually very sparse, and the computational complexity of the matrix vector multiplication is linear with respect to the vector length, which grows as 4*Nb*+<sup>1</sup> assuming no reduction by symmetry is applied.

Since practical calculations are limited to a finite order, *N*, the impurity Green function can be represented by *N* coefficients of the Chebyshev polynomials expansion. It is worthwhile to mention that the expansion of the Green function in terms of orthogonal polynomials is independent of the method for obtaining the coefficients of the polynomials. Instead of employing exact diagonalization, a recent study shows that one can obtain the expansion coefficients by representing the quantum states by matrix products [48].

#### **4. Feature Vectors and the Labels for the Machine Learning Algorithm**

Our strategy is to train a neural network with a large set of variables for the host, i.e., the bath sites, the impurity interaction and the impurity chemical potential. The impurity solver calculates the impurity Green function for a given bath Green function, local impurity interaction and chemical potential that is a total of 2 + 2*Nb* variables for the input.

The output is the impurity Green function which can be represented by *N* coefficients of the Chebyshev polynomials expansion. Using the above method, the spectral function is effectively represented in terms of *N* coefficients. It allows us to naturally employ the supervised learning method by identifying the 2 + 2*Nb* variables as the input feature vectors and the *N* variables as the output labels.

While the kernel polynomial method grows exponentially with the number of sites, the end result is represented by a finite number of coefficients which presumably does not scale exponentially with the number of sites. Once the neural network is properly trained, we can use it to predict the impurity Green function without involving a calculation which scales exponentially.

#### **5. Results**

We generated 5000 samples for training the neural networks with randomly chosen parameters which are drawn uniformly from the range listed as follows:

$$\begin{aligned} t\_{i,\uparrow} &= t\_{i,\downarrow} = [0, 1.5], \\ \epsilon\_{i,\uparrow} &= \epsilon\_{i,\downarrow} = [-5, 5], \\ \mathcal{U} &= [0, 10], \\ \epsilon\_f &= [-2.5, 2.5]. \end{aligned} \tag{13}$$

We assume that the electron bath has a symmetric density of states. That is, *ti* = *ti*<sup>+</sup>*Nb*/2 and *<sup>i</sup>* = −*i*+*Nb*/2 for *i* = 1 to *Nb*/2 and *Nb* even. This further reduces the number of variables in the feature vector to *Nb* + 2.

Before embarking on training the neural network, we would like to have some understanding of the range of Chebyshev coefficients. For this purpose, we randomly pick 32 samples and plot the coefficients in Figure 2. There are two prominent features of the coefficients: 1. There are clear oscillations and the coefficients do not decrease monotonically. 2. For all cases shown here, the coefficients essentially vanish for orders around 200 and higher. Due to these two reasons, we decide to train the neural network for the coefficients between 0 and 255.

**Figure 2.** *Cont.*

**Figure 2.** The coefficients of the Chebyshev polynomial expansion for 32 randomly chosen parameter sets of the finite cluster. Only the first 256 coefficients are shown, as higher order terms are vanishingly small. Only the coefficients directly calculated from the kernel polynomial method (KPM) are shown here. The coefficients obtained from the neural network match very closely with the ones from the KPM and would not be visible by laying them on the same plot, and thus they are omitted. The spectral functions and the corresponding parameters are presented in Figures 3 and 4 respectively. We demonstrate the quality of the coefficients obtained from the neural network in Figure 5. The magnitude of the the last five coefficients in each panel is smaller than 10<sup>−</sup>5.

With the above approximations, the task of solving the Anderson impurity model boils down to mapping a vector containing *Nb* + 2 variables to a vector containing *N* coefficients. For the particular case, we study we choose *Nb* = 6 and *N* = 256. Machine learning algorithms can thus be naturally applied to this mapping.

We set up an independent dense neural network for each coefficient. The neural network has 14 layers. The input layer contains *Nin* = *Nb* + 2 units, and the output layer contains the expansion coefficient for one specific order. The 12 hidden layers have the following number of units: 2*Nin*, 2*Nin*, 4*Nin*, 4*Nin*, 8*Nin*, 8*Nin*, 8*Nin*, 8*Nin*, 4*Nin*, 4*Nin*, 2*Nin*, 2*Nin*.

As we consider a total of 256 orders, we train 256 independent neural networks. Considering the coefficients at different orders separately may lose some information contained in the correlations between them. While it is possible to predict a few coefficients by only one neural network, we do not obtain a good prediction using a single neural network for all 256 coefficients without an elaborated fine tuning. Therefore, here instead of searching for a optimal number of coefficients to be predicted for one neural network, we consider each coefficient independently.

In Figure 3, we show the spectral functions corresponding to the same 32 samples used in Figure 2. Both the results from the direct numerical calculation based on the Lanczos method and recurrence relations and those predicted by the neural network are plotted. They basically overlap each other. There is a slight difference for the range of energies where the spectral function is nearly zero. This is perhaps due to the incomplete cancellation among the expansion terms at different orders. An improvement may be attainable if we consider the correlations of the coefficients for different orders. The input parameters of each of the 32 samples are plotted in Figure 4.

**Figure 3.** The spectral function, *A*(*ω*), is plotted for the same 32 randomly chosen parameter sets used in Figure 2. Both the results from the KPM and from the neural network are shown. They match each other very closely and visually overlap. A closer inspection reveals that there are slight oscillations in the spectral function when the weights are very small. This may be due to the inexact cancellations of different orders in the coefficients generated by the neural network method. In general, these oscillations are rather small and only appear when the spectral weight drops to near zero.

**Figure 4.** The input parameters of the 32 samples used in Figures 2 and 3; *i* = 1 corresponds to *U*; *i* = 2 corresponds to *<sup>f</sup>* ; *i* = 3, 4, 5 correspond to 1, 2, 3; and *i* = 6, 7, 8 correspond to *t*1, *t*2, *t*3.

Evidence of the capability of the neural network approach can be seen in Figure 5 where we plot the comparison of the first 32 expansion coefficients obtained by the direct numerical calculation and the neural network prediction. We find that two methods give very close results for the 1000 testing samples we consider. Please be reminded that 5000 samples are used for training. All the 1000 testing samples exhibit a linear trend. This clearly shows that a neural network is capable of providing a good prediction. There were no exceptional outliers among the 1000 testing samples we tested.

**Figure 5.** Comparison of the first 32 coefficients as computed by the KPM and the neural network method; 1000 samples are plotted in each figure. The figures are ordered from left to right and top to bottom from the order 0−th to the order 31−th.

#### **6. Conclusions**

We demonstrate that the supervised machine learning approach, specifically the neural network method, can be utilized as a solver for small quantum interacting clusters. This could be potentially useful for the statistical DMFT or the typical medium theory for which a large number of impurity problems have to be solved to perform the disorder averaging [29–33,36]. The main strategy is to devise a finite number of variables as the feature vector and the label for the supervised machine learning method. In line with the exact diagonalization method for the single impurity Anderson model, the feature vector is represented by the hoppings and the energies of the lattice model. The output we consider, the spectral function, is represented in terms of Chebyshev polynomials with a damping kernel. The labels are then the coefficients of the expansion. By comparing the coefficients directly calculated by the Lanczos method and the recurrence relation with the ones by the neural network, we find the agreement between these two methods is very good. Notably, among the 1000 samples we tested, there is no exceptional outlier. They all have a good agreement with the results of the numerical method.

For a simple impurity problem, the present method may not have an obvious benefit, as a rather large pool of samples have to be generated for training purposes. The situation is completely different for the study of disorder models, such as those being studied by the typical medium theory, where the present method has a clear advantage. Once the neural network is trained, the calculations are computationally inexpensive. For systems in which disorder averaging is required, this method can beat most if not all other numerical methods in term of efficiency. Moreover, the present approach is rather easy to be generalized for more complicated models, such as the few impurities models required in the dynamical cluster approximation. In addition, this method can be easily adapted to the matrix product basis proposed for the kernel polynomial expansion [48,53].

The range of parameters we choose in the present study covers the metal insulator transition of the Hubbard model within the single impurity dynamical mean field theory. An interesting question is the validity of the trained neural networks for the parameters well outside the range of the training data. We expect that the results could deteriorate without additional training data.

The size of the bath we choose in the present study is *Nb* = 6. This choice is somewhat arbitrary, for more accurate results for the single impurity problem a larger number of bath sites is preferred. The computational time required to generate the training set scales exponentially with the number of bath sites; however the computational time required to train the neural networks only scales as the third power with the number of bath sites. A larger number of bath sites should pose no problem for the present approach.

The ideas presented in this paper are rather generic. They can be generalized to the solutions by different solvers. For example, this approach can be adapted to Quantum Monte Carlo results as long as they can be represented in some kind of series expansion [54,55]. Our approach can also be adapted to predict the the coefficients of the coupled cluster theory [19,20].

**Author Contributions:** Conceptualization, K.-M.T.; Data curation, K.-M.T.; Formal analysis, J.M.; Funding acquisition, K.-M.T. and J.M.; Investigation, K.-M.T.; Methodology, K.-M.T., N.W., S.K. and Y.Z.; Project administration, J.M.; Supervision, J.M.; Writing—original draft, K.-M.T.; Writing—review & editing, J.M., N.W., S.K. and Y.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is funded by the NSF Materials Theory grant DMR-1728457. Additional support was provided by NSF-EP- SCoR award OIA-1541079 (N.W. and K.-M.T.) and by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0017861 (Y.Z., J.M.).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Correspondence and requests for data should be addressed to K.-M.T.

**Acknowledgments:** This work used the high performance computational resources provided by the Louisiana Optical Network Initiative (http://www.loni.org) and HPC@LSU computing, accessed from 1 January 2019 to 1 January 2022.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

