Modeling of the metal–insulator transition temperature in alio-valently doped VO₂ through symbolic regression

Vanadium dioxide (VO₂), well-known for its insulator–metal transition (IMT),^1–4 experiences a sharp alteration in both electrical conductivity and optical transmittance near room temperature.^3,4 This transition involves shifting from a high-temperature metallic state with a rutile structure to a low-temperature insulating state exhibiting a monoclinic structure.^3,4 Among materials considered to exhibit Mott physics as the basis for electronic transitions, VO₂ stands out due to its ultrafast transformation dynamics, making it a model candidate for low-power electronics applications and has been explored in various modern electronic devices^5,6 and optical switches.⁷ Understanding the mechanistic origin of the metal-to-insulator transition is crucial for precise control over the transition. Aliovalent element doping is an effective approach to regulate carrier concentration. Alteration of the charge state and local environment, provides a means to engineer and manipulate electronic and structural properties with continuous transition points.^8,9 Among various strategies, aliovalent doping with tungsten (W)^1,10,11 has been widely studied in the context of suppressing the IMT temperature.

To study the dynamic responses of such systems, one must gain insights into the origins of the intertwined relationship between structural deformation and complex electronic changes associated with the introduction of impurities into the lattice. Over the years, there have been many incumbent theories emphasizing mechanistic insights into the origin of the IMT in VO₂ from theoretical and experimental endeavors. Two major dimensions of theory exist in understanding the underlying mechanism: the Peierls transition, with emphasis on structural change associated with the dimerization of the V atoms along the rutile crystal axis and consequently opening of the gap, or the Mott transition⁴ to an insulating state, where the gap opens due to the strong Coulomb repulsion between the localized V atom 3d orbital and related dynamical effects. However, despite existing contrasting theories,^4,10 one primary commonality is the local structural changes due to the insertion of dopants in the lattice and relevant electronic changes associated with the changes in the occupancy of the d orbitals, allowing control over IMT. The incorporation of locally inserted W ions causes pileup of itinerant electrons on the neighboring V-sites, which works as a destabilizing agent.¹¹ 

While the complexity associated with experimental procedures limits the scope in these systems, quantum calculation methods, such as density functional theory (DFT), have been an extensively used tool for modeling the behavior theoretically. However, when confronted with strongly correlated systems, such as VO₂, the efficacy of existing DFT techniques in accurately predicting electronic structures encounters challenges.¹² The crux of the issue lies in the inadequacy of DFT’s exchange-correlation energy functional, which is often approximated as semi-local and static. This sometimes leads to results having large departures from experimental insights, making these models difficult to implement for gaining insight on the general trends, particularly when dealing with impure crystals.¹³ 

Decades of experimental effort have yielded data encoding crucial insights into the thermal control of transitions in these systems. These data are a potential seed for modern artificial intelligence (AI)/machine learning (ML)^14,15 algorithms to uncover complex relationships. However, most widely used AI/ML approaches demand much larger data sets and often produce black-box models, complicating the understanding of the independent contributions from different input parameters and, thus, posing interpretability challenges. Additionally, the limited quantity of data achievable from experiments is associated with dopant-driven IMT in VO₂ is primarily due to procedural complexities, impeding the development of traditional AI/ML models. Consequently, only a few empirical linear models are used.¹ A promising approach for bridging the gap between physical reasoning and data-driven methodologies is symbolic regression¹⁶ (SR).¹⁷ SR can identify nonlinear analytical expressions linking a target property to key input parameters, even with relatively small data sets. The absence of a predefined functional form in symbolic regression necessitates the use of optimization algorithms differing from conventional analytic or numerical methods. However, these expressions may provide valuable insights into the underlying processes governing the property.

In this work, we develop analytical models that capture the trends in the IMT in VO₂ influenced by different dopant control parameters through the exploration of the symbolic search space. Our data set, consisting of experimental data compiled from an extensive literature survey over the past two decades, focuses on W dopant-driven IMT of VO₂ at various concentrations. We focused on W as a dopant since it is perhaps the most widely studied impurity in VO₂ and has broad relevance to development of thermochromic materials for energy-efficient windows.¹⁸ In recent years, similar compounds have gained importance in the design of low-voltage switchable artificial neurons and oscillators for neuromorphic computing.¹⁹ We also contrast these studies against other dopants wherever possible from literature data. We observe a dual trend in the data set: one trend aligns with the prediction by the free energy-based analytical model developed by Adler et al.,^20,21 while a second trend has not been noted earlier. Drawing inspiration from the Adler model, we develop a generalized analytical expression capturing the trend in IMT against different dopant concentrations in the symbolic space. Our model captures both the trend and variability in the experimental observations and ties them to the fundamentals of semiconductor physics. Additionally, we analyze other dopant candidates and their influence toward IMT in VO₂. Our analysis reveals qualitative understanding of the preference of dopants for tuning the properties of these systems.

The symbolic regression workflow utilized in this study is illustrated in Fig. 1(a). It involves four main steps for developing a complete symbolic regression model. Initially, the data set is analyzed to identify existing trends. Data that is extremely noisy or clean based on presumed correlations can adversely affect both the accuracy of the model and the insights gained from the data. Therefore, trends with prominent segregation in features are separated, and a symbolic model is developed for each trend. The symbolic regression section of the workflow employs the PySR package,²⁴ which implements an evolutionary algorithm-based multi-population approach. The main loop of PySR operates on each population independently and utilizes a classic evolutionary algorithm based on tournament selection for individual. The goal is to identify simple analytic expressions for accurate and interpretable models. Additionally, the PySR regression package implements a simulated annealing-based approach to accept or reject mutations, which helps retain diversity and aids in convergence for a given population. During training, multiple symbolic searches are conducted with random initialization. After training, several symbolic models capturing trends in the data with various degrees of complexity are obtained. Appropriate selection of the simplest model (least complexity) based on the physics behind the captured data is crucial (see the model selection section in Methods). If a satisfactory model is obtained, we stop the symbolic search.

Here, $x d$ represents the dopant concentration $( x d ≪ 1)$⁠. There exists a point at which the free energy of the metallic phase becomes lower than that of the semiconducting phase, leading to a first-order transition.

For an expression $E$⁠, the total loss term is $l( E)$⁠, l_pred (E) is the prediction loss, and the complexity is $C( E)$⁠. With a suitable choice of “parsimony” (see note 1 in the supplementary material), extremely complex equations are avoided during the fitting process itself. After removing equations based on their complexity and MAE error, we select the final candidate based on physical criteria. The model should be monotonously decreasing and exhibit no discontinuities or sharp changes in gradients. Additionally, the model should not have any singular points at lower dopant concentrations.

The data set used for the development of the symbolic model, as depicted in Fig. 1(c), comprises different experimental data points showcasing the change in IMT of VO₂ for tungsten (W) doping up to a doping concentration of 5%. A total of 56 points are in the data set collected across a wide range of literature.^11,25–35 As can be seen in the data set, there is some variability that we have characterized. However, we carefully chose data sets from thin films of W-doped VO₂ in almost all cases to minimize microstructural effects due to bulk or Nano confined crystals. While the thin film deposition methods did vary between references such as sputtering, pulsed laser deposition, etc., we have limited to data from such films. We hope this is a more representative sub-group of experimental data despite differences in the deposition method. A closer look at the data set reveals a dual nature of transition temperature change [Fig. 2(a)]. We additionally validate the trends mathematically (refer to note 2 in the supplementary material). One trend (“trend 1”) depicts relatively low depression in IMT temperature with increasing W concentration [red points in Fig. 2(a)], while the rest of the points indicate a second trend (“trend 2”) with relatively high depression in temperature but noisy in nature. Additional data^36–39 include the rate of change (RC) (°C/at. %) and IMT temperature for various dopant species [Figs. S1(a) and S1(b) in the supplementary material].

