Carrier-induced Antisymmetric-symmetric Tendencies of Spin Stiffness in Zigzag Graphene Nanoribbons

TL;DR

This paper investigates how doping influences the spin stiffness in zigzag graphene nanoribbons, revealing antisymmetric-symmetric tendencies and their dependence on ribbon width, using the generalized Bloch theorem and spin spiral formations.

Contribution

It introduces a method to analyze doping effects on spin stiffness in nanoribbons using spin spiral formations and the generalized Bloch theorem, highlighting antisymmetric-symmetric tendencies.

Findings

High spin stiffness confirmed by previous supercell calculations.

Doping induces antisymmetric-symmetric tendencies in spin stiffness.

Ribbon width affects the spin stiffness dependence on doping.

Abstract

The generalized Bloch theorem was applied to calculate the spin stiffness and to consider its tendencies when introducing the doping in zigzag graphene nanoribbons. To reach the intentions, two different flat spin spiral formations were constructed by fixing the ferromagnetic and antiferromagnetic spin arrangements at the two different edges by applying a constraint scheme method. A spin stiffness was then calculated by means of a quadratic order function, which maps the total energy difference in the self-consistent calculations onto the Heisenberg Hamiltonian. We found a very high spin stiffness, as predicted previously by the supercell calculation. We also showed that the antisymmetric-symmetric tendencies of spin stiffness are induced by the hole-electron doping. The dependence of ribbon widths of zigzag graphene nanoribbon on the spin stiffness was also provided with similar tendencies when the doping is taken into account.

Tables

Table 1. Table 1: The spin stiffness dependence of the ribbon width for the nondoped case.

Width	$D_{\text{FM}}$ (meVÅ2)	$D_{\text{AFM}}$ (meVÅ2)
6	2465	3320
10	2966	3583
12	3175	3745

Equations

$${\mathbf{M}}_i(t)$$

$$⟨\mathbf{r}|{\psi }_{\nu \mathbf{k}}⟩={\psi }_{\nu \mathbf{k}}\left(\mathbf{r}\right)$$

$$E$$

$$J_{\mathbf{q}}^{\mu \nu }=-\sum _m(1-{\delta }^{0m}{\delta }^{\mu \nu })J^{\mu (\nu +m)}{e}^{-i\mathbf{q}({\mathbf{R}}_{\nu +m}-{\mathbf{R}}_{\mu })}.$$

$$J_{\mathbf{q}}^{\mu \nu }=J_{-\mathbf{q}}^{\nu \mu },\text{Re}[J_{\mathbf{q}}^{\mu \nu }]=\text{Re}[J_{\mathbf{q}}^{\nu \mu }],\text{Im}[J_{\mathbf{q}}^{\mu \nu }]=-\text{Im}[J_{\mathbf{q}}^{\nu \mu }].$$

$$E=E_0+\frac{1}{2}\sum _{\mu \nu }{M}_{\mu }{M}_{\nu }\left\{(1-{\theta }_{\mu }^2/2-{\theta }_{\nu }^2/2)J_0^{\mu \nu }+{\theta }_{\mu }{\theta }_{\nu }\text{Re}[J_{\mathbf{q}}^{\mu \nu }]\right\}.$$

$$\text{Re}[{\tilde{J}}_{\mathbf{q}}^{\mu \nu }]={\left[\frac{1}{M_{\mu }{M}_{\nu }}\frac{{\partial }^{2}E}{\partial {\theta }_{\mu }\partial {\theta }_{\nu }}\right]}_{\theta =0},$$

$$\text{Re}[{\tilde{J}}_{\mathbf{q}}^{\mu \nu }]=J_{\mathbf{q}}^{\mu \nu }-{\delta }^{\mu \nu }\sum _{\alpha }\frac{M_{\alpha }}{M_{\mu }}{J}_0^{\alpha \mu }.$$

$$i\hslash \frac{d{\hat{\text{M}}}_i}{dt}=\left[{\hat{\text{M}}}_i,\hat{H}\right],$$

$$\left[{\hat{\text{M}}}_i^{\alpha },{\hat{\text{M}}}_j^{\beta }\right]=i{\mu }_B{\delta }_{ij}{ϵ}_{\alpha \beta \gamma }{\hat{\text{M}}}_i^{\gamma }.$$

