Disentangling morphology and conductance in amorphous graphene

Department of Physics, McGill University, Montréal, Québec, Canada

Understanding the role of microscopic structure in determining the macroscopic properties of a material is an important overarching goal in chemistry, physics, and molecular engineering. Of particular interest is establishing structure-function relationships for materials classified as “amorphous”. Although the adjective “amorphous” may suggest that the structure is simply a random arrangement of atoms, it is well known that the properties of real amorphous materials sensitively depend on the preparation route (deposition rates, temperature of substrate, etc.) implying configurational tunability at the atomic-level. The rich chemical and conformational landscapes of these materials pose many challenges to theory and experiment, and drive innovation in both ¹.

Disorder inherent in amorphous materials is difficult to control and characterize experimentally and it is hard to replicate in a simulation ¹. Amorphous graphene, or amorphous monolayer carbon (AMC) has recently emerged as a uniquely suitable model system for studying structure-function relationships in amorphous materials ². Thanks to the two-dimensional configuration which can be visualized using microscopy and to the conjugated-carbon composition that is relatively easy to model computationally the correspondence between simulation and experiment can be established precisely. In recent years AMCs have garnered significant research interest regarding the effect of disorder on the electronic ^{3, 4, 5}, thermal ^{6, 7}, structural ^{8, 9, 10}, and mechanical properties ¹¹. However, it is the novel synthetic protocol by Tian et al. that turned AMCs into a class of amorphous materials resolved along a well characterized morphological axis ¹².

Our focus here is on electrical conductance in AMCs, a property which was shown to have remarkable sensitivity to morphology. It was reported that a small change in synthetic conditions of AMC films led to a dramatic billion-fold increase in the sheet conductance ¹². Building a mathematical relationship between morphology and conductance is a way to gain insight into this instability. From the computational perspective, however, modeling charge transport in such systems is extremely challenging because precise atomistic AMC structures need to be analyzed at the mesoscale, the scale of experimental measurements. Since AMC is a covalently-bonded network of conjugated carbons and cannot be meaningfully fragmented, a good model must generate and analyze high quality atomistically resolved mesoscale samples. By anchoring computational results in experimental data, this paper aims to unravel the morphology-conductance relationship in AMCs by using machine-learning-augmented simulation methodologies ^{13, 14, 15}, thus demonstrating that computational modeling can be a helpful tool in amorphous materials research.

For the purposes of our analysis we generate three ensembles of morphologically-distinct mesoscale AMCs: two to be compared directly with experimentally reported data, and one to be used to demonstrate an ambiguity in the currently accepted AMC morphology classification method. Samples of the three morphological classes are shown in Figure 1. Each mesoscopic AMC has the dimensions of roughly 40 nm $\times$ 40 nm and was generated using a sampling approach based on generative machine learning: the Morphological Autoregressive Protocol (MAP) ^{13, 14}. This approach takes advantage of decaying structural correlations characteristic of amorphous materials to sample large-scale conformations by extrapolating from smaller (order of correlation length) samples. The training samples in this work were produced by a Monte-Carlo bond switching protocol driven by ML energies from the C-GAP-17 potential ¹⁵. In order reduce the error bars on computational predictions to acceptable levels more than 200 samples were generated in each ensemble. For full technical details on dataset generation and the MAP implementation, see Section S1 of Supporting Information.

To classify the morphology we use the same descriptors as Tian et al.: (i) the pair correlation function $g(r)$ averaged over each ensemble (see Figure 2a); (ii) the short- and medium-range order parameters, $\rho_{\text{sites}}$ and $\log\eta_{\text{MRO}}$ defined below, (iii) the distribution of different types of carbon rings, and (iv) the area distribution of crystalline inclusions. The order parameter $\rho_{\text{sites}}$ is the number of connected clusters of undistorted carbon atoms per unit area in a given structure. A carbon atom is considered undistorted if the bond lengths and bond angles it forms with its immediate neighbours deviate by under 10% from those in pristine graphene. This order parameter therefore measures the persistence of short-range order in a AMC sample. The medium-range order parameter $\log\eta_{\text{MRO}}$ order parameter is related to the integral of the pair correlation function in the medium range distance between 4 Å  and 12 Å. For technical details on all four descriptors, see Section S2 of Supporting Information. A discussion of our ensembles’ bond length and bond angle distribution can also be found in Section S3 therein.

