Tiago de Sousa Araújo
Cassiano
a,
Geraldo Magela
e Silva
*a and
Pedro Henrique de Oliveira
Neto
ab
aInstitute of Physics, University of Brasília, 70919-970, Brasília, Brazil. E-mail: magela@unb.br
bInternational Center of Physics, University of Brasília, 70919-970, Brasília, Brazil
First published on 26th February 2024
Organic thermoelectric devices allow the conversion of heat into electricity in a sustainable way, making them strong candidates to solve the present energy crisis. In this context, integrating graphene nanoribbons (GNRs) into thermoelectrics holds great potential for addressing this challenge. The development of a physical description of charge carriers under thermal influence is a paramount step toward this objective. However, to this day, the effects of temperature on charged quasiparticles hosted on GNRs remain elusive. In this work, we propose an adaptation to the long-established Su–Schrieffer–Heeger (SSH) model Hamiltonian to accommodate thermal effects on GNRs. The results show that random lattice deformations can significantly alter polarons’ and bipolarons’ localization profiles. Moreover, we report a thermally-induced re-balance of carrier stability. As temperature increases, the probability of observing bipolarons decays in favor of the formation of two independent polarons. The results are especially relevant to Seeback-based thermoelectrics because they rely on temperature gradients. With the thermal stability of charge carriers, local thermal differences could regulate GNR-based currents with quasiparticle population control.
Since the isolation of graphene about 20 years ago, 2D materials have been driving a revolution in materials science. Graphene derivatives play a special role in this development. Most of their appeal resorts to the emergence of unique quantum phenomena after simple structural/environmental changes.3–10 Nowadays, this attribute is the core of novel optoelectronic devices. Concerning semiconducting apparatuses, graphene nanoribbons (GNRs) stand out as the flagship choice in ultra-narrow applications. Among the many edge shapes,3,11 the armchair type (AGNR) has an encouraging prospect, already integrating novel devices.5,11,12 Given this scenario, using GNRs in next-generation thermoelectric devices is a natural path to proceed.13
Charge transport is fundamental for the description of the thermoelectric effect. In organic materials, the mechanism is driven by charged quasiparticles. They are polarized defects coupled with phonon clouds resulting from the collective response of the system.14,15 The polaron, a half-spin structure of charge e, is usually the main carrier during transport. Besides its importance, other high-order quasiparticles can emerge during the transport mechanism. In a highly doped configuration,16,17 two independent polarons can merge together, recombining into a new carrier known as a bipolaron.18 This bosonic quasiparticle has twice the polaron's charge. Naturally, it is an active carrier with its own transport and magnetic properties. Importantly, bipolarons also participate in photonic mechanisms19 because the scattering process may induce the formation of excitonic quasiparticles.16,18,20 Therefore, the stability of polarons and bipolarons reverberates in electronic and optical performance.
Importantly, due to the coupled nature of a charged quasiparticle, its shape and stability directly impact the mobility and effective mass. Therefore, characterizing charge carriers allows to predict performance capabilities and provides new insights into transport mechanisms. Such a strategy had previously led to important results for conjugated polymers,21–24 GNRs,6,25–27 and molecular crystals.28,29
However, describing the charge carriers is not enough to understand thermoelectric effects. A realistic physical description requires an explicit consideration of thermal response. Substantial effort was directed to developing adequate strategies.30 The lattice can be highly sensitive to thermal influence, making the use of hybrid Hamiltonians, with explicit electronic and lattice interactions, a desirable approach.30 The matter is especially relevant when considering systems with a strong electron–phonon coupling. In a hybrid Hamiltonian approach, one can model thermal effects as random distortions of the lattice bonds.31–33 In this context, the Su–Schrieffer–Heeger hybrid Hamiltonian15,34 is a robust framework to investigate thermal mechanisms in quantum systems, being the backbone of important studies.31,32,35–37
Progress was made through this route, which helped uncover the functioning of intricate phenomena in highly confined materials. For instance, with a mixed 1D Hamiltonian, Troisi recovered the long-established mobility power law dependence with temperature in rubrene crystals.28,29 Other studies reported a thermal-assisted reduction of electron–hole pair binding energy,32 the hampering of polaron stability,33,36 and shortening of quasiparticle recombination37 in conjugated polymers. These works demonstrate that the theoretical investigation of thermal properties holds a latent potential to elucidate and provide insights in nanoelectronics and photonics. However, besides the importance of GNRs, the impact of thermal influence over its charge carriers remains an open topic.
In this work, we simulate the formation of polarons and bipolarons on armchair graphene nanoribbons under several temperature regimes. The model consists of an adaptation of the SSH hybrid Hamiltonian. The proposed modification aims to consider thermal effects as bond distortions. We selected four armchair nanoribbons with different widths to implement the study. Bipolaron stability is determined. The results show that raising the temperature can progressively hamper bipolaron cohesion, favoring the formation of independent polaron pairs instead. Additionally, we report an increase of 25% of polarons for a temperature variation of 150 K. The results could benefit future thermoelectric applications based on the Seeback effect while providing a phyisically accurate description of charge carriers under thermal influence.
H = Hlatt + Helec | (1) |
The former is modeled under the Harmonic approximation, reading8,38,39
(2) |
Here, K is the harmonic constant, M is the site mass, Pi is the conjugated momentum of the i-th site, and ηi,j is the deformation of the bond connecting sites i and j. The bracket at the lower sum index indicates a summation over neighboring bonds.
The electronic part is treated through a tight-binding formalism, which models the π-electron hopping as38
(3) |
The electronic and lattice interactions are coupled through a first-order expansion of the hopping integral in the pristine configuration (t0):15,21,34
ti,j = (t0 − αηi,j). | (4) |
The lattice coupling to the hopping integral is determined by the electron–phonon coupling α.
Following the Ehrenfest Theorem, the lattice solution can be obtained by solving the Euler–Lagrange equations of the expected value of the Lagrangian.8,21,38 In the case of stationary configurations, when no kinetic interaction is considered (Pi = 0), the minimization condition simplifies to:
(5) |
Here, is the expected value of the potential energy. The state is a Slater determinant formed by the occupied π-electron orbitals. After explicit manipulation, this expression turns into38
(6) |
(7) |
Since the Hamiltonian explicitly depends on the lattice configuration, a self-consistent algorithm is applied to solve the hybrid model.8 The procedure begins by choosing an appropriate guess for the bond distortion set {ηi,j} (which is usually taken as 0). Then, Helec can be numerically built and diagonalized, yielding the corresponding eigenvectors and eigenvalues. The new electronic solution is then applied to recalculate {ηi,j} according to eqn (6). Then, the newly calculated lattice configuration is compared with the previous set. This involves calculating the root mean squared error between the two solutions. If the difference is below a preset threshold, the algorithm stops. Otherwise, the procedure repeats by using the updated {ηi,j} in the initial step.
Charge states are obtained by removing the electrons from the material. This means extracting spin–orbitals from the neutral state Slater determinant. For instance, a positive polaron is created after removing an electron occupying the highest occupied molecular orbital (HOMO). In turn, a positive bipolaron emerges after removing two electrons from the HOMO. The corresponding carriers with a negative charge are obtainable by adding electrons to the unoccupied molecular orbital (LUMO). Importantly, the positive and negative quasiparticles are bound to exhibit the same properties due to the electron–hole symmetry.
The parameters of the model Hamiltonian are based on previous studies on AGNRs and similar conjugated systems. That is t0 = 2.7 eV,4,40K = 21 eV Å−2,34,41 and α = 5.2 eV Å−1.26,39,42,43
Xki,j(T) = Xi,j + Δki,j(T), | (8) |
However, it is always possible to recognize ηi,j into Xi,j because they differ only by a constant. In other words,
Xi,j = L0 + ηi,j, | (9) |
Because of that, the lattice changes due to thermal effects can be thought as corrections to the bond length distortion, turning each ηi,j into
ηki,j(T) = ηi,j + Δki,j(T). | (10) |
The Δki,j(T) is sampled according to the Boltzmann probability distribution:30,32,33
(11) |
Here, kb is the Boltzmann constant and x is the distortion random variable. For each ensemble member, the sampling is done on all existing bonds. The process consists of two steps: (1) for each bond, an x value is sampled from the distribution in eqn (11) through the Box–Muller algorithm; (2) the distortions are implemented in the lattice, rendering a new set {η}. Then, Helec is diagonalized using the updated lattice, yielding the corrected eigenvalues and eigenvectors of this particular configuration. Finally, ensemble statistics will provide the physical description of the GNRs. As a side remark, we recall that the normalization of the distribution is irrelevant for the sampling because it is a common factor for all x at a fixed T.
Fig. 1(d) and (e) illustrate the sampling procedure. In Fig. 1(d), the converged bond length Xi,j connects the sites i and j. For a finite temperature T, Xi,j is corrected by a random distortion Δi,j that is sampled from the distribution in eqn (11). The result is a bond distorted from its equilibrium length, as shown in Fig. 1(e). This process is repeated for all unique bonds, giving one ensemble member as the end result. By repeating the procedure, an ensemble of GNR configurations is generated. Their collective response will translate into the physical description of the material under thermal effects.
All cases were simulated by constructing ensembles with 5000 members. This number was found to be adequate after verifying the statistical robustness of averaging processes. Fig. S1 (ESI†) illustrates this test. The variation of the total energy in the neutral, polaron and bipolaron states is presented in Fig. S1(a), (c) and (d) (ESI†). The solid lines represent the energy of a given ensemble element, while the dashed lines are the corresponding cumulative averages. As can be seen, the solid curves display an oscillating behavior. On the other hand, the dashed lines are stable, exhibiting no significant changes when the ensemble has more than 1000 elements. Therefore, any succeeding quantity that relies solely on the total energies should be equally stable.
We emphasize that a small ensemble can lead to an underrepresented thermal response. To see this, we present the zoomed versions of Fig. S1(a), (c) and (d) (ESI†), respectively, in Fig. S1(b), (d) and (e) (ESI†). As can be seen, the dashed curves are not straight lines anymore. For a small ensemble, the cumulative mean is not stable.
In the following discussion, the k index will be suppressed when its presence is not relevant. We emphasize that ordering the ensemble geometries has no physical meaning and is only important for the sake of clarity and practical reasons. Finally, all the simulated results obtained using the traditional model (with no explicit consideration of thermal effects) will be referred to as the “0 K” case.
In addition, we recall that in general, thermal effects influence the electrons through changes in the electronic occupation. Since electrons are fermionic particles, the Fermi–Dirac distribution determines the temperature dependence of the occupation. This effect was not included because all the simulated GNRs are intrinsic semiconductors, with a gap of ≈1–3 eV. In such cases, a significant modification in the electronic occupation is only visible for thermal baths at high temperatures ≈5000 K. Since we are concerned with near-ambient temperatures, the occupation change is negligible. Naturally, the situation can drastically change when dealing with semimetals, metals, and highly doped semiconductors. For those materials, even a low-temperature reservoir can thermally excite electrons to the conduction band. The impact of thermal effects on the electronic occupation depends on the actual density of charge carriers. In a high-density charge carrier regime, the system has many particles at the conductance (or valence) band. In such cases, occupation changes become relevant because there is a large number of particles to fill the vacant states. On the other hand, the intra-band excitation, although present, has a minimal impact. This is because it alters the occupation of only a few electrons (or holes).
Note that this does not mean that electrons are insensitive to thermal effects in the model. We recall that the electronic Hamiltonian has a contribution from the lattice through the electron–phonon coupling. Therefore, the thermal distortions in the sites’ positions reverberate due to the electronic eigenstates and eigenvalues.
BE = FEbip − 2FEpol, | (12) |
FEbip = SCFbip − SCFneutr and |
FEpol = SCFpol − SCFneutr. | (13) |
Here, SCFbip/pol is the self-consistent energy obtained from the stationary algorithm for the bipolaron/polaron state. In a similar manner, SCFneutr is the corresponding energy for the neutral configuration. If BE > 0, the formation of two independent polarons is energetically more favorable than a single bipolaron. In other words, bipolaron formation would occur only with the aid of an external agent. When BE < 0, the opposite becomes true: bipolarons are expected to form spontaneously.
It is important to note that SCFbip, SCFpol, and SCFneutr are distributions when considering thermal effects. Each state is represented by an ensemble of 5000 separate calculations. For this reason, eqn (12) needs to be adapted. First, it is important to realize that each ensemble member shares the same probability of being observable. In addition, BE should be a distribution that obeys eqn (12). Therefore, BE calculation is done by the sampling values obtained from FEbip and FEpol and applying them to eqn (12). However, although independent in principle, the polaron and bipolaron formation energies cannot be sampled separately for BE calculation. During each selection process, the same SCFneutr sampled energy was used to calculate both FEbip and FEpol. That way, we avoid the consideration of artificial recombinations involving ground-state lattice reorganization.
(14) |
It is also suggested to look at the order parameter heatmap of ensemble members. Fig. 2(d)–(f) correspond to 0 K, 95 K and 350 K, respectively. For all cases, dark tones indicate a bond expansion, while light colors represent a contraction. Bonds in red display no distortion. When temperature is not considered, the 4AGNR bonds form an alternating pattern, intercalating between contraction and expansion along the width lines. This is the dimerization effect, which is widely reported for GNRs8,38 and other conjugated materials such as polymers.15,34 As the temperature increases, this pattern is progressively destroyed. Some bonds that were once contracting at 0 K start expanding now and vice versa. Another possibility is for the bond to display no distortion at all, indicated by the color red. The results show that thermal effects can significantly change the lattice. Therefore, any mechanism that rely, to some extent, on lattice properties is expected to be also influenced. This is especially relevant for quasiparticle formation that usually depends on the lattice symmetry for the generation of stable carriers. Therefore, it is expected that carriers created at high-temperature regimes are less stable than the 0 K case.
One may notice that the profile for 197 K is considerably closer to the 350 K configuration than it is from the 95 K configuration. In other words, the thermal response in a high T regime is expected to be less sensitive to temperature variations. This trend will be visible in all the following results.
Fig. S3–S5 (ESI†) present the same panel of Fig. 3 for the 6, 7, and 9AGNR. As in the 4AGNR case, these nanoribbon bonds display an increasing degree of expansion or shrinking as the temperature rises. Moreover, wider nanoribbons present bond length distributions with peaks closer to 0 Å. Because of that, thermal distortions in wider nanoribbons are more prone to change the bonding state (if it is expanding or shrinking). Therefore, the result foreshadows an enhanced sensitivity of wider nanoribbons to dimerization changes.
To further explore this indication, we proceed to analyze the bond distributions more systematically. Here, we present Fig. 3(a)–(d). They are the corresponding heatmaps of distortion distribution of all bond types combined for the 4, 6, 7, and 9AGNR. The yellow tones indicate a high probability of occurrence, while purple represents a low probability. A white line is drawn at η = 0 Å to serve as a guide. Each different nanoribbon has at least two high-probability regions. The first one corresponds to the expanded bonds (η > 0) and the other represents contraction (η < 0). The position of these zones along the η axis depends on the nanoribbon width. For instance, when T ≈ 95 K, the 4AGNR heatmap has an orange region near η = 0.06 Å, while the 6AGNR's corresponding peak is at η = 0.09 Å. In all cases, increasing the temperature seems to spread the high-probability regions over the η-axis, a response already observed in Fig. 2(a)–(c).
When T is high enough, the peaks that once were related to contracted bonds begin to overlap with the expansion peaks for all GNRs. As a result, the bond distortion sign begins to flip in some configurations. To better visualize the effect, we present DR as a function of the temperature for the 4 (blue), 6 (red), 7 (yellow), and 9AGNR (green) as shown in Fig. 3(e). Qualitatively, all curves share the same behavior: as T increases, DR decays monotonically. This is an expected outcome. It results from the growing overlap between the distributions’ peaks associated with contracted and expanded bonds.
Some nanoribbons are more prone to dimerization degradation than others. At 95 K, the 4AGNR has DR ≈ 96%, bearing almost no change compared with the 0 K solution. On the other hand, the 9AGNR has a DR of ≈62% for the same temperature. This is a result of the initial position of the η peaks shown in Fig. 2(a)–(d). The peaks closer to 0 Å will cross this value at lower temperatures, providing a non-zero probability to the flipping of the distortion sign of some bonds. The wider nanoribbons present bonds closer to this region. Therefore, the results suggest that width extension enhances the temperature changes in the nanoribbon structural stability. Since quasiparticle formation depends on lattice cohesion, charge carriers hosted on wider nanoribbons are expected to be less stable.
Although wider GNRs suffer more from temperature effects, the relation is not monotonic. For instance, the 6AGNR suffers a stronger dimerization degradation than the 7AGNR. That is due to the grouping of armchair families.10,46 Members belonging to the same family are expected to group together, making width effects to be non-monotonous. Because of that, the 4AGNR and 7AGNR have to display a similar response. The same applies to the 6AGNR and the 9AGNR. A clear example of such behavior is the dependence of the AGNR energy bandgap on the nanoribbon's width, as widely reported the in literature.10,46
Although temperature significantly impacts the dimerization of individual members of the ensemble, we observe no effects when considering the mean values of η. Fig. S1 (ESI†) displays this calculation for the 4AGNR at 300 K. Four ensemble individuals are displayed in Fig. S1(a) (ESI†). The average of them, along with the remaining 4996 units of the set, yields the heatmap shown in Fig. S1(b) (ESI†). This profile is equivalent to the one shown in Fig. 2(d). The explanation for this lies in the Gaussian form of the probability distribution used to sample the distortions. The average of Δki,j(T) is zero for any T because eqn (11) is symmetric. Therefore, summing out the contribution from each ensemble member will just set each bond distortion to its original value from the stationary algorithm, ηi,j.
Provided the neutral state's general response, we proceed to analyze the effects on the charged states. Fig. 4(a) presents the charge density (ρ) for the 7AGNR. Four heatmaps are displayed, corresponding to individuals from ensembles of different temperatures. Hot colors signal high charge accumulation, while cold tones indicate a low local density. As can be seen, the carrier morphology suffers substantial changes due to temperature. The symmetry patterns present in the 0 K polaron are often eroded due to the random distortions. For instance, when T = 0 K, the amount of charge on the left side equals the right side. However, this symmetry is lost when considering the 197 K case. Here, the right side has slightly more charge. Changes in the quasiparticle size are observed too. The polarized region at 350 K extends over ≈30 Å, while the 0 K structure has a size of ≈25 Å. Similar asymmetries can be spotted for the other nanoribbons (not shown). These findings are consistent with previous studies on conjugated polymers, in which the polaronic picture of the carrier is challenged with increasing temperature.33
A complete overview of a charged quasiparticle requires examination of the nanoribbon lattice deformation profile. The order parameter heatmap is displayed in Fig. 4(b). At 0 K, a distinguishable local distortion is present at the center of the heatmap, representing the polaron local lattice deformation. Due to the coupling of carriers, the deformation size, shape, and magnitude directly impact the resulting charge accumulation.6,26,43,47 Importantly, the lattice disarray can induce asymmetries in the charge profile. Previous works regarding charge profiles on conjugated polymers indicate that cohesion loss can potentially hinder the transport capability.48 Here, we extend this reasoning by concluding that the thermally-induced random distortions may reduce the carrier's mobility. This is an expected outcome, corroborated by widely reported mobility degradation with increasing temperature in molecular crystals.28,31,49 Confirmation of this response is achievable through dynamical simulations.
We remark that the profiles displayed in Fig. 4 cannot be interpreted as the average behavior of the nanoribbon in that temperature regime. These profiles are snapshots of individual ensemble members. Therefore, a comparison between these structures holds no physical significance. Only the collective response of the ensemble set has this meaning. The heatmaps are presented to illustrate how distorted the 0 K polaron picture becomes when temperature effects emerge.
Now, we move to quantitative analysis. Fig. 5(a)–(d) display the formation energy distributions of the 4AGNR at 100, 150, 200, and 350 K. The curves in red refer to the polaron state, while the ones in blue represent the bipolaron response. Regardless of the temperature magnitude, each state distribution exhibits a Gaussian-like shape with a well-defined peak. Moreover, the peak position of the polaron case is always at a lower energy than that of the bipolaron. This is an expected outcome, reflecting that the system requires more energy to create a bipolaron than a polaron. That should be the case because the additional charge requires a more profound lattice reorganization to form a localized polarized structure. However, as the temperature rises, an interesting effect begins to unfold: the distributions overlap. In other words, certain polaron states, assisted by thermal effects, are energetically more costly to form than bipolaronic configurations. This nonintuitive setting suggests that quasiparticle stability can be widely affected by thermally-induced distortions. The overlapping feature is present in all simulated GNRs.
Fig. 5 4AGNR probability density distribution of the polaron (red) and bipolaron (blue) formation energies. (a)–(d) Refer to, respectively, T = 95, 197, 299, and 350 cases. |
The variants shown in Fig. 5 for the 6, 7, and 9AGNR are presented in Fig. S6–S8 (ESI†). Qualitatively, all nanoribbons display the same behavior: the polaron and bipolaron distributions are monomodal curves that begin to overlap when increasing the temperature. However, close analysis shows that the overlap degree depends on the width. This can be verified by comparing Fig. 5(d) with Fig. S8(d) (ESI†). We recall that the dispersion of the distribution is associated with the material thermal response. The broader the curve, the more intense the change induced by the distortions. Therefore, the result evidences that wider nanoribbons are more sensitive to thermal influence.
For some configurations, EFbip < EFpol. Then, one should expect that the bipolaron stability distribution (see eqn (12)) may have a non-zero probability of being positive. In other words, a non-zero temperature favors the destabilization of bipolarons, effectively facilitating polaron formation. To verify this reasoning, we present the bipolaron stability heatmaps in Fig. 6(a)–(d) for the 4, 6, 7, and 9AGNR, respectively. The colors indicate the probability density magnitude: the zones in blue have a high probability density. In turn, the BE regions in yellow have a low probability density. For each heatmap, a red line is drawn at BE = 0 eV to serve as a reference. Near this line, all distributions display a single strong blue signal. Then, the density progressively drops for BE values away from the central region. This behavior is characteristic of a monomodal distribution. The FE curves shown in Fig. 5 display a similar response.
Fig. 6 (a)–(d) depict the bipolaron binding energy distributions heatmaps for the 4, 6, 7, and 9AGNR, respectively. The regions in blue present high probability density, while those in yellow have a low probability. (e) displays the percentage of probability of stable polarons as a function of the temperature. This quantity is obtained by integrating each BE distribution in the region BE > 0 eV. The colors of the curves indicate the nanoribbon width, following the same convention shown in Fig. 3. |
Note that the center of the distributions is not at BE = 0 eV. In the absence of thermal effects, BE equals −0.36, −0.19, −0.23, and −0.10 eV for the 4, 6, 7, and 9AGNR, respectively. Therefore, bipolarons are more stable than two independent polarons when T = 0 K. As the temperature rises, the distribution expands along the y-axis, making high BE energy states accessible. A consequence of this feature is that the distributions begin to cross the red line, exhibiting a non-zero probability of occurring positive BE configurations.
That thermally-induced behavior has the potential to represent significant implications for GNR-based nanoelectronics. We recall that bipolarons are active agents in photonics. They are involved in excitonic formations, influencing the emission mechanism. Since temperature destabilizes bipolaron formation, thermal effects will reverberate due to the optical properties of GNRs. If bipolarons are prone to split into polarons, biexciton formation might be hampered. As a result, characteristic biexciton peaks at emission spectra are expected to drop with temperature increase. These findings could be explored in thermoelectric applications. Seeback-based gadgets make use of spatial temperature gradients to function. Then, the results suggest that the ratio of polaron/bipolaron populations will be different at the device's edges. Such configuration could be rationalized as a spin filter tool because polarons are fermionic carriers and bipolarons are bosons.
It is important to note that the current model does not include electron–electron repulsion terms. While these terms are pertinent in the context of collision-induced recombination50,51 and magnetic processes,52–54 their incorporation beyond these scenarios does not yield substantial enhancements for the model. The effect is to cause minor adjustments to the gap magnitude and quasiparticle profile. In the present case, the quasiparticles are isolated. Therefore, this addition can be avoided.
Moreover, the rising probability density for positive BE is bound to change charge transport efficiency. Bipolarons often present different transport properties compared with polarons.6,44,47 Our results show that temperature favors the formation of polarons. Therefore, the polaron population is expected to grow with T. As a result, the nanoribbon charge transport mechanism could be characterized by polarons and their properties even in regimes of strong doping. This feature can affect future thermoelectric applications in which transport efficiency (i.e. the current magnitude) is temperature-dependent.
The degree of bipolaron destabilization depends on the nanoribbon width. To investigate this in detail, Fig. 6(e) presents the probability of obtaining two stable polarons as a function of the temperature. The curves were obtained by integrating the probability density curve inside the positive BE region. The color of the curves indicates the nanoribbon width, following the same convention shown in Fig. 3(e). As expected, increasing the temperature raises the probability of positive BE for all AGNRs. However, the nanoribbon width has a profound impact on the probability value. For instance, the 4AGNR has a 26% chance of having a configuration with two stable polarons at 95 K. When T = 350 K, this probability goes to 45%. At the same temperature, in the case of the 9AGNR, the probability escalates from 45% to 49%. Another point worth mentioning is that due to the monomodal character of BE distribution, the probability of positive BE cannot be greater than 50%. The reason is that the probability density also spreads towards the negative BE region. Consequently, at least half of the distribution will still have BE < 0. This interesting result shows that one cannot indefinitely favor the bipolaron split through temperature increments.
When it comes to charged states, thermal effects induce significant changes in the formation energies of polarons and bipolarons. As an outcome of ensemble statistics, these quantities become distributions. Their dispersion increases with rising temperature. As a result, an interesting result emerges for high T: bipolaron and polaron FE probability density distributions overlap. This means that certain configurations display a lower bipolaron formation energy compared with a single polaron. Further investigation of these cases led to the calculation of bipolaron binding energy distribution. The results show that temperature promotes the occurrence of states with BE > 0 eV, effectively favoring the formation of polarons instead of bipolarons. In contrast to 0 K simulations in which bipolarons are stable, here we report that almost 50% of the ensemble states contain stable polarons for high T. These unexpected results paint an encouraging prospect to envision thermoelectric devices that control the population of charge carriers. Moreover, the results illustrate the importance of including thermal effects to represent realistic scenarios.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ma01181j |
This journal is © The Royal Society of Chemistry 2024 |