Home | About | Journals | Submit | Contact Us | Français |

**|**Scientific Reports**|**PMC5353753

Formats

Article sections

Authors

Related links

Sci Rep. 2017; 7: 44044.

Published online 2017 March 16. doi: 10.1038/srep44044

PMCID: PMC5353753

Received 2016 September 19; Accepted 2017 January 31.

Copyright © 2017, The Author(s)

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

We study numerically the voltage-induced breakdown of a Mott insulating phase in a system of charged classical particles with long-range interactions. At half-filling on a square lattice this system exhibits Mott localization in the form of a checkerboard pattern. We find universal scaling behavior of the current at the dynamic Mott insulator-metal transition and calculate scaling exponents corresponding to the transition. Our results are in agreement, up to a difference in universality class, with recent experimental evidence of a dynamic Mott transition in a system of interacting superconducting vortices.

Materials exhibiting electric field-driven (dynamic) Mott metal-insulator transition (MIT) have a high potential for replacing semiconductors due the unique property of controllable energy gap, making them extremely promising for future low-energy electronics. While the physical mechanism behind dynamic MITs in most experimental systems is still unclear^{1}^{,2}, several theories have been proposed, including avalanche breakdown^{3}^{,4} and Schwinger-Landau-Zener tunneling^{5}^{,6}^{,7}^{,8} associated with parity-time symmetry-breaking^{9}.

In this Letter we investigate a classical system of long-range interacting charged particles experiencing voltage-induced Mott MIT near half-filling on a square lattice at small temperatures. At sufficiently low applied voltage across the system, particles arrange themselves in a checkerboard pattern. Strong inter-particle interaction impedes any motion at exactly half-filling, forming a Mott insulator. This state can be broken by either increasing temperature (thermodynamic transition) or by applying sufficiently strong external electric field or voltage, causing a dielectric breakdown characterized by finite conductivity. We observe scaling behavior at the dynamic Mott MIT with differential conductivity of the system being a universal function of , where *V* is the applied voltage, |*f*−*f*_{c}| is deviation from commensurate particle density, *f*_{c}=0.5, *ε* is some scaling exponent, and *V*_{c} is the critical amplitude of applied voltage inducing the transition.

We study a lattice gas model with long-range Coulomb interactions on a two-dimensional square lattice, with energy of the system given by the expression

where *n*_{i}=0, 1 represents the particle occupation number of site *i*, and *r*_{ij} is the distance between sites *i* and *j*. At low temperatures *T*→0 and half-filling *f*=0.5 the system exhibits Mott localization and displays a corresponding checkerboard charge order pattern (periodic boundary conditions make our system a compact surface of a 3D torus, thus Earnshaw’s theorem stating that a stable configuration of free charges interacting according to the Coulomb law is impossible, does not hold)^{10}. To realize the dynamic Mott transition we apply an electric field along the *x*-direction,

Here *V* is the electric potential and *x*_{i} is the *x*-coordinate of the site *x*. In the remainder of this section we will describe how we simulated this model.

We consider a square lattice of linear dimension *L* with periodic boundary conditions and take into account the long-range nature of Coulomb interaction by employing the Ewald summation method^{11}. In the two-dimensional case, the Ewald sum is split into a constant energy term and two rapidly converging sums over real and reciprocal space, correspondingly:

where **n** and **m** are integer vectors.

We performed Monte Carlo simulation with the heat-bath local update algorithm. At each computational step one randomly chosen particle is proposed to move to one of its neighboring sites. The acceptance probability is , where Δ*E* is the corresponding change in energy. Since we aim at relating our study with experimental results of ref. ^{1} obtained in the conditions of thermally activated dynamics, the explicit inclusion of quantum tunneling is not necessary. To promote particle conductivity, the electric field applied along the *x* axis of equation (2) is modeled by including in Δ*E* a lowering (raise) by *V* if the suggested move is to the right (left). Note that due to periodic boundary conditions, the total energy of a particle configuration is only defined up to a multiple of *VL*. Because the Monte Carlo simulation is only dependent on energy *differences*, however, this does not pose a problem.