We start our morphology analysis with the most disordered ensemble, see Figure 1a for sample visualization, and the green line in Figure 2a for evidence of fast-decaying pair correlations (note that unlike in Ref 12, we normalise our pair correlation functions such that $g(r)\xrightarrow[r\rightarrow\infty]{}1$). As shown in Figure 2b, the green set of points corresponding to this ensemble, it is positioned in the upper left corner of the $\log\eta_{\text{MRO}}$ vs. $\rho_{\text{sites}}$ plot – close to the position of AMC-500 ensemble reported by Tian et al. ¹², and its distribution of rings shown Figure 2c (in green) is close to the experimental AMC-500 ensemble as well. Due to its similarity in morphological metrics to the AMC-500 from Ref. 12, we will refer to it as “simulated AMC-500”, or sAMC-500 henceforth.

The next ensemble we discuss is the most grahene-like ensemble shown in Figure 1c. It is much more ordered as it contains large sections of almost perfect hexagonal order (graphene-like hexagonal motifs are highlighted in lime green in Figures 1a-c). This ensemble displays ring statistics that are close to the AMC-300 samples in Tian et al. ¹² and it is found in a similar region of the $\log\eta_{\text{MRO}}$ vs. $\rho_{\text{sites}}$ space. We will therefore refer to this ensemble as “simulated AMC-300”, or sAMC-300.

The third ensemble, visualized in Figure 1b, is somewhat difficult to classify. Similarly to Tian et al.’s AMC-400, our third set of structures exhibits a degree-of-disorder greater than that of the sAMC-300 samples, but lower than that of those in the sAMC-500 as inferred from observations of disorder in Figures 1a-c, and the intermediate decay of correlations in Figure 2a. It is found in a position similar to AMC-400 in $(\log\eta_{\text{MRO}},\rho_{\text{sites}})$ space. However, the ring statistics (Figure 2c in purple) show a balance of crystalline hexagons (6-c) vs. non-crystalline hexagons (6-i) that is reversed with respect to the experimentally generated AMC-400 samples (see Figure 2h in Ref. 12): in our simulated structures, the crystalline hexagons 6-c are much more common than non-crystalline 6-i hexagons, while the opposite is true in the experimental samples. Because of this mismatch, we will refer to this simulated ensemble as sAMC-q400 (for “quasi-400”). This ensemble highlights the possibility of non-unique classification of morphologies based on position in the ($\log\eta_{\text{MRO}}$, $\rho_{\text{sites}}$) space.

It is widely accepted that the lack of periodicity in amorphous materials leads to spatially localized electronic states ^{16, 17, 18}. This phenomenon prevents charge carriers from travelling coherently across the material. Instead, charge transport in amorphous phases is modelled within a variable range hopping (VRH) picture, wherein charge carriers incoherently travel between localised sites by exchanging energy with a bath of molecular vibrations ¹⁹. While the billionfold enhancement in conductance displayed by the high-crystallinity AMC samples in reference 12 might lead one to suspect that a coherent band transport mechanism may dominate, transport measurements show that their resistance-temperature relation obeys the two-dimensional Mott law^{12, 19, 2}, which is typical of VRH in systems with a slowly varying density of states (DOS) ^{19, 20, 21}. This is strong evidence that the VRH picture remains valid, even in for such ordered AMCs.

The hopping sites in the VRH picture are usually taken to be spatially localised eigenstates of the electronic hamiltonian ²¹. Since AMCs are $sp^{2}$-hybridized conjugated carbon systems ², it may be safely assumed that only the $\pi$-network electrons contribute to conduction. We therefore model the electronic structure of mesoscopic AMC fragments using an all-atom tight-binding hamiltonian $\mathcal{H}=\sum_{\langle i,j\rangle}t_{ij}(|\varphi_{i}\rangle\langle\varphi_{j}|+|\varphi_{j}\rangle\langle\varphi_{i}|)$, where $|\varphi_{i}\rangle$ denotes the $2p_{z}$ orbital centered on the $i$th carbon in the AMC sample. The sum is carried out over all nearest-neighbour pairs $\langle i,j\rangle$, and the semi-empirical parametrisation of the hopping elements $t_{ij}$ is adapted from prior works ^{22, 23}, see Section S4 of Supporting Information for details.