With the analyzed experimental data, we further assess the predictability of Adler models. As seen in Fig. 2(a), the Adler model agrees relatively well with “trend 1,” indicating a relatively low depression in IMT with increased W doping concentration. However, the model disregards the rest of the points predominantly observed in experiments. The occurrence of predominant trends can be attributed to various factors. For instance, the extent of crystallinity in VO₂ differs between monocrystal samples and polycrystalline films and powders,⁴ leading to different rates of dispersion in IMT. Additionally, the emergence of intermediate phases,⁴⁰ such as the M2 phase identified owing to doping, can also significantly contribute to the variability in experimental observations. The preparation of the sample, strain in the lattice, and the mechanism for dopant-induced phase change profoundly affect IMT, as previously observed in the VO₂ phase diagram.⁴¹ The effect of W doping involves more than just the injection of electrons; it also introduces lattice strain, making IMT a complex phenomenon resulting from competitive interactions among electrons, orbitals, and lattices. Although the Adler model initially captures specific and more idealistic trends in the data (similar to single crystal), it exhibits shortcomings in its idealistic assumptions that prevent it from representing more realistic scenarios in IMT trends. The additional complexity associated with the experiments surpasses the scope of the Adler model, which is based on ideal theoretical assumptions. While the Adler model does not account for these additional factors, it may serve as a valuable baseline model, providing a direct link between IMT and semiconductor physics.

The nature of the dopant has a profound effect on dictating the IMT in VO₂. While tungsten is the most widely used dopant to modulate the IMT in VO₂, other candidates have also been explored. Here, we discuss 29 other dopants and their corresponding effects on IMT and the rate of change (RC) (reduction) (°C/at. %) utilizing an experimental data set, as depicted in Fig. 4(a), collected from different pieces of the literature.^36–39 The bubble sizes [Fig. 4(a)] indicate the dopant concentration in at. %. It is intuitive that upon doping, the dopant modifies the local VO₂ structure either electronically, structurally, or in combination. In substitutional doping, replacement of the transition metal cation V⁴⁺ often occurs in a different formal oxidation state. Based on the valency of the dopant, we can categorize the doping in three different ways: i.e., hole doping (valency of dopant less than +4), neutral doping (valency of dopant equal to +4), and electron doping (valency of dopant greater than +4). For example, W dopant with a +6 formal oxidation state provides additional electrons per dopant atom to VO₂ due to a charge transfer between W⁶⁺ and V⁴⁺ ions, causing the reduction of V⁴⁺ to V³⁺ ions.⁴⁵ In contrast, a dopant like Ti with a +4-oxidation state is unlikely to alter the electronic structure significantly due to the lack of charge transfer between Ti⁴⁺ and V⁴⁺ ions. The overall effect of the oxidation state of the dopant on IMT can be seen in Fig. 4(b). For hole doping (e.g., Cu³⁺, Fe³⁺, and F⁻), the IMT is altered by either increment or reduction, but the degree of modulation is comparatively low. For neutral doping, the IMT modulation is very small, and the mean IMT remains very close to that of the undoped system. However, electron doping significantly alters the IMT, primarily reducing it. A similar trend can be observed for the rate of change (RC) (reduction) in Fig. 4(c). The predominant RC can be observed for electron doping. To understand the qualitative effect of the electronic changes and structural changes, we compared the SR2 model for W doping in VO₂ with neutral doping of VO₂ with Ti doping data [Fig. 4(d)] collected from different literature studies.^46,47

The introduction of Ti⁴⁺ induces negligible changes in the valence band of VO₂. Without a donor or acceptor nature of doping, the mechanism of the regulation of Tc by Ti⁴⁺ doping is predominantly associated with the local structure perturbations induced.⁴⁵ At the same time, W replacement destabilizes the monoclinic phase through the disruption of the dimeric V⁴⁺–V⁴⁺ bond to form W⁶⁺–V³⁺ and V⁴⁺–V³⁺ pairs, thus lowering the energy barrier for the transition to the rutile structure. While both Ti⁴⁺ and W⁶⁺ have similar ionic radii⁴⁸ [Fig. 4(d)], it is difficult to discretize the exact electronic and structural contribution during the transition and compare. For dopants such as W, it is evident that imminent changes in the electronic structure lead to structural destabilization, altering the phase stability of the phases. This further corroborates the development of a symbolic model with the Adler model as a basis, which stems from electronic theory.

