Structural Relaxation and Aging Scaling in the Coulomb and Bose Glass Models
Abstract
We employ Monte Carlo simulations to study the relaxation properties of the twodimensional Coulomb glass in disordered semiconductors and the threedimensional Bose glass in typeII superconductors in the presence of extended linear defects. We investigate the effects of adding nonzero random onsite energies from different distributions on the properties of the correlationinduced Coulomb gap in the density of states (DOS) and on the nonequilibrium aging kinetics highlighted by the density autocorrelation functions. We also probe the sensitivity of the system’s equilibrium and nonequilibrium relaxation properties to instantaneous changes in the density of charge carriers in the Coulomb glass or flux lines in the Bose glass.
1 Introduction
Disordered semiconductors near a metalinsulator transition have been the focus of experimental and theoretical investigations for the past two decades Efros1984 (); Efros1985 (); Pollak2013 (). Experimental analysis of the conductivity of a twodimensional silicon sample revealed nontrivial nonequilibrium relaxation properties: Jaroszynski and Popović drove the system out of equilibrium by altering the gate voltage, which is equivalent to suddenly changing the charge carrier density. They monitored the conductivity relaxation and discovered that the relaxation rate of the sample depended on the waiting time popovic20071 (); popovic20072 () before performing any measurements; this served as evidence for physical aging in such systems and therefore as an indication of their complex relaxation dynamics.
Our understanding of experimental results and theoretical studies investigating such strongly correlated disordered materials stems from the Coulomb glass model that was introduced by Efros and Shklovskii Efros1984 (); Efros1985 (); Efros1975 (). This model consists of charge carriers localized at random positions and interacting via strong longrange repulsive forces. Charge carriers are exposed to pinning defect sites and they reach thermal equilibrium by rearranging themselves via variablerange hopping. The strong anticorrelations originating from the longrange repulsive interactions cause the density of states to vanish near the chemical potential, creating a soft “Coulomb gap” which in turn affects the carrier transport and equivalently the system’s conductivity. The correlationinduced soft Coulomb gap was confirmed by electron tunneling experiments in doped semiconductors Massey1995 (); Butko2000 (). In addition to the correlationinduced equilibrium properties, this model displays complex nonequilibrium relaxation dynamics that has motivated careful theoretical investigations. Monte Carlo simulations performed by Grempel et al. Grempel2004 (); Kolton2005 () discovered that the autocorrelation function measurements for different waiting times can be collapsed onto a master curve for various temperatures; an evidence of aging effects and specifically dynamical scaling Henkel2010 (). This theoretical study motivated Ovadyahu to investigate the emergence of aging effects in indiumoxide films after applying a “quenchcooling” protocol where the sample’s temperature is lowered and the system then relaxes at the new temperature for a certain waiting time Ovadyahu2006 (). The gate voltage is then abruptly changed exciting the sample and the conductance is henceforth measured to confirm that the system’s relaxation displays aging effects. Similar experimental work has been carried out by Grenet and Delahaye to study the outofequilibrium properties and specifically aging in disordered insulators at low temperatures, focusing on insulating granular aluminum thin films. They concluded that the sample’s relaxation indeed depends on the waiting times since the quenchcooling application Grenet2010 ().
Disordered typeII superconductors similarly require deeper understanding due to the many competing energy scales and thus rich thermodynamic phases and transport properties Blatter1994 (). The most desirable property of typeII superconductors is the zero dissipation in current flow, hence one has to prevent vortex motion due to externally applied currents to eliminate the resulting Ohmic resistance. Optimal pinning mechanisms are utilized in these materials to prevent flux flow in the presence of applied currents. Among the different pinning configurations used, columnar defects have experimentally proven to be more efficient than uncorrelated pointlike disorder Civale1991 (). At low temperatures, magnetic flux lines become attached to the entire length of these linear pinning sites. This results in the emergence of the strongly pinned “Bose glass” phase FisherWeichman1989 (); Kwok1992 (); Nelson1992 (); Lyuksyutov1992 (); Nelson1993 () due to the localization of flux lines to these defect sites.
Physical aging has been observed in some experiments involving superconducting materials. Du et al. observed that the voltage response to an externally applied current in a superconducting 2HNbSe sample depended on the pulse duration, and that was an evidence of the occurrence of physical aging in disordered superconducting materials Du2007 (). On the numerical front, Pleimling and Täuber employed an elastic line model and Monte Carlo simulations to study magnetic flux lines in typeII superconductors analyzing their nonequilibrium relaxation properties, starting from somewhat artificial initial conditions where straight flux lines were randomly placed in the sample Pleimling2011 (); Pleimling2015 (). The resulting complex aging features (with identical initial conditions and parameter values) were subsequently confirmed in a very different microscopic representation of the nonequilibrium vortex kinetics through Langevin molecular dynamics Dobramysl2013 (). Thereafter, Refs. Assi2015 (); Assi2016 () utilized these Langevin molecular dynamics simulations to investigate the effects of flux line density quenches on the lines’ transverse fluctuations and found that the aging scaling exponents for these transverse fluctuations depended on the choice of initial conditions Dobramysl2013 (); Assi2015 (); Assi2016 (). The structural rearrangements of these vortex lines could not be probed in the computationally accessible time scales.
The Bose glass phase in typeII superconductors and the Coulomb glass model in semiconductors share many similarities Nelson1992 (); Nelson1993 (); Dai1995 (); Tauber1995 (). At low temperatures where the magnetic vortex lines are straight and at weak magnetic fields (, with being the matching field), the Bose glass can be mapped to the Coulomb glass by considering the magnetic flux lines localized to columnar defects as particles pinned to point defects in a twodimensional plane Nelson1993 (); Dai1995 (). The interacting flux lines in the Bose glass “hop” from their occupied defect sites to surrounding unoccupied sites. This leads to the formation of doublekinks or “superkink” excitations which are tonguelike doublekinks, each extending from one pinning site to another with a comparable energy Nelson1992 (); Nelson1993 (). Therefore, flux lines are not required to hop to nearest neighboring defect sites, which is similar to the variablerange hopping feature in the Coulomb glass. The phononmediated distancedependent tunneling term in the Coulomb glass is replaced with an elastic energy term in the Boltzmann factor reducing the probability for longdistance double kinks in the Bose glass. Due to the remarkable similarity between the Coulomb and Bose glasses and the resulting mapping, one expects both systems to display similar properties. Analogous to the Coulomb glass, the Bose glass is composed of longrange interacting “particles” and spatial disorder, and hence this model’s density of states was confirmed to similarly display a soft gap Dai1995 (); Tauber1995 (), and the nonequilibrium relaxation dynamics display rich aging and scaling properties Shimer2010 (); Shimer2014 (). Numerical and analytical studies by Somoza et al. investigated the density of states of twodimensional systems with logarithmic interactions at zero and finite temperatures Somoza2015 (). They confirmed that their results from the two approaches are in perfect agreement at zero temperature. Furthermore, they observed that at finite temperatures the density of states is highly dependent on the system’s temperature and can be scaled with this characteristic energy scale.
We employ Monte Carlo simulations to investigate the effects of random onsite energies and carrier density quenches on the density of states properties and aging scaling behavior in the Coulomb and Bose glass models. Our paper is organized as follows: The next section provides a brief theoretical background on the Coulomb glass model in disordered semiconductors and the Bose glass model in typeII superconductors in the presence of extended linear defects. It also discusses the Coulomb gap that forms in the density of states near the chemical potential due to the interactioninduced correlations. Section 3 presents an overview of the model we employ and the Monte Carlo algorithm we utilize to carry out our numerical investigations. Furthermore, it describes the protocol we perform, specifies the material parameter values we use, and presents the quantities we measure to analyze the system’s equilibrium properties and nonequilibrium and aging properties. Section 4 demonstrates the effects of adding nonzero random onsite energies on the equilibrium and nonequilibrium relaxation properties in the Coulomb and Bose glass models. Section 5 analyzes the effects of sudden changes in the density of charge carriers in the Coulomb glass model or magnetic flux lines in the Bose glass model on the density of states properties and the nonequilibrium relaxation dynamics. We finally summarize our work in the concluding section 6.
2 Theoretical Description
In this section, we describe the Coulomb glass model and the properties we analyze throughout our investigation. We thereafter move to adapt the Coulomb glass model to the Bose glass in typeII superconductors in the presence of extended linear defects.
2.1 Coulomb Glass Model
The Coulomb glass model was introduced by Efros and Shklovskii in 1974 Efros1975 () to describe localized charge carriers exposed to pinning defect centers in amorphous or doped crystalline semiconductors Efros1984 (); Efros1985 (); Efros1975 (). Since this model assumes that the localization length of charge carriers is of the same order or smaller than the mean separation between acceptor and donor sites, the carriers are confined to the randomly distributed sites. Furthermore, the carriers rearrange themselves to minimize the total interaction energy and thus reach thermal equilibrium at low temperatures. The carrier distribution rearrangements occur by variablerange hopping which is influenced by the phononassisted quantum tunneling between the acceptor and donor sites and thermallyinduced jumps over energy barriers.
The system’s Hamiltonian is given by Efros1984 (); Efros1985 (); Efros1975 ()
(1) 
where is the electron’s charge, is the dielectric constant, and and respectively denote the position vector and occupancy of defect site ( with being the total number of defect sites), and the sites are randomly distributed on a square lattice Yu1999 (). The site occupancy is restricted to or due to the strong intrasite repulsions. represents the th site’s bare energy, and thus the first term in (1) quantifies random site energies. The second term captures the repulsive Coulomb interactions. A uniform relative charge density is introduced to maintain global charge neutrality, where is equal to the total carrier density per site or the system’s filling fraction. When we consider cases with all onsite energies , the system’s Hamiltonian (1) is characterized by particlehole symmetry: and (where is any deviation from the halffilled system) are equivalent.
The resulting energy difference for a hop between two sites and is
(2) 
where is the distance between the two sites, and the interacting site energies are given by
(3) 
It is worth noting that if we replace the sites occupation numbers with the Ising spin variables , the Coulomb glass model would map into a randomsite, randomfield antiferromagnetic Ising model with longrange exchange interactions Davies1982 ().
2.2 Bose Glass Adaptation
The Coulomb glass model above can be adapted to describe the lowtemperature properties of magnetic flux lines in typeII superconductors in the presence of strong correlated columnar disorder Nelson1993 (); Dai1995 (); Tauber1995 (). These magnetic vortices become localized at the material’s extended linear defects in the lowtemperature Bose glass phase and subsequently thermal transverse fluctuations are suppressed. Therefore, the threedimensional system becomes essentially twodimensional by considering the crosssection of the superconductor where a pinning site is the location of the intersection of the superconductor’s “slice” transverse to the magnetic field direction and a columnar defect. The charge carriers’ role in the Coulomb glass model is played by the magnetic flux lines pinned to columnar defects in the Bose glass model.
The mutual repulsive Coulomb interaction between two occupied sites is now replaced with a modified Bessel function , which is essentially a longrange logarithmic potential screened on the scale of the London penetration length . The Bose glass effective Hamiltonian is then
(4) 
where represents the system’s energy scale with the magnetic flux quantum . In the current study, we only consider a dilute lowmagnetic field regime where all distances between sites are and thus .
Therefore, the energy difference for hops between two sites and in the Bose glass model is
(5) 
with the interacting site energies
(6) 
where (similar to the Coulomb glass model) and respectively denote the th site random onsite energy and occupation number, and is the total vortex density per site; corresponds to with being the matching field.
2.3 Coulomb Gap Properties
Efros and Shklovskii argued that the strong repulsive Coulomb interactions between charge carriers cause the system’s density of states to display a soft Coulomb gap in equilibrium Efros1984 (). Efros and Shklovskii’s meanfield arguments yield that the density of states vanishes at the chemical potential , which separates the lowenergy filled states and the higherenergy empty states. For repulsive interactions of the form , they predict that the density of states follows a power law
(7) 
for Efros1984 (); Efros1985 (); Efros1975 (); Davies1982 (); Grunewald1982 (); Xue1988 (); Grannan1993 (); Menashe2000 (); Menashe2001 (); Surer2009 (); Mobius1992 (), with a positive gap exponent . Their meanfield calculations, specifically applicable for wide randomsite energy distributions Efros1979 (), predict that the gap exponent is Efros1985 (); Tauber1995 (). These meanfield arguments assume singleparticle hops and that all sites are evenly distributed within the energy interval , whereas in fact multiparticle hops play an important role in the system and sites often cluster Davies1982 (). Therefore, represents a lower bound for the exponent Efros2011 (), for which Möbius, Richter, and Drittler obtained deviations from the predicted meanfield results Mobius1992 ().
The existence of the correlationinduced Coulomb gap in the density of states significantly impedes the charge carrier mobility Efros1975 () effectively lowering the sample’s conductivity: The gap exponent alters the temperature dependence of the conductivity and the system’s variablerange hopping length. In dimensions, the conductivity scales as where in the thermallyactivated transport regime at low temperatures . For a constant DOS, and thus one recovers the Mott variablerange hopping exponent in three dimensions and in two dimensions. However, in the presence of the Coulomb gap, the nonzero gap exponent yields a larger : in three and two dimensions for , confirming that the Coulomb gap effectively lowers the lowtemperature conductivity.
Meanfield arguments predict an exponential suppression of the density of states in the Bose glass phase, corresponding to an infinite gap exponent value. Therefore, the Coulomb gap in the Bose glass model strongly suppresses flux creep and impedes flux transport resulting in the very favorable effect of a reduced resistivity which scales as (with the driving currents ) in the thermally assisted creep regime Tauber1995 ().
3 Model and Algorithm Description
We employ Monte Carlo (MC) simulations to investigate the Coulomb glass model in disordered semiconductors and the Bose glass model in typeII superconductors with correlated columnar defects. This section specifies the MC algorithm we have utilized to study both models, the system parameters we considered, and the quantities we measured.
3.1 Model Description
Monte Carlo simulations are utilized to study the Coulomb and Bose glass models in dimensions. The simulation system is initiated by randomly placing pinning defect sites within a square lattice with periodic boundary conditions and charge carriers/flux lines that occupy sites, where the site occupancy is restricted to or . An example of the initial layout in one of our simulations is displayed in Fig. (a)a.
Then, charge carriers/flux lines attempt to hop from their corresponding occupied sites to unoccupied sites. The transition rate from site to site is determined by
(8) 
where denotes a microscopic time scale, represents the spatial extension of the charge carrier wave function or thermal wandering of the magnetic flux lines, and given by Eq. (2) for the Coulomb glass model or (5) for the Bose glass model is the difference between the energies of the two sites in question. Therefore, the success of hops is determined by two factors in (8): an exponential term capturing the strongly distancedependent phononmediated tunneling in semiconductors and vortex superkink proliferation in typeII superconductors, and a Metropolis term describing thermallyactivated jumps through energy barriers. These hops allow the system in Fig. (a)a to relax to a new configuration shown in Fig. (b)b after MC time steps (the MC algorithm is outlined in the following section).
3.2 Monte Carlo Algorithm
One Monte Carlo time step (MCS) in our simulations occurs as follows:

