**5. Numerical Results**

In this section we present numerical results obtained with the BITLLES simulator (see Section 4) that demonstrate the ability of the Bohmian conditional wavefunction approach to provide dynamics information for both Markovian and non-Markovian scenarios. We simulate a two-terminal electron device whose active region is a graphene sheet contacted to the outer by two (ohmic) contacts. Graphene is a 2*D* material that has attracted a lot of attention recently because of its high electron mobility. It is a gapless material with linear energy band, which differs from the parabolic energy bands of traditional semiconductors. In graphene, the conduction and valence bands coincide at an energy point known as the Dirac point. Thus, the dynamics of electrons is no longer governed by an (effective mass) Schrödinger equation but by the Dirac equation, allowing transport from the valence to the conduction band (and vice versa) through Klein tunneling. A Bohmian conditional bispinor (instead of a conditional scalar wavefunction) is used to describe electrons inside the device. The change from a wavefunction to a bispinor does not imply any conceptual difficulty but just a mere increment of the computational cost. More details can be found in Appendix C.

In particular, we want to simulate electron transport in graphene at very high frequencies (THz) taking into account the electromagnetic environment of the electron device. Typically, nanoscale devices are small enough to assume that, even at THz frequencies, the electric field is much more relevant than the magnetic field. Therefore, only the Gauss law (first Maxwell's equations) is enforced to be fulfilled in a self-consistent way (i.e., taking into account the actual charge distribution in the active region). However, the environment of nanoscale devices is commonly a metallic element of macroscopic dimensions. In there, the magnetic and electric fields become both relevant, acting as active (detecting or emitting) THz antennas. For the typical electromagnetic modes propagating in the metals, the magnetic and electric fields are translated into the language of currents and voltages and the whole antenna is modeled as a part of an electric circuit. In this work, the graphene device interacts with an environment that is modeled by a Resistor (R) and a capacitor (C) connected in series through ideal cables (see the schematic plots in Figure 5a–c).

The active region of the graphene device is simulated with the Bohmian conditional wavefunction approach explained in the previous section, while the RC circuit is simulated using a time-dependent finite-difference method. We consider the system plus environment to be in equilibrium. Specifically, the self-consistent procedure to ge<sup>t</sup> the current is as follows: an initial (at time *t* = 0) zero voltage is applied at the source (*VS*(0) = 0) and drain (*VD*(0) = 0) contacts of the graphene active region. At room temperature this situation yields a non-zero current from Equation (20) (i.e., *I*(0) = 0) because of thermal noise. Such current *I*(0) enters the RC circuit and leads to a new voltage *VS*(*dt*) = 0 at the next time step *dt* (where *dt* represents the time step that defines the interaction between the RC circuit and the quantum device which was set to *dt* = 0.5fs). The new source *VS*(*dt*) = 0 and fixed drain *VD*(*dt*) = 0 voltages now lead to a new value of the current *I*(*dt*) = 0 in (20) which is different from zero not only because of thermal noise but also because there is now a net bias ( *VD*(*dt*) − *VS*(*dt*) = 0). This new current *I*(*dt*) is used (in the RC circuit) to ge<sup>t</sup> a new *VS*(2*dt*) that is introduced back in the device to obtain *I*(2*dt*) and so on. Importantly, as the system and environment are in equilibrium, the expectation value of *I*(*t*) is zero at any time, i.e., *I*(*t*) = 0 ∀*t*.

We consider three different environments (with different values of the capacitance). In Figure 5a we plot the total (particle plus displacement) electrical current at the end of the active region when *R* = 0 and *C* = ∞. The same information is shown in Figure 5b,c for two different values of the capacitance *C* = 2.6 × 10−<sup>17</sup> F and *C* = 1.3 × 10−<sup>17</sup> F. In all cases the value of the resistance is *R* = 187 Ω, and we assumed the current *I*(*t*) to be positive when it goes from drain to source.

The effect of the RC circuit is, mainly, to attenuate the current fluctuations, which are originated due to thermal noise. This can be seen by comparing Figure 5a with Figure 5b,c. The smaller the capacitance the smaller the current fluctuations. This can be explained as follows: when the net current is positive, the capacitor in the source starts to be charged and so the voltage at the source increases trying to counteract the initially positive current. Therefore, the smaller the capacitance the faster the RC circuit reacts to a charge imbalance.

**Figure 5.** Total (particle plus displacement) electrical current *I*(*t*) evaluated at the ammeter as a function of time for a graphene device connected to three different RC circuits with *R* = 187 Ω. The values of the capacitances are: (**a**) *C* = <sup>∞</sup>, (**b**) *C* = 2.6 × 10−<sup>17</sup> F and (**c**) *C* = 1.3 × 10−<sup>17</sup> F.

In Figure 6 we plot the total (particle plus displacement) current–current correlations as a function of the observation time *τ* for the three scenarios in Figure 5. Correlations at very small observation times provide information of the variance of the current, which, as explained above, is reduced as the value of the capacitance is increased. Numerical simulations (not shown here) exhibit that the role of the resistor *R* is less evident because the active region itself has a much larger (than *R* = 187 Ω) associated resistance. Numerically the distinction between Markovian and non-Markovian dynamics boils down to the comparison of time correlations as defined in Equations (9) and (13). Since there is no net bias applied to the graphene device (i.e., it is in equilibrium), an ensemble average of the current (over an infinite set of trajectories like the one depicted in Figure 5) yields *I*(*t*) = 0 ∀*t*. Time correlation functions computed in Equation (9) are thus zero by construction, i.e., *I*(*t*)*I*(*t* + *τ*) = 0 ∀*t*, *τ*. Therefore, the non-Markovian dynamics occurring at very high-frequencies (below the ps time-scale in Figure 6 expressly shows) fixes the correlation time of the environment at *tD* ∼ ps. Although all three

values of the capacitance *C* in Figure 6 yield the same order of magnitude for *tD* ∼ ps, it seems also true that the smaller the value of the capacitance, the smaller *tD*.

Current–current correlations shown in Figure 6 can be better understood by assessing the transit time of electrons. For a velocity of roughly 10<sup>6</sup> m/s inside an active region of *L* = 40 nm length is roughly *τT* = *L*/*vx* = 0.04 ps. Positive correlations correspond to transmitted electrons traveling from drain to source (as well as electrons traversing the device from source to drain). While 0 < *t* < *τT* electrons are transiting inside the active region, such electrons provide always a positive (or negative) current as seen in expression (20). In other words, if we have a positive current at time *t* because electrons are traveling from drain to source, we can expect also a positive current at times *t* satisfying *t* < *t* < *t* + *τT*. The negative correlations belong to electrons that are being reflected. They enter in the active region with a positive (negative) velocity and, after some time *τR* inside the device, they are reflected and have negative (positive) velocities until they leave the device after spending roughly 2*τR* in the active region. Thus, during the time *τR* < *t* < 2*τR* which will be different for each electron depending on the time when they are reflected, we can expect negative correlations. Interestingly, during the 4 ps simulation the number of Bohmian trajectories reflected are double in the black (*C* = ∞) simulation than in the red one (*C* = 1.3 × 10<sup>17</sup> F). This can be explained in a similar way as we explained the reduction of the current fluctuations. The fluctuations of the electrical current imply also fluctuations of the charge inside the active region, which are translated (through the Gauss law) into fluctuations of the potential profile. Thus, the larger the noisy current, the larger the noisy internal potential profile. This implies a larger probability of being reflected by the Klein tunneling phenomenon. Therefore, if one aims at describing the dynamics of nanoscale devices with a time-resolution *τ* that is comparable to (or goes beyond) the electron transit time *τT*, a non-Markovian approach is necessary. This is so because the total current *I*(*t*) (which has contributions from the displacement and the particle currents) shows correlations at times that are smaller than the electron transit time.

**Figure 6.** Total current–current correlation as a function of time for the three different experiments in Figure 5. The zero is indicated by a dashed line to show the tendency of the total current, understood as a property of the environment, to vanish at long times *τ*. Zero autocorrrelation implies an independence between *I*(*t*) and *I*(*t* + *τ*) which is typical for Markovian scenarios. This is not true for the short *τ* considered here which are the representatives of the non-Markovian dynamics.
