A A method to solve the generalized eigenvalue problem

A Numerical Formulation to Calculate the Conductance of Mesoscopic Conductors Using Singular Value Decomposition

\abst

We present a new formulation to calculate the electric conductance of mesoscopic conductors by utilizing the singular value decomposition, which is a mathematical technique to manipulate matrices. Our formulation is useful in treating conductors with rather complicated atomic structures, for which naive recursion formula is cumbersome. It also has an advantage in scaling up the calculation by using parallel computation, which potentially allows us the real-scale calculations at the atomic level. On the other hands, the effects of electron-electron interactions are hard to be treated within this framework, since it depends crucially on the linearity of the system. In this paper, we study graphene nanoribbons with external leads for a simple example. \kwordTransport, Conductance, Landauer Formula, Graphene, Nano-constrictions

1 Introduction

Numerical calculation is one of the most powerful tools to study mesoscopic systems. Especially, the effects of impurities and randomness on the electric conduction have been intensively studied in the light of the Anderson localization [1, 2, 3, 4, 5], and the recursive calculation method of the scattering states of electrons has been developed [2, 3, 6, 7]. These numerical approaches are strengthened by the Landauer formula[8], by which one can reduce the calculation of the conductance to the calculation of the scattering matrix of the conductor [9, 10]. The so-called recursive Green function method has been established as a standard to calculate mesoscopic conductance [11]. This method has been extended to multi-terminal geometries in order to calculate the Hall conductance[5]. Until now many sophisticated numerical algorithms have been developed to be applied to molecular or nano-devices [12, 13, 14, 15, 16, 17, 18].