We work with 40 nm$\,\times\,$40 nm atomistically resolved samples, and thus do not resort to coarse-graining, applying periodic boundary conditions, modeling only partially resolved structures, or any other kind of commonly employed simplifications. While the enormous size of the hamiltonian matrices and the large number of samples in our ensembles preclude us from resolving the full eigenspectrum using standard numerical routines, doing so is not necessary for this problem. Since thermally activated conductance is expected to be mediated by a relatively narrow band of thermally accessible states, we partially diagonalise each hamiltonian using the Lanczos algorithm ²⁴, which is well-suited to the tight-binding hamiltonians’ sparse structure (Section S5 of Supporting Information). We thereby obtain the subset of molecular orbitals (MOs) whose energies lie within $4k_{\text{B}}T$ above the chemical potential $\mu$ at room temperature ($T=300$ K), under different gating conditions. We focus on three regimes: (i) the ungated regime in which $\mu=\epsilon_{F}$, the sample’s Fermi level at half-filling; (ii) the regime where a strong negative gate voltage is applied and $\mu=\epsilon_{0}$, the smallest eigenvalue of $\mathcal{H}$; and (iii) the regime where a strong positive gate voltage is applied and $\mu=\epsilon_{N}-4k_{\text{B}}T\,$ where $\epsilon_{N}$ is the greatest eigenvalue of $\mathcal{H}$.

It is expected that solutions obtained from such a low level of theory may possess some degree of artefacts. In this case, we noticed that some of the eigenstates we obtain are delocalised in ways that we do not expect to withstand the effects of decoherence and localization due to factors that were not included in our calculation but are certainly present in the real system. Some of the factors that are not included in our model are the electron-electron and electron-phonon interactions, as well as the interactions with substrate, and distortions such as ripples or buckling that AMC sheets are expected to undergo at ambient conditions ^{25, 26}, which would disrupt the extent of $sp^{2}$ conjugation and contribute to electronic localisation ²⁷. In our case, the delocalized eigenstates/molecular orbitals (MOs) which we regarded as pathological typically have several disjoint pockets of high electronic density separated by distances that can span 10’s of nanometres – almost the entire AMC structure’s length (see Figs. S3 and S4 in Supporting Information). This behavior was more common among MOs found near $\epsilon_{F}$, whereas the eigenstates at the band edges expectedly tended to be localised, with very few exceptions. In order to model VRH of charges we have made the choice to construct conduction space in which we artificially localize the pathological MOs into multiple disjoint charge-hopping sites.

To construct the state space for the VRH network, we have developed a procedure based on $k$-means clustering to extract hopping sites from the MOs supported by each AMC fragment, (see section S7 of Supporting Information). This approach partitions each (n^th) MO $|\psi_{n}\rangle$ into a set of $m_{n}$ localised states $\{|s_{n,i}\rangle\}_{i=1,\ldots,m_{n}}$, from which site positions $\{\mathbf{R}_{n,i}\}_{i=1,\ldots,m_{n}}$ and delocalisation radii $\{a_{n,i}\}_{i=1\ldots,m_{n}}$ can be extracted. MOs that are already localised are left unchanged by this procedure. As the last correction mechanism for numerical artifacts in the electronic structure modeling, the occasional sites whose effective area $\pi a_{n,i}^{2}$ exceeded the area of the largest crystalline inclusion found in the ensemble (Fig. 2d) were considered numerically artificial and removed from our calculation. We tabulate the maximum allowed site radii $a_{max}$ for the three AMC ensembles in Table 1. We note that the experimental crystallite sizes reported in reference 12 are much smaller than the ones we obtained from our simulated structures. This may be because the AMC microscopy images processed by Tian et al are much smaller ($\sim$ 5 nm $\times$ 5 nm) than the sAMC samples we use in our analysis ($\sim$ 40 nm $\times$ 40 nm).

Having constructed the VRH space, we apply percolation theory to estimate AMC conductances. Percolation theory has been very successfully applied as an analytical framework to estimate the VRH conductance of various classes of disordered semiconductors whose density of states (DOS) has a simple closed-form expression (e.g. a Gaussian distribution) ^{20, 28, 29, 30, 31, 32, 33, 34, 35, 36}. Practically, percolation-based approaches are attractive because they do not suffer from stability or convergence issues associated with the usual techniques of VRH simulations like the solution of a transport master equation ^{37, 38}, or Monte Carlo sampling of hopping trajectories ^{39, 40}. Seeing as the DOS profiles for the three morphological classes we consider are not known analytically, we developed a numerical implementation of percolation of theory, which yields an ensemble-averaged estimate for charge conductance.

We define the hopping rate $\omega_{ij}$ between sites $i$ and $j$ using a Miller-Abrahams expression which we modified to accommodate for variable site radii: $\omega_{ij}\sim\,e^{-\xi_{ij}}$, where the dimensionless quantity $\xi_{ij}$ can be thought of as an effective distance between sites $i$ and $j$ and depends on the sites’ positions, energies, and radii (see equation (S5) for $\omega_{ij}$ in Supporting Information). Each structure admits a critical distance $\xi_{c}$ at which a cluster of sites obeying $\xi_{ij}\leq\xi_{c}$ percolates the sample by connecting its right edge to its left edge. The randomness inherent to the AMC fragments, as well as their finite size, will lead $\xi_{c}$ to fluctuate from fragment to fragment. Following an approach similar to Rodin and Fogler ³⁵, we estimate each ensemble’s conductance $\sigma$ as follows:

where $P(\xi)$ is the probability of having a percolating cluster through sites obeying $\xi_{ij}\leq\xi$ in a given ensemble, $\omega_{0}=1$ fs^-1 is the escape frequency, $q_{e}$ is the elementary charge, $k_{\text{B}}$ is Boltzmann’s constant, and $T=300$ K denotes temperature. See Section S8 of Supporting Information for technical details.

Next, we summarize and discuss the results of our modeling. First, we reproduce experimentally observed conductances reported in Ref. 12 for sAMC-500 and sAMC-300 ensembles. We then discuss the deviation of the predicted conductance for the sAMC-q400 ensemble relative to the experimentally characterized in Ref. 12 AMC-400. Our results imply that a unique map between the $(\log\eta_{\text{MRO}},\rho_{\text{sites}})$ space and conductance, in contradiction to previous assertions (e.g., Figure 4a in Ref. 12), does not exist. We conclude with the discussion of the possibility of controlling conductance in AMCs by applying gate voltage and in effect modifying the characteristics of the charge transport pathways while keeping the AMC sample unchanged.

We performed our conductance calculations at $T=300\,$K, under different gating conditions. The results are plotted in Figure 3a, and a summary is tabulated in Table 1. First, we focus on the conductance in the ungated case in which charge transport is carried by states close to the Fermi energy ($\epsilon_{F}$) of the AMC, the middle set of data-points in Figure 3a. The two ensembles that we may directly compare to experimental data are sAMC-300 shown in orange and sAMC-500 shown in green. We observe a conductance gap of $5$ orders of magnitude between the two in our calculations which matches well the gap observed in Ref.12, see Figure 3d therein¹¹1Note that all data plotted Figure 3d of Ref. 12 were collected at room temperature ⁴¹. Its abscissa thus corresponds to the substrate temperature during the AMC films’ growth, which controls their degree of disorder, rather than the temperature of the films during the resistivity measurements..

Of particular interest is the high conductance of the sAMC-q400 ensemble which is found to be on par with the high, relative to other AMCs, conductance of sAMC-300. The sAMC-q400 ensemble falls close to the experimental AMC-400 ensemble in the $(\log\eta_{\text{MRO}},\rho_{\text{sites}})$ space with the center close to (-1.3,0.3) point compared to approximately (-1.25,0.35) in the experimental case. However, AMC-400 was experimentally and theoretically shown to be a perfect insulator. Previous modeling approaches assign a smooth map from the $(\log\eta_{\text{MRO}},\rho_{\text{sites}})$ space to conductance with vanishing values in this particular region of the $(\log\eta_{\text{MRO}},\rho_{\text{sites}})$ space, and would therefore struggle to reconcile the sharp difference between the behaviors of AMC-400 and sAMC-q400. Our result, on the other hand, demonstrates that the morphology-conductance relationship in AMC is more complex than previously thought. The morphological difference between sAMC-q400 and AMC-400 is visible in ring-distribution statistics (Figure 2c).

Finally, we discuss the possibility to weakly decouple the conductance from morphology in AMCs. The idea behind this is that molecular orbitals in different regions of the spectrum tend to have differing morphological characteristics. We previously discussed the emergence of edge-states around the middle of the spectrum and bulk-localization of MOs towards the edges of the “band” in amorphous graphene nanoflakes ⁴². Similar behavior is observed in mesoscale samples of AMCs as well, and by applying gate voltage we, in effect, modify morphological characteristics of the conducting VRH networks while keeping the overall atomistic morphology fixed.

We note that the effects of applying a gate voltage in our simplified modeling reflect only the characteristics of the MOs that carry the current under the different conditions - the changes to resistance of contacts for example are neglected. Under this idealized assumption, we observe that conductance depends rather strongly on gating and that more disordered AMCs show larger sensitivity. To understand the origin of this effect we need to quantify the morphological changes, i.e., the structural metamorphosis, that conducting networks undergo when a gate voltage is applied.

The connection between the morphology and conductance can be clarified by focusing on the structure of only those regions of AMC samples on which the electronic states that dominate charge transport are localized. To do so we define the crystallinity $\chi$ of a given VRH state $|\psi\rangle$ as follows:

where $\mathcal{C}$ corresponds to the set of crystalline atoms – i.e. atoms belonging to a crystalline hexagon (previously referred to as 6-c). In words, $\chi(|\psi\rangle)$ corresponds to the aggregated density of $|\psi\rangle$ which lies on the crystalline regions of a given AMC sample.

Figures 3b-d show the crystallinity of each VRH site belonging to a percolating cluster against the fraction of crystalline atoms $\phi_{c}$ in the relevant AMC sample. Focusing first on the gated regimes depicted in Figures 3b-c, we find that conducting sites at the band edges (i.e. far from $\epsilon_{F}$), tend to be localized preferentially on crystalline inclusions since they tend to lie above the diagonal and thereby satisfy $\chi>\phi_{c}$, meaning that they are more crystalline than the rest of the structure. We infer from this that charge transport in these regimes is mostly carried by hopping from one crystalline site to another, giving a partial rationalization to an intuitive and common modeling assumption.

The picture changes when we move to the ungated case shown in Figure 3d. There, we find that unlike the states at the band’s edges, the mid-spectrum eigenstates produce VRH sites which predominantly satisfy $\chi<\phi_{c}$ (Fig.  3d), i.e., that are preferentially found on defects in all three ensembles. Charge hopping in the ungated regime therefore takes place predominantly over the disordered regions of AMC. Within the framework of VRH, the bias of the conduction networks towards disordered regions puts the propensity of defects in competition with the localization of the electronic states. The interplay between these two properties gives rise to the predicted conduction trends, and originates in the dependence of the inter-site hopping rates $\omega_{ij}$ on the sites’ radii and site-site distances (equations (S5) and (S6) in Supporting Information). For instance, our calculations show that sAMC-q400 exhibits ungated conductance close to sAMC-300 in spite of being much more disordered overall. This happens because it strikes a balance between the two competing effects: it contains more extensive disordered regions to increase the density of hopping sites relative to sAMC-300, while also retaining enough structural order to give rise to sites with reasonably large radii.

Figure 3 also highlights the contrast between the crystallinity distributions of conducting sites in the different sAMC ensembles, especially in the strong gating regimes (Figs. 3b-c). In this respect, the ordered ensembles sAMC-300 (orange points in Figs. 3b-c) and sAMC-q400 (purple) are qualitatively very similar: their low- and high-energy conducting sites are clustered near the top of Figures 3b-c – indicating their highly crystalline character – while the extremal energy sites in sAMC-500 (green) follow a much more uniform crystallinity distribution. We discuss such qualitative differences in the crystallinity distributions of our three ensembles’ conducting sites in greater detail in Section S9 of the Supporting Information, but comment on them briefly here. First, we note that the crystallinity of a structure’s conducting sites is strongly influenced by the crystallinity of the energy eigenstates from which they originate (see Figure S5 in the Supporting Information). Starting with sAMC-500, its structures are akin to continuous random networks (CRNs), which are statistically homogeneous ⁴³. The MOs supported by sAMC-500 structures will therefore tend to evenly sample a diverse set local atomic environments ^{44, 45}. These MOs therefore exhibit a correspondingly diverse set of crystallinity values (cf. Fig. S5 in the Supporting Information), which in turn leads to the more uniform nature of the crystallinity distributions of their corresponding hopping sites (Figs 3b-c). Conversely, sAMC-300 and sAMC-q400 structures feature extended crystalline domains, and universally obey $\phi_{c}>0.5$. These two ensembles are therefore much more similar to defected graphene than to a CRN. Defected semiconductors are known to exhibit defect-localised states near $\epsilon_{F}$, while the states deep within the occupied and virtual manifolds will retain a highly crystalline character ^{21, 46}. This is entirely consistent with what we observe in Figure S5, and therefore explain the very high crystallinity of their hopping sites.

Under our simplifying assumptions, the gate-voltage modulation charge transport in AMCs is owed to the metamorphosis of the conducting states along the electronic energy spectrum. This ready and reversible tunability of AMC’s conductance – taken together with the low heat conductance which is owed to its inherently disordered bonding network ⁷ – make is potentially an attractive candidate for thermoelectric applications. This observation is likely to remain at the conceptual level until effective strategies are invented to increase conductivities in AMCs to levels sought after in thermoelectric materials, i.e. $\sim 10$ S cm^-1 47. In the past, nitrogen doping has been shown enhance the electrical conductance of AMC by an order of magnitude ⁴⁸ and sets a precedent for future developments in this space.

In this work we have combined deep learning-enhanced simulation techniques with percolation theory to model charge conductances in three morphologically distinct mesoscale AMCs. We have overcome the challenges of modeling electronic conductance in mesoscale atomistically resolved covalently bonded networks of conjugated carbons by developing a custom partial diagonalisation procedure based on the Lanczos algorithm and adapted the percolation theory calculation of charge conductance to the pecularities of the AMC system. Our protocol is noteworthy in that is avoids the artifacts that may arise from applying periodic boundary conditions to aperiodic amorphous structures (a common practice), or from oversimplifying assumptions regarding transport mechanism. We reproduce the reported dependence of charge conductance on morphology and discuss the ambiguous relationship between incomplete/partial measures of morphology and conductance. By conducting a crystallinity analysis of the conducting sites we show that they metamorphose in response to gate voltage from being localized on crystallites at band edges to being localized on defects around the Fermi energy. Inspired by this observation, we explore the possibility to decouple the AMC morphology from electronic conductance by applying a gate voltage and comment on the potential of thermoelectric tunability in amorphous conductors.

The following files are available free of charge.

• •

  Supporting Information: Detailed description of our methods for (i) modeling AMCs using MAP, (ii) characterising AMC morphologies, (iii) electronic structure calculations, (iv) defining the VRH sites, and (v) computing AMC conductance from numerical percolation theory; Discussion of the bond length and bond angle distributions in our simulated structures; Comparison of our method for diagonalising large tight-binding Hamiltonians with an exact benchmark; Discussion of the crystallinity distributions of tight-binding eigenstates in different of our sAMC structures’ energy spectrum. Repositories containing the code used in this work are also included. (PDF).

We thank Zakariya El-Machachi and Volker Deringer for sharing amorphous graphene dataset they published in Ref.15 and Alessandro Troisi for helpful tips on convergence issues in charge conductance calculations. Funding from NSERC Discovery Grant, IVADO, and AI4Design program of the National Research Council of Canada is gratefully acknowledged. Computations were performed on Calcul Quebec supercomputers (Narval) funded by Canada Foundation for Innovation (CFI).