**6. Conclusions and Final Remarks**

Theoretical approaches to open quantum systems that rely on the manipulation of state vectors instead of a reduced density matrix have well known computational advantages. Two major benefits are the substantial reduction of the dimensionality of the involved mathematical objects and the preservation of complete positivity [18]. However, substituting density matrices by state vectors constitutes also an attempt to achieve a more detailed description of the dynamics of open quantum systems [6,19]. It is well recognized, for example, that the continuous measurement of an open quantum system with associated Markovian dynamics can be described by means of a SSE (see Table 2 O4). The conditional state solution to such an equation over some time interval can be linked to a "quantum trajectory" [12,19] of one property of the environment. Thus, the conditional state can be interpreted as the state of the open system evolving while its environment is under continuous monitoring. This is true in general for Markovian systems, no matter whether or not the environment is being actually measured (i.e., it is valid for both Figure 1a,b). This fact is of grea<sup>t</sup> importance for designing and experimentally implementing feedback control in open quantum systems [35]. If this interpretation could also be applied to non-Markovian SSEs [33,37], then this would be very significant for quantum technologies, especially in condensed matter environments (e.g., electron devices), which are typically non-Markovian [6].

**Table 2.** Validity of Bohmian vs. orthodox conditional states to provide dynamic information of open quantum system depending on the relation between the environment decoherence time *tD* and the observation period *τ*. Here (un)measured refers to unmeasured and measured indistinctively.


Unfortunately, for non-Markovian conditions, the above interpretation is only possible for the rather exotic scenario where the environment is being continuously monitored and the system is strongly coupled to it. As no correlation between the system and the environment can build up, the evolved system is kept in a pure state. This is the well-known quantum Zeno regime [65,66], under which conditional states can be trivially used to describe the frozen properties of the system (see Table 2 O1). Without the explicit consideration of the measurement process (as in Figure 1a), however, the postulates of the orthodox theory restrict the amount of dynamical information that can be extracted from state vectors (see Table 2 O2). In most general conditions, for *τ* > 0 and non-Markovian dynamics, while conditional states can be used to reconstruct the reduced density matrix, they cannot be used to evaluate time-correlations (see Table 2 O3) [20,23]. This is not only true when the environment is being measured (as in Figure 1b), but also when it is not measured (as in Figure 1a).

Therefore, we turned to a nonorthodox approach: the Bohmian interpretation of quantum mechanics. The basic element of the Bohmian theory (as in other quantum theories without observers) is that the intrinsic properties of quantum systems do not depend on whether the system is being measured or not. Such ontological change is, nevertheless, fully compatible with the predictions of orthodox quantum mechanics because a measurement-independent reality of quantum objects is not in contradiction with non-local and contextual quantum phenomena. Yet, the ontological nature of the trajectories in Bohmian mechanics introduces the possibility of evaluating dynamic properties in terms of conditional wavefunctions for Markovian and non-Markovian dynamics, no matter whether the environment is being actually measured or not (see Table 2, B1–B4 and Figure 7a,b).

In summary, the Bohmian conditional states lend themselves as a rigorous theoretical tool to evaluate static and dynamic properties of open quantum systems in terms of state vectors without the need of reconstructing a reduced density matrix. Formally, the price to be paid is that for developing a SSE-like approach based on Bohmian mechanics one needs to evaluate both the trajectories of the environment and of the system see (d) Table 1. Nonetheless, we have seen that this additional computational cost can be substantially reduced in practical situations. For THz electron devices (see Section 5), for example, we showed that invoking current and charge conservation one can easily ge<sup>t</sup> rid of the evaluation of the environment trajectories. This reduces substantially the computational cost associated to the Bohmian conditional wave function approach (as shown (e) in Table 1). Let us also notice that here we have always assumed that the positions of the environment are the variables that the states of the system are conditioned to. However, it can be shown that the mathematical equivalence of the SSEs with state vectors conditioned to other "beables" of the environment (different from the positions) is also possible. It requires using a generalized modal interpretation of quantum phenomena, instead of the Bohmian theory. A review on the modal interpretation can be found in [67,68].

**Figure 7.** (**a**) Figure depicting the Unmeasured Bohmian approach in which the computation of any property (electric current in an electron device) is independent of the measuring apparatus. (**b**) The continuous measurement of the electric current through an ammeter (measuring apparatus) can be also described in Bohmian mechanics by including the degrees of freedom of the measuring apparatus.

As an example of the practical utility of the Bohmian conditional states, we have introduced a time-dependent quantum Monte Carlo algorithm, called BITLLES, to describe electron transport in open quantum systems. We have simulated a graphene electron device coupled to an RC circuit and computed its current–current correlations up to the THz regime where non-Markovian effects are relevant. The resulting simulation technique allows to describe not only DC and AC device's characteristics but also noise and fluctuations. Therefore, BITLLES extends to the quantum regime the computational capabilities that the Monte Carlo solution of the Boltzmann transport equation has been offering for decades for semi-classical devices.

**Author Contributions:** Conceptualization, D.P., G.A. and X.O. ; methodology, D.P., G.A. and X.O.; software, E.C. and X.O.; validation, D.P., G.A., E.C.and X.O.; formal analysis, D.P., G.A. and X.O.; investigation, D.P., G.A. and X.O.; resources, D.P., G.A. and X.O.; data curation, D.P., E.C., G.A., and X.O.; writing–original draft preparation, D.P., G.A. and X.O.; writing–review and editing, D.P., E.C., G.A. and X.O.; visualization, D.P., E.C., G.A. and X.O.; supervision, G.A. and X.O.; project administration, G.A. and X.O.; funding acquisition, X.O and G.A.

**Funding:** We acknowledge financial support from Spain's Ministerio de Ciencia, Innovación y Universidades under Grant No. RTI2018-097876-B-C21 (MCIU/AEI/FEDER, UE), the European Union's Horizon 2020 research and innovation programme under gran<sup>t</sup> agreemen<sup>t</sup> No Graphene Core2 785219 and under the Marie Skodowska-Curie gran<sup>t</sup> agreemen<sup>t</sup> No 765426 (TeraApps). G.A. also acknowledges financial support from the European Unions Horizon 2020 research and innovation programme under the Marie Skodowska-Curie Grant Agreement No. 752822, the Spanish Ministerio de Economa y Competitividad (Project No. CTQ2016-76423-P), and the Generalitat de Catalunya (Project No. 2017 SGR 348).

**Acknowledgments:** The authors are grateful to Howard Wiseman and Bassano Vacchini for enlightening comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A. D'Espagnat Distinction between "Proper" and "Improper" Mixtures**

An alternative explanation on the difficulties of state vectors in describing open quantum system comes from the distinction between "proper" and "improper" mixtures by D'Espagnat.


D'Espagnat claims that the ignorance interpretation of the "proper" mixture cannot be given to the "improper" mixture. The D'Espagnat's argumen<sup>t</sup> is as follows. Let us assume a pure global system (inclding the open system and the environment) described by Equation (1). Then, if we accept that the physical state of the system is given by |*ψq*(*t*), we have to accept that the system-plus-environment is in the physical state |*q*⊗|*ψq*(*t*) with probability *<sup>P</sup>*(*q*). The ignorance interpretation will then erroneously conclude that the the global system is in a mixed state, not in a pure state as assumed in Equation (1). The error is assuming that the system is in a well-defined state that we (the humans) ignore it. This is simply not true. D'Espagnat results shows that a conditional state cannot be a description of an open system with all the static and dynamic information that we can ge<sup>t</sup> from the open system.

In addition, let us notice that the conclusion of D'Espagnat applies to any open quantum system without distinguishing between Markovian or non-Markovian scenarios. However, indeed, there is no contradiction between the D'Espagnat conclusion and the attempt of the SSE of using pure states to describe Markovian open quantum systems for static and dynamic properties. Both are right. D'Espagnat discussion is a formal (fundamental) discussion about conditional states, while the discussion about Markovian scenarios is a practical discussion about simplifying approximation when extracting information of the system at large *τ*.

Finally, let us notice that the D'Espagnat conclusions do not apply to Bohmian mechanics because the Bohmian definition of a quantum system involves a wave function plus trajectories. The conditional state *ψ*˜ *<sup>Q</sup><sup>i</sup>*(*t*)(*<sup>x</sup>*, *t*) = <sup>Ψ</sup>(*<sup>x</sup>*, *Qi* (*t*), *t*), together with the environment and system trajectory *Qi* (*t*) and *X<sup>i</sup>* (*t*), contains all the (static and dynamic) information of the open system in this *i*-th experiment. An ensemble over all experiments prepared with the same global wavefunction <sup>Ψ</sup>(*<sup>x</sup>*, *q*, *t*) requires an ensemble of different environment and system trajectories *Qi* (*t*) and *X<sup>i</sup>* (*t*) for *i* = 1, 2, ..., *M* with *M* → ∞.

### **Appendix B. Orthodox and Bohmian Reduced Density Matrices**

The orthodox and Bohmian definitions of a quantum state are different. The first uses only a wave function, while the second uses the same wave function plus trajectories. It is well known that both reproduce the same ensemble values by construction. Here, we want to discuss how the orthodox density matrix (without trajectories) can be described by the Bohmian theory with trajectories.

We consider a system plus environment defined by a Hilbert space H that can be decomposed as H = H*x* ⊗ H*q* where *x* is the collective positions of particles of the system while *q* are the collective positions of the particles of the environment.The expectation value of any observable *Ox* of the system can be computed as *Ox* = Ψ|*O*<sup>ˆ</sup> *x* ⊗ ˆ*Iq*|Ψ with ˆ*Iq* the identity operator for the environment. We describe the typical orthodox procedure to define the orthodox reduced density matrix by tracing out all degrees of freedom of the environment:

$$
\rho(\mathbf{x}, \mathbf{x}', t) = \int dq \Psi^\*(\mathbf{x}', q, t) \Psi(\mathbf{x}, q, t) \tag{A1}
$$

From Equation (A1) the mean value of the observable *Ox* can be computed as,

$$
\langle O\_{\mathbf{x}} \rangle = \int d\mathbf{x} \left( O\_{\mathbf{x}} \rho(\mathbf{x}, \mathbf{x}', t) |\_{\mathbf{x}' = \mathbf{x}} \right) \tag{A2}
$$

In this appendix, we want to describe Equations (A1) and (A2) in terms of the Bohmian conditional wavefunctions and trajectories described in the text. The conditonal wavefunction associated to the system during the *i*-th experiment conditioned on a particular value of the environment *Qi* (*t*) is defined as *ψQ<sup>i</sup>*(*t*)(*<sup>x</sup>*, *t*) = <sup>Ψ</sup>(*<sup>x</sup>*, *Qi* (*t*), *t*), being <sup>Ψ</sup>(*<sup>x</sup>*, *q*, *t*) = *<sup>x</sup>*, *q*|Ψ the position representation of the global state. We start from the general expression for the ensemble value in the position representation as,

$$
\langle O\_{\mathbf{x}} \rangle = \int d\mathbf{x} \int d\boldsymbol{q} \,\,\Psi^\*(\mathbf{x}, \boldsymbol{q}, t) O\_{\mathbf{x}} \Psi(\mathbf{x}, \boldsymbol{q}, t) \tag{A3}
$$

Multiplying and dividing by <sup>Ψ</sup><sup>∗</sup>(*<sup>x</sup>*, *q*, *t*) we get,

$$\begin{split} \langle \mathcal{O}\_{\mathbf{x}} \rangle &= \quad \int d\mathbf{x} \int dq \, \frac{|\Psi(\mathbf{x}, q, t)|^{2} O\_{\mathbf{x}} \Psi(\mathbf{x}, q, t)}{\Psi^{\*}(\mathbf{x}, q, t)} \\ &= \quad \frac{1}{M} \sum\_{i=1}^{M} \int d\mathbf{x} \int dq \frac{\delta[(\mathbf{x} - \mathbf{X}^{i}(t)) \delta[(q - \mathbf{Q}^{i}(t)) \mathcal{O}\_{\mathbf{x}} \Psi(\mathbf{x}, q, t)}{\Psi^{\*}(\mathbf{x}, q, t)}] \\ &= \quad \frac{1}{M} \sum\_{i=1}^{M} \int d\mathbf{x} \frac{\mathcal{O}\_{\mathbf{x}} \Psi(\mathbf{x}, \mathcal{Q}^{i}(t), t)}{\Psi^{\*}(\mathbf{x}, \mathcal{Q}^{i}(t), t)} \delta(\mathbf{x} - \mathcal{X}^{i}(t)) \end{split} \tag{A4}$$

where we have used the quantum equilibrium condition |Ψ(*<sup>x</sup>*, *q*, *t*)|<sup>2</sup> = 1*M* ∑*Mi*=<sup>1</sup> *<sup>δ</sup>*[(*x* − *<sup>X</sup><sup>i</sup>*(*t*)]*δ*[(*q* − *Q<sup>i</sup>*(*t*)] with *M* → ∞. Now, we multiply and divide by <sup>Ψ</sup><sup>∗</sup>(*<sup>x</sup>*, *Q<sup>i</sup>*(*t*), *t*) to get,

$$\begin{split} \langle O\_x \rangle &= \frac{1}{M} \sum\_{i=1}^{M} \int dx \frac{O\_x \Psi^\*(x, Q^i(t), t) \Psi(x, Q^i(t), t)}{|\Psi(x, Q^i(t), t)|^2} \delta(x - X^i(t)) \\ &= \int dx \left[ O\_x \sum\_{i=1}^{M} P\_i \tilde{\Psi}^{i\*}(x', t) \tilde{\Psi}^i(x, t) \right]\_{x' = x = X^i(t)} \end{split} \tag{A5}$$

where *Pi* = 1/*M* can be interpreted as the probability associated to each *i* = 1, 2, ..., *M* experiment and we have defined:

$$\Psi^i(\mathbf{x}, t) \equiv \frac{\Psi(\mathbf{x}, Q^i(t), t)}{\Psi(X^i(t), Q^i(t), t)}. \tag{A6}$$

Now, once we arrive at Equation (A5), one can be tempted to define a type of Bohmian reduced density matrix in terms of the conditional wavefunctions for *i* = 1, 2, ..., *M* experiments as,

$$\rho\_{\mathbf{w}}(\mathbf{x}^{\prime},\mathbf{x},t) = \sum\_{i=1}^{M} P\_{i} \,\tilde{\Psi}^{i\*}(\mathbf{x}^{\prime},t) \tilde{\Psi}^{i}(\mathbf{x},t) = \sum\_{i=1}^{M} P\_{i} \,\frac{\Psi^{\*}(\mathbf{x}^{\prime},Q^{i}(t),t)}{\Psi^{\*}(X^{i}(t),Q^{i}(t),t)} \frac{\Psi(\mathbf{x},Q^{i}(t),t)}{\Psi(X^{i}(t),Q^{i}(t),t)},\tag{A7}$$

where we have arbitrarily eliminated the role of the trajectories. However, strictly speaking Equation (A1) is not equal to Equation (A7). If we include all *i* = 1, 2, ..., *M* experiments in the computation of (A7), there are trajectories *Q<sup>i</sup>*(*t*) and *Qj*(*t*) that at the particular time *t* can be represented by the same conditional wavefunction *ψ*˜*<sup>i</sup>*(*<sup>x</sup>*, *t*) = *ψ*˜*j*(*<sup>x</sup>*, *t*) if *Q<sup>i</sup>*(*t*) = *Qj*(*t*). Such over-summation due to the repetition of the same trajectories is not present in (A1).

To simplify the subsequent discussion, let us assume that *q* is one degree of freedom in a 1D space. Let us cut such 1D space into small intervals of length Δ*q*. Each interval is defined as *j* Δ*q* < *q* < (*j* + 1) Δ*q* and it is labelled by the index *j*. Then, we can define *Gj*(*t*) as the number of positions *Q<sup>i</sup>*(*t*) that are inside the *j*-interval at time *t* as:

$$G^{j}(t) = \sum\_{i=1}^{M} \int\_{j\Delta q}^{(j+1)\Delta q} \delta[q - Q^{i}(t)] dq \tag{A8}$$

With this definition, assuming that Δ*q* is so small that all *Q<sup>i</sup>*(*t*) inside the interval and all the corresponding Bohmian conditional wave functions <sup>Ψ</sup>(*<sup>x</sup>*, *Qj*(*t*), *t*), system positions *<sup>X</sup><sup>i</sup>*(*t*), and probabilities *Pi* are almost equivalent, and given by *qj*, <sup>Ψ</sup>(*<sup>x</sup>*, *qj*, *t*), *xj* and *Pj* respectively, we can change the sum over *i* = 1, ..., *M* experiments into a sum over *j* = ..., −1, 0, 1, ... spatial intervals to rewrite Equation (A7) as:

$$\rho\_{\mathbf{w}}(\mathbf{x}',\mathbf{x},t) = \sum\_{j=-\infty}^{j=+\infty} \mathbf{G}^j(\mathbf{t}) \mathbf{P}\_j \frac{\mathbf{F}^\*(\mathbf{x}',q^j,t)}{\mathbf{F}^\*(\mathbf{x}',q^j,t)} \frac{\mathbf{F}(\mathbf{x},q^j,t)}{\mathbf{F}(\mathbf{x}',q^j,t)} \approx \int dq^j \mathbf{N}^j(t) \mathbf{\varPsi}^\*(\mathbf{x}',q^j,t) \mathbf{\varPsi}(\mathbf{x},q^j,t) \tag{A9}$$

where *Nj*(*t*) = *Gj*(*t*)*Pj* / Ψ<sup>∗</sup>(*xj*, *qj*, *<sup>t</sup>*)Ψ(*xj*, *qj*, *<sup>t</sup>*). So, finally, a proper normalization of the Bohmian conditional states allows us to arrive to Equation (A1) from Equation (A7). Such normalization is already discussed in quation (11) in the text. The moral of the mathematical developments of this appendix is that open systems are more naturally described in terms of density matrix than in terms of conditional states when using the orthodox theory, while the contrary happens when using the Bohmian theory. Because of the additional variables of the Bohmian theory, the conditional states are a natural Bohmian tool to describe open systems.

### **Appendix C. Equations of Motion for Single-Electron Conditional States in Graphene**

As said in the text, graphene dynamics are given by the Dirac equation and not by the usual Schrödinger one. The presence of the Dirac equation on the description of the dynamics of electrons in graphene is not due to any relativistic correction but to the presence of a linear energy-momentum dispersion (in fact, the graphene Fermi velocity *vf* = 10<sup>6</sup> *m*/*s* is faster than the electron velocity in typical parabolic band materials but still some orders of magnitude slower than the speed of light). Thus, the conditional wavefunction associated to the electron is no longer a scalar but a bispinor. In particular, the initial bispinor is defined (located outside of the active region) as:

$$
\begin{pmatrix}
\psi\_1(x, z, t) \\
\psi\_2(x, z, t)
\end{pmatrix} = \begin{pmatrix}
1 \\
s e^{i\theta\_{\vec{k}\_\xi}}
\end{pmatrix} \Psi\_{\mathcal{S}}(x, z, t),
\tag{A10}
$$

where <sup>Ψ</sup>*g*(*<sup>x</sup>*, *z*, *t*) is a Gaussian function with central momentum *k c* = (*kx*,*c*, *kz*,*<sup>c</sup>*), *s* = 1 (*s* = −1) if the electron is in the conduction band (valence band) and *θkc* = *atan*(*kz*,*<sup>c</sup>*/*kx*,*<sup>c</sup>*). The wave packet can be considered as a Bohmian conditional wavefunction for the electron, a unique tool of Bohmian mechanics that allows to tackle the many-body and measurement problems in a computationally efficient way [25,62]. The two components are solution of the mentioned Dirac equation:

$$i\hbar \frac{\partial}{\partial t} \begin{pmatrix} \psi\_1\\ \psi\_2 \end{pmatrix} = \begin{pmatrix} V(\mathbf{x}, z, t) & -i\hbar \boldsymbol{\upsilon}\_f \frac{\partial}{\partial \mathbf{x}} - \hbar \boldsymbol{\upsilon}\_f \frac{\partial}{\partial z} \\\ -i\hbar \boldsymbol{\upsilon}\_f \frac{\partial}{\partial \mathbf{x}} + \hbar \boldsymbol{\upsilon}\_f \frac{\partial}{\partial \mathbf{z}} & V(\mathbf{x}, z, t) \end{pmatrix} \begin{pmatrix} \psi\_1\\ \psi\_2 \end{pmatrix} = \quad -i\hbar \boldsymbol{\upsilon}\_f \left( \vec{\boldsymbol{\sigma}} \cdot \vec{\nabla} + V \right) \begin{pmatrix} \psi\_1\\ \psi\_2 \end{pmatrix}, \quad \text{(A11)}$$

where *vf* = 10<sup>6</sup> *m*/*s* is the mentioned Fermi velocity and *<sup>V</sup>*(*<sup>x</sup>*, *z*, *t*) is the electrostatic potential. *σ* are the Pauli matrices:

$$
\vec{\sigma} = (\sigma\_x, \sigma\_z) = \left( \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}, \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \right). \tag{A12}
$$

Usually, in the literature, one finds *σz* as *<sup>σ</sup>y*, however, since we defined the graphene plane as the *XZ* one, the notation here is different. Then, we can obtain a continuity equation for the Dirac equation and then we can easily identify the Bohmian velocities of electrons as [44]

$$
\vec{v}(\vec{r},t) = \frac{v\_f \psi(\vec{r},t)^\dagger \vec{v} \psi(\vec{r},t)}{|\psi(\vec{r},t)|^2}. \tag{A13}
$$

By time integrating (A13) we can obtain the quantum Bohmian trajectories. The initial positions of the trajectories must be distributed according to the modulus square of the initial wavefunction, i.e., satisfying the quantum equilibrium hypothesis and thus certifying the same empirical results as the orthodox theory [44]. All this formalism was introduced in the BITLLES simulator in order to correctly model graphene and other linear band structure materials.
