Charge carrier nonadiabatic dynamics in non-metal doped graphitic carbon nitride

Two-dimensional (2D) graphitic carbon nitride (GCN) based on heptazine units has become one of the most promising metal-free polymeric semiconductors in photocatalysis research.^1–5 The narrow bandgap of around 2.7 eV, low cost, chemical and thermal stability, and unique surface and electronic properties make GCN an ideal material in photocatalytic applications.^6–8 At present, GCN is used mainly in hydrogen production via water splitting, photodegradation of pollutants, and several other photocatalytic oxidation and reduction reactions.^1,2,4,8–12

Although GCN has a wide range of potential applications, poor specific surface area, low visible light absorbance, rapid electron–hole recombination, and low charge separation efficiency limit the photocatalytic performance of the material.^3,13,14 Therefore, the modification of GCN is one of the hot research topics in the field of catalysis. In the past, several methods have been developed to overcome the inherent limitations of GCN and improve the efficiency of the material, such as chemical doping modification, heterostructured nanocomposite fabrication, and noble metal deposition.^15–18 Element doping has been considered as one of the most easy and effective means to tune the electronic structure of the material. A few experimental and theoretical works reported on defect modification in g-C₃N₄ have shown to improve the light absorption, accelerate the separation and transport of charges, and prolong the carrier lifetime.^15,19–21 In consideration of the metal-free nature of the modified GCN, non-metal doping has attracted more attention, as also explored in this work.

Non-metal elements, such as O, N, C, P, and B, have been used to improve the photocatalytic performance of GCN through the modification of electronic structures and morphology of the material.^15,22–30 For example, work by Qiu et al. demonstrates higher photocatalytic efficiency in O-doped GCN due to increased separation of electron–hole pairs and improved visible light absorption.²⁴ Thaweesak et al. reported enhanced photocatalytic water splitting in B-doped GCN due to decreased bandgap and decreased photoluminescence (PL) intensity.²⁸ However, a detailed theoretical charge carrier dynamics study on the doped GCN systems is still lacking. Quantum dynamics simulations would provide a detailed and atomistic level understanding of the effect of the doped elements on the electronic properties and charge dynamics of GCN.

Motivated by both experimental and theoretical works,^15,22,24,28 we carry out ab initio nonadiabatic (NA) molecular dynamics (MD) simulations combined with time-dependent density functional theory (TDDFT)³¹ to investigate the photoinduced charge trapping and recombination dynamics in the doped GCN. Here, we consider three doped systems: O doped at N (O–N–GCN), B doped at N (B–N–GCN), and B doped at C (B–C–GCN). The simulations show that the presence of trap states in doped GCN facilitates rapid charge recombination by providing additional nonradiative relaxation pathways as compared to the pristine system. Moreover, the decreased energy gap between the trap and the band edge is responsible for better visible light absorption. This work focuses on the excited state dynamics, leading to charge separation and recombination that are part of the photocatalytic processes. A full evaluation of the photocatalytic capabilities requires a comprehensive study of the micro-reaction processes involved.

The electronic structure and adiabatic molecular dynamics (MD) trajectories are obtained using the Vienna Ab initio Simulation Package (VASP)^32–34 and the projector augmented wave method. The exchange–correlation energy is calculated with the Perdew–Burke–Ernzerhof (PBE) functional³⁵ and ultrasoft pseudopotentials,³⁶ unless otherwise specified. Since PBE tends to underestimate bandgaps, we employ the HSE06 functional^37,38 for initial electronic structure calculations. The van der Waals interactions are described via optB86b-vdW functional employed in the vdW-DF method.³⁹ A 900 eV kinetic energy cutoff³⁹ and a 3 × 3 × 1 Γ-centered k-point mesh are used. The high energy cutoff allows us to obtain accurate NA couplings (NAC) and provides an accurate evaluation of the vdW-DF corrections. To obtain an accurate electronic structure, a denser 9 × 9 × 1 k-point mesh is employed with the PBE functional. All systems are simulated using $3$ × $3$ unit cell with lattice parameters as given in Ref. 40 with 42 atoms. A 20 Å vacuum layer is introduced in the z-direction to eliminate interactions between the layers. The optimized structures of systems under investigation are shown in Fig. 1. After optimization, the structures are heated to 300 K using uniform velocity rescaling. Then, 5 ps adiabatic MD trajectories are obtained with 1 fs time step in the microcanonical ensemble (NVE) and further used for NAC calculations using the CA-NAC package.⁴¹ A 3 × 3 × 1 k-points grid is used to obtain the wavefunctions required for the NAC calculations. The NACs are computed at the Γ-point only since the defect states are localized and exhibit weak k-point dependence. Using a single k-point reduces the extent of interaction between defect periodic images, since defect states corresponding to multiple k-points are not considered.

