Effective Mass of Quasiparticles in Armchair Graphene Nanoribbons

Armchair graphene nanoribbons (AGNRs) may present intrinsic semiconducting bandgaps, being of potential interest in developing new organic-based optoelectronic devices. The induction of a bandgap in AGNRs results from quantum confinement effects, which reduce charge mobility. In this sense, quasiparticles’ effective mass becomes relevant for the understanding of charge transport in these systems. In the present work, we theoretically investigate the drift of different quasiparticle species in AGNRs employing a 2D generalization of the Su-Schrieffer-Heeger Hamiltonian. Remarkably, our findings reveal that the effective mass strongly depends on the nanoribbon width and its value can reach 60 times the mass of one electron for narrow lattices. Such underlying property for quasiparticles, within the framework of gap tuning engineering in AGNRs, impact the design of their electronic devices.

, and charge ± e associated with a local lattice distortion 28 . Bipolarons, in turn, present two narrower intragap energy levels and stronger lattice distortion than polarons, charge ± 2e and are spinless quasiparticles 28 . The lattice distortions produced by either quasiparticle result in the observed larger effective masses responsible for reduced charge mobility. In this sense, the interplay between the carrier's effective mass and the different properties of AGNRs is a crucial aspect that should be understood to promote the enhancement of graphene-based devices figures of merit.
Herein, we study the drift of charge carriers in AGNRs to phenomenologically characterize their effective masses (m eff ). By means of a 2D generalization of the Su-Schrieffer-Heeger (SSH) model 29,30 , along with a Stokes dissipation model, we numerically investigate the dynamics of polarons and bipolarons in these systems. In the scope of our approach, we determine terminal velocities and effective masses of the charge carriers for AGNRs with different widths. Our findings show that effective mass strongly depends on carrier type and ribbon width, varying up to two orders of magnitude. Importantly, different carbon-based systems (or even inorganic-based nanomaterials 31 ) have different electronic structures and may present distinct responses when it comes to transport properties, in the sense their symmetry and doping level may alter these properties substantially 32 . Zigzag GNR does not present an energy gap with which polarons and bipolarons are usually associated 3,4 . In this way, they are not considered here.

