Introduction

Colloidal semiconductor nanoparticles (NPs) are singularly promising nanoscale building blocks for fabricating mesoscale materials that exhibit emergent collective properties. There is a growing interest to use NPs for numerous optoelectronic applications1,2, including third generation solar cells3,4 light emitting diodes5, and field effect transistors (FET)6,7.

Electron wavefunctions are localized on the individual NPs. This “quantum confinement” makes the electronic parameters tunable with the NP size, and thus makes the NP solids a very versatile platform for applications8.

However, the very same quantum confinement also suppresses the transport between NPs, and thus drives NP solids insulating. As a result, without the application of specific transport-boosting fabrication steps, the electron mobility in NP solids is often in the range of 10−2 − 10−1 cm2/Vs9,10. These mobilities are typically measured in FET arrangements. This is orders of mangitude below the mobilities that would be acceptable for electronic applications. Therefore, increasing the mobility and transport in NP solids is one of the central challenges on the way to realize the promise of NP solids.

Various experimental groups managed to boost the mobility by enhancing the inter-NP transition rate with a variety of methods, including: ligand engineering11,12,13, band-alignment engineering14,15, chemical-doping9,16, photodoping17, metal-NP substitution18, epitaxial attachment of NPs19,20, and atomic layer deposition methods21. Encouragingly, these efforts recently translated into notable progress, as NP solids were reported to exhibit band-like, temperature-insensitive mobilities, with values exceeding 10 cm2/Vs at room temperatures9,10. It is important to note that the claims of band-like transport are controversial some experiments reported data that can be interpreted as evidence for band-like transport. One of these is the relative temperature independence of the observed mobilities, in contrast to hopping insulators where an activated temperature dependence is expected. However, the absolute values of the mobilities remain relatively low compared to most metals, and this makes conservative commentators stop short of identifying this transport as metallic9,10.

On the theoretical front, there have been efforts from several groups to understand electronic transport in NP films and solids. Density functional theory (DFT)-based ab initio calculations of the energy levels of a single NP alone are already limited to only hundreds of atoms for higher-reliability methods, and a few thousands for more approximate methods by prohibitive CPU times. These translate to diameters less than 2–3 nm, whereas experimental NP diameters often exceed 5–6 nm. Next, the accurate computation of the NP-NP transition rates would require the simulation of two NPs. And even if this calculation is completed, it does not address that the NP-NP transport is not metallic but insulating; the disorder of the parameters from NP to NP; and finally the defects of the NP solids. In total, ab initio descriptions alone are very far from being capable of describing transport in NP solids. Cleary, there is a pressing need for developing mesoscopic transport simulations that somehow integrate ab initio calculations.

Shklovskii et al. have developed transport calculations for a NP array in a FET geometry, where they focused on the effects of the Coulomb interaction22. The interplay of transport and Coulomb interactions was studied in Refs. 23 and 24, albeit on very small samples. Over the last few years, our group developed the Hierarchical Nanoparticle Transport Simulator (HINTS) platform that starts with an ab initio calculation of the energetics of individual nanoparticles, then forms a NP solid of several hundred NPs, and finally simulates transport across this NP solid by a Kinetic Monte Carlo method24,25. HINTS can simulate 500–2,000 nanoparticles. A reassuring validation of HINTS emerged from simulating the dependence of the mobility of PbSe NP layers as a function of the NP diameter. The results in24,25 closely tracked the experimental results of Liu et al., who studied the electron mobility of PbSe layers in a FET geometry26. More recently, we studied commensuration effects in bilayer NP solids27.

However, these theoretical efforts only considered NP solids with homogeneous disorder: the NPs were arranged either in a close-packed glassy/jammed structure, or on an ordered superlattice (SL) with disorder only in the NP size. In contrast, representative scanning electron microscope (SEM) images, like in Fig. 1, taken of NP solids with millions of NPs, conspicuously reveal that typical NP solids are also characterized by disorder on much larger length scales. These defects, often on the µm length scale, have sizes well beyond the capabilities of any published technique, including HINTS. Therefore, there is a need for transport simulation methods that are capable of capturing mesoand macro-scale defects and their effect on transport.

Figure 1
figure 1

Common planar defects in PbSe nanoparticle superlattices. (a) Grain map of a typical region of an PbSe NP superlattice film showing the location of several types of planar defects. Blue, yellow, and red lines denote wide-angle grain boundaries, twin planes, and more complex, unclassified planar defects and defect clusters, respectively. The image is a montage of fifteen low-magnification, high-resolution SEM images. Voids, step edges, vacancies, and other types of non-planar defects are also visible in the montage. (b) Higher-magnification secondary electron image of a bicrystalline region of an epi-SL film with two (100)SL-oriented SL grains meeting at a twin plane. (c) Image of another region of the same epi-SL film showing multiple grain boundaries between (100)SL- and (01¯1)SL-oriented grains, as well as several other planar defects.