Hundred geometries from the initial 2 ps MD trajectories are used as the initial condition for the nonadiabatic molecular dynamics (NAMD) and 100 stochastic hop sequences are sampled for each initial condition. To make the calculations computationally affordable, the NA Hamiltonian is iterated multiple times to perform the long 1000 ps ab initio NAMD simulations. The main assumption involved in iterating the Hamiltonian is that the short simulations sample the motions that are relevant for the NA dynamics. The evolution of the ab initio energy levels (shown in Fig. S1) demonstrates that the energy levels fluctuate rapidly and that multiple periods of the large amplitude fluctuations are sampled by the 5 ps trajectories. 5 ps correspond to the period of an oscillation with a frequency of about 6 cm⁻¹. The system is composed of light elements, with bond stretching and bending frequencies of around 1000 cm⁻¹ and above. The only motions that are not sampled by the 5 ps trajectories are very slow acoustic modes. Sampling of such modes would require both longer trajectories and larger simulation cells. However, the NAC is proportional to the nuclear velocity, and higher frequency modes have larger velocities at a particular temperature. Therefore, higher frequency modes should contribute most to the NAC. The assumption is tested by splitting the NAMD Hamiltonian into shorter parts and performing the calculations based on each part separately, as discussed below.

The volume per charge carrier is smaller in the simulation cell than in most experiments due to computational limitations on the system size. However, there is only one electron and one hole in the simulation. Higher-order states, such as trions or bi-excitons, are not considered, and Auger-type charge–charge scattering that would occur in experiments at high carrier concentrations is not included in the simulation. Even though the charges are confined to small parts of the polymer, only processes occurring in the low carrier limit are considered, representing effectively the situation encountered in most experiments. By not allowing carriers to extend beyond the simulation cell, we put them close to each other. Therefore, we describe situations in which a carrier has already diffused to the defect and the free carrier with which the trapped carrier recombines is nearby. Long-range carrier transport requires a different type of simulation involving the modeling of carrier diffusion.

The simulations are performed using the mixed quantum–classical approach implementing NAMD within TDDFT in the Kohn–Sham (KS)^42,43 framework. A higher level of theory could provide a more rigorous description. In particular, it is desirable to use hybrid DFT functionals that produce better bandgaps and to include excitonic effects through either linear-response TDDFT wavefunctions or GW and Bethe–Salpeter theories. Unfortunately, both hybrid functionals, and linear-response TDDFT or Bethe–Salpeter theories greatly increase the computational cost. Pure DFT functionals, such as PBE used here, tend to over-delocalize electrons, overestimating overlap of electronic charge densities. Both hybrid functionals and excitonic effects should localize the electronic states. Over-delocalization of the electronic states should produce larger NACs. Underestimation of the energy gaps by PBE should also produce larger NACs. Therefore, the current level of theory likely produces dynamics that are too fast. A detailed description of the methodology adopted in this study can be found in other Ref. 44 and 45.

The NAMD simulations are carried out using the decoherence-induced surface hopping (DISH) method⁴⁶ implemented in the PYXAID software,^45,47 under the classical path approximation (CPA).⁴⁸ The CPA assumes that thermal atomic fluctuations are more significant than geometry changes induced by photoexcitation. In order to test the validity of the CPA, we used the delta-SCF method to compute changes in the excited and ground-state geometries and compared these changes to thermal fluctuations. In particular, we promoted an electron from the KS orbitals corresponding to the valence band maximum (VBM) to the KS orbitals corresponding to the conduction band minimum (CBM). We also computed the geometry of the electron trap state. i.e., we obtained the delta-SCF geometry with an electron promoted from the VBM to the localized electron trap state. The results reported in Table S1 of the supplementary material show that the thermal atomic fluctuations are larger in magnitude than the geometry changes between the ground and excited states, justifying the CPA. DISH incorporates the loss of coherence within the electronic subsystem induced by coupling to quantum phonons.⁴⁹ Decoherence effects become important while studying the slow charge recombination process taking place across a large energy gap.^50–52 The decoherence time is estimated by computing the pure-dephasing time using the optical response theory formalism.^53–55 The method has been extensively applied to study the photoinduced excited-state dynamics in a broad range of periodic materials.^44,56–60