Choose a random occupied defect site .

Choose a random unoccupied defect site with probability .

Pick a random number , and attempt a charge carrier hopping from to . The probability of the transition success is given by the Metropolis term, , i.e., move the charge carrier from defect site to site if . If the hopping attempt fails, return to the first step.
One MCS consists of N iterations of the outlined steps.
3.3 Simulation Protocol
The first system we study as a continuation to the work in Refs. Shimer2010 (); Shimer2014 () is the initial hightemperature quench with completely random initial conditions. The simulation lattice with the corresponding number of defect sites and carriers is prepared, and the system relaxes for a sufficiently long initial relaxation time for both models.
The second system we consider involves an abrupt change in the number of charge carriers, where carriers are randomly added or removed from the sample. This serves as an attempt to bridge our computational work with experimental research and to better understand the system’s sensitivity to the choice of initial conditions. For this part, we let the system relax for a long time beyond any microscopic time scale . At MCS, the number of charge carriers is instantaneously increased from to or lowered from to . Then, we let the resulting system relax up to a waiting time measured after the quench. A snapshot of the system is taken at and the twotime carrier density correlation function defined in Section 3.5 is measured at later times . A similar study of flux line density quenches is performed to analyze structural rearrangements of flux lines in the Bose glass phase in disordered typeII superconductors.
3.4 System Parameters
All simulation distances are measured in units of the average distance between neighboring defect sites . Simulation energies in the Coulomb glass model are measured in units of the typical Coulomb energy , whereas energies in the Bose glass model are in units of . Furthermore, all simulation times are measured in units of Monte Carlo time steps (MCS). The spatial extension of the charge carrier wave function is chosen to be Grempel2004 (); Kolton2005 () following the careful investigation of different values of in Ref. Shimer2010 (). In the Bose glass model, the London penetration length is set to .
Earlier work by Efros et al. Efros2011 () focused on random site energies from a uniform distribution of certain width, while here we investigate different energy distributions to conclude their effects on the system’s properties. We consider cases with zero random onsite energies () and others where the onsite energies are taken from either a normalized Gaussian or a flat distribution of zero mean and different widths .
All other simulation parameters are chosen in accordance with the findings in Ref. Shimer2010 (). The simulation lattice length is , since the study of various lengths proved the absence of measurable finitesize effects. The system’s ambient temperature is set to , since the temperature range that yields reasonable dynamics is : the system’s kinetics slowed down too much below , and no aging was observed above . The filling fraction is kept within the range since the dynamics in systems with filling fractions or freezes out in the accessible simulation time scales. Periodic boundary conditions are utilized, and each simulation configuration shown was analyzed with to independent simulation runs.
3.5 Measured Quantities
The density of states (DOS) , with the interacting site energies given by (3) for the Coulomb glass model and (6) for the Bose glass model, is measured by finding the number of sites that are located in each energy bin, with the bin size chosen as . We also study the soft gap that forms in the DOS and the speed of its formation over time given different simulation configurations. Meanfield calculations predict the DOS to display a power law near the chemical potential , which directs us to measure the gap exponent and compare it to the meanfield predicted .
The nonequilibrium relaxation dynamics and response to the hightemperature quench and the charge carrier (or magnetic flux line) density abrupt changes are analyzed using the normalized twotime carrier density autocorrelation function
(9) 
where is the occupation number of defect site at time , is the waiting time, and is the filling fraction. Furthermore, denotes an average over all sites, followed by an average over several thousand independent disorder realizations and thermal noise.
4 Random Onsite Energy Effects
Via Monte Carlo simulations, we have studied the emergence of the soft gap in the density of states in the Coulomb and Bose glass models in two dimensions. We compare the gap’s properties and nonequilibrium relaxation dynamics of systems without random onsite energies to those with onsite energies from different distributions. Section 4.1 discusses the effects of adding nonzero onsite energies on the density of states and the Coulomb gap’s characteristics in the Coulomb and Bose glass models, while Section 4.2 highlights the ensuing nonequilibrium relaxation dynamics and resulting aging scaling exponents. In this section, there are random pinning sites available to charge carriers/flux lines.
4.1 Coulomb Gap Properties
4.1.1 Coulomb Glass Model
We first investigate the Coulomb gap formation in the Coulomb glass model under random initial conditions in the absence of random onsite energies, where a fixed number of charge carriers is exposed to randomly distributed defect sites. At the beginning of our simulations before any hops occur, the density of states has a maximum at the chemical potential . At MCS, a local minimum in the DOS already starts to form as a clear indication for the occupation of some lowenergy states due to the charge carriers hopping between sites. The DOS in Fig. 2 is almost totally suppressed at the chemical potential and the Coulomb gap is pronounced. This confirms the occupation of the lowenergy states by the present charge carriers, which is pictured by the peak located at , and the separation from the unoccupied higherenergy states at .
We have monitored the evolution of the DOS in time in the absence of onsite energies to gauge the effect of the inclusion of such energies on the speed of the gap’s formation. To this end, we studied systems without onsite energies and others with random onsite energies chosen from either a normalized Gaussian or a flat distribution of different widths. Fig. 3 shows the gap in the DOS starting to form within MCS, and the value of decreasing as time progresses until (an almost) total suppression is observed in our twodimensional Coulomb glass model. The threedimensional system was observed in Refs. Shimer2010 (); Shimer2014 () to display a faster total suppression of the Coulomb gap than the twodimensional system.
Fig. 3 confirms that the inclusion of random onsite energies from a Gaussian distribution of different widths (here and are utilized) does not influence the speed of the gap’s formation and it does not affect the suppression of the DOS at the chemical potential at long times. On the other hand, one observes that onsite energies from a flat distribution of zero mean and width might be slowing down the suppression of the DOS at , where reaches a plateau value after MCS of our accessible time scales as seen in Fig. 3. It is possible that the inclusion of this energy distribution might have caused the density of states to freeze out in our smallsized Coulomb glass system for the accessible time scales considered here.
Comparing Fig. 4 to Fig. 2 without onsite energies at the same time ( MCS) asserts that onsite energies influence the arrangement of charge carriers and their pinning behavior to defect sites, thereby eliminating the symmetry of the DOS around for halffilling and more importantly delaying the time when total suppression of the DOS at is reached.
To better understand the effect of nonzero random onsite energies on the Coulomb gap, we have studied the asymptotic behavior of the density of states near the chemical potential. As discussed earlier, the DOS near the chemical potential diminishes via a power law and meanfield arguments compute the gap exponent as for the Coulomb repulsive interactions of the form . This implies that the meanfield gap exponent for the Coulomb glass model in two dimensions is . We study the asymptotic behavior of the density of states near the chemical potential for the cases characterized by the absence and later presence of different onsite energies and observe that a power law behavior is displayed in the DOS near ; an example is shown in Fig. 5.
As indicated in Fig. 5, the computed effective gap exponent for the system without random onsite energies is found to be , which markedly differs from the meanfield prediction . As discussed in Section 2, it is expected for the exponent computed in these simulations to be different from that predicted by meanfield arguments. Furthermore, the meanfield prediction pertains to systems with strong disorder (wide random site energy distributions) Efros1979 (). We have measured the gap exponent in cases with different onsite energy distributions, and the values we found are summarized in Table 1.
Note that in order to estimate the characteristic error associated with the slope of the linear fit of our data, thus accounting for varying fluctuation sizes in different portions of the data, we divide the dataset into multiple regions. We first use a simple linear regression to obtain a line of best fit for the entire data range, producing the major line. We then divide the dataset into equalsized regions along the xaxis and regressively fit a straight line to each region and note the slope of the line and standard deviation in the slope, two numbers that define an interval for the slope. This gives us two slope values associated with each region, i.e., the upper and lower bound of the interval. We then compare all the slope values obtained to the slope of the major line and report the largest difference as the characteristic error associated with the linear fit to our data.
The flat random onsite energy distribution has a much more pronounced influence on the value of the computed exponent than the Gaussian energy distributions. Including random onsite energies from this flat distribution decreases the gap exponent and renders it closer to the meanfield predictions, as expected Efros1979 (). Sharp energy distributions are less effective due to the interactioninduced correlations, and one obtains results that are very similar to the case when (we confirmed this with flat and Gaussian distributions of zero mean and width ). In contrast, broader distributions with variances on the same scale as the correlationinduced Coulomb gap’s width have a more pronounced effect on the density of states around the chemical potential. A lower effective gap exponent than the case without onsite energies implies a lower , which in turn translates to an increased conductivity and thus charge carrier transport since conductivity scales as .
(a)  (b) Gaussian  (c) Gaussian  (d) Flat  