Full size image

We performed one step in this direction previously by extending our HINTS method to include percolative effects into homogeneously disordered NP solids25. This simulation captured physics on the longer length scales of percolative clusters. Our main message was that a metal–insulator transition (MIT) occurs when a percolating path of metallicconnected NPs develops across the entire sample. We described this MIT as a Quantum Percolation Transition. However, this work still did not incorporate planar defects.

Simulation methods

To answer the above needs, in this paper we report our work that boosted the capability of our HINTS platform by introducing additional hierarchical layers to capture the effect of planar defects on the transport in NP solids. First, we used HINTS to individually model a NP superlattice (SL) with one planar defect that was either a generic grain boundary or a twin plane. Second, we simulated transport across a large number of such single-defect SLs, and determined the distribution of the mobilities of the single-grain-boundary NP SLs and that of the single twin-plane NP SLs. We also determined the distribution of the mobilities of undefected NP SLs with only homogeneous NP disorder. Third, to reach a simulation scale approaching the scale of the NP solids in the experiments, we built a resistor network where the individual resistor values were taken from the three mobility distributions with predetermined probabilities. Motivated by our previous work25, we determined the resistance of the entire resistor network by changing the fraction of undefected NP SLs within the network. Finally, we analyzed our results by a finite size scaling method.

We call this boosted HINTS platform the TRIDENS: the “TRansport In DEfected Nanoparticle Solids” Simulator. With TRIDENS, we are capable of capturing the physics from atomistic length scales up to the scale of NP solids in the experiments by integrating the simulations on several hierarchical layers. The complete hierarchical structure of TRIDENS is presented below.

  1. (1)

    The energy levels of individual PbSe NPs are determined by adapting a k·p calculation within the NP diameter range of 5–7 nm. The valence band and conduction band values have been validated via comparison to optical experiments28. Here we focus on PbSe NPs because they are of considerable interest for solar applications due to their large Bohr exciton radius and small direct bulk bandgap29, and exhibit the possibly game-changing multipleexciton generation (MEG)30. For these reasons, PbSe NPs are often thought to have strong promise for solar applications.

  2. (2)

    The electron–electron interaction is included on the level of an on-site, or self-charging energy expression, EC, defined as:

    $$ E_{c} = n\left( {\sum^{0} + \frac{(n - 1)}{2}}\sum \right) $$
    (1)

    where n is the number of electrons on the NP. Σ0 is the self-charging energy of loading the first electron onto a neutral NP. Σ is the extra energy it takes to load each additional electron onto the NP due to repulsive Coulomb interaction with the (n-1) electrons already on the NP, as well as the interaction with the induced image charge.

This self-charging energy can be calculated by a variety of methods, including the semi-empirical pseudopotential configuration interaction method of Zunger and coworkers31 and the single NP empirical-perturbative hybrid calculations of Delerue32. In this paper we report results with the latter approach. In this approach

$$ \sum^{0} = \frac{{q^{2} }}{{8\pi \varepsilon_{0} R}}\left( {\frac{1}{{\varepsilon_{solid} }} - \frac{1}{{\varepsilon_{NP} }}} \right) + 0.47\frac{{q^{2} }}{{4\pi \varepsilon_{0} \varepsilon_{NP} R}}\left( {\frac{{\varepsilon_{NP} - \varepsilon_{solid} }}{{\varepsilon_{NP} + \varepsilon_{solid} }}} \right) $$
(2)

and

$$ \sum = \frac{{q^{2} }}{{4\pi \varepsilon_{0} R}}\left( {\frac{1}{{\varepsilon_{solid} }} + 0.79\frac{1}{{\varepsilon_{NP} }}} \right) $$
(3)

For the dielectric constant inside the NP, we assume that it equals the bulk high frequency dielectric constant of PbSe, taken to be 22.0. To model the dielectric constant of the medium surrounding the NP, we account for both the organic ligand shell of the NP as well as the presence of neighboring NPs. We assume that the ligands themselves have a dielectric constant of 2.0. The dielectric constant of the entire solid is then calculated using the Maxwell–Garnett (MG) effective medium approximation:

$$ \varepsilon_{solid} = \varepsilon_{ligand} \left( {\frac{{\varepsilon_{NP} (1 + \kappa f) - \varepsilon_{ligand} (\kappa f - \kappa )}}{{\varepsilon_{ligand} (\kappa f - \kappa ) + \varepsilon_{NP} (1 - \kappa f)}}} \right) $$
(4)

where κ is 2 for spherical NPs, and f is the filling factor.

We note that the long range part of the Coulomb interaction can be easily included into the calculation. The long range interactions change the nature of transport from activated hopping to Efros-Shklovski type variable range hopping. However, many experiments show that while this Efros-Shklovskii hopping dominates at low temperatures, as the temperature is raised past 150–200 K, the transport becomes dominated by nearest neighbor hopping33 Since our work focuses on temperature ranges around ambient room temperature, representing the Coulomb interactions with the on-site term only is appropriate.

  1. (3)

    We modelled the electron transitions between neighboring NPs via a Miller–Abrahams phonon-assisted hopping mechanism:

    $$\begin{aligned} \Gamma_{i \to j} =\left\{\begin{array}{ll} {{vg_{ij} \beta_{ij} \exp \frac{{ - \Delta E_{ij} }}{{\kappa_{b} T}}}} &\quad {if\Delta E_{ij} > 0} \\ {{vg_{ij} \beta_{ij} }}&\quad {if - \Delta E_{ij} < = 0} \\ \end{array} \right. \end{aligned}$$
    (5)

    where ν is an attempt frequency, chosen to be 1012 s−1, gij is the product of the initial density of states on NPi and the final density of states on NPj, and βij is the tunneling amplitude. βij is calculated using the WKB approximation as:

    $$ \beta_{ij} = \exp - \left( {2\Delta x\sqrt {\frac{{2m^{*} (E_{vac} - E_{ij} )}}{{h^{2} }}} } \right) $$
    (6)

Here ∆x is the NP-NP surface-surface separation distance. m∗ is the effective mass of the electrons in the tunneling medium, approximated as 0.05me, the effective mass of electrons in bulk PbSe. It is noted that m∗ was estimated to be 0.3me in NPs26. Evac is the vacuum energy level, set to be zero as all other energy levels are defined relative to the vacuum. Eij is the tunneling energy, taken to be the average of the initial and final states of the hopping transition: Eij = (Ei + Ej)/2, where Ei is the energy level of NPi.

  1. (4)

    To reach the length scale of hundreds of nanometers, we generated triclinic NP superlattices, of PbSe NPs with a 20 × 20x2 geometry, inspired by the 2D channel geometries of FETs used in transport experiments34. The triclinic unit cell was described with lattice constants a1 = a2 = a3 = 6.9 nm and angles α = β = γ = 99◦. The average NP diameter was 6.0 nm. Size and location disorder were introduced by assigning the NPs a diameter and lattice displacement vector according to Gaussian distributions of widths σ(diameter) = 0.4 nm and σ(location) = 0.3 nm respectively. See Fig. 2a for an example of one of these SLs. In our undefected SLs, these parameters yield a hβiji ' 0.015 ± 0.02.

Figure 2
figure 2

Top down views of the three types of the simulated NP SLs: (a) An undefected NP SL, characterized by the NPs having only size and location disorder; (b) A NP SL, containing a twin plane, as denoted by the dashed line, also with NP size and location disorder; and (c) An NP SL, containing a grain boundary, also with NP size and location disorder. NP color corresponds to NP diameter, as indicated in the colorbar on the left.

Full size image

Layers (1)–(4) are the main constituents of our HINTS platform. HINTS is suitable for capturing the effects of homogeneous disorder, i.e. disorder associated with the size and location of the NPs that varies from site to site of the NP superlattice, but does not involve planar defects. Next, we describe the additional layers of the TRIDENS that enable us to access length scales well beyond the reach of HINTS.

  1. (5)

    As a first step, we introduced a single planar defect into each generated NP superlattice (SL). We carefully analyzed SEM images of NP solids with millions of NPs of the type of Fig. 1, and determined the predominant types of defects and their statistics, such as the lengths and densities of the planar defects. Based on the SEM image analysis, the most relevant and oft-occurring planar defects were twin planes and grain boundaries.

  2. (a)

    Twin planes: the NP superlattice is mirrored across a boundary plane, creating two crystallites, or grains,which have reflected in-plane unit cells. Twinned grains can also be related by a 180◦ rotation normal to the twin plane. The twin plane is always parallel to a possible crystal face (but not any planes of complete symmetry, e.g. it is distinct from all space group symmetries), and thus requires the two grains to share some NP lattice sites along the boundary plane. See Fig. 1b and Fig. 2b. The high symmetry nature of twin planes makes them a low disorder defect, compared to the more highly disordered grain boundaries discussed below. In our generated samples the twin planes were created in a (100) in-plane oriented SL, and the orientation of the twin plane itself was randomly selected on a sample-by-sample basis from all possible crystal planes which would span the entire NP simulation SL in the x-direction (in order to bisect the SL). As an example, one such boundary orientation is that of a (01¯2)/(02¯1) twin boundary, where (01¯2) is the orientiation of the boundary plane in grain 1, and (02¯1) is the orientiation of the boundary plane in the mirrored grain 2.

  3. (b)

    Grain boundaries: the NP superlattice is again fractured by a boundary plane. However, unlike with twinplanes, the superlattice is not mirrored across the boundary plane. There are two main types of grain boundaries. 1) Tilt grain boundaries, where the in-plane SL orientation is the same in the two grains, but they are spatially rotated in-plane relative to each other. The angle of rotation can be divided into “low angle” and “high angle” regimes, where the higher the angle of rotation, the more disordered the grain boundary (with large areas of poor fit). 2) Twist grain boundaries, where the in-plane superlattice orientations of the two grains are different (rotation occurs along an axis perpendicular to the boundary plane). Such grain boundaries will result in two crystallites/grains with different in-plane superlattice orientations (e.g. a boundary between a (100)SL in-plane orientation and a (101)SL in-plane orientation).

