Tuning phononic and electronic contributions of thermoelectric in defected S-shape graphene nanoribbons

Thermoelectrics as a way to use waste heat, is essential in electronic industries, but its low performance at operational temperatures makes it inappropriate in practical applications. Tailoring graphene can change its properties. In this work, we are interested in studying the transport properties of S-shape graphene structures with the single vacancy (SV) and double vacancy (DV) models. The structures are composed of a chiral part, which is an armchair graphene nanoribbon, and two zigzag graphene ribbons. We investigate the changes in the figure of merit by means of the Seebeck coefficient, electronic conductance, and electronic and phononic conductances with the vacancies in different device sizes. The transport properties of the system are studied by using the non-equilibrium Green’s function method, so that the related Hamiltonians (dynamical matrices) are obtained from the tight-binding (force constant) model. The maximum figure of merit (ZT) obtains for the DVs in all lengths. Physical properties of such a system can be tuned by controlling various parameters such as the location and the type of the defects, and the device size. Our findings show that lengthening the structure can reduce phononic contribution, and single vacancies than double vacancies can better distinguish between electronic thermal conductance behavior and electronic conductance one. Namely, vacancy engineering can significantly increase thermoelectric performance. In the large devices, the SVs can increase the ZT up to 2.5 times.

www.nature.com/scientificreports/ the dimensions decrease, quantum confinement effects become important 17 , which can help maximize ZT 18 by manipulating physical properties. This can be accompanied by phonon scattering due to nanostructure boundaries 19 . In 2D materials, especially hexagonal structures, the edge geometry of a ribbon provides a degree of freedom to tailor its physical properties 16,20,21 . Introducing defects, doping, and applying mechanical strain can also alter the physical properties of graphene [22][23][24][25] .
Since the shape and geometry of the nanodevices are important in tuning physical properties at the nanoscale, we are interested in studying the S-shape graphene structures with three different lengths. S-shape structure is a mix-up of zigzag and armchair edge geometries, which can help to tune physical properties. In an S-shape graphene nanoribbon (GNR), electronic contributions can be significantly altered due to the quantum confinement and edge effects 26 . In this work, the temperature is considered 350 K, close to what is to be controlled in processor units 27,28 . To evaluate the thermoelectric performance, the figure of merit, a dimensionless parameter, is investigated. The figure of merit can be calculated as ZT(µ, T) = gS 2 κ e +κ ph T , with g as electronic conductance, S as Seebeck coefficient, κ e as electronic thermal conductance, κ ph as phononic thermal conductance, and T as absolute temperature. These parameters are individually plotted for each of the studied structures.
The two experimentally observed vacancies, single vacancy (SV) and divacancy (DV) 29 , are introduced, and their impact on both electronic and phononic contributions related to thermoelectric performance is studied.
We have used the non-equilibrium Green's function (NEGF) method to calculate the interested quantities. Hamiltonians are obtained from the tight-binding (TB) approach by considering up to third nearest-neighbor (3NN) interactions. To be more accurate 30 , overlap integrals are taken into account. For phononic thermal conductance, force tensor matrices obtained via the force constant (FC) model by considering up to 4NN.
Vacancies are introduced and named as indicated in Fig. 1, e.g., a single vacancy located at the eleventh atomic position in the armchair direction and the fifteenth one in the zigzag direction is identified by SV-11-15; the number of atomic positions is also presented for each direction.
The article is arranged as follows; in the next section, we will describe the model, with a brief introduction on the TB, and FC formulations. Results and discussion are in section three. In the last section, we conclude our study.

Model and methods
In this section, we describe a system consisting of left and right contacts and a central device connected to them.
To start, a schematic of the structure of system is presented in Fig. 1. The system is divided into three parts with black boxes; the device section is an S-shape graphene structure, the right and left contacts are two semiinfinite ZGNRs with a width of 12 atoms. The grey dashed boxes show the unit cells in contacts. Also, the first, second, and third nearest neighbors are displayed with concentric circles in the left contact, so that magenta dashed circle shows the first nearest-neighbor domain, the cyan dashed circle shows the 2NN domain, and red dashed circle indicates the 3NN domain. Vacancies are identified by their position in the armchair and the Three concentric circles show the first, second, and third nearest-neighbor domains, respectively, by magenta, cyan, and red colors for a selected atom. Vacancies were considered in different places in the hatched area (the green arrow shows this direction). Numbers in the device were used for this purpose as described in the text. Vacancies are identified with the general form of VT-m-n-or, in which the VT is the vacancy type, here it can be SV or DV, m and n are atomic positions in the armchair and zigzag directions, respectively, and the last part indicates respective orientation to the nanoribbon axis. For single vacancies, it is omitted, but for divacancies, the relative orientation of the hypothetical line between two removed atoms determines the last part. If a DV is perpendicular to the ribbon axis (or parallel to the ribbon width), it is indicated with 'pr' , it is indicated with 'or' . To be more precise, the cyan box in Fig. 1 indicates a divacancy DV-7-8-or, in which its first atom (as numbers, from left to right) is in the seventh atomic position in the armchair direction, and the eighth atomic position in the zigzag direction, this divacancy is oriented respect to the ribbon axis which is indicated by "or" in the name of the DV. The direction, in which vacancies move in the structure is marked by the green arrow. The dashed lines are bonds that are affected throughout the study. The hatched area shows the zone where the vacancies are introduced.
To employ the NEGF method for electronic and phononic contributions, matrices that describe electron and phonon energies and their interaction with nth nearest neighbors, are essential. To form matrices for electrons (i.e., Hamiltonians) in the tight-binding approach, the unit cell should be defined (as depicted in Fig. 1 with dashed gray rectangles). In the non-orthogonal tight-binding approach, the Hamiltonian of the system, its elements, and the elements of overlap matrix are as 7 : where ε i is the on-site energy and t i,j , and S i,j are the interatomic and overlap parameters, respectively. There are several sub and superscripts that LC means the left contact, RC means the right contact, and D is the device. As indicated in Fig. 1, H LC 0,0 is the Hamiltonian of the unit cell 0 in the left contact, and H LC −1,0 is the coupling Hamiltonian between the unit cell number − 1 and 0 in the left contact.
Hopping and overlap parameters are presented in Table 1 as reported in Ref. 31 . The electronic energy dispersion for a periodic system, like the left contact, can be obtained by solving the eigenvalue problem 32 : where H k and S k are given by: Overlap parameters (eV) www.nature.com/scientificreports/ in which k and a are the wave vector and the lattice constant, respectively. The transmission probability (for electrons T e , and for phonons it is shown by T ph ) can be calculated using Green's function method 33 ; details are provided in the Supplementary Materials. By having transmission probability, one can calculate the electronic conductance g(µ, T) , the Seebeck coefficient S(µ, T) , and the electronic thermal conductance κ e (µ, T) as 11,34 : here e is the elementary charge, and L n is given by: which its numerical form is as follows: with h as plank constant and k B as the Boltzmann constant. This is the discrete form of integral. The summation is over the whole energy range. By considering l as total steps, the integration element for numerical integration in the rectangular method is The secular equation for phonons, which derives from Newton's second law, is: in which, U is the matrix containing the vibrational amplitude of all atoms, ω is the angular frequency, and D is the dynamical matrix: where M i is the mass of the ith atom, and K i,j represents 3 × 3 force tensor between the ith and jth atoms: with θ ij as the angle between the ith and the jth atom. The unitary matrix U θ i,j is defined by the rotation matrix in a plane as: also K 0 i,j is given by: where ϕ r , ϕ t i , ϕ t o are force constant parameters in the radial, in-plane, and out of plain directions of the jth atom, respectively. To be more clear about these matrices, e.g., for the D D , which represents the dynamical matrix of the device section, regarding Eq. (11), and to write what each atom feels (or when i = j ), one must consider all 4 NN effects in the summation, including atoms in the neighboring unit cells, i.e., For coupling terms like the D D−RC elements, the interaction between the first atom of the device and the first atom of the right neighbor (i.e., diagonal elements) is already accounted, so one can safely set this to zero. Force constants 35 , and other essential parameters are presented in Table 1. www.nature.com/scientificreports/ Phononic band structure can be obtained by solving the following eigenvalue problem 36 : with r i,j = r i − r j as the distance between the ith and jth atoms, and k as the wave vector. To calculate the phononic density of states or vDOS, with v as vibrational, one can use vDOS = − 2ω π Im[Trace(G(ω))], within Green's function method 37,38 , or by using Gaussian smearing of the Dirac Delta 7 : with n as band index of the phonon, q as wave vector, and η is a small positive number. By having the transmission function, one can obtain vibrational conductance as 11 : When L D is much shorter than the phonon mean free path (MFP), phononic transport is considered ballistic 39 . The MFP is also a function of the width of the GNR 24,39 . Here, we assume phononic transport is ballistic, because the devices are small enough, and the width of the system is relatively small.
In the next section, we investigate the length of the structure and defect location effect on the electronic and phononic contributions in the ZT formula for the studied structures.

Results and discussion
We have studied three structures with different lengths to find out what contributions are played a major role in determining thermoelectric performance in S-shape graphene ribbons and seeking any meaningful vacancy place impact on it. First, an S-shape device with L D ≈ 34.43 Å is studied. For this configuration, we have moved vacancies in the shaded area between two green dash-dotted lines in the green arrow direction indicated in Fig. 1. All terms in the ZT formula are plotted by taking the ratio between the defected and pristine values. We call this hereafter the ratio. Also, we will mention κ e and g as electronic terms.
The 30° chiral part of the system is a 10-AGNR with a finite length. Electronic band structures and phononic dispersions plotted in Fig. 2 are for an infinite length of GNRs. A 10-AGNR with infinite length is a semiconductor, as evidenced by the electronic band structure plotted in Fig. 2a, with an energy gap ~ 1.1 eV. The Fermi energy is zero. The band structure is not symmetric respect to the Fermi energy, which is because of the inclusion of the 3NN (with overlap) TB model.
The left part of panels (b) and (d) in Fig. 2 show the phononic dispersion and the vDOS for the 10-AGNR and the 12-ZGNR, respectively. The acoustic bands are located at low-frequencies and they usually possess much higher group velocities compared to those of the optical ones, so they contribute mostly to the thermal transport. Therefore, low-frequency bands play a dominant role in thermal conductance. Comparing low-frequency phonon bands of the AGNR with ZGNR, shows that phononic bands in the AGNR are less dispersive respect to the dω. www.nature.com/scientificreports/ ZGNR, therefore, since the phononic transmission coefficient is equal to the sum of phonon modes, dispersive phonon bands lead to larger values of the transmission coefficients. The above discussion suggests that a ZGNR can act as good thermal conductors, while the AGNR is a better candidate for thermoelectric applications. Figure 2c shows electronic dispersion for the 12-ZGNR, with no energy gap. Therefore, it is a meal, which is consistent with the fact that all zigzag nanoribbons are metal in simple TB model 40,41 . Metallicity can strongly reduce the ZT value in ZGNRs 42 . The vibrational density of states, vDOS, are plotted in the right part of panels (b) and (d) in Fig. 2. Also, for low-frequencies, the AGNR has the lower vDOS values than the ZGNR (check the black dotted line on Fig. 2e). The thermal conductance of the AGNR is lower than that of the ZGNR, which is due to the lower phonon density of states at low-frequencies. The above results are consistent with Refs. 31,43,44 . As shown in Fig. 3, all electronic and phononic terms, together with the Seebeck coefficient and ZT, show a symmetric behavior respect to the middle of the device (purple dashed line in Fig. 1). Figure 3a shows the maximum ZT (magenta squares) as a function of SV locations. The chemical potential for each ZT max is indicated by numbers near the symbols. In this spectrum, the maximum figure of merit, ≃ 0.12 , is for SV-10-12 (or SV-13-19) at µ = − 1.54 eV.  www.nature.com/scientificreports/ According to Fig. 3b, SVs can alter the ZT max (magenta circles) up to 1.18 times, and Seebeck coefficient square (black stars) up to 4 times greater than the pristine structure. As shown in Fig. 3c, by increasing κ e and g and reducing S 2 , vacancies can also degrade the ZT. Electronic and phononic terms can increase or decrease by varying the place of the defect. Inclusion of an SV, decreases the phononic term, and changes it to ∼ 0.8 of its pristine case, and the electronic ones fluctuate between 1.5 and 0.1 times of its defect-free structure. The κ e is shown by the square symbol, g by the circle symbol, and κ ph by stars. The SV-9-11 reduces the ZT to 0.71 times of the pristine structure. This reduction is because of increasing κ e and g to 1.48 times, accompanied by reducing S 2 to 0.48 times compared to the pristine case.
The highest ZT value with a DV, ≃ 0.21 , in L D ≈ 34.43 Å , is achieved for the DV-6-8-pr (DV-16-23-pr) (Fig. 3d) at µ = − 1.52 eV . In the presence of the DVs, S 2 and the ZT max are about 4.3, and 2.1 times greater than the pristine structure (Fig. 3e), respectively. Also, we should mention that the ZT for this vacancy can decrease by 0.52 times. The behavior of electronic and phononic terms, κ e , g , and κ ph are shown in Fig. 3f. A meaningful change in the figure of merit occurs whenever electronic terms take apart from each other. Since double vacancies can reduce available conduction channels for phonons, they affect κ ph more than single vacancies, as evidenced in Fig. 3f. By recalling the ZT formula, higher g , i.e., lower suppression on g in comparison to κ e , can help to achieve the higher ZT and vice versa. In this length, all terms are affected comparably large in the case of the DVs than the SVs.
One can see that the SV-7-8 and SV-16-23 are located in the chiral section of the system. Figure 3c shows that electronic terms between these two are symmetric with respect to the center of the device section, and out of that, their behavior is different.
Because of the symmetric behavior discussed earlier, one can only move vacancies up to the middle of the device.
In this step, we can track the changes of the above quantities for a longer length of the device. The second structure has L D ≈ 47.96 Å, which we did not address it here. The results for this length are presented in the Supplementary Materials. The third structure has L D ≈ 59.03 Å. As Fig. 4a shows, the ZT max occurs for a single vacancy in the middle (SV-18-24) with µ = 0.99 eV . Figure 4b indicates that the presence of an SV in the middle of the device can increase the ZT value by 2.5 times with respect to the pristine case. It also reduces the ZT to 0.6 times of the device with L D ≈ 34.43 Å in the presence of SVs. The ZT also shows higher performance in comparison to the shorter case. Figure 4c shows a slightly downward trend in the phononic term. As the SVs get closer to the middle of the device, the ratio of the phononic term, reduces. The phononic term shows a reduction of ~ 0.8 times of the clean case, which is close to that of the shorter length of the device. In the first system, electronic terms follow each other closely, but this correlation tends to demolish in longer device, as it can be seen for electronic conductance and electronic thermal conductance for the SV-18-24.
The DVs have a slightly poor effect on the thermoelectric performance in comparison to the case of L D ≈ 34.43 Å . The ZT has its highest value with ≃ 0.21 at the chemical potential of − 1.52 eV for the DV-6-8-pr. However, this double vacancy in the shorter length also gives the highest ZT. Phonon thermal conductance reduction is almost negligible compared to the shorter length, but DVs induce stronger fluctuations in electronic terms, as evidenced in Fig. 4f. The ZT has a smaller value for this length in the presence of DVs.
By comparing panels (c) and (f) of Fig. 4, and the same ones in Fig. 3, one can conclude that in the case of SVs, as the length increase, the figure of merit also rises, and κ e and g become more independent. For DVs, lengthening device reduces the ZT, but it decouples electronic terms. In the longer length, the fluctuation of chemical potential corresponds to the ZT max , becomes smaller.
Vibrational local DOS (vLDOS) for the first and third lengths are shown in Fig. 5. Regarding Fig. 3c,f, the lowest thermal conductance is for the SV-10-12 (or SV- [13][14][15][16][17][18][19], which is placed on an SV with the highest vLDOS; the same is true for the DV-8-11-pr (or DV-14-20-pr) (Fig. 5a). Increasing the device length causes the system to experience the lower vLDOS (Fig. 5b). In the long devices, the higher vLDOSs are related to the edge atoms, suggesting edge defects can induce a more substantial effect on κ ph .
We also studied electronic and phononic transport properties for ZGNR with the exact width of the contacts, and an AGNR similar in width to the chiral section of the device, together with the two S-shape GNR lengths studied here. As it is depicted in Fig. 6a, electronic conductance for large chemical potentials (> 1.5 eV) is almost close for 10-AGNR (blue line) and 12-ZGNR (green line), and two S-shape graphene structure lengths (shorter length is indicated by red and the longer length is indicated by cyan), for the structure with a longer length, the behavior of the electronic conductance for values of µ close to zero is almost similar to that of the 10-AGNR.
The behavior for electronic thermal conductance is in a similar trend for g , as depicted in Fig. 6b. Seebeck coefficient is plotted in Fig. 6c versus chemical potential, which for the shorter length, it has smaller values than the longer length. Also, the peaks between − 1 and 1 for chemical potentials are opposite, which shows the change of charge carriers 11 , and closeness to the AGNR behavior. The κ ph decreases by lengthening the system (check the results for L D ≈ 47.96 Å in the Supplementary Materials, which confirms this trend), which is the attribution of anharmonicity of phonon modes in chiral and zigzag parts (Fig. 6d).
As the length of the chiral part increase, the AGNR characteristics become stronger. The Seebeck coefficient enhances with a bandgap 45 , so one should be aware of these tradeoffs between phononic and electronic terms and the length impact when designing an efficient thermoelectric structure. Phonon mismatch between the AGNR and the ZGNR parts becomes stronger as the AGNR characteristics become dominant by lengthening the chiral section. Moreover, vibrational modes occupy a narrower frequency range in comparison to the ZGNR sections (Fig. 2b,d), which limits thermal conductance 43 .
We choose the system with L D ≈ 59.03 Å due to the strong increase of its ZT in the presence of SV-18-24. Here, we study the transmission coefficient for this system, both for electrons and phonons. According to Fig. 7a, there is a small peak close to zero energy (cyan line) in the transmission spectrum of the system caused by the SV-18-24, which shows the close behavior of transmission coefficient of the defected S-shape GNR with the www.nature.com/scientificreports/ 12-ZGNR (green line). Although, the behavior of the pristine S-shape GNR is likely dominated by its chiral section (red line). The transmission spectrum of the pristine system has a semiconducting bandgap, similar to the 10-AGNR gap (blue line). As one can see in Fig. 7b, low-frequency phonons are affected more than high-frequency phonons in the pristine S-shape GNR (red line), compared to the 10-AGNR (blue line), and 12-ZGNR (green line). Besides, the transmission coefficient of the bent systems shows more suppression in the range of 0 to 1.5 × 10 14 Rad/s. This is, in general, a good change since low-frequency phonons are known to be more responsible in thermal conductance 46 .
Additionally, the impact of temperature for various chemical potentials on the ZT is plotted in Fig. 8 for the SV-18-24. As evidenced in this figure, the maximum ZT of 0.32 can be achieved at 1000 K for µ ≈ 0.98eV.  www.nature.com/scientificreports/

Conclusions
In summary, we have studied the behavior of electronic and phononic contributions in S-shape graphene structures with different lengths by including single and double vacancies in different locations. We have investigated the electronic and phononic transport properties of the system by using the NEGF method and the tight-binding approach by considering the 3NN for electronic contribution, and the force constant model with 4NN, for the phononic part. Our numerical results show the symmetric behavior of the terms involved in the ZT, in the presence of defects respect to half of the device. The maximum of ZT is obtained for the DVs in the shorter length of the device. Also, the device with a longer length has a ZT of < 0.21. By increasing the length of the chiral part of the system, the ZT can be enhanced, not only by reducing the vibrational contribution, but also by separating the electronic terms.
Single vacancies can detach the electronic thermal conductance and electronic conductance better than double vacancies; namely, the SVs can magnify the ZT up to 2.5 times in the studied structures. Detaching κ e and g from each other is a way to improve the ZT performance.
Lengthening the system causes the chiral section (AGNR) characteristics to become dominant. One can tune the electronic thermal conductance and electronic conductance properties of the system by altering the parameters, such as the type and the location of the vacancy defects and the device size.