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 unclear 1,2 , several theories have been proposed, including avalanche breakdown 3,4 and Schwinger-Landau-Zener tunneling 5-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.

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 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 ≈ .
. f [0 485, 0 515]. 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.

Results
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 . 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, c c c c and as a function of 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 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 − > . ∼ V V 0 02 c 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  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: . We have performed scaling analysis of the I(V, f) data in the form: c c c c 1/ 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 > . ∼ V 0 26 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. Conclusions 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.