$${\hat{\text{M}}}_i(t)={\mu }_B{\int }_{V_i}{d}^3\mathbf{r}{\psi }^{†}(\mathbf{r},t)\hat{\mathbf{\sigma }}\psi (\mathbf{r},t),$$

$$\hslash \sin {\theta }_i\frac{d{\phi }_i^{\mathbf{q}}}{dt}={\mu }_B\sum _{j(\ne i)}{J}_{ij}{M}_j[\sin {\theta }_i\cos {\theta }_j-\cos {\theta }_i\sin {\theta }_j\cos ({\phi }_i^{\mathbf{q}}-{\phi }_j^{\mathbf{q}})],$$

$$\hslash \frac{d{\theta }_i}{dt}=-{\mu }_B\sum _{j(\ne i)}{J}_{ij}{M}_j\sin {\theta }_i\sin ({\phi }_i^{\mathbf{q}}-{\phi }_j^{\mathbf{q}}),$$

$$\sqrt{M_{\mu }}{\theta }_{\mu }\hslash {\omega }_{\mathbf{q}}={\mu }_B\sum _{\nu }\sqrt{M_{\mu }{M}_{\nu }}\text{Re}[{\tilde{J}}_{\mathbf{q}}^{\mu \nu }]\sqrt{M_{\nu }}{\theta }_{\nu }.$$

$$\text{det}({\delta }_{\mu \nu }\hslash {\omega }_{\mathbf{q}}-{\mu }_B\sqrt{M_{\mu }{M}_{\nu }}\text{Re}[{\tilde{J}}_{\mathbf{q}}^{\mu \nu }])=0.$$

$$\hslash {\omega }_{\mathbf{q}}$$

$$\hslash {\omega }_{\mathbf{q}1}={\mu }_{B}M\text{Re}[{\tilde{J}}_{\mathbf{q}}^{11}],\hslash {\omega }_{\mathbf{q}2}={\mu }_{B}M\text{Re}[{\tilde{J}}_{\mathbf{q}2}^{22}].$$

$$E=E_0+M^2(1-{\theta }^2)J_0^{11}+M^2{\theta }^2\text{Re}[J_{\mathbf{q}}^{11}].$$

$$\text{Re}[{\tilde{J}}_{\mathbf{q}}^{11}]=\frac{1}{2M^2}\frac{{\partial }^{2}E}{\partial {\theta }^2}.$$

$$\hslash {\omega }_{\mathbf{q}}=\frac{{\mu }_B}{M}\frac{\Delta E(\mathbf{q},\theta )}{{\sin }^2\theta },$$

Full text

Keywords: spin-wave excitations, spin stiffness, graphene nanoribbon

1 Introduction

The exploration of graphene nanoribbons (GNR) is of interest for the spintronics applications. Present studies confirmed that both zigzag graphene nanoribbons (ZGNR) and armchair graphene nanoribbons (AGNR) are the feasible candidates for the spintronics application by applying the electric field [1, 2, 3, 4], fixing the angle between two different magnetic moment of C atoms at the different edges [5], or introducing the doping [6, 7, 8, 9, 10]. The rare property that has not been fully elucidated is the spin-wave excitations in the GNR. Edwards and Katsnelson predicted that the $sp$ electron systems should have higher spin stiffness than that of 3$d$ transition metals [11]. They analyzed by comparing the dominant contribution between the Stoner excitations and the spin-wave excitations. Previous studies considered the spin-wave excitations using the supercell calculation [12] or the Hubbard model approach [13, 14, 15] in the zigzag graphene nanoribbons (ZGNR).