In our generated samples, we simulated grain boundaries with a combination of tilt and twist mismatching. Specifically, the boundary plane separated grains of (100)SL in-plane orientation and (01¯1)SL in-plane orientation respectively, with the relative in-plane spatial orientation of the two grains depending on the angle of the boundary plane (chosen at random on a sample-by-sample basis). The boundary plane was always limited to angles which would span the entire NP simulation SL in the x-direction (in order to bisect the SL). This results in grain boundaries which are much more extensively disordered than twin planes, particularly when the boundary plane results in a high-tilt grain boundary. See Figs. 1c and 2c. Hereafter, we will refer to our specific combination of boundary mismatching as simply a “grain boundary”.

In total, we generated 30,000 NP superlattices, containing either one or two grains, where we varied the disorder of the NP diameters (see color code in Fig. 2), the on-site NP location disorder, and the orientation of the planar defects as viewed out-of-plane. Of these 30,000 NP SLs, 10,000 NP SLs had no planar defects, the next 10,000 NP SLs contained one twin plane, and the last 10,000 NP SLs contained a grain boundary.

Figure 1 shows other types of defects as well. We determined that point vacancies have minimal effect on the mobilities in the insulating phase. One can also see tears/rips/voids/cracks in the SEM image. NP superlattice fabrication technologies will be ready for technical application when they can minimize or eliminate such disruptive tears. For these reasons, we did not model either of these defects.

  1. (6)

    Next, we determined the electron mobility across each of the 3 × 10,000 defected NP SLs. To do this, each NP SL was populated with electrons, randomly placing them on NPs, using the Mersenne Twister, until a predetermined electron density was reached. The chance of an electron being placed on any particular NP was uniform, independent of electron occupation and NP parameters. Data was only taken well after the system achieved equilibrium. A small voltage was applied across the sample to induce transport in the linear I-V regime, with periodic boundary conditions. Finally, the electron transport was simulated by evolving time via a kinetic Monte Carlo (KMC) algorithm. The sodetermined mobilities of the 3 × 10,000 NP SLs were used to create the mobility distributions for the homogeneously disordered NP SLs, the twin-plane-supporting NP SLs, and the grain-boundary-supporting NP SLs. The first class of NP SLs will also be referred to as “undefected NP SLs”, the latter two classes as “defected NP SLs”.

  2. (7)

    To simulate NP solids on mesoscopic length scales of the order of 10 µm or longer, we generated a network of thousands of these NP SLs, some with planar defects, some without planar defects. In the present work, we filled the network with these NP SLs at Random. Each defected NP SL classical resistor network, with resistors chosen at random from the distributions determined in step 6. Which distribution the resistors were chosen from was also randomly determined, according to a pramaeter that describes the fraction of defected resistors Random numbers were generated using standard numpy libraries. Each resistor represents an NP SL with an L = 20 length planar defect. One notes that in Fig. 1 many of the planar defects are considerably longer than L = 20. Representing planar defects with longer lengths is possible in TRIDENS by placing defected NP SLs correlated along the lines of the network. Such longer range defect-correlations were not pursued in the present work, but will be included in future work. The mobilities of the individual NP SLs were drawn from the distributions determined in step (6).

With these preparations, the mobility of the overall NP solid was determined by treating this NP SL network as a resistor network. We used the Laplacian method of F. Y. Wu et al. to calculate the overall resistance across the entire network35. The electrodes were modeled as equipotential metallic strips spanning the entire length of the sample edge and, thus making them equivalent to a single node on a resistor network. These electrodes were coupled to the sample by contact resistors that were chosen according to the same rules as the bulk resistors.

  1. (8)

    Having determined the overall mobility of the network of defected NP SLs, we adapted finite size scaling methods to analyze whether this resistor network model built from defected NP SLs had a phase transition, or a crossover, and if so, what are the properties of this transition. To this end, we repeated step (7) for resistor networks of various sizes, including 32 × 32, 64 × 64, and 128 × 128. As detailed below, our finite size scaling found a percolation transition that separates a low mobility insulator from a high mobility insulator. We used finite size scaling to determine the critical properties of this transition, including the critical point, the critical exponents and the universal scaling function.

