Influence of nuclear quantum effects on the electronic properties of amorphous carbon.

We carry out quantum simulations to study the physical properties of diamond-like amorphous carbon by coupling first-principles molecular dynamics with a quantum thermostat, and we analyze multiple samples representative of different defective sites present in the disordered network. We show that quantum vibronic coupling is critical in determining the electronic properties of the system, in particular its electronic and mobility gaps, while it has a moderate influence on the structural properties. We find that despite localized electronic states near the Fermi level, the quantum nature of the nuclear motion leads to a renormalization of the electronic gap surprisingly similar to that found in crystalline diamond. We also discuss the notable influence of nuclear quantum effects on band-like and variable-hopping mechanisms contributing to electrical conduction. Our calculations indicate that methods often used to evaluate electron-phonon coupling in ordered solids are inaccurate to study the electronic and transport properties of amorphous semiconductors composed of light atoms.

Amorphous semiconductors are interesting materials for countless applications including thermoelectric (1, 2), solid-state memory (3), and photovoltaic devices (4, 5) and for the development of optical fiber technologies (6). In particular, diamond-like amorphous films (DLCs) (7–9), i.e., amorphous carbon materials with densities close to that of crystalline diamond, have attracted widespread attention in the last few decades, because they are easy to synthesize, environmentally friendly, and mechanically hard. In addition, they exhibit tunable band gaps, which make them promising components of thin-film transistors, solar cells, and microelectromechanical devices (5, 9). However, the presence of defect-like states and low carrier mobilities remain challenging factors for optimizing the performance of electronic devices (5, 9) based on amorphous carbon.

Diamond-like films have been extensively investigated experimentally, using photoelectron spectroscopy (10, 11) and infrared (IR) (12, 13), Raman (11, 13, 14), and near-edge X-ray absorption fine-structure (NEXAFS) (13) spectroscopies, and, in most cases, the interpretation of experiments relies on the availability of atomistic structural models accounting for the complexity of disorder (15). Since the development of first-principles molecular dynamics (FPMD) (16), much progress has been reported in deriving structural models of amorphous carbon and DLCs from computer simulations based on first principles (17–22), as well as using empirical potentials (23–29), with a focus on structural and vibrational properties and, to a lesser extent, electronic properties (15). However, none of these studies has so far accounted for the quantum motion of the nuclei.

Nuclear quantum effects (NQEs) are known to be significant in determining the electronic properties of crystalline diamond and molecular forms of carbon (30) and are expected to play a key role also in the case of DLCs. Based on previous heat-capacity measurements (31), our estimate for the Debye temperature of DLC films is 1,400 K, which indicates that NQEs may play an important role in determining the electronic properties of DLCs, even at room temperature. Recently, using first-principles molecular dynamics simulations coupled with a quantum thermostat (32), we investigated the effect of quantum vibronic coupling on the electronic gaps of carbon systems and indeed found a substantial impact not only on the band gap of diamond but also on a representative structural model of disordered diamond-like films (30).

Here we present detailed predictions of the influence of nuclear quantum effects on the physical properties of diamond-like films without the presence of any element lighter than carbon (e.g., hydrogen). We investigate several structural models representative of various defective structures present in the amorphous network (21). We show that quantum vibronic coupling substantially affects the electronic properties of DLCs, in particular, electronic and mobility gaps; however, it has a moderate influence on the structural properties of the films. Our results point to the critical role of nuclear quantum effects in determining the electronic transport properties and ultimately the electrical conductivity of DLCs, even when hydrogen is absent.

In our simulations we included NQEs by coupling FPMD with a quantum thermostat, following the same protocol recently adopted to study crystalline diamond (30). We refer to these calculations as “quantum” simulations. We also performed FPMD simulations with a stochastic velocity rescaling thermostat (33), hence neglecting NQEs, and we refer to these calculations as “classical” simulations.