The aim of this paper is to calculate the spin stiffness and to investigate the tendencies of spin stiffness with respect to the hole-electron doping in the ZGNR by using the generalized Bloch theorem (GBT) within the linear combination of pseudo-atomic orbitals (LCPAO) method. For calculating the spin stiffness, we are given a choice to involve the Heisenberg model or not. In fact, some authors found the reliable spin stiffness in some magnetic systems without using the Heisenberg model [16, 17], i.e., the spin stiffness is achieved by fitting the total energy difference excluding the magnetic moment. However, for the most cases, such as in 3$d$ transition metals, the Heisenberg model should be used to get the reliable results [18, 19, 20, 21]. For this situation, Yazyev and Katsnelson proved that the Heisenberg model should be included [12]. When they ignored the Heisenberg model, they found the spin stiffness of 320 meVÅ2. Based on Ref. [11], this value seems less reliable if comparing to the spin stiffness in 3$d$ transition metals. For the comparison, we take fcc nickel with the magnetic moment per atom about 0.6 $\mu_{B}$. The magnetic moment of carbon atom at the edge is about 0.28 $\mu_{B}$, two times smaller than that of fcc nickel. The DFT calculations showed that the spin stiffness of fcc nickel is larger than 700 meVÅ2 [18, 19, 20, 21], so the spin stiffness of ZGNR should be much larger than that of fcc nickel. Yazyev and Katsnelson then obtained a high value of spin stiffness of 2100 meVÅ2 after considering the Heisenberg model [12], which is a more reliable value. The important point is that this high value is related to the allowable length of spin correlation at the low temperature, at which the ZGNR-based devices can work. This means that controlling the spin stiffness, in this case, has a vital role to develop ZGNR-based spintronics applications.

In this paper, we include the Heisenberg model to calculate the spin stiffness based on the above reason. To employ the GBT, we construct two flat spiral configurations with the constraint scheme method to fix the direction of the magnetic moments of C atoms at two different edges. The frozen magnon method is applied to calculate the magnon energy of the system without considering the spin-orbit interaction. We also show that the magnitude of the magnetic moment is drastically reduced for the high spiral vector. Therefore, the calculation of spin stiffness can only be performed in a set of the lower spiral vectors than that of 3$d$ transition metals, as carried out in our previous paper [22]. We obtain the high spin stiffness for the ferromagnetic and antiferromagnetic spin arrangements at two different edges and find that the tendencies of spin stiffness become different when introducing the doping. The ferromagnetic configuration achieves the antisymmetric tendency while the antiferromagnetic one yields the symmetric tendency.

Note that the calculated spin stiffness via DFT calculation strongly needs verification from the experiment. In the case of the low magnetic moment, such as fcc nickel and ZGNR, the Stoner excitations cannot be neglected. In 3$d$ transition metals, the calculated spin stiffnesses are in good agreement with experiments only in bcc iron and fcc cobalt, whereas the overestimated spin stiffnesses occur in fcc nickel. This situation may be hold for the ZGNR case [12]. This means that the calculated spin stiffness of ZGNR via DFT calculation could be either underestimated or overestimated with the available experiment. In this paper, the rest of the discussions is organized as follows. The model of ZGNR and the mathematical expression to calculate the spin stiffness will be discussed thoroughly in section 2. We also give the available considerations to simplify the calculation based on the distance of two magnetic carbon atoms at the edges. The calculated spin stiffness and its tendencies when the doping is taken into account will be given in section 3. We also vary the ribbon width of ZGNR to see the tendencies of spin stiffness. We show that the spin stiffness is proportional to the ribbon width of ZGNR. We summarize our discussion based on the results in section 4.

2 Model of Ferromagnetic and Antiferromagnetic Zigzag Graphene Nanoribbons

To consider the spin-wave excitations using the GBT, the magnetic moment of atoms should be rotated under the direction of a spiral vector $\mathbf{q}$ to establish a spin spiral configuration [23], which is given by

[TABLE]

This magnetic moment is characterized by two angles, i.e., the cone angle $\theta$ and the azimuthal angle $\varphi^{\mathbf{q}}$. To consider the spin-wave excitations, the cone angle $\theta$ remains unchanged while the azimuthal angle $\varphi^{\mathbf{q}}$ rotates in the direction of spiral vector $\mathbf{q}$ defined as $\varphi_{i}^{\mathbf{q}}(t)=\mathbf{q}\cdot\mathbf{R}_{i}+\omega_{\mathbf{q}}t$, where ${R}_{i}$ is the lattice vector and $\omega_{\mathbf{q}}$ is the frequency of magnon. The first-principles calculation in this paper was carried out by using the OpenMX code [24], which applies the LCPAO method and the norm-conserving pseudopotentials [25]. The implementation of GBT in the OpenMX code was conducted by writing the Bloch wavefunctions on the atomic site $\tau_{i}$ in terms of an LCPAO as [22]

[TABLE]

where the pseudo-atomic orbital (PAO) $\phi_{i\alpha}$ is constructed by the confinement scheme method [26, 27]. For applying the GBT, we first build the model of ferromagnetic and antiferromagnetic configurations at two different edges of ZGNR, as shown in Figs. 1(a) and 1(b), where the structural optimization was carried out by using the nonmagnetic state. The flat spiral configuration will be then generated by fixing the cone angle $\theta=\pi/2$ with the applied penalty functional in the constraint scheme method. The penalty functional works if the direction of the magnetic moment deviates away from the initial $\theta$ in the self-consistent calculation. This functional will force the deviated direction of the magnetic moment back to the initial condition, for the detailed explanations on the constraint scheme method one can see Refs. [28, 29, 30]. So, the total energy and the magnetic moment will then be evaluated self-consistently with a fixed $\theta$, where $\varphi$ rotates in the $\mathbf{q}$ direction.

We assume that the total energy of the system can be written in terms of Heisenberg Hamiltonian as

[TABLE]

where $N$ defines the number of unit cells and $E_{0}$ is the ground state energy. To consider the periodicity in our system, we rewrite the lattice vector ${R}_{i}$ with the double index notation ${R}_{m+\mu}=T_{m}+\tau_{\mu}$ as the linear combination of the translation vector $T_{m}$ and atomic site $\tau_{\mu}$ [31], see Fig. 1(c). Since the frozen magnon method works in the reciprocal space, it is convenient to transform the exchange coupling constant $J_{ij}$ in the real space to the reciprocal space by means of the Fourier transform

[TABLE]

From Eq. (9), it can be derived the following symmetry properties of the exchange coupling constant in the reciprocal space

[TABLE]

Applying the small deviation of $\theta$ ($\theta<<$), the total energy in Eq. (8) then becomes [32]

[TABLE]

Since we deal with the small $\theta_{\mu}$, we can derive the exchange coupling constant in terms of second derivative of the total energy with respect to $\theta_{\mu}$

[TABLE]

where we define

[TABLE]

The next step is to formulate the eigenvalue equation to calculate the magnon energy. We consider the quantum version of the Heisenberg motion for the operator of magnetic moments.

[TABLE]

with the commutation relations

[TABLE]

In this case, the Greek and Roman indexes refer to the Cartesian coordinates and the lattice sites. Note that the commutation relation in Eq. (15) is referred to the definition of the operator of magnetic moment at the lattice site $i$

[TABLE]

where $V_{i}$ is the crystal volume and $\hat{\bm{\sigma}}$ denotes the operator of Pauli matrices. Extracting Eq. (14) and Eq. (15) in terms of the components of operator of magnetic moment, we obtain the following equations

[TABLE]

as well as the time-independent of the magnitude of magnetic moment $dM_{i}/dt=0$. By applying the small $\theta$ and the Fourier transform to $J_{ij}$ in Eq. (9), those equations above yield the eigenvalue problem

[TABLE]

We can deduce that Eq. (19) corresponds to the eigenvector $\sqrt{M_{\mu}}\theta_{\mu}$ and the eigenvalue $\hbar\omega_{\mathbf{q}}$ obtained by diagonalizing the matrix $\mu_{B}\sqrt{M_{\mu}M_{\nu}}\textrm{Re}[\tilde{J}_{\mathbf{q}}^{\mu\nu}]$. Thus, to obtain the magnon energy, we have to solve the following eigenvalue problem

[TABLE]

As directly seen in Figs. 1(a) and 1(b), the two models of ZGNR contain two magnetic carbon atoms (two sublattices) in the unit cell with the same magnetic moment $M_{1}=M_{2}=M$. In this case, the ferromagnetic configuration indicates the initial angles $\theta_{1}=\theta_{2}=\pi/2$ and $\varphi_{1}^{0}=\varphi_{2}^{0}=0$, while the antiferromagnetic one indicates the initial angles $\theta_{1}=\theta_{2}=\pi/2$ and $\varphi_{1}^{0}=0$, $\varphi_{2}^{0}=\pi$. Consequently, the eigenvalues of the magnon energy attained from Eq. (20) are given by