Results
The quasiparticles present local lattice distortions that accumulate charge. Figure 1(a,b) depict the time evolution of charge density in 6-AGNR, where hot colors represent such charge accumulation. In Fig. 1(a), the charge density profile represents a polaron moving under the influence of an external electric field applied in the vertical direction. Similarly, Fig. 1(b) shows a bipolaron, subject to the same electric field strength. Both quasiparticles respond to the applied field. However, the polaron experiences stronger acceleration. The slower response to the electric field as presented for the bipolaron is so in spite of the its charge (2e), which results in twice as much force applied on bipolarons when compared to polarons. Such behavior goes to show that the extra force is not enough to balance the increased inertia of a bipolaron. Indeed, bipolarons carry along a stronger distortion of the nanoribbon lattice. Comparison between both distributions of charge density demonstrates that the polaron's charge is distributed over 40 sites in the vertical direction, whereas the bipolarons' charge spread over less than 30 sites. As such, the polaron is more delocalized than the bipolaron. The combination of more charge confined in shorter  regions results in more significant lattice deformation, which is responsible for the increased inertia observed for bipolarons. Therefore, Fig. 1 illustrates that increased localization will take a toll on charge mobility.
By taking the center of the charge density distribution and registering its position, it is possible to obtain the time evolution of the charge carrier's position in the nanoribbon. Figure 2(a) shows the polaron displacement, in the vertical direction, as a function of time. Each color represents the behavior observed in different AGNR families. Red curves refer to 4-AGNR, green to 5-AGNR, and blue to 6-AGNR. In Fig. 2(a) it can be seen that the displacement curves approach linear behavior, denoting that charge carriers reach a terminal velocity. Similar results are seen in Fig. 2(b) for bipolarons. Note that a curve representing the 3p + 2 family is absent in Fig. 2(b) since bipolarons are not stable in lattices belonging to this family 33 . Regarding polarons, it is worth to mention that this quasiparticles are present only in thinner AGNRs for the 3p + 2 family 33 , and because of that we consider as representative systems the ribbons 5-AGNR and 8-AGNR in which polarons can be formed 33 . Despite similar qualitative behavior, the two quasiparticles show different characteristic times for reaching terminal velocity. For polarons, this time corresponds roughly to 1 ps, around half of the simulation time. In the case of bipolarons, an extra 0.5 ps is necessary. A simple Stokes dissipation model is used to describe the observed carrier motion. This model considers a dissipation term proportional to the first power of the velocity: F d = −bv. Therefore, we can obtain terminal velocity v t when the dissipation term equals the force produced by the electric field v t = qE/b, where q stands for the carrier's charge and E the electric field. Under these conditions, the displacement of the charge center as a function of time (t) is given by where m eff stands for the effective mass of the carrier and b is the drag coefficient. By fitting the displacement curves, changing the values of m eff and b, we can evaluate both polarons' and bipolarons' effective masses as well as the drag coefficients. Note that meff is a fitting parameter. In this sens, we are able to evaluate both effective masses and drag parameter by using its value. Importantly, the main goal of the present work is proposing a route to evaluate the effective masses of quasi-particles in graphene nanoribbons. This quantity is fundamental to experimental approaches to the evaluation of charge mobilities, as well as the understanding of underlying properties of charge transport in graphene systems. Regarding the model used for the determination of the effective mass, as mentioned above, we chose Stoke's model for its simplicity and its phenomenological aspect. A vital feature of this model comes from observation of terminal velocity and thus the necessity of a dissipation term.
Through the procedure described above, it is possible to understand how the ribbon width affects charge carriers inertia. Figure 3(a) shows the polaron effective mass in units of electron mass (m e ) as a function of nanoribbon's width. The blue, red, and green lines correspond to the 3p, 3p + 1 and 3p + 2 families, respectively. In all cases, the effective mass correlates inversely with the ribbon width, ranging from 0.2 to 14 m e . As expected, effective masses are lower for the 3p + 2 family, given its quasi-metallic nature. The remaining two families, however, show no clear ordering, with both curves intersecting each other. Analogous conclusions hold for bipolarons, as shown in Fig. 3(b). www.nature.com/scientificreports www.nature.com/scientificreports/ In contrast to polarons, bipolaron effective masses are considerably larger, ranging from around 10 to near 60 m e in the case of the 3p + 1 family and from 5 to 15 m e in the case of the 3p family. These results reflect the above mentioned more considerable inertia observed in the bipolaron's response to the electric field, responsible for the longer times needed for terminal velocity to be reached. Due to the differences in effective mass, the required electric field intensity required to move both charge carriers differ considerably. The insets of Fig. 3 show representative results of the evaluation of effective mass under different electric fields. One can note that in the case of polarons, the obtained masses are practically field independent. The same situation occurs for bipolarons until fields of around 3000 V/cm. After this critical field strength, the effective masses increase, and the fits to the model are not appropriate, indicating that the model may be not valid in high electric field regimes. This behavior takes place due to strong electron-phonon interaction in GNRs. Below a critical field strength, electrons and phonons are strongly coupled. Above this critical limit, electrons are decoupled from the lattice assuming supersonic velocities. Therefore, the lattice distortions and electron starts to move disconnectedly, and the kinematic model is not valid anymore.
Finally, the interplay between effective masses and ribbon width can be understood from a microscopic perspective by analyzing the charge distribution for quasiparticles in different AGNRsy. Figure 4(a,b) show these distributions for 4-AGNR and 6-AGNR, considering a lattice containing a polaron (left panels) and a bipolaron (right panels). It is clear that as ribbon width increases, the polaron becomes more delocalized, as can be seen by the hotter colors in smaller widths. Despite increasing the ribbon width, the charge tends to concentrate laterally. The same qualitative behavior takes place for lattices containing a bipolaron. The combination of these two underlying effects for the net charge localization leads to an increase in the local lattice deformations associated with the presence of charge for narrower AGNRs. Contrarily, for wider AGNRs, the interplay of these two effects decreases the local distortions that are interacting with the charge. Moreover, in Fig. 4 it is possible to note that lattices containing a bipolaron presents a higher degree of charge localization (that are represented by the signatures in red). Bipolarons quasiparticle have a similar extension to the polaron, approximately 30 Å. Since polarons and bipolarons are composite quasiparticles in which the local lattice distortions are coupled to an additional charge, both evolve in time together during the transport of these quasiparticles. Therefore, the higher the degree of distortion more lattice energy should be transferred between neighboring sites to accomplish the polaron/ bipolaron transport. Consequently, this mechanism for charge transport increases the effective mass of more localized charge carriers.

Methodology
The model Hamiltonian employed here is given by H = H latt + H elec , where the first and second terms govern the lattice and electronic degrees of freedom, respectively. By employing a harmonic approximation 30 , we treat the lattice dynamics classically. In this sense, its Hamiltonian assume the following form where P i is the momentum of the i-th site with mass M, and K is the force constant associated with the σ bond 30 . The electronic Hamiltonian, in turn, describes the π-electrons dynamics according to the equation below,  www.nature.com/scientificreports www.nature.com/scientificreports/ The summation runs over π-electrons in neighboring i and j sites with spin s (see Fig. 5). † C i s , and C i,s denote the creation and annihilation of an electron in states denoted by their subscript indices. To consider an external electric field ( → E ), we use a vector potential according to → = − E t c A t ( ) (1/ ) ( ). The exponentials come from the Peierls substitution method 34 . The unit vector r i j , points from j site to i site. Finally, the parameter γ in H elec is defined as γ ≡ ea c  , where a is the lattice parameter, e the fundamental charge, and c the speed of light. The term t i,j is the hopping integral, which couples the π-electrons to the lattice according to In Eq. 4, α is the electron-phonon coupling constant and η i,j is the relative displacement of the lattice sites from their equilibrium positions.
The dynamics calculation starts from an arbitrary initial set of coordinates {η i,j }, that is necessary to solve the electronic part of our model Hamiltonian initially. As a consequence, this procedure leads to an eigenvalue-eigenvector equation for the electronic component of the system, where the eigenvalues are E k and the eigenvectors are ψ k,s (i,t = 0). These quantities can be related as follows: k k s i j l s i j l s i j l s , ,, where i, j, j′ and j″ are neighboring sites.
To solve the classical component or our model, that describes the lattice structure, we turn to the Euler-Lagrange equation 23 . From the solution of the electronic part, we evaluate the expectation value of the wave function 〈Ψ|L|Ψ〉. This equation leads to: where www.nature.com/scientificreports www.nature.com/scientificreports/ couples the electronic and lattice degrees of freedom. The primed sum means that only the occupied states are considered.
The solution of the Euler-Lagrange equation with P i = 0 leads to a new set of coordinates {η i,j } that is used to recalculate the electronic Hamiltonian. This process is repeated iteratively until they reach the convergence criteria. As a result, this self-consistent procedure yields the ground state geometry that considers the interdependence between charge and lattice.
After achieving the convergence criteria, the time evolution of the initial state can be accomplished using the full Euler-Lagrange equation 23 . The time evolution of the electronic part is governed employing the time-dependent Schrödinger equation. To do so, we expand the wave function ψ k,s (t) in the basis of eigenstates of the electronic , at a given time t. Therefore, the wave function in time t + dt can be expressed as The dynamics of the electronic structure is carried out by using Eq. 8, that is evaluated numerically and then employed to the calculation of the expectation value of a new Lagrangian 23 . The Euler-Lagrange equation leads to a Newtonian type expression that takes into account the neighboring bonds: The applied electric field is turned on adiabatically, to avoid numerical error, in the following scheme: where t f is the total time of simulation and τ is the time needed until the electric field reaches its full strength.
To avoid edge effects, we consider periodic boundary conditions in the vertical direction of the nanoribbon, where the field is applied. Here, we use the notation NxM-AGNR, where N and M represent the number of sites on the horizontal and vertical directions of the nanoribbon, respectively. As all systems considered have M = 70, for the sake of simplicity, we use the notation N-AGNR. In the studied cases, N vary within the interval of 4-12.

Conclusions
In summary, charge carrier dynamics in AGNRs under the influence of an external electric field were analyzed employing a 2D generalization of the SSH Hamiltonian. AGNRs with widths ranging from N = 4 to N = 12 were studied. Results point to the polaron and bipolaron formation on such systems, where these quasiparticles respond differently to the external electric field, being the inertia of bipolarons larger. Eventually, however, both quasiparticles stop accelerating under the electric fields, moving afterward with constant velocity. Making use of a Stokes dissipation model, we were able to determine the effective mass of the charge carriers for several ribbon widths. It is shown that the effective mass of these quasiparticles varies drastically depending on two aspects, the system's width and the particular kind of quasiparticle present in the system. The effective mass for polarons presented values from 0.31 m e to 14.7 m e . In the case of bipolarons, the effective mass had values between 4 m e and 60 m e .