Effects of Coulomb Blockade on the Charge Transport through the Topological States of Finite Armchair Graphene Nanoribbons and Heterostructures

In this study, we investigate the charge transport properties of semiconducting armchair graphene nanoribbons (AGNRs) and heterostructures through their topological states (TSs), with a specific focus on the Coulomb blockade region. Our approach employs a two-site Hubbard model that takes into account both intra- and inter-site Coulomb interactions. Using this model, we calculate the electron thermoelectric coefficients and tunneling currents of serially coupled TSs (SCTSs). In the linear response regime, we analyze the electrical conductance (Ge), Seebeck coefficient (S), and electron thermal conductance (κe) of finite AGNRs. Our results reveal that at low temperatures, the Seebeck coefficient is more sensitive to many-body spectra than electrical conductance. Furthermore, we observe that the optimized S at high temperatures is less sensitive to electron Coulomb interactions than Ge and κe. In the nonlinear response regime, we observe a tunneling current with negative differential conductance through the SCTSs of finite AGNRs. This current is generated by electron inter-site Coulomb interactions rather than intra-site Coulomb interactions. Additionally, we observe current rectification behavior in asymmetrical junction systems of SCTSs of AGNRs. Notably, we also uncover the remarkable current rectification behavior of SCTSs of 9-7-9 AGNR heterostructure in the Pauli spin blockade configuration. Overall, our study provides valuable insights into the charge transport properties of TSs in finite AGNRs and heterostructures. We emphasize the importance of considering electron–electron interactions in understanding the behavior of these materials.


Introduction
The study of two-dimensional (2D) materials has gained significant attention since the discovery of graphene in 2004 [1][2][3][4]. Although graphene has limited applications in optical, semiconductor, and thermoelectric devices due to its gapless semi-metal nature, other 2D materials, such as MoS 2 and WSe 2 , have shown potential for use in transistors and optoelectronics due to their direct band gaps. Low electron mobility and high contact resistance are two of the main challenges that need to be overcome to improve the performance of two-dimensional (2D) electronic devices, such as those made of MoS 2 and WSe 2 [5]. In recent years, researchers have explored the use of heterostructures made of 2D materials to improve the performance of quantum devices with ballistic transport [6].
To implement high-power output quantum devices, it is crucial to reduce the high contact resistance of 2D materials caused by contacted metallic electrodes [6]. On the other hand, certain low-power output quantum devices, such as single-electron transistors [7], single-photon emitters [8][9][10], spin-current conversion devices [11], single-quantum-dot heat engines [12], and solid-state quantum bits [13][14][15][16][17][18][19][20][21], require high contact resistances. Hence, 2D material nanostructures with small dielectric constants may have promising applications in low-power quantum devices. The exploration of electronic structures with deep energy levels that are well-separated from the band states in 2D material nanostructures is crucial for their development [22][23][24][25][26][27][28] since these quantum devices are operated on the basis of few discrete states that are well-separated from other continuous states to reduce charge transport density. A variety of graphene nanoribbons (GNRs) have been extensively studied, including armchair GNRs (AGNRs) [22], AGNR heterostructures [23], cove-edged zigzag GNRs [24], and chevron GNRs [25]. These GNRs have demonstrated the ability to exhibit diverse electronic topological states (TSs) based on their width, edge shape, and end terminations. The controllable manipulation of topological invariants in materials is a highly pursued research area. For instance, the utilization of electric fields and lattice strains has been explored to modulate TSs [26][27][28].
Significant progress has been made in the fabrication of graphene nanoribbons (GNRs) and heterostructures using the bottom-up synthesis technique [13][14][15][16][17][18][19][20][21]. The topological states (TSs) of the end zigzag edges of semiconducting AGNRs and interface states of AGNR heterostructures are well-separated from the conduction and valence subbands [19,20]. Therefore, these TSs may have potential for realizing low-power devices operating at room temperature. Due to graphene's small dielectric constant, electron Coulomb interactions are expected to be significant in the TSs of finite AGNRs and heterostructures [29]. Investigating the effects of electron Coulomb interactions on low-power devices made of GNRs is desirable [30][31][32][33]. However, to date, there has been a lack of theoretical and experimental analysis concerning the Coulomb blockade effect in charge transport through serially coupled TSs (SCTSs) of AGNRs and AGNR heterostructures [34][35][36][37][38]. While various systems involving GNRs have exhibited topological phases [22][23][24][25][26][27][28], the utilization of AGNRs and AGNR heterostructures holds particular advantages for the realization of low-power quantum devices and circuits. This advantage arises from their symmetrical end-edge structures, which can be easily connected to line-contacted electrodes using current bottom-up synthesis techniques [38].
This study aims to investigate the charge transport mechanisms of SCTSs in two distinct systems: the end zigzag edge states of finite AGNRs and the topologically protected interface states of AGNR heterostructures coupled to leads, as shown in Figure 1a,b, respectively. A two-site model is employed to accurately represent the electrical conductance spectra resulting from charge transport through their SCTSs. A two-site Hubbard model with intra-and inter-site Coulomb interactions and Green's function techniques are used to reveal the effects of Coulomb blockade on the charge transport of SCTS. In the linear response regime, we calculate the electrical conductance (G e ), Seebeck coefficient (S), and electron thermal conductance (κ e ) of AGNRs' SCTSs. The Seebeck coefficient shows greater sensitivity to many-body spectra than electrical conductance at low temperatures. At high temperatures, electron Coulomb interactions do not affect the optimized values of S. In the nonlinear response regime, we observe a tunneling current exhibiting negative differential conductance (NDC). NDC arises due to inter-site electron Coulomb interactions in the absence of bias-dependent orbital offset. Moreover, in the asymmetrical tunneling junction systems of finite AGNRs, we observe current rectification behavior due to intersite electron Coulomb interactions. Since the wave functions of TSs are far away from the contacted electrodes, 9-7-9 AGNR heterostructures can be readily set up in the Pauli spin blockade configuration. As a result, the tunneling current in a certain applied bias direction is strongly suppressed due to SCTSs that are highly occupied by two electron triplet states. This current rectification feature in the PSB configuration is very useful for spin-current conversion applications [11].  Figure 1. (a) Schematic diagram of an AGNR connected to electrodes. The tunneling rate of electrons between the left (right) electrode and the leftmost (rightmost) atoms of the AGNR is denoted by Γ L (Γ R ), respectively. The temperature of the left (L) and right (R) electrodes is represented by T L and T R , respectively. The charge density of the localized zigzag edge state with energy ε = 9.087 meV is shown for an AGNR with size characterized by (N z , N a ) = (7,36). (b) Schematic lattice structure of AGNR heterostructures formed by three segments. Each segment is characterized by (N z = 9(7), N a = 12). Note that N a = 12 corresponds to three unit cells (u.c.). The charge density of the topological state with energy ε = 0.12156 eV is plotted for a 9/7/9 AGNR heterostructure. The radius of the circle represents the intensity of the charge density.

Calculation Methods
To model the transport properties of finite AGNRs and heterostructures connected to the electrodes shown in Figure 1, it is a good approximation to employ a tight-binding model with one p z orbital per atomic site to describe the electronic states [39][40][41]. The Hamiltonian of the nano-junction system depicted in Figure 1, including two different AGNR structures, can be written as H = H 0 + H AGNR [42], where The first two terms of Equation (1) describe the free electrons in the left (L) and right (R) electrodes. a † k (b † k ) creates an electron with wave number k and energy k in the left (right) electrode. V L k, ,j=1 (V R k, ,j=N z (N a ) ) describes the coupling between the left (right) lead with its adjacent atom in the -th row. The Hamiltonian for AGNRs can be expressed as: where E ,j is the on-site energy for the p z orbital in the -th row and j-th column. Here, the spin-orbit interaction is neglected. d † ,j (d ,j ) creates (destroys) one electron at the atom site labeled by ( , j ) where and j, respectively, are the row and column indices as illustrated in Figure 1. t ( ,j),( ,j ) describes the electron hopping energy from site ( , j ) to site ( , j ). The tight-binding parameters used for AGNRs is E ,j = 0 for the on-site energy and t ( ,j),( ,j ) = t ppπ = 2.7 eV for the nearest neighbor hopping strength.
To study the transport properties of an AGNR junction connected to electrodes, it is convenient to use the Keldysh Green function technique [42]. In the linear response regime, the electrical conductance (G e ), Seebeck coefficient (S), and electron thermal conductance (κ e ) are given by G e = e 2 L 0 , S = −L 1 /(eTL 0 ), and κ e = 1 T (L 2 − L 2 1 /L 0 ), respectively. Here, L n (n = 0, 1, 2) is defined as The Fermi distribution function of electrodes at equilibrium temperature T and chemical potential µ is given by f (ε) = 1/(exp (ε−µ)/k B T + 1). The transmission coefficient T LR (ε), as shown in Equation (3), is a critical factor in electron transport between the left (L) and right (R) electrodes. The numerical code can be used to calculate T LR (ε), which is given by Here, Γ L (ε) and Γ R (ε) represent the tunneling rate (in energy units) at the left and right leads, respectively. Additionally, G r (ε) and G a (ε) correspond to the retarded and advanced Green functions of the AGNR, respectively [43][44][45].

Charge Transport through a Finite Armchair Graphene Nanoribbon (AGNR)
The electronic behavior of AGNRs is determined by their widths, which follow the rule N z = 3m + 2, N z = 3m + 1, and N z = 3m, where m is an integer. AGNRs exhibit semiconducting behavior for N z = 3m + 1 and N z = 3m, resulting in semiconducting phases for widths such as N z = 7, N z = 9, and N z = 13 with corresponding band gaps of 1.26 eV, 0.948 eV, and 0.714 eV, respectively, in the absence of electron Coulomb interactions [40,41]. Notably, AGNRs with N z = 3m + 2 maintain their semiconducting phases, as determined by first-principle calculations in references [22,28]. However, the oneband tight-binding model still captures the main physics for the cases of N z = 3m and N z = 3m + 1, which is the primary focus of this article. To investigate charge transport through AGNRs with line contacted electrodes, as depicted in Figure 1a, we present the calculated electron conductance (G e ) as a function of µ at zero temperature for various N z values of finite AGNRs with N a = 40 in Figure 2. The magnitude of G e is given in units of the quantum conductance of G 0 = 2e 2 h = 1/(12.9 KΩ) = 77.5 µS. For the case of N z = 7 in Figure 2a, we observe a peculiar peak labeled by Σ 0 at the charge neutrality point (CNP) (µ = 0). In contrast, the spectra around the CNP are vanishingly small for N z = 9, as seen in Figure 2b. The uniform peak separation in Figure 2c corresponds to the linear dispersion of AGNRs with a metallic phase at N z = 11. In particular, Σ 0 is split into two peaks around the CNP for N z = 13 in Figure 2d. The electrical conductance near the CNP is attributed to electron transport through the end zigzag edges of finite AGNRs. The wave functions of the left and right zigzag edges of AGNRs decay along the armchair directions, and their overlaps are vanishingly small in infinitely long AGNRs. When N a is finite, the zero mode energy levels are lifted due to the coupling between these two localized wave functions which can be regarded as a SCTS. The magnitude of the lifted energy level can be determined by the electron hopping strength t LR between the left and right zigzag edges [20]. In Figure 2, we have |t LR | = 5.323 meV, 0.1 meV, and 39.67 meV for N z = 7, N z = 9, and N z = 13, respectively. Since |t LR |/Γ t 1, the peak of Σ 0 is degraded in the case of N z = 9. The results of Figure 2 indicate that t LR is significantly affected by the AGNR width, which is either N z = 3m + 1 or N z = 3m. In appendix A, we discuss the electronic structures of 9-7 AGNR superlattices and charge transport through the SCTSs of 9-7-9 AGNR heterostructures, which were proposed as useful quantum bits [20].
To better understand the transport behavior of the SCTSs, Figure 3a shows the calculated electrical conductance (G e ) as a function of µ for various N a at Γ t = 90 meV, T = 0 K and N z = 13. As N a increases for a fixed N z = 13, t LR decreases, with values of 39.67 meV, 22.29 meV, 12.59 meV, and 7.13 meV for N a = 40, N a = 48, N a = 56, and N a = 64, respectively. As seen in Figure 3a, ε HO and ε LU , which are the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), respectively, can be resolved only for N a = 40. For N a = 64, the magnitude of electrical conductance is smaller than 0.4 G 0 . To clarify the effect of contacted electrodes on the transport of the SCTS, we plot the calculated G e for various Γ t values for N a = 64 and N z = 13 in Figure 3b. The electrical conductance reaches the value of one quantum conductance G 0 at Γ t = 27 meV. When Γ t is smaller than 27 meV, the peaks ε HO and ε LU can be distinguished. The results shown in Figure 3 indicate that resolving the peaks of ε HO and ε LU is challenging, as it requires proper modulation of the t LR /Γ t ratio. (a) Electrical conductance G e as a function of chemical potential µ for AGNRs of different lengths, while keeping the tunneling rate fixed at Γ t = 90 meV and the temperature at T = 0 K. (b) Electrical conductance G e as a function of chemical potential µ for AGNRs with N a = 64 (L a = 6.674 nm), while varying the tunneling rate Γ t at T = 0K.
Up to now, it remains unclear about the thermoelectric properties of SCTSs [19][20][21]. In Figure 4, we present the calculated electrical conductance G e , Seebeck coefficient S, power factor PF = S 2 G e , and thermoelectric figure of merit ZT as functions of chemical potential µ for various temperatures for AGNRs with N z = 13 and N a = 64. The spectra of G e exhibit a symmetrical function of µ, while S displays a bipolar behavior with respect to µ. At T = 48 K, the maximum S h(e),max = ±2.753. Notably, S can be well described by S = µ/T at room temperature (T = 324 K), indicating that the thermoelectric properties of zigzag edge states of AGNRs with N z = 13 and N a = 64 resemble those of a single localized state at room temperature [46][47][48]. In Figure 4c, the maximum power factors (PF h(e),max ) are 0.76 and 0.339 for T = 48 K and T = 324 K, respectively. At T = 324 K, PF h(e),max occurs at µ = ±66.68 meV, with the ratio of µ/k B T = ±2. 46. In ref. [46], the authors predicted that the optimized PF of a single energy level without energy level broadening satisfies µ/k B T = ±2.4. In Figure 4d, we calculated ZT using the equation ZT = S 2 G e T/(κ e + κ ph ).
Here, we used a phonon thermal conductance value of κ ph = F s * π 2 k 2 B T 3h , where we assumed a scattering factor of F s = 0.1. For the remainder of this article, we will use the same value of In a previous study [48], the authors proposed introducing defects in GNRs to increase phonon scattering while maintaining the electron transmission coefficient of topological states. More recently, we demonstrated [43] that zigzag edge states in AGNRs with vacancies remained robust. The maximum ZT h(e),max = 1.966 at T = 324 K occurs at µ = ±68.

Effective Hamiltonian and Formula for Tunneling Current
It is challenging to calculate the tunneling current through SCTSs in the Coulomb blockade region using either first-principle methods [22] or tight-binding models [44]. Hence, we introduce a two-site Hubbard model to clarify the thermoelectric coefficients and tunneling current of SCTSs in the Coulomb blockade region. The Hamiltonian of the two-site Hubbard model is given by where E j represents the spin-independent energy level of the TSs, U j = U L(R) = U 0 and U j, = U LR = U 1 denote the intra-site and inter-site Coulomb interactions, respectively, and n j,σ = c † j,σ c j,σ . U 0 and U 1 are calculated using 1 with the dielectric constant 0 = 4, and U cc = 4 eV at i = j. U cc arises from the two-electron occupation in each p z orbital. Ψ L(R) (r i ) is the wave functions of TSs [29].
Based on the effective Hamiltonian given by Equation (4), we derive the tunneling current through the SCTSs coupled to the electrodes after tedious algebra [49]. We obtain the expression of tunneling current leaving from the left (right) electrode as The Fermi distribution function for the α-th electrode is denoted by is the chemical potential of the left (right) electrode with an applied bias of V bias /2 (−V bias /2), such that µ L(R) = µ ± eV bias /2. The transmission coefficient for charge transport through the SCTSs, denoted by T 2−site LR (ε), has a closed-form expression given by where L = ε − E L + iΓ e,L and R = ε − E R + iΓ e,R . Γ e,L(R) is an effective tunneling rate for the left (right) TS coupled to the left (right) electrode. The transmission coefficient given by Equation (6) consists of eight terms, each corresponding to one of the eight possible configurations of the SCTS that an electron with spin σ from the left electrode may encounter. The probabilities of these configurations are as follows: where the intra-site and inter-site two-particle correlation functions are denoted by n ,−σ n ,σ and n ,σ n j,σ ( n ,−σ n j,σ ), respectively. n ,−σ n ,σ n j,σ is the three-particle correlation function. These correlation functions can be solved self-consistently, and it should be noted that probability conservation is satisfied by ∑ m P m = 1.
The expression for T 2−site LR (ε) goes beyond not only mean field theory [38], but also our previous work [49], where we only considered one-particle occupation numbers and intra-site two-particle correlation functions. Theoretical analysis has shown that intersite two-particle correlation functions play a significant role when inter-site Coulomb interactions are large (U 1 > t LR ) [50]. To obtain T 2−site RL (ε), we can simply exchange the indices in Equation (6). Based on Equation (6), the expression for the electrical conductance at zero temperature for a 2-site model without Coulomb interactions is given by: Using Equation (7), we calculated G e for different values of t LR at Γ e,L = Γ e,R = Γ e,t = 22.3 meV in Figure 5a. The considered t LR values correspond to N a = 40, N a = 48, N a = 56, and N a = 64. Similarly, we calculated G e for various Γ e,t values at t LR = 7.13 meV in Figure 5b. It is worth noting that the curves shown in Figure 5a,b are identical to the curves of Figure 3a,b, respectively. These results in Figure 5 illustrate that charge transport through the SCTS can be well explained with a two-site model that introduces an effective tunneling rate Γ e,t . In the following discussion, we will consider charge transport through the SCTSs in the presence of electron Coulomb interactions, both in the linear and nonlinear regimes.

Linear Response Regime
In this subsection, we examine the effects of the Coulomb blockade on the charge transport through the end zigzag edges of AGNRs in the linear and nonlinear response regimes. We utilize Equations (3) and (6) to calculate the correlation functions (CF), electrical conductance G e , Seebeck coefficient S, and electron thermal conductance κ e of the zigzag edge states as functions of µ at T = 1.2 K and Γ e,t = 1 meV, and present the results in Figure 6. Figure 6a shows the electron single particle occupation number N j,σ = n j,σ , which displays four major plateaus corresponding to the values of 1/4, 1/2, 3/4, and 1 (or total electron number N T = ∑ j,σ N j,σ corresponding to 1, 2, 3, and 4). These plateaus arise due to electron Coulomb interactions. In contrast to the two-electron singlet state in each site n L,−σ n L,σ and n R,−σ n R,σ , which occur at a large µ value near 197 meV, the inter-site triplet state n R,σ n L,σ and the inter-site singlet state n R,−σ n L,σ appear near µ ≈ 42 meV. Figure 6b presents the spectra of G e , which can reveal electron hopping strength (t LR = 7.13 meV), inter-site Coulomb interactions (U 1 = 42 meV), and intra-site Coulomb interactions (U 0 = 155 meV) due to low temperature and weak Γ e,t . The asym-metrical electrical conductances labeled by ε HO and ε LU indicate that the probability of P 1 depends on the µ of the electrodes. In the Coulomb gap region (N j,σ = 0.5), there are two tiny peaks labeled by ε 4(7) and ε 2(5) corresponding to P 4(7) and P 2 (P 5 ) channels. Because they are off-resonant channels as E L = E R , their electrical conductances are quite small. In the Pauli spin blockade configuration, they are in resonant channels, which will be discussed later. Figure 6c shows the spectra of Seebeck coefficient S, which can more clearly reveal the two tiny spectra of G e (ε 4(7) and ε 2(5) ) because S = − ∂G e (µ) ∂µ at low temperatures. Figure 6d displays the spectra structure of electron thermal conductance κ e , which is quite similar to G e spectra at low temperatures, whereas they are quite different at high temperatures (see results of Figures 7 and 8). Note that κ e is in units of κ 0 = 0.62 nW/K. To clarify the contact properties in the presence of Coulomb interactions, we then calculated several thermoelectric quantities, including G e , S, κ e , PF, ρ = κ e /(TG e ), and ZT, as functions of µ for various Γ e,t at a temperature of 48K, and plotted the results in Figure 7. Consistent with the results in Figure 5, the maximum G e still occurred at the condition of Γ e,t = t LR , which is not affected by the intra-site and inter-site Coulomb interactions. At T = 48 K, the ε HO and ε LU peaks are washed out in the spectra of G e , whereas the Coulomb gap between ε 2 and ε 3 arising from intra-site Coulomb interactions remains. Although the spectra of S become more complex in the presence of electron Coulomb interactions, it does not significantly affect the maximum S h(e),max . Comparing with the results of S h(e),max at T = 48 K shown in Figure 4, the maximum S h(e),max has a value of ±2.684 at Γ e,t = t LR , which is almost the same as S h(e),max of Figure 4. The spectra of κ e show two main structures, each with three peaks, and an extra peak appears in the charge blockade region (as seen in the G e spectra) as the temperature increases. The maximum PF does not occur at Γ e,t = t LR , as seen in Figure 7d. The Lorenz number of charge transport through the zigzag edges does not satisfy the Wiedemann-Franz law L z = κ e /(TG e ) = π 2 k 2 B 3e 2 in Figure 7e, indicating that κ e and G e may not be relevant thermoelectric quantities in discrete energy level systems. Figure 7f shows that the maximum ZT values prefer smaller Γ e,t and ρ = κ e /(TG e ) < L z .
To effectively apply heat engines [12], it is crucial to clarify the thermoelectric quantities over a wide temperature range. In Figure 8, we depict the behavior of G e , S, κ e , PF, ρ, and ZT as functions of µ for various temperatures at Γ e,t = t LR . It can be observed that the magnitude of G e reduces with increasing temperature. When the temperature is at 200 K and 324 K, the inter-site Coulomb interactions responsible for the structure in G e spectra are washed out. In Figure 8b, the sophisticated spectra of S become an N-shaped curve at room temperature. In Figure 8c, the highest κ e occurs at the middle Coulomb gap when the temperature is at 240 K and 324 K. This behavior can be understood using the expression for electron thermal conductance κ e = 1 T (L 2 − L 2 1 /L 0 ) = L 2 T − S 2 G e T. When µ is located at the middle Coulomb gap, the electron-hole balance requires S to be close to zero. Therefore, κ e is dominated by L 2 T , which generally has a significant contribution at high temperatures based on the thermionic procedure where µ does not align with resonant channels. This explains why κ e reaches its maximum value in the middle Coulomb gap. The maximum power factors are PF = 0.476 and PF = 0.274 for T = 48 K and T = 324 K, respectively. These values are lower than the maximum PF values in Figure 4c, attributed to the reduction in G e resulting from the Coulomb blockade. As seen in Figure 8e, ρ = κ e /(TG e ) displays temperature-dependent behavior. Clearly, the reduction in G e also suppresses the maximum ZT values. Nevertheless, the maximum ZT h(e),max = 1.6 at T = 324 K in the Coulomb blockade region still reaches 80% of ZT h(e),max = 1.966 shown in Figure 4, without considering electron Coulomb interactions.

Nonlinear Response Regime
Recent experimental studies have reported negative differential conductance (NDC) in charge transport through the AGNR heterostructure tunneling junction [37,38]. However, the underlying mechanism for NDC at low bias remains unclear [37,38]. In this study, we propose a novel mechanism for NDC resulting from the inter-zigzag edge Coulomb interactions U 1 (or inter-site Coulomb interactions). By utilizing Equations (5) and (6), we calculate the tunneling current as a function of applied bias V bias for different temperatures, as shown in Figure 9a,b. We considered U 1 values of 42 meV and 0 meV for Figure 9a,b, respectively. In the forward (reversed) bias, the tunneling current is determined by the transmission coefficient of T 2−site LR (ε) (T 2−site RL (ε)). As seen in Figure 9a, at low temperatures (T = 12 K), the tunneling current in the low bias region decreases with increasing V bias . This intriguing behavior demonstrates the NDC phenomenon, which is absent at higher temperatures.
In Figure 9b, the NDC behavior disappears when we artificially set the inter-site Coulomb interactions to zero (U 1 = 0) and only take into account the intra-site Coulomb interactions (U 0 = 155 meV). The plateau of the tunneling current is attributed to the intrasite Coulomb interaction of U 0 , as U 1 is zero. It is worth noting that in Figure 9a,b, we do not consider the effect of bias-dependent orbital offset. Therefore, such NDC characteristics are not a result of off-resonant channels introduced by applied bias [38]. To clarify the NDC behavior shown in Figure 9a, we added two curves in Figure 9a, which were calculated considering the P 1 and P 3 channels of Equation (6), respectively. The P 3 channel plays an important role for |V bias | < 50 mV. However, it should be noted that the tunneling current is dominated by the P 1 channel when |V bias | > 70 mV. From the curves of P 1 and P 3 , we understand that the NDC behavior of the tunneling current is due to the reduction in the probability of the P 3 channel. In Figure A4, we provide the single-particle occupation number and two-site two-particle correlation functions, such as the inter-site triplet states n R,σ n L,σ ( n L,σ n R,σ ) and singlet states n R,−σ n L,σ ( n L,−σ n R,σ ), which determine the probabilities of the resonant channels P 1 and P 3 in T 2−site LR (ε) (T 2−site RL (ε)). The electronic transport behavior can be significantly impacted by the properties of the contact between AGNRs and the electrodes [19,20,[51][52][53][54]. Therefore, it is important to investigate the effect of the parameter Γ e,t on the NDC behavior. In Figure 10a, we present the tunneling current for different values of Γ e,t at T = 12 K, where the variation of Γ e,t corresponds to finite AGNRs or AGNR heterostructures with different contact properties [53]. To analyze the effect of the applied bias on the orbital offset, we adopted the method of Ref. [55], where the bias-dependent energy level of TSs is given by E L(R) = η L(R) eV bias in Figure 10b. As shown in Figure 10a, the tunneling current increases as Γ e,t increases, whereas the NDC region is reduced. Figure 10b shows that the tunneling current is significantly suppressed in the region of |V bias | > 50 mV when the bias-dependent orbital offset is taken into account. We observed that a second NDC region appears at high applied bias for η = 0.16. Up until now, our calculations assumed a symmetrical contact junction with Γ e,L = Γ e,R . However, in Figure 11, we present tunneling current results for Γ e,L = Γ e,R . In addition to NDC, we also observed current rectification behavior. The maximum current values in the low forward and reversed applied bias are denoted as J F,1 and J R,1 , respectively. The ratio of |J R,1 /J F,1 | increases with increasing Γ e,R when Γ e,L = 3 meV. It reaches 1.91 for Γ e,R = 12 meV. As discussed in Figure 9, J F,1 and J R,1 are primarily attributed to the P 3 channel of Equation (6). The probability weight of P 3 is primarily determined by the single-particle occupation number N R,σ = n R,σ (N L,σ = n L,σ ) and the two-site singlet state n R,−σ n L,σ ( n L,−σ n R,σ ). When |V bias | < 150 mV, the intra-site two-particle occupation number n R,−σ n R,σ ( n L,−σ n L,σ ) and two-site three-particle correlation functions are negligible. In the case of Γ e,L = Γ e,R , we have N L,σ (V bias ) = N R,σ (−V bias ) (see Figure A5). In contrast, if Γ e,L < Γ e,R , then N L,σ (V bias ) = N R,σ (−V bias ). Meanwhile, the two-site singlet state (2-site-S) in the reversed bias region is smaller than that in the forward bias region. Therefore, we observe that |J R,1 | is larger than J F,1 . Note that if the inter-site Coulomb interaction U 1 is set to zero, the current rectification is significantly reduced as Γ e,L = Γ e,R . This suggests that the inter-site Coulomb interactions play a crucial role in determining the behavior of negative differential conductance (NDC) and current rectification. Note that the nonlinear current-voltage characteristics (NDC) and current rectification behavior discussed for a finite AGNR can also be observed in an AGNR heterostructure.

Tunneling Current through TSs of 9-7-9 AGNR Heterostructures under Pauli Spin Blockade Configuration
Because the topological states (TSs) of 9-7-9 AGNR heterostructures are located at the interfaces between 9-7 junctions, they can be far away from the contacted electrodes based on the design (refer to Appendix A or Figure 1b). Therefore, it is possible to arrange the two gate electrodes in a way that allows for the modulation of the energy levels of TSs and setting them in the Pauli spin blockade (PSB) configuration. We present the calculated tunneling current through the TSs of 9-7-9 AGNR heterostructures in the PSB configuration as a function of V bias for various Γ e,L = Γ e,R = Γ e,t values at T = 12 K in Figure 12, where each AGNR segment has eight unit cells. Such TSs have U 0 = 125 meV, U 1 = 49 meV and t LR = 7.54 meV. As seen in Figure 12, a remarkable current rectification behavior is observed in the PSB configuration at the symmetrical tunneling rate Γ e,L = Γ e,R . The ratio of |J R,max /J F,max | are 4.242, 2.908, 2.477 and 2.286 for Γ e,t = t LR /3, Γ e,t = 2t LR /3, Γ e,t = t LR and Γ e,t = 4t LR /3, respectively. To understand the rectification behavior shown in Figure 12, we first note that under PSB, the channel P 2 in Equation (6) behaves as a resonant channel in the forward applied bias, and the tunneling current is mainly contributed by P 2 . We therefore focus on the transmission coefficient of P 2 in the PSB configuration, which can be expressed as: Here, the probability of the resonant channel is given by P PSB = N R,σ − n R,σ n L,σ − n R,−σ n R,σ + n R,−σ n R,σ n L,σ (P PSB = N R,σ − n L,σ n R,σ − n L,−σ n R,σ + n L,−σ n L,σ n R,σ ) in the forward (reversed) applied bias. In the small applied bias region |V bias | ≤ 50 mV, the resonant channels are almost unaffected by bias-dependent orbital offset, and P PSB is the main factor determining the magnitude of the tunneling current. For |V bias | ≤ 150 mV, only single occupation number and inter-site singlet and triplet states are important. We show these correlation functions in Figure 13 to clarify the bias-dependent P PSB .
As seen in Figure 13a, charge transport in the PSB region is a two-electron process because of the total electron number 1.5 ≤ N T ≤ 2. In the forward applied bias, N R,σ and N L,σ are in discharging and charging processes, respectively, and their values reach 0.5 at V bias = 150 mV. Due to the large probability of inter-site triplet state (see the curve of 2-site-T), the P PSB is degraded. In contrast, N R,σ and N L,σ in the reversed applied bias are in the charging and discharging processes, respectively. The highly enhanced N R,σ and suppressed inter-site correlation functions (see 2-site-T and 2-site-S) result in a large P PSB . This explains the significant current rectification observed in Figure 12a. With increasing Γ e,t , the magnitude of 2-site-T is degraded for V bias > 0, leading to an enhancement of J F,max and a reduction in |J R,max /J F,max | in Figure 12. In addition, we observe an interesting behavior of phase transformation. The merge together 2-site-T and 2-site-S curves in the reversed applied bias region is splitting into two curves in the forward applied bias. Finally, we examine SCTSs in the PSB configuration for weak coupling strength of t LR values. In Figure 14, we present (a) single-particle occupation numbers N j,σ , (b) inter-site two-particle correlation functions, and (c) tunneling current of 9-7-9 AGNR heterostructures as functions of µ at T = 12 K, t LR = 0.88 meV, Γ e,L = t LR /3, and Γ e,R = 3t LR . To achieve the small coupling parameter of t LR = 0.88 meV, we use 12 u.c. segments to form 9-7-9 AGNR heterostructures. We can disregard the bias-dependent orbital offset between TSs, as the wave functions of TSs are distant from the electrodes and the channel length between the electrodes is significantly extended [55]. In Figure 14a, the total electron number N T is constrained within the range of 1.5 < N T ≤ 2. As observed in Figure 14b, the 2-site-T curve reaches its maximum value of 0.486, while the 2-site-S curve approaches nearly zero in the forward applied bias. When compared to the results in Figure 13, it is apparent that the forward applied bias and weak t LR lead to high occupancy of inter-site triplet states in SCTSs. Consequently, the SCTSs occupied by triplet states generate an extremely small current (J F ) as shown in Figure 14c. The significant current rectification demonstrated in the PSB configuration, as depicted in Figure 14, holds great value for spin-current conversion devices [11]. (a) Single-particle occupation numbers N j,σ , (b) inter-site two-particle correlation functions (CF), and (c) tunneling current of 9-7-9 AGNR heterostructures as functions of µ for a small t LR = 0.88 meV at T = 12 K, Γ e,L = t LR /3, and Γ e,R = 3t LR . We realize the small t LR = 0.88 meV using 12 u.c. segments to form 9-7-9 AGNR heterostructures. Their intra-site and inter-site Coulomb interactions are U 0 = 111 meV and U 1 = 36.97 meV. The energy levels of TSs are E L = −U 1 and E R = −U 0 .

Conclusions
In summary, this study provides an in-depth investigation of the charge transport properties of two distinct graphene nanoribbon (GNR) structures: the end zigzag edges of AGNRs and the topological states of 9-7-9 AGNR heterostructures. Our findings demonstrate that 9-7-9 AGNR heterostructures with deep topological states (TSs) exhibit superior characteristics, such as wave functions of TSs far away from the electrodes, making them more suitable for low-power quantum devices. The two-site model with effective tunneling rates provides an excellent description of the electrical conductance spectra of the serially coupled topological states (SCTSs). Additionally, we analyzed the Coulomb blockade effect on the thermoelectric coefficients and tunneling current of the zigzag edge states using the two-site Hubbard model. Our results show that electron Coulomb interactions have a more significant impact on electrical conductance than on the Seebeck coefficient. The Lorenz number of the zigzag edge states does not satisfy the Wiedemann-Franz law due to their localized characteristic.
We observed negative differential conductance (NDC) in the nonlinear response regime of the tunneling current, attributed to inter-zigzag edge Coulomb interactions. Additionally, we observed current rectification behavior in the end zigzag edges of AGNRs when asymmetrical contacted electrode junctions were present. The tunneling current through SCTSs formed by 9-7-9 AGNR heterostructures in the Pauli spin blockade configuration exhibits remarkable current rectification behavior, even in AGNR heterostructures with symmetrical contacted electrodes. For a weak coupling parameter t LR , the forward tunneling current is almost blocked due to the high occupation of SCTSs by two-electron triplet states. This property is highly useful in spin-current conversion devices. Overall, our study provides valuable insights into the charge transport properties of TSs in finite AGNRs and heterostructures, emphasizing the importance of considering electron-electron interactions in understanding their behavior.

Data Availability Statement:
The data presented in this study are available upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. 9-7-9 AGNR Heterostructures
The results presented in Figure 2b indicate that the wave function of the left zigzag edge state has only a small overlap with that of the right zigzag edge state for the case of N z = 9. Consequently, it is difficult for electrons from the electrodes to be transported through these localized states. Recently, interesting topological states (TSs) have been identified in the electronic structures of 9-7 AGNR heterostructures. In these structures, the wave function of TS accumulates at the interface between the 9-7 junction [17][18][19][20]. However, the charge transport through these serially coupled TSs (SCTSs) has not been clarified [38,51,52]. In the Appendix A, we not only further clarify their electronic structures but also investigate how the contacted electrodes influence the charge transport through the SCTSs.
As depicted in Figure 1b, the 9-7-9 AGNRs are formed by three u.c. segments. The calculated eigenvalues of the AGNR heterostructures for various periods are shown in Figure A1. As seen in the case of (a), two peaks labeled by ε HO and ε LU with a splitting energy of |2t LR | = 0.24312 eV between E c and E v are observed. Here, E c and E v denote the conduction subband minimum and valence subband maximum, respectively. The charge density distribution of ε LU = 0.12156 eV is presented in Figure 1b, which shows that the charge densities are very dilute for lattices close to the end zigzag edges. Such behavior indicates that ε HO and ε LU have weak coupling strengths with the electrodes. As each 7-AGNR segment provides two energy levels, there are 4, 6, 8, and 10 energy levels within the band gap for (b), (c), (d), and (e), respectively, as shown in the results. When the number of 7-AGNR segments approaches infinity, the minibands are formed. To reveal such interesting electronic structures formed by TSs, we show the calculated electronic structures of 9-7 AGNR superlattices (SLs) in Figure A2, where (a) and (b) consider different 9-AGNR segment lengths at fixed three u.c. 7-AGNR segments. Two minibands are formed near the charge neutrality point (CNP), the tiny gaps in (a) and (b) are 76 meV and 68 meV, respectively. We find that the electronic dispersions of minibands can be well described by the Su-Schrieffer-Heeger (SSH) model with a closed expression of E SSH (k) = ± t 2 1 + t 2 2 − 2t 1 t 2 cos(k π/L) [16,19,54]. Using E SSH (k), we can determine t 1 = 0.167 eV and t 2 = 0.129 eV in (a) and t 1 = 0.118 eV and t 2 = 0.084 eV in (b). As seen in Figure A2, the miniband widths become narrow as the length of the 9-AGNR segment increases. According to the SSH model, the metal, semiconductor, and insulator phases are determined by the relationships between t 1 and t 2 . Once t 1 = t 2 , the band gap is vanishingly small.
The charge densities of the topological states (TSs) in 9-7-9 AGNR heterostructures are distinct from those of the end zigzag edge states found in finite AGNRs. The wave functions of the TSs are located far away from the line-contacted electrodes, which allows for the observation of charge transport through the TSs of 9-7-9 heterostructures even at large coupling strengths between the zigzag edge carbon atoms and the electrodes (i.e., large Γ t values). We calculated the electrical conductance of 6-unit cell (u.c.) 9-7-9 AGNRs for various Γ t values, and the results are presented in Figure A3. The splitting energy between the highest occupied (ε HO ) and lowest unoccupied (ε LU ) states is |2t LR | = 44.4 meV, and their magnitudes and widths increase with increasing Γ t . Notably, even when Γ t = 2.7 eV, the two peaks of ε HO and ε LU are still clearly resolved. The excellent agreement between the results of Figure A3 and Equation (7) (not shown here) is due to the TSs' localized wave functions. Therefore, a 2-site Hubbard model is also suitable for describing charge transport through the TSs of 9-7-9 AGNR heterostructures.  Figure A1. Eigenvalues of AGNR heterostructures with varying numbers of 7-AGNR segments. 9-7-9 AGNR heterostructures are formed by 3 u.c segments. (a-e) correspond to one to five 7-AGNR segments, respectively, with an energy range of |E| ≤ 2 eV.  Figure A2. Electronic structures of 9-7 AGNR superlattices (SLs). (a,b) correspond to 9-AGNR segments with 4 u.c. and 6 u.c., respectively, with an energy range of |E| ≤ 1.2 eV. Note that the 7-AGNR segment with 3 u.c. is fixed. L is the length of the super unit cell. We have L = 7 u.c. and L = 9 u. c. in (a,b), respectively.

Appendix B. Correlation Functions of the End Zigzag Edges of AGNRs
In Figure 9, we illustrated the tunneling current of the end zigzag edge states as a function of V bias for various temperatures. The behavior of the tunneling current is determined by the resonant channels and their probabilities. Since E L and E R do not have a bias-dependent offset, the NDC behavior arises from the bias-dependent probabilities. Therefore, it is essential to understand the behavior of P n in Equation (6). Figure A4 displays the single-particle occupation number N L(R),σ = n L(R),σ and inter-site twoparticle correlation functions (CF) as functions of the applied bias (V bias ), with and without inter-site Coulomb interactions. Due to the symmetrical structure of the junction system, we have N L,σ (V bias ) = N R,σ (−V bias ). The 2-site singlet (2-site-S) state n R,−σ n L,σ and 2-site triplet (2-site-T) state n R,σ n L,σ (or n L,−σ n R,σ and n L,σ n R,σ ) also exhibit this symmetric character with respect to V bias . We analyzed their behavior in the case of a forward bias situation (V bias > 0). N L(R),σ has a finite value at zero bias (V bias = 0 mV) since E L and E R are below µ. As V bias increases, N L,σ increases, but N R,σ decreases. Eventually, N R,σ will saturate. The charge filling of E R by the left electrode is attributed to the resonant channel between E L and E R . Since N R,σ is indirectly charged by the left electrode, in the forward applied bias range, N R,σ is smaller than N L,σ . We found that the probability of 2-site-S is larger than that of 2-site-T at finite bias. The bias-dependent probabilities of P 1 and P 3 channels, which provide resonant energy levels, are plotted in Figure A4a. The decline behavior of the probability of P 3 explains the behavior of the tunneling current with NDC shown in Figure 9. To understand the current rectification behavior observed in Figure 11, we present N L(R),σ and the inter-site two-particle correlation functions (2-site-S and 2-site-T) for various Γ e,R values at Γ e,L = 3 meV in Figure A5. As shown in Figure A5, the condition of N L,σ (V bias ) = N R,σ (−V bias ) is lifted when the junctions have asymmetrical contacted properties. P 3 channels dominate the behavior of the tunneling current at low bias. Therefore, the forward maximum current (J F,1 ) and the reversed maximum current (J R,1 ) can be determined by P 3 = N R,−σ − n R,−σ n L,σ and P 3 = N L,−σ − n L,−σ n R,σ , respectively. In the low bias region, we have N L,σ ≈ N R,σ and n R,−σ n L,σ > n L,−σ n R,σ , which explains why J R,1 is larger than J F,1 and the current rectification behavior is observed in Figure 11.