where $j$ and $k$ runs over the sites of top and bottom layers of graphene. The retarded Green’s function is calculated in energy space with $G^r(E) = ( E I - H_C - \Sigma^r)^{-1}$ where $\Sigma^r$ is the sum of total self energies due to the four leads. Actually, both $H_C$ and $\Sigma^r$ are block diagonal with respect to the layer index since there is no direct electronic coupling. As a result, $G^{>,<,r}$ is also block diagonal. The lead self energies are calculated by standard iterative algorithms of surface Green’s functions.

The polarization functions $\Pi^{>,<,r}$ are calculated according to the time-domain formulas and then Fourier transformed to frequency domain. This appears very fast computationally, but it brings about numerical instability for large systems. $4 \times 4$ system with about 200 atoms is probably the largest system we can obtain reliable result from. We have used 23000 fast Fourier transform points with a spacing $\hbar \Delta \omega =
22$meV. The nonequilibrium information is incorporated through the lead temperatures and chemical potentials by the Keldysh equation, $G^< = G^r \Sigma^< G^a$, at the very beginning when calculating the polarization function $\Pi^{>,<}$, and not through somewhat ad hoc procedure such as Doppler shifting 