Experimental methods

Materials

Lead oxide (PbO, 99.999%), oleic acid (OA, technical grade, 90%), diphenylphosphine (DPP, 98%), 1octadecene (ODE, 90%), ethylene glycol (EG, 99.8%, anhydrous), acetonitrile (99.99%, anhydrous), hexanes (≥ 99%, anhydrous), toluene (99.8%, anhydrous), and (3-mercaptopropyl)trimethoxysilane (3-MPTMS, 95%) were purchased from Sigma Aldrich and used as received. Trioctylphosphine (TOP, technical grade, > 90%) and selenium (99.99%) were acquired from Fluka and mixed for 24 h to form a 1 M TOP-Se stock solution. Ethylenediamine (EDA, > 98.0%, anhydrous) was purchased from TCI and mixed with acetonitrile in a 1:1 volume ratio to make a 7.5 M EDA stock solution, this is a slight modification to a published procedure34.

Quantum dot synthesis

In this experimental section we adopt the alternative terminology of “quantum dots” to refer to the nanoparticles, to accommodate alternative terminologies preferred by different communities. PbSe QDs were synthesized and purified air-free using a slight modification of a published procedure34. Briefly, PbO (1.50 g), OA (5.00 g), and ODE (10.00 g) were mixed and degassed in a three-neck round-bottom flask at room temperature. The mixture was heated to 120 °C under vacuum to form dissolved Pb(OA)2 and dry the solution. After 1 h at 120 °C, the Pb(OA)2 solution was heated to 180 °C under argon flow and 9.5 mL of a 1 M solution of TOP-Se containing 200 µL of DPP was rapidly injected into this hot solution. An immediate darkening of the solution was observed, and the QDs were grown for 105 s at ∼160 °C. The reaction was quenched with a liquid nitrogen bath and injection of 10 mL of anhydrous hexanes. QD purification and SL fabrication were performed in glove boxes with < 0.5 ppm O2 content. The QDs were purified by two rounds of precipitation/redispersion using acetonitrile/toluene and stored as a powder in the glove box.

Substrate preparation

Following and slightly modifying the procedure seen in34, a single-side polished Si substrate was cleaned using 10 minutes of sonication in acetone, Millipore water, and then isopropanol, followed by drying in a stream of flowing air. The cleaned substrate was immersed in a 100 mM solution of 3-MTPMS in toluene for 1 hour to functionalize its native SiOx surface for improved QD film adhesion, then rinsed with neat toluene and dried in flowing air.

Superlattice fabrication, electron microscopy imaging

Quantum dot superlattice films were fabricated and imaged using (modified) published procedures34. An oleate-capped superlattice was prepared in a glovebox (<2 ppm O2) by drop casting 60 µL of 20 g/L dispersion of PbSe QDs in hexanes onto 7 mL of ethylene glycol (EG) in a Teflon well (3.5 × 5 × 1 cm). After depositing the QD solution, the well was immediately covered with a glass lid. The hexane evaporated over 30 minutes, resulting in a smooth, dry QD film floating on the EG surface. The glass lid was then removed and 0.1 mL of a 7.5 M solution of EDA in acetonitrile was slowly injected (5–10 sec) into the EG under the QD film using a 1 mL syringe. After 30 seconds of exposure to EDA, the resulting epi-SL film was stamp transferred to the Si substrate using a vacuum wand, rinsed vigorously with acetonitrile, and dried under flowing N2.

Scanning electron microscopy (SEM) imaging was performed on an FEI Magellan 400 XHR SEM operating at 10 kV. Grain maps were produced by stitching together fifteen 6144 × 4415 pixel images acquired at 50,000 × magnification, providing the ability to resolve individual QDs in the sample. Image stitching was performed in Adobe Photoshop. Grain boundaries, twin planes, and other planar defects were then located by eye and drawn in manually.

Superlattice samples for TEM analysis were prepared by stamping QD films from the EG surface onto holey carbon TEM grids without a carbon film coating. The use of TEM grids free of a carbon film was critical for high-quality secondary electron imaging (SEI) in the TEM. SE imaging was performed on a JEOL JEM-2800 TEM operating at 200 kV using a probe size of 1.0 nm.

Results and discussion

TRIDENS simulations

We laid the foundation of our simulation by carrying out steps (1)–(4), as done in a standard HINTS study. To carry out step (5), we generated 3 × 10,000 defected NP SLs by starting with homogeneously disordered but undefected 20 × 20 × 2 NP SLs whose shape broadly corresponded to FET geometries, and then inserted a twin plane planar defect into 10,000 NP SLs, and a grain boundary planar defect into another 10,000 NP SLs. The latter sometimes involved removing a few NPs to keep the shape of the NP SLs largely unchanged.