Before describing our results, we discuss the choice of our simulation protocol and its advantages. Within the Born–Oppenheimer approximation, several methods are available in the literature for the study of vibronic coupling and the calculation of electron–phonon interaction in solids and molecules from first principles; these methods include the many-body perturbation theory (MBPT) (34) and the frozen-phonon method (35), which incorporates NQEs using the harmonic approximation. The latter, however, is not justified when describing amorphous materials (36) that usually exhibit anharmonic potential energy surfaces and strong phonon–phonon coupling. Anharmonic potenial energy surfaces may be conveniently sampled using molecular dynamics (MD) simulations, and several atomistic models of amorphous carbon have been derived using first-principles MD (17–22) and snapshots utilized to investigate the electronic properties of the system. However, as mentioned above, FPMD does not include nuclear quantum effects, and, in addition, it describes the distribution of vibrational states in an approximate, classical manner (37), instead of including the proper quantum (Bose–Einstein) statistics. Hence, FPMD is expected to yield accurate results only at temperatures close to or higher than the Debye temperature. Below the Debye T one could use path-integral FPMD simulations (38–40), which constitute an accurate framework to describe NQEs; however, they are computationally rather demanding and they would be prohibitive, from a computational standpoint, for the study of an amorphous solid for which the sampling of several configurations is required. In contrast, a quantum thermostat based on a colored noise generalized Langevin equation is computationally much less demanding than straightforward path integral simulations. Hence, we chose to use a quantum thermostat whose performance against path integral FPMD simulations was previously tested in the case of crystalline diamond (30). We now turn to the descriptions of our results.

Results

Structural Properties.

We carried out quantum and classical simulations starting from seven selected configurations generated in ref. 21 using FPMD and cells containing 216 C atoms. We computed the electronic gap between highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO). In the absence of any electron–phonon interaction, i.e., when the positions of the C atoms are at rest, these samples exhibit a wide variety of static HOMO-LUMO gaps ranging from 0.84 to 3.63 eV when the Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional (41, 42) is used in density functional theory (DFT) calculations (Table 1). Hence, they represent an adequate set of specimens to investigate the electronic properties out of the chemical configurational space of amorphous carbon at a density close to that of the crystal. Although the PBE functional is expected to underestimate electronic gaps, it is a useful tool to study trends in various samples and to discuss the effect of the nuclear motion on the properties of DLCs.

Table 1.

Summary of HOMO-LUMO, renormalization values (ZPRs), and mobility gaps computed for several diamond-like amorphous carbon samples of density 3.25 g/cm3

	Static gap, eV		ZPR, eV
	HOMO-LUMO	Mobility	HOMO-LUMO	Mobility
1)	0.844	4.274	–0.093	–0.228
2)	1.203	4.746	–0.100	–0.473
3)	1.930	4.411	–0.155	–0.309
4)	2.511	4.600	–0.234	–0.407
5)	2.786	4.359	–0.238	–0.371
6)	3.149	4.488	–0.370	–0.375
7)	3.630	4.920	–0.396	–0.803

ZPR values are computed by extrapolating the quantum simulation results using the Viña model fit (46).

We start by considering the structural and bonding properties of the various samples. We computed pair correlation functions (PCFs) between C atoms and maximally localized Wannier function centers (43) (denoted by X) for crystalline diamond as well as amorphous C configurations obtained from classical and quantum FPMD simulations at 100 and 500 K (Fig. 1). The positions of Wannier function centers give insight into bonding charges in the system. Following previous studies (44, 45), we define two C sites as bonded if they share a Wannier function center (X) within a given cutoff distance, chosen to correspond to the first minimum of the C–X pair correlation function of diamond.

Fig. 1.

(A and B) The pair correlation function between C atoms and Wannier function centers (X) computed for crystalline diamond with a cell with 216 atoms and for seven a-C samples studied using classical FPMD (C) and FPMD coupled with a quantum thermostat (Q) at 100 and 500 K, respectively. (C and D) The percentage of threefold-coordinated C atoms as a function of the static HOMO-LUMO gap (E) for the same seven a-C samples at 100 and 500 K, respectively.

As expected, when considering NQEs, the distribution of bonds in both crystalline and amorphous C is substantially broadened, with classical and quantum distributions exhibiting more marked differences in the case of the crystal.

The samples considered here contain threefold-coordinated (3c) and fourfold-coordinated (4c) sites; the latter are sp3 bonded, while the former can be either sp2 or sp3 bonded, depending on the angles between bonds. In other words, 3c sites are not always sharing angles of 120∘, and their first neighbors do not belong to the same plane, as in graphite; rather, the great majority of the 3c sites are diamond-like sites with a dangling bond. However, a low concentration of sp2 sites is present in the sample with the lowest HOMO-LUMO gap. The percentage of 3c sites obtained from classical and quantum simulations at 100 and 500 K is shown in Fig. 1, respectively, as a function of the static HOMO-LUMO gap, which uniquely defines a sample. The fraction of 3c atoms does not appear to be correlated with the static HOMO-LUMO gap. We find that the inclusion of NQEs increases the fraction of 3c sites by a relatively small amount, at most at 100 K, and as expected, the fraction of 3c sites slightly increases at higher T, while the effect of nuclear quantum motion decreases. We conclude that the overall influence of NQEs on the structural properties of DLCs is rather weak. Such an influence is instead important on the electronic properties of the system, as we show below.

Electronic Properties.

Fig. 2 shows the electronic density of states (EDOS) at 100 K, computed using 700 snapshots obtained from classical and quantum simulations for three representative amorphous carbon (a-C) samples with static HOMO-LUMO gaps of 0.84, 1.93, and 3.63 eV, respectively. (For comparison, the fundamental gap of diamond computed at the same level of theory is 4.19 eV.) The results for four additional samples are reported in . We find that the electronic gap decreases due to NQEs, with the valence band maximum (VBM) (or HOMO) being higher and the conduction band minimum (CBM) (or LUMO) lower in energy than in classical simulations. The position of defect states in Fig. 2 is instead largely unaffected by NQEs, although their distribution is broadened.

Fig. 2.

EDOS at 100 K computed with classical FPMD (Classical) and FPMD coupled with a quantum thermostat (Quantum) for three representative samples with static HOMO-LUMO gap of (A) 0.84 eV, (B) 1.93 eV, and (C) 3.63 eV. Shaded regions represent occupied states (valence band) and unshaded regions represent unoccupied states (conduction bands).

The electronic gap renormalization due to electron–phonon interaction for the same three samples in Fig. 2 is shown in Fig. 3 (see for four additional samples). We find that quantum vibronic coupling leads to a renormalization of 80 to 400 meV even at 0 K (corresponding to 8 to 12% of the static gap) and that the results of classical simulations underestimate the gap renormalization obtained with quantum simulations: for example, by 100 to 350 meV at ∼300 K.

Fig. 3.

(A–C) Values of the HOMO-LUMO gap renormalization due to electron–phonon interaction as a function of temperature for three a-C samples with HOMO-LUMO gaps of (A) 0.84 eV, (B) 1.93 eV, and (C) 3.63 eV, respectively. The orange and blue symbols are the results obtained with first-principles molecular dynamics simulations performed with 1) a classical and 2) a quantum thermostat, respectively. The blue (orange) lines are obtained from a Viña (linear) model fit (46) of the quantum (classical) results, while the red lines are the results of FPH calculations. (D) Absolute value of the ZPR of the HOMO-LUMO gap of the seven a-C samples as a function of the static HOMO-LUMO gap (Table 1). Dashed blue and red lines are linear fit of the quantum simulations and FPH results, respectively (Table 1).

Fig. 3 also compares the results of quantum simulations for the zero-point renormalization (ZPR) with those of the frozen-phonon harmonic (FPH) method, which includes NQEs, although only within the harmonic approximation. We find relatively good agreement for samples with a low concentration of defect states and large discrepancies for the other samples, showing that overall, the frozen-phonon results are not robust.

To understand, in detail, the origin of the differences between FPH and quantum simulation results, we computed the overlap integrals between the electronic orbitals of a given snapshot extracted from FPMD simulations and those of their closest local minimum determined by a structural optimization at 0 K. We chose all states within a given energy window (2 eV) from the Fermi level (corresponding to 700 and 500 snapshots at 100 and 500 K, respectively) and obtained an ensemble average of all overlap matrices at each temperature (for the results of quantum trajectories obtained at 100 K, see for the same samples shown in Fig. 3). We found that for configurations with significant mixing between orbitals due to vibronic coupling, FPH and quantum simulation results show poor agreement, due to higher-order electron–phonon coupling effects (i.e., the electronic orbital energy is not quadratic as a function of ionic coordinates). There are instead samples with minimal orbital mixing and for which we find good agreement between FPH and quantum simulation results (Fig. 3). However, the absence of orbital mixing is a necessary but not sufficient condition to ensure that FPH and quantum simulations agree. In fact, we identified one sample () that does not exhibit any significant vibronic mixing between HOMO (LUMO) and other occupied (empty) orbitals, and for which FPH and quantum simulation results show poor agreement. In this case, the discrepancy arises from the presence of strong phonon–phonon scattering affecting the HOMO-LUMO gap of the system. Hence, there are two main physical origins for the inadequacy of a FPH treatment of vibronic coupling in DLCs: the neglect of higher orders in the description of electron–phonon coupling and/or the inaccuracy in describing strong phonon–phonon couplings. Interestingly, for all samples considered here both quantum simulations and FPH calculations show an increase of the ZPR as a function of increasing static gap, although the renormalized values are different in the two cases.

Orbital Localization and Mobility Gap.

We now turn to analyze the impact of NQEs on the localization of the orbitals near the Fermi level and on the mobility gap (47, 48). As an illustrative example, we focus on the sample with the lowest HOMO-LUMO gap, which shows a number of defect-like states within the mobility gap. In Fig. 4, we show the joint probability distribution of electronic energy and inverse participation ratio (IPR) for the orbitals near the Fermi level; the IPR of the ith orbital is calculated as , with a higher IPR value indicating a more localized single-particle wave function. The mobility gap is computed as the energy difference between occupied and empty states with an IPR less than a certain threshold (49), chosen to be 0.0012. At this value we found a sharp transition of the IPR as a function of energy (). We emphasize that only molecular dynamics simulations may be utilized for the calculations of mobility gaps; the strong vibronic coupling at mobility edges would make the use of the frozen-phonon method inaccurate and most likely would lead to numerical instabilities when computing second derivatives using finite differences.

Fig. 4.

Joint probability distribution of electronic energies (Kohn–Sham eigenvalues) and IPRs for an a-C sample with static HOMO-LUMO gap of 0.84 eV. A–C and D–F show the results of classical and quantum first-principles molecular dynamics simulations, respectively; Left (A and D), Center (B and E), and Right (C and F) display results at 100, 300, and 500 K, respectively. Black stars represent the IPR and electronic energies for the configuration that is the closest local minimum to the a-C sample and the black dashed line shows the chosen IPR cutoff for the calculation of the mobility gap (main text).

As expected, we find that most of the electronic states are more delocalized in the presence of NQEs than in classical simulations. However, there exists a much smaller proportion of states whose localization is instead increased. The effects of nuclear quantum fluctuations are prominent at 100 K, still notable at 300 K, and becoming less important at 500 K.

Interestingly, our computed mobility edges show marked differences when obtained with classical and quantum simulations. The electron–phonon renormalizations of the mobility gaps are given in Fig. 5 and Table 1. While classical simulations show a variability of mobility gaps between 4.27 and 4.92 eV at T = 0, among the different samples, quantum simulations yield smaller and more uniform values between 4 and 4.2 eV. Additionally, the behavior of mobility gaps as a function of temperature is qualitatively and quantitatively different in classical and quantum simulations. In amorphous semiconductors, electronic transport occurs via both band and variable-hopping conduction. The component of the conductivity arising from band conduction is determined by an activation energy that is given by the mobility edge: , where denotes the mobility gap at temperature T. Hence, our results indicate that the influence of NQEs on band conduction is expected to be substantial at temperatures way above room temperature (Fig. 5). In addition, NQEs will likely influence variable-hopping transport mechanisms as well. Although we cannot make quantitative predictions, we can draw interesting conclusions from our results on EDOS and orbital localization (Figs. 2 and 4): The valence mobility edge moves closer to the defect-like states in quantum simulations and the distribution of defect-like states broadens, pointing at a more favorable hopping mechanism between tail states at the quantum level of theory, due to an increased overlap between states. Overall, our results show that NQEs enhance both conduction and hopping transport mechanisms in DLCs. We note that, irrespective of whether a quantum or classical description is adopted, there is no correlation between mobility and HOMO-LUMO gaps and hence the latter should not be used to make predictions, even qualitative ones, about transport properties.

Fig. 5.

(A) The static mobility gaps, E (black squares) and quantum mobility gaps at 0 K (red asterisks) for different amorphous C samples shown as a function of the respective static HOMO-LUMO gaps (E). The quantum mobility gap at 0 K is calculated as the sum of E and the mobility gap renormalization (ZPR) (Table 1). (B–D) The temperature dependence of the electron–phonon renormalization of the mobility gaps obtained from classical and quantum first-principles molecular dynamics simulations for three amorphous C samples. The static HOMO-LUMO gaps of these samples are 0.84, 2.51, and 3.63 eV, respectively, whereas the static mobility gaps are 4.27, 4.60, and 4.92 eV, respectively. The symbols represent the simulation results and the solid lines represent the simulated results fitted with linear (classical) and Viña (quantum) models. See for mobility gap renormalization results for four additional samples.

Discussion

We investigated the influence of nuclear quantum effects on the structural and electronic properties of diamond-like amorphous carbon, using quantum simulations and analyzing several representative configurations with various defect sites. We found that the influence of the quantum motion of the nuclei is minor on structural properties but substantial on electronic properties, including both fundamental and mobility gaps.

Despite the presence of localized electronic states near the Fermi level, introduced by disorder, the renormalization of the electronic gap due to quantum vibronic coupling is similar to that found in crystalline diamond; namely, it varies from –100 to –400 meV, depending on the specific configuration, and is comparable in magnitude to that of the crystal [–605 meV (30), when computed at the same level of theory]. This finding is consistent with the study by Atta-Fynn et al. (50), where for amorphous silicon, large electron–phonon coupling was observed for localized states although the study did not include NQEs. Overall, the zero-point quantum motion of the nuclei renormalizes the HOMO-LUMO gap by 8 to 12% of the static gap.

We also calculated the static mobility gaps and their electron–phonon renormalizations, amounting to 5 to 16% of the gap value; we obtained ZPR values that lie in a range of –230 to –800 meV. Surprisingly, for some samples, ZPR values are even larger than that obtained for the indirect gap of crystalline diamond. We found that the description of the nuclear motion obtained in classical FPMD simulations is inadequate to properly evaluate the effect of vibronic coupling on fundamental and mobility gaps. We showed that, even at room temperature, neglecting NQEs would greatly underestimate the electrical conductivity of diamond-like carbons, possibly up to one to two orders of magnitude, although we could not make quantitative predictions on doped samples. Our results also indicate that frozen-phonon and many-body perturbative methods based on the harmonic approximation should not be used for amorphous carbon, due to the inaccurate description of electron–phonon and phonon–phonon couplings in the systems. Although for some specific configurations, the FPH method may yield a qualitative correct description.

Work is in progress to address the impact of nuclear quantum effects on the electronic and transport properties of amorphous carbon with densities lower than diamond and with various degrees of hydrogenation. In addition, work is ongoing to improve the level of electronic structure theory, by utilizing GW corrections with a transfer learning-accelerated GW approach, which is necessary due to the large number of configurations required to perform accurate statistical averages.

Materials and Methods

We started our first-principles simulations from configurations obtained in a previous first-principles MD study (21). We used the Qbox code (51) for classical FPMD simulations and utilized a stochastic velocity rescaling thermostat (33). In the case of quantum simulations we adopted a colored noise generalized Langevin equation thermostat also known as a quantum thermostat (QT) (32). We used Qbox coupled (30) to the i-PI code, where the i-PI driver (52) moves the nuclear coordinate and Qbox serves as a DFT engine, to compute forces from DFT. We also utilized the i-PI–Qbox coupling scheme for FPH calculations, where the output of the calculations was postprocessed using the pyEPFD package (see section S5 of ref. 30) for the computation of the renormalized band gaps. The PyEPFD package is currently not available to the public and will be released in near future.

For all calculations, we used the PBE exchange correlation functional (41, 42) and norm-conserving pseudopotentials (53) with 50 Ry kinetic energy cutoff.

We chose an MD time step of 0.5 fs and we performed 20-ps-long simulations at 100 K and 15-ps-long simulations at all higher temperatures. We calculated all equilibrium properties (PCFs, HOMO-LUMO gaps, IPRs, mobility gaps) by sampling the trajectories after discarding the first 2.5 ps of our simulations.

We computed the Wannier function centers at 25-fs intervals for 100- and 500-K simulations using the simultaneous diagonalization algorithm implemented in the Qbox code (54). We postprocessed trajectory files using the TRAVIS code (55, 56) to compute the radial distribution functions. For all simulations, the Kohn–Sham Hamiltonian was diagonalized and eigenvalues were computed every 4 fs and HOMO-LUMO gaps were computed as averages.

Following the protocol described by Herbstein (57), the Debye temperature of DLC is estimated by fitting the previously measured heat capacities (31) with the following equation: , where C, N, k, and Θ denote the measured heat capacities, Avogadro’s number, Boltzman’s constant, and the Debye temperature, respectively.