Controlled Coherent Coupling in a Quantum Dot Molecule Revealed by Ultrafast Four-Wave Mixing Spectroscopy

Semiconductor quantum dot molecules are considered promising candidates for quantum technological applications due to their wide tunability of optical properties and coverage of different energy scales associated with charge and spin physics. While previous works have studied the tunnel-coupling of the different excitonic charge complexes shared by the two quantum dots by conventional optical spectroscopy, we here report on the first demonstration of a coherently controlled interdot tunnel-coupling focusing on the quantum coherence of the optically active trion transitions. We employ ultrafast four-wave mixing spectroscopy to resonantly generate a quantum coherence in one trion complex, transfer it to and probe it in another trion configuration. With the help of theoretical modeling on different levels of complexity, we give an instructive explanation of the underlying coupling mechanism and dynamical processes.

This document contains: S1. Four-wave mixing experiment S1A. Experimental method and device S1B. Coherence dynamics S1C. Rabi rotations S1D. Identification of the neutral exciton-biexciton complex S1E. Power dependence of PL spectra S1F. Fitting 2D FWM spectra S2. Theory S2A. Semi-empirical model S2B. Calculation of 2D spectra for ultra-short pulses S2C. Dynamical simulations S2D. Inhomogeneous dephasing S2E. k · p model S1. Four-wave mixing experiment In this section we present details on the performed experiment and used sample (Sec. S1A).
We show standard characterization measurements of the coherence (Sec. S1B) and lightmatter coupling (Sec. S1C). Further, we identify the neutral exciton transition (Sec. S1D), demonstrate the power dependence of the photoluminescence (PL) spectra (Sec. S2E), and describe our fitting procedure for the 2D spectra (Sec. S1F).

S1A. Experimental method and device
We employ the heterodyne spectral interferometry technique introduced in Ref. 1 to detect four-wave mixing (FWM) signals from the quantum dot molecule (QDM). Its two-beam version was described in detail in the Supplementary Material of Ref. 2 , while its extension to a three-beams configuration was presented in the Supplementary Material of Ref. 3 . We use a Ti:Sapphire femto-second oscillator tuned into resonance with bright excitonic (trion) transitions of the QDM around 930 nm (1333 meV). In order to avoid driving spectrally nearby transitions, for instance biexcitons, the pulses are spectrally shaped (with a passive, diffraction grating-based pulse shaper) increasing their duration to ≈ 0.4 ps (full width at half maximum of the intensity spectrum). The main beam is then split into the reference pulse E R and three excitation pulses E 1,2,3 . In the latter, the consecutive pulses in the train are distinctly phase-shifted using acousto-optic modulators driven at the radio-frequencies  Figure S1: Schematic picture of the sample structure. Above a GaAs/AlAs distributed Bragg structure the InGaAs quantum dot molecule (QDM) layer is sandwiched between charge-doped layers, which are connected to a bias source. For details on the sample design and processing we refer to Ref. 5 . The optical in-and out-coupling is improved by a circular Bragg grating, also called a bulls-eye photonic structure, on the top surface. The laser pulses used for the optical excitation in the four-wave mixing experiment are focused on the center of the grating, while the reference beam hits the free surface in the vicinity.
(Ω 1 , Ω 2 , Ω 3 ) = (80, 79, 79.77) MHz. Using mechanical delay stages, they then acquire respective delays τ 12 and τ 23 , between the first two and last two arriving beams, respectively. E 1,2,3 are recombined into the same spatial mode and propagate unidirectionally. Using an external microscope objective, E 1,2,3 are focused on the sample surface, impinging the central part of the circular Bragg grating, while E R is placed on the flat surface around 15 µm apart ( Fig. S1). The signal is collected in the reflection direction. The reflection from the QDM, which also contains the FWM signal, is interfered with E R and heterodyned at the frequency Ω 3 + Ω 1 − Ω 1 in the same manner as explained in previous publications. [1][2][3][4] Note, that in this work we apply the degenerate two-pulse FWM scheme, i.e., τ 23 = 0.
The sample is kept in a He-flow cold-finger cryostat at the temperature of 7 K. To tune the spectral separation between quantum levels in the QDM, we vary the gate voltage ∆V    Figure S2 shows an exemplary coherence dynamics measurement by varying the delay τ 12 in the FWM experiment. While the neutral exciton line only decays due to dephasing processes, the trion transitions exhibit pronounced oscillations. These involved dynamics originate from the coherent coupling between the different states and lead to the off-diagonal peaks in the 2D FWM spectra (see e.g. Fig. S7).

