Graphene Nanoribbon Based Thermoelectrics: Controllable Self‐ Doping and Long‐Range Disorder

Control of both the regularity of a material ensemble and nanoscale architecture provides unique opportunities to develop novel thermoelectric applications based on 2D materials. As an example, the authors explore the electronic and thermal properties of functionalized graphene nanoribbons (GNRs) in the single‐sheet and helical architectures using multiscale simulations. The results suggest that appropriate functionalization enables precise tuning of the doping density in a planar donor/acceptor GNR ensemble without the need to introduce an explicit dopant, which is critical to the optimization of power factor. In addition, the self‐interaction between turns of a GNR may induce long‐range disorder along the helical axis, which suppresses the thermal contribution from phonons with long wavelengths, leading to anomalous length independent phonon thermal transport in the quasi‐1D system.


Electrical transport property
With the electronic structures predicted by DFT simulations, the transport coefficients were obtained by semi-classical Boltzmann theory with energy dependent relaxation time 5,6 using the BoltzTraP package (modified from v1.2.5). 7 For the functionalized graphene sheets, the hole/electron transport coefficients were calculated using their donor/acceptor GNR components, because the BoltzTraP package does not account for band crossing, and the carrier concentrations of the entire sheets do not reflect the local charge distributions.
The electrical conductivity (σ), the electric part of the thermal conductivity (κ e ), and the Seebeck coefficient (S) were computed by 5,7 σ = e 2 n dk where e is the electronic charge, T is the temperature, µ is the chemical potential, v nk = (1/h) k nk is the group velocity in the n th band at k, nk is the energy eigenvalue, f ( nk ) is the Fermi-Dirac function, τ e ( nk ) and τ p ( nk ) are the electron and phonon relaxation time respectively. All simulations were performed at 300 K, including bands in the energy range of 0.4 Ry around the fermi level.
Since acoustic phonon scattering is the dominant scattering mechanism in GNRs, 8 the relaxation time for electron could be estimated by deformation potential theory, 5,6 where k B is the Boltzmann constant, D ac is the deformation potential coupling constant, ρ m is the mass density, v ph is the sound velocity, A is the area of the sample, and k δ( k − k )/A is the density of states. Since the relevant orbitals are localized in the conjugating area, the effective area A was modified accordingly. v ph = 2.1 × 10 6 cm/s and D ac = 16 eV 9-11 of pristine graphene were applied to all samples.
For very clean electronic systems, the Wiedemann-Franz ratio could be reduced substantially because the e-e scattering event contributes to the thermal relaxation rate but not the charge relaxation rate. Thus, the phonon relaxation time is approximated as 12 τ p = τ e 1 + τ e /τ ee th (5) where τ ee th is the relaxation time of the thermal current due to e-e collisions. According to the previous work that estimates the τ ee th based on quasiparticle lifetime, 12 the value doesn't change too much in the carrier concentration range that is relevant to our system (10 12 -10 13 cm −2 ), therefore we set τ /τ ee th = 0.1 ps in all simulations.

Phonon transport property
The phonon contribution to thermal conductivity (κ p ) of the functionalized graphene sheet and the helix GNRs were calculated by classical molecular dynamics (MD) simulations, as implemented in the LAMMPS package (v2014). 13 All simulations were carried out at 300 K with a time step of 0.05 fs.
For the functionalized graphene sheet, equilibrium molecular dynamics (EMD) based on fluctuation-dissipation theorem was chosen to reduce finite-size effects. κ p of the donor and acceptor domains were estimated by those of the corresponding graphene sheets composed of single type of domains.
According to the Einstein relation, 14,15 κ p was computed from the fluctuations of the heat current, where T is the temperature, V is the system volume set to the product of the sheet area times the interplanar distance (3.35Å) of graphene. [R α (t) − R α (0)] 2 is the mean square displacement of the integrated microscopic heat flux along the α direction with where i and r iα are the energy and position of the i th atom. The microscopic heat flux were recorded for 2×10 7 MD steps (1 ns) after a 100 ps equilibration period to reach convergence.
10 separate simulations with different initial conditions were averaged for each system. Since the difference between the κ p calculated using a 2 × 6 supercell (about 50 × 25Å) and those with larger supercells are less than 5% in our test case, the 2 × 6 supercell is adopted in all other samples.
Even though the second-generation reactive empirical bond order (REBO) potential 16 have been proved to be sufficient for graphene systems, 5 the more expensive adaptive intermolecular reactive empirical bond order (AIREBO) potential 17 is necessary for our cases, which appropriately describes the interactions between conjugating ligands. In order to account for the long-wavelength phonon contributions, the calculated κ p were rescaled using graphene as a reference to the previous theoretical work that accounts for the long-wavelength phonon contribution by the Klemens Model. 5 For the GNR helixes, reverse nonequilibrium molecular dynamics (RNEMD) simulations were performed with the Muller-Plathe algorithm which imposes the heat flux to form temperature gradient. 18 To study the length dependence of κ p , various simulation boxes with periodic boundary condition were used, which lengths along the heat transfer direction were set to twice of the corresponding GNR helixes. Each simulation box is divided into 50 slabs along the heat transfer direction. After a 100 ps equilibration period, momentum exchange between the coldest and the hottest slabs was performed every 10 fs. Our convergence tests confirm that this exchange frequency is sufficiently low such that κ p is independent of the exchange frequency. The nonequilibrium steady states can be achieved within 1 ns simulation time for the longest GNR helix. The heat flux was then computed from the exchanged kinetic energy between two selected layers, where m h/c and v h/c are the atomic masses and the velocities of the hot and cold atoms respectively, t swap is the entire swap time, N swap is the number of swaps. The temperature profile was calculated as where T i is the temperature of the i th slab, N is the number of atoms inside the slab, and p j is the momentum of the j th atom. Finally, κ p was obtained by where A is the cross section of heat transfer. For the combined helix and H-SiNW systems without constraint, A = π(r +s/2) 2 , where r is the radius of the helix and s is the interplanar distance (3.35Å) of graphene. For the combined helix and H-SiNW systems with frozen H-SiNW cores, A = s × w p /sinθ, where w p is the width of the associated planar GNR, and θ = arctan( lc 2πr ) with lc the distance between neighboring turns. The rescaled κ p of a helix is defined as κ p /(sinθ) 2 , which is just a convenient way of excluding the bare geometric factor for direct comparison with the associated planar GNR. The time intervals for calculating κ p were set to be sufficiently long such that both the heat flux and temperature profile are independent of the time interval.
We note that while the NEMD approach simulates the temperature gradient directly and introduces various degree of boundary scattering by modification of simulation box size, it does not account for the contributions from long-wavelength phonons. Nevertheless, NEMD simulation has been successfully applied to illustrate the power law dependence of thermal conductivity on system length in the range of about 10 2 -10 3 nm 20,21 and the transition between ballistic and diffusive transport regioms 20,22 in quasi-1D materials including CNTs 20 and SiNWs, 21 which are consistent with both the predictions from theoretical analysis [23][24][25][26][27] and the violation of Fouriers' Law observed in experiments. 28,29 This indicates that the scattering of long-wavelength phonon transport may be similar to that of short-wavelength phonons, and that the NEMD method may capture the essence of the length dependent scattering phenomena, which enables us to qualitatively investigate the effects of self-interactions in GNR helixes.
Finally, the phonon spectra were computed by the Fix-Phonon package, 30 with the dynamical matrix constructed by observing the displacements of atoms during MD simulations.
Convergence is reached with a simulation period of 400 ps and a length of about 330Å.

Figure of merit
The power factor was calculated by where σ is the electrical conductivity, and S is the Seebeck coefficient. The figure of merit ZT of each material was obtained by where κ p and κ e are the phonon and electron contribution to the thermal conductivity respectively. Given the parameters from both donor and acceptor materials, the figure of merit of the entire system was evaluated as B. Comparison of electronic structures between an entire graphene sheet and its GNR components As illustrated by the prototype system gnraH-w8/gnraCN-w8 (notations as explained in

C. Doping densities of functionalized graphene sheets with various passivations and domain widths
In order to illustrate the possibility to control the doping density of functionalized graphene sheets, we considered a series of prototype systems with various ligands and domain widths, and the electronic structures of the ones exhibiting ground state charge transfer characteristics are summarized in Figure S4 and Figure S5.
The impact of domain width was studied using the gnraH/gnraCN systems. When both the donor and acceptor GNRs are in category-I (gnraX-w6/w12/w18) with comparatively large bandgaps, the energy level alignments tend to be type-I or type-II, and thus the doping densities are zero. When both the donor and acceptor GNRs are in category-II (gnraX-w8/w14/w20) with comparatively small bandgaps, the differences in chemical end groups of the ligands introduce sufficient potential drops inside the graphene sheets to form a type-III energy level alignment and therefore incur ground-state charge transfer ( Figure S4, 1 st and 2 nd rows). In these cases, the bandstructures of all samples are similar to each other, and  Figure S4, 3 rd rows). To gain extra tunability of the functionalized graphene sheets, the donor GNRs gnraH were replaced by gnraOH with electron donating chemical groups. As expected, a higher degree of charge transfer is reached, leading to relatively high doping densities ( Figure S5).

D. Electrical transport properties of GNRs with various widths and passivations
Because of the various quantum confinement, graphene and GNRs exhibit different electronic structures and electron relaxation time as shown in Figure S6, which result in distinct electrical transport properties ( Figure S7). Note that our calculated density of states and relaxation time for graphene is consistent with previous simulations. 5 The electrical transport properties of a GNR depend on its category, width, and passivation, and are especially sensitive to the doping density (carrier concentration). The Seebeck  and power factor (S 2 σ) of our prototype GNR systems in the doping density ranges that could be achieved in our proposed single sheet architecture ( Figure 2) are summarized in Figure S7 and Figure S8, and the following discussions focus on these relevant carrier concentration ranges.
For all of our prototype GNRs, the seebeck coefficient peaks near zero carrier concentration owing to the large derivative of DOS, and then decays sharply with increasing doping density. The S of GNRs in category-I is larger than that in category-II, while both decrease with increasing width (Figure S7 (a)(b)). In contrast, electrical conductivity increases with increasing width, with its optimal values achieved at relatively high doping densities. Notice  that σ has a flatter profile near the optimal value in category-II GNRs than that in category-I GNRs, which may favor the applications that require high stability of σ ( Figure S7 (c)(d)).
The electron contribution to thermal conductivity has a comparably weak dependence on GNR width, and increases with increasing doping density ( Figure S7 (e)(f)). While κ e is generally considered to be much smaller than the phonon contribution in realistic system due to defect scattering, our results suggest that for extremely clean GNRs with acoustic phonon scattering serving as the dominant scattering mechanism, κ e could be as large as E. Comparison of phonon thermal transport between combined GNR-SiNW system and its components In order to confirm that the thermal conductivities calculated by our MD simulations reflect the intrinsic properties of GNRs rather than the influences of GNR-SiNW interfaces, we compared the κ p of combined GNR-SiNW systems and their components in both helix and planar architectures. As shown in Figure S9, with the helix architecture, the κ p of both isolated SiNW and combined helix-SiNW system without constraint increase with increasing GNR length following the power law relation, while κ p of the GNR helix wrapping around the frozen SiNW is almost independent to the GNR length. (here the cross sections of the SiNW and the GNR helix are set to that of the combine system for direct comparison.) Therefore, while the contribution from the GNR helix is more important in the short-length limit, κ p of the combined system is mainly attributed to the transport in SiNW at relatively long lengths.
By contrast, in the planar architecture where the planar GNRs lying on top of Si slabs, κ p of the combined systems are similar to those of the associated isolated GNRs, which increases with increasing GNR length following the power law relation ( Figure S10, the cross sections of the combine systems are set to those of the corresponding planar GNRs for direct comparison). That means the scattering induced by GNR-Si interfaces is weak, and should not be the main reason for the abnormal thermal transport properties observed in GNR helixes.
To further explore the influence of architecture on heat transport, we compared the phonon density of states (DOS) of the isolated planar GNR to that of the combined GNR helix-SiNW system ( Figure S11). While the peaks of the GNR helix-SiNW system is generally broader than that of the planar GNR due to its larger geometric inhomogeneity, the peak positions are quite similar in both cases, implying that the scattering of the phonon modes rather than the mode distribution is responsible for the distinct thermal transport behaviors in different architectures.
Section F. Length dependent thermal conductivity of

GNR helix
While the additional phonon scattering mechanism caused by the long range disorder affects the thermal conductivity only slightly for short GNRs, the reduction of κ p in helix compared to planar architecture becomes non-negligible for GNRs longer than 100 nm as illustrated in Figure S12. Here the κ p of our prototype GNR is estimated by extrapolation from the fitting to power law relation using the MD simulation data for short GNRs as shown in Figure 4(b).
We assume the power law relation is still valid for GNRs with lengths up to 5000 nm at 300 K, based on a previous analytical study on CNT systems. 25 Section G. Phonon DOS of helix-SiNW systems with various side chains For obtaining strong self-interaction effects by tuning the GNR edge passivations, the ligands need to be chosen appropriately such that they are either chemically bonded or sufficiently close to each other as shown in Figure 4(f). Otherwise, the influence of the ligands would be weak because no new phonon transport channel is generated at the ligand-ligand interfaces ( Figure S13).

Section H. Comparison of electronic properties between planar and helical GNRs
As shown in Figure S14, the surface curvature of the helical GNR weakly affects the electronic structure, leading to a slight decrease of both the electrical conductivity and electron contribution to thermal conductivity, while a slight increase of the Seebeck coefficient compared to the associated planar GNR.   Figure S11: Phonon density of states of the isolated SiNW, planar GNR and combined GNR helix-SiNW system. Figure S12: Length dependent reduction of rescaled κ p in GNR helixes compared to that of planar GNRs.