The projected density of states (PDOS) of all three defect GCN systems using the HSE functional are presented in Fig. 2. The corresponding PBE PDOS are shown in Fig. S2. Here, we have shown the spin-polarized PDOS for two of the defect systems, which have an odd number of electrons. As evident from the figure, the defect systems have an additional state close to the band edges, which are labeled as trap states, which is absent in the pristine system. For further calculations, we consider the spin component that has additional states compared to the pristine since this would be involved in the dynamics, limiting the nonradiative recombination in the system. For the B-doped-at N system, there is a midgap state, which serves as an electron trap since the system is electron deficient. For spin-up PDOS of O-doped-at N [shown in Fig. 2(c)], although the trap state is very close to the CBM, it would act as a hole trap state due to the half-filled occupancy of that state, resulting from the extra electron in this system due to the doped electron-rich oxygen. For spin-down PDOS of B-doped-at C, the trap state is close to VBM; however, it acts as an electron trap state since it is empty. The intrinsic bandgap (CBM–VBM) of the defect systems is similar to the pristine system, as expected. However, there are smaller energy gaps in the systems with defect states due to the presence of additional trap states. The presence of an additional smaller gap would be responsible for better absorption of visible light in the higher wavelength ranges, as mentioned in the literature.²³ For an even better understanding of the electronic structure of the systems studied, we have also plotted the electronic band structures of all three systems in Fig. S3 of the supplementary material. The data show that the minimum VBM–CBM bandgap lies at the Γ k-point, which is used for the NAMD simulations.

The corresponding charge densities of the band edge states (VBM and CBM) and the defect states of all the systems under investigation using the HSE functional are shown in Fig. 3, plotted using VESTA.⁶¹ As evident from the figure, the defect states are more localized near the tri-s-triazine rings containing the doped atoms as compared to the delocalized band edge states. Therefore, the VBM and CBM wavefunctions have a larger overlap as compared to that of the overlap of the defect and band edge wavefunctions, leading to stronger NAC for the band edge states. In all systems, the VBM is dominated by N atoms, but the CBM is spread more evenly between the C and N atoms, and the doped atoms have negligible contribution in the PDOS. However, NAC also depends on the energy gap, and since the defect states are closer to the band edge states than the band edges themselves, NACs between the band edges is expected to be lesser than other pair of states.

The nonradiative recombination lifetime depends on the energy gap, NAC, and pure-dephasing time (listed in Table I for the systems under study). The intrinsic bandgap of the doped systems is quite similar to that of the pristine; however, the energy gap between the trap and band edge states, scaled NAC, and pure-dephasing time could lead to different nonradiative recombination lifetimes.

Relaxation of electrons and holes inside the conduction and valence bands involves electronic states from multiple k-points. The relaxation through dense manifolds of states is fast, sub-picosecond, compared to the charge recombination that takes place across a large bandgap. Our simulations, therefore, consider that the charges have already relaxed to the corresponding band edges, eliminating the need to consider NAC between multiple k-points. Thereafter, NAMD is performed considering all relevant electronic configurations constructed from the band edge and trap state orbitals (shown pictorially in Fig. S4). The DISH simulation results (shown in Fig. 4) are fitted to exponential functions, with the fitting constants reported in Table II. In particular, the decay of a state population is fitted to $Pt=Aexp−t/τ$⁠. The timescale of a population rise is fitted to $Pt=B[1−exp−t/τ]$⁠. A = B = 1 for the initial and final states, while these constants are fitted to be between 0 and 1 for the trap states. In the cases when the population dynamics is not complete, in particular, in Fig. 4(c), we have used the short-term linear approximations to the exponential decay, e.g., $Pt=exp−t/τ≈1−t/τ$⁠.

The pristine system has slower nonradiative recombination dynamics compared to the doped systems since the charges recombine across a larger energy gap between the band edges. The presence of trap states leads to shorter nonradiative recombination lifetimes in the doped systems as compared to the pristine system. This is because the presence of trap states provides additional nonradiative pathways for charge recombination. Once the charge is trapped after photoexcitation, the charges rapidly recombine due to the smaller energy gap between the trap state and the band edge state (VBM/CBM) as compared to charge recombination across the band edges in pristine. This result is also in agreement with the shorter lifetime of charge carriers in doped GCN as compared to the pristine.^27,30 To validate the assumption involved in iterating the NA Hamiltonian multiple times to generate long NAMD trajectories, we performed additional calculations by dividing the full ab initio trajectories into shorter, 800 fs parts, and performing the NAMD calculations by iterating the NAMD Hamiltonian of each of the shorter parts (shown in Fig. S5 of the supplementary material). The timescales and errors reported in Table S2 demonstrate that the results are robust and are comparable to the original results reported here.

For systems containing hole traps (O-doped) and electron traps (B-doped), charge trapping timescale is considered as the time involved in trapping of the hole from the VBM to the hole trap state (VBM–defect) and trapping of the photoexcited electron from the CBM to the electron trap state (CBM–defect), respectively, prior to recombination. This charge separation is absent in the pristine system. In O–N–GCN and B–C–GCN, we do not observe a substantial charge trapping time since the trap states are closer in energy to the ground state, and the NAC for these pair of states is large. As a result, once the trap state is populated, it rapidly decays to the ground state. In comparison, for B–N–GCN, the trap state is in the middle of the bandgap, and the NAC from this midgap state to the ground state (VBM–defect) is comparatively small. Consequently, the population stays in the trap state for some time, allowing the population amplitude to grow. In order to test the influence of decoherence on the trap state populations, we have repeated the calculations using the fewest-switches surface hopping (FSSH) method (shown in Fig. S6). FSSH shows qualitatively similar results as obtained with DISH, although the FSSH timescales are much faster, as expected in the absence of decoherence effects. Overall, the results would favor the use of B–C–GCN in photocatalytic applications involving suitable oxidation and reduction reactions over the other two doped systems due to the charge separation and slower charge recombination.

To establish that the excited state lifetime of charge carriers is limited by nonradiative rather than radiative decay, radiative charge recombination lifetimes are calculated by taking inverse of the Einstein coefficient for spontaneous emission, A₂₁, between states 1 and 2 given by the following formula:

Here, v is the transition frequency, and e, ɛ₀, m_e, and c are the fundamental constants. The state degeneracies g₁ and g₂ are considered to be 1 for this system. The dipole strength, D₁₂, is calculated from the transition dipole moment between two KS states (here, CBM and VBM) that are computed using the plane-wave coefficients of the pseudo-wavefunction.⁶² The oscillator strength, f₁₂, is related to the dipole strength, D₁₂, as

The calculated radiative recombination lifetime for the single optimized geometry of all the systems is given in Table III. The emission lifetime is suppressed in O–N–GCN and B–C–GCN, with most suppressed in the latter. Moreover, these lifetimes that occur over a timescale of several nanoseconds are much slower compared to the nonradiative recombination times. Hence, nonradiative electron–hole recombination constitutes the main process responsible for excited state decay. The radiative lifetime is shortest for the B-doped-at N system since the replacement of N with B leads to a larger difference in the negative and positive charges, resulting in its higher transition dipole moment, which is inversely related to the radiative lifetime. Moreover, we have also calculated the radiative lifetimes using five different geometries sampled from the MD trajectories at the ambient temperature (see Table S3). In comparison, the radiative decay times obtained for the sampled MD geometries are similar to those computed for the optimized geometries.

Increased charge separation time or the charge trapping time and longer nonradiative and radiative recombination lifetimes are desirable for a better photocatalytic performance of a material since then, charges would be available for a longer period to participate in the photocatalytic oxidation–reduction reactions. Looking at the timescales for these charge carrier dynamics processes, B–C–GCN is the most suitable candidate, among three doped systems studied here, to find better applications in photocatalysis.

In summary, we have performed a time-domain ab initio study of photoexcitation charge carrier dynamics in three doped GCN by combining NAMD with real-time TDDFT. The study reveals how the electron-rich and electron-deficient elements doped at different atomic sites influence the electronic structure, charge trapping, and nonradiative and radiative recombination in GCN. The doped elements introduce a midgap state or a state close to the band edges that serves as either a hole or an electron trap center. The presence of these trap states and charge localization leads to charge separation in doped systems, which gives an edge in useful applications as compared to the pristine. These states also improve the visible light absorption range of GCN due to smaller energy gaps between consecutive states. Furthermore, the electron–hole recombination lifetime of the doped systems becomes notably shorter, as also observed in experiments due to additional nonradiative relaxation pathways. Slower radiative recombination lifetimes, compared to the nonradiative, validate that the nonradiative electron–hole recombination is the limiting factor for excited state decay. The narrow bandgap, efficient charge separation, and longer charge recombination indicate B–C–GCN as the more efficient photocatalyst, among three systems studied here. Thus, the results show that boron doping at carbon can promote the photocatalytic activity of GCN. This study provides an atomistic perspective toward understanding the charge carrier dynamics in doped GCN and further guides the design of highly efficient doped GCN-based photocatalysts for improved practical applications.

See the supplementary material for the PDOS and the band structure of three doped GCN structures obtained using the PBE functional, ab initio energy levels, active space configurations, justification of CPA, and additional NAMD results using the DISH and FSSH methods.

S.A. and O.V.P. acknowledge the support from the U.S. Department of Energy (Grant No. DE-SC0014429). D.J.T. acknowledges the support from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (Grant No. TG-CHE190008). A.S.V. acknowledges support from the Basic Research Program of the HSE University.

The authors have no conflicts of interest to disclose.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.