Spectroscopy of hot-electron pair emission from a driven quantum dot

On-demand emission of individual electrons for the implementation of flying qubits and quantum electron-optics experiments requires precise knowledge and tunability of emission times and energies. Crucially, for confined electron sources such as driven quantum dots, the effect of local Coulomb interaction on these emission properties needs to be understood, in particular if multiple particles are emitted close in time or near-simultaneously. This paper theoretically analyzes electron-pair emission from an ac driven quantum dot, detailing the competing effects of the electron-electron interaction, the time-dependent potential forming the quantum dot, and of the quantum-state properties such as degeneracy on the emission times and energies. We complement a numerical analysis of the coherent Schr\"odinger evolution of two particles in a driven potential with a master-equation description for strongly interacting electrons tunneling stochastically into a weakly coupled conductor. This captures a broad range of different influences on the emitted particles and thereby guides the development of single-electron sources higher control over two-particle emission properties.


I. INTRODUCTION
The invention of on-demand single-electron sources [1][2][3] has opened up for recent research on electron-based flying qubits [4] and more generally on quantum optics experiments based on the controlled emission of single electrons [3,5,6].This endeavour relies not only on precisely timing the electron emission, but also on controlling the energy of the emitted particles.Indeed, recent experiments have manipulated electrons in an energy-selective manner [7,8] and time-resolved measurements [7,9] and tomographic techniques [10][11][12][13][14] have revealed information about the emission characteristics of different sources.
Reliable two-qubit gates within the electron-optical flying qubit platform require further, crucial progress in manipulating two-electron quantum states [3,6,8,10,12,13,[15][16][17].The most important distinguishing feature of electronic devices in this context is the possibly strong Coulomb interaction between electrons.This can play a role during particle propagation along wave guides (implemented via quantum Hall edge states) [18,19], or when particles are brought to collide in a tunable manner at quantum-point-contact based beam splitters [8,[15][16][17]20].In contrast, the study of the interplay with local Coulomb interaction, when electrons are emitted from a periodically driven quantum dot [1,2,[21][22][23], has until now mostly been limited to single-level quantum dots with constant, energy-independent tunnel barriers [24][25][26].Previous analyses of the emission energy have emphasized the role of the time-and energydependence of the potential-well forming the quantum dot [10,14,27].However, local Coulomb interaction is expected to strongly affect the emission times and energies of emitted electrons [28][29][30][31][32][33], which in turn are important for propagation velocities and further processing Gates on top of a conductor time-dependently modulate the potential landscape [Eq.(1)] to eject pairs of electrons.The positions x d,L , x d,C , and x d,R indicate the boundaries and center of the quantum-dot region, x f,L and x f,R set the regions, where a filter would measure time and energy at which electrons propagate through it.(b) Sketch of the emission process: the potential is driven such that the first electron can tunnel out, until only one electron remains in the quantum dot.Further driving allows also the second electron to tunnel out until the potential well has basically disappeared.
of flying qubits.
The present paper thus theoretically analyses the emission times and energies of two-electron emission processes from a quantum dot into a conductor separated by a gatetunable barrier [34], as sketched in Fig. 1(a).A particular focus of this study lies on the effect of local Coulomb interaction on such processes.
Typically, the emission into the conductor takes place far above the Fermi energy, and is hence referred to as hot-electron emission.The precise emission time and energy is controlled by a time-dependently modulated potential of the dot and its exit tunnel barrier, as sketched in Fig. 1(b).Electrons escape the dot once any of the addition energies rises to a point above the conduction band bottom at which the exit barrier becomes sufficiently transparent.The emission process and the separation of the particles in time and energy [7,8], as pictured in Fig. 1(b), is hence determined by a complex interplay between the modulated potential shape and the two-particle dot state dynamics due to quantum interference as well as Coulomb interaction.
To get hold of the various competing effects, we approach the problem with two complementary methods.First, we model the two-particle problem by numerically solving and analyzing the fully coherent two-particle Schrödinger dynamics in a 1 dimensional (1d) potential landscape with a dip between two barriers acting as the dot [Fig.1].This captures how different addition energies for the two emitted electrons due to charging energy and single-particle level splitting lead to different emission times, since the time-dependent potential modulation causes these energies to reach the onset of finite exit-barrier transparency at subsequent times.The simulation, however, also exposes how this straightforward time-energy separation can be nontrivially modified by the possibility of Coulomb interaction transferring energy between the particles during the first emission event.Similar interaction-based energy transfer between particles scattered at a barrier has been identified in Ref. [19].
Second, the emission times themselves can be enhanced or delayed differently for the two particles, thereby further altering the emission time-energy relations.Key underlying mechanisms include single-particle state degeneracies leading to a faster first emission [24,28,35,36], as well as time-dependent dot and exit barrier variations between the emission events shifting the required transition energy of the second particle relatively to the first.A suitable experiment to study these mechanisms is to measure the emission energies in sweeps of the dot potential ramp speed, as the latter changes the rate of energy increase relatively to the exit barrier escape rate.Simulating such sweeps using the full, interacting twoparticle Schrödinger dynamics is challenging.Results may also depend on many experimentally inaccessible details of the potential landscape, and the numerical cost of repeating the full simulation for many different ramp speeds is rather high.We tackle these challenges by complementing the two-particle simulation in 1d with a time-dependent quantum master equation for an effective two-orbital quantum dot, capturing in a simple way the interplay between Coulomb interaction, spatial dependence of single-electron states with dot-internal dynamics, and energy-dependent dot-conductor couplings.Unlike the coherent 1d Schrödinger evolution, this only applies to the weak dot-conductor coupling regime, gov-erned by electron emission processes due to stochastic tunneling.The crucial advantage is, however, that the much smaller state-and parameter space enables us to systematically classify and compare the competing effects impacting the emission times and energies as a function of the dot potential ramp speed.This approach in particular accounts for degeneracy-enhanced emission rates, and for classical ensemble-averaging which is typically performed in experiments as well.
We implement our investigation strategy by organizing the remainder of this paper as follows.The model and theoretical background of the coherent two-particle Schrödinger dynamics within the time-dependently modulated, 1d potential landscape is introduced in Sec.II; Sec.III then discusses the emission times and energies obtained from the corresponding numerical simulations.The quantum master equation approach for the dynamics of the simplified quantum-dot model is set up in Sec.IV; the resulting effects in the energy-time spectroscopy of the emission process are detailed in Sec.V. We highlight the main insights and open questions emerging from the two complementary approaches of describing two-particle emission times and energies in the concluding section VI.This paper also contains an appendix with details on the numerical simulation of Sec.II and Sec.III and on the quantum master equation approach.Throughout the manuscript, we set ℏ = k B = |e| = 1.

II. TWO-PARTICLE SIMULATION IN 1D
Experimental realizations of periodically driven quantum-dot electron pumps typically rely on electrostatic confinement of two-dimensional electron gases via tunable gates, as sketched in Fig. 1(a).A full theory of two-electron emission from such dots would need to determine the time evolution of a many-body electron system with Coulomb interaction in an inhomogeneous and time-dependent potential.Some approximations hence need to be made to make the problem manageable.
In a first step, Secs.II and III, we make two approximations: First of all, we decide to treat a two-particle problem, which is motivated by the fact that the interplay between the emitted 'hot' electrons and the Fermi sea in the contacts is rather weak.Furthermore, in order to keep the numerical simulation tractable, we model the system in one dimension.

A. Hamiltonian with time-dependent potential for hot-electron emission
The potential landscape of the dot and its environment we consider is sketched in Fig. 1(a).With V 0 representing the conduction band bottom high above the Fermi sea, we envision the three gates defining the quantum dot to have a Gaussian-like effect on the potential landscape: The latter approximation is appropriate if emission takes place during some small fraction of a periodic drive signal [10] with an amplitude A ∼ V 0 , such that the only relevant transition energies E lie close to the emission point, We set up the system Hamiltonian by discretizing the potential landscape V (x ≥ 0, t) [Eq.(1)] into R+1 points defining R equidistant real space intervals δx = L/R on some large, but finite length L ≫ L d .Importantly, we explicitly model only the dynamics of the two particles initially occupying the dot.This relies on the assumption that for hot-electron emission with sufficient energy splitting between empty conductance and completely filled valence band, the latter appears to be nearly chargeless to the emitted electrons, apart from some overall screening effect.Hence fixing the total particle number to 2, we obtain a Hilbert space spanned by the (R + 1)(2R + 1) orthonormal, anti-symmetric states , rσ⟩ in discrete positions x r = rδx, r ∈ {0, . . ., R} with spin-z projection σ =↑, ↓.The particles are created from the vacuum |Vac⟩ by the fermionic field operators Ψ † rσ obeying the usual anti-commutation rules where N rσ = Ψ † rσ Ψ rσ are the occupation numbers at positions x r , and the kinetic energy The coupling λ = 1/(2m eff δx2 ) with effective electron mass m eff arises when discretizing the kinetic term Ψ † rσ where k e is the bare value and κ accounts for material properties screening.Note that for a spin-up and spindown electron in the same position x r , we regularize the Coulomb potential heuristically by setting the denominator in Eq. (2c) to δx/2.This is appropriate if the discretization step δx is small enough for the two electrons to reside in a state with a typical separation larger than δx prior to emission.Only in situations irrelevant to our analysis may our δx choice still be insufficient, such as collisions of counter-propagating particles.While we here put forward one -experimentally relevant-choice for the confinement and interaction potentials, modifications of this could be envisioned in the future.One extension would be to implement a Thomas-Fermi model for the effect of a background electronic density on the screening of interactions or via image-charge screening [37], as used for example in Ref. [15].Furthermore, one could optimize the confinement potential for the emission, similarly to how it has recently been shown for the process of loading the dot [38].

B. Unitary two-particle evolution and spectroscopy
We analyze the emission process, which is induced by ramping up V d,C (t) until the two electrons have enough energy to pass through the exit barrier at x d,R into the environment channel, at a potential energy V 0 .We start by setting the initial state |Φ(t = 0)⟩ to the ground state |GS⟩ of H 2P (t = 0), for which the potential dip V d,2 (0) ≫ E d C is deep enough for the dot to be stably occupied by both electrons, given a typical dot charging energy See, e.g., Ref. [39][40][41] for details of the loading process.
The ground state and its subsequent evolution due to the gradually lifted potential dip V d,2 (t) → 0 are obtained numerically using the Hamiltonian given in Eq. ( 2).We have developed a GPU implementation of a direct, norm-preserving leap-frog solver for a large, but bounded system, see Appendix A. This offers non-prohibitive run times even for boundaries far enough away that unphysical reflections do not affect the particle emission.In particular, it enables us to use the kernel polynomial method [42] to compute the full emission energy distribution as a function of time.
With the time-dependent two-particle state |Φ(t)⟩ at hand, we track the emission of the two electrons from the dot with the expectation values ⟨O⟩(t) = ⟨Φ(t)|O|Φ(t)⟩ of various observables O. Namely, next to the Coulomb energy H Coul and the spatial charge density |Ψ 1 (r, t)| 2 = σ=↑,↓ ⟨Φ(t)|N rσ |Φ(t)⟩/δx, we extract the charge and kinetic energy in the dot (l = d) as well as in a region away from the dot representing the energy filter (l = f), see Fig. 1(b).Concretely, we obtain the charge N rσ and the kinetic energy H l kin by modifying Eq. (2b) by replacing in all sums the lower limit 0 by R l 1 and R in the upper limit by R l 2 ; the indices R l 2 > R l 1 delimit the spatial ranges [R l 1 δx, R l 2 δx] on which these regions are defined.The time-dependent charge and kinetic energy in the filter region as function of time are the main quantities of interest in the analysis of Sec.III.Their properties are further supported by the time-and energy resolved emission distribution extracted from the spectral density in the filter, Technically, this distribution is analogous to spectral decompositions in Green's function approaches, see e.g.Ref. [43].From a practical point of view, it connects to the experimentally broadened distribution functions of Refs.[7,10], see also the averaged emission distribution (13) of Sec.IV.However, in contrast to the stochastic effects discussed later in Secs.IV and V, Eq. ( 4) only emerges from quantum effects and two-particle interaction.Note that this quantity should not be confused with quasi probabilities, such as for example the Wigner function.

III. TWO-PARTICLE DYNAMICS -RESULTS
We now discuss the emission dynamics from the twoparticle simulation described in Sec.II, analyzing the time-dependent expectation values of the charge N d in the dot, of the charge N f and kinetic energy E f kin in the filter region, as well as of the energy-and time-resolved filter region spectral density [Eq.( 4)].We compare the cases of vanishing (κ = 0) and finite (κ ̸ = 0) Coulomb repulsion strength, both for opposite spins ↑, ↓ initially occupying dot orbitals of equal single-particle energy, and for two ↑ spins initially residing in different, energy-split dot orbitals due to Pauli exclusion.
Videos of the time-dependent particle density |Ψ| 2 along the 1d potential landscape for the different scenarios can be found in the supplementary material [44].
In the following, we focus on the situation where particles are emitted separately in time, such as found in recent experiments [7,8].This regime furthermore enables to clearly attribute which of the effects described in the following stem from Coulomb interaction, and which are merely related to energy splittings or to time-dependent potential variations.To achieve comparable situations independently of the absence or presence of Coulomb interaction, we choose for the simulation in Fig. 2, a smaller effective mass m eff than the literature values known, e.g., for InAs or GaAs.However, the interaction-related effects identified below remain equally important for larger effective masses m eff , see Appendix B.

A. Effects of single-particle level splitting
We start by analyzing the case of fully screened Coulomb interaction (κ = 0) between two spin-↑ particles with different orbital energies in the dot before emission, see Fig. 2(a,b).Due to this energy splitting, the timedependent V (x, t) modulation causes the dot to emit the first particle clearly separated from the second, as shown by the successive dot particle number reduction from two over one to zero (black line in panel (a)).The charge and energy expectation values N f , E f kin (blue and red line) accordingly indicate one particle after another to arrive in the filter region, with a time difference to emission set by the propagation velocity of the coherent wave packet [44].Notably, the second emitted particle, once fully inside the filter region (N f → 1), has an overall higher energy than the first one.We attribute this to the fact that by lifting the Gaussian potential dip in Eq. ( 1) with the prefactor is also slightly raised, as demonstrated in the wave propagation animations [44]; the second particle emitted at a later time therefore needs to be lifted to a higher energy before it can escape [7].
While the left panel visualizes the main results to be presented here, this trend is also confirmed by the behavior of the spectral density.The spectral density ϕ 2P (t, E) in Fig. 2(b) is the time-dependent representation of the traveling wave packets in terms of the discrete energy eigenstates of the filter region, i.e., of a one-dimensional problem with finite length x f,R − x f,L = 1350nm; it hence reveals how the particle energies are distributed as a function of time, where the peak at E = 0 means that no particle is present in the filter region.The plot exhibits two time intervals with non-zero weight at finite ener- gies, corresponding to the two emitted particles just as in Fig. 2(a).The discrete values with nonvanishing density stem from the finite size of the filter region; their broadening stems from the truncation scheme that we employ here.Interestingly, we observe a steep yet finite E − tslope of the spectral density, comparable to the potential ramp speed v r (see the turquoise dashed lines).This reflects the experimentally observed effect of the potential drive increasing the energy of a wave packet during its emission, a so-called energy-time chirp [10].
Since the high-energy wave function components also propagate faster through the filter, the spectral density at lower energies remains finite for the longest time duration.In fact, these slowly traveling components of the wave function are also the main reason for the nonvanishing N f in between the two peaks in (a), and for an N f even slightly larger than 1 in the second peak.Furthermore note that during periods in which the wave packets are only partly localized in the filter region (N f < 1), there is also partial overlap with low-energy filter states that do not exhibit any significant spectral density once the particle is fully inside the filter region.This is highlighted by the horizontal, turquoise-dotted line in Fig. 2(f), and creates the visual appearance of voids in the spectral density at low energy.
In Fig. 2(c,d), we contrast the above case of two equal spins subject to Pauli exclusion against the degenerate situation with two opposite spins initially occupying the same dot orbitals.In the absence of Coulomb interaction, κ = 0, this implies equal addition energies for the two particles, and thus simultaneous emission.The filter-region charge and kinetic energy accordingly exhibit only one peak, where the energy is approximately twice as large compared to the second peak of the spin-split case [Fig.2(a,b)].The spectral density in panel (d) furthermore features much more closely spaced lines once both particles are localized well within the filter region.The reason for this is that the spectral density is composed of all possible pairs of (not equally spaced) singleparticle energies stemming from the contributions of the two emitted particles.

B. Energy transfer via Coulomb interaction
Having covered the non-interacting limit, we proceed with Fig. 2(e-h) showing the corresponding dynamics in the presence of a finite Coulomb interaction between the two particles, κ ̸ = 0. Focussing on the relative emission energy of the first and second particle, we see that the first emitted particle carries significantly more energy than the second.This is the opposite outcome to the κ = 0 case [Fig.2(a-d)], and largely independent of the particle spin orientations, which generally play a less prominent role, c.f., Fig. 2(e-h).
The weak sensitivity to the particles' spins is intuitively clear since the Coulomb repulsion keeps the particles in the dot apart, such that fermionic anti-bunching becomes less of a factor.The main feature we highlight here is, however, the fact that the first emitted particle now always has more energy than the second.We attribute this to the combination of our system being 1d, and a Coulomb-interaction mediated energy exchange between the particles at the first emission event.Namely, if both electrons are inside a one-dimensional dot in which they cannot move about each other, they are bound to repel each other into orbital configurations higher in energy both in response to Coulomb or exchange interaction.However, once the first electron starts to leave, the two-particle Coulomb interaction allows the remaining electron inside the dot to relax to a lower energy by transferring the energy difference as additional kinetic energy to the escaping electron.This is the key difference between a mere single-particle level splitting related to Pauli exclusion, and the transition energy difference introduced by Coulomb repulsion.In fact, a close comparison between panels (e,f) and (g,h) reveals that the E f kin difference between first and second electron is even slightly higher for opposite spins without Pauli exclusion.This suggests an even more efficient energy exchange, possibly due to the fact that the particles can approach each other even more without fermionic anti-bunching.
These features described above are confirmed by the properties of the spectral density, plotted in panels (f) and (h).The higher-energy contributions of the first emitted particle compared to the second are clearly visible, whereas the second emission strongly resembles the non-interacting, spin-split case, shown in panel (b).The higher kinetic energy of the first emitted particle goes along with a higher propagation velocity, which can be seen in the smaller time-window in which the high-energy states are occupied in the filter region.
Another Coulomb-repulsion related effect on the first emitted particle is simply the static, 1/r-dependence of the potential originating from the second particle residing in the dot.This features most prominently when V (x, t) is chosen to be a sharp potential well with a flat bottom shifted up in time.In this case, the first emit-ted wave packet remains sharper over time and travels faster than the second, since the first particle effectively runs down the 1/r potential while the second diffuses out into a flat potential landscape, see box-potential animations in Ref. [44].This also highlights the more general fact that the potential landscape needs to fall off into the environment in order to see well confined wave packets2 .
Additional smaller features occur in the spectral function of particles in the filter region, which are expected to derive from the precise realization of the potential landscape and the resulting complex coherent two-particle wave-packet dynamics.Attributing these features to specific physical mechanisms is hindered by the large available parameter space and the numerical cost of sweeping this space.This includes variations of the ramp speed v r as an experimentally feasible method to expose nontrivial deviations from the above identified time-energy chirp.In the following, we therefore switch to the complementary quantum master equation description of an effective, two-orbital quantum-dot, in which these effects can be studied more systematically.

IV. MASTER EQUATION FOR EFFECTIVE DOT MODEL
In addition to highlighting the importance of dotinternal Coulomb scattering, the numerical analysis in Sec.III shows more generally that the emission times and energies are dependent on a complex interplay between interaction, interference and the spatial structure of the potential landscape.However, isolating these features within the full two-particle model is tricky and timeconsuming due to the many different parameter interdependencies, and the analysis would be tied to a 1d setting.We therefore approximate the emission dynamics with a simpler, fermionic two-orbital dot with local Coulomb interaction and a tunnel-coupled reservoir.Given sufficient spatial confinement for two-particle emission, this model is appropriate for zero-to three dimensional quantum dots.It furthermore allows for a full many-body master-equation study that can to a significant extent be carried out analytically.This analysis yields the average emission energies at a given emission time but does not provide information about the subsequent particle propagation, which is instead accessible via the complementary numerical approach of Secs.II and III and the related movies in Ref. [44].

A. Quantum dot coupled to reservoir
Our model contains two fermionic single-particle states ℓ = R,L with local Coulomb charging energy in a time-dependently driven potential coupled to a reservoir [Fig.3(a)].In the corresponding Hamiltonian H tot = H(t) + H res + H tun (t), the dot dynamics are given by (5) This includes the driven onsite-energy ϵ R (t) = ϵ(t) = ϵ 0 − A cos (Ωt) with offset ϵ 0 and amplitude A, a singleparticle level splitting ∆ϵ = ϵ L (t) − ϵ R (t), an internal transition amplitude τ ≥ 0, and the spatially independent, but possibly time-dependent onsite interaction strength U (t) > 0. While the latter cannot describe dotinternal Coulomb scattering as seen for the extended 1d model [Sec.II B], the time-dependence does reflect how possible dot size variations during the driving cycle can affect the Coulomb energy, as detailed in Sec.V A. The occupation number operators ) and annihilation (d ℓ ) operators for single-particle state ℓ = R,L.
Note that despite the suggestive notation, we do not yet specify the physical nature of these states: ℓ may, at this point, not only refer to two localized levels ℓ = R,L, but also to a spin projection (ℓ =↑, ↓) in a spinful setup, or to two orbitals (ℓ = s,p) in a strongly spin-polarized system.
The effectively non-interacting Hamiltonian H res = k,ν ϵ kν c † kν c kν of the electronic reservoir describes a quasi-continuum of fermionic modes with wave number k and all discrete quantum numbers ν necessary to characterize these modes.The respective creation and annihilation operators are denoted by c † kν , c kν .The charge conserving coupling between these modes and the dot state is governed by The coupling amplitudes τ kν (t) account for timedependently varying couplings due to the driven potential landscape; the relative amplitude X RL kν models that depending on the setup, one dot mode ℓ may couple more strongly to the environment mode kν than the other dot mode.Averaged over all environment modes as carefully prescribed in Sec.IV B, the coupling amplitudes determine the typical rates (7) for tunneling in or out of dot mode ℓ = R,L with Kronecker deltas δ ℓR , δ ℓL .Importantly, the energydependent barrier transparency Γ ℓ (E), beyond the common wideband limit [45], should capture that the potential V (x, t) [Eq.( 1)] induces emission once V d,C (t) is small enough for a dot transition energy E to exceed the potential energy outside the dot.We achieve this with a barrier height V 0 > µ above which Γ ℓ (E) increases smoothly from an uncoupled dot with Γ ℓ (µ The modified sigmoid ) implements the ramp-up from exactly 0 to Γ ℓ using an analytic function in the finite interval 3 ).We address this further below, see Sec.V C, to analyze the effect of dot potential modulations on the tunnel barrier, which we have already identified for the unitary two-particle evolution above as the cause of the different energy peak heights in Fig. 2(a).

B. Dot state dynamics
We describe the time-dependent dot state with a master equation for the reduced density matrix ρ(t), i.e., with the reservoir modes traced out as detailed, e.g., in [46,47].Indeed, the emission of electrons far above the Fermi energy of the conductors is equivalent to a situation with infinitely large bias, for which a weakcoupling master equation approach is applicable [48][49][50].The broadening in the currents calculated in this section hence stem from the stochastic nature of the tunneling.
With the above in mind, we use the approach of Ref. [51] to derive a Lindblad master equation for the dot tunnel-coupled to the bath.Care needs to be taken since the tunneling (6) couples a single environment mode kν to both dot states ℓ = L, R for X RL kν ̸ = 0.This introduces up to four independent decay channels per particle-and hole process corresponding to the two possible orthogonal hole-like operators αd R + βd L and the analogous particle-like fields.As we have not yet tied ℓ to any specific single-particle basis, and since the shape of H(t) is invariant under any unitary transform (after properly regauging the fields to ensure τ ≥ 0 and by removing any constant energy shift) of the single-particle basis [Eq.( 5)], we demand that each of the maximally 4 orthogonal channels in H tun are proportional to only one of the four dot operators d ℓ , d † ℓ .The corresponding master equation for the density operator ρ(t) then reads with Lamb shift Λ [Eq.(C1)] and Lindblad operators The states |i⟩, |j⟩ are the four instantaneous many-body energy eigenstates of the dot Hamiltonian H(t), and E ij = E i − E j the corresponding energy differences, where our notation suppresses their time-dependence for better readability.The operator d ηℓ = δ η+ d † ℓ + δ η− d ℓ combines the corresponding dot creation and annihilation operators, and the Fermi function f for η = +, respectively 1 − f for η = −.The master equation ( 10) relies on the instantaneous-time approximation, in which the Lindblad operators (11) only depend on the parameters at the current emission time t [51][52][53].We expect that this leads to reliable results for driving that is limited by AΩ ≪ (E b − V 0 ) 2 [50,54].For fast driving of tunablebarrier dots, see also Ref. [8,27,55].
As discussed above, we represent the hot electron setting by taking the T → 0 limit of the Fermi function, i.e., f η (x) → Θ(η [µ − x]) with Heaviside function Θ(x).Moreover, in the lowest-order coupling approximation assumed here, the environment-induced Lamb shift Λ only modifies the unitary dynamics in the single-particle sector of the local Hamiltonian H.The splitting ∆ϵ and coupling τ are hence generally shifted for coherent L-R rotations, but the jump rates (11) are only affected by the bare energies.Also, our discussion focuses on a tunnel coupling Γ/U = 0.001 much smaller than the range |E b − V 0 | ∼ 0.05U on which the barrier changes its opacity [Eq.( 8)], so that |Λ| ∼ Γ ≪ |E b − V 0 |.Thus, while we do account for the Lamb shift Λ in solving Eq. ( 10), it is in fact irrelevant for the emission dynamics in all cases considered in this paper, i.e., both for ∆ϵ = τ = 0 and for |∆ϵ|, |τ | ≫ Γ.More detailed derivations of the master equation are given in Appendix C.

C. Time-resolved emission energy
The observables of interest for our master-equation based analysis are the ensemble-averaged time-dependent charge current and energy current carried by particle transfer into the environment.Local charge conservation [H(t), N ] = 0 and charge conserving tunneling implies that the particle current into the environment equals the particle current out of the dot.Furthermore, in the weak coupling limit, no energy is stored in the barriers [56], ⟨H tun ⟩ = 0, so that the energy flow into the reservoir is the flow out of the dot minus the energy exchange with the work source.We hence write In the regime studied here, the energy flow is concomitant with net particle transfer; hence the kernel W I,E involves only those state transitions of the master equation (10) with changing dot occupation, as detailed in Appendix D.
We start the dynamics (10) from double occupation, , and an initial dot potential ϵ(t = 0) = ϵ 0 − A, sufficiently below the energy where the exit barrier becomes opaque to keep the electrons stable in the dot, 2(ϵ 0 − A) + U ≪ V 0 .The driven transition energies E ij eventually reach the opaque-barrier interval Γ ℓ (E > V 0 ) > 0 causing a particle to be emitted.The resulting particle current I N (t) is a measure for the emission probability at time t, and the ratio E mit (t) = I E (t)/I N (t) is associated to the energy emitted per particle.We combine the above quantities to define the time-and energy-resolved emission distribution Unlike experimental data in, e.g., [7,10,15] or the twoparticle spectral density (4), this distribution is sharp in E because the energy E mit (t) per emitted particle is already defined in terms of both classical and quantum averaging.The observed energy spread in experiments [9,10,15] is considered to stem to a large extent from the noise and timing jitter of pump and detector drive signals, obscuring the quantum uncertainty.The reconstructed Wigner distribution [10] is likely of mixed states, in other words, a classical ensemble of different emitted states, rather than of a pure state.Our simplified, low-energy model in Eq. ( 13) captures how an uncertainty in emission time results in a classical uncertainty in emission energy; quantum mechanical energy smearing is hence not at the focus of the present paper, but could instead potentially be obtained from the full two-particle approach previously described in Sec.II.
We consider a central question to be how the two emitted particles differ in the distribution ϕ(t, E) -and in particular in their average emission times and energiesas a function of the system-and driving-parameters, with special focus on the potential ramp speed ∂ t ϵ(t).Assuming two particles released within two separated time intervals, as sketched in Fig. 3(b), with the first transition from double to single occupied dot and the second from single occupied to empty dot indicated by x ∈ {1, 2}, we obtain these averages by integrating over ϕ(t, E): 14) Sweeps of the offset potential ϵ 0 , within a range fulfilling V 0 ≤ ϵ 0 + U + A, modify the potential ramp speed ∂ t ϵ(t) in the time in which the addition energies are in the opaque-barrier interval.This modifies the emission energies and emission times t x (as well as the ramp speed at the time of emission).We sketch in Fig. 3(c) a parametric time-energy plot of (t x , E mit,x − ϵ 0 ) for the two-particle emissions, indicated by stars for one specific choice of ϵ 0 .When sweeping ϵ 0 , we expect the emission energies to lie on the indicated lines, extrapolating with the cosineshape of the addition energies, time-dependent via ϵ(t).
To connect the obtained trajectories to available experimental data [7,8,57], we compare the two emission curves in two different ways.First, we fix the offset ϵ 0 , meaning we compare the emission times and energies of the first and second particle for the same driving protocol : Compare to the horizontal distance betwen the stars shown as an example in Fig. 3(c).If a particle with energy E was instantaneously emitted as soon as E > V 0 with constant V 0 , we would not expect any difference between first and second particle, ∆E mit → 0. We, however, show in Sec.V that ∆E mit ̸ = 0 arises due to the finite emission time ∼ 1/Γ ℓ competing with the driving parameters Ω, ϵ 0 determining the ϵ(t)-ramp speed at emission, and with dot-internal transitions ∼ 1/τ .Second, we compare different driving protocols with different ϵ 0 , ϵ ′ 0 > ϵ 0 , black vertical line in Fig. 3(c), to obtain the apparent transition energy difference between first and second particle: For well separated emission events, we intuitively expect ∆E tr to be given by the difference between the largest possible double-to-single transition energy and the largest available single-to-zero transition energy after the first emission event.Based on the dot Hamiltonian Eq. ( 5), this difference is given by the sum of interaction strength and splitting between the two singleparticle states, ∆E tr ≈ ∆E ref tr = U + 2 |τ | 2 + ∆ϵ 2 /4, see Eq. (C2).While we indeed find this to be mostly the case, there are, however, deviations from ∆E tr ≈ ∆E ref tr that we further explore in Sec.V D.

V. TIME-RESOLVED EMISSION SPECTROSCOPY -RESULTS
To discuss the time-resolved emission spectroscopy obtained from the master equation, we first analyze the difference in emission energy of the two particles ∆E mit and in particular its dependence on the ramp speed v r,1 = AΩ sin(Ωt 1 ), as shown in Fig. 4. To stay in line with previous experiments [7,15], we tune v r,1 by sweeping the offset ϵ 0 determining at which phase, and thus at which slope of the cosine-shaped ϵ(t) the particles are emitted.This entails two features visible in all results shown in the following subsections: first, all graphs begin at a finite v r,1 , determined by the minimum offset requirement ϵ 0 + A ≥ V 0 for the driving cycle with amplitude A to emit both particles from the dot.Second, the ∆E mit -lines bend downwards when approaching this minimum ramp speed from above.This stems from the stronger curvature of the cosine-shaped driving potential around the turning point ϵ(t) ≈ ϵ 0 + A ≈ V 0 ; it causes the two ramp speeds to differ significantly, v r,1 ≫ v r,2 = AΩ sin(Ωt 2 ), and thus the first particle to be emitted at a higher energy, ∆E mit < 0.
In the following we provide a detailed analysis of how emission times and energies are impacted by local Coulomb interaction effects as compared to effetcs due to the potential landscape creating the dot (here visible as energy-and dot-state-dependent couplings).We address the impact that quantum-state degeneracies have enabled by Coulomb interaction, Sec.V A, the effect of coupling asymmetry and internal dynamics (as they could realistically arise also from local Coulomb repulsion, see Sec.III), Sec.V B, compare to the effect that results from driving-dependent tunnel couplings, Sec.V C, and show the effect of a time-varying local Coulomb interaction, Sec.V D.

A. Interplay between interaction and degeneracy
We start by studying the influence of degenerate orbitals and of the coupling asymmetry, see results presented in Fig. 4(a).We, therefore, first set ∆ϵ = τ = 0 in the Hamiltonian (5).In this case of degenerate single-particle levels, as is commonly the case for spindegeneracies for example, the first emitted electron can be either the one occupying level R or level L, leading to an increased rate Γ R + Γ L for the emission of the first particle [24].Since the second particle is emitted more slowly than the first, it is also emitted at a higher energy due to the continuous lifting of the energy level by the driving potential.This effect of the level degeneracy can be seen in Fig. 4(a), where ∆E mit is shown for 3 different situations with zero detuning, ∆ϵ = 0.The red and the orange-dashed lines show this degeneracy-induced effect, demonstrating a ∆E mit that increases with increasing ramp speed v r,1 .This effect is stronger for the red line, where the equal rates Γ R = Γ L are smaller than for the orange line.Therefore the time that passes until the emission of the second particle is larger such that the level has been shifted to a higher energy.The blue-dashed line shows that this degeneracy effect indeed yields similar results to having one of the levels being coupled more strongly than the other.In contrast, as soon as the dotorbital degeneracy of two equally coupled levels is lifted, ∆ϵ ̸ = 0, the degeneracy-induced delay of the second emission vanishes, as shown by the black dashed-dotted line that is close to zero and almost independent of the driving speed.
This effect due to the degeneracy of single-particle energy levels can only be observed in the presence of strong Coulomb interaction, separating the emission of the two particles from the dot.Note that the impact of degeneracy on decay rates [58] of a quantum dot has been observed in various experimental settings [35,59,60] and further schemes have been theoretically proposed to read out such differences in decay rates [31,61].

B. Coupling asymmetry and internal dynamics
One possible factor governing the emission spectrum in the 1D two-particle dynamics discussed in Sec.III is the dot-internal repulsion due to Pauli exclusion and Coulomb interaction, pushing one particle further away from the tunnel barrier than the other.The resulting asymmetric coupling to the environment can be reflected within the master equation approach by setting Γ L ̸ = Γ R .Calculations which illustrate the result of an (effective) coupling asymmetry are shown in Fig. 4(b).
We have previously seen in Fig. 4(a) that ∆E mit grows with increasing v r,1 for degenerate levels (∆ϵ = 0) or if the higher-lying state was at least equally if not stronger coupled to the environment.If we instead assume a much more weakly coupled high-lying state (∆ϵ > 0, Γ L ≪ Γ R ), the slower tunneling rate in the first emission event leads to a significantly higher emission energy than for the second electron emitted from ϵ R (t) with rate Γ R .As long as the two emission events remain clearly separated in time, this leads to exactly the opposite situation with a negative ∆E mit < 0 becoming more negative with growing v r,1 , as shown by the green-dashed line in Fig. 4(b).Note however, that if v r,1 is increased even further, approaching or exceeding Γ L , the lower addition energy ϵ R + U may cross V 0 and induce emission from state L before the higher-lying particle from state R could leave the dot.This can result in almost simultaneous emission and a ∆E mit approaching positive values again, as indicated by the positively curving green-dashed line in Fig. 4(b) for v r,1 /(AΩ) → 1.
A relevant situation where such an asymmetric tunnel coupling, and hence a negative ∆E mit (v r,1 )-slope can be realized is when two localized orbitals can be occupied in the dot, and tunneling to the environment from one of them can occur only through coherent coupling to the other.Concretely, we realize this by setting Γ L = 0 and turning on a finite τ > 0 in the Hamiltonian (5).The ratio |∆ϵ/τ | then determines to what degree the single-particle eigenstates of the dot are localized in orbitals L and R, interpolating between perfect (anti-)bonding states (|∆ϵ/τ | → 0) and perfectly localized states (|∆ϵ/τ | → ∞).With state 2 decoupled from the bath (Γ L ), a splitting ∆ϵ = 0 results in equally coupled bonding-and anti-bonding states, whereas a positive splitting ∆ϵ > 0 yields a higher-lying state with weaker effective environment coupling, see Appendix C. Figure 4(b) accordingly shows a very weak v r,1 -dependence of ∆E mit for ∆ϵ = 0 (turquoise solid line) due to the nearly symmetric environment coupling, and a clearly negative ∆E mit (v r,1 )-slope due to the strong couplingasymmetry for finite splitting ∆ϵ ∼ τ > 0 (violet-dashed and blue-dotted line).In particular, the effect gets larger with larger splittings ∆ϵ, localizing the single-particle eigenstates more strongly into orbital L and R.
Finally, it is interesting to compare the above described effect leading to higher emission energy of the first emitted particle (∆E mit < 0) in the weak tunneling regime to the mechanism causing ∆E mit < 0 in the coherent two-particle simulation [Fig.2].In the latter case, the first emitted particle attains additional energy because energy is transferred from the particle remaining in the dot due to a dot-internal rearrangement process.This is physically distinct from the above described effect of asymmetric coupling, which can exist even in the absence of Coulomb interaction.

C. Impact of time-modulated barrier potential
Until here, we have kept the coupling Γ ℓ constant in time, while only the energy levels were modulated.However, a major cause of finite ∆E mit between the two emitted particles, which we have already identified from the coherent two-particle simulation, Sec.III, can be the effect of the dot potential modulation on the tunnel barrier.We can systematically analyze this within the master equation framework by letting the bare coupling strengths Γ ℓ (E, t) explicitly depend on the timedependent potential ϵ(t).Our approach is to model a situation where the energy window in which the coupling is modulated from zero to Γ ℓ , is shifted upwards during the driving by a total amount ∆E b .We therefore add to V 0 and E b in Eq. ( 8) the expression interpolating with the smooth sigmoid-like function ( 9) between an initial (V 0 , E b ) and final The influence of the varying exit barrier ( 17) on the ramp-speed dependent emission-energy difference ∆E mit is illustrated in Fig. 5.We find that if the total barrierheight shift ∆E b takes place well between the two emission events on the scale of the typical tunneling time, the barrier shift expectedly leads to a shift of ∆E mit (v r,1 ) in equal direction, sgn (∆E b ) = sgn (∆E mit ), as illustrated by the light grey solid, red and dark grey dashed lines in Fig. 5.The precise magnitude |∆E mit |, however, depends on how fast the barrier height is ramped compared to both the emission time 1/Γ ℓ and the dot potential ramp time (E b − V 0 )/v r,1 itself.This becomes even more relevant if the barrier height is ramped instead within a dot potential interval [ϵ V , ϵ b ] close to, or even containing a particle emission event.In this case, comparable time and energy scales ε(t) ∼ v r,1 ∼ ∆E b ϵ V −ϵ b v r,1 can also change the slope of ∆E mit (v r,1 ) compared to the case of constant V 0 , E b .For example, the light grey dashed line in Fig. 5 for ∆E b = ϵ V − ϵ b = 0.25U and ϵ V = V 0 − U corresponds to a case in which the onset of the opaque region of the barrier V 0 is essentially raised together with the addition energy ϵ(t) + U of the first particle once the latter has reached V 0 .The imminent particle emission is thereby continuously deferred for ϵ V < ϵ(t) < ϵ b , and v r,1 affects ∆E mit more strongly nonlinearly, via the combination of ϵ(t) itself and V 0 (ϵ(t)), E b (ϵ(t)).Apart from the opposite sign of barrier height shift, a similar nonlinearity is also seen in the dashed-dotted black line in Fig. 5 for ∆E b = −2U and ϵ V − ϵ b = 2U .

D. Apparent transition energy changes due to time-dependent charging energy
We finish our analysis by studying how the apparent equal-time transition energy difference ∆E tr between first and second emitted particle [Eq.( 16)] is affected by the dot parameters.The effects analyzed up to here are found to not significantly affect ∆E tr .This changes when the time-dependent dot potential ϵ(t) also influences the charging energy U , via the changing spatial confinement in the dot's potential landscape.A similar effect could in FIG. 6. Apparent transition energy difference ∆Etr [Eq.( 16)] as function of the ramp speed vr,1 for a potential-dependent interaction strength U (ϵ(t)) as given in Eq. ( 18).We fix ΓR = ΓL = 500Ω, ϵU = E b − U + 0.01U = 34.01U, δϵU = 0.29U .For the light grey solid line underneath the blue dashed line, we set ∆ϵ = 0.5U, τ = 0.75U , whereas ∆ϵ = τ = 0 in all other graphs.The remaining parameters equal those stated in the caption of Fig. 4.
practice also arise via a modulated single-particle splitting ∆ϵ, but we here only consider U for simplicity and concreteness.We model this in analogy to Eq. ( 17) with the smooth transition where S(ϵ, ϵ U , ϵ U + δϵ U ) (9) establishes a smooth transition between the initial (U ) and final (U + ∆U ) interaction strength within the potential interval [ϵ U , ϵ U + δϵ U ].In Fig. 6, we compare ∆E tr (v r,1 ) with constant U to cases in which the interaction strength either disappears (∆U = −U ) or doubles when ϵ is ramped up through [ϵ U , ϵ U + δϵ U ]. First note that the expected reference energy difference ∆E ref tr = U + 2 |τ | 2 + ∆ϵ 2 /4 includes both the interaction strength U for ϵ < ϵ U and the singleparticle orbital splitting due to ∆ϵ, τ .The light-grey solid and blue dashed lines in Fig. 6 show results for constant U , where we choose ∆ϵ = τ = 0 (blue) and ∆ϵ = 0.5U , τ = 0.75U .These lines clearly show that the apparent transition energy difference closely approaches this expected difference ∆E tr → ∆E ref tr .Instead, a finite ∆U ̸ = 0, quantifying the timedependent shift of the interaction energy, can affect ∆E tr if the interaction strength changes prior to or during the first emission event: A U -shift prior to emission only leads to a constant, v r,1 -independent shift to ∆E tr → ∆E ref tr + ∆U , as intuitively expected.However, if the interaction strength changes closely to the likely emission time on the scale of the typical emission rate, the effect depends on how fast ϵ, and hence U (ϵ) are shifted relatively to this emission rate.The resulting v r,1 -dependence in this case -with ϵ U + U ≈ V 0 purposely coinciding with the addition energy of the first particle reaching V 0 -is illustrated by the dark-grey dashed/black dashed-dotted line in Fig. 6.Note that we consider not only decreasing (∆U = −U ), but also increasing interaction strength (∆U = +U ): while ∆U < 0 reflects the physically more intuitive scenario in which a rising ϵ reduces the spatial confinement, one could also imagine a special potential shape in which this confinement tightens at least within a limited ϵ interval, meaning ∆U > 0. As apparent from the concrete example ∆U = +U in Fig. 6, ∆E tr − ∆E ref tr is then accordingly shifted in positive direction, albeit with smaller magnitude than the negative bending for ∆U = −U .This different magnitude occurs because a shift ∆U < 0 close to the emission time lowers the addition energy and hence further defers the emission time, thereby providing more time during which the emission energy also decreases.

VI. CONCLUSION AND OUTLOOK
This manuscript has provided a detailed analysis of the energy-and time-dependent emission of a pair of hot electrons from a driven quantum-dot potential.To capture a broad range of physical mechanisms impacting the emission process, we have employed two complementary models and ensuing approaches: a numerical two-particle simulation in a time-modulated 1d potential including Coulomb interaction, and an effective master equation description of a two-orbital quantum dot with time-dependent energy levels.Both methods have generally revealed the impact of the driven, dotinternal dynamics and of the possibly time-dependent dot exit-barrier height on the emission spectrum.A specific insight of the coherent two-particle simulation is how dot-internal energy transfer via Coulomb interaction can influence the emission process -especially in an effectively one-dimensional geometry.By contrast, the master equation analysis has demonstrated how the emission energies are impacted by dot-level modulation speed, dot-degeneracy effects typical of zero-or higherthan-one dimensional systems, and by orbital dependent dot-environment coupling asymmetries.
With the theory results obtained here, a systematic comparison with experiments is enabled.This is thanks to the largely isolated discussion of the above stated physical mechanisms that are furthermore to a good extent separately accessible/tunable in experiment.
For very fast dot level modulation, we expect the master equation approach used here to break down and we anticipate that nonadiabatic effects can become important.The treatment of these effects is however beyond the scope of the current manuscript and left open to future studies.H2020 research and innovation program under grant agreement No. 862683.Funding from the Knut and Alice Wallenberg Foundation through the Academy Fellow program (J.Sp.and J.Sc.), from the Danish National Research Foundation (J.Sc.), and the Danish Council for Independent Research Natural Sciences (J.Sc.) is also gratefully acknowledged.

APPENDIX Appendix A: Details of two-particle dynamics simulation
This appendix details how we perform the two-particle simulation outlined in Sec.II on a GPU.Generally, we need to overcome two numerical challenges.First, the two-particle interaction complicates the common approximation of unbounded systems by finite systems with an absorbing, non-Hermitian boundary Hamiltonian [62,63].As pointed out in Sec.II, our solution to this is to resort to a reflecting boundary in a large system, such that unphysical boundary reflections happen far enough away from the dot to not interfere with the emission process.Second, the time-dependent potential significantly increases the numerical cost of the conventional Lanczos method [64,65] of solving Schrödingers equation for time-independent Hamiltonians, since the Krylov subspace truncation must eventually be repeated after a few potential update.We thus instead choose a direct leap-frog [66] solver as discussed below.
More precisely, we start by describing how we map out and index the Hilbert space [Sec.A 1] and Hamiltonian [Sec.A 2].We then detail the leap-frog time evolution scheme as well as the calculation of expectation values [Sec.A 3].This appendix finishes with an overview of the simulation steps [Sec.A 4].The full code implemented in C++ and in the Nvidia CUDA progamming language is provided in Ref. [44].

State representation on triangular 2d grid
We start by considering the Hilbert space basis of occupation number states |rσ, rσ⟩ in discrete positions x r = rδx, r ∈ {0, . . ., R} with spin-z projection σ =↑, ↓ and fixed total occupation number N = 2.The full Hilbert space of two particles in 1d can be arranged in a 2d grid of projections of the full state |Φ⟩ onto these N S linearly independent basis states |i⟩ := |r i σ i , r ′ i σ ′ i ⟩.However, due to the anti-symmetry of these states under exchange of spin and position, this does not form a full square, but rather a triangular arrangement.The grid used in our simulations is illustrated in Fig. 7: The blue-shaded boxes in this grid number the i = 0, . . ., N S − 1 projections onto the basis states |i⟩.The second particle in these states with quantum number pair r ′ σ ′ is tied to the discrete horizontal grid FIG. 7. Two-dimensional triangle-grid mapping of the twoparticle state |Φ⟩ with(a)/without(b) spin degree of freedom σ in a potential landscape with R + 1 sites r.The-blue shaded boxes with number i = 0, . . ., NS − 1 indicate the corresponding projection ⟨i|Φ⟩ onto the i-basis state, given by |i⟩ = |riσi, r ′ i σi⟩ for the spinful, and by |ri, r ′ i ⟩ for the spinless case.The smaller black-crossed boxes represent empty, zero-valued grid points acting as padding elements in memory.The horizontal and vertical axes display the grid coordinates (A1) corresponding to the state numbers i.
coordinate n x = 0, 1, . . ., where (r ′ σ ′ )-index pairs with larger r ′ are located farther to the right, and spin σ ′ =↓ is always to the right of σ ′ =↑ for equal r ′ .The vertical coordinate n y = 0, 1, . . .likewise indicates the singleparticle state of the first particle with quantum numbers rσ, where larger r are further below and spin ↓ is below ↑ at equal r.Without loss of generality, we arrange the triangular grid to ensure r ≤ r ′ , i.e., in upper-right format.The grid points with rσ = r ′ σ ′ are forbidden by Pauli's principle, but are nevertheless considered part of the triangle: together with two additional rows (one per spin-projection) at the top, and one additional column to the right of the triangle (smaller black-crossed boxes in Fig. 7), they are regarded as constant, zero-weight grid points.Our numerical implementation allocates these points in memory merely to act as padding elements for sparse matrix-vector multiplications, see Appendix A 2.
Considering only the state-representing grid points, the R + 1 possible r-values and two spin-projections yield, by application of the Gauss sum rule, one grid point/vector components for each of the N S = (2R + 1)(2R + 1 + 1)/2 = (R + 1)(2R + 1) linearly independent basis states.We enumerate these non-empty points from i = 0 to i = N S − 1 by starting from the top non-empty grid row, advancing from left to right, and then by jumping to the left of the next row, skipping all empty grid points.The top-left non-empty grid point thereby corresponds to i = 0, and the bottom right to i = N S − 1.Our GPU-based numerical approach associates GPU threads directly to the index i, but not to the quantum numbers r i σ i , r ′ i σ ′ i which can be inferred from the grid-location (n y i , n x i ).Each thread hence first determines n y i , n x i from i. Carefully contemplating the above defined grid point arrangement and again making use of the Gauss sum rule, we derive and use where σ = 0/1 is synonymous to σ =↑ / ↓.The memory address of the state vector component at grid coordinate (n i y , n i x ) is given by the basis state index i -shifted by all preceding empty padding elements-to ensure sufficient data locality, see full code [44] for details.
Note that while the square root in Eq. ( A1) is in principle computationally demanding, our implementation nevertheless calculates it on-the-fly in every GPU thread without any severe performance impact.This works because our simulation mostly relies on sparse matrixvector multiplication, whose performance is in any case primarily memory-bandwidth limited.To efficiently use the available bandwidth, we compute and store the vector components only in single precision (32 bit) instead of double precision (64 bit), as our comparisons between the two have not revealed any appreciable inaccuracies for the former.Moreover, due to the absence of any Zeeman field and spin-orbit coupling, the equal-spin case effectively reduces to a spinless system.Neglecting the quantum number σ altogether for this case, we obtain a triangle grid with approximately 4 times fewer points compared to the spinful case, N S → ÑS = R(R+1) 2 , see Fig. 7(b).Moreover, we then only add a single row of top-padding empty points, and no padding column to the right.Grid and particle coordinate also become directly related, i.e., we set r i = n y i and r ′ i = n x i in Eq. (A1) without dividing by the spin factor 2.

Implementing the sparse Hamiltonian matrix
The absent Zeeman and spin-orbit coupling terms in our model mean that any Hamiltonian appearing in our simulation can always be gauged to have a fully real matrix representation with respect to the above defined occupation number basis.Moreover, the only non-diagonal elements within this matrix appear for the kinetic energy H kin [Eq.(2b)], yielding a sparse matrix containing nonzero elements only in the main diagonal and in maximally 4 sub-diagonals, i.e., 2 subdiagonals per particle/grid dimension.Exploiting this simple structure, we implement the sparse matrix-vector multiplication of any Hamiltonian with a state as a CUDA kernel acting on an array that stores the grid-ordered vector components according to the previous Sec.A 1.An actual calculation is of course only performed for output-array elements corresponding to non-empty grid points associated with a basis state index i, skipping the computation for constantly zero-valued padding elements.These padding elements are nevertheless important for any input array, as explained below.
The main-diagonal elements of the Hamiltonian matrix, i.e., those elements relating grid points i to themselves, consist of three different contributions.One stems from the onsite term 4λ of H kin [Eq.(2b)] for two particles, which is independent of time t and state i.The second contribution comes from the Coulomb potential (2c), which only requires a single division by δx|r i −r ′ i | per index i, and is hence always computed on-the-fly.This is in contrast to third contribution, given by the diagonal elements of H pot (t).The latter namely equal the sum of two single-particle potentials V (δxr i , t) + V (δxr ′ i , t) [Eq.( 1)], meaning that while several computationally expensive exponential functions are required to evaluate V (x, t), the same value for V enters many different main-diagonal elements of H pot (t).The potential V (x, t) is thus pre-computed only once for each of the R + 1 different spatial coordinates x = rδx, and stored in the GPU memory.The subsequent calculation of the main-diagonal matrix elements of H pot (t) and of the full two-particle Hamiltonian H 2P (t) then simply reads these values.
The only contribution to the sub-diagonals of the Hamiltonian matrix come from the bilinear coupling terms ∼ λ in H kin [Eq. (2b)].As such, they relate state vector components associated with different indices i ′ ̸ = i, both in horizontal (n x ) and vertical (n y ) grid direction, see full code at [44].Care needs to be taken in determining the correct fermionic exchange sign in the spinful case, and in correctly truncating the coupling at any edge of the non-empty grid.Points on these edges namely coincide with a particle at one of the ends of the 1d potential landscape, where the coupling beyond the edge no longer contributes to H kin .The above mentioned zero-value padding elements are added precisely to make any branching logic for these edge cases in the application of H kin unnecessary.This helps not only with code readability, but also with avoiding divergent GPU threads.

Time evolution of state, expectation values, and spectral density
Having defined the state space and set up the Hamiltonian matrix, we can now describe how we perform the Schrödinger time evolution (3) and obtain the expectation values for the charge and energy in the dot and in the filter region.
For a time-independent Hamiltonian, the Schrödinger equation ( 3) can be solved efficiently by truncating the evolution to the relevant Krylov subspace, such as in the often employed time-dependent Lanczos approach [64,65].However, the modulated potential V (x, t) entering our time-dependent two-particle Hamiltonian H 2P (t) [Eq.( 2)] reduces the efficiency and increases the complexity of this method, because the truncation must be repeated regularly, with possibly V (x, t)dependent truncation cutoffs.Our numerical implementation avoids these complications by instead solving the Schrödinger equation (3) directly, demanding a spatial(δx)-and temporal(δt) resolution to reliably capture the desired maximum kinetic energy E max , i.e., E max ≤ 2 π 2 2m eff δx 2 (lowest energy of two opposite-spin electrons in a box of size δx) and E max ≤ 2π/δt (using E = ω = 2π/∆t).Given kinetic energies E ∼ meV in a potential landscape of typical size L ∼ 10 4 nm, we choose a leap-frog solver suitable for the quite large number of basis states N S ∼ 10 6 required by these criteria.This method saves memory bandwidth by allocating only few temporary resources per time step, and in particular preserves the state norm without any explicit renormalization [66]: As the latter would be computationally expensive due to the necessary overflow protection for the given N S , this represents a crucial advantage.To check whether our discretization steps are actually small enough for our results to converge, we also run a few higher-resolution test-runs for comparison.
The leap-frog scheme begins by initializing the system in the ground state |GS⟩ of H 2P (0), i.e., with respect to the initial potential configuration V (x, t = 0).We achieve this by first seeding a two-particle state |Seed⟩ that has only non-zero vector components inside the dot region, 0 < rδx, r ′ δx < L d , and then by evolving this state along the imaginary time axis to project out any excited state: Each renormalization by the inverse square-root norm 1/ ⟨Φ( t)|Φ( t)⟩ is computed on the GPU using the cuBLAS Level-1 API from the standard CUDA library [44].
To prepare the actual real-time leap-frog evolution, we first divide the state into real and imaginary part.Due to measurably better performance, the latter are not stored in a block format, but instead in an interleaved fashion.This means that we do not distribute real and imaginary parts into separate grids, but rather split each state component [|Φ(t = 0)⟩] i and each empty point in the grid shown in Fig. 7 into a pair of successively stored singleprecision floating point numbers, the first corresponding to the real part and the second to the imaginary part; the empty padding elements are represented by two constant zero-values.The crucial aspect of the leap-frog method is now to introduce a half-time step between real and imaginary part.This means that the initial-state components where we have implemented a custom CUDA kernel to perform the sparse matrix-vector multiplication 1 − iH 2P (0) δt 2 |x⟩, and exploited that |GS⟩ is already normalized and has purely real components.Using the likewise purely real matrix representation of H 2P (t), the leap-frog time evolution then proceeds according to [66] Re where the CUDA kernel for the sparse matrix-vector product [H 2P (t)δt] |Φ⟩ reads the single-particle potential V (x, t) [Eq.( 1)] pre-generated for time t from GPU memory prior to evaluating the time step (A5).As already highlighted above, the key and name-giving feature of the leap-frog scheme is the initialization (A4) with half-time separation δt/2 between real and imaginary part.This stabilizes the state norm during the time evolution (A5) -and hence renders explicit renormalizations unnecessary-if that norm and the expectation value are slightly redefined [66]: Note that the interleaved memory storage of real and imaginary parts enables us to compute the norm (A6b) on the GPU with a single call of the cuBLAS Level-1 scalar-product routine, thereby preventing numerical overflow while maintaining a high runtime performance.
A particularly important expectation value for our simulation is the time-dependent spectral density ϕ 2P (t, E) [Eq.( 4)].We evaluate this density by expanding the H f kin -dependence of δ(E − H f kin ) in terms of a finite number N CP of Chebyshev polynomials of first kind [42]: with energy scale E CP chosen such that all eigenvalues of H f kin /E CP lie within the open interval (−1, 1), with the recursively defined Chebyshev vectors and with the Chebyshev expansion coefficients The Jackson kernel J n suppresses unphysical oscillations in the Chebyshev polynomial expansion of δ(E − x) due to the Dirac δ being sharp in x.Note again that the norm and vector overlaps in Eq. (A7) are defined by Eq. (A6), which, due to the employed leap-frog scheme, slightly differs from the usual definition.
To compute the sum over n in Eq. (A7) for a fixed time t on the GPU, we first obtain all ⟨Φ(t)| p n H f kin ECP |Φ(t)⟩ and buffer them in a single vector of size N CP ; all coefficients K n (E/E CP ) for all considered ratios E/E CP are stored into a matrix.Performing the sum in Eq. (A7) thereby reduces to a single call of the cuBLAS dense matrix-vector product of the coefficient matrix with the Chebyshev polynomial vector.Note, also, that we do not call a scalar-product routine for each individual polynomial ⟨Φ(t)| p n H f kin ECP |Φ(t)⟩ entering Eq. (A7).We instead employ a custom CUDA kernel for the matrixvector product H f kin |Φ(t)⟩ to recursively compute a subset of all N CP required Chebyshev vectors (limited by the GPU memory capacity) according to Eq. (A8), and first buffer these vectors as columns of a large matrix.Then, we calculate the given subset of Chebyshev polynomials with a single cuBLAS dense vector-matrix multiplication of ⟨Φ(t)| (⟨Φ(t − δt)| for the imaginary part [Eq.(A6c)]) with this matrix.We thereby benefit from cuBLASspecific optimizations of the dense vector-matrix product, which sizably reduces the time to evaluate Eq. (A7) compared to calculating all scalar products individually (∼ 10 − 20%, depending on vector size).

Overview of simulation steps
Having detailed the individual aspects of the simulation, we conclude Appendix A with an overview of the steps performed, and the system parameters used to obtain the data presented in Fig. 2  2. Compute the initial potential landscape V (x, 0) [Eq.( 1)] and initialize the system state into the ground state |GS⟩ of the initial Hamiltonian H 2P (0) [Eq.( 2)] with the scheme described in and around Eqs. (A2)-(A3).
3. Starting from |GS⟩, prepare the initial state of the leap-frog time evolution scheme according to Eq. (A4).
4. Perform a number N test GS of time evolution steps (A5) with a constant Hamiltonian H 2P (0) to let the system state settle to the actual initial simulation state |Φ(0)⟩, and confirm that H 2P (0) does not further evolve this state |Φ(0)⟩ in time apart from a global phase.For all simulations shown in Fig. 2 and Fig. 8, the system is discretized in R + 1 = 2 11 = 2048 sites.The total simulation time t end = 307.2ps is divided into numerical time steps of length δt = 2.5 × 10 −5 ps, and we compute all observables for 2 9 = 512 equidistant time points within [0, t end ], each using N CP = 2 17 = 131072 Chebyshev polynomials to generate ϕ 2P .We use N GS = 5×2 16 = 327680 imaginary time steps δ t = 1.25×10This vector evolves according to [Eq.(10)] where the elements W y,z = [W] y,z of the transition matrix W are obtained from an expansion of the master equation (10) with the help of Eqs.(C2)-(C4): as well as the following (square-root) rates: where we use As in the main text, we denote energy differences as E ij = E i − E j .
In the following, for completeness, we describe how to evaluate the Lamb shift.Note, however, that for the parameter regimes considered here, with Γ ≪ U, δE b , the Lamb shift has no visible impact on the presented results and could therefore equally be neglected.To calculate these, we assume the here relevant zerotemperature limit and model the barrier transparency as in Eqs. ( 8), (9) for energies E > µ; to enable the dot to be recharged for E < µ, we also model the barrier as reopening for decreasing energies below µ, choosing for convenience the same relative energy range δE b = |E b − V 0 | on which the barrier becomes transparent: Physically, the details of this barrier reopening is irrelevant for the emission process, apart from the mere fact that it brings the dot back to its original doubly occupied state at the end of the ϵ(t) driving cycle.However, Eq. (C14) formally enters F +ℓ just as Eq. ( 8) enters F −ℓ : It is typical for such principal value (P) integrals (C15)-(C18) to have a logarithmic divergence for infinite bandwidth |D| → ∞.We stress, however, that the shifted splitting and coupling [Eq.(C12)] only involve differences of these integrals in which all divergent terms cancel out.Furthermore, the poles E χ , E dχ avoided by principal value integration are either outside the finite support of the respective Γ ℓ (E) or fully enclosed by the integration interval; their contribution thus mostly cancels out due to the sign switch.Indeed, as confirmed by our numerical results, the effective contribution to τ , ∆ε is then ∼ Γ = ℓ Γ ℓ .Note that while we evaluate the functions (C15)-(C18) numerically as described below, the above analytical expressions already highlight a few special cases for Λ.First, we consider τ = 0, which means that the eigenvectors |±⟩ are the bare dot states |L/R⟩.We then have C ℓ + C ℓ − = C L + C L + = 0, which according to Eqs. (C12),(C13) results in τ = τ = 0. Analogously, ∆ϵ = 0 and τ > 0 implies perfectly L-R (anti-)symmetric states |C L χ | = |C R χ | and therefore ∆ε = ∆ϵ = 0 according to Eq. (C12).This, importantly, proves that Λ cannot rotate the single-particle eigenbasis if the lat-ter is given by the dot states L,R or by the corresponding bonding and antibonding states.In case of τ = 0, Eq. (C13) furthermore predicts the proportionality (∆ε − ∆ϵ) ∼ (Λ ++ − Λ −− ) ∼ (Γ L − Γ R ), meaning that the Lamb shift then scales with the environment coupling asymmetry between the two bare dot states.For the setup studied in Fig. 5 and for all curves of Fig. 6 except the light grey one, the Lamb shift therefore vanishes exactly, ∆ε = ∆ϵ, τ = τ .
We emphasize again that all system parameters entering Eqs.(C2)-(C10) can in principle be time-dependent in the way discussed in Sec.IV.We therefore need to solve the master equation (C6) numerically, which we achieve with a standard fourth-order Runge-Kutta integration scheme [44].We use N t = 2 22 = 4194304 time steps over a single period 2π/Ω to create a single current trace as in Fig. 3(a), and substitute the result into Eq.( 14) to obtain the emission times and energies (t x , E mit,x ) for x = 1, 2. This procedure is repeated for N sw = 2 9 = 512 equidistant ϵ 0 -points in the closed range [0, 35U ] to determine the energy-time relations E mit,1 (t) and E mit,2 (t) as sketched in Fig. 3(b).Note that equidistant ϵ 0 -steps do not translate into equidistant time steps for these E mit,x (t)-curves.To determine the apparent transition energy ∆E tr [Eq.( 16)], we hence first pick E mit,1 at a given time t, and then sample E mit,2 by interpolating between E mit,2 (t < ) and E mit,2 (t > ) with t < and t > being the closest times before and after t, so that t < ≤ t ≤ t > .
Finally, to obtain the Lamb shift Λ, we use Simpson's 1/3 rule with N Simp = 2 17 + 1 = 131073 equidistant points to evaluate the principal value integrals (C15)-(C18).Importantly, if any resonance lies within the integration interval, we implement the principal-value integration to achieve stable convergence by placing these points symmetrically around this resonance.Since the integrals need to be obtained for a large range of E χ or E dχ in many different driving protocols, we reduce the number of required computations by precomputing the integrals for N Int = 2 23 = 8388608 equidistant energy points in a range [−2.5A, 2.5A].The matrix elements Eq. (C13) are then obtained by linearly interpolating between these points.A further reduction of the integral computation time is achieved by likewise linearly interpolating between N Γ = 2 18 = 262144 precomputed points of the Γ(E) functions.We have checked for convergence by individually increasing the above stated number of points and comparing the results.

FIG. 1 .
FIG. 1.(a) Sketch of a tunable-barrier quantum-dot setup.Gates on top of a conductor time-dependently modulate the potential landscape [Eq.(1)] to eject pairs of electrons.The positions x d,L , x d,C , and x d,R indicate the boundaries and center of the quantum-dot region, x f,L and x f,R set the regions, where a filter would measure time and energy at which electrons propagate through it.(b) Sketch of the emission process: the potential is driven such that the first electron can tunnel out, until only one electron remains in the quantum dot.Further driving allows also the second electron to tunnel out until the potential well has basically disappeared.

.
The left Gaussian peak defines the left quantum dot 'wall' of height V d,L at the reference coordinate x d,L ≡ 0 with width σ d,L , and the right peak represents the exit barrier of height V d,R and width σ d,R at x d,R .The maximum of the resulting emission barrier with respect to the potential in the conductor V 0 at the initial time t 0 defines a characteristic energy scale of the potential and is indicated by δE b in Fig.1.The central, inverted Gaussian with constant width σ d,C but time-dependent depth V d,C (t) at x d,C > x d,L = 0 establishes a tunable potential well.At initial time t 0 ≡ 0, we choose V d,C (0) > 0 such that a dot forms inside this potential dip, and confines two particles within the typical spatial dot range1 [0, L d = 5σ d,C ].The dip potential V d,C (t)is then raised until the particles are emitted into the flat channel right to the exit barrier, where we linearize

FIG. 2 .
FIG. 2. Energy-and time dependent emission spectroscopy.Upper row: Vanishing Coulomb interaction, κ = 0, for equal spins, ↑, ↑, in panels (a) and (b) as well as for degenerate spin states ↑, ↓ in panels (e) and (f).Bottom row: Finite Coulomb interaction κ = 0.125, for equal spin states ↑, ↑ in (c) and (d) and for degenerate spin states ↑, ↓ in (g) and (h).We show the number of electrons in the dot N d , and the average kinetic energy and particle density in the filter region, E f kin and N f as function of time (left columns, a,c,e,g); the right columns (b,d,f,h) show the time-resolved energy-density, ϕ2P(t, E), where the discrete contributions stem from the discrete energies of the finite filter region.The typical energy scale δE b is the numerically determined difference between the initial (t = 0) exit barrier height and the conductance band potential bottom V0 [Fig.1(a)].The red stars mark the time-energy peaks of the average E f kin , the turquoise dashed and dotted lines, respectively, indicate the ramp speed vr and the lowest level with sizable density once the wave package is fully located in the filter region.We set t end = 307.2ps, m eff = 0.016me, L = 12000 nm, x d,L = 0, x d,C = 115 nm, x d,R = 255 nm, L d = 300 nm, x f,L = 1000 nm, x f,R = 2350 nm, and a potential landscape V (x, t) [Eq.(1)] with vr ≈ 0.39 meV/ps, σ d,L = 90 nm, σ d,C = 60 nm, σ d,R = 50 nm, V d,L = 2.2V0, V d,C (0) = 0.8V0, V d,R = 0.006V0, δE b ≈ 0.0057V0, where V0 = 150 meV.With these parameters, the initial single-particle energy splitting of electrons with equal spin in the lowest energy states in the potential well ((a),(b),(e),(f)) is 4.69meV ≈ 5.5δE b and the initial charging energy in the interacting cases ((e),(f),(g),(h)) is ⟨H Coul ⟩(t = 0) = 8.35meV ≈ 9.83δE b .Appendix A 4 provides all parameters related to the numerical implementation.

FIG. 3 .
FIG. 3. (a) Sketch of the discrete single-particle energy levels ϵ L/R , their coupling τ , as well as two-particle interaction U in the quantum dot.Electrons can be emitted from the dot far above the Fermi level µ due to the time-dependent driving Vt and energy-dependent coupling Γ(E) to the reservoir.(b) Example for the dot occupation number and the emitted charge and energy currents as function of time due to the modulated potential.(c) Example for emission energies of the two electrons as function of emission time.Stars indicate energy-time pairs of one driving protocol on a line extrapolated by the time-dependent addition energies, see Sec.V D for definitions of ∆Emit and ∆Etr.

FIG. 4 .
FIG. 4.Emission energy differences ∆Emit [Eq.(15)] as function of the ramp speed vr,1.(a) shows the case of zero level-coupling τ = 0.In (b) we set for all lines ΓR = 500Ω.The black dashed-dotted lines in panels (a,b) are identical for reference.The time-independent parameters in both panels are U = 1000Γ, A = 34U , V0 = 35U , E b = 35.05U .Appendix C states all parameters relevant for the numerical implementation.

FIG. 8 .
FIG.8.Energy-and time dependent emission spectroscopy for effective mass m eff = 0.016me as in main text (a,b) compared to the corresponding quantities for m eff = 0.067me as in GaAs (c,d).We assume two equal ↑-spins, κ = 0.125, and all other parameters as stated in the caption of Fig.2and at the end of Appendix A 4.