Exact two-body quantum dynamics of an electron-hole pair in semiconductor coupled quantum wells: a time-dependent approach

We simulate the time-dependent coherent dynamics of a spatially indirect exciton (an electron-hole pair with the two particles confined in different layers) in a GaAs coupled quantum well system. We use a unitary wave-packet propagation method taking into account in full the four degrees of freedom of the two particles in a two-dimensional system, including both the long-range Coulomb attraction and arbitrary two-dimensional electrostatic potentials affecting the electron and/or the hole separately. The method has been implemented for massively parallel architectures to cope with the huge numerical problem, showing good scaling properties and allowing evolution for tens of picoseconds. We have investigated both transient time phenomena and asymptotic time transmission and reflection coefficients for potential profiles consisting of i) extended barriers and wells and ii) a single-slit geometry. We found clear signatures of the internal two-body dynamics, with transient phenomena in the picosecond time-scale which might be revealed by optical spectroscopy. Exact results have been compared with mean-field approaches which, neglecting dynamical correlations by construction, turn out to be inadequate to describe the electron-hole pair evolution in realistic experimental conditions.


I. INTRODUCTION
Electronic bilayer systems have exposed a huge amount of new physics driven by inter-layer Coulomb interactions. A few examples are fractional quantum Hall effect states in semiconductor 1,2 and graphene bilayers, 3 phase transitions in Quantum Hall ferromagnets, 4 complex Wigner crystal ordering 5 . Electron-hole bilayers gained importance in their own right. Excitons -bound electronhole quasi-particles with a bosonic character-have long being predicted to undergo quantum condensation. 6,7 Recently, signatures of condensation have been found in systems of spatially indirect excitons (IXs), electron-hole pairs optically excited in semiconductor coupled quantum well systems, with the two charges confined in different layers by a static electric field. [8][9][10] In a different perspective, IXs are at the heart of a new class of opto-electronic devices, made possible thanks to the small electron-hole overlap which extends their intrinsic lifetime from nanoseconds 11 to microseconds. 8,12 Indeed, although IXs are neutral excitations, they carry a large finite electric dipole which can be used to drive the evolution of IXs in the coupled quantum well (CQW) planes by electric field gradients generated, e.g., by metallic gates. Finally, since a bias normal to the QW-plane can control the overlap of the pair, IX recombination can be induced at arbitrary time, thereby 'measuring' the result of the evolution. IX gases have been exploited to demonstrate several functionalities, such as fast data storage, 13,14 acceleration with electrostatic ramps 12,15 and interdigital devices, 16 field effect transistors. 17 Furthermore, trapping of single IXs has been recently demonstrated, 18 opening the way to single IX electronics. This requires the development of theoret-ical concepts to describe the evolution of IX wave-packets in complex electrostatic fields.
Scattering of a composite quantum objects with internal degrees of freedom 19 (DoFs) like an IX in the presence of an electrostatic field gradient, 20 is an important topic in its own, with applications in molecular and nuclear physics (see, e.g., Refs. 21-24, and references therein). Indeed, in the presence of a scattering potential, energy can be transferred between the centerof-mass (CM) kinetic energy and internal excitations, which may strongly influence the transmission and reflection probabilities. However, due to the difficulty to evolve the quantum equations of motion for several DoFs with open boundaries, exact calculations are often limited to idealized situations, such as collinear scattering, purely one-dimensional (1D) systems, and/or very simple potential profiles. 19,22,[25][26][27][28] Indeed, the numerical approach scales exponentially with the number of DoFs. For realistic situations, such as the 3D problem of colliding molecules with complex inter-molecular interactions, mean-field methods have been applied. 29 IXs in CQWs are an important system from this point of view, since in principle their evolution could be probed by accurate time-dependent optical means, not only at asymptotic times, but also during the scattering event. Recently, we used an idealized 1D model to study the evolution of IX wave-packets under the action of external fields, taking explicitly into account the internal structure of the electron-hole pair. 30 Our model allowed to investigate different regimes/potential profiles where inter-particle Coulomb correlations may lead to internal excitations or even dissociation, as a result of scattering. However, such 1D calculations are too simplistic to be applied to realistic CQW quasi-two-dimensional (2D) systems, where electrons and holes evolve in a complex 3D structure.
In this paper we report time-dependent simulations of the coherent dynamics of a single IX wave-packet in a semiconductor CQW structure with complex in-plane electrostatic potentials, taking into account the 3D structure of realistic devices through a 2D+1D effective model. We used a unitary wave-packet propagation which includes the long-range electron-hole Coulomb interaction and arbitrary scattering potentials acting on the electron and the hole separately. The large numerical problem to simulate exactly the present four-DoF system has been tackled by the Fourier split-step method implemented on a massively parallel architecture. This allowed propagation for tens of picoseconds in typical potential landscapes. We have investigated potential profiles consisting of extended barriers/wells, and a single slit geometry, finding genuine signatures of the two-body dynamics in realistic experimental conditions, with transient phenomena in the picosecond time-scale. Our results could be directly compared with time resolved optical spectroscopy. We have compared the unitary evolution results with a mean-field approach at different levels of approximation, the so-called rigid exciton (RIX) model and the time-dependent Hartree (TDH) method. The comparison shows that a mean-field approach is in general inadequate to describe the electron-hole pair evolution in realistic samples, thereby showing that IX dynamics in CQW might be a particularly interesting system to investigate correlation effects and to test theoretical modeling.
In Sec. II we define our Hamiltonian description of an IX in a typical semiconductor CQW (II A) and we provide details on the full (II B) and mean-field (II C) wave-packet propagation methods. Initial conditions are discussed in Sec. II D. In Sec. III we calculate the freeexciton properties of our model (III A), while numerical details of wave-packet propagation are discussed in Sec. III B. Results are summarized for scattering potentials which are weak (Sec. III C) or strong (in Sec. III D) with respect to internal excitations, and for a single-slit geometry (Sec. III E). Section IV discusses, in particular, the predictivity of the different approaches. A formal derivation of the mean-field propagation scheme is provided in the Supplemental Material. 46

II. THEORETICAL APPROACH
A. The electron-hole Hamiltonian Our reference system is sketched in Fig. 1(a). A symmetric GaAs CQW structure grown along z is embedded in a Al x Ga 1−x As matrix. A vertical electric field F z along the growth direction separates electrons and holes in different layers. 31 While typical CQW confinement energies amount from tens to hundreds of meV, in-plane (xy) potential landscapes generated by metallic gates, as well as kinetic energies which can be impressed upon IXs, are in the meV range. Therefore, we factorize the in-plane and vertical component of the IX 3D wave function (here and throughout we use wave function in place of envelope function) as where r i and z i are the 2D xy coordinate and z coordinate of the two particles, respectively. ζ e (z e ) and ζ h (z h ) are calculated explicitly for a given structure and field. ζ e is calculated from the 1D effective-mass equation with the Hamiltonian where W e (z e ) and m z e (z e ) are the band edge and the material-dependent effective mass, respectively, of conduction electrons in the CQW structure. No spindependent term is considered. Indeed, the validity of Eqs. (1) and (3) is limited to samples with narrow CQWs which are the typical heterostructures used in the experiments we are addressing to [12][13][14][15][16][17][18] . Here, in-plane Coulomb binding energy and scattering potentials are in the few meV range, while vertical confinement energies is at least one order of magnitude larger.
Equation (2) with the Hamiltonian (3) is called a Ben-Daniel Duke problem. 32 E ze and ζ e (z e ) are obtained from (2) on a homogeneous real-space grid by a finite difference approach, taking into account the material-dependent effective mass. E z h and ζ h are computed similarly, with parameters appropriate to the valence band electrons. In this case m z h is the effective mass given by the heavyhole diagonal mass tensor, related to the Luttinger's parameters, 33 γ 1 and γ 2 , by An example of these calculations is reported in Fig. 1(b) showing that the two carriers are well localized in either wells by the external bias. This justifies the separability of the electron and hole wave functions in the z direction assumed in Eq. (1).
To simulate the electron-hole quantum dynamics we need to include accurately the mutual interaction. While in a ideal 2D electron system this is described by the bare r −1 Coulomb potential, in a typical CQW heterostructure carriers are delocalized in the wells over a length which is comparable to the carrier separation, given by the effective Bohr radius. Therefore, the effective interaction is modified at short range, up to a distance comparable to the well width. To account for this effect, we consider electrons and holes in the ground state of the confinement potential (a reasonable assumption due to the large energy gaps in the growth direction) and we consider the effective 2D interaction U C as the mean value of the Coulomb interaction over the wave functions in the growth direction, Here, r is the electron-hole in-plane distance and r the relative permittivity of the well material. Therefore, the full 3D problem has been mapped into an effective 2D model. The two-body wave function Ψ is then propagated in time according to the 2D Hamiltonian where is the free IX Hamiltonian. Here, m h is the in-plane component of the strongly anisotropic mass tensor of GaAs, It is actually numerically convenient and more transparent to work in the CM and relative coordinate system where M = m e + m h is the in-plane exciton mass and m = [m −1 e + (m h ) −1 ] −1 is the in-plane reduced effective mass. In this representation, the free IX Hamiltonian separates as and the free in-plane wave function can be factorized as As we shall discuss later, the CM energy is typically of the order of tenths of meV, much less than the internal excitation energy, which is of the order of several meV. Therefore, the grid parameters related to CM and relative dynamics are quite different and can be optimized separately. On the contrary, using electron and hole coordinates requires to use (almost) the same grid, and the same accuracy would be reached at a greater computational cost. 30 In this representation, the inclusion of one-body external potentials, U ext = U e (r e ) + U h (r h ) = U ext (R, r) in the full Hamiltonian couples the CM and relative coordinates, and removes the separability of the wave function, Eq. (14). Therefore, to propagate Ψ(R, r; t) we need to deal with a four DoFs propagation scheme. We neglect the interaction of IXs with environmental degrees of freedom, like phonons of the semiconductor lattice. In typical experiments E CM is in the tenths of meV range, well below the optical phonon energy in GaAs, while low temperature strongly suppresses acoustic phonon population. In our simulation, the scattering time is in the order of tens of picoseconds, well below the LA-phonon assisted relaxation of IXs in these devices, which is in the nanoseconds range. 8

B. Full numerical propagation
The quantum propagation of the IX is obtained through the numerical application of the evolution operator, U(t + ∆ t ; t), between two consecutive times t and t + ∆ t as Our numerical solution relies on the Fourier split step (FSS) approach. 34,35 This unitary method, numerically exact as ∆ t → 0, is based on the Suzuki-Trotter factorization 36 of the evolution operator in the product of two exponential operators, containing the kinetic or the potential operators, respectively, each diagonal either in direct or in reciprocal space (see Appendix in Ref. 30 for details): where T = P 2 /(2M ) + p 2 /(2m) and U tot = U C + U ext are the total kinetic and potential energy operators of the system (P and p are CM and relative in plane linear momenta, respectively). Hence, at each time step the IX wave function must be switched from position to momentum representation, and vice versa, through Fourier transformation, F, according to The use of the Fast Fourier Transform (FFT) algorithm, therefore, results in high computational efficiency with respect to other methods, particularly those based on finite difference discretized Hamiltonian, as the Crank-Nicolson method. The coupling between CM and relative DoFs introduced by the one-body external potentials U ext (R, r) requires the numerically full propagation of all four DoFs of the IX for a sufficiently long time for the scattering process to conclude on a sufficiently dense and extended real-space grid. For the energy scales into play here, this turns out to be a demanding numerical task, both to store the complex-valued IX wave function in memory, and to numerically compute the application of the evolution operators to it. A code exploiting massive parallelization has been developed to cope with these issues. The four-dimensional domain is discretized in a grid of about 4295 × 10 6 points. Typically, the computation of a single 40 fs time step takes 10.5 seconds, and one Gb of memory per core is used on a 256-core run (with 16 cores in two 2.4 GHz Haswell Xeon processors per sharedmemory node), corresponding to a speedup of about 43 with respect to a serial run.

C. Mean-field propagation
In semiconductor physics, when dealing with excitons in weak potentials, it is often justified to apply the socalled rigid exciton model (RIX), 37 which consists in the assumption that the IX is frozen into its relative motion ground state, φ 0 (r); hence, only the quantum evolution of the CM component, χ(R), is taken into account. The internal DoFs are integrated out, leading to an effective potential Therefore, the IX moves as a rigid object with coordinates R, its dynamics being determined by the external potential averaged on the relative-motion ground state. For example, IX wave function localization in weak traps can be calculated in this way. 38 The CM effective evolution operator is then applied to the CM wave function, χ(R, t), and the numerical evolution can be obtained again by the FSS method. Obviously, the RIX model requires a much lower computational effort with respect to full propagation, since only the two DoFs wave function χ(R) needs to be propagated, which can be readily obtained on a standard personal computer. Clearly, every effect of the internal dynamics is neglected by the RIX approximation.
The RIX model is the lowest order example of a more general mean-field strategy, also known as the time dependent Hartree (TDH) method in atomic and molecular scattering. 29,39 Within this approach, at each time t the global wave function of the composite object is assumed to be factorized into a CM and a relative wave function, The evolution of χ(R; t) and φ(r; t) is determined by an effective potential representing the expectation value of the external potential on the relative and CM wave function, respectively, at that specific time t, i.e. (see Ref. 39, and Suppl. Mat. 46 ) and with This is a mean-field model, since the evolution of the CM wave function is determined by the average field generated by the relative wave function, and vice versa, at each time step. Again, the numerical propagation is performed using the FSS method. The FSS algorithm needs to be applied twice at each time step, independently of the CM and relative wave functions, the two evolutions being coupled through the effective potentials, Eqs. (54), and (55). Clearly, the RIX model consists in assuming a rigid φ(r; t) = φ 0 (r), and only the χ(R; t) component needs to be propagated.
Note that, even for a stationary external potential U ext (R, r), the mean-field propagation method requires to evolve the two components of the wave function under time-dependent potentials. This increases substantially the computational cost of the simulation with respect to the RIX model, since the propagator needs to be calculated at each time step rather than only once. However, the TDH approach is still far less demanding than the full evolution.

D. Initial state
In order to start a time dependent simulation, we need to choose a proper initial state. In typical CQW systems, an IX thermalizes and relaxes to the ground state φ 0 (r) of the free exciton Hamiltonian H 0 within nanoseconds from photogeneration, due to active scattering mechanisms (acoustic phonons). We do not consider in our simulations this transient, which is short compared to the IX photo-recombination lifetime. 8 Therefore we initialize the IX wave function as The CM wave function is chosen as the minimum uncertainty wave-packet centered at the initial CM position (X 0 , Y 0 ), having widths σ X and σ Y , and propagating with an average CM wave vector where E CM is the most probable CM kinetic energy at t = 0 and θ identifies the initial propagation direction with respect to a properly defined normal incidence direction. In the following simulation we take, as initial dispersion, σ X = σ Y = 80 nm unless otherwise specified. This value is of the same order of magnitude of the confinement length of IX traps 18 . This is a reasonable compromise between a sufficiently narrow momentum distribution (hence a delocalized exciton) and a spatially localized IX. The CM momentum can be controlled, e.g., through the application of acceleration ramps. Localization of the initial state can be controlled by electrostatic traps. See Refs. 12-18, and 49 for typical parameters, comparable to those used in our simulations. Note that, while the CM energy clearly affects transmission and reflection coefficients, due to the linearity of the equations the momentum distribution does not affect much the dynamics, as long as the momentum dispersion is not too broad.
We remark that, even if the initial state is factorized, in the full propagation the wave function is correlated, evolving under the influence of both the electron-hole interaction and the external potential, with no specific a priori decomposition. On the contrary, factorization for the IX wave function is assumed in the RIX and TDH approximations at any intermediate time t.

III. RESULTS
In the following we investigate the IX dynamics in a 8nm/4nm/8nm GaAs/Al 0.33 Ga 0.67 As/GaAs CQW system. 7,8,40 Band parameters are indicated in Tab. I. An homogeneous (i.e. constant along the xy planes) electric field F = 1mV/nm, generated by the application of a bias voltage between the top and the back gates of the sample, is assumed (see Fig. 1

).
A. Free IX: effective Coulomb interaction and relative motion eigenvalue problem In Fig. 1(b) we show the square modulus of the calculated electron and hole wave function components (ζ e , ζ e ), together with the band profile along the growth axis of the heterostructure. Clearly, the electric field localizes the electron and the hole in different layers.
ζ e , ζ h are used to compute the effective electron-hole interaction U C (r) (see Eq. (5)) along the CQW planes, shown in Fig. 2(a). At large electron-hole distances U C amounts to the bare Coulomb interaction. At small distances, the Coulomb divergence is removed due to the separation of the electron and the hole in different layers. Bound electron-hole states and energy levels E n are calculated from H rel with a finite difference approach on a 2D uniform square grid. The grid density is the same as for the wave packet evolution (see Sec. III B). The lowest bound states are shown in Fig. 2 These parameters satisfy the Nyquist criterion 44 for the spatial frequency sampling far beyond the considered energy ranges. Another criterion in order to apply the FSS method with no ambiguity requires that the phase exponent U ext ∆ t /(2π ) 1. If this condition is not fulfilled, the external potential strength is not well defined, because the class of mod (2π /∆ t )−defined potentials has more than one element. In this case the method gives the same results with different potentials, which is physically inconsistent.
In order to estimate the IX localization during the time dependent simulations and, in the asymptotic times, the transmission and reflection probabilities, we define three regions of the CM space: (i) the reflection region A, i.e., the subspace with vanishing external potential, U e = U h = 0, where the IX wave function is initially localized; (ii) the potential region B; (iii) the transmission region C, i.e., the subspace with vanishing external potential which can be reached from A only by crossing the potential region. Note that these regions refer to the CM DoFs only, while in the four DoF model the external potential U ext depends also on the relative coordinate r. Thus, in order to compare consistently the different methods, we define the above regions according to the potential of Eq. (18). Specifically, we choose as boundaries of the region B the Y positions where the external effective potential drops to 5% of its maximum value. The related coefficients are then defined as integrals, over the corresponding regions • of the CM-part of the wave function square modulus, |χ(R; t)| 2 , for the mean field approximations; • of the CM marginal probability, for the full propagation.
These integrals are then normalized to the whole domain, so that the A,B,C-coefficients are defined in the interval [0, 1]. An example of the three regions is showed in Fig. 3. Clearly, the three coefficients evolve in time. Coefficients A and C at asymptotic times set to a constant which are the reflection and transmission probabilities, respectively. Below we shall investigate the time-dependent dynamics of a wave-packet, prepared as described in Sec. II D, scattering against two classes of potentials, a uniform, infinitely long barrier or well, and a barrier with a slit, which mimics the potential generated by a split gate, as typically realized in 2D heterostructures. A few comments are in order: -The bias and/or gating potential drops are much smaller than the confinement energies in the quantum wells of the structure. Therefore, wave functions ζ e , ζ h are not distorted by local variations of the external potential, which is thus uninfluential on the vertical localization of the electron and the hole. In other words, the effective interaction U C can be considered independent of space (R) and time; -In CQW systems U e and U h are in general different, with opposite sign. A metallic gate on top of the structure, for example, generates an electrostatic potential which is opposite for electrons and holes, and it is slightly different in strength between the two layers, due to the different distance from the gate. A simple capacitor model 30 shows that the difference is typically of a few meV and comparable to the generated in-plane voltage drop. 45 -We shall investigate external potentials which are short ranged and vanish exactly outside the scattering region. According to the energy conservation law, excitations to higher internal IX levels, up to the dissociation threshold, are allowed only inside the scattering region B; in the asymptotic regions (A, C), the CM energy (which in our simulations is always much smaller than the lowest internal excitation threshold, see Tab. II) is not sufficient to excite the internal dynamics, and the IX can only be transmitted/reflected in the internal ground state.
Two regimes can be identified, with the strength of the scattering potential being weak or strong with respect to the internal motion excitation energies. Below we shall investigate separately these two regimes for uniform, infinitely long well/barriers potentials, chosen separately for the two particles, but with a common width L y , C. Weak external potential In Fig. 3 we show a simulation for a weak potential well for the hole, U h,0 = −1.0 meV, U e,0 = 0, and L y = 40 nm. Note that the external potential strength is smaller than the lowest internal excitation energy, E 1 − E 0 = 2.35 meV (see Tab. II). The IX is initialized with a kinetic energy E CM = 0.2 meV with normal incidence (θ = 0) to the well. Figure 3(a) shows that the IX is completely transmitted as a bound state in the internal ground state, as in a single-particle scattering. This is confirmed in Fig. 3(b) where we plot the time evolution of the A,B,Ccoefficients, showing that scattering is over in ∼ 20 ps, and the wave-packet is completely transmitted. Furthermore, the full propagation and the RIX and TDH propagations give indistinguishable results (we did not plot the TDH calculation for clarity). This proves that in this regime i) the wave function can be factorized into the product Ψ(R, r; t) = χ(R; t)φ(r; t) and ii) it remains in the ground state φ(r; t) = φ 0 (r)e −iE0t/ of the internal DoFs, the fundamental assumption in the RIX approximation, during the whole propagation.

D. Strong external potential
We next investigate scattering of an IX against external potentials with an energy scale comparable to the IX internal excitations.

Electron well
We first consider scattering of a IX with an external potential consisting of a square well applied to the electron, U e,0 = −3.0 meV, U h,0 = 0, and L y = 40 nm (see Fig. 4(a)). The initial CM kinetic energy of the IX is set to E CM = 0.4 meV. Now the external potential intensity is slightly below the dissociation energy, −E 0 = 3.63 meV, but it is sufficient to excite the IX to higher energy internal states. The dissociation phenomena, not possible here due to energy conservation, has been analyzed elsewhere for a simpler 1D geometry 30 .
To highlight the role of internal excitations, in Fig. 5 we show the projections on the first internal eigenstates j = 1s, 2p y , 2s at t = 10 ps. The p x -state projection is not shown, since it cannot be excited due to potential symmetry reasons. We see that the excitation of the internal DoF takes place especially at the edge of the well, where the change in potential energy is abrupt. ρ j CM (R; t), for j = 1s, vanish as t 20 ps, i.e. when the scattering process is almost concluded.
In Fig. 4(b-d) we show the A,B,C-coefficients at normal incidence θ = 0, and θ = π/6 and θ = π/4 incidence. In all these cases, scattering is over in ∼ 15 ÷ 20 ps, the scattering time being larger for larger incident angles. The transmitted wave-packet and the amplitude in the potential region are smaller as the angle increases. This is in agreement with the lower momentum of the exciton in the direction normal to the potential well, if scattering is not at a resonance energy.
In this regime, RIX or TDH substantially overestimate the transmission obtained from full propagation, while RIX and TDH give very similar results between each other. Note however, that the trend with incident angle is similar for the three methods and that in the potential region C all methods give very similar results. Therefore, even though the CM localization is similar during scattering, the internal DoFs have a strong effect on the transmission and reflection at asymptotic times.
We also verified by explicit simulations that transmission with θ > 0 coincides with that obtained at normal incidence but with a kinetic energy which corresponds to the normal component of CM wave vector, which for the present case is 0.3 meV at θ = π/6, and 0.2 meV at θ = π/4. This holds true both in the full propagation and in the mean-field methods. This is clearly to be expected in   the RIX approximation, which is effectively a one-particle problem, since the evolving CM wave function is separable into an X and a Y part in a translationally invariant potential. For full and TDH propagations, where the Coulomb interaction couples x and y directions, note that for an external potential which is invariant in one planar direction, say, with respect to x e , x h , the total potential is invariant with respect to X, U C (x, y)+U ext (Y, y). Therefore, the X Fourier components of the full wave function is not scattered by the potential, and separates as Ψ(R, r) = e iK X X Ψ ⊥ (Y, x, y). Accordingly, the dynamics corresponds to the free-particle one in the X direction.

Hole well
We next consider the equivalent scattering problem, but with the potential well applied to the hole, with U h,0 = −3.0 meV, U e,0 = 0. Selected results are shown in Fig. 6(b-d). The kinetic energy and scattering angles are identical as in Fig. 4(a). Due to the different effective masses of electrons and holes, however, this problem is not equivalent to the previous one.
Several differences with respect to the electron case (see Fig. 4) can be recognized. First, the component in the potential well during the scattering is larger and, accordingly, scattering times are longer. This is expected, due to the heavier mass of the hole. However, the localization in the potential well (B-coefficient) is largely underestimated by the mean-field calculations, contrary to the previous case. Second, when the normal component of the exciton momentum is changed, there is no definite trend between the full calculation and mean-field calculations: at normal incidence the transmission coefficient is lower for the full propagation with respect to the RIX and TDH approaches. At θ = π/6, Fig. 6(c), the C-coefficient is lower in the full propagation only until t ≈ 20 ps, but the asymptotic value is larger for the full calculation than in the mean-field approximations. Finally, at θ = π/4, Fig. 6(d), the full C-coefficient is always larger than the RIX value, both during scattering and at asymptotic times. This should be ascribed to the activation of a resonant transmission channel. This coupling is able to excite the IX into a superposition of higher internal states and it breaks the separability into CM and relative motion parts. This is why even the TDH method is not able to reproduce the full results, and merely mimics the RIX approximation.
A peculiar behavior is shown in Fig. 6(b) which exhibits a plateau in the C-coefficient during the scattering within the full dynamics. This indicates a non uniform transmission of the IX wave function. In a semiclassical picture, the hole is trapped into the well while the electron is partially transmitted. Due to Coulomb interaction, however, the transmitted electron inverts its motion, the IX bounds again and the pair is finally transmitted. Therefore, the transmitted wave-packet splits into an advanced and a delayed IX. Clearly, this requires the excitation of the internal modes, and this behavior is not reproduced by the RIX model. 30 To confirm this interpretation, we calculate the classical internal oscillation period, (see Suppl. Mat. 46 ) and estimate it in τ ≈ 5.5 ps, which is indeed comparable to the time duration of the transmission plateau in Fig. 6(b).

Symmetric potential barrier/well
For completeness, we next discuss a symmetric potential with opposite sign for the two particles, U h,0 = 3.0 meV and U e,0 = −3.0 meV. The potential profile is sketched in Fig. 7(a). This is a somehow special situation, because the average potential is zero.
We consider θ = 0 and θ = π/4 at the CM kinetic energy of E CM = 0.4 meV. For both incident angles, the transmission coefficient for the full propagation is small,  due to the repulsive barrier felt by the hole. Most interestingly, it is substantially smaller than the one obtained from the mean-field methods, for which the average of the external potential over the internal DoFs makes the hole barrier smooth, thus favoring transmission. Note also that for θ = π/4 the TDH result deviates substantially from the RIX calculation. Nevertheless, it is still far from the full calculation.
To show that this is not an accidental situation, we have calculated the transmission coefficient as a function of the width of the external potential (Fig. 7(d)). Note that the behavior of the transmission coefficient is not monotonous with the the width of the external potential, L y , as a possible consequence of resonant transmission for the CM DoFs alone. In particular, both for the full and the mean-field propagation, the transmission coefficient is larger for L y = 28 nm and L y = 52 nm than for L y = 40 nm. Still, the transmission is systematically overestimated by the mean-field approach.

E. Single slit potentials
We finally investigate evolution through a single slit potential. This consists of an aperture of width ∆ in an otherwise infinitely long barrier/well potential similar to those investigated in the previous sections. In Figs. 8 we summarize results for a slit with ∆ = 240 nm and a slightly asymmetric electron/hole potential U e = 6 meV, U h = −3 meV, and L y = 20 nm. These values exclude tunneling through the barrier, while satisfying the condition U ext ∆ t /(2π ) 1. The IX is initialized with a CM kinetic energy E CM = 0.4 meV, moving towards the mid-point of the slit with normal incidence, and an initial width of the CM minimum uncertainty wave packet σ X,Y = 160 nm. These parameters produce several diffraction lobes in the transmitted wave-packet in the RIX approximation, shown in 8(a). Moreover they avoid the IX wave from spreading too much before reaching the slit. 47 Figure 8(b) shows a snapshot of the CM marginal probability of Eq. (28) at t = 36 ps. A comparison with the equivalent RIX calculation in panel (a) helps to identify several important features. First, while the reflected part of the CM wave-packet is very similar in the two calculations, the diffraction lobes are almost suppressed in the full calculation. Second and most interesting, part of the CM wave-packet propagates as edge states along the barrier, far from the aperture. In semiclassical terms, this corresponds to a IX, with the hole trapped inside the well and the electron trapped on either side of the barrier by the electron-hole attraction. Therefore, this is a genuine correlation effect which cannot be reproduced by a mean-field approach, and indeed it is completely absent in Fig. 7(a).
The evolution of the A,B,C-coefficients are shown in Fig. 8(d) for full propagation and the RIX model (the TDH approximation almost coincides with the RIX one and it is not shown). Note that here the potential region for integration has been extended to ten times the external potential range (200 nm, see Fig. 8(a)) to measure accurately the asymptotic coefficients, since part of the wave function remains trapped at the edges of the potential, as we discussed above. These coefficients show a qualitative agreement between full and mean-field methods. This is because in the region where the wave function is large (the slit) the potential is vanishing. Therefore, the average external potential contribution is weak, and we are in a regime similar to Sec. III C. However, at asymptotic times, the full calculation shows a trans-mission coefficient which is smaller than in the RIX calculation, the difference being the fraction of the wave function propagating along the edges of the external potential. Therefore, even if in a weak potential regime, correlation effects are exposed in the full calculation but not in the mean-field propagation.
In Fig. 8(c) we plot |Ψ(R, r = 0)| 2 for the full calculation. On the one hand, there is very little difference with respect to the CM marginal probability of panel (b), indicating that the r = 0 contribution is by far the largest in the relative coordinate average which provides the CM marginal probability. On the other hand, this is proportional to the IX optical recombination probability and shows that an optical luminescence experiment with sub-µm resolution would be able to probe the wave function with accuracy. 48 We have repeated similar calculations, but with the electron and the hole potential exchanged, U e = −3.0 meV, U h = +6.0 meV. We observed, in this case, no qualitative difference with respect to the previous one. While full and mean-field calculations agree overall, there is a substantial part of the CM wave function which propagates as a bound IX along the potential edges, which is not captured by the mean-field calculation.
In Fig. 8 we show a somehow special situation, where transmission and reflection probabilities are almost equal. Therefore, we show for completeness in Fig. 9 the evolution of the A,B,C-coefficients when the slit potential is 'open', i.e., most part of the CM wave-packet is transmitted, or 'closed', i.e., the CM wave-packet is almost fully reflected. Results are shown for the same potentials used in Fig. 8. Again, we repeated the calculation interchanging the particle potential, finding no quantitative difference with the original one. Interestingly, in all cases a similar fraction of the CM wave function propagates along the edges of the potential barrier and well. This phenomenon seems thus to be mostly related to the Bohr radius of the exciton, rather than the slit aperture: a larger (smaller) ∆ determines the asymptotic value of the B-coefficient, at a given σ, only through the fact that a greater (smaller) part of the IX wave packet shall hit the edges of the slit.

IV. CONCLUSIONS
Scattering of composite particles is an important issue in several fields beyond semiconductors, such as molecular and nuclear scattering. 21,22 However, theoretical methods are relatively little developed, due to the numerical complexity. Here we have developed a numerical scheme to approach the particularly hard problem of a Coulomb bound complex scattering against arbitrary potentials, and we have applied the method to the specific case of IXs confined in CQW systems. On the one hand, this system offers the unique possibility to probe scattering of a composite particle with detailed optical means. On the other hand, we have shown here that full numer-ical propagation can be obtained for this system, due to the limited number of DoFs allowed by the low dimensionality ensuing from CQW quantum confinement. This allows both to compare with perspective experiments and to test approximate and numerically simpler methods, such as mean-field methods.
For the present case we have shown that mean-field methods are predictive only for external potentials which, while coupling the CM and relative dynamics, are sufficiently weak, i.e., with an energy scale which is much smaller than the first excitation gap for the internal motion of the IX. In this case, the CM and relative DoF can be factorized at any time of the evolution.
The RIX and the TDH approximations, however, are inadequate in predicting the transmission and reflection coefficients for stronger potentials, i.e., when the strength of the external potential is comparable to internal levels spacing. In this situation the scattering potential partially excites (locally) the IX to higher internal levels, and the wave function cannot be factorized during evolution in the potential region. In such a case, the asymptotic IX transmission, as computed from the full propagation, can be larger or smaller than the transmission extracted from mean-field methods, depending on the specific parameters of the system.
We have identified other signatures of the internal dynamics of the IX, coming into play during scattering, which cannot be captured by any mean-field calculation. For example, when plotted against time, the transmission through a well may exhibit plateaux before scattering is completed, which are absent in the simulations performed within RIX and TDH approximations. Moreover, in the single slit scattering problem, part of the wave-packet is not transmitted or reflected, but propagates along the edges. Again, this requires higher internal level excitations in the potential region through coupling of relative and CM DoFs which cannot be captured by a mean-field approach. This genuine correlation effect could be detected in optical experiments. It is also interesting to note that the RIX and TDH approximations basically coincide in almost all numerical simulations, indicating that, besides keeping the CM and relative subsystems decoupled, a mean field approach also smoothens the external potentials in such a way that it cannot exchange energy with the internal motion DoFs -which is de facto always in its ground state -but only with the CM one.
We finally note that, even for the case of only four DoFs, full wave-packet quantum propagation is a demanding task which required the development of a massively parallel code. On the one hand, exact calculations as the present one may serve as a severe benchmark for approximate, beyond mean-field methods which might be less computationally intensive. On the other hand, the implemented FSS method is suitable to treat also timedependent external potentials, which might be produced in heterostructures for driving single carriers or bound IX wave packets, such as surface acoustic waves produced by interdigital devices, at small additional numerical cost.  Supplemental Material to "Exact two-body quantum dynamics of an electron-hole pair in semiconductor coupled quantum wells: a time-dependent approach" where are single subsystem Hamiltonians, and T j (j) and U j (j) are the kinetic energy and potential energies, respectively, for the j − th subsystem. U (1, 2) is some potential which couples the coordinates for the two subsystems.
Consider Eq. (37). We choose the phases of A(1, t) and B(2, t) such that the scalar products of these functions with their time derivatives, computed at the same time, vanish. This is equivalent, as dt → 0, to A(1; t + dt)|A(1; t) = B(1; t + dt)|B(1; t) = 1, ∀t ∈ I. By inserting these equations into Eq.
To summarize, the quantum evolution of a system composed of subsystems 1 and 2 have been separated in the quantum evolutions each subsystem separately, each with a time-dependent Hamiltonians. For, say, subsystem 1, an effective potential u eff (1; t) arises, which can be written as the expectation value of the coupling potential averaged on the wave function of subsystem 2. Equivalently for particle 2, with particle coordinates exchanged. In this sense, this is a mean-field-like approach.
The Fourier split step method (FSS) for each of the two subsystems can then be applied to obtain the global evolution. In the present case, the systems 1 and 2 are representing the CM and relative motion coordinates, respectively; in this case, Eqs. (48) and (49) read as Eqs. (21) and (22) in the main text: