Coherent manipulation of charge qubits in double quantum dots

The coherent time evolution of electrons in double quantum dots induced by fast bias-voltage switches is studied theoretically. As it was shown experimentally, such driven double quantum dots are potential devices for controlled manipulation of charge qubits. By numerically solving a quantum master equation we obtain the energy- and time-resolved electron transfer through the device which resembles the measured data. The observed oscillations are found to depend on the level offset of the two dots during the manipulation and, most surprisingly, also the on initialization stage. By means of an analytical expression, obtained from a large-bias model, we can understand the prominent features of these oscillations seen in both the experimental data and the numerical results. These findings strengthen the common interpretation in terms of a coherent transfer of electrons between the dots.


Introduction
Coherent control of nanoscale devices is one of the main topics of current research on electron transport in systems on the nanometer scale. A manifest realization of control is given by pump-probe schemes known from molecular physics [1,2]. In the context of electron transport the "pump" and "probe" steps consist in switching to and from a regime, where transport is largely blocked [3]. This provides an effective decoupling of the electronic system from the reservoirs and allows for a coherent evolution between pump and probe triggers. Such a setup has been successfully used to coherently control charge [4] and spin qubits [5,6] in double quantum dots (DQDs). The theoretical description of these experiments is a very demanding task, since a fully time-resolved calculations are necessary. Although several formalisms exist to address this issue, only very few numerical schemes are available in this context. Typically, the numerical approaches to time-dependent electron transport with arbitrary driving rely either on equations of motion for non-equilibrium Green functions (NEGF) [7][8][9][10][11][12] or generalized quantum master equations (QME) for the reduced density matrix [13,14]. Such calculations are very helpful to gain a deeper understanding of the mechanisms relevant for coherent control, since time-resolved quantities, such as occupations and currents, are readily accessible.
In this article we concentrate on the experiment of Fujisawa and coworkers on coherent manipulation of charge states in double quantum dots [4]. One of the hallmarks of the experiment is the observation of oscillations of the so-called number of pulseinduced tunneling electrons as a function of pulse length (see also Sec. 1.1). These oscillations are commonly associated with coherent tunneling processes between the two quantum dots during the manipulation stage. While this picture qualitatively explains the main features of the experiments, it neglects the influence of the initialization and the measurement stages on the coherent time-evolution. As we will show, consistently considering the whole pump-probe sequence provides even stronger evidence for actual coherent control. Moreover, the new insights may help to gain a deeper understanding of the dynamics in nanoscale devices.
First of all we briefly introduce the concept of charge qubits in double quantum dots and review the main results of the coherent manipulation experiment. The theoretical model and the relevant tools used in this work are presented in Sec. 2. Numerical results and a detailed analysis using an analytically solvable model are given in Sec. 3. The article concludes with a summary in Sec. 4.

Charge Qubits in Double Quantum Dots
From quantum information theory it is well known that in principle any two-level system may be used as a single bit of (quantum) information, which is then called qubit [15]. In laterally coupled quantum dots one may find such two-level systems by considering charge states denoted by (N , N r ), i. e., systems with N electrons in the left and N r electrons in the right dot. An electron tunneling from left to right corresponds to the sequence (N , N r ) → (N +1, N r ) → (N , N r +1) → (N , N r ). The probability for such a tunneling event is largest if the two charge states, | = (N +1, N r ) and |r = (N , N r +1) have a vanishing energy difference ε [16,17]. In the vicinity of such a resonance the coupled quantum dots can be described by a two-level system, which is characterized by the energy difference ε and the interdot tunnel coupling T c [4,17]. Correspondingly, the associated qubit is called charge qubit.
The Hamiltonian of the two coupled charge-states in the DQD, i. e. two coupled energy levels, reads The operators c † n (c n ) create (annihilate) an electron in the left (n = ) or right (n = r) dot, respectively. The interdot charging energy U suppresses double occupancy of the DQD. The time-dependence of the energies ε and ε r may be given by external gatevoltages. Figure 1 shows the corresponding energy scheme for a charge qubit additionally coupled to a source and a drain reservoir. In order to use a quantum system as a qubit, it has to fulfill at least the following three conditions [18]: (i) initialization of the qubit into a well defined state, (ii) application of unitary operations (quantum gates) and (iii) readout of the qubit state. Using a pump-probe scheme with rectangular bias-voltage pulses, these three steps have successfully been implemented for a charge qubit in a double quantum dot [4]. The scheme reminds in many ways of the usual pump-probe experiments with atoms or molecules [19][20][21]. In the present case the raising edge of the pulse triggers the coherent dynamics (pump), while the trailing edge starts the measurement (probe). Figure 1. Pump-probe scheme for the coherent manipulation of a charge qubit [4].
. e) Current (orange line) as a function of time. The stationary value J 0 is taken before the manipulation phase, t < 0, and for long times, t → ∞. The number of transferred electrons N (yellow-shaded area), without the stationary portion (gray-shaded), has contributions from the manipulation and the measurement phase.
The main idea of the experiment consists in suddenly switching between the Coulomb blockade regime and a transport regime by using a bias-voltage pulse. In the former case the DQD is effectively isolated from the reservoirs since sequential tunneling is strongly suppressed [22]. This provides the possibility to coherently control the charge qubit [3]. The transport regime is used to initialize the system and to readout the charge state after manipulation. The whole sequence is schematically shown in Fig. 1. Before and after the pulse, i. e. during the initialization and measurement phase, the sourcedrain voltage is V SD = V 0 and transport through the dot is possible since both charge states are within the transport window defined by the source-drain voltage. During the pulse, i. e. during the manipulation phase, the source-drain voltage is switched to V SD = V 1 and no transport is possible.
In the experiment, the time-averaged current is measured as a function of pulse length τ and energy difference ε ≡ ε r − ε . The latter is tuned by applying a suitable gate-voltage V g to the right quantum dot (not shown in Fig. 1). Notice, that due to capacitive couplings [4], the energy difference ε for a given gate-voltage V g also changes with the source-drain voltage, which is why we use ε 0 and ε 1 for high and low source-drain voltage, respectively, in Figs. 1a-c. The pulse is repeated with a repetition rate f rep = 100 MHz and by using a lock-in technique the pulse-induced current J p is obtained, which does not contain the asymptotic (stationary) current of the initialization and measurement phase. With J p one finally gets the number of pulseinduced tunneling electrons N = J p /ef rep [4]. It is then argued that this quantity is equivalent to the occupation of the right dot at the end of the pulse [4]. As we will show in Sec. 3.2, this assumption is not always valid. The experimental results show pronounced oscillations of N as a function of pulse length. Hereby, the frequency and the amplitude of the oscillations depend on the energy difference between the charge states. These oscillations are interpreted as a signature of coherent tunneling between the charge states (Rabi oscillations). Apart from the oscillations there are two other noteworthy features in the experimental results: for a fixed pulse length the function N is asymmetric around ε 1 = 0 and, in particular for ε 0 = 0, the number of pulse-induced electrons takes negative values. It has been argued that these features are artifacts of incomplete initialization or imperfect pulse shapes [4]. In the remaining part of the article we will demonstrate, that this must not be the case and the additional features instead strengthen the view of a coherent manipulation.

Theoretical Description
After specifying the Hamiltonian we briefly present the time-dependent description of the DQD using a quantum master equation. Furthermore we discuss a model which allows for deriving analytical expressions and facilitates a detailed analysis. Specifically, we consider the large-bias limit and use a description of the dynamics, which has originally been developed for photon-assisted transport through DQDs [23].

Setup
Apart from the DQD Hamiltonian [Eq. (1)], the total Hamiltonian also contains terms describing the electron reservoirs and the tunnel coupling, The reservoirs are described as usual by non-interacting electrons, where {b † αk } and {b αk } are the electron creation and annihilation operators for the αreservoir state k, respectively. The reservoir single-particle energies have the general form ε αk (t) = ε 0 αk + ∆ α (t) with the ∆ α accounting for a time-dependent bias. The tunneling Hamiltonian for the linear setup of two quantum dots in series, which are each coupled to a single reservoir, reads explicitly with T αk,n denoting the coupling matrix element between QD n = , r and the k-th mode of the respective reservoir α = L, R. In the second line of the tunneling Hamiltonian we have introduced the abbreviations B n , which will simplify the notation later on. Further, making the wide-band assumption renders the tunnel-coupling elements independent of the reservoir state, T αk,n = T α,n , which yields a flat spectral density,

Time-local Quantum Master Equation (TLQME)
The time-evolution of the density operator characterizing the total system is determined by the Liouville-von Neumann equation, where H is given by Eq.
(2). Since we are only interested in the properties of the DQD, we can trace out the reservoir degrees of freedom and thus obtain the reduced density operator σ = Tr res . In order to get an equation of motion for σ, we employ the time-convolutionless projection operator technique [24][25][26], which yields a time-local differential equation for σ [13], Here and in the following we set = 1. The TLQME is supplemented by the following auxiliary operators [13] Λ (x) These auxiliary operators are determined by the free DQD propagator U S (t, t ) and the reservoir correlation functions C αm (t, t ). The former is given in terms of the DQD where T is the usual time-ordering prescription. The correlation functions describe the influence of the reservoir on the system dynamics. In the current context they are given by [13] The Fermi function f α (ε) characterizes the initial state of reservoir α, and the spectral density Γ αmm has been defined in Eq. (6). In order to include the energy-level broadening induced by the back-action of the reservoirs on the DQD, we modify the correlation functions by multiplying with an exponential decay factor [27][28][29]. Consequently, we have C (xy) [28]. Although the equation of motion for the reduced density matrix is local in time, it is still very demanding to solve it numerically. An efficient method to propagate σ is the so-called auxiliary-mode expansion technique, which is based on a decomposition of the Fermi function [13]. As a consequence the operators Λ (x) andΛ (x) are also decomposed into a sum and for each term an individual equation of motion can be derived. The details of this procedure are given in Appendix A. The auxiliary-mode expansion has successfully been used in the context of time-nonlocal and time-local quantum master equations [13,14,30] and has recently been applied to non-equilibrium Green functions [12].
It remains to discuss the calculation for the time-resolved electric current through the tunneling barriers in a way consistent with the QME [31]. The electric current through the tunneling barrier α is given by the rate of change of the particle number in reservoir α [27,32], Plugging the formal solution of the Liouville-von Neumann equation [Eq. (7)] in the interaction representation into Eq. (11) and keeping only terms up to 2nd order in H tun , one finds within the time-local approximation for the time-dependent current [13], This expression is consistent with the TLQME, which is also of second order in H tun . Note, that in order to calculate the current one only needs the auxiliary operators Λ αm and the reduced density operator σ.

Markovian Approximation in the Large-Bias Limit (LBL)
Considering the experimental situation, one is lead to an even simpler description of the dynamics in the DQD. Firstly, for a large inter-dot interaction strength U only the following three states are relevant: |0 , | = c † |0 and |r = c † r |0 [23]. These correspond to an empty DQD, one electron occupying the left dot and one electron occupying the right dot, respectively. These three states are used as a basis for the reduced density operator in the following considerations. The second simplification arises due to the large bias in the experiment. In this case, using the Markovian limit of the QME [Eq. (8)] is well justified [23,33].
For the initialization and measurement phase (cf. Fig. 1) one obtains the following differential equations for the components of the reduced density matrix [23,34] Obviously, the coupling to the reservoirs introduces transitions between the charge states |0 , | and |r . The first two equations describe the evolution of the chargestate occupations, which change due to tunneling between the two dots and because of tunneling from and to the reservoirs. The remaining equation yields the dynamics of the coherences σ r = σ * r . Tunneling out of the DQD leads to loss of coherence. In the manipulation phase an effective decoupling of the DQD from the reservoirs is achieved by switching the source-drain voltage. An electron can still enter the DQD, but leaving the system is strongly suppressed. In contrast to the initialization phase, where tunneling out of the DQD leads to a vanishing coherence, here the dynamics stays approximately coherent. Therefore, we assume the following equations of motion during the voltage pulse, The rate γ has been introduced to account for additional processes leading to decoherence, such as background-charge fluctuations or back-action of the reservoirs on the DQD. In principle, these processes are also present in the other stages, but there they are dominated by the transport from source to drain. The expressions for the time-resolved currents are very simple and may, for instance, be extracted from Eqs. (13) and (14). They are explicitly given by in the initialization/measurement and manipulation phase, respectively [23,33].

Results and Discussion
In the following we will discuss a DQD system with parameters based on the experimental values [4]. We summarize all relevant quantities in Table 1. As sketched in Fig. 1, we consider a perfect rectangular pulse with duration τ and assume that the DQD is an stationary state at t = 0, when the manipulation phase starts. In the numerical calculation, at t = 0 the source-drain voltage switches from V 0 to V 1 and at t = τ it switches back. The level of the left dot is fixed at ε = −eV 0 /2. In order to get an energy-resolved picture we can shift the right level ε r = −V 0 /2 + eV g by changing the gate voltage V g at the right dot, similar to the experiment [4]. Due to capacitive parameter value capacitive level offset δε = 30 µeV interdot tunnel-coupling T c = 4.5 µeV tunnel rates Γ L = Γ R = 30 µeV source-drain voltages V 0 = 650 µV V 1 = 0 µV electron temperature 1/β = 10 µeV couplings there is an additional offset δε during the pulse and the energy of the right dot is ε r = −V 0 /2 + eV g + δε. Thus, the level offset ε = ε r − ε is ε 0 = eV g and ε 1 = eV g + δε.
In principle it is sufficient to specify either ε 0 or ε 1 , since this fixes automatically the other one. However, for convenience we will use both notations in the following.

Numerical Solution of the TLQME
We have investigated the pump-probe scheme presented in Sec. 1.1 using the TLQME given by Eqs. (8) and (9). The operators are represented in the basis {|0 , | , |r }. The resulting system of differential equations is propagated by means of an auxiliary-mode expansion described in Appendix A. In order to calculate the number of pulse-induced tunneling electrons the numerically determined current [Eq. (12)] was integrated over time and the stationary contribution was subtracted. Therefore, the following integral has to be calculated The dependence on ε accounts for both situations ε = ε 0 and ε = ε 1 , which have a fixed relation ε 1 = ε 0 + δε as explained above. J 0 refers to the stationary current reached at the end of the initialization and the measurement phase. The last term has been added to correct for the different stationary current during the manipulation phase (cf. Fig. 1e). Figure 2a shows the current J 0 at the end of the initialization phase as a function of the level offset ε 0 . This current is maximal when the energies of the two charge states coincide (ε 0 = 0). Also shown is the analytical result for large source-drain voltages [34,35]. In Fig. 2b the number of pulse-induced electrons N is shown as a function of pulse length τ and energy difference ε 0/1 . Hereby, the function N (ε 1 , τ ) is asymmetric around ε 1 = 0. Moreover, one can clearly observe an oscillatory behavior of N . The frequency of the oscillations increases with increasing values of ε 1 . This is explicitly shown in Fig. 2c, where the function N (ε, τ ) is shown for two energy differences, ε 0 = 0 and ε 1 = 0, respectively. For the latter case, the oscillations have the largest amplitude. Following the hypothesis that the number of tunneling electrons reflects just the occupation of the right dot at the end of the voltage pulse, one may naturally interpret the oscillations seen in Fig. 2b as Rabi oscillations between the two charge states [4], with the frequency given by Consequently, the maxima of N should appear at pulse lengths τ n = (2n+1)π/Ω with n = 0, 1, . . ., whereby the frequency Ω depends on ε 1 according to Eq. (17). These positions are indicated by white dashed lines in Fig. 2b. Around ε 1 = 0 the maxima of N are indeed found at the expected positions. However, for large ε 1 the maxima are clearly shifted compared to τ n . It is interesting to compare the numerical data for N (ε, τ ) with the instantaneous occupation of the right dot σ rr (t). This is shown in Fig. 3 for two energies ε. For ε 1 = 0 (ε 0 = −30 µeV) both quantities are almost identical to each other. In the  case ε 1 = 30 µeV (ε 0 = 0) one observes considerable deviations between N and σ rr . In contrast to the transferred electrons N , the occupation σ rr does have maxima at the times expected from the Rabi oscillations. Furthermore, N also takes negative values. Altogether, this is in clear contradiction to the assumption that the number of transferred electrons is just the occupation of the right dot. In summary, the numerical results and the experimental observations are qualitatively very similar. In particular, the described features of the behavior of N can also be seen in the experimental results. However, the reasons for the peculiar features, the asymmetry of N around ε 1 = 0, the negative values of N for large ε 1 and the apparent shift of the maxima, cannot be explained by the numerical investigations alone. Therefore, in the remaining part of this section we will consider an analytically solvable model based on the Markovian quantum-master equations of Sec. 2.3.

Model Calculations in the LBL
Starting point of the analysis is the initialization phase described by Eq. (13). Its stationary state provides the input for the manipulation phase. The stationary solution of Eq. (13) directly gives the the initial state within the considered pump-probe scheme. Settingσ = 0 yields the stationary populations σ , σ rr and coherences σ r The stationary state is in general a mixed state. Considering the case of strong coupling to the reservoirs, Γ T c , gives Obviously, in this case the initialization leads to an almost perfect localization of an electron in the left quantum dot. However, it is important to notice that at the same time finite coherences are unavoidable, which will have consequences for the manipulation phase.
For the manipulation phase it is convenient to introduce the following combinations of matrix elements of the density matrix σ, The respective equations of motion [Eqs. (14)] can be solved using a Laplace transformation [36], which is done in Appendix B. Under the assumption Γ L = Γ R = Γ one finds for the total occupation and for the difference Together, these two expressions allow for a calculation of the occupations in the DQD. Moreover, the current through the right barrier is given by for 0 < t ≤ τ , e Γ σ rr (t) for τ < t .
Hereby, the stationary current is J 0 = e Γ σ rr (0) = e Γ σ rr (∞). The typical time dependence of J R (t) for γ = 0 is shown in Fig. 1e. Thus, we are ready to get the number of pulse-induced tunneling electrons according to Eq. (16). This total number has contributions from the manipulation and measurement phase, i. e. N (ε, τ ) = N 1 (ε, τ ) + N 0 (ε, τ ), which are according to Eq. (22) The first contribution can explicitly be obtained from Eq. (21). For the case γ = 0 and s(0) ≈ 1, one finds s(t) ≈ 1 and consequently For the second contribution we use again a Laplace transform, cf. Appendix B. The result of this procedure is the following expression for the number of pulse-induced electrons in the measurement phase, Obviously, the value of N 0 depends in a non-trivial way on all elements of the reduced density matrix. For weak inter-dot couplings, T c Γ, we can use Eq. (19) for the initial density matrix. As in Eq. (19) we keep only terms linear in T c /Γ. Thereby we obtain, together with Eq. (24), a compact expression for the number of pulse-induced electrons This is a central result of this work as it shows quantitatively the differences between N and σ rr and, as we discuss in the following, explains the main features seen in the experiment. Firstly, Eq. (26) implies that for large energy differences, ε 0 T c , the number of pulse-induced tunneling electrons indeed corresponds to the occupation of the right dot at the end of the pulse. Therefore, it is a good measure for the occupation only for non-resonant initialization as already seen in Fig. 3. In this case the second term in Eq. (26) becomes very small. However, if the energy difference ε 0 , which is relevant in the initialization and in the measurement phase, vanishes, this term cannot be neglected. The corrections may lead to N taking negative values, which is shown in Fig. 4. Secondly, for energy differences ε 1 close to zero one finds N (ε, τ ) ≈ σ rr (τ ) . In this case, N is mainly determined by the population difference w(τ ), which is explicitly given by Eq. (21b). An important consequence of this dependence is the occurrence of an asymmetric behavior of w(t) as a function of ε 1 for a finite real part u(0) of the initial coherences. For example, considering times t n = (2n+1)π/Ω with n = 0, 1, . . . and assuming for the moment γ = 0, Eq. (21b) yields Obviously, this expression is only symmetric with respect to the energy difference ε 1 if u(0) = 0. In general, one finds w (t n ; ε 1 ) = w (t n ; −ε 1 ). In summary, the analysis within the described model provides a good explanation for the main features seen in the numerical results. These features result from an additional dependence of N on all the matrix elements of the reduced density matrix. Finally, we will briefly discuss damping of the oscillations. To this end, the analytical result given by Eq. (23) and the numerically obtained results from Fig. 3 are  shown together in Fig. 4. Without additional damping processes, i. e. γ = 0 in Eqs. (14), one finds for both energy differences undamped oscillations (dashed lines in Fig. 4). Hereby, the frequency and the positions of the maxima are in good agreement with the numerical results. For ε 1 = 30 µeV the number of pulse-induced electrons, N , takes negative values. An almost perfect agreement with the numerical results can be achieved by introducing decoherence into the manipulation phase (γ > 0). Taking γ = 0.46 µeV (full lines in Fig. 4) yields a very good description of the pulse-length dependence of N . Obviously, decoherence is a necessary ingredient to obtain a consistent picture. In the numerical calculations, the decoherence originates from the finite source-drain voltage and the broadening of the energy levels. This leads to a non-vanishing probability of the electron leaving the DQD [37]. In the experiment other sources of decoherence exist and these will typically dominate the damping of coherent effects. The most important processes in this regard are interaction of electrons with phonons, background-charge fluctuations and cotunneling [4]. From the energy-difference dependence of the damping rate, which is extracted from N (ε, τ ), one may infer about the nature of the relevant decoherence processes [3,4].

Conclusions
In summary, we have theoretically investigated a pump-probe scheme realized in a recent experiment on the coherent manipulation of charge qubits in double quantum dots [4]. To this end we have numerically simulated the pump-probe scheme using a time-local quantum master equation. The equations are solved by the auxiliary-mode expansion technique described in Appendix A, which, in general, provides a flexible and viable method to study time-resolved electron transport. The numerical results for the number of pulse-induced electrons N (ε, τ ) show good qualitative agreement with the experimental results. In particular, the main feature seen in the experiment, i. e. clear oscillations of N as a function of pulse length τ , are also observed. Moreover, two other characteristics of the experimental results are seen in the numerical data: the asymmetry of N around ε 1 = 0 and the occurrence of negative values. To adress these, so far unexplained, features we considered the Markovian limit of the QME in the respective pump-probe stages. The resulting equations can be solved analytically and lead to the main result of the article, namely the expression for the transferred electrons N . It turns out, that the value of N depends in a non-trivial way on all elements of the reduced density matrix. Only for large initial energy offsets, ε 0 T c , the number of pulse-induced tunneling electrons corresponds to the occupation of the right dot at the end of the pulse. For small values of ε 0 , larger differences between the occupation and the number of pulse-induced electrons are expected, which explains the occurrence of negative values. Further analysis shows that the asymmetry can be attributed to unavoidable initial coherences resulting from the initialization stage. Thus, the presented findings strengthen the common interpretation of the observed features in terms of coherent manipulation.
The expression for N (ε, τ ) may readily be tested against the experimental results. An interesting question for further studies would be the analysis of the decoherence rate as a function of energy difference. To this end, one could extent the numerical description by including phonons [38,39] or background-charge fluctuations. The decoherence rate may then be extracted from the numerical data using the analytical expression for the pulse-induced tunneling electrons from the LBL model.

Appendix A. Auxiliary Mode Propagation
In order to perform the energy integration in Eqs. (10) we expand the Fermi function f (ε) as a finite sum over simple poles with χ ± p = µ±x p /β and Im x p > 0. Here, µ is the chemical potential and β the inverse electron temperature. Instead of using the Matsubara expansion [40], with poles x p = iπ(2p−1), we use a partial fraction decomposition of the Fermi function [41], which converges much faster than the standard Matsubara expansion. In this case the poles x p = ±2 √ z p are given by the eigenvalues z p of an N F ×N F matrix [41]. The poles are arranged such that all poles χ + p (χ − p ) are in the upper (lower) complex plane. Employing the expansion given by Eq. (A.1), one can evaluate the energy integrals necessary to compute the reservoir correlation functions given by Eqs. (10) by contour integration. Thereby, the integrals in Eqs. (10) become (finite) sums of the residues. One gets for t > t with the auxiliary modes for reservoir α given by Hereby ∆ α (t) is due to the time-dependent single-particle energies ε αk (t) of the reservoir Hamiltonian [Eq. (3)]. The last term containing Γ αmm leads to the correct broadening of the energy levels as discussed in Sec. 2.2.
With the expansion (A.2) one can write the the auxiliary operators [Eq. (9)] in terms of the partial correlation functions. One readily obtains with the "partial" auxiliary operators An analogous expression holds forΛ

Manipulation Phase
Assuming Γ L = Γ R = Γ renders the equation for the total occupation s independent from the others,ṡ(t) = −[γ+2Γ]s(t) + 2Γ, with the solution given by Eq. (21a). The other components w(t), u(t), and v(t) are given by coupled differential equations which are solved by applying a Laplace transform. The resulting linear system reads  where we have used the Rabi frequency Ω from Eq. (17).

Measurement Phase
The Laplace transformation of Eqs. (13) leads to the following linear system  This system may be solved for the componentsσ (s),σ rr (s),σ r (s) andσ r (s). In order to calculate the number of pulse-induced tunneling electrons [Eq. (25)], one just needs σ rr (s). The stationary value of σ rr for t → ∞ can be calculated through the following limit, Using the expressions for the stationary density matrix [Eqs. (18)] finally yields Eq. (25).