Topological interface states and nonlinear thermoelectric performance in armchair graphene nanoribbon heterostructures

Abstract

We investigate the emergence and topological nature of interface states (IFs) in N-AGNR/(N − 2)-AGNR/N-AGNR heterostructure (AGNRH) segments lacking translational symmetry, focusing on their relation to the end states (ESs) of the constituent armchair graphene nanoribbon (AGNR) segments. For AGNRs with R₁-type unit cells, the ES numbers under a longitudinal electric field follow the relations N = N_A(B) × 6 + 1 and N = N_A(B) × 6 + 3, whereas R₂-type unit cells exhibit (N_A(B) + 1) ESs. The subscripts A and B denote the chirality types of the ESs. The Stark effect lifts ES degeneracy and enables clear spectral separation between ESs and IFs. Using a real-space bulk boundary perturbation approach, we show that opposite-chirality states hybridize through junction-site perturbations and may shift out of the bulk gap. The number and chirality of IFs in symmetric AGNRHs are determined by the difference between the ESs of the outer and central segments, N_O and N_C, according to N_IF,β = |N_O,B(A) − N_C,A(B)|, where β labels the chirality. Depending on whether N_O > N_C or N_C > N_O, the resulting IFs acquire B- or A-chirality, respectively. Calculated transmission spectra reveal that AGNRHs host a topological double quantum dot (TDQD) when IFs originate from the ESs of the central AGNR segment. Using an Anderson model with effective intra-dot and inter-dot Coulomb interactions, we derive an analytical expression for the tunneling current through the TDQD via a closed-form transmission coefficient. Thermoelectric analysis shows that TDQDs yield enhanced nonlinear power output in the electron-dilute and hole-dilute charge states, with Coulomb blockade suppressing thermal current but not thermal voltage. The thermal power output of the TDQD is significantly enhanced by nonlinear effects, even under strong electron Coulomb interactions.

---

1. Introduction

Since the groundbreaking discovery of graphene in 2004,¹ extensive experimental and theoretical studies have focused on graphene nanoribbons (GNRs).^2–13 Advances in bottom-up fabrication now allow atomically precise GNRs with diverse geometries.^2–13 Among these, armchair GNRs (AGNRs) exhibit width-dependent electronic structures,^2–6 enabling tunable semiconducting phases suitable for quantum device applications. Of particular interest are zero-dimensional topological states (0D-TSs) associated with the interface states (IFs) of AGNR heterostructures (AGNRHs), which have been observed in scanning tunneling microscopy measurements of the local electronic density.^7–11 These 0D-TSs emerge within the mid-gap region of semiconducting AGNRHs,^10,11 providing a promising route toward atomically precise topological quantum dot (TQD) devices.^14–20 Their robustness allows the design of double TQDs and TQD arrays with controllable electron hopping and Coulomb interactions between TSs.¹¹

Topological states in GNRs were originally predicted using the Zak number Z₂ of bulk GNRs,²¹ which requires time-reversal and translational symmetries.^21–25 However, experimentally synthesized 9–7–9 and 7–9–7 AGNRHs lack translational symmetry, as illustrated in Fig. 1(a). Because the wave-function decay length of the TS is very short, 9–7–9 AGNRH segments do not exhibit significant size effects on the TSs,¹¹ which explains why predictions based on the Zak number remain valid for 9–7–9 and 7–9–7 AGNRH segments.²⁴ For wider AGNRHs, such as 15–13–15, 21–19–21, and 27–25–27 segments, the wave-function decay lengths of TSs increase with width, making Zak number calculations insufficient. This motivates a real-space analysis of the relationship between ESs of wider AGNR segments possessing multiple end states (ESs) and IFs in wider AGNRHs, which is critical for designing two-dimensional TS-based crystals for novel quantum devices.²⁶

Fig. 1 (a) Schematic illustration of a 9_w–7_x–9_y armchair graphene nanoribbon heterostructure (AGNRH) composed of three AGNR segments. The size of the AGNRH is characterized by (N, M), where N and M denote the row and column numbers. R₁ and R₂ represent the unit cells (u.c.) of the 9-AGNR and 7-AGNR segments, respectively. (b) A and B sublattice sites are indicated by white and blue colors. The inter-AGNR electron hopping strengths t_es connect adjacent atoms between the 9-AGNR and 7-AGNR segments. (c) AGNRH segment with zigzag edge terminations coupled to left and right electrodes with equilibrium temperatures T_L and T_R. Γ_L(R) denote the tunneling rates for electrons tunneling from the left (right) electrode into the adjacent atoms at the zigzag edges.

A real-space approach involves solving the Schrodinger equation to obtain eigenvalues and eigenfunctions of AGNRHs. Exact analytical solutions exist only for triangular or rectangular GNR structures within the one-band tight-binding model.^27–32 Numerical methods can alternatively determine the number of IFs in AGNR heterojunctions.^33–43 However, distinguishing the locations and counts of IFs in numerical calculations is challenging due to their zero-energy modes. To overcome this, we propose applying a longitudinal electric field to AGNRH segments, as illustrated in Fig. 1(b). The resulting Stark effect not only resolves the number of IFs but also clarifies their spatial locations.

This study has two main objectives. First, we aim to elucidate the correlation between the ESs of isolated AGNR segments and IFs in N-AGNR/(N − 2)-AGNR/N-AGNR heterostructures using a real-space formulation. While only 9–7–9 and 7–9–7 AGNRHs have been experimentally realized via bottom-up synthesis,¹¹ we extend our analysis to semiconducting AGNRs in the N = 9, 15, 21, and 27 families. To explore the evolution of terminal states of isolated AGNRs into IFs of AGNRHs, we introduce a tunable inter-AGNR hopping parameter t_es at the junction sites (Fig. 1(b)), following the bulk-boundary perturbation method.⁴³ When t_es = 0, the AGNRH decouples into three independent segments; varying t_es effectively modifies the boundary conditions of each segment. Sublattices are labeled A and B in Fig. 1(b) to emphasize the chiral symmetry of the structure.²⁴

Second, we analyze charge transport through IFs in the Coulomb blockade regime when zigzag-terminated ends of the AGNRH are connected to electrodes (Fig. 1(c)). Tunneling current measurements probe the spatial distribution of localized states, providing essential information for designing one- and two-dimensional 0D-TS crystals. We also tackle the challenging calculation of nonlinear thermoelectric power in topological double quantum dots (TDQDs) formed by AGNRH segments. By employing an analytical solution for tunneling current in the Coulomb blockade regime, we evaluate the influence of Coulomb interactions on both thermal voltage and thermal current, providing a comprehensive understanding of TDQD thermoelectric performance.

2. Calculation methodology

The system Hamiltonian, illustrated in Fig. 1(c), is written as: H = H₀ + H_GNR, where H₀ describes the Hamiltonian of the electrodes and H_GNR represents the AGNRH. The Hamiltonian of H₀ is given bywhere the first two terms describe the free electrons in the left and right electrodes. The operators a^†_k,σ (b^†_k,σ) create an electron with momentum k and spin σ in the left (right) electrode, each with energy ε_k. The terms and describe the coupling between the electrode and its adjacent atoms in the -th row of the AGNRH. The electronic state of the GNR is described using a tight-binding model with one p_z orbital per carbon atom.^44–46 The AGNRH Hamiltonian H_GNR is written as:where denotes the on-site energy for the p_z orbital at row and column j. Spin–orbit interaction is neglected in this model. Graphene possesses exceptionally weak spin orbit coupling and negligible hyperfine interaction, owing to the dominant presence of ¹²C atoms with zero nuclear spin.^47–51

The operators and create and annihilate electrons at the site , respectively. The hopping integral describes electron transfer between sites and . The tight-binding parameters used for the AGNRH are for all sites and eV for nearest -neighbor hopping. A perturbative hopping term t_es is introduced to probe the interaction between pairs of states with opposite chirality, as illustrated in Fig. 1(b). The effect of a longitudinal electric field F_y is incorporated through an additional potential U_y = eF_yy added to , where F_y = V_y/L_a, with V_y being the applied bias and L_a the length of the AGNRH.

To investigate the electron transport through the AGNRH junction, the transmission coefficient is evaluated using , where Γ_L(R)(ε) denote the tunneling rates of the left (right) electrode, and G^r(a)(ε) are the retarded (advanced) Green functions of the AGNRH.^52–56 In the tight-binding formulation, Γ_α(ε) and Green's functions are matrices. The expression for Γ_L(R)(ε) is derived from the imaginary part of the self-energies, denoted as , and is given by . For simplicity, we adopt the wide-band limitation and assume Γ_L(R)(ε) to be energy-independent, denoted simply as Γ_L(R).

3. Results and discussion

3.1. End states of AGNR segments and interface states of AGNRH segments

Because the end states (ESs) of an AGNR segment possess zero-energy modes, it is difficult to resolve their degeneracies through numerical calculations alone. Meanwhile, to experimentally reveal the number of ESs, we apply a Stark effect induced by an external electric field to lift the zero-energy modes and thereby determine the number of ESs in AGNR segments with an R₁ unit cell (u.c.), as shown in Fig. 1(a). Fig. 2 displays the energy levels of N-AGNR segments under a uniform electric field applied along the armchair direction (defined as the y-direction). The conduction- and valence subband edge states, E_c and E_v, exhibit red Stark shifts, while the localized ESs exhibit blue Stark shifts. These localized states show a linear dependence on the applied voltage V_y. The subband gap Δ_g = E_c − E_v decreases with increasing ribbon width for a family of AGNR, as illustrated in Fig. 2.

Fig. 2 Energy levels of armchair graphene nanoribbon (AGNR) segments as functions of the applied voltage V_y. (a) 13-AGNR, (b) 15-AGNR, (c) 19-AGNR, (d) 21-AGNR, (e) 25-AGNR, and (f) 27-AGNR segments. All segments have a length specified by M = 96, and their widths are characterized by N.

Furthermore, the number of ESs differs among AGNR segments of various widths. We find that 13-AGNR and 15-AGNR segments possess two left ESs (N_A = 2) and two right ESs (N_B = 2). The energies of the left ESs are labeled as Σ_c1 and Σ_c2, while those of the right ESs are labeled as Σ_v1 and Σ_v2. The subscripts c and v denote states above and below the charge-neutral point (CNP), respectively. These multiple ESs contrast with those in narrower segments such as 7-AGNR and 9-AGNR, which each exhibit only one left ES and one right ES (N_A = N_B = 1). For wider AGNRs, such as 19-AGNR and 21-AGNR, we obtain N_A = N_B = 3, and for 25-AGNR and 27-AGNR, we obtain N_A = N_B = 4. Based on the numerical results in Fig. 2, the ES numbers follow the rule N = 6N_A(B) + 1 and N = 6N_A(B) + 3. The former corresponds to semiconducting AGNRs with N = 3p + 1 and the latter to semiconducting AGNRs with N = 3p, where p is an integer.

The Stark shift of ES energy levels in AGNR segments can also be applied to AGNRH segments. For the N-AGNR/(N − 2)-AGNR/N-AGNR heterostructures considered in this work, we focus on the semiconducting family with N = 9, 15, 21, and 27. The outer AGNR segments in these heterostructures belong to the N = 3p family, while the central segment belongs to the N = 3p + 1 family. Fig. 3 presents the energy spectra of these AGNRHs as functions of the applied voltage V_y.

Fig. 3 Energy levels of armchair graphene nanoribbon heterostructure (AGNRH) segments as functions of the applied voltage V_y for different widths. (a) 9–7–9 AGNRH, (b) 15–13–15 AGNRH, (c) 21–19–21 AGNRH, and (d) 27–25–27 AGNRH. Here, we adopt w = y = 4 for the outer AGNR segments with R₁-type unit cells (u.c.) and x = 16 for the central AGNR segment with an R₂-type unit cell (u.c.).

As shown in Fig. 3(a), the 9–7–9 AGNRH exhibits two low-energy modes above the CNP and two below it. The modes labeled Σ_c,1 and Σ_v,1 correspond to the left and right ESs of the 9–7–9 AGNRH segment. The modes Σ_IF,c and Σ_IF,v represent the left and right topological symmetry-protected interface states (IFs) at the 9–7 and 7–9 AGNR heterojunctions. At V_y = 0.18 V, the energy levels are Σ_c,1 = 0.087 eV, Σ_v,1 = −0.087 eV, Σ_IF,c = 0.062 eV, and Σ_IF,v = −0.062 eV.

For the wider AGNRH structures shown in Fig. 3(b)–(d), only one left IF and one right IF appear, even though the AGNR segments contain multiple ESs. For example, in Fig. 3(b), the 15–13–15 AGNRH at V_y = 0.18 V exhibits six in-gap states: four ESs with energies Σ_c,2 = 0.08932 eV, Σ_v,2 = −0.08932 eV, Σ_c,1 = 0.08483 eV, and Σ_v,1 = −0.08483 eV, and two IFs with energies Σ_IF,c = 0.05972 eV and Σ_IF,v = −0.05972 eV.

Similarly, the 21–19–21 AGNRH in Fig. 3(c) exhibits eight in-gap states at V_y = 0.18 V: six ESs with energies Σ_c,3 = 0.08966 eV, Σ_v,3 = −0.08966 eV, Σ_c,2 = 0.08823 eV, Σ_v,2 = −0.08823 eV, Σ_c,1 = 0.0824 eV, and Σ_v,1 = −0.0824 eV, along with two IFs at Σ_IF,c = 0.05728 eV and Σ_IF,v = −0.05728 eV. The 27–25–27 AGNRH contains ten in-gap states (eight ESs and two IFs). Compared with the in-gap states of the isolated N-AGNR segments shown in Fig. 2, the N-AGNR/(N − 2)-AGNR/N-AGNR heterostructures exhibit two additional in-gap states arising from the IFs. Thus, the Stark effect reveals not only the number of ESs in AGNR segments but also the number of IFs in AGNRH structures.

3.2. Bulk boundary perturbation method

Although the number of IFs in AGNR heterojunctions can be determined using the winding number Z, its evaluation requires momentum–space calculations.²⁴ However, this approach is not directly suitable for finite AGNRH segments. To elucidate the correlation between the ESs of AGNR segments and the IFs of AGNRH segments, we analyze the energy levels of the four AGNRH structures shown in Fig. 3(a–d) as functions of the inter-AGNR coupling parameter t_es. In Fig. 4(a), eight in-gap energy levels appear at t_es = 0 because the 9–7–9 AGNRH is decoupled into isolated left-, central-, and right-AGNR segments. Each 9-AGNR segment with an R₁-type unit cell hosts N_L(R),A = 1 and N_L(R),B = 1 ESs. The 7-AGNR segment, which has an R₂-type unit cell, possesses N_C,A = 2 and N_C,B = 2 ESs. One of the N_C,A (N_C,B) states is an additional zero-energy mode with an extremely localized wave-function at sublattice-A(B), denoted as Ψ_7,A (Ψ_7,B). Because the 9-AGNR and 7-AGNR segments have finite lengths of 4 and 6 unit cells, respectively, these in-gap ESs hybridize and form bonding and antibonding energy levels. However, the coupling between Ψ_7,A and Ψ_7,B remains very weak.

Fig. 4 Energy levels of four AGNRH segments with different widths as functions of the inter-AGNR electron hopping strength t_es. (a) 9₄–7₆–9₄ AGNRH segment, (b) 15₄–13₆–15₄ AGNRH segment, (c) 21₄–19₆–21₄ AGNRH segment and (d) 27₄–25₆–27₄ AGNRH segment. The parameter t_es is illustrated in Fig. 1(b).

As t_es increases from 0 to the pristine hopping value t = 2.7 eV, four of the in-gap states exhibit linear t_es dependence, forming the curves labeled Σ_AB,c and Σ_AB,v. These levels merge into the bulk bands at t_es = 0.45t. The origin of Σ_AB,c and Σ_AB,v is the interaction between pairs of ESs with opposite chirality through weak boundary perturbation (see also Fig. 5). The curves labeled Σ_IF,c and Σ_IF,v correspond to the bonding and antibonding combinations of the robust ESs (Φ_7,A and Φ_7,B) in the 7-AGNR segment. These states are weakly dependent on t_es. Specifically, we find Σ_IF,c = 0.04552 eV (Σ_IF,v = −0.04552 eV) at t_es = 0, and Σ_IF,c = 0.02286 eV (Σ_IF,v = −0.02286 eV) at t_es = t_ppπ = 2.7 eV. The evolution of these levels with t_es demonstrates that the robust ESs of the 7-AGNR segment give rise to the interface states of the 9–7–9 AGNRH. The terminal states of the 9–7–9 AGNRH, labeled Σ_A,c and Σ_B,v near the CNP, are insensitive to t_es.

Fig. 5 (a–c) Probability densities corresponding to Σ_AB,c2 = 0.24355 eV, Σ_AB,c1 = 0.20466 eV and Σ_IF,c = 0.12651 eV in the 15₄–13₆–9₄ AGNRH segment with t_es = 0.1t. (d–f) Probability densities corresponding to Σ_AB,c2 = 0.47859 eV, Σ_AB,c1 = 0.32919 eV and Σ_IF,c = 0.12399 eV in the 15₄–13₆–15₄ AGNRH segment with t_es = 0.2t. (g–i) Probability densities of Σ_IF,c in the 15₄–13₆–9₄ AGNRH segment for various values of t_es: (g) Σ_IF,c = 0.10343 eV at t_es = 0.5 t, (h) Σ_IF,c = 0.07926 eV at t_es = 0.8t, and (i) Σ_IF,c = 0.06531 eV at t_es = 1t.

For the 15–13–15 AGNRH shown in Fig. 4(b), we obtain 14 in-gap levels at t_es = 0. The 15-AGNR segments possess N_L,A(B) = 2 and N_R,A(B) = 2, whereas the 13-AGNR segment has N_C,A(B) = 3. As t_es increases to t = 2.7 eV, eight ESs move out of the bulk gap. The levels Σ_IF,c and Σ_IF,v originate from the bonding and antibonding combinations of the robust ESs of the 13-AGNR segment (Φ_13,A and Φ_13,B). In this case, the robust ESs exhibit a significantly longer decay length than those in the 9–7–9 AGNRH. In Fig. 4(c) and (d), 12 ESs and 16 ESs leave the bulk gap in the 21–19–21 and 27–25–27 AGNRHs, respectively. The 21–19–21 AGNRH retains two IFs (forming Σ_IF,c and Σ_IF,v) and six ESs within the gap, while the 27–25–27 AGNRH contains two IFs and eight ESs. These results lead to the general rule for N-AGNR/(N − 2)-AGNR/N-AGNR heterostructures: N_IF,L,A = |N_C,A − N_L,B| and N_IF,R,B = |N_C,B − N_R,A|. This relation shows that the number and type of IFs in a symmetric N_out-AGNR/N_cen-AGNR/N_out-AGNR heterostructure with ΔN = (N_out − N_cen)/2 = 1 are determined by the ES-number difference and the chirality of the ESs in the AGNR segments with the maximal ES number. As demonstrated in Appendix A.1, this rule is also valid for ΔN ≠ 1, for example, ΔN = 3, 7, 9.