Next, we executed step (6) by determining the mobility distribution of the defected NP SLs. The mobility distribution for the homogeneously-disordered, undefected NP SLs is shown in Fig. 3a. The mobility distribution of the twin-plane NP SLs is shown in Fig. 3b. Finally, the mobility distribution of the grain-boundary NP SLs is shown in Fig. 3c. All three distributions were approximately normal, and could be well characterized by a mean and a standard deviation. The mobility of the undefected NP SLs was 0.42 ± 0.1 cm2/Vs, the mobility of the NP SLs containing twin planes was 0.16 ± 0.06 cm2/Vs, and the mobility of the NP SLs containing grain boundaries was 0.09 ± 0.05 cm2/Vs, as shown.

Figure 3
figure 3

Mobility distributions of the 3 types of NP SLs: (a) Undefected NP SLs, the NPs having size and location disorder; (b) NP SLs, each containing a twin plane, also with NP size and location disorder; (c) NP SLs, each containing a grain boundary, also with NP size and location disorder. Displayed in each panel is the average mobility and standard deviation of Gaussians fitted to the distributions. All simulations were performed at a temperature of T = 300 K.

Full size image

We then performed step (7) by assembling a resistive network whose individual links had mobilities selected from the above determined mobility distributions. To identify the paradigmatic aspects of the behavior of the mobility of the NP solid, we selected the links from the highest mobility undefected NP SL distribution with a probability p, and from the lowest mobility grain boundary NP SL distribution with a probability (1 − p). We used this p, the fraction of the high mobility undefected NP SLs/links as the control parameter of our study. Initially we expected that when the probability (1 − p) of the low mobility defected NP SLs becomes small, then the electrons will be able to “flow around” the low mobility links through the high mobility links. Put differently, the high mobility links will be able to approximately short out the low mobility links. Had this expectation been true, then the mobility of the NP solid should have exhibited a saturation as p approached 1.

Figure 4 shows the evolution of the mobility of NP solids with p. Visibly, our initial expectation was not confirmed as the mobility did not show a saturation as the fraction p of the undefected NP SLs approached 1.0. The high mobility links did not “short out” the low mobility links. A possible explanation is that the mobilities of the different NP SLs were not different enough for such a short-out. We checked the robustness of this result, and found the same characteristic non-saturating shape in other 2D geometries, as well as in 3D bulk networks.

Figure 4
figure 4

Mobility of a 32 × 32 resistor network whose resistors/links are chosen from two distributions, the low mobility grain boundary NP SLs mobility distribution and the high mobility undefected NP SLs mobility distribution. The fraction p of the high mobility NP SLs sweeps the 0 to 1 region. The error b ars are smaller than the data points.

Full size image

Bimodal percolation transition: Next, we investigated whether there is a percolative critical behavior as the p fraction of high mobility links is varied. The simple case of a resistor network, where the links either have a finite (electronic) mobility µ with probability p, or a non-conductive zero mobility with probability (1 − p), has been extensively analyzed. The conductivity of such a resistor network exhibits a critical behavior across the percolation critical point pc with a power law dependence µ ∝ (p − pc)t, where the critical exponent t > 1 is universal, and pc is the percolation threshold.

Remarkably, the closely related bimodal problem of the links of a network having a high conductivity σ(high) with probability p, or a low but finite conductivity σ(low) with probability (1 − p) has been rarely analyzed. Efros and Shklovskii (ES) established the broad framework for the analysis, when they made the analogy between this bimodal distribution problem and the problem of how the critical behavior of a spin system gets modified by the presence of a symmetry breaking magnetic field36. They hypothesized a power law critical behavior for the network conductivity, where the universal scaling function at the critical point p = pc is anchored by the ratio of the high and low conductivities. However, they did not determine either the critical exponents, or the universal scaling function.

Finite size scaling: Next, we attempt to adapt the ES bimodal framework to describe our TRIDENS-simulated results. The finite size L of the simulated samples makes it necessary to analyze the results by finite size scaling. Normally this is handled by the introduction of a scaling function with a single variable: the ratio of the sample size L to the correlation length ξ = ξ0P−ν,where \(p = \frac{{\left| {p - pc} \right|}}{pc}\) that smoothes over the non-analytic critical behavior.

However, for the present, bimodal mobility distribution problem ES argued that the ratio of conductivities plays the role similar to an external magnetic field in a critical magnetic system: \(\frac{{\sigma_{low} }}{{\sigma_{high} }} = h\), and already smoothes over the critical behavior. Therefore, the model needs to be analyzed by a two variable finite size scaling form, where the ratio L/ξ is a second factor that smoothes the critical behavior:

$$ \mu \left( {P,L,h} \right) \, = P^{ - \alpha } \mu \left( {hP^{ - \Delta } ,LP^{\nu } } \right) $$
(7)