Recently, many people have strong interests in the transport properties of graphene [19, 20, 21, 22]. The conductance through the graphene nano-wires has been studied, taking into account the effects of impurities, sample edges, and so on [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. The experimental studies are also developing[35, 36]. In addition, there exist various intriguing proposals of the new devices. The effects of the top gate on the transport in graphene have been studied and several new structures are proposed [37, 38, 39, 40, 41]. A method to control the transport phenomena in graphene using strain has been proposed [42, 43, 44, 45] and experimental studies are going on[46, 47]. Together with the nanopatterning of graphene [48], these kinds of technologies are promising for the future application of graphene.

In this paper we introduce a new formulation to calculate numerically the electric conductance of the mesoscopic conductors within the tight-binding approximation. We utilize the so-called singular value decomposition (SVD) in deriving our formalism. The SVD is a mathematical technique similar to the matrix diagonalization. However, by the SVD, we can treat even the non-square matrices. We use the SVD to manipulate the wave-function basis of the conductor. Actually, we can omit, using the SVD, some degrees of freedom, which are not relevant for the transport phenomena, thus reducing the dimensions of the matrices in the calculation. In addition, this procedure is rather easily parallelized, as we will see in Sec. 6. Therefore, it is promising for the future application to the numerical simulation of the mesoscopic systems in the realistic size.

On the other hands, a major price we have to pay for the above advantages is that our formulation is not useful in treating the electron-electron interactions. As one can imagine, once the electron-electron interactions are switched on, the omitted degrees of freedom become relevant to the calculation, thus spoiling the advantages of our formulation. We however consider that the large scale calculation of the transport coefficients is worth pursuing even without the electron-electron interactions, since it gives us the information of the mesoscopic systems, which can be directly compared to the experiments.

The rest of this paper is organized as follows: In Sec. 2, we describe the basic formalism. In Sec. 3, we study the ideal wires, which are used as the model for the external leads. In Sec. 4, we introduce the Landauer formula and derive the equation for the electric conductance. In Sec. 5, we apply our formalism to the graphene nanoribbons with external leads. In Sec. 6, remaining problems and possible future studies are discussed. A prospect of the application of parallel computation to our formulation is also discussed. Supplemental issues are given in Appendices.

2 Basic Formalism

2.1 Singular value decomposition

In developing our formulation, we especially utilized a mathematical technique called “singular value decomposition (SVD)”[49, 50]. Using this, we can decompose an arbitrary complex matrix as

(1)

where and are and unitary matrices, respectively, and , , and are block matrices obtained by partitioning and . ( is the Hermite conjugate of the matrix or vector .) The integer is the matrix rank of , i.e., , and ’s are positive numbers called the singular values of . In the following, we write , , and , where ’s and ’s, respectively, are and dimensional columnar vectors.

The following mathematical features are easily confirmed: We define the range of , , and the null space of , , as

where is the set of complex numbers. Then, is spanned by and is spanned by , as one can see from Eq. (1). It is also useful to see that the left null space of , namely a set of vectors satisfying , is spanned by , which is equivalent to . For convenience, we define an operator, which extract from the SVD of the matrix , i.e., Eq. (1), as,

(2)

The most widely known application of the SVD may be the so-called Moore-Penrose quasi-inverse matrix, which we denote by in this paper. The matrix is given by

(3)

Using this, we can obtain an approximate solution of the equation as . If , this solution is exact (but not unique since we have freedom to add arbitrary superposition of to ). Useful numerical algorithm to calculate SVD is described in Ref. \citenPress:2007:NRE:1403886.

2.2 Schrödinger equation in a recurrence formula

Here we introduce our model of the mesoscopic conductor. The system is described by the tight-binding Hamiltonian,

(4)

where is the annihilation operator of an electron at the -th site and the summation is over all the bonds. We suppress the spin index throughout this paper.

From the Heisenberg form,

(5)

the probability conservation relation is expressed as,

(6)

where is the probability density at the -th site, is the set of sites connected to the -th site and stands for the probability current flowing from the site to . We write the field operator of the electron as

(7)

and, then, is the wave function of the system.

An example of the conductor is depicted in Fig. 1. Each circle corresponds to an atomic site of the tight binding model. The whole system is decomposed into several blocks, numbered by , , and in the figure. We divide the wave function of the -th block into three components, , and , each of which is defined as a columnar vector. The vector () represents the wave function of the sites shown by gray (black) circles in Fig. 1, which are located at the left (right) edge of the -th block and are connected by the bonds to the sites in the -th (-th) block. The vector represents the wave function of the sites in the -th block, which are not included in either or . Here we assume that there is no overlap in the elements of and . If some elements are overlapped, we can extend the area of the -th block until the overlap is dissolved. Such extension is always possible unless we introduce the long-range hopping.

Our aim is to re-express the Schrödinger equation of the whole system into a set of equations, which relate to . In order to carry this out, the vectors, ’s, should be truncated out (or integrated out) by some means. We will show that this process can be readily performed by employing the SVD and the Moore-Penrose quasi-inverse matrices. [49]

Figure 1: The -th, , -th block of the conductor are shown. Gray, black and white circles correspond to , and , respectively.

Let us denote the dimensions of , and by , and , respectively, and, then, the total number of sites in the -th block is given by . The Schrödinger equation within the -th block is given by

(8)

where is the energy, is the identity matrix of order and is a zero vector whose dimension is . Here is the Hamiltonian of sites within the -th block, and () is the hopping matrix elements between the -th and the -th (the -th and the -th) block. Since the Hamiltonian is Hermitian, the relation, , holds. The current flowing from the -th block to the -th block is given by

(9)

where we have introduced the matrix . Since, in the equilibrium, the probability current conserves at all the boundaries between the blocks, the current does not depend on .

Rewriting the Eq. (8), so that only is on the left hand side, yields

(10)

where , and are block matrices composing the matrix , whose dimensions are , and , respectively. The solvability condition of the Eq. (10) with respect to is expressed as

(11)

From this, the relation between and is obtained in the following way: The Eq. (11) means that the vector does not belong to the left null space of . This is expressed, by introducing (see Eq. (2) for the definition), as

(12)

Conversely, if this holds, we obtain the solution for as . (Note that this is not the unique solution.)

Introducing the partitioning of , and (into upper rows, middle rows, and lower rows) as

(13)

the Eq. (12) is rewritten as

(14)

where

(15)

Rearranging these equations, we obtain the relation between and as

(16)

where

(17)

In case of , we should put and omit all the terms with symbol.

2.3 Boundary condition and transfer matrix

Let us write as . When the whole system is composed of blocks, the transport propertied of the system is described by the vectors, . The vector includes and includes . The wave functions and are located out of the blocks and we assume that these give the boundary conditions for the conductor.

Our next goal is to derive the equation which directly relates the components, and . We rewrite whole equations into the form where

(18)

As discussed in the previous section, the solvability condition of with respect to is given by , from which we can derive an equation satisfied by .

We can see that satisfies the condition , where . We introduce the partitioning , where and are block matrices, whose columns correspond respectively to the components, and , of in Eq. (18). Then, we obtain the relation between and as

(19)

Denoting and by and , respectively, and introducing the wave functions of the left lead and the right lead , we obtain the transfer matrices which connect the left and the right leads as

(20)

From this we can calculate the scattering matrix of the conductor.

3 A Generalized Ideal Wire

Here we study the ideal wires as a model for the leads. We assume that the wire is made of a periodic reputation of the set of atoms. In that case, all the matrices () do not depend on and we denote them simply by ().

We assume that the eigenstates of the lead satisfy the relation [11]

(21)

where is a complex number. From the Eq. (16), the equation determining and the corresponding eigenvector is given by

(22)

Here we call the transfer eigenvalue of the wire and the corresponding eignevector, , the transfer eigenvector (or eigenstate). We should note that and are not necessarily square matrices. In that sense, this is a generalized eigenvalue problem. A method to solve this equation is given in the AppendixA.

There are several types of eigenmodes: the mode with corresponds to a propagating mode (or channel), and that with to an evanescent mode. When () the wave function grows (decays) as increases. Some evanescent modes also have imaginary part in and they oscillate as well as grow or decay.

Here we may assume the following mathematical features:

  1. If is a transfer eigenvalue, also is. In case of , .

  2. The current matrix (see Eq. (9)) is diagonalized as

    (23)

    by the transfer eigenstates. (Some note is given in the Appendix B.)

Here we note that 1) seems clear if the system has a time-reversal symmetry. However, similar feature may also exist even without time-reversal symmetry.

In the preceding works, it has sometimes been assumed that the direction parallel and perpendicular to the wire are separable, such as the case of a simple wire made of a square lattice. In that case, we can define the transverse modes rather easily and the property expressed in 2) is clearly satisfied. Our statement is the generalization of such a case. It seems that 2) is quite reasonable from the physical point of view even in generalized situations, although we do not go into details in this paper.

4 Landauer Formula and Conductance

We consider the situation depicted in Fig. 2. The scatterer is connected to the left and the right leads. Two leads are not necessarily identical to each other. However each of them should be have a completely periodic structure, so that the channels are well-defined. We assume that the two leads are semi-infinitely long.

Figure 2: Geometry of the conductor is schematically shown. The shaded area is the scatterer connected to the left and the right leads. Vertical dashed lines are boundaries of the blocks. Symbols with “ev.”  indicate the evanescent modes. The numbers of each modes are also indicated.

As depicted in Fig.2, we consider the right- and the left-moving modes and the evanescent modes decaying with distance from the scatterer. We assume that there are () channels and () evanescent modes in the left (right) lead. (Note that one channel corresponds to a pair of the right- and the left-moving mode.)

Now we calculate the conductance of this conductor. To do this, we utilize the so-called Landauer formula. First we introduce matrix , which relates the probability currents of the incoming waves to that of the outgoing waves, as

(24)

where , etc. represent the probability currents carried by the propagating modes; “ ”  and “”  indicate respectively the left and the right lead, and “”  and “”  the right-moving and the left-moving mode.

We denote the wave functions at the left and the right lead by and , respectively. Each function is given as a superposition of the right-moving, the left-moving and the evanescent modes as

(25)

where and are constants. The vectors and represent the eigenvectors of the right-moving and the left-moving modes, respectively, in the left (right) lead. The vector represents the wave function of the evanescent modes in the left (right) lead.

Here we introduce the matrix representation of a basis of the wave function as . Then the Eqs. (25) are rewritten as

(26)

where , etc. are columnar vectors. From the Eq. (20), we obtain

(27)

This can be rewritten in a matrix form as

(28)

where

(29)

Here we introduce (see Eq. (2)). Then, the solvability condition of the Eq. (28) with respect to and is given by

(30)

which relates the amplitudes of the incoming and the outgoing waves. Note that, if there are no evanescent modes, should be the identity matrix of order .

Writing , where is the left (right) columns of , and introducing , , Eq. (30) becomes . Using the Moore-Penrose inverse of and defining , we reach the expression

(31)

From the argument of the Sec. 3 and the Eq. (23), the total probability current is given by the sum of the independent contributions of each channel. Let us denote the contribution of the -th “in (out)”  channel by , the -th components of and are given by

Using Eq. (31), the probability current carried by the -th outgoing channel is given by

(32)

Here we separated the summation into the diagonal and the off-diagonal terms with respect to and . We assume that electrons are injected to the modes from the reservoir incoherently, and then the cross term (the last line of Eq. (32)) vanishes after time averaging. As a result we obtain the following relation,

(33)

Then the elements of the scattering matrix is given by

(34)

Using the Landauer formula we can calculate the conductance as [10]

(35)

In actual calculations, the wave functions should be normalized appropriately. However, the normalization is not relevant to the results, since the present formulation depends only on the ratio of the incoming and the outgoing current. Therefore we can adopt any normalization for our convenience.

5 Application to Graphene Constrictions

We apply the present formulation to the graphene constrictions of the several forms. The character of the graphene is taken into account by assuming the honeycomb lattice structure. The hopping is restricted to the nearest neighbors for simplicity.

Figure 3: The geometries of the graphene nanoribbons studied in this paper are shown. I, II and II have constrictions in the middle of the conductor, though the lengths of the constricted regions differ. The constricted region is armchair nanoribbon in all cases as one can see from the magnification. Both leads are modeled by square lattices.
Figure 4: The energy dependence of the conductance of the wires depicted in Fig. 3. (a), (b) and (c) correspond to I, II, and III of Fig. 3, respectively. The vertical axis is the conductance including spin degree of freedom. The step-like function, shown in three figures, is the conductance of an ideal armchair wire. (d) is the magnification near the charge neutrality point . The peaks are those of I and II, as shown in the figure. III also shows a small peak at , however it is non-distinctive.

We treat the graphene wires with a short (I), a middle (II), and a long (III) constriction as depicted in Fig. 3. The narrowest parts of the wires are composed of armchair nanoribbons. The left and the right leads are assumed to be made of two-dimensional square lattices laid parallel to the graphene surface, whose lattice spacing is (: the minimum C-C spacing in the graphene). The directions of the square lattices are intentionally rotated from the wire axis. The transfer integrals within the leads () and those between the graphene and the leads () are set as and . Only the nearest neighbor hopping is assumed within the leads, and the hopping between the leads and the graphene is assumed only between the sites, whose horizontal separation (parallel to the graphene surface) is smaller than .

The conductance of the wires at zero temperature are calculated as a function of the energy and the results are shown in Fig. 4 (a) (c). The apparent step like function (indicated in blue) is the conductance of an ideal armchair wire. As one can see, the conductance curves of the constrictions consist of spiky peaks, which may arise from the various resonances of the conduction electrons. The resemblance to the blue steps are not obvious in all cases. However at some energies the conductance approaches the perfect transmission, especially in the lowest step .

Near zero energy, some peaks are seen in (a) and (b) (The magnification is shown in (d)). These may come from the transmission through the evanescent modes decaying in the constricted region. This explanation is plausible since such peaks are not significant in the case of the longest constriction (c), though a tiny peak still exists.

The conductance at a finite temperature is calculated from

(36)

where is the Fermi distribution function and is the conductance at the bias and the temperature . In Fig. 5, we have shown the results of the constriction I for and . At , we can see several step-like structures, however their relation to the steps of the ideal wire is not clear at present. More intensive study is required to clarify the conductance quantization in graphene constrictions.

Figure 5: Conductance of Fig.3 (a) after the thermal average. The thin green curve, thick red curve and dashed purple curve show the result of , and , respectively.

We note on the particle-hole symmetry of these systems. It has been known that the graphene nanoribbons show particle-hole symmetry if the hopping is limited to the nearest neighbors (NN’s). The same property holds for the square lattices. With the particle-hole symmetry, the conductance is symmetric with respect to the energy inversion . The particle-hole symmetry in the graphene is broken when the next nearest neighbor (NNN) hopping is switched on. In the present calculation, the graphene and the leads separately hold particle-hole symmetry. However, the coupling of the graphene to the leads effectively induces the NNN hopping in the graphene, which may probably break the particle-hole symmetry. The slightly asymmetric behaviors, as seen in Figs. 4 and 5, may originate from such a mechanism.

6 Discussion

In this paper we have introduced a new formulation to calculate the conductance of the mesoscopic conductors. The key point of our formulation is to omit, using the SVD, the degrees of freedom in the scatterer irrelevant for the transport phenomena. This procedure may be available for other scattering problems, although, until now, the SVD does not seem to be used so often in the studies of physics.

In the present study we have not treated the random potential, the magnetic field, the spin orbit coupling, and the lattice distortion. Inclusion of these into our formalism may be straight forward. They will provide useful information in the various fields, such as spintronics or topological insulators[51].

It is also important to treat the multi-terminal geometry[5]. For example, the Hall conductance belongs to this category. We consider that the present treatment can be extended to such situations by introducing appropriate partitioning of the sample, which, however, needs some more theoretical efforts.

Finally, we point out a possibility of applying parallel computation to our formulation. Here we show that the dimension of the large matrix appearing in the Eq. (18) can be reduced by the following method. Let us take a part of the equation,

(37)
(38)

We rewrite this by introducing block matrices as,

(39)

Then, we obtain the solvability condition for as

(40)

Rewiring the blocks of the matrix by

(41)

we result in

(42)

Now is deleted and the matrix dimension is reduced.

By using the above transformation we can reduce the number of blocks in the large matrix in the Eq. (18) by one. Applying this method repeatedly, the matrix dimension is remarkably reduced. Furthermore, this process can be carried out for several different ’s simultaneously, which allows us to speed up the calculation by using parallel computation.

Using the above procedure, we can enlarge the system size drastically, which potentially allows us the real-scale calculations at the atomic level. As we have pointed out before, our formulation is not useful for the systems with the electron-electron interaction. However we consider that the advantages of our formulation overcome such a disadvantage.

7 Summary

We have proposed a new formulation to calculate the electric conductance of the mesoscopic conductors, which has a potential advantage in treating huge systems using parallel computation. A simple example of the calculation is shown for graphene nanoribbons with external leads.

Acknowledgement

The author is grateful to H. Yoshioka, A. Kanda and H. Tomori for valuable discussions. Especially, he owes the argument on the particle-hole symmetry in Sec. 5 largely to H. Yoshioka.

This work was supported by JSPS KAKENHI Grant Numbers 24540392, 22540329.

Appendix A A method to solve the generalized eigenvalue problem

We present a simple method to solve the equation of the following type,

(43)

where and are matrices ( in general). A number, , and a vector, , should be determined. In solving Eq. (43), it is advantageous to formulate the problem on the basis of the eigenvalue problem, although some attention is needed. The calculation consists of the following two steps.

a.1 Removal of the intersection of the null spaces of and

First we note the case where the intersection of the null spaces of the matrices and is not empty. Let us suppose that a vector is in the null spaces of and simultaneously, namely, satisfies and . Then, is a solution of Eq. (43) irrespective of the value . In such case, we cannot determine the eigenvectors other than from Eq. (43), since even if and are square matrices, the characteristic polynomial vanishes for any and, then, the eigenvalues are not determined. In order to avoid this situation we need to remove the intersection of the null spaces of and from the basis of the vectors. This process is executed as follows (see Chap. 12.4 of Ref. \citenGolub:2012wp).

First we consider a block matrix made from and and introduce its SVD as follows ,

(44)

where ’s are the singular values of , and . Here the columns of form a complete orthonormal basis of the -dimensional space, and the columns of form the basis of the null space of , which is actually the intersection of the null spaces of and . Therefore if we restrict the basis of   to , the intersection of the null spaces of and is removed from the solution space of .

Here we set

(45)

If we divide into upper rows and lower rows and put

(46)

the relations and hold and the Eq. (43) is rewritten as

(47)

Note that . At this stage, and are not necessarily square matrices.

a.2 Making matrices and square

Next we make the matrices and square by truncating some of the null spaces of or . We put the dimensions of and to be and also define [52]. The SVD of and are given as