4.1.2 Bose Glass Model
We have similarly monitored the density of states in the Bose glass model in two dimensions and analyzed the Coulomb gap’s properties that characterize this system. A very pronounced difference with the Coulomb glass model is that the soft gap forms quicker and the DOS is suppressed faster, i.e., in shorter time, in about MCS. The total suppression of the DOS at for the Bose glass can be observed in Fig. 6.
The Coulomb gap forms much faster and is broader than that displayed by systems with Coulomb repulsive interactions. Analyzing the evolution of in time confirms this observation and the direct effect of the addition of random onsite energies from various distributions. Comparing the speed of formation of the Coulomb gap in systems with logarithmic interaction displayed in Fig. 7 to that in systems with Coulomb interactions in Fig. 3 proves that the logarithmic repulsive interactions drive the gap to form much faster.
Aside from the faster formation of the gap, the DOS is also strongly suppressed for systems without onsite energies and energies from flat and Gaussian distributions alike. Therefore, one concludes that the inclusion of onsite energies from different distributions does not change the overall properties and speed of formation of the gap in the Bose glass model, unlike our findings in the twodimensional Coulomb glass model where energies from a flat distribution affect the suppression of . To further investigate the effects of random onsite energies on the density of states in the Bose glass model, we have checked for the power law fit and the gap exponent in this case. Similar to the Coulomb glass model in Section 4.1.1, we observe that follows a power law in the vicinity of the chemical potential .
However, the gap exponents computed in the Bose glass are large when compared to those in the Coulomb glass. One has to keep in mind that the suppressed motion within the Bose glass model leads to a slow relaxation and thus a large gap exponent. Calculating the gap exponent from the power law best fit allows us to study the influence of random onsite energies on the density of states. We have computed the Coulomb gap exponents from the power law best fit in the Bose glass system in the absence of any onsite energies and in the presence of random energies from a flat distribution of width and a Gaussian distribution of widths or . The summary of our findings is presented in Table 2.
(a)  (b) Gaussian  (c) Gaussian  (d) Flat  