where µ(x,y) is the universal finite size scaling function, and α,ν and ∆ are critical exponents. The analysis is more tractable if the singular P dependence is absorbed by factoring out \(hp^{ - \alpha /\Delta }\) from µ, leaving us with:

$$ \mu \left( {P,L,h} \right) \, = h^{ - \alpha /\Delta } \mu^{0} \left( {hP^{ - \Delta } ,LP^{\nu } } \right) $$
(8)

Since µ0(x,y) is a two-variable function, the full testing of the finite size scaling hypothesis would require a quite extensive computational effort. Therefore, we narrowed our analysis of the finite size scaling assumption to the first variable, hP−∆, while keeping the second variable, LPν, constant. As the lattice sizes were varied from L = 32 to L = 128, keeping LPν constant required the appropriate modification of P. We then chose the h values so that the critical regime on either side of the critical point pc was well sampled.

Figure 5 shows the scaled mobility hα/∆ as a function of h−1P∆ for a fixed value of LPν, for three lattice sizes, varying from L = 32 to L = 128. Reassuringly, we were able to achieve very good data collapse within the (− 0.1,0.1) critical regime around the critical point at P = 0, which remained acceptable out to (− 0.2,0.2). Using the literature values of v = 4/3 and pc = 0.5, the best collapse was reached with exponents α =  − 0.99 ± 0.02 and ∆ = 2.02 ± 0.02.

Figure 5
figure 5

Scaled data of µhα/∆ vs. h−1P∆ for 3 different lattice sizes. The product LPν is held constant.

Full size image

The success of the finite size scaling shows that the Efros Shklovskii analogy to critical spin systems in an external magnetic field is indeed appropriate for this bimodal percolation problem: as the fraction p of the high mobility undefected NP SLs increases, one can think of the evolution of the overall mobility as a modified percolation transition, rounded by the finiteness of the mobility of the low mobility NP SLs. As far as the authors know, this is the first report of the critical exponents and the universal scaling function of the bimodal distribution resistor network problem.

The following points are worth making. Figure 5 shows that the overall network electron mobility, or conductivity, displays a marked transition from a low mobility insulator behavior when the high mobility NP SLs do not percolate yet, to a high mobility insulator behavior once the high mobility NP SLs percolated. Of course, both of these regimes are insulators, so while the geometry of the NP solid undergoes a genuine percolation transition, the conductive properties exhibit only a low mobility insulator-to-high mobility insulator transport crossover, not a genuine phase transition.

We have studied a version of this problem recently on the level of our HINTS code25, where the NPs were connected with either low mobility activated insulating links, or high mobility, non-activated metallic links. In that version of the problem, the underlying percolation transition of the metallic links of the NP solid drove a genuine metal–insulator-transition (MIT), as the percolation of the metallic links created a genuine metallic phase. We conceptualized that MIT as a Quantum Percolation Transition, and adapted the ES bimodal mobility distribution percolation model for its description, as at any fixed temperature that version of the problem also consisted of a bimodal mobility distribution with low mobility links and high mobility links. Whether the high mobility phase is a metal or a high mobility insulator can be identified from the temperature dependence of its conductivity. In the absence of the present detailed finite-size scaling study, in25 we developed a simple, mean-field model form for the scaling function that described the mobility’s evolution from the low mobility insulator to the high mobility metal. With the notation of the present paper, the mean field exponents were α =  − 1 and ∆ = 1. This enabled us to create the dynamical phase diagram of the model on the electron filling—disorder plane, where the MIT separated the insulating phase with activated conductivity from the metallic phase with non-activated conductivity.

The present bimodal TRIDENS study scales up our previous bimodal HINTS work to much larger length scales. The key distinction is that the building blocks of the HINTS network were the individual NPs, whereas in TRIDENS the building blocks are the NP SLs with around a thousand NPs (in the present work, with 800 NPs). Further, in HINTS the origin of the bimodal mobility distribution was the presence or absence of metallic links between individual NPs, whereas here the origin of the bimodal mobility distribution is the presence or absence of an planar defect across the individual NP SLs. Obviously, the HINTS transport modeling that tracks individual electrons transitiong between individual NPs is more detailed than the resistor network of the present TRIDENS work. Nevertheless, since the building blocks of both the bimodal HINTS and the bimodal TRIDENS are low mobility links and high mobility links, we expected that the same ES bimodal percolation model with the same universal scaling function and exponents will capture the critical behavior of the bimodal TRIDENS results. The success of the data collapse with the ES finite size scaling form validated this expectation. While, of course several different analytic forms can be fitted to the universal scaling function that emerged in Fig. 5, nevertheless it was reassuring that in particular the mean-field function of the bimodal HINTS study:

