  Vertically stacked two-dimensional (2D) van der Waals (vdW) heterostructures
based on graphene electrodes represent a promising architecture for
next-generation electronic devices. However, their first-principles
characterizations have been so far mostly limited to the equilibrium state due
to the limitation of the standard non-equilibrium Green's function approach. To
overcome these challenges, we introduce a non-equilibrium first-principles
calculation method based on the recently developed multi-space
constrained-search density functional formalism and apply it to
graphene/few-layer hexagonal boron nitride (hBN)/graphene field-effect
transistors. Our explicit finite-voltage first-principles calculations show
that the previously reported negative differential resistance (NDR)
current-bias voltage characteristics can be produced not only from the
gating-induced mismatch between two graphene Dirac cones but from the
bias-dependent energetic shift of defect levels. Specifically, for a carbon
atom substituted for a nitrogen atom (C$_N$) within inner hBN layers, the
increase of bias voltage is found to induce a self-consistent electron filling
of in-gap C$_N$ states, which leads to changes in voltage drop profiles and
symmetric NDR characteristics. On the other hand, with a C$_N$ placed on outer
interfacial hBN layers, we find that due to the pinning of C$_N$ levels to
nearby graphene states voltage drop profiles become bias-independent and NDR
peaks disappear. Revealing hitherto undiscussed non-equilibrium behaviors of
atomic defect states and their critical impact on device characteristics, our
work points towards future directions for the computational design of 2D vdW
devices