Onsite energies from a Gaussian distribution of different widths show a much weaker effect on the gap exponent than the flat distribution, similar to our findings in the Coulomb glass model. The introduction of random onsite energies from the utilized flat distribution increases the effective gap exponent rendering it closer to the infinite value that meanfield arguments predict in the Bose glass phase.
4.2 NonEquilibrium Relaxation Dynamics
In addition to our investigations of the effects of random onsite energies on the density of states, the Coulomb gap formation, and the effective gap exponent, we have analyzed the corresponding effects on the nonequilibrium relaxation properties of the systems studied following a quench from a completely uncorrelated, hightemperature initial state. Section 4.2.1 presents our findings for the Coulomb glass model in disordered semiconductors, and Section 4.2.2 summarizes our results for the Bose glass model in typeII superconductors in the presence of extended linear defects.
4.2.1 Coulomb Glass Model
We have investigated the nonequilibrium relaxation dynamics in the Coulomb glass model by measuring the twotime carrier density autocorrelation function and its properties.
In the absence of random onsite energies, the system shows slow dynamics and breaking of time translation invariance due to its dependence on the waiting time, as seen in Fig. 8. The slow relaxation of this system with Coulomb repulsive interactions is apparent in Fig. 8 since the system has not reached a steady state even after MCS. The system displays two different relaxation regimes, similar to those shown in Fig. 9 for the system with Gaussian onsite energies.
The first regime is characterized by a fast relaxation produced by the displacements of charge carriers around their initial positions and the carrier clusters’ group motion. On the other hand, the second regime displays a slower relaxation that is dominated by the divisions of charge carrier clusters and the resulting relaxation of the system towards a steady state. The slow dynamics and breaking of time translation invariance are two characteristics of physical aging, while the third property to be considered here is dynamical scaling Henkel2010 ().
The carrier density autocorrelation function in the aging scaling regime obeys the general scaling form
(10) 
where is an aging scaling exponent. When in the case of full aging, the scaling function is a power law , where is the autocorrelation exponent and is the dynamic exponent Henkel2010 ().
Investigating the dynamical scaling property in this regime, we observe that scaling collapse is obtained with the exponents and for the twodimensional Coulomb glass model without onsite energies, see Fig. (a)a. The goal of this section is to compute these scaling exponents after the addition of onsite energies from different distributions and compare them to the above case. We focus on energies from flat and Gaussian distributions, and we present the measured scaling exponents in Table 3.
(a)  (b) Gaussian  (c) Gaussian  (d) Flat  

From these results, we conclude that the addition of random onsite energies from a Gaussian distribution of different widths does not noticeably affect the aging scaling exponents. On the other hand, onsite energies from a flat distribution of width highly influence both exponents, see Fig. (b)b. The introduction of broad onsite energy distributions of widths on the same scale as the correlationinduced Coulomb gap’s width drives the system to display a faster relaxation dynamics, as confirmed by the higher .
4.2.2 Bose Glass Model
A similar study of the nonequilibrium relaxation properties of the Bose glass model in two dimensions has been performed, where we have analyzed the effects of random onsite energies on the density autocorrelation function and its scaling form and exponents.
Starting with the Bose glass model with zero onsite energies, slow dynamics is also a signature of such a system and the autocorrelation function displays a tworegime relaxation similar to that of the Coulomb glass model in Fig. 9. Investigating the universality of the scaling exponents, we study the full aging scaling forms for the various onsite energy distributions.
To monitor the effects of onsite energies on the twodimensional Bose glass system’s nonequilibrium relaxation properties, we similarly investigate any discrepancies present between the scaling exponents in the case with zero onsite energies and those in the cases with different energy distributions. The results are summarized in Table 4.
(a)  (b) Gaussian  (c) Gaussian  (d) Flat  

Our results in Table 4 confirm that adding random onsite energies from a Gaussian distribution of different widths has a minute effect on the scaling exponents and , see Fig 11, which implies that this addition hardly modifies the relaxation dynamics of the Bose glass system with essentially logarithmic interactions. On the other hand, it is evident that a flat onsite energy distribution of width has a pronounced effect on the nonequilibrium relaxation properties of this system. Both aging scaling exponents and assume considerably larger values as a response to the introduction of nonzero random onsite energies from a flat distribution, which hence imposes a faster relaxation on the system analogous to our results for the Coulomb glass model.
5 Density Quench Effects
Previous studies addressed the Coulomb gap and nonequilibrium relaxation properties under random initial conditions Grempel2004 (); Kolton2005 (); Shimer2010 (); Shimer2014 (). It is worthwhile to analyze the effect of abrupt changes in the density of charge carriers in the Coulomb glass model or analogously flux lines in typeII superconductors with linear defects on the Coulomb gap’s properties and aging scaling exponents to probe the system’s sensitivity to sudden changes in parameters. These quenches are experimentally achieved by quickly switching the gate voltage for semiconductors or the external magnetic field in typeII superconductors. We have utilized Monte Carlo simulations to study the effects of density quenches on the Coulomb gap in Section 5.1 and on the aging scaling behavior and exponents in Section 5.2, all in the absence of random onsite energies ().
5.1 Coulomb Gap Properties
5.1.1 Coulomb Glass Model
In the Coulomb glass model in two dimensions, we investigate the effects of suddenly increasing the filling fraction from to (upquench) or alternatively decreasing the filling fraction from to (downquench) on the Coulomb gap and its evolution in time.
We first keep the density of charge carriers fixed with the filling fraction and let the system relax in time for MCS until the Coulomb gap forms at the chemical potential where for halffilled systems (this is when the simulation clock is reset to ). After confirming that the density of states displays a relaxed Coulomb gap, we perform the density quench. Discussing first the case of the density upquench and measuring the DOS, we have to keep in mind that the system is not at halffilling anymore but at , which implies that we can expect a shift in the DOS as compared to the case with a fixed density. To test this expectation, we follow the DOS evolution in time.
At the moment of the charge carrier density upquench, Fig. (a)a shows that the peak on the left denoting the lowerenergy states reduces while the opposite peak increases. This occurs since new charge carriers are abruptly introduced to the system which in turn increases the typical Coulomb repulsion, and thus an enhanced peak on the higherenergy states initially emerges. As the system relaxes away from this initial arrangement and towards equilibrium, the new charge carriers start exploring their surroundings for energetically more favorable pinning sites. This results in a change of the shape of the DOS through a flip in the asymmetry, compare Fig. (a)a with Fig. (b)b.
Since after the density quench, this disappearance of the symmetry that characterizes the halffilled system is expected. In the case of an abrupt increase in away from halffilling, the equilibrium chemical potential is shifted to a new value .
The case of density downquench is similarly performed, where the filling fraction is abruptly lowered from to . At the moment of the quench at , the peak at the lower energies is enhanced, which is expected due to the reduction in the number of charge carriers which implies a decrease in the typical Coulomb repulsion. A similar but reversed asymmetry to that in Fig. (b)b therefore appears at MCS after the downquench. As the system relaxes, a smaller number of charge carriers occupies the original pinning defect sites, and hence the asymmetry is reversed to show an increased peak for the unoccupied states.
Another property that can be utilized to investigate the effects of instantaneous density changes is the power law fit for the density of states in the vicinity of the chemical potential. We compute the effective gap exponent at the same time for the cases when the filling fraction is kept fixed at half or suddenly raised/lowered away from halffilling.
(a) Fixed density  (b) Upquench  (c) Downquench  

From Table 5, one obtains the effective gap exponents for the three systems with different final filling fractions. It is worth noting that since the Coulomb gap exponent is an equilibrium property, one expects to obtain similar results when either starting from random initial conditions with filling fraction or abruptly changing the filling fraction from to and thereafter letting the system relax towards equilibrium.
5.1.2 Bose Glass Model
Having established that the DOS configuration is affected by sudden changes in the charge carrier density in the Coulomb glass model, we aim to carry out the same investigation for the Bose glass model with logarithmic interactions. This study is significant because it analyzes the effects of abrupt changes in the external magnetic field on the spatial rearrangements of magnetic vortex lines in typeII superconductors, while our earlier work in Ref. Assi2015 (); Assi2016 () focused on the flux lines height fluctuations and the structural relaxation time regime remained inaccessible.
Similar to the study above, we start with the Bose glass model with a fixed filling fraction , and we obtain the density of states in Fig. 6, where a Coulomb gap appears and the DOS is broad and totally suppressed at the chemical potential , which implies that the system here is relaxed beyond microscopic time scales. In the case of suddenly increasing the system’s filling fraction from to , one notices the same behavior in the density of states in the Bose glass where the system starts with an enhanced peak at the higherenergy states, similar to Fig. (a)a. Hence, the symmetry that is witnessed with a fixed number of magnetic flux lines at halffilling is broken due to the initial introduction of new magnetic flux lines into the system with a constant number of pinning centers. As the system relaxes to MCS, the newlyadded flux lines have explored the surrounding landscape and become pinned to the defect sites, and hence a peak on occupied sites is now enhanced. When we consider the consequences of abruptly reducing the number of flux lines due to changing the filling fraction from to , we observe the same reversed asymmetry that is witnessed in the downquench case in the Coulomb glass system.
(a) Fixed density  (b) Upquench  (c) Downquench  

We studied the effects sudden changes in the number of flux lines have on the effective gap exponent in the Bose glass model. The results are summarized in Table 6. The effective gap exponents for the three different final filling fractions in Table 6 are almost equal, which implies that this equilibrium property is similar in systems with the studied filling fractions , , and in the Bose glass model.
5.2 NonEquilibrium Relaxation Dynamics
In addition to studying the effects of charge carrier/flux line density quenches on the density of states and the soft Coulomb gap, we analyzed the effects of sudden addition/removal of carriers on the nonequilibrium relaxation properties and the aging scaling exponents: Section 5.2.1 presents the results of this study in the Coulomb glass model and Section 5.2.2 summarizes our findings for the Bose glass model, all in the absence of random onsite energies ().
5.2.1 Coulomb Glass Model
In the twodimensional Coulomb glass model, we relax the system for MCS (this is when the simulation clock is reset to ) to then quench the density of charge carriers and measure the resulting twotime carrier density autocorrelation function for different waiting times measured after the quench.
In the case where the density is kept fixed at , the dynamics in Fig. 13 becomes less influenced by the choice of waiting times due to the long initial relaxation time that this finite system undergoes. The density autocorrelation function of this relaxed system with a fixed does not acquire the dynamical scaling property, regardless of the choice of scaling exponents, possibly due to the finitesize effects that emerge at the time scales considered. When the charge carrier density suddenly changes, we observe different features to the case with a fixed filling fraction .
For both cases when (a) the charge carrier density abruptly increases from to and (b) the density abruptly decreases from to , the dynamics is again highly dependent on the waiting times chosen as can be seen in Fig. (a)a. Furthermore, the case of density downquench shows a similar relaxation dynamics to that of the upquench case in Fig. (a)a.
Dynamical scaling is not present in the relaxed system with a fixed carrier density, while a similar system with random initial conditions displays dynamical scaling with the scaling exponents and . Investigating the dynamical scaling property after the influence of sudden changes in the density of charge carriers reveals that aging scaling is again a property of the system. Due to the similar relaxational dynamics in the cases of density upquench and downquench, we confirm that the dynamical aging scaling exponents and are equal in both cases. This dynamical scaling behavior is displayed in Fig. (b)b, from which we obtain the values and .
5.2.2 Bose Glass Model
The relaxational dynamics of the Bose glass model in two dimensions was also analyzed with a similar goal of studying the effects of sudden changes in the density of magnetic flux lines.
When the density of flux lines is kept fixed at and the system is relaxed for a sufficiently long initial time MCS, the relaxation of the autocorrelation function is very weakly dependent on the waiting times chosen, see Fig. 15.
On the other hand, when the flux line density is suddenly increased from to or decreased from to , we confirm that the system is again dependent on the waiting times . Moreover, both systems similarly display an initial fast relaxation regime and a slower later regime, see Fig. 16.
Dynamical scaling is not present in the relaxed system with a fixed flux line density, while a similar system with random initial conditions displays dynamical scaling with the scaling exponents and . In situations with abrupt increases or decreases in the magnetic flux line density, the following values for these scaling exponents are found: and .
The equal aging scaling exponents and in the Bose glass systems when the magnetic flux line density is abruptly increased or decreased implies that this system undergoes similar relaxational dynamics when new flux lines are added to the sample or after some flux lines are removed from the sample due to particlehole symmetry. This observation was also made in the Coulomb glass model with sudden changes in the charge carrier density.
6 Conclusions
We have employed the Coulomb glass model in disordered semiconductors and adapted it to the Bose glass model in typeII superconductors in the presence of extended linear defects to investigate the density of states and the nonequilibrium relaxation properties in the twodimensional Coulomb and Bose glass systems via Monte Carlo simulations.
Earlier investigations focused on these systems with zero onsite energies Grempel2004 (); Kolton2005 (); Shimer2010 (); Shimer2014 (), thus one goal of this study is to analyze the effects of adding random onsite energies from different distributions into the Coulomb and Bose glass systems with a fixed charge carrier/flux line density. We conclude that adding onsite energies from a Gaussian distribution of zero mean and of various widths does not noticeably affect the speed of formation of the Coulomb gap or the effective gap exponent in both models. On the other hand, a flat distribution of a similar width to the one utilized for the Gaussian distribution causes the suppression of the density of states at the chemical potential to become slower in the Coulomb glass model (keeping in mind the possibility that the flat onsite energy distribution might have caused the density of states to freeze out in our smallsized Coulomb glass system in the accessible simulation time scales), while not affecting this feature in the Bose glass model. Furthermore, adding random onsite energies from this flat distribution causes the gap exponent’s effective value to approach the meanfield predictions in both models. Such broad onsite energy distributions of widths on the same scale as the correlationinduced Coulomb gap width affect the density of states near the chemical potential in these two models.
The nonequilibrium relaxation dynamics shows the effects of nonzero random onsite energies as well: The Gaussiandistributed onsite energies do not change the aging scaling exponents, while a wide flat distribution has a pronounced effect on these exponents in both systems. Including random onsite energies from a flat distribution into the Coulomb and Bose glass models causes their relaxation to become faster. This change in the aging scaling exponents implies that these exponents are not universal and in fact depend on various microscopic system parameters Shimer2010 (); Shimer2014 ().
Another goal of this study is to analyze the effects of abrupt changes in the density of charge carriers in the Coulomb glass model or magnetic flux lines in the Bose glass model to gauge the systems’ response to sudden perturbations. In the cases of density upquench and downquench in both models, the system first responds to the quench by shifting the whole DOS curve to higher (upquench) or lower (downquench) energies and displaying an asymmetry in the density of states around its minimum value. When the modified system reaches equilibrium, the larger peak shifts to the opposite side of the new chemical potential due to the relaxation of the collective system with the newlyadded carriers or the relaxation of the same system after removing some carriers. Since the Coulomb gap exponent is an equilibrium property, one should of course obtain the same exponent in a system with a specific filling fraction regardless of the utilized initial conditions.
We analyzed the effects of density quenches on the nonequilibrium relaxation dynamics and the aging properties in the Coulomb and Bose glass models. Systems with a fixed density relax for a long time until time translation invariance is displayed by the carrier density autocorrelation function. We observe that dynamical scaling is not a property of finite Coulomb and Bose glass systems at very long times when the density is fixed. On the other hand, when the density is suddenly raised or lowered, the system displays dynamical scaling and aging scaling exponents that are equal in the cases of density upquench and downquench in both models. The equal aging scaling exponents when the density of charge carriers or flux lines abruptly changes confirms that the corresponding systems undergo similar relaxational dynamics as a response to the addition of new carriers or the removal of an equivalent number of carriers due to particlehole symmetry. Within error bars, the aging scaling exponents in the Coulomb and Bose glass models with random initial conditions and after sudden changes in the carrier density coincide, suggesting that structural rearrangements of charge carriers/flux lines are governed by the lowenergy equilibrium site energy distributions inside the Coulomb gap.
Acknowledgments
We gladly acknowledge stimulating and helpful discussions with Boris Shklovskii, Matt Shimer, and Shadi Esmaeili. This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, under Grant No. DEFG0209ER46613.
Authors Contribution Statement
Hiba Assi carried out the numerical work, analyzed the data, and wrote the first manuscript draft. Harshwardhan Chaturvedi assisted in the computational part and coding. Michel Pleimling coadvised the research, contributed to data analysis and physical discussions, and assisted in finalizing the paper. Uwe C. Täuber proposed the research topic, supervised the research progress and data analysis, and edited the first manuscript draft.
References
 (1) B. I. Shklovskii and A. L. Efros, Electronic Properties of Doped Semiconductors (Springer, New York, 1984).
 (2) A. L. Efros and M. Pollak, ElectronElectron Interactions in Disordered Systems (North Holland, Amsterdam, 1985).
 (3) M. Pollak, M. Ortuño, and A. Frydman, The Electron Glass (Cambridge University Press, New York, 2013).
 (4) J. Jaroszynski and D. Popović, Phys. Rev. Lett. 99, 046405 (2007).
 (5) J. Jaroszynski and D. Popović, Phys. Rev. Lett. 99, 216401 (2007).
 (6) A. L. Efros and B. I. Shklovskii, J. Phys. C: Solid State Phys. 8, L49 (1975).
 (7) J. G. Massey and M. Lee, Phys. Rev. Lett. 75, 4266 (1995).
 (8) V. Yu. Butko, J. F. DiTusa, and P. W. Adams, Phys. Rev. Lett. 84, 1543 (2000).
 (9) D. R. Grempel, Europhys. Lett. 66, 854 (2004).
 (10) A. B. Kolton, D. R. Grempel, and D. Domínguez, Phys. Rev. B 71, 024206 (2005).
 (11) M. Henkel and M. Pleimling, Nonequilibrium Phase Transitions, Volume 2 – Ageing and Dynamical Scaling far from Equilibrium (Springer, Heidelberg, 2010).
 (12) Z. Ovadyahu, Phys. Rev. B 73, 214204 (2006).
 (13) T. Grenet and J. Delahaye, Eur. Phys. J. B 76, 229 (2010).
 (14) For a now classic review, see: G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, Rev. Mod. Phys. 66, 1125 (1994).
 (15) L. Civale, A. D. Marwick, T. K. Worthington, M. A. Kirk, J. R. Thompson, L. KrusinElbaum, Y. Sun, J. R. Clem, and F. Holtzberg, Phys. Rev. Lett. 67, 648 (1991).
 (16) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
 (17) W. K. Kwok, S. Fleshler, U. Welp, V. M. Vinokur, J. Downey, G. W. Crabtree, and M. M. Miller, Phys. Rev. Lett 69, 3370 (1992).
 (18) D. R. Nelson and V. M. Vinokur, Phys. Rev. Lett. 68, 2398 (1992).
 (19) I. F. Lyuksyutov, Europhys. Lett. 200, 273 (1992).
 (20) D. R. Nelson and V. M. Vinokur, Phys. Rev. B 48, 13060 (1993).
 (21) X. Du, G. Li, E. Y. Andrei, M. Greenblatt, and P. Shuk, Nature Phys. 3, 111 (2007).
 (22) M. Pleimling and U. C. Täuber, Phys. Rev. B 84, 174509 (2011).
 (23) M. Pleimling and U. C. Täuber, J. Stat. Mech. (2015) P09010.
 (24) U. Dobramysl, H. Assi, M. Pleimling, and U. C. Täuber, Eur. Phys. J. B 86, 228 (2013).
 (25) H. Assi, H. Chaturvedi, U. Dobramysl, M. Pleimling, and U. C. Täuber, Phys. Rev. E 92, 052124 (2015).
 (26) H. Assi, H. Chaturvedi, U. Dobramysl, M. Pleimling, U. C. Täuber, to appear in Molecular Simulation (2016) [arXiv:1509.02227].
 (27) U. C. Täuber, H. Dai, D. R. Nelson, and C. M. Lieber, Phys. Rev. Lett. 74, 5132 (1995).
 (28) U. C. Täuber and D. R. Nelson, Phys. Rev. B 52, 16106 (1995).
 (29) M. T. Shimer, U. C. Täuber, and M. Pleimling, Europhys. Lett. 91, 67005 (2010).
 (30) M. T. Shimer, U. C. Täuber, and M. Pleimling, Phys. Rev. E 90, 032111 (2014).
 (31) A. M. Somoza, M. Ortuño, T. I. Baturina, and V. M. Vinokur, Phys. Rev. B 92, 064201 (2015).
 (32) C. C. Yu, Phys. Rev. Lett. 82, 4074 (1999).
 (33) J. H. Davies, P. A. Lee, and T. M. Rice, Phys. Rev. Lett. 49, 758 (1982).
 (34) M. Grünewald, B. Pohlmann, L. Schweitzer, and D. Würtz, J. Phys. C: Solid State Phys. 15, L1153 (1982).
 (35) W. Xue and P. A. Lee, Phys. Rev. B 38, 9093 (1988).
 (36) E. R. Grannan and C. C. Yu, Phys. Rev. Lett. 71, 3335 (1993).
 (37) D. Menashe, O. Biham, B. D. Laikhtman, and A. L. Efros, Europhys. Lett. 52, 94 (2000).
 (38) D. Menashe, O. Biham, B. D. Laikhtman, and A. L. Efros, Phys. Rev. B 64, 115209 (2001).
 (39) B. Surer, H. G. Katzgraber, G. T. Zimanyi, B. A. Allgood, and G. Blatter, Phys. Rev. Lett. 102, 067205 (2009).
 (40) S. D. Baranovskii, A. L. Efros, B. L. Gelmont, and B. I. Shklovskii, J. Phys. C: Solid State Phys. 12, 1023 (1979).
 (41) A. Möbius, M. Richter, and B. Drittler, Phys. Rev. B 45, 11568 (1992).
 (42) A. L. Efros, B. Skinner, and B. I. Shklovskii, Phys. Rev. B 84, 064204 (2011).