Next, we consider the case of Fig. 4(b) as an example to illustrate the charge (probability) density distribution of ESs and IFs at different t_es values, and to clarify why some ESs shift out of the bulk gap, leaving only two IFs in the N-AGNR/(N − 2)-AGNR/N-AGNR structure. As shown in Fig. 5(a–c), the AGNRH behaves as three weakly coupled segments at t_es = 0.1t. In Fig. 5(a), the probability density of Σ_AB,c2 is confined at the interface sites, reflecting interactions between sublattice-A and sublattice-B states. The probability density of Σ_AB,c1, shown in Fig. 5(b), also reflects opposite-chirality interactions but is primarily distributed over the outer 15-AGNR segments. By contrast, the probability density of Σ_IF,c is localized entirely within the central 13-AGNR segment (Fig. 5(c)).

At t_es = 0.2t, the distributions of Σ_AB,c2, Σ_AB,c1, and Σ_IF,c remain nearly unchanged. To clarify the effect of t_es on Σ_IF,c, we plot the probability densities for t_es = 0.5t, 0.8t, and 1.0t in Fig. 5(g–i). As t_es increases, the probability density of Σ_IF,c gradually extends into the outer 15-AGNR segments while conserving normalization. Based on Fig. 4(b) and 5, we conclude that the interface states of the N-AGNR/(N − 2)-AGNR/N-AGNR heterostructure originate from the ESs of the narrow (N − 2)-AGNR segment, which adopts the R₂-type unit cell. It is important to note that the IF states exhibit a direction-dependent decay length. Similar direction-dependent decay behavior has also been reported for the two-dimensional topological interface states in HgTe/CdTe heterostructures.⁵⁷

3.3. Interface states of AGNRH segments functioning as a single TDQD

The bulk gap of the 9–7–9 AGNRH segment is larger than that of the 15–13–15, 21–19–21, and 27–25–27 AGNRH segments, which is advantageous for suppressing thermal noise. However, the decay length of its IF wave functions is too short to allow the fabrication of gate electrodes. Increasing the width of AGNRH segments leads to longer decay lengths of IF wave functions, as reflected by the increased energy separation between Σ_IF,c and Σ_IF,v shown in Fig. 4. Because the bulk gaps of the 21–19–21 and 27–25–27 AGNRH segments are smaller than that of the 15–13–15 AGNRH segment, we focus on the transport properties of the interface states in the 15–13–15 AGNRH segment. The calculated transmission coefficients for this structure are shown in Fig. 6.

Fig. 6 Transmission coefficient of the 15_w–13_x–15_y AGNRH segment with zigzag-terminated edge structures coupled to electrodes. (a) of the 15₄–13₆–15₄ AGNRH segment for different values of t_es at Γ_L = Γ_R = Γ = 2.7 eV. (b) of the 15₄–13₆–15₄ AGNRH segment with a staggered sublattice potential δ, applied to the outer AGNR segments, for Γ = 2.7 eV and t_es = t. (c) of the 15₄–13₆–15₄ AGNRH segment with t_es = t and δ = 0 for different values of Γ, corresponding to various electrode materials. (d) of the 15_w–13₆–15_w AGNRH segment for different values of w, with t_es = t, Γ = 2.7 eV and δ = 0. A parameter of γ = 90 meV is used in panels (b) and (c).

Fig. 6(a) presents the transmission coefficient of a 15₄–13₆–15₄ AGNRH with zigzag edge terminations for different values of t_es. As t_es decreases, the spectrum splits into two resonant peaks, each reaching a maximum value of one. These two peaks clearly demonstrate that the 15–13–15 AGNRH segment hosts two nondegenerate resonant channels. The peak positions shift away from the CNP as t_es decreases. This behavior arises because reducing t_es effectively increases the barrier height experienced by the IFs, strengthening their confinement. Consequently, the peak widths also become narrower. In the context of t_es ≠ t, one may consider employing STM techniques to manipulate the electron-hopping strengths at the interface sites.¹¹

Fig. 6(b) shows for different staggered potentials δ_A = −δ_B = δ induced by the two-dimensional substrates supporting the outer 15-AGNR segments. Increasing δ enhances quantum confinement, and once the 15-AGNR segments behave as potential barriers, strong backward electron scattering is expected. However, since strain engineering can be used to minimize the staggered potential,^58–62 we neglect substrate-induced staggered potentials in the following discussion.

In Fig. 6(c), we show for the 15₄–13₆–15₄ AGNRH segment for various coupling strengths Γ, corresponding to different electrode materials. While the resonant peak positions remain unchanged, the peak widths vary according to the coupling between the electrodes and the edge atoms (i.e., variation of Γ_L = Γ_R = Γ). These results further highlight the topologically protected and symmetric localization characteristics of the IFs.

In Fig. 6(a)–(c), the length of the 15₄–13₆–15₄ segment is fixed. To examine size effects, we also calculate for 15_w–13₆–15_y AGNRHs with symmetrically varied outer segment lengths (w = y). Because the IFs are strongly localized, increasing w reduces the effective coupling strength between the electrodes and the IFs. As a result, the resonant peaks become narrower with increasing 15-AGNR segment length. Overall, the results in Fig. 6 demonstrate that the IFs in a 15–13–15 AGNRH function as a topological double quantum dot (TDQD). Each quantum dot hosts a single energy level that is well isolated from bulk states by an effective bulk gap on the order of one electron-volt (see Fig. 3(b)).

3.4. Nonlinear thermoelectric power of TDQD

The transmission coefficients shown in Fig. 6 were calculated within a single-particle framework. However, evaluating in the Coulomb blockade regime is extremely challenging.^63–68 When intra-site Coulomb interactions are included in the Hamiltonian of eqn (2), current theoretical approaches remain limited to mean-field treatments, i.e., one-particle approximations.³⁶ Since our interest lies in the low-energy modes near the CNP, we employ an extended Anderson model that incorporates effective intra-dot and inter-dot Coulomb interactions to investigate tunneling through the TDQD. The effective Hamiltonian is given by H_eff = H_SD + H_TDQD, where H_SD describes the source and drain electrodes and H_TDQD represents the TDQD:Here, E_j is the spin-independent energy level of dot j, and U_j = U_L(R) = U₀ and denote the intra-dot and inter-dot Coulomb interactions, respectively. The operator n_j,σ = c^†_j,σc_j,σ is the number operator, and t_x is the interdot hopping amplitude. For 15–13–15 AGNRH segments, both end states and interface states emerge within the middle gap and are well separated from the conduction and valence subbands by a band gap of approximately 1 eV [Fig. 3(b)]. This large separation suppresses optical-phonon assisted transport due to the phonon bottleneck effect;⁶⁹ therefore, thermal noise effects arising from bound-to-continuum states can be safely neglected in eqn (3). Noting that the phonon mean free path in graphene nanostructures can be reduced to 10 nm.^70–76 Since the lengths of the AGNRH segments considered in this work are smaller than this scale, electron-acoustic phonon scattering can be safely ignored.

The Green-function technique provides a powerful framework for calculating tunneling currents in strongly correlated nanostructures under nonequilibrium conditions.^77–80 Using the equation-of-motion method, the tunneling current from the left (right) electrode through the TDQD with electron–electron interactions is given by

The Fermi–Dirac distribution of electrode α is defined as f_α(ε) = 1/(exp[(ε − µ_α)/k_BT_α] + 1),with chemical potentials µ_L(R) = µ ± eV_a/2 under an applied bias of +V_a/2 and −V_a/2. For convenience, we set the Fermi energy µ = 0. The temperature bias is defined as ΔT = T_L − T_R. A closed-form expression for is given in ref. 80:

Here, ε_L = ε − E_L + iΓ_e,L and ε_R = ε − E_R + iΓ_e,R; Γ_e,L(R) is the tunneling rate determined by the dot-electrode coupling. The eight terms above correspond to the eight possible TDQD occupation configurations encountered by an incoming spin–σ electron. The configuration probabilities C_m are expressed as

C₁ = 1 − N_L,σ − N_R,σ − N_R,−σ + 〈n_R,σn_L,σ〉 + 〈n_R,−σn_L,σ〉 + 〈n_R,−σn_R,σ〉 − 〈n_R,−σn_R,σn_L,σ〉

C₂ = N_R,σ − 〈n_R,σn_L,σ〉 − 〈n_R,−σn_R,σ〉 + 〈n_R,−σn_R,σn_L,σ〉

C₃ = N_R,−σ − 〈n_R,−σn_L,σ〉 − 〈n_R,−σn_R,σ〉 + 〈n_R,−σn_R,σn_L,σ〉

C₄ = 〈n_R,−σn_R,σ〉 − 〈n_R,−σn_R,σn_L,σ〉

C₅ = N_L,σ − 〈n_R,σn_L,σ〉 − 〈n_R,−σn_L,σ〉 + 〈n_R,−σn_R,σn_L,σ〉

C₆ = 〈n_R,σn_L,σ〉 − 〈n_R,−σn_R,σn_L,σ〉

C₇ = 〈n_R,−σn_L,σ〉 − 〈n_R,−σn_R,σn_L,σ〉

C₈ = 〈n_R,−σn_R,σn_L,σ〉,

Here, N_L,σ is the single-particle occupation at site L. Two-particle and three-particle correlation functions are also included, and the full set is solved self-consistently so that , ensuring probability conservation. The reverse-bias current is obtained by exchanging the indices of . Note that the transmission formula is valid only for temperatures above the Kondo temperature.^81–83

Although many-body effects in voltage-driven transport have been widely explored, their impact on electrical power generation in nanostructures driven by a temperature bias-both in quasi-equilibrium and far-from-equilibrium regimes-remains insufficiently understood.^84–106 Here, we investigate these many-body effects for thermoelectric generators formed by the topological states (TSs) of AGNRHs. In realistic experimental conditions, the left (right) gate voltage not only tunes the energy level of the left (right) TQD but also affects the right (left) TQD; therefore, the left and right TQD energy levels are modulated as E_L = −(η₁eV_L,g + (1 − η₁)eV_R,g) + η_LeV_a and E_R = −((1 − η₂)eV_L,g + η₂V_R,g) + η_ReV_a, where η₁ and η₂ represent the effects of intra- and interdot Coulomb interactions¹⁰⁷ and η_L(R) (|η_L(R)| < 1/2) describe the orbital shifts induced by V_a (see Fig. (3)). We take U₀ = 90 meV and U₁ = 25 meV, consistent with their weak dependence on AGNR lengths.⁵⁶ In this study, we use η₁ = η₂ = 0.8 and a small orbital offset η_L = −η_R = 0.1. The tunneling rates Γ_e,L and Γ_e,R and the interdot hopping t_x are treated as tunable parameters controlled by the 13- and 15-AGNR segment lengths.

To analyze nonlinear temperature-driven power generation, we define the optimal electrical power as Ω_op = −J_op(V_op, ΔT) × V_op(ΔT), where J_op(V_op, ΔT) and V_op(ΔT) denote the optimized thermal current and thermal voltage, respectively.^93,103 The thermal voltage V_op(ΔT) is obtained by maximizing the electrical power Ω = −J_L(R)(V_a, ΔT) × V_a with respect to V_a.¹⁰³ This optimization is highly nontrivial in the Coulomb blockade regime under nonequilibrium conditions; thus, many theoretical studies avoid correlated transport and power optimization altogether.^91–104 The electrode temperatures are expressed as T_L = T₀ ± ΔT/2 and T_R = T₀ ∓ ΔT/2, where T₀ = (T_L + T_R)/2 is the average temperature.

3.4.1 Quasi-equilibrium and out-of-equilibrium scenarios. Fig. 7(a)–(c) show the charge stability diagram (N_t), electrical conductance (G_e), and Seebeck coefficient (S) as functions of the two gate voltages V_L,g and V_R,g at T₀ = 12 K. In Fig. 7(a), the charge stability diagram gives the total charge number . Nine distinct charge configurations (N_L, N_R) are observed, corresponding to TDQD occupancies of zero to four electrons. These results are consistent with experimental charge stability diagrams reported for other serial DQD systems at low temperature.^107–111 In Fig. 7(b), the electrical conductance exhibits eight peaks aligned with the charge-transition lines in Fig. 7(a). These peaks indicate that the TDQD provides eight resonant transport channels under weak coupling (Γ_e,t = t_x = 1 meV). The first peak (P₁) corresponds to the lowest-energy TDQD state, while the last peak (P₈) corresponds to the highest-energy state.⁸⁰ The Seebeck coefficient in Fig. 7(c) exhibits a bipolar behavior, with negative (positive) values associated with electron (hole-like) transport. The S spectra thus encode both the thermal transport characteristics and the charge stability structure of the TDQD. The quantities G_e and S are expressed in units of G₀ = 2e²/h = 77.5 µS and k_B/e = 86.25 µV K⁻¹, respectively. These results represent the quasi-equilibrium regime of TDQD thermoelectric transport.

Fig. 7 (a) Charge stability diagram (N_t), (b) electrical conductance (G_e), and (c) Seebeck coefficient (S) as functions of gate voltages V_L,g and V_R,g at T₀ = 12 K with ΔT = 0. (d) Total occupation number N_t, (e) optimized thermal current J_op(ΔT), and (f) optimized thermal voltage V_op(ΔT) at T₀ = 24 K with ΔT = 12 K. A weak interdot coupling of t_x = Γ_e,t = 1 meV is used.

To examine the out-of-equilibrium regime, Fig. 7(d)–(f) present the total charge number N_t, optimized thermal current J_op, and optimized thermal voltage V_op at average temperature T₀ = 24 K with a thermal bias ΔT = 12 K. In Fig. 7(d), the charge-transition region between (0, 1) and (1, 1) is broader than that between (1, 0) and (1, 1), reflecting the asymmetric electrode temperatures (T_L = 30 K, T_R = 18 K). Unlike G_e, the optimized thermal current J_op [Fig. 7(e)] displays a bipolar dependence on the symmetric gate voltage V_L,g = V_R,g = V_g, with vanishingly small values at the conductance peaks. This indicates that maximum J_op does not occur when TDQD energy levels align with the Fermi energy. The optimized thermal voltage V_op [Fig. 7(f)] also shows bipolar behavior but with broader features than the Seebeck coefficient. Although J_op and V_op exhibit opposite Coulomb oscillations with respect to V_g, the optimized output power Ω_op remains positive.

3.4.2 Comparison between noninteracting and interacting cases. Many previous studies predicted the power output of nanoscale thermoelectric generators without considering electron–electron Coulomb interactions.^88–105 To clarify the role of Coulomb interactions and temperature bias on V_op and J_op at room temperature, Fig. 8 shows V_op, J_op, and Ω_op as functions of gate voltage V_g and temperature bias ΔT at T₀ = 288 K. V_g is varied from −100 mV to 0, and ΔT from 0 to 240 K. Fig. 8(a)–(c) correspond to the noninteracting case, and Fig. 8(d)–(f) to the interacting case. As seen in Fig. 8(a) and (d), V_op differs only slightly, mainly near V_g ≥ −10 mV: without interactions, V_op(V_g = 0) = 0 for all ΔT due to electron–hole symmetry, whereas interactions lift this symmetry, giving small but finite V_op at V_g = 0. In contrast, J_op shows large differences: the noninteracting maximum J_op,max = 0.566 at V_g = −25 mV and ΔT = 240 K [Fig. 8(b)] is significantly overestimated compared with the interacting value J_op,max = 0.418 at V_g = −32 mV [Fig. 8(e)]. In the (0, 0) configuration, the Coulomb blockade suppresses J_op,max but has little effect on V_op,max. The behavior of V_op resembles that of the Seebeck coefficient S, which is largely insensitive to electron Coulomb interactions.⁸⁶

Fig. 8 (a–c) Noninteracting case: optimized thermal voltage V_op, thermal current J_op, and electrical power Ω_op as functions of symmetric gate voltage V_L,g = V_R,g = V_g and temperature bias ΔT at T₀ = 288 K. (d–f) Interacting case under the same conditions. Other physical parameters are identical to those in Fig. 7. Units are J₀ = 0.773 nA and Ω₀ = 77.3 pW.

To further clarify the relationship between Ω_op and ΔT shown in Fig. 8(f) and 9(a)–(c) present S_op = V_op/ΔT, J_op, and Ω_op as functions of ΔT for V_g = −32, −46, and −60 mV. For V_g = −60 mV, Ω_op reaches its maximum value at k_BΔT = 20 meV. When the temperature bias (ΔT) approaches zero, the Seebeck coefficient (S_op = V_g/T₀) (for V_g = 60 mV) depends only on the separation between the TDQD energy levels and the Fermi energy (µ = 0), as well as the average temperature (T₀). This indicates that the Seebeck coefficient of the TDQD in the weak-coupling limit (t_x = Γ_et → 0) can function as an ultra-sensitive thermal detector.¹¹² As shown in Fig. 9(a), S_op is only weakly dependent on ΔT for V_g = −32 mV, which also explains why V_op exhibits a linear dependence on ΔT for V_g > – 24 mV in Fig. 8(d). When the TDQD energy levels lie far from the Fermi energy (V_g = −60 mV), S_op remains linear in ΔT, and this linearity persists up to k_BΔT = 20 meV. It is noteworthy that is not a perturbative parameter, yet the nonlinear correction to V_op is restricted to a quadratic term, with no observable cubic or higher-order contributions which is an unexpected behavior. This finding is consistent with the recent analytical results reported in ref. 105, where a single noninteracting quantum dot with one energy level was considered. In Fig. 9(b), we find that J_op remains a linear function of ΔT even at V_g = −60 mV. Based on the trends identified in Fig. 9(a), Ω_op = −J_opV_op for V_g = −46 and −60 mV, with Ω_op,max deviating from the expected (ΔT)² scaling due to nonlinearities in V_op.

Fig. 9 (a) Optimized thermal voltage S_op = V_op/ΔT, (b) thermal current J_op(ΔT), and (c) electrical power output Ω_op(ΔT) = −J_op(ΔT) × V_op(ΔT) as functions of ΔT for V_g = −32, −36, and −60 mV at T₀ = 288 K. Other physical parameters are the same as in Fig. 8. The units are J₀ = 0.773 nA and Ω₀ = 77.3 pW.

3.4.3 Electrical power rectification. Finally, we examine the nonlinear thermoelectric properties of the TDQD under asymmetric conditions, V_L,g ≠ V_R,g and Γ_e,L ≠ Γ_e,R. Fig. 10(a)–(d) show 2D plots of total occupation number N_t, optimized thermal voltage V_op(ΔT), thermal current J_op(ΔT), and electrical power Ω_op(ΔT) as functions of temperature bias ΔT and interdot hopping t_x, for parameters Γ_e,L = t_x, Γ_e,R = t_x/10, E_L = 2 meV, E_R = 8 meV, and T₀ = 48 K. In Fig. 10(a), N_t decreases with forward temperature bias (ΔT > 0) and increases with backward bias (ΔT < 0) at weak hopping (t_x = 1 meV). N_t is primarily determined by N_L, as the left dot energy level is close to the Fermi energy (µ = 0). At stronger coupling (t_x = 10 meV), the broadening of E_L due to Γ_e,L = t_x makes N_L nearly independent of bias direction, explaining the reduced sensitivity of N_t to ΔT sign.

Fig. 10 (a) Total occupation number N_t, (b) optimized thermal voltage |V_op(ΔT, T₀)|, (c) optimized thermal current |J_op(ΔT, T₀)|, and (d) optimized power output Ω_op(ΔT, T₀) = −J_op(ΔT, T₀) × V_op(ΔT, T₀) as functions of temperature bias ΔT and interdot hopping t_x at k_BT₀ = 4 meV, E_L = 2 meV, E_R = 8 meV, Γ_e,L = t_x, and Γ_e,R = 0.1t_x. Other physical parameters are the same as in Fig. 8. Units are J₀ = 0.773 nA and Ω₀ = 77.3 pW.

The optimized thermal voltage V_op is less sensitive to the direction of ΔT (Fig. 10(b)). In contrast, J_op and Ω_op exhibit pronounced rectification (Fig. 10(c) and (d)), providing strong evidence that electron Coulomb interactions govern power rectification; asymmetrical structures alone cannot induce this effect in noninteracting nanostructures. The asymmetry of J_op and Ω_op can be understood from the first term of the transmission coefficient:

andFor ΔT > 0, the relevant probability P₁ ≈ 1 − N_L,σ − N_R,σ − N_R,−σ, while for ΔT < 0, P₁ ≈ 1 − N_R,σ − N_L,σ − N_L,−σ. Since N_L(ΔT < 0) ≫ N_L(ΔT > 0), this accounts for the observed asymmetry in J_op and Ω_op. Rectification of J_op(ΔT) and Ω_op(ΔT) requires TDQDs with strong electron Coulomb interactions, asymmetric gate voltages, and asymmetric 15-AGNR segment lengths.

4. Conclusion

In this study, we systematically clarified the number and types of interface states (IFs) that emerge in N-AGNR/(N − 2)-AGNR/N-AGNR heterostructure (AGNRH) segments lacking translational symmetry, based on their real-space geometric structure. We established a direct correspondence between the end states (ESs) of individual AGNR segments and the topological IFs of the resulting AGNRHs. For AGNR segments with R₁-type unit cells, the numbers of ESs were determined under a uniform longitudinal electric field for N = 13, 15, 19, 21, 25, and 27. These ESs satisfy the relations N = N_A(B) × 6 + 1 and N = N_A(B) × 6 + 3, where N_A(B) denotes the number of ESs with A-(B-)chirality. For AGNR segments with R₂-type unit cells, the total number of ESs increases to N_A(B) + 1.

The Stark effect induced by the applied electric field lifts the degeneracy of ESs in AGNRH segments and enables clear spectral distinction between ESs and IFs. The real-space bulk boundary perturbation approach further shows that states with opposite chirality can hybridize through junction-site perturbations and shift out of the bulk gap. We demonstrated that the number and chirality of IFs, N_IF,β, in symmetric AGNRHs composed of N = 3p and N = 3p + 1 AGNR segments are fully determined by the ESs belonging to the outer and central AGNR components. Denoting by N_O and N_C the numbers of ESs of the outer and central AGNR segments, respectively, the IF number satisfies N_IF,β = |N_O,B(A) − N_C,A(B)|, where β specifies the chirality of the IFs. When N_O,B > N_C,A, the IFs exhibit B-chirality, while they acquire A-chirality when N_C,A > N_O,B. For instance, the IFs of a 15–13–15 AGNRH segment originate from the ESs of the central 13-AGNR segment with R₂-type unit cells, whereas the IFs of a 27–13–27 AGNRH segment arise from the ESs of the outer 27-AGNR segments.

Using the calculated transmission coefficient for N_out-AGNR/N_cen-AGNR/N_out-AGNR segments with zigzag edge structures coupled to electrodes, we demonstrated that IFs act as a topological double quantum dot (TDQD) when the IFs are formed by the ESs of the central AGNR segment. By employing an Anderson model incorporating effective intradot and interdot Coulomb interactions, we derived an analytical expression for the tunneling current through the TDQD via the transmission coefficient . The thermoelectric performance of AGNRH-based TDQDs was further analyzed in the context of their potential application as graphene-nanoribbon power generators.¹¹³

The nonlinear thermoelectric power output of TDQDs exhibits several notable features. (a) The optimized power output is favored in either the electron-dilute (0, 0) or hole-dilute (2, 2) charge configurations. In the (0, 0) configuration, Coulomb blockade strongly suppresses J_op,max while leaving V_op,max largely unaffected. (b) Nonlinear temperature bias induces only a quadratic correction to the thermal voltage V_op, with no cubic or higher-order terms observed. (c) Even in the presence of strong electron Coulomb interactions, the thermal power output remains highly enhanced under nonlinear temperature bias. (d) Direction-dependent power output arises from strong electron correlations and the asymmetry in outer AGNR lengths. Owing to their well-isolated energy levels located deep within the bulk gap of 15–13–15 AGNRHs, TDQDs represent promising candidates for high-temperature thermoelectric power generation.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that supports the finding of this study are available within the article.

A. Appendices

A.1. Interface states of wide-AGNR/narrow-AGNR/wide-AGNR AGNRH segments

For N-AGNR/(N − 2)-AGNR/N-AGNR AGNRH segments, the width difference ΔN = (N_out − N_cen)/2 equals ΔN = 1, where N_out and N_cen denote the atom widths of the outer and central AGNR segments, respectively. In Fig. 4, our analysis focused on the case ΔN = 1. In this appendix, we examine how increasing ΔN affects both the number and chirality of interface states (IFs) in AGNRHs.

Fig. 11(a–d) display the energy spectra of four AGNRH segments as functions of t_es: (a) 27₄–25₁₀–27₄, (b) 27₄–21₁₀–27₄, (c) 27₄–13₁₀–27₄, and (d) 27₄–9₁₀–27₄. These correspond to ΔN = 1, 3, 7, and 9, respectively, with a fixed outer width of N_out = 27. In all cases, the outer AGNRs have an R₁-type unit cell, whereas the central AGNR adopts an R₂-type unit cell.

Fig. 11 Energy levels of four AGNRH segments with varying central widths as functions of the inter-AGNR electron hopping strength t_es: (a) 27₄–25₁₀–27₄, (b) 27₄–21₁₀–27₄, (c) 27₄–13₁₀–27₄, and (d) 27₄–9₁₀–27₄ AGNRH segments.

In Fig. 11(a), one interface state of sublattice-A (sublattice-B) arises at the 27-AGNR/25-AGNR (25-AGNR/27-AGNR) junction. The resulting pair of IF levels, Σ_IF,c and Σ_IF,v, originates from their mutual coupling. Compared with the 27₄–25₆–27₄ case in Fig. 4(d), the energy separation between Σ_IF,c and Σ_IF,v in Fig. 11(a) is smaller due to the longer 25₁₀ central segment.

In Fig. 11(b), eight in-gap states appear, all attributable to the end states of the 27–21–27 AGNRH system. Since both 27-AGNR and 21-AGNR segments have N_A(B) = 4, the number of IFs is zero, consistent with the rule that IFs are determined by the difference in ES counts.

In Fig. 11(c), a single interface state with sublattice-B(A) character is present at each 27-AGNR/13-AGNR (13-AGNR/27-AGNR) junction because the 13-AGNR segment has N_A(B) = 3. These IFs originate from the ESs of the outer AGNRs; consequently, their mutual overlap is extremely small, whereas the IFs strongly couple to the terminal states of the 27–13–27 AGNRH segment. Their evolution with respect to t_es is therefore markedly different from IFs originating from the central AGNR segment.

In Fig. 11(d), two IFs appear at the 27-AGNR/9-AGNR (9-AGNR/27-AGNR) junctions because the 9-AGNR segment possesses N_A(B) = 2.

A particularly large energy separation between Σ_IF,c2 and Σ_IF,v2 is seen in Fig. 11(d). To understand its origin, Fig. 12 presents the energy levels of 27_w–9₆–27_w AGNRHs for four different outer lengths w = 6, 7, 8, and 9. As w increases, the splitting between Σ_IF,c2 and Σ_IF,v2 decreases significantly. This trend indicates that these two levels arise from IFs coupled to end states of the 27-AGNR segments that have long decay lengths. In contrast, the splitting between Σ_IF,c1 and Σ_IF,v1 shows only weak dependence on w, implying that these states originate from IFs associated with short-decay-length end states of the 27-AGNR segments.

Fig. 12 Energy levels of 27_w–9₆–27_w AGNRH segments as functions of the inter-AGNR hopping strength t_es for different outer 27-AGNR lengths: (a) 27₆–9₆–27₆, (b) 27₇–9₆–27₇, (c) 27₈–9₆–27₈, and (d) 27₉–9₆–27₉ AGNRH segments.

Fig. 13 shows the probability densities of representative IF levels specifically, the pair Σ_IF,c2,a = 102.2335 meV and Σ_IF,c2,b = 102.2232 meV, the level Σ_IF,c1 = 15.4212 meV, and the terminal-state level Σ_A,c = 0.539 meV corresponding to Fig. 12(a). In Fig. 13(a) and (b), the spatial distributions of Σ_IF,c2,a and Σ_IF,c2,b are nearly indistinguishable; furthermore, their amplitudes inside the central 9-AGNR segment are extremely weak, indicating that transport through Σ_IF,c2 is strongly suppressed.

Fig. 13 Probability densities of selected energy levels in the 27₆–9₆–27₆ AGNRH segment with t_es = t: (a) Σ_IF,c2,a = 102.2335 meV, (b) Σ_IF,c2,b = 102.2232 meV, (c) Σ_IF,c1 = 15.4212 meV, and (d) Σ_A,c = 0.539 meV.

In contrast, Σ_IF,c1 exhibits a pronounced probability density within the central 9-AGNR segment, making it a more favorable transport channel. As shown in Fig. 13(c), Σ_IF,c1 effectively behaves as a set of four topological dots (TDs). Fig. 13(d) demonstrates that the state Σ_A,c is primarily associated with the terminal states of the 27–9–27 AGNRH structure.

Acknowledgements

This work was supported by the National Science and Technology Council, Taiwan under Contract No. MOST 107-2112-M-008-023MY2.

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Electric Field Effect in Atomically Thin Carbon Films, Science, 2004, 306, 666.
2. J. Cai, P. Ruffieux, R. Jaafar, M. Bieri, T. Braun, S. Blankenburg, M. Muoth, A. P. Seitsonen, M. Saleh, X. Feng, K. Mullen and R. Fasel, Atomically precise bottom-up fabrication of graphene nanoribbons, Nature, 2010, 466, 470.
3. Y. C. Chen, T. Cao, C. Chen, Z. Pedramraz, D. Haberer, D. G. de Oteyza, F. R. Fischer, S. G. Louie and M. F. Crommie, Molecular bandgap engineering of bottom-up synthesized graphene nanoribbon heterojunctions, Nat. Nanotechnol., 2015, 10, 156.
4. S. Y. Wang, L. Talirz, C. A. Pignedoli, X. L. Feng, K. Mullen, R. Fasel and P. Ruffieux, Giant edge state splitting at atomically precise graphene zigzag edges, Nat. Commun., 2015, 7, 11507.
5. J. P. Llinas, A. Fairbrother, B. G. Borin, W. Shi, K. Lee, S. Wu, B. Y. Choi, R. Braganza, J. Lear and N. Kau, Short-channel field-effect transistors with 9-atom and 13-atom wide graphene nanoribbons, Nat. Commun., 2017, 8, 633.
6. N. Merino-Diez, A. Garcia-Lekue, E. Carbonell-Sanroma, J. C. Li, M. Corso, L. Colazzo, F. Sedona, D. Sanchez-Portal, J. I. Pascual and D. G. de Oteyza, Width-Dependent Band Gap in Armchair Graphene Nanoribbons Reveals Fermi Level Pinning on Au(111), ACS Nano, 2017, 11, 11661.
7. O. Groning, S. Wang, X. Yao, C. A. Pignedoli, G. B. Barin, C. Daniels, A. Cupo, V. Meunier, X. Feng and A. Narita, et al., Engineering of robust topological quantum phases in graphene nanoribbons, Nature, 2018, 560, 209.
8. D. J. Rizzo, G. Veber, T. Cao, C. Bronner, T. Chen, F. Zhao, H. Rodriguez, S. G. Louie, M. F. Crommie and F. R. Fischer, Topological band engineering of graphene nanoribbons, Nature, 2018, 560, 204.
9. L. H. Yan and P. Liljeroth, Engineered electronic states in atomically precise artificial lattices and graphene nanoribbons, Adv. Phys.: X, 2019, 4, 1651672.
10. D. J. Rizzo, G. Veber, J. W. Jiang, R. McCurdy, T. C. C. Bronner, T. Chen, S. G. Louie, F. R. Fischer and M. F. Crommie, Inducing metallicity in graphene nanoribbons via zero-mode superlattices, Science, 2020, 369, 1597.
11. D. J. Rizzo, J. W. Jiang, D. Joshi, G. Veber, C. Bronner, R. A. Durr, P. H. Jacobse, T. Cao, A. Kalayjian, H. Rodriguez, P. Butler, T. Chen, S. G. Louie, F. R. Fischer and M. F. Crommie, Rationally designed topological quantum dots in bottom-up graphene nanoribbons, ACS Nano, 2021, 15, 20633.
12. Q. Sun, Y. Yan, X. L. Yao, K. Mullen, A. Narita, R. Fasel and P. Ruffieux, Evolution of the topological energy band in graphene Nanoribbons, J. Phys. Chem. Lett., 2021, 12, 8679.
13. S. T Song, Y. Teng, W. C. Tang, Z. Xu, Y. Y. He, J. W. Ruan, T. Kojima, W. P. Hu, F. J. Giessibl, H. Sakaguchi, S. G. Louie and J. Lu, Janus graphene nanoribbons with localized states on a single zigzag edge, Nature, 2025, 637, 580.
14. E. Leobandung, L. J. Guo, Y. Wang and S. Y. Chou, Observation of quantum effects and Coulomb blockade in silicon quantum-dot transistors at temperatures over 100 K, Appl. Phys. Lett., 1995, 67, 938.
15. G. D. Mahan and J. O. Sofo, The best thermoelectric, Proc. Natl. Acad. Sci. U. S. A., 1996, 93, 7436.
16. D. Loss and D. P. DiVincenzo, Quantum computation with quantum dots, Phys. Rev. A, 1998, 57, 120.
17. K. Ono, D. G. Austing, Y. Tokura and S. Tarucha, Current Rectification by Pauli Exclusion in a Weakly Coupled Double Quantum Dot System, Science, 2002, 297, 1313.
18. D. P. DiVincenzo, Double quantum dot as a quantum bit, Science, 2005, 309, 2173.
19. H. M. Wang, H. S. Wang, C. X. Ma, L. X. Chen, C. X. Jiang, C. Chen, X. M. Xie, A. P. Li and X. R. Wang, Graphene nanoribbons for quantum electronics, Nat. Rev. Phys., 2021, 3, 791.
20. J. Zhang, L. Qian, G. B. Barin, P. P. Chen, K. Mullen, P. Ruffieux, R. Fasel, J. Zhang, M. Calame and M. L Perrin, Double quantum dots in atomically-precise graphene nanoribbons, Mater. Quantum Technol., 2023, 3, 036201.
21. T. Cao, F. Z. Zhao and S. G. Louie, Topological Phases in Graphene Nanoribbons: Junction States, Spin Centers and Quantum Spin Chains, Phys. Rev. Lett., 2017, 119, 076401.
22. K. S. Lin and M. Y. Chou, Topological Properties of Gapped Graphene Nanoribbons with Spatial Symmetries, Nano Lett., 2018, 18, 7254.
23. J. W. Rhim, J. H. Bardarson and R. J. Slager, Unified bulk-boundary correspondence for band insulators, Phys. Rev. B, 2018, 97, 115143.
24. J. W. Jiang and S. G. Louie, Topology Classification using Chiral Symmetry and Spin Correlations in Graphene Nanoribbons, Nano Lett., 2021, 21, 197.
25. C. H. Lu and E. Y. T Li, A new graph theory to unravel the bulk-boundary correspondence of graphene nanoribbons, Carbon, 2024, 230, 119624.
26. Y. W. Gu, Z. Qiu and K. Mullen, Nanographenes and Graphene Nanoribbons as Multitalents of Present and Future Materials Science, J. Am. Chem. Soc., 2022, 144, 11499.
27. L. Maalysheva and A. Onipko, Spectra of π Electrons in Graphene as a Macromolecule, Phys. Rev. Lett., 2008, 100, 186806.
28. A. Onipko, Spectra of π electrons in graphene as an aletrnant macromoleculae and its specific features in quantum conductance, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 245412.
29. A. V. Nikolaev, A. V. Bibikov, A. V. Avdeenkov, I. V. Bodrenko and E. V. Tkalya, Electronic and transport properties of rectangular graphene macromolecules and zigzag carbon nanotubes of finite length, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 045418.
30. A. V. Rozhkov and F. Nori, Exact wave functions for an electron on a graphene triangular quantum dot, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 155401.
31. P. Potasz, A. D. Guclu and P. Hawrylak, Zero-energy states in triangular and trapezoidal graphene structures, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 033403.
32. A. Talkachov and E. Babaev, Wave functions and edge states in rectangular honeycomb lattices revisited: Nanoflakes, armchair and zigzag nanoribbons, and nanotubes, Phys. Rev. B, 2023, 107, 045419.
33. M. Pilar Lopez-Sancho and M. Carmen Munoz, Topologically protected edge and confined states in finite armchair graphene nanoribbons and their junctions, Phys. Rev. B, 2021, 104, 245402.
34. M. Ohfuchi and S. Sato, Energetics and magnetism of topological graphene nanoribbons, J. Appl. Phys., 2021, 129, 064305.
35. F. Z. Zhao, T. Cao and S. G. Louie, Topological Phases in Graphene Nanoribbons Tuned by Electric Fields, Phys. Rev. Lett., 2021, 127, 166401.
36. M. J. J. Mangnus, F. R. Fischer, M. F. Crommie, I. Swart and P. H. Jacobse, Charge transport in topological graphene nanoribbons and nanoribbon heterostructures, Phys. Rev. B, 2022, 105, 115424.
37. S. Traverso, M. Sassetti and N. T. Ziani, Role of the edges in a quasicrystalline Haldane model, Phys. Rev. B, 2022, 106, 125428.
38. N. V. Tepliakov, J. Lischner, E. Kaxiras, A. A. Mostofi and M. Pizzochero, Unveiling and Manipulating Hidden Symmetries in Graphene Nanoribbons, Phys. Rev. Lett., 2023, 130, 026401.
39. A. H. Huang, S. S. Ke, J. H. Guan, J. Li and W. K. Lou, Strain-induced topological phase transition in graphene nanoribbons, Phys. Rev. B, 2024, 109, 045408.
40. D. M. T. Kuo, Topological states in finite graphene nanoribbons tuned by electric fields, J. Phys.: Condens. Matter, 2025, 37, 085304.
41. J. Ostmeyer, L. Razmadze, E. Berkowitz, T. Luu and U. G. Meibner, Effective theory for graphene nanoribbons with junctions, Phys. Rev. B, 2024, 109, 195135.
42. H. Abdelsalam, D. Corona, R. B. Payod, M. A. S. Sakr, O. H. Abd-Elkader, Q. F Zhang and V. A. Saroka, Topological Junction States in Graphene Nanoribbons: A Route to Topological Chemistry, Nano Lett., 2025, 25, 10594.
43. M. Guzman, D. Bartolo and D. Carpentier, Geometry and Topology Tango in Ordered and Amorphous Chiral Matter, SciPost Phys., 2022, 12, 038.
44. K. Nakada, M. Fujita, G. Dresselhaus and M. S. Dresselhaus, Edge state in graphene ribbons: Nanometer size effect and edge shape dependence, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 17954.
45. K. Wakabayashi, M. Fujita, H. Ajiki and M. Sigrist, Electronic and magnetic properties of nanographite ribbons, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 8271.
46. K. Wakabayashi, K. Sasaki, T. Nakanishi and T. Enoki, Electronic states of graphene nanoribbons and analytical solutions, Sci. Technol. Adv. Mater., 2010, 11, 054504.
47. M. T. Allen, J. Martin and A. Yacoby, Gate-defined quantum confinement in suspended bilayer graphene, Nat. Commun., 2012, 3, 934.
48. S. Das Sarma, S. Adam, E. H. Hwang and E. Rossi, Electronic transport in two-dimensional graphene, Rev. Mod. Phys., 2011, 83, 407.
49. V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea and A. H. Castro Neto, Electron-Electron Interactions in Graphene: Current Status and Perspectives, Rev. Mod. Phys., 2012, 84, 1067.
50. C. Volk, C. Neumann, S. Kazarski, S. Fringes, S. Engels, F. Haupt, A. Mueller and C. Stampfer, Probing relaxation times in graphene quantum dots, Nat. Commun., 2013, 4, 1753.
51. M. Brotons-Gisbert, B. Artur, S. Kumar, R. Picard, R. Proux, M. Gray, K. S. Burch, K. Watanabe, T. Taniguchi and B. D. Gerardot, Coulomb blockade in an atomically thin quantum dot coupled to a tunable Fermi reservoir, Nat. Nanotechnol., 2019, 14, 442.
52. Q. F. Sun and X. C. Xie, CT-Invariant Quantum Spin Hall Effect in Ferromagnetic Graphene, Phys. Rev. Lett., 2010, 104, 066805.
53. D. M. T. Kuo and Y. C. Chang, Contact Effects on Thermoelectric Properties of Textured Graphene Nanoribbons, Nanomaterials, 2022, 12, 3357.
54. D. M. T. Kuo, Thermal rectification through the topological states of asymmetrical length armchair graphene nanoribbons heterostructures with vacancies, Nanotechnology, 2023, 34, 505401.
55. D. M. T. Kuo, Charge transport through the multiple end zigzag edge states of armchair graphene nanoribbons and heterojunctions, RSC Adv., 2024, 14, 20113.
56. D. M. T. Kuo, Room-Temperature Pauli Spin Blockade and Current Rectification in 15-13-15 Armchair Graphene Nanoribbon Heterostructures, Nanoscale, 2025, 17, 18920.
57. Y. C. Chang, J. N. Schulman, G. Bastard and Y. Guldner, Effects of quasi-interface states in HgTe-CdTe superlattices, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 2557.
58. G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly and J. van den Brink, Substrate-induced band gap in graphene on hexagonal boron nitride:: Ab initio density functional calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 073103.
59. S. Y. Zhou, G. H. Gweon, A. V. Fedorov, P. N. First, W. A. De Heer, D.-H. Lee, F. Guinea, A. H. Castro Neto and A. Lanzara, Towards wafer-size graphene layers by atmospheric pressure graphitization of silicon carbide, Nat. Mater., 2009, 6, 770.
60. Z. H. Qiao, S. Y. A. Yang, B. Wang, Y. Q. Yao and Q. Niu, Spin-polarized and valley helical edge modes in graphene nanoribbons, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 035431.
61. C. P. Lu, G. H. Li, K. Watanabe, T. Taniguchi and E. Y. Andreil, MoS2: Choice Substrate for Accessing and Tuning the Electronic Properties of Graphene, Phys. Rev. Lett., 2014, 113, 156804.
62. J. Zhao, P. Ji, Y. Q. Li, R. Li, K. M. Zhang, H. Tian, K. C. Yu, B. Bian, L. H. Hao, X. Xiao, W. Griffin, N. Dudeck, R. Morol, L. Mal and W. A. de Heerl, Ultrahigh mobility semiconducting epitaxial graphene on silicon carbide, Nature, 2024, 625, 60.
63. F. Sols, F. Guinea and A. H. Castro Neto, Coulomb Blockade in Graphene Nanoribbons, Phys. Rev. Lett., 2007, 99, 166803.
64. K. Todd, H. T. Chou, S. Amasha and D. Goldhaber-Gordon, Quantum Dot behavior in Graphene Nanoconstrictions, Nano Lett., 2009, 9, 416.
65. K. Cernevics, O. V. Yazyev and M. Pizzochero, Electronic transport across quantum dots in graphene nanoribbons: Toward built-in gap-tunable metal-semiconductor-metal heterojunctions, Phys. Rev. B, 2020, 102, 201406(R).
66. M. E. Abbassi, M. L. Perrin, G. B. Barin, S. Sangtarash, J. Overbeck, O. Braun, C. J. Lambert, Q. Sun, T. Prechtl, A. Narita, K. Mullen, P. Ruffieux, H. Sadeghi, R. Fasel and M. Calame, Controlled Quantum Dot Formation in Atomically Engineered Graphene Nanoribbon Field-Effect Transistors, ACS Nano, 2020, 14, 5754.
67. W. H. Huang, O. Braun, D. I. Indolese, G. B. Barin, G. Gandus, M. Stiefel, A. Olziersky, K. Mullen, M. Luisier, D. Passerone, P. Ruffieux, C. Schonenberger, K. Watanabe, T. Taniguchi, R. Fasel, J. Zhang, M. Calame and M. L. Perrin, Edge Contacts to Atomically Precise Graphene Nanoribbons, ACS Nano, 2023, 17, 18706.
68. J. Zhang, L. Qian, G. B. Barin, A. H. S. Daaoub, P. P. Chen, K. Mullen, S. Sangtarash, P. Ruffieux, R. Fasel, H. Sadeghi, J. Zhang, M. Calame and M. L. Perrin, Contacting individual graphene nanoribbons using carbon nanotube electrodes, Nat. Electron., 2023, 6, 572.
69. J. W. Fang, P. Li, J. Y. Zhang, C. C. Qin, T. Debnath, W. Huang and R. F. Chen, Stimulating Phonon Bottleneck Effect in Organic Semiconductors by Charge-Transfer-Mediated J-Aggregation, J. Am. Chem. Soc., 2024, 146, 961.
70. Y. M. Zuev, W. Chang and P. Kim, Thermoelectric and magneto thermoelectric transport measurements of graphene, Phys. Rev. Lett., 2009, 102, 096807.
71. P. Wei, W. Z. Bao, Y. Pu, C. N. Lau and J. Shi, Anomalous thermoelectric transport of Dirac particles in graphene, Phys. Rev. Lett., 2009, 102, 166808.
72. Y. Xu, Z. Y. Li and W. H. Duan, Thermal and thermoelectric properties of graphene, Small, 2014, 10, 2182.
73. Y. H. Wang, J. Yang, L. Y. Wang, K. Du, Q. Yin and Q. J. Yin, Polypyrrole/graphene/polyaniline ternary nanocomposite with high thermoelectric power factor, ACS Appl. Mater. Interfaces, 2017, 9, 20124.
74. M. Li, D. L. Cortie, J. X. Liu, D. H. Yu, S. M. K. N. Islam, L. L. Zhao, D. R. G. Mitchell, R. A. Mole, M. B. Cortie and S. X. Dou, Ultra-high thermoelectric performance in graphene incorporated Cu2Se: Role of mismatching phonon modes, Nano Energy, 2018, 53, 993.
75. S. Ghosh, S. Harish, M. Ohtaki and B. B. Saha, Thermoelectric figure of merit enhancement in cement composites with graphene and transition metal oxides, Mater. Today Energy, 2020, 18, 100492.
76. Y. Y. Wang, D. R. Chen, J. K. Wu, T. H. Wang, C. S. Chuang, S. Y. Huang, W. P. Hsieh, M. Hofmann, Y. H. Chang and Y. P. Hsieh, Two-dimensional mechano-thermoelectric heterojunctions for self-powered Strain Sensors, Nano Lett., 2021, 21, 6990.
77. D. M. T. Kuo, S. Y. Shiau and Y. C. Chang, Theory of spin blockade, charge ratchet effect, and thermoelectrical behavior in serially coupled quantum dot system, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 245303.
78. C. C. Chen, Y. C. Chang and D. M. T. Kuo, Quantum interference and electron correlation in charge transport through triangular quantum dot molecules, Phys. Chem. Chem. Phys., 2015, 17, 6606.
79. C. C. Chen, D. M. T. Kuo and Y. C. Chang, Quantum interference and structure-dependent orbital-filling effects on the thermoelectric properties of quantum dot molecules, Phys. Chem. Chem. Phys., 2015, 17, 19386.
80. D. M. T. Kuo, Temperature-stable tunneling current in serial double quantum dots: insights from nonequilibrium Green's functions and Pauli spin blockade, Phys. Chem. Chem. Phys., 2025, 27, 5238.
81. V. Madhavan, W. Chen, T. Jamneala, M. F. Crommie and N. S. Wingreen, Tunneling into a Single Magnetic Atom: Spectroscopic Evidence of the Kondo Resonance, Science, 1998, 280, 567.
82. D. Goldhaber-Gordon, H. Shtrikman, D. Mahalu, D. Abusch-Magder, U. Meirav and M. A. Kastner, Kondo effect in a single-electron transistor, Nature, 1998, 391, 156.
83. D. Goldhaber-Gordon, J. Gores, M. A. Kastner, H. Shtrikman, D. Mahalu and U. Meirav, From the Kondo Regime to the Mixed-Valence Regime in a Single-Electron Transistor, Phys. Rev. Lett., 1998, 81, 5225.
84. G. Chen, M. S. Dresselhaus, G. Dresselhaus, J. P. Fleurial and T. Caillat, Recent developments in thermoelectric materials, Int. Mater. Rev., 2003, 48, 45.
85. J. N. Hu, X. L. Ruan and Y. P. Chen, Thermal Conductivity and Thermal Rectification in Graphene Nanoribbons: A Molecular Dynamics Study, Nano Lett., 2009, 9, 2730.
86. D. M. T. Kuo and Y. C. Chang, Thermoelectric and thermal rectification properties of quantum dot junctions, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 205321.
87. Z. W. Tan, J. S. Wang and C. K. Gan, First-Principles Study of Heat Transport Properties of Graphene Nanoribbons, Nano Lett., 2011, 11, 214.
88. S. Hershfield, K. A. Muttalib and B. J. Nartowt, Non-linear thermoelectric transport: A class of nano-devices for high efficiency and large power output, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 085426.
89. Y. Xu, Z. Gan and S. C. Zhang, Enhanced thermoelectric performance and anomalous Seebeck effects in topological insulators, Phys. Rev. Lett., 2014, 112, 226801.
90. P. H. Chang, M. S. Bahramy, N. Nagaosa and B. K. Nikolic, Giant thermoelectric effect in graphene-based yopological insulators with heavy adatoms and nanopores, Nano Lett., 2014, 14, 3779.
91. J. Azema, P. Lombardo and A. M. Dare, Conditions for requiring nonlinear thermoelectric transport theory in nanodevices, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 205437.
92. R. S. Whitney, Finding the quantum thermoelectric with maximal efficiency and minimal entropy production at given power output, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 115425.
93. J. H. Jiang and Y. Imry, Enhancing Thermoelectric Performance Using Nonlinear Transport Effects, Phys. Rev. Appl., 2017, 7, 064001.
94. V. Talbo, J. Saint-Martin, S. Retailleau and P. Dollfus, Non-linear effects and thermoelectric efficiency of quantum dot-based single-electron transistors, Sci. Rep., 2017, 7, 14783.
95. D. M. T. Kuo, C. C. Chen and Y. C. Chang, Large enhancement in thermoelectric efficiency of quantum dot junctions due to increase of level degeneracy, Phys. Rev. B, 2017, 95, 075432.
96. P. A. Erdman, F. Mazza, R. Bosisio, G. Benenti, R. Fazio and F. Taddei, Thermoelectric properties of an interacting quantum dot based heat engine, Phys. Rev. B, 2017, 95, 245432.
97. P. A. Erdman, J. T. Peltonen, B. Bhandari, B. Dutta, H. Courtois, R. Fazio, F. Taddei and J. P. Pekola, Nonlinear thermovoltage in a single-electron transistor, Phys. Rev. B, 2019, 99, 165405.
98. N. Taniquchi, Quantum control of nonlinear thermoelectricity at the nanoscale, Phys. Rev. B, 2020, 101, 115404.
99. G. Marchegiani, A. Braggio and F. Giazotto, Nonlinear Thermoelectricity with Electron-Hole Symmetric Systems, Phys. Rev. Lett., 2020, 124, 106801.
100. B. Bhandari, P. A. Erdman, R. Fazio, E. Paladino and F. Taddei, Thermal rectification through a nonlinear quantum resonator, Phys. Rev. B, 2021, 103, 155434.
101. S. Mishra and C. Benjamin, Reaching the Van den Broeck limit in linear response and the Whitney limit in nonlinear response in edge mode quantum thermoelectrics and refrigeration, Phys. Rev. B, 2024, 110, 235430.
102. P. Senapati and P. Parida, Charge and spin thermoelectric transport in benzene-based molecular nano-junctions: a quantum many-body study, Nanoscale, 2024, 16, 2574.
103. J. Fast, H. Lundstrom, S. Dorsch, L. Samuelson, A. Burke, P. Samuelsson and H. Linke, Geometric Symmetry Breaking and Nonlinearity Can Increase Thermoelectric Power, Phys. Rev. Lett., 2024, 133, 116302.
104. J. Balduque and R. Sanchez, Scattering Theory of Thermal and Bipolar Thermoelectric Diodes, Phys. Rev. Lett., 2025, 134, 186301.
105. R. J. Hartig, I. Grosu and I. Tifrea, Thermoelectric response of nanoscale devices in the nonlinear regime, Phys. E, 2025, 171, 116236.
106. N. Sobrino, Thermoelectric properties of interacting double quantum dots, Phys. Rev. B, 2025, 112, 235101.
107. D. Das Sarma, X. Wang and S. Yang, Hubbard model description of silicon spin qubits: Charge stability diagram and tunnel coupling in Si double quantum dots, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 235314.
108. N. S. Lai, W. H. Lim, C. H. Yang, F. A. Zwanenburg, W. A. Coish, F. Qassemi, A. Morello and A. S. Dzurak, Pauli spin blockade in a highly tunable silicon double quantum dot, Sci. Rep., 2011, 1, 110.
109. D. M. Zajac, A. J. Sigillito, M. Russ, F. Borjans, J. M. Taylor, G. Burkard and J. R. Petta, Resonantly driven CNOT gate for electron spins, Science, 2018, 359, 439.
110. A. C. Johnson, J. R. Petta, C. M. Marcus, M. P. Hanson and A. C. Gossard, Singlet-triplet spin blockade and charge sensing in a few-electron double quantum dot, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 165308.
111. M. Eich, R. Pisoni, A. Pally, H. Overweg, A. Kurzmann, Y. J. Lee, P. Rickhaus, K. Watanabe, T. Taniguchi, K. Ensslin and T. Ihn, Coupled Quantum Dots in Bilayer Graphene, Nano Lett., 2018, 18, 5042.
112. M. Preeti, N. Natthapon, H. A. Eric and L. Heiner, A Nanoscale Standard for the Seebeck Coefficient, Nano Lett., 2011, 11, 4679.
113. M. L. Perrin, Nanoscale Power Plants: Turning Heat into Power with Graphene Ribbons, Swiss National Science Foundation, 2024.

---