To understand the correlation between dopant characteristics and the relative effectiveness of tunability in IMT, we use a set of relevant periodic table and electronic structure-based features for dopant elements as depicted in Fig. 5. The features are based on specific periodic table-based attributes such as group, column number (atomic number), Mendeleev number, atomic weight, row, valency, ionic radius, and element block (s, p, d, and f), as well as electronic properties. We select the last two orbital attributes of the electronic configuration of the element, such as principal quantum number (n1 and n2), types of orbital (o1_s, o2_s, etc., for s and p orbitals), and number of electrons (e1 and e2) in them. Through feature-driven correlation analysis, we shed light on the attributes crucial for the choice of dopants. The feature importance is computed using SHAP (SHapley Additive exPlanations)⁴⁹ [Figs. 5(a) and 5(b)] using random forest (RF) models that are fitted to predict IMT and RC on the constructed feature space. For fitting the random forest models, a fivefold cross-validation strategy was used, and SHAP was applied to each fold separately and the results were combined afterward. As depicted in Fig. 5(a), the most important feature for predicting IMT is e1, indicating that elements with partially filled last s orbitals (low feature values) (e.g., Nb and Mo) are likely to impact IMT negatively (reduce it) compared to dopants having filled s orbitals. Apart from dopant concentrations (at. %), we see that valency can impact IMT, particularly for the elements with high valency [Fig. 5(a)] are likely to reduce it. For RC, we see a different trend, as valency has the most importance, with high-valent dopants likely to have a positive RC (reduce IMT). Other important features include atomic weight and atomic number. Additionally, we see that ionic radius impacts IMT, with elements having a low ionic radius likely to have a positive RC contribution (reduce IMT) (e.g., Boron) compared to elements with a high ionic radius (e.g., Tb). Overall normalized feature importance is given in Fig. 5(c). It is also imperative to remember that some of these features are correlated among themselves. Hence, to present the degree of correlation between features, we also showcase a dendrogram in Fig. 5(c). The height of the dendrogram represents the similarity among the features. For example, one can start from the extreme bottom left side of the dendrogram with two connected features (block_d and O2_d), which are practically the same. As one moves upward, it gets connected with o1_s. The d block elements share some degree of similarity in the electronic structure of the last s orbital. Overall we can see that valency, ionic radius, atomic number, and dopant concentration are among the crucial factors that dictate IMT in VO₂.

Modulating the insulator–metal transition (IMT) temperature through aliovalent doping is an effective method for regulating carrier concentration, offering a robust means to adjust device performance. While tungsten (W) is widely explored for this purpose, other dopants also influence the IMT in VO₂. Given the complexity of experiments and limited experimental data, we have developed symbolic models to capture IMT trends in VO₂ under various dopant control parameters. The symbolic models, namely, SR1, SR2, and SR3, encompass the dual trend observed in the data set. Specifically, SR2 and SR3, along with a correlation parameter a, effectively capture the variability in experimental conditions and transition mechanisms, aligning with the principles of semiconductor physics. These models provide valuable alternatives to experiments, with the a parameter requiring significantly fewer trials for determination, facilitating extrapolation to other dopant concentrations. Furthermore, our analysis of different dopant candidates sheds light on specific dopant features critical for tuning the properties of VO₂ systems, aiding in informed dopant selection for desired device characteristics. Our approach may be effective for exploring other phase change materials by modeling the control parameters and their effect on the phase space.

The supplementary material consists of data sets used to develop the SR model and related analysis.

S.R. acknowledges AFOSR (No. FA9550-23-1-0215) for support. Work performed at the Center for Nanoscale Materials, a U.S. Department of Energy Office of Science User Facility, was supported by the U.S. DOE, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. This material is based on work supported by the DOE, Office of Science, BES Data, Artificial Intelligence, and Machine Learning at DOE Scientific User Facilities program (Digital Twin). This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231.

The authors have no conflicts to disclose.

S. Banik: Formal analysis (lead); Writing – original draft (lead). S. V. Shriram: Data curation (equal); Writing – review & editing (supporting). S. Ramanathan: Conceptualization (supporting); Writing – review & editing (supporting). S. K. R. S. Sankaranarayanan: Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.