In what follows we will assume that all particles have a unit charge, and take the lattice spacing as a unit of distance, making *V* also a measure of applied electric field in dimensionless units. To calculate the current generated during simulations, we count the number of particles crossing the *x*=0 line per single Monte Carlo sweep, where one sweep is defined as *L*^{2} proposed moves. Note that the current measured this way is limited by the number of particles present in the system and, therefore, has an unphysical upper bound. This fact limits the validity of our approach to studying particle conductivity at low voltages. To what extent such saturation effects influence the results can be probed by checking the acceptance rate of moves in the *x* direction.

Each complete Monte Carlo simulation has been performed at a fixed overall particle density *f* in two stages. First, the system’s thermal equilibrium state was reached by annealing in the absence of an external applied voltage at temperature *T*=0.04, followed by incremental increases in the applied electric field by *dV*=1/600 and measurements of particle conductivity at each voltage.

Our results were obtained for lattice size *L*=36 with 628 to 668 particles, corresponding to the range of densities . Each data point presents an average over 2.88 million Monte Carlo sweeps. We also studied the system at smaller sizes *L*=12, 24 which gave similar results, however, for clarity we will only present the data for the largest lattice size *L*=36. Differential conductivity data, *dI*/*dV*, were obtained from *IV* curves by a five-point stencil.

From the calculated *dI*/*dV* curves, we find the critical voltage *V*_{c}=0.233±0.005 near half-filling, see Fig. 1, below which the system behaves as an insulator, and above which it is conducting with significant non-zero current flowing through the system. Simulations revealed that in the immediate vicinity of *f*=0.5 (region of absent data in Fig. 1), particle current is mostly generated by an avalanche-like motion of melted clusters in the checkerboard arrangement and not by excitation of individual particles driving the dynamic MIT. We found that current, as a function of particle density *f*, experiences discontinuity at *f*=0.5, i.e. . A study of this collective effect lies outside the scope of the present Letter and will be considered elsewhere.

To show that the field-driven MIT considered here is a phase transition, we study the behavior of charge current generated in the system by applied transverse voltage. As expected for a phase transition, we observe power-law scaling of the current both as a function of applied voltage,

and as a function of particle density,

see Figs 2 and and3,3, correspondingly.

(**a**) Current, *I*, as a function of applied voltage, *V*, near the dynamic Mott MIT for a range of densities near *f*_{c}, *f*_{c}−*f*<0.015. Red line represents the best fit of data in the form (4) for *f*=0.499. **...**

In Fig. 2a we plot the *IV* curves for *f*<0.5, revealing a sharp increase in measured current above the critical voltage *V*_{c}, which becomes progressively more pronounces as particle density *f* approaches *f*_{c}. Nonlinear regression analysis was performed for the *I( f*_{c}, *V*) data^{12} to achieve the best fit to the function (4) and resulted in the scaling exponent *β*=0.5±0.1 and critical voltage *V*_{c}=0.238±0.002, which is in full agreement with *V*_{c} determined based on the *dI*/*dV* curves from Fig. (1). The critical region corresponding to the power-law fitting (4) was determined based on the extent of the linear range of *I( f*_{c}, *V*−*V*_{c}) in double-logarithmic coordinates, see Fig. 2b. The nearest to *V*_{c} data point seems to be largely affected by fluctuations of the measured current near the transition, and, together with the data at does not belong to the critical regime of the Mott MIT transition. In fact, we observe a noticeable deviation of the scaling exponent *β* from the 0.5 value at *V*≥0.26, which can be attributed to the onset of finite system size effects.

In the regime of fixed voltage, at *V*=0.237 near the critical voltage value, we find the power-law scaling of *I( f*) in the form of equation (5) in the vicinity of *f*_{c}=0.5, see Fig. 2a, with the critical exponent 1/*δ*=0.5±0.1. Figure 2b reveals the extent of the critical regime with power-law scaling. The nearest to *f*_{c}=0.5 data points appear to be affected by the avalanche physics, while at the power-law scaling starts to deviate from the investigated critical one, given by Eq. (5).

Our main result comes from analysis of the whole set of data around the critical point, {*V*_{c}, *f*_{c}}, where we observe universal scaling behavior of the measured current in the following form:

*F*_{±}(*x*) are the scaling functions in the metallic (*V*>*V*_{c}) and insulating (*V*<*V*_{c}) phases, correspondingly. It follows from the scaling relations (4)–(5) that , and .

We have performed scaling analysis of the *I( V, f*) data in the form:

The scaling parameters were obtained by maximizing the non-adjusted coefficient of determination, *R*^{2}, for the non-linear fit model (7), which resulted in 1/*δ*=0.5, *V*_{c}=0.24 and *ε*=1.0, see Fig. 4a. The above values are in full agreement with the separate analysis based on Eqs (4) and (5). We found that the scaling relation *ε*=(*δβ*)^{−1}, cf. Eqs (6) and (7), holds with remarkable accuracy.

Due to the current-saturation effect at caused by finite system size, scaling of the upper branch (*V*>*V*_{c}) in Fig. 4a is considerably less accurate. In fact, for *V*>0.26 one obtains universal scaling behavior with *ε*=0.66, see Fig. 4b, which explains poor scaling in the transient voltage range between 0.24 and 0.26.

Following ref. ^{1}, we have analyzed differential conductance data, *dI*/*dV*(*V, f*). We have also observed universal scaling behavior in the form