S1B. Coherence dynamics
The beat dynamics in the FWM signal is suppressed when only one of the transitions is optically driven. For such a measurement the spectrally integrated FWM amplitude of the neutral exciton is shown in Fig. S3. By fitting the long-time decay with a single exponential function, we deduce a dephasing time of 200 ps. This dephasing rate is used for the simulations in Fig. S10. The initial rapid drop of the FWM amplitude is most likely stemming from a phonon-induced dephasing process 6 and is not further investigated here. delay τ 12 (ps) Figure S3: Spectrally integrated coherence dynamics as a function of τ 12 of the neutral exciton line, measured far from the hole tunneling resonance (blue points). The orange points represent the noise level. The long-time decay is fitted with a single exponential with the dephasing time 200 ps. The initial faster decay is due to phonon-induced dephasing.

S1C. Rabi rotations
To characterize the strength of the light-matter coupling of the QDM in Fig. S4 we plot the spectrally integrated FWM amplitude (for τ 12 = 0) as a function of the applied laser pulse amplitude. We find the expected Rabi rotation behavior, 7 which shows that a pulse area of π/2 is found for a pulse amplitude around √ P 1 = 0.7 √ µW. For all measurement in this work we restrict ourselves to small excitation powers in the linear regime of the Rabi rotation, i.e., √ P 1 < 0.05 √ µW. S1D. Identification of the neutral exciton-biexciton complex As stated in the main text, the bright spectral line in Fig. 3(a) stems from a neutral exciton transition. The standard proof of this type of transition is the detection of a coupled biexciton state. 2 We do this through the 2D FWM spectrum in Fig. S5, where the biexciton appears as off-diagonal peak with a binding energy of around 2.8 meV, which is a typical value for InGaAs QDs. 8

S1E. Power dependence of PL spectra
We have seen in the main text, that the bias ranges where the hole and trion resonances appear differ between the presented PL and FWM measurements. In Fig. S6(a) we show the same data as in Fig. 1(a) in the main text for an optical excitation power of P = 100 nW and in (b) the same measurement for P = 400 nW. We find the same spectral features of line-shifts and avoided-crossings but on different bias intervals. The range is not only shifted but also scaled by almost a factor of 2. This shows that the bias-dependence of the different state energies strongly depends on the specific optical excitation conditions. We also note that even more different scenarios are observed upon additional white light illumination (not shown). Varying optical excitation, changes the free carrier density and thus also the energetic level structure in the p-i-n diode and also the alignment of the quantized levels in the QDM itself. Therefore, because already the general optical excitation schemes in PL and FWM are off-resonant and resonant, respectively, and consequently entirely different, the two experiments cannot be quantitatively compared regarding their bias range properly.

S1F. Fitting 2D FWM spectra
To efficiently remove the prominent impact of the neutral exciton transition from the 2D spectra depicted in Fig. S7 we fit each spectrum with a sum of 2D Lorentzians. While for most studied bias values we use one diagonal peak for the exciton and two diagonal and two off-diagonal ones for the trion transitions, for the bias values ∆V = −307, −312, and −320 mV we have to consider one additional diagonal and two additional off-diagonal peaks.
This appearance of additional peaks in the bias range directly reflects the development of an avoided-crossing in the linear spectra (see Fig. 3(a) in the main text).
Our goal is to determine the peak ratio in the 2D spectra consistently in the measured and the simulated data. Therefore, we simply focus on four parts marked by colored areas (dashed squares) in Fig. S7(o). The signal amplitude summed over the red shaded areas gives the diagonal contributions and the sum over the green shaded areas the off-diagonal ones. This procedure can be easily automated to retrieve the peak ratios (green areas/red areas) for many bias values in the simulation, i.e., many 2D spectra. To determine the experimental peak ratios we use the fitted 2D spectra to reduce the uncertainties stemming from different discretization or spectral oscillations due to insufficiently long delay ranges.

S2. Theory
Here, we present the details of the semi-empirical theory for the coupled hole and trion system (Sec. S2A). This model yields transition energies and amplitudes for dipole-allowed transitions that are used for determining the FWM spectrum from the analytical formulas in the ultra-short pulse limit (Sec. S2B) or from dynamical simulations (Sec. S2C). Further, we discuss the influence of an inhomogeneous dephasing on the 2D spectra (Sec. S2D). Finally, we present the k · p model and the spectra calculated from the system's eigenstates and eigenenergies obtained in this way (Sec. S2E).

S2A. Semi-empirical model
In the following electron and hole destruction (creation) operators are given by a The relevant single particle states of the QDM are given by the single hole states forming the ground states, which are characterized by the QD number and the spin orientation: |1 ↑ , |1 ↓ , |2 ↑ , |2 ↓ . The relevant trions forming the excited states emerge by the creation of an additional electron-hole pair, with the electron always in QD1. These three-particle complexes now possess two holes, whose spin coupling results in new eigenstates that can be sorted into singlet S (m,n) (m, n ∈ {1, 2}) and triplet states |T α (α ∈ {+, −, 0}) defined by To get the full trion wave function we have to add the spin orientation σ ∈ {↑, ↓} of the remaining electron to the states such that we get σ, S (m,n) and |σ, T α . Overall we have 4 ground/hole states and 12 excited/trion states.
The ground/hole state Hamiltonian reads: where V 2 is the energy difference between the ground states in the two QDs, F = edE z is the dipole energy associated with the charge displacement in the axial electric field E z , where d is the effective distance between the QDs and e is the elementary charge, and t is the tunnel coupling.
The excited/trion state Hamiltonian reads (dropping the conserved electron spin index σ on the right-hand side): Here, J is the hole exchange splitting, V c is the Coulomb energy of exciton dissociation, and V 2 is the hole energy difference between the QDs in the presence of the exciton (taking the Coulomb correlation energy into account). The total trion Hamiltonian is The Hamiltonian for the electron-hole short-range exchange interaction is most conveniently written in the second quantization form, reading The considered values of the system parameters are V c = 20 meV, V 2 = 0 (this parameter merely sets the reference for the electric field; with this choice the F = 0 is fixed at the hole resonance), V 2 = 0.785 meV, δ = 0.21 meV, J = 0. With these parameters, the model is equivalent to the one used in Ref. 5 The dipole operators for σ + and σ − polarized light are, respectively, and which determines the optically active transitions in the system.

S2B. Calculation of 2D spectra for ultra-short pulses
Using the matrix of transition amplitudes between the hole and trion states M iu , we can derive the two-pulse FWM spectrum in a general analytical way in the limit of ultra-short laser pulses.
The diagonalized Hamiltonian is of the form where |h i and |T j are the hole and electron eigenstates, respectively, with the corresponding eigenenergies E i and E u , and M is one of the operators in Eq. (S6) or Eq. (S7), depending on the polarization of light.
In the interaction picture this reads Before any optical excitation the system state is a mixture of hole states represented by a density matrix ρ (0) . The first optical excitation happens at t = −τ 12 and we assume that the laser pulse is much faster than any internal dynamics, i.e., dephasing, tunneling, or decay.
The density matrix immediately after this pulse is then given by where h.c. contains the phase +ϕ 1 , which is irrelevant for the final FWM signal and can therefore directly be neglected. After introducing the pulse spectrum, the FWM-relevant density matrix after the first pulse reads Before the interaction with the second pulse, the system can undergo dephasing and decay dynamics, such that the density matrix immediately before the next pulse is given by where γ iu is the dephasing rate for a given coherence.
The second pulse acts via a second order process, which leads to the density matrix When considering again only contributions with the correct phase dependence for the FWM signal (∼ e 2iϕ 2 ) we arrive at Finally, taking into account that the emission resulting from a coherence |T v h k | is proportional to M * kv e (−iω kv −γ kv )t , we arrive at the FWM signal Assuming that the spectrum of the laser pulses is much broader than the energy differences between the relevant transitions, we can setŝ(ω − ω mn ) ≈ 1. We further assume that the initial state is a statistical mixture of the hole states, i.e., where we take the thermal distribution for the probabilities p n . With this the FWM signal can be simplified to From this we directly find the 2D FWM spectrum where f iu (ω) is broadened by γ iu but also takes the spectral broadening by a finite spectral resolution in the detection process into account.  Figure S9: Simulated 2D FWM spectrum using Eq. (S17) for different external electric fields as labeled in each panel.
A few examples for the analytically calculated 2D spectra are shown in Fig. S9 for different bias values as given in the plot. We assume that the spectral broadening is dominated by the instrumental resolution and use the same value as in Fig. 2 in the main text, which roughly agrees with the spectrometer resolution in the experiment. We clearly see that the off-diagonal peaks only appear when the states are efficiently tunnel-coupled in the vicinity of the avoided-crossings in the bias scan (Fig. 2 or S10).

S2C. Dynamical simulations
To simulate the FWM signal while taking a non-vanishing pulse duration into account we calculate the system's dynamics by solving the Lindblad equation for several combinations of pulse phases e iϕ 1 and e iϕ 2 . 9, 10 We then phase-filter the optically active coherences p in the density matrix by numerically integrating ϕ 1 ,ϕ 2 pe −i(2ϕ 2 −ϕ 1 ) dϕ 1 dϕ 2 . For a given pulse delay τ 12 we then determine the FWM spectrum by Fourier transforming the filtered coherence with respect to the real time t after the last arriving laser pulse. The bias scan of the FWM spectrum calculated dynamically is depicted in Fig. S10. We consider a pulse duration

S2D. Inhomogeneous dephasing
In general, inhomogeneous dephasing of a single emitter system originates from fluctuations of the transition energy which happens on timescales slower than a single repetition of the FWM experiment but faster than the repetition rate of the laser. Formally, this situation corresponds to an ensemble measurement. In order to include such slow fluctuations, we assume the transition frequencies ω kv , ω nu in Eq. (S16) to be Gaussian random variables with mean values ω kv , ω nu and standard deviation σ, which are uncorrelated if k = n or v = u. Upon averaging over the fluctuations, we arrive at where the dephasing factor is f deph (t, τ 12 ) = e −σ 2 (t 2 +τ 2 12 ) , k = n or v = u, e −σ 2 (t−τ 12 ) 2 , k = n and v = u.
The impact of the inhomogeneous dephasing is shown in Fig. S11, where we show the 2D in (a) and σ = 0.1 meV in (b). When presenting the results, we subtract the "background" formed by the power-law tails of the strong diagonal peaks in the spectral area of the offdiagonal ones. We clearly see that the previously sharp peaks get additionally broadened predominantly into the diagonal direction in the 2D spectrum, as well known from several previous works. 11 In Fig. S11(c) we show peak ratios of the form A off-diag / A diag as in the main text and increase the inhomogeneous broadening as given in the plot. The values of A off-diag and A diag are determined by integrating the simulated signal over the areas shown in Fig. S7(o), with A off-diag representing the difference of the counting results for the full signal and for the signal with off-diagonal contributions artificially switched off (when only the tails of the strong diagonal peaks contribute). As discussed in the main text, the additional dephasing significantly reduces the visibility of the off-diagonal peaks representing the coherent coupling between the different trion states.

S2E. k·p model
Our k · p model describes a vertically aligned QDM formed by InGaAs in a GaAs matrix.
Both QDs are modeled as lens-shaped with the geometry of the bottom QD being set to 18.0 nm in diameter and 4.0 nm in height, while for the top QD it is 22.0 nm and 4.4 nm, respectively. The dots are separated by a distance of 8 nm (base to base) and they are placed on 0.6 nm thick wetting layers (WLs). Note, that all the values are subject to discretization with a mesh adjusted to 0.565325 nm, which is the lattice constant of GaAs at low temperatures. The maximum In content in the dots and in the WL is 45% with a gradual decrease of indium at the edges, which simulates material intermixing. This is done by processing the composition profile with a Gaussian blur with a standard deviation of 0.6 nm. We calculate single-particle electron and hole states using an eight-band k · p model with the envelope function approximation. 12,13 The strain which arises due to the lattice mismatch of InAs and GaAs is determined by the continuous elasticity approach. 14 The piezoelectric field is included up to the second order in polarization, 15,16 using parameters from Ref. 17 Details of the model are given in Ref. 18 with further improvements of Ref. 19 The trion states are calculated using the configuration-interaction (CI) method, where we consider a basis of 4 electron and 8 hole single-particle states. The electric field is accounted for at the stage of the exact diagonalization of the CI Hamiltonian 20 (with the matrix elements calculated in the basis of single-particle states obtained at zero electric field). The electron-hole exchange is also included at the CI stage. The relevant matrix elements are calculated by mesoscopic expansion of the Coulomb interaction up to second order as in Ref. 21 We then estimate the atomic integrals using rescaled hydrogen-like wave functions 22 with the effective value of the dielectric constant for local (on-site) integrals being equal to 4. The single-and three-particle eigenstates are then used to compute the amplitudes of the optical transitions. 23 The PL and FWM spectra obtained from the k · p modelling are shown in Fig. S12.
Qualitatively, in spite of the much reacher basis, the k · p approach yields the same features in the PL, as the phenomenological model, including the hole and trion resonances as well as the characteristic crossing of transition lines. The FWM spectrum also has exactly the same  Figure S12: The luminescence (a) and FWM (b) spectra simulated using the system states obtained from the k · p model. structure, showing essentially two horizontal branches interrupted by avoided-crossings and with intensities decaying in the opposite directions. Quantitatively, the energy scale (which is determined by the Coulomb correlation energies of the three-particle system) is larger in our k · p simulation because of a smaller inter-dot distance, which turned out to be necessary to guarantee the stability of the relevant states in the absence of the AlGaAs barrier in our simplified geometry of the model. One additional feature in the k · p -based PL spectrum is the line at about 1343 meV that crosses all other transition lines without any interaction (avoided-crossing). This line belongs to a transition between states with an excited hole, i.e., electron-hole recombination in the presence of a spectator hole in an excited state (in the p-shell). It is visible in our PL simulation, where all the trion states are assumed to be equally occupied. This is consistent with the measurements shown in Fig. S6, where a similar line becomes pronounced at higher excitation powers. There is no corresponding line in the FWM spectrum, since this is a resonantly excited coherent measurement in which only the lowest hole states are involved.