H =- \frac{1}{A} \int_{0}^{+\infty} \frac{d\omega}{2\pi} \hbar \omega {\rm Tr} (D^> \Pi^<_1 - D^< \Pi^>_1 ). \end{aligned}$$ Here $D^{>,<}$ is the greater/lesser Green’s function for the scalar photon, which is calculated from the Keldysh equation, $D^{<,>} = D^r \Pi^{<,>} D^a$, and retarded Green’s function is obtained by solving the Dyson equation, $D^r = v + v \Pi^r D^r$ ($v$ is the bare Coulomb potential with matrix element $1/(4\pi\epsilon_0 r_{ij})$ between the tight-binding sites $i$ and $j$ of a distance $r_{ij}$). $\Pi^{<,>} = \Pi^{<,>}_1 + \Pi^{<,>}_2$ is block-diagonal and is obtained with the random phase approximation (RPA). $A$ is the area of the graphene. Further computational details will be presented in Appendix \[app-d\].

Comparing with the calculation based on FE, the NEGF method provides a rigorous way to extend the drift effects into nanoscale ballistic systems without any phenomenological assumptions. Due to the smallness of the sample, the transport is ballistic with the electric current at the bottom driven layer given by $I_1 = (4 e/h) (\mu_{1R} - \mu_{1L})$, independent of the sample width or length. We do not observe Coulomb drag effect [@narozhny2016coulomb], as the electric current in the top unbiased layer is very close to 0, while the thermal current 