Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 44044.
Published online 2017 March 16. doi:  10.1038/srep44044
PMCID: PMC5353753

Universality and critical behavior of the dynamical Mott transition in a system with long-range interactions


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 unclear1,2, several theories have been proposed, including avalanche breakdown3,4 and Schwinger-Landau-Zener tunneling5,6,7,8 associated with parity-time symmetry-breaking9.

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 An external file that holds a picture, illustration, etc.
Object name is srep44044-m1.jpg, where V is the applied voltage, |f  fc| is deviation from commensurate particle density, fc = 0.5, ε is some scaling exponent, and Vc is the critical amplitude of applied voltage inducing the transition.

Model and simulation details

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

An external file that holds a picture, illustration, etc.
Object name is srep44044-m2.jpg

where ni = 0, 1 represents the particle occupation number of site i, and rij 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,

An external file that holds a picture, illustration, etc.
Object name is srep44044-m3.jpg

Here V is the electric potential and xi 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 method11. 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:

An external file that holds a picture, illustration, etc.
Object name is srep44044-m4.jpg

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 An external file that holds a picture, illustration, etc.
Object name is srep44044-m5.jpg, 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 L2 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 An external file that holds a picture, illustration, etc.
Object name is srep44044-m6.jpg. 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 Vc = 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. An external file that holds a picture, illustration, etc.
Object name is srep44044-m7.jpg. A study of this collective effect lies outside the scope of the present Letter and will be considered elsewhere.

Figure 1
dI/dV curves in the vicinity of half-filling.

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,

An external file that holds a picture, illustration, etc.
Object name is srep44044-m8.jpg

and as a function of particle density,

An external file that holds a picture, illustration, etc.
Object name is srep44044-m9.jpg

see Figs 2 and and3,3, correspondingly.

Figure 2
(a) Current, I, as a function of applied voltage, V, near the dynamic Mott MIT for a range of densities near fc, fc  f < 0.015. Red line represents the best fit of data in the form (4) for f = 0.499. ...
Figure 3
(a) Current, I, as a function of particle particle density, f, for a range of applied voltages, V. Red lines show the best fit of I( Vc, f) curves, Eq. (5) in the critical regime near fc. (b) Current as a function of deviation from critical particle density, ...

In Fig. 2a we plot the IV curves for f < 0.5, revealing a sharp increase in measured current above the critical voltage Vc, which becomes progressively more pronounces as particle density f approaches fc. Nonlinear regression analysis was performed for the I( fc, V) data12 to achieve the best fit to the function (4) and resulted in the scaling exponent β = 0.5 ± 0.1 and critical voltage Vc = 0.238 ± 0.002, which is in full agreement with Vc 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( fc, V  Vc) in double-logarithmic coordinates, see Fig. 2b. The nearest to Vc data point seems to be largely affected by fluctuations of the measured current near the transition, and, together with the data at An external file that holds a picture, illustration, etc.
Object name is srep44044-m10.jpg 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 fc = 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 fc = 0.5 data points appear to be affected by the avalanche physics, while at An external file that holds a picture, illustration, etc.
Object name is srep44044-m11.jpg 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, {Vc, fc}, where we observe universal scaling behavior of the measured current in the following form:

An external file that holds a picture, illustration, etc.
Object name is srep44044-m12.jpg

F±(x) are the scaling functions in the metallic (V > Vc) and insulating (V < Vc) phases, correspondingly. It follows from the scaling relations (4)–(5) that An external file that holds a picture, illustration, etc.
Object name is srep44044-m13.jpg, and An external file that holds a picture, illustration, etc.
Object name is srep44044-m14.jpg.

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

An external file that holds a picture, illustration, etc.
Object name is srep44044-m15.jpg

The scaling parameters were obtained by maximizing the non-adjusted coefficient of determination, R2, for the non-linear fit model (7), which resulted in 1/δ = 0.5, Vc = 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.

Figure 4
(a) Universal scaling of current according to Eq. (6). Lower branch corresponds to applied voltages below critical, V < Vc, while the upper branch shows a small range of data at voltage directly above critical, 0.24 < ...

Due to the current-saturation effect at An external file that holds a picture, illustration, etc.
Object name is srep44044-m16.jpg caused by finite system size, scaling of the upper branch (V > Vc) 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

An external file that holds a picture, illustration, etc.
Object name is srep44044-m17.jpg

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

Figure 5
Scaling of the differential conductivity dI/dV around the dynamical critical point at Vc = 0.238 and fc = 0.5.
An external file that holds a picture, illustration, etc.
Object name is srep44044-m18.jpg

which for ε = 1, 1/δ = 0.5 and β = 0.5 yields ε = 1.5. Note that only scaling of differential conductance on the lower branch (V < Vc) 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 > Vc. 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 Vc. 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 array1, 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 sites13, 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 by14 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 example15,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.

Additional Information

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