The critical parameters determined from the best fit to Eq. (8) were found to be *V*_{c}=0.238 and *ε*′=1.5, see Fig. 5. It follows from Eq. (6) that the scaling exponent *ε*′ must satisfy the following relation:

which for *ε*=1, 1/*δ*=0.5 and *β*=0.5 yields *ε*′=1.5. Note that only scaling of differential conductance on the lower branch (*V*<*V*_{c}) was performed due to the saturation effect at high voltages mentioned earlier, which limited the amount of scaling-suitable data points with reasonably small errors at *V*>*V*_{c}. Additionally, noise levels in *dI*/*dV* data are significantly higher than in the raw *I( V, f*) data due to numerically calculated first derivative.

In a simple model of interacting classical particles with long-range interactions, we have observed universality of critical behavior near the transition between the insulating Mott state with checkerboard order and the conducting liquid-like state above a critical value of applied voltage *V*_{c}. The main critical exponents, *β*=0.5, 1/*δ*=0.5 and *ε*=1.5 satisfy to a remarkable accuracy scaling relations corresponding to a scaling form of the current and differential conductance around the critical point.

A few comments are in order. First, our exponents differ from the ones in a dynamical MIT in a Josephson junction array^{1}, where *ε*=2/3 at density *f*=1 and *ε*=1/2 at *f*=0.5. We attribute the difference in universality classes of these two transitions to the form of interaction potential. While classical particles considered in the present Letter interact via 1/*r* Coulomb interaction, the interaction potential between vortices in Josephson junction arrays is logarithmic. An experimental realization of true Coulomb interactions is possible by taking a two-dimensional electron gas (2DEG) and applying an external periodic potential to mimic the square lattice. Second, although the checkerboard charge-ordered system is unstable toward the electronic phase separation at any deviation from the half-filling of the lattice sites^{13}, one can expect that the critical dynamics may not be much affected by the specific charge configuration as long as the state we are interested in can be still identified as an insulator. This conclusion is supported by the results by^{14} that demonstrated the same critical BKT behavior for both regular and random systems. Yet, effects of the phase separation are important and call for further studies.

Finally, it is worth mentioning that so far most of the works on a delocalization transition have focused on disordered systems, see, for example^{15}^{,16}^{,17}^{,18}. The transitions investigated there are fundamentally different as their main mechanism is depinning, rather than the reduction of the Mott gap. The dynamical Mott transition in clean systems, as first observed in ref. ^{1}, demands a much firmer theoretical understanding.

**How to cite this article**: Rademaker, L. *et al*. Universality and critical behavior of the dynamical Mott transition in a system with long-range interactions. *Sci. Rep.*
**7**, 44044; doi: 10.1038/srep44044 (2017).

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank Ivar Martin and Lusine Khachatryan for fruitful discussions. L.R. is supported by the Dutch Science Foundation (NWO) through a Rubicon grant and the National Science Foundation under Grant No. PHY11-25915 and Grant No. NSF-KITP-17-019. V.V. and A.G. are supported by the US Department of Energy, Office of Science, Materials Sciences and Engineering Division.

The authors declare no competing financial interests.

**Author Contributions** A.G. and V.V. conceived the work, L.R. did the Monte Carlo simulations, A.G. instigated the research, A.G. and V.V. analyzed the data. All authors contributed to the interpretation of the results and writing the manuscript.

- Poccia N. et al. . Critical behavior at a dynamic vortex insulator-to-metal transition. Science 349, 1202–1205 (2015). [PubMed]
- Li J., Aron C., Kotliar G. & Han J. E. Electric-Field-Driven Resistive Switching in the Dissipative Hubbard Model. Physical Review Letters 114, 226403 (2015). [PubMed]
- Guiot V. et al. . Avalanche breakdown in GaTa4Se8-xTex narrow-gap Mott insulators. Nature Communications 4, 1722 (2013). [PubMed]
- Stoliar P. et al. . Universal Electric-Field-Driven Resistive Transition in Narrow-Gap Mott Insulators. Advanced Materials 25, 3222–3226 (2013). [PubMed]
- Oka T. Nonlinear doublon production in a Mott insulator: Landau-Dykhne method applied to an integrable model. Physical Review B 86, 075148 (2012).
- Oka T. & Aoki H. Dielectric breakdown in a Mott insulator: Many-body Schwinger-Landau-Zener mechanism studied with a generalized Bethe ansatz. Physical Review B 81, 033103 (2010).
- Oka T., Arita R. & Aoki H. Breakdown of a Mott Insulator: A Nonadiabatic Tunneling Mechanism. Physical Review Letters 91, 066406 (2003). [PubMed]
- Schwinger J. On Gauge Invariance and Vacuum Polarization. Physical Review 82, 664–679 (1951).
- Tripathi V., Galda A., Barman H. & Vinokur V. M. Parity-time symmetry-breaking mechanism of dynamic Mott transitions in dissipative systems. Phys. Rev. B 94, 041104 (2016).
- Rademaker L., Pramudya Y., Zaanen J. & Dobrosavljević V. Influence of long-range interactions on charge ordering phenomena on a square lattice. Physical Review E 88, 032121 (2013). [PubMed]
- Toukmaji A. Y. & Board J. A. Ewald summation techniques in perspective: A survey. Computer Physics Communications 95, 73–92 (1996).
- Marković N., Christiansen C., Mack A. M., Huber W. H. & Goldman A. M. Superconductor-insulator transition in two dimensions. Physical Review B 60, 4320–4328 (1999).
- Kagan M. Yu., Kugel K. I. & Khomskii D. I. Phase Separation in Systems with Charge Ordering. Journal of Experimental and Theoretical Physics 120, 415–423 (2001).
- Ortuño M., Somoza A. M., Vinokur V. M. & Baturina T. I. Electronic transport in two-dimensional high dielectric constant nanosystems. Scientific Reports 5, 9667 (2015). [PubMed]
- Middleton A. A. & Wingreen N. S. Collective transport in arrays of small metallic dots. Physical Review Letters 71, 3198–3201 (1993). [PubMed]
- Ladieu F., Sanquer M. & Bouchaud J. P. Depinning transition in Mott-Anderson insulators. Physical Review B 53, 973–976 (1996). [PubMed]
- Altshuler B. L., Kravtsov V. E., Lerner I. V. & Aleiner I. L. Jumps in Current-Voltage Characteristics in Disordered Films. Physical Review Letters 102, 176803–4 (2009). [PubMed]
- Ovadia M., Sacépé B. & Shahar D. Electron-Phonon Decoupling in Disordered Insulators. Physical Review Letters 102, 176802–4 (2009). [PubMed]

Articles from Scientific Reports are provided here courtesy of **Nature Publishing Group**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |