Quantum magnetism of topologically-designed graphene nanoribbons

Graphene nanoribbon (GNR) is a planer strip of graphene with extraordinary electronic and physical propertiesYazyev (2010); Son et al. (2006a); Yazyev (2013). It is first introduced as a theoretical model to study the size effect in grapheneNakada et al. (1996). Recent advance in bottom-up techniques have allowed production of atomically precise GNRs with armchair, zigzag and other sophistated edgesKimouche et al. (2015); Ruffieux et al. (2016). GNRs have potential applications in the next-generation nanoelectronics, and have been extensively investigated.

GNRs exhibit electronic properties that are not present in grapheneRuffieux et al. (2016). In the simplest form, the edges can be either armchair or zigzag. Tight binding calculations predict that zigzag GNRs are always metallic, while armchair ones can be insulating or metallic. The graphene features two inequivalent Dirac points, which are characterized by $\pm\pi$ Berry phasesCastro Neto et al. (2009). In the presence of zigzag edges, the two Dirac points are projected to different momenta of the one-dimensional (1D) Brillouin zone, and flat bands connecting them appear. The low-energy states on the flat bands are localized at the edge. The large density of states at the Fermi energy is sensitive to the electron-electron interactions, and the nanoribbon has an instability to magnetic ordering. Remarkable edge magnetism develops in the weak-coupling regime $U/t<2$, where bulk magnetic order is absentSorella and Tosatti (1992); Paiva et al. (2005); Meng et al. (2010); Sorella et al. (2012). The magnetic moments form long-ranged ferromagnetic order along the zigzag edge while they are antialigned at the opposite edgeFeldner et al. (2011); Hikihara et al. (2003); Feldner et al. (2010); Luitz et al. (2011); Schmidt and Loss (2010); Golor et al. (2013); Hagymási and Legeza (2016). The opposite spin polarization along the zigzag edges opens a band gap, and the system becomes a semiconductor. Recent experiments provide direct evidence of the edge magnetismJoly et al. (2010); Magda et al. (2014). In particular, it is found that the magnetic order on zigzag edges can be stable even at room temperatureMagda et al. (2014).

The gap values of the armchair GNRs depend on width. When the number of dimer lines across the ribbon is $N=3p+2$ with $p$ a positive integer, the armchair GNRs are metallic. Otherwise, they are insulating with the gaps inversely proportional to the widthsSon et al. (2006b); Brey and Fertig (2006); Ezawa (2006). Recently, the topological properties of the insulating armchair GNRs are investigated, which demonstrates that different widths and end terminations lead to distinct topological phasesCao et al. (2017); Lin and Chou (2018). The topology is characterized by a $\mathbb{Z}_{2}$ invariant, and is manifested by localized boundary state between two segments with different topologies. Later experiments confirming the predictions are reportedGroning et al. (2018); Rizzo et al. (2018). Furthermore, the junction states are successfully used as building block to engineer the famous Su-Schrieff-Heeger topological modelGroning et al. (2018); Rizzo et al. (2018); Su et al. (1979).

The armchair GNRs have no edge states, and the critical interaction to antiferromagnetism is similar to that of bulk graphene. In the topologically-designed GNRs, there are zigzag-edge segments at the boundary, which suggests the possibility of magnetic ordering at weak interactions. While the mean-field (MF) calculation finds a critical interaction smaller than that of bulk system ($U_{c}=2.23t$)Sorella and Tosatti (1992), it is still larger than the realistic value in graphene, and supports that the experimentally realized structures are nonmagnetic. Since the quantum fluctuation is strong in two dimensions, it is very necessary to perform an exact numerical study of the effect of interactions, and validate the absence of the magnetism.

In the manuscript, we study quantum magnetism of topologically-designed GNRs based on the Hubbard models using exact numerical simulations. We first derive a two-band model for the low-energy topological bands, and study the corresponding 1D Hubbard model using DMRG and DQMC methods. It is found the spin correlations decay quickly with the distance, and the local moment is extrapolated to zero in the presence of symmetry-breaking terms. The results show that the two-band Hubbard chain is nonmagnetic. We then include the Hubbard interaction to the topological-designed GNRs, and calculate the spin correlations and the local moments in the presence of symmetry-breaking terms. We find the magnetic order develops from a critical interaction, which is much larger than the realistic value in graphene. We conclude the experimentally relevant GNRs is nonmagnetic, which is consistent with the experimental results.

We consider a nearest-neighbor (NN) tight-binding Hamiltonian on a topologically-designed GNR with the geometry shown in Fig. 1,

where $c_{j\sigma}$ and $c^{{\dagger}}_{j\sigma}$ are the annihilation and creation operators at site $j$ with spin $\sigma=\uparrow,\downarrow$. The hopping amplitudes between the NN sites $l$ and $j$ are $t$, which is about $3~{}eV$ in graphene.

The geometry in Fig. 1(b) is obtained by adding small segments with zigzag edges periodically on the pristine armchair GNR with the width $N=7$ ($7$-AGNR). The resulting geometry is a superlattice of alternating unit cells of $7$-AGNR and $9$-AGNR. The topological property of the $7(9)$-AGNR composed of such unit cells [see shaded region in Fig. 1(b)] is described by a $Z_{2}$ invariant: $Z_{2}=\frac{1\pm(-1)^{\lfloor\frac{W}{3}\rfloor+\lfloor\frac{W+1}{2}\rfloor}}{2}$, where $\lfloor x\rfloor$ is the floor function and the sign is negative (positive) for $W=7(9)$Cao et al. (2017). The $Z_{2}$ value is $0(1)$ for $7(9)$-AGNR, and the topological and trivial unit cells alternate along the chain. Thus the localized boundary states between adjacent unit cells may form low-energy dispersing bands in the gap.

Performing Fourier transformation, the Hamiltonian in Eq.(1) becomes $H_{0}=\sum_{\bf k}\psi^{\dagger}_{\bf k}{\cal H}({\bf k})\psi_{\bf k}$, where the basis is $\psi_{\bf k}=(c_{1,{\bf k}},c_{2,{\bf k}},...,c_{N_{s},{\bf k}})^{T}$ ($N_{s}$ is the number of sites in a unit cell), and ${\cal H}({\bf k})$ is the $N_{s}$-by-$N_{s}$ Hamiltonian matrix in the momentum space. By diagonalizing ${\cal H}({\bf k})$, the band structures are directly obtained, and are shown in Fig. 2. Figure 2(a) shows the band structure of the pristine $7$-AGNR for comparison, which has been known to be a semiconductor with the gap size about $0.5t$. By including small segments, two new dispersing bands appear in the gap[see Fig. 2(b)], and their wave functions mainly distribute on the zigzag sites of the small segments.

Since the two low-energy bands in Fig. 2(b) are isolated from other bands, an effective model for them can be constructed using maximally localized Wannier functionsMarzari et al. (2012), which is a 1D chain containing hoppings of different ranges. The main contribution is from the NN ones, and the corresponding hopping amplitudes are $t_{1}=0.1493t$ and $t_{2}=-0.1824t$, respectively. It is found that the low-energy two bands are fitted quite well only using the above NN hoppings (see Fig. 2). Thus the designed GNR realizes an effective 1D tight-binding model analog to the famous Su-Schrieffer-Heeger one.

It has been known that the pristine $7$-AGNR supports a topological phase. The topological invariant depends on the shape of its termination, and the zigzag termination yields a nontrivial topological invariant. It is also desirable to know the topological class of the above topologically-designed GNR. The 1D topological property is characterized by the Berry phaseGuo and Shen (2011); Guo et al. (2012); Wang and Zhang (2014); Resta (1994); Xiao et al. (2010),

where $k$ varying from $0$ to $2\pi$ is the wave vector, and $\psi_{k}$ is the periodic part of the Bloch wave function corresponding to $k$. The Berry phase of the Hamiltonian Eq.(1) on the geometry shown in Fig. 2 has a nontrivial value $\pi$, corresponding to which a pair of zero boundary modes appear under open boundary condition.

To study the magnetic property of the topologically-designed GNR, we consider the Hubbard interaction, which writes as,

where $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is the number operator of fermion and $U$ is the strength of the on-site interaction. In the following discussions, we include the Hubbard term into both the Hamiltonian in Eq.(1) and the two-band effective model, and study the induced magnetism.

The interacting Hamiltonian is solved numerically by means of the DQMC and DMRG methodsWhite et al. (1989); Hirsch (1985); White (1992). In the DQMC approach, one decouples the on-site interaction term through the introduction of an auxiliary Hubbard-Stratonovich field (HSF). The fermions are integrated out analytically, and then the integral over the HSF is performed stochastically. The only errors are those associated with the statistical sampling, the finite spatial lattice and inverse temperature discretization. All are well-controlled in the sense that they can be systematically reduced as needed, and further eliminated by appropriate extrapolations. The DMRG method is most efficient for the ground state of 1D systemsBauer et al. (2011). We use it to solve the effective model directly at zero temperature. In the following calculations, we focus on half-filled bands, when the DQMC method is free of the infamous ‘minus-sign problem’.

To study the magnetic behavior, we measure the local moment $m^{z}_{j}=\langle S^{z}_{j}\rangle$ with $S^{z}_{j}=\frac{1}{2}(n_{j\uparrow}-n_{j\downarrow})$. Since the original Hamiltonian preserves the $SU(2)$ symmetry, a symmetry-breaking term should be included to induce the magnetism. We break the symmetry of the magnetic phase to $z-$axis by adding an alternating Zeeman term $H_{B}=\sum_{j}B_{j}(n_{j\uparrow}-n_{j\downarrow})$ with $B_{j}=\pm B$ depending on the sublattices. The value of the local moment is extrapolated to the limit $B=0$, and a nonzero value marks the existence of the magnetism. The equal-time spin correlation function is also calculated, which is given by

where the spin raising and lowing operators are $S^{+}_{\bf j}=c^{\dagger}_{\bf j\uparrow}c_{\bf j\downarrow},S^{-}_{\bf j}=c^{\dagger}_{\bf j\downarrow}c_{\bf j\uparrow}$, respectively.