$$ \mu^{0} \left( {hP^{ - \Delta } ,LP^{\nu } \to \, \infty } \right) \, = \, 1/\left( {1 \, + P/h} \right) $$
(9)

was also consistent with it. Further, the α =  − 0.99 exponent of the TRIDENS scaling is approximately equal to the α =  − 1 mean field value within the margin of error. We noted that there was a difference regarding the ∆ exponent: TRIDENS gave a ∆ = 2.02, whereas in the mean field theory ∆ = 1. However, such differences occur typically between mean-field and numerically determined exponents. All in all, the substantial correspondence between our HINTS and TRIDENS works demonstrated that the ES bimodal percolation model is a good, quantitative description of how the underlying percolation transition of the NP solid drives a low-mobility insulator-to-high-mobility insulator transport crossover.

For completeness we mention that we implemented TRIDENS with randomly selecting high or low mobility NP SLs for the links. This corresponds to planar defects with a length of tens of NPs. However, the sample in Fig. 1 has many defects that are much longer. Such long defects can be modelled by selecting defected SLs along lines of links in TRIDENS, in a correlated manner. Such correlated TRIDENS models will be pursued in a future work.

Conclusions

Transport in nanoparticle solids must be simulated on extremely large length scales, corresponding to millions of NPs, because NP solids exhibit spatial structures on several length scales, from the subtleties of individual NPs through the sensitive modeling of inter-NP transitions and through transport across homogeneously disordered SLs all the way to transport in NP SLs with large planar defects. Single-level computational methods are manifestly unable to span these length scales. This is why in our previous work we developed the multi-level HINTS method that was capable of simulating transport across NP solids with up to a thousand NPs. However, even HINTS is unable to capture the effect of planar defects on transport in NP solids of the size of tens of microns.

In this paper, we reported the development of the TRIDENS method that adds three further hierarchial layers on top of the HINTS method. In TRIDENS, we first introduced planar defects into individual NP SLs that comprised the order of about a thousand NPs. Then we used HINTS to simulate the transport across these defected NP SLs. We performed these HINTS transport simulations for tens of thousands of defected NP SLs, and constructed the distribution of the NP SL mobilities with planar defects. Second, the defected NP SLs were assembled into a resistor network with more than 104 NP SLs, thus representing about 107 individual NPs. This translated to length scales of tens of microns, approaching the experimental scales for NP solids. Third, and finally, the TRIDENS results were analyzed by finite size scaling to explore whether the percolation transition, separating the phase where the low-mobility-defected NP SLs percolate from the phase where the high-mobility-undefected NP SLs percolate, drives a low-mobility-insulator-to-high-mobility-insulator transport crossover that can be extrapolated to genuinely macroscopic length scales.

Our extensive TRIDENS analysis generated the following results. On the level of individual NP SLs, we found that the average of the mobility for undefected NP SLs was 0.42 ± 0.1 cm2/Vs, for twin-plane-defected NP SLs 0.16 ± 0.06 cm2/Vs, and for grain-boundary-defected NP SLs 0.09 ± 0.05 cm2/Vs. On average, grain boundary defects hinder transport about twice as much as twin planes. This result makes sense, as grain boundaries are more disruptive to lattice periodicity than twin planes, and transport across the grain boundaries involves longer hops between more distant NPs, whereas transport across twin planes proceeds across many NPs shared by the grains on the two sides of the twin plane, and thus it involves regular hop lengths. It is noteworthy that the introduction of planar defects into NP SLs reduced their mobility by a factor of up to 5. On one hand, this is a substantial suppression of the mobilities that drives a transport crossover, and thus demonstrates the imperative of minimizing the density of planar defects in NP solids to help their suitability for applications. On the other hand, this is not a qualitative, order-ofmagnitude suppression of the transport that indicate a Metal–Insulator-Transition: those are driven by the loss of phase coherence.

On the highest, resistor network-level analysis of TRIDENS, we observed that the introduction of the planar defects immediately started to reduce the network mobility. This finding suggests that even small concentrations of planar defects are not shorted out in NP solids, and thus every reduction of the density of planar defects will lead to further improvements of the transport in NP solids. Among the planar defects, the elimination of grain boundaries pays more dividends than that of twin planes.

For the theoretical description, we adapted the Efros-Shklovskii bimodal mobility distribution percolation model. We performed a finite size scaling analysis of the TRIDENS network mobilities. We demonstrated that increasing the density of the undefected NP SLs drives an underlying, structural/geometric percolation transition in the NP solid, which in turn drives a low-mobility-insulator-to-high-mobility-insulator transport crossover. We demonstrated that our adaptation of the ES bimodal theory’s two-variable scaling function is an effective tool to quantitatively characterize this low-mobility-insulator-to-high-mobility-insulator transport crossover. For context, we discussed the analogies with the Quantum Percolation Transition we developed in our earlier, MIT-focused work25.