[TABLE]

where $\textrm{Re}[\tilde{J}_{\mathbf{q}}^{11}]$ is the exchange interaction between C atoms in one edge, $\textrm{Re}[\tilde{J}_{\mathbf{q}}^{22}]$ is the exchange interaction between C atoms in the other edge, and $\textrm{Re}[\tilde{J}_{\mathbf{q}}^{12}]$ is the exchange interaction between C atom in one edge and C atom in the other edge. Since the distance of two magnetic carbon atoms at the different edges in the unit cell is much larger than the length of periodic unit cell for the two configurations in Figs. 1(a) and 1(b), we are allowable to approximate Eq. (21) to obtain two magnon energies

[TABLE]

Since two magnetic C atoms, in this case, are equivalent, $\textrm{Re}[\tilde{J}_{\mathbf{q}}^{11}]$ and $\textrm{Re}[\tilde{J}_{\mathbf{q}}^{22}]$ become equivalent, too. Thus, we can reduce the computation by performing one calculation to obtain the magnon energy, as carried out by Essenberg $et$ $al$ [31]. In this case, we make a different approach with theirs for obtaining the final expression of magnon energy. Realizing $\theta_{1}=\theta_{2}=\theta$, Eq. (11) then becomes

[TABLE]

We then differentiate Eq. (23) with respect to $\theta$ to obtain

[TABLE]

By replacing the left-hand side using Eq. (22) and replacing the second derivative of total energy in terms of the total energy difference at the right-hand side, Eq. (24) is cast into

[TABLE]

where $\Delta E(\mathbf{q},\theta)=E(\mathbf{q},\theta)-E(\mathbf{0},\theta)$ is the total energy difference. Note that, one can get Eq. (24) by employing the L’H$\hat{\textrm{o}}$pital’s theorem at the right-hand side of Eq. (25). Although the derivation employs the small $\theta$, the similar magnon dispersion relation formulated in Eq. (25) can also hold for the larger $\theta$, for example see Refs. [33, 34, 35].

3 Results

The computation was performed by a $90\times 1\times 1$ $k$ point sampling in the Brillouin zone with the generalized gradient approximation (GGA) for the exchange-correlation potential [36]. For the basis sets, we used two valence orbitals $s$ and two valence orbitals $p$ for C atoms with the cutoff radius of 4.0 a.u. (atomic unit). Meanwhile, we used two valence orbitals $s$ and one valence orbital $p$ for H atoms with the cutoff radius of 6.0 a.u. We also used the cutoff energy of 150 Ryd to ensure the convergence of self-consistent calculations. In addition, we set the experimental lattice constant of graphite of 2.46 Å for the length of the periodic unit cell ($x$ axis) and fixed the length of the nonperiodic cell larger than 25 Å to create a sufficient vacuum condition. For the first discussion, we focused on the ribbon width $N=10$. Since the self-consistent calculations using the frozen magnon approach must be carried out for each $\mathbf{q}$, at which the magnetic moment is nearly constant, we set the nonzero $\mathbf{q}$ very close to $\mathbf{q}=0$. This is due to the drastic reduction of the magnetic moment as $\mathbf{q}$ increases, as shown in Fig. 2(a).

After calculating the magnon energy for each $\mathbf{q}$ using Eq. (25), we fit the data by using the quadratic dispersion $E=D\mathbf{q}^{2}$ to obtain the spin stiffness $D$. Figure 2(b) shows the spectrum of magnon energy of 10-ZGNR for the ferromagnetic configuration for the nondoped case. Our calculations yield $D=2966$ meVÅ2 for the ferromagnetic configuration and $D=3583$ meVÅ2 for the antiferromagnetic one for the nondoped case. Note that the antiferromagnetic configuration is more stable than the ferromagnetic one. We see that the spin stiffness of the antiferromagnetic configuration is larger than the ferromagnetic one due to the higher magnon energy of the antiferromagnetic one. This means that the spin stiffness of the most stable configuration in 10-ZGNR should have the highest spin stiffness. We also find that the magnetic moment of the antiferromagnetic configuration is always higher than that of the ferromagnetic one. This means that there is a close relationship between the spin stiffness and the magnetic moment. We will see more clearly this relationship when introducing the doping or varying the ribbon width in the next discussion.

To introduce the hole-electron doping as the uniform background charge, we employ the Fermi level shift method to keep the system neutral. We find the different tendencies for the ferromagnetic and antiferromagnetic configurations for the hole-electron doping. As shown in Fig. 3, for the ferromagnetic configuration, as the hole (electron) doping increases, the spin stiffness increases (reduces). This situation is similar to our last paper [22], which used the fourth-order fit to obtain the spin stiffness. The difference between the quadratic dispersion and the fourth-order dispersion lies in the fitting error. The fitting error of quadratic dispersion is usually higher than that of the fourth-order dispersion. However, the quadratic dispersion, in this case, is sufficient to produce the curve of magnon energy, as shown in Fig. 2(b). For the antiferromagnetic case, on the contrary, the situation is very different. As the hole (electron) doping increases, the spin stiffness reduces (reduces). Therefore, by introducing the hole-electron doping, the ferromagnetic and antiferromagnetic configurations create the antisymmetric and symmetric tendencies, respectively. These tendencies also hold for the magnetic moments, which almost follow the tendencies of spin stiffness. Note that the magnetism in ZGNR disappears if the introduced doping is sufficiently high.

The last discussion is to see the dependence of spin stiffness on the ribbon width. As shown in Fig. 4, we see that the spin stiffness is proportional to the ribbon width. As the ribbon width increases, the spin stiffness increases, see table 1. This means that the required energy becomes higher when the ribbon width increases to make the spin-wave excitations. We also show that all the tendencies for 10-ZGNR also hold for the other ribbon widths when introducing the hole-electron doping. The same tendencies between the spin stiffness and the magnetic moment are due to the mapping of the energy of Heisenberg Hamiltonian onto the total energy difference in the self-consistent calculation. When the exchange constant $J_{ij}$ remains unchanged, the energy of Heisenberg Hamiltonian is proportional to the magnetic moment. At the same situation, the spin stiffness proportionally depends on the total energy difference in the self-consistent calculation. Therefore, the spin stiffness and the magnetic moment should have the same tendency.

4 Conclusions

We obtain the very high spin stiffness for $N$-ZGNR compared to 3$d$ transition metals, as predicted previously by Edwards and Katsnelson [11]. Compared to Yazyev and Katsnelson’s result by the supercell calculation [12], our calculated spin stiffnesses of both the ferromagnetic and antiferromagnetic configurations are higher but still in the same order. We also show that the antiferromagnetic configuration, which is more stable than the ferromagnetic one, has a higher spin stiffness than that of the ferromagnetic one. This implies that the most stable configuration in $N$-ZGNR should have the highest spin stiffness. It is also shown that the spin stiffness increases as the ribbon width increases. It suggests that the distance of two magnetic carbon atoms at the different edges in the unit cell determines the high or low spin stiffness.

When the doping is taken into account, we find the different tendencies of spin stiffness. For the ferromagnetic (antiferromagnetic) configuration, as the hole doping increases, the spin stiffness increases (reduces). Contrarily, the spin stiffness reduces (reduces) as the electron doping increases for the ferromagnetic (antiferromagnetic) configuration. It seems that introducing the doping can control the spin stiffness of $N$-ZGNR.

Acknowledgments

The computations were partly performed using ISSP supercomputers at the University of Tokyo, while the rest of the computations was performed at the Universitas Negeri Jakarta. This work was supported by Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research on Innovative Area, “Nano Spin Conversion Science” (Grant Nos. 15H01015 and 17H05180). It was also supported by a JSPS Grant-in-Aid for Scientific Research on Innovative Area, “Discrete Geometric Analysis for Material Design” (Grant No. 18H04481). It was partially supported by a JSPS Grant-in-Aid on Scientific Research (Grant No. 16K04875)