Metal-Insulator Transition in Nanoparticle Solids: Insights from Kinetic Monte Carlo Simulations

Progress has been rapid in increasing the efficiency of energy conversion in nanoparticles. However, extraction of the photo-generated charge carriers remains challenging. Encouragingly, the charge mobility has been improved recently by driving nanoparticle (NP) films across the metal-insulator transition (MIT). To simulate MIT in NP films, we developed a hierarchical Kinetic Monte Carlo transport model. Electrons transfer between neighboring NPs via activated hopping when the NP energies differ by more than an overlap energy, but transfer by a non-activated quantum delocalization, if the NP energies are closer than the overlap energy. As the overlap energy increases, emerging percolating clusters support a metallic transport across the entire film. We simulated the evolution of the temperature-dependent electron mobility. We analyzed our data in terms of two candidate models of the MIT: (a) as a Quantum Critical Transition, signaled by an effective gap going to zero; and (b) as a Quantum Percolation Transition, where a sample-spanning metallic percolation path is formed as the fraction of the hopping bonds in the transport paths is going to zero. We found that the Quantum Percolation Transition theory provides a better description of the MIT. We also observed an anomalously low gap region next to the MIT. We discuss the relevance of our results in the light of recent experimental measurements.

Finally, all of the events that may have changed are recalculated. In our previous work, 2 a residence time algorithm was used, neglecting updating events that may have changed due to executed transitions. Simulation is stopped when the measured physical observable reached a steady state value within a user-defined threshold.
We start the simulation by randomly placing charges on NPs with predefined density then we switch on the KMC algorithm. The mobility is measured as: harvested charges × L z t × total number of carriers × F ext where L z is the length of the simulation box in the conducting direction. As stated above, we stop the simulation once the mobility reaches steady state: typically millions of time steps are needed to reach convergence. We used periodic boundary conditions in all three Cartesian directions and the number of harvested charges were measured by counting the net number of electrons crossing the z = L z plane.
In the current work, we extend the previous method 2 by introducing a metallic hopping channel. Concretely, we calculate the energy difference associated with the event. If the energy difference is greater than a preselected Overlap Energy (OE), then the regular thermal hopping is used. Otherwise a metallic, non-activated form is used.
For the regular hopping, we use Miller-Abrahams (MA) thermally assisted nearestneighbor hopping or tunnelling. Other approaches include the Marcus theory of electron transfer, 3 which also takes into account nuclear relaxation effects after the hopping, and the model developed by Nelson and Chandler that closely resembles Marcus theory. 4 The validity and differences of these approaches have been analyzed in detail. 5,6 The attempt frequency ν and the metallic transition prefactor ν' are assumed to be size and ligand independent and they set the time scale of the simulations. We chose ν in order to qualitatively match the order of magnitude mobitilies measured by the Matt Law group. 7 ν', however, is chosen such that the metallic transition is at least twice as fast as the regular hopping for a typical charge transfer path.
Here, ν' is assumed to be temperature independent, but in principle it could have temperature dependence similar to the mobility of the parent bulk material. ∆V is the potential difference induced by external field between source and target NPs. β is tunnelling amplitude and we evaluate it in the WKB approximation: Here ∆x is the NP-NP surface-to-surface distance, which in practice is chosen to be twice the ligand length. m * is the effective mass of the tunnelling medium, which also depends on the effective mass of the barrier. Here we approximated m * with the effective masses of electrons and holes in bulk PbSe. 8 An alternative approach is to use the Bardeen formula of tunnelling. 9 Here we refer everything to the vacuum level E vac which is thus set to zero in all simulations. If the NP solid was embedded in a matrix this would represent the conduction band minimum of the embedding matrix. E tunnelling is the tunnelling energy.
It is not immediately clear what energy should be used for E tunnelling . In the spirit of the thermally assisted hopping approach of Chandler & Nelson, E sp was defined as an average of the energies of initial and final states of the hopping: E tunnelling = (E sp a + E sp b )/2, where E sp is the single particle energy.
The energy difference ∆E ab in Eq.1 is the barrier for hopping, which can be written as: where the first term on the RHS is the difference in single particle energies of the initial and final states of the hopping: We used the energies from kṗ perturbation theory as obtained by Kang&Wise. 10 We then applied a rigid shift to align the infinite diameter limit of the conduction band edge to the work function of bulk PbSe. 11 The second term on RHS of Eq. 3 is the contribution from the external voltage V : where L z is the length of the NP solid in the conducting direction and z b and z a are the z position of the center of the NPs. Care is exercised to make sure simulations are always in linear I-V regime. In particular, we set the external voltage so that |E F ab | = 0.1kT (@30K).
As mentioned in the main text, the disorder of the NP energies did not exactly average to zero in our samples, and thus generated an internal bias field. We eliminated this bias by always taking the pairwise average of the currents with a forward and a backward applied voltage.
Finally, the last term is due to the on-site Coulomb interaction: where we introduced self energy, or the on-site charging energy: Σ 0 is the energy that needs to be paid upon the load of the first charge onto the neutral NP, while Σ is the energy it takes to load each additional charge. Both of them can be written in the form of Σ(d particle ) = q 2 /2C(d particle ), where is C is the self-capacitance of the NP. Some groups also include here the mutual capacitance of the array further decreasing the self-energy. 12 The capacitance can be taken to be proportional to the diameter d. This is the approach we followed in our previous work and X C was chosen according to the work of Zunger. 13 In this work, instead, we use Delerue's model, 14 which provides a semi-analytic form for Σ and for Σ 0 : We used the Maxwell-Garnett (MG) effective medium approximation 15 to compute the dielectric constant of the entire NP solid. According to MG, the dielectric constant of the solid can be approximated as where κ is 2 for spherical NPs and f is the filling factor.
Determining the dielectric constant of NPs is a field on its own. 16 Having calculated the total energy difference ∆E ab , we then compare it against a preselected Overlap Energy (OE). If ∆E ab is greater than OE, the hopping is considered to be an activated one and rate is calculated as shown in Eq.1.
Having defined the energetics of the NPs we can now discuss the transition, or hopping rates. The probability of an electron transferring from the initial NP a to the NP b is: where i/j denote kinetic energy levels, g i /g j are their degeneracy and f i /f j is the Fermi occupation function, n is the number of electrons on the respective nanoparticle.
Since we are assuming that the NP solid is weakly charged and the average charge per nanoparticle is less than the band degeneracy, the Fermi occupation functions can be replaced by their zero temperature limit and the sum over bands is limited to the first states: Following earlier works and assuming that we our NP solid is not operating in the extremely confined size regime, the degeneracy g of the band edge states comes from the valley-degeneracy of bulk PbSe which is eight, including spin. Later work of Delerue showed that there is minor split of these states due to intervalley coupling. 20 The above described KMC framework was coded in Cython. Random numbers were generated with the Mersenne Twist Pseudorandom Number Generator. 21 We also need to define our NP solid. The most common approach is to build a cubic lattice of NPs. This may already well capture the physics as was shown in several cases. 2,9,22 However, (i) in experiments highly monodisperse NPs tend to form denser fcc or hexagonal lattices (ii) the mobility may be underestimated due to the underestimation of the number of nearest neighbors, (iii) since most colloidal NPs have a substantial size dispersion cubic lattices of size-disordered NPs, with NPs clamped to ideal lattice positions, may contain unrealistic voids.
In order to investigate more realistic lattice disordered NP solids we set up random closed packed models by using an event-driven Molecular Dynamics code. 23,24 In either case, NP diameters (d) were drawn from a Gaussian distribution with an average diameter µ and standard deviation of σ: In order to sufficiently capture disorder effects we averaged over one hundred different random NP lattices. Error bars in our calculations represent the standard deviation of the mean. Table S1 summarizes the key parameters of the methodology used in our study.
Lengths and widths of the samples were chosen to sufficiently capture disorder effects. In particular, in our earlier study we found that it is necessary to include at least 50 sites in the conducting direction. 2 We typically studied systems with a thickness of 4 layers. The lateral transverse direction was not particularly wide because generating the sample for larger sizes was a limiting factor for the code used to generate samples with high-packing density. We often used 4 NP wide samples. We chose the dimensions of our size disordered models accordingly. However, it was increasingly hard to generate reliable packing of NP solids and run the KMC simulations while keeping the density the same thus we halved the distance in the conducting direction for the charge transport calculations in the paper. However, in order to explore the size dependence of the percolation probability, we also generated 100 samples of NP solids with 800 NPs and determined the percolation probability in both sets of samples with a Union-Find algorithm. The results are shown in Fig. S1. Figure S1: Percolation probability, averaged over 100 samples with 400 NPs, and with 800 NPs.
When comparing Fig. S1 to Fig. 1 in the paper with 400 NPs, we see that the threshold for the onset of percolation showed only minimal size dependence. This onset of the percolation transition was the key parameter used in critically analyzing the type of the MIT: see the dotted line in Figure 5, and therefore is one indicator that the size dependence of the conclusions of our analysis may be limited.
In the last column of the Table S1 we also showed the packing density of the size-, and lattice-disordered models. Let us recall that in the case of monodisperse spheres cubic lattice is the loosest packing with packing density of ≈ 0.52. The densest monodisperse packing (with coordination number of 12) has a packing density of ≈ 0.74. Packing densities in between are associated with loose or the mentioned close random packing. The simulations presented in the paper were obtained for samples with an average packing density of 0.64 and 100 random samples.