Radial and Angular Electron Ejection Patterns of Two-Electron Quantum Dots in THz Fields

: In this work, we develop and apply an ab-initio method to calculate the joint radial- and-angular electron distributions following the interaction of two-electron spherical quantum dots (QD) with intense terahertz pulses of subpicosecond duration. By applying the method to two QDs of different size, we could investigate two particular ionization mechanisms: the direct and the sequential two-photon double ionization. According to our results, the two ionization mechanisms show some similarity in the angular distribution patterns, whereas the corresponding radial distributions are distinctly different, associated with their joint kinetic energy spectrum. We also discuss the time-evolution of the ionization process in the context of the different nature of the interaction of the QD with the external radiation and the electron–electron correlation interactions.


Introduction
The study of the optical properties of semiconductor quantum dots (QDs), wires, and wells has attracted a lot of attention in the domains of fundamental theory and applications for quantum information processing [1], solar energy harvesters, optoelectronics, and digital imaging [2]. Quantum confined structures (QCS) exhibit a variety of enhanced optical properties compared with atoms, molecules, and bulk materials. This is due to the high degree of flexibility in QCS design where it is possible to artificially tailor their transport and optical properties at will. Nonlinear spectroscopy has proven to be a powerful and widely used tool that can not only to probe the electronic structure and the electron-hole dynamics in quantum systems but also has important technological applications [3,4]. For example, the nonlinearities in QCS have many applications in laser technology to create materials with very large values of second-and third-order susceptibilities, which are essential for harmonic generation radiation in the mid-infrared regime [5][6][7][8][9][10].
In the past, most of the attention was paid to the non-linear properties of coupled dots (molecular quantum dots) in strong laser fields (see [11], and references therein). In contrast, the non-linear properties of two-electron QDs in laser fields is a largely unexplored scientific area. The main reason is that in QDs the electrons are confined within the QD potential and therefore they strongly interact, thus complicating the representation of the electronic structure. The existing theoretical approaches rely on the description of a considerably simplified system by ignoring part of the interelectronic correlations (essential states approximation), an approach which may be not valid when the applied laser field is strong (i.e., when the internal forces acting on each electron are of the order of the laser-electron interaction). In the latter case, the electronic QD continuum and multiply-excited states play a decisive role in the system's dynamics. In this context, it is well known that the theoretical investigations (and associated numerical and computational implementation) of nonlinear processes (e.g., the above threshold ionization and high-order harmonic generation (HOHG)) is a demanding problem both from the viewpoint of theory and implementation (computationally).
In the present work, we apply an ab-initio electronic structure method to fully describe the double ionization of two-electron semiconductor QDs, following the interaction with a terahertz (THz) laser field. The ab-initio calculation method using a numerical local basis has been thoroughly tested and applied to describe non-linear interactions of ultrashort fields with atomic and the molecular systems [12][13][14][15][16][17][18][19] throughout the years, while its adaptation and application to the case of QDs is described in detail in [20,21]. Here, we extend our theory to include the calculation of the joint radial and angular distributions of the ejected electrons, a subject that is theoretically and computationally highly non-trivial, mainly because the configuration space of the electrons extends to thousands of times the size of the residual target (here, the QD).
The method is based on a non-perturbative solution of the time dependent Schrödinger equation (TDSE) for the QD interacting with a linearly-polarized THz pulse. The model at hand describes a spherical two-electron quantum dot built from a narrow bad gap semiconductor crystal of approximate radius R Q . Within the effective mass approximation, the crystal potential is taken into account via the electrons effective mass, m * e , and the crystal's dielectric constant, κ. In this case, the TDSE that describes the QD is given by: i ∂ ∂t ψ(r 1 , r 2 , t) = Ĥ m * e ,k (r 1 , r 2 ) +D m * e ,k (r 1 , r 2 , t) ψ(r 1 , r 2 , t), where H is the field-free QD two-electron Hamiltonian, ψ is the time-dependent wavefunction, and D describes the interaction of the QD with the radiation field. In the above, the electron's coordinates are denoted as (r 1 , r 2 ) . The field-free QD Hamiltonian is modeled through a Gaussian potential as, whereĥ m * e (r i ) is the one-electron QD's confinement Hamiltonian. The above expressions constitute the fundamental equations of the present formulation. The details of calculating the electronic structure, the time-dependent ψ(r 1 , r 2 , t), and the joint radial and angular distributions are discussed in the next section.

Calculational Method of the Two-Electron QD Dynamics
In the following, the formulation for the ab-initio calculation of the radial and angular distributions of the two ejected electrons is developed. Briefly, we calculate the two-electron wavefunction at the end of the THz pulse and then extract the desired information either by projecting onto the asymptotic two-electron wavefunction which describes two uncorrelated electrons or by evaluating the squared modulus of the two-electron wavefunction.
As the method of solving for the electronic structure of the quantum dot and its dynamics under the interaction with a laser field is discussed in detail in [20,21], we make a short discourse of this and the equations necessary for the current theory are presented.

Calculation of the QD Electronic Structure
We assume the electrons confined in a spherically symmetric Gaussian potential well of depth V 0 and width parameter β, where β is related to the quantum dot radius by R Q = ln 2/β. In this case, the one-electron QD confinement Hamiltonian is modeled as, The model also takes into account the inter-electronic interaction via the electrostatic Coulomb potential. The equations used to represent the quantum dot-laser system are presented in scaled atomic units (s.a.u), which is discussed in [20]. We thus replace the set of the material constants by Q = (m * e , k, β, V 0 ) and denote the one-electron Hamiltonian byĥ q and the two-electron QD Hamiltonian in Equation (2) as, H Q . Therefore, the eigenvalue problem to be solved is the time-independent Schrödinger equation (TISE), where Φ NΛ (r 1 , r 2 ) are the eigenstates ofĤ Q (two-electron field free wavefunctions) and the index Λ is shorthand for (L, S, M L , M S ), representing the total orbital and spin angular momenta, and their projections onto the z-axis, respectively. While the notation is kept general, only the singlet symmetry case is considered where (S, M S ) = (0, 0). The system is assumed to be confined in a sphere of radius R, much larger than the QD's size, R >> R Q . Within our particular approach, we implement the so-called fixed boundary conditions, which require the wavefunction to strictly vanish at the origin and the boundaries, r = R. Consequently, all QD's states (including the continuum) become discretized, allowing the bound and continuum states to be treated in the same way, all subject to the unity normalization. For the solution of equation (4), a configuration interaction (CI) method is employed where the quantum dot eigenstates are expanded on an uncorrelated two-electron basis φ Λ n 1 l 1 n 2 l 2 (r 1 , r 2 ), formed by Slater determinants of angular momentum coupled one-electron states (configurations) found from the one-electron TISE. More specifically, the uncorrelated two-electron basis functions are the eigensolutions of the Hamiltonian Equation (2) but without the 1/|r 1 − r 2 | interaction term. The singlet antisymmetric uncorrelated basis with sharp angular momentum values, (L, M L ), respectively, are as below, where Y LM L l 1 l 2 (Ω 1 , Ω 2 ) are the bipolar spherical harmonics, containing the angular momentum coupling coefficients (Clebsch-Gordon coefficients) and A 12 is the anti-symmetrization operator which acts to exchange the coordinates of the two electrons. The purely radial functions, P nl (r), are calculated as solutions of the radial one-electron eigenvalue QD Hamiltonian, The radial functions are numerically solved for by expanding on a basis of B-Splines with excellent properties for representing continuum states [22? ]. Having completed these calculations, the correlated two-electron basis in terms of the uncorrelated basis may be written as, where the expansion coefficients ν NΛ n 1 l 1 n 2 l 2 are called the configuration CI coefficients. Substituting (7) into the TISE and projecting over φ Λ n 1 l 1 n 2 l 2 (r 1 , r 2 ) converts it into a matrix equation, which upon diagonalization retrieves the eigenenergies E NΛ and the CI coefficients ν NΛ n 1 l 1 n 2 l 2 . The energies of the ground state, single ionization threshold, and double ionization threshold are denoted by E 0 , E 1 , and E 2 , respectively, where E 2 ≡ 0 [sketches in Figure 1]. If a two-electron eigenstate has an energy E NΛ < E 1 , then the state is of bound character. States with energies E 1 < E NΛ < E 2 represent singly ionized QD with an ejected electron (including resonance states), and states with E NΛ > E 2 may represent both singly and doubly ionized QDs with one-electron ejected or two-electrons ejected, respectively.

Quantum Dot Laser Field Dynamics
Having solved for the electronic structure of the two-electron quantum dot, its dynamics under the influence of an intense and ultra-short linearly polarized laser pulse will be solved for. This amounts to solving the time dependent Schrödinger equation (TDSE) for the quantum dot-field system (Equation (1)). We use a semiclassical representation for the laser field in the Coulomb gauge and the long wavelength approximation for the field, whereẑ is a unit vector (not an operator) along the z-axis, representing the linear polarization of the laser pulse along this axis. p 1 and p 2 are the electron momenta and A(t) is the electromagnetic potential field of the laser pulse, E(t) = −Ȧ(t) in the Coulomb gauge (E(t) is the electric amplitude), given by, where E 0 is the field amplitude and ω is the carrier frequency. For a realistic ultrashort laser pulse, we use a squared sinusoidal envelope, which satisfies the requirements that the envelope varies slowly with respect to the carrier period, and rises and falls to zero. τ p is the laser pulse duration, related to the field period (T 0 = 2π/ω) by τ p = n c T 0 , where n c is the number of field cycles in the pulse. The solution for Equation (1) is found by an expansion of the time-dependent wavefunction on the two-electron eigenstates with time dependent coefficients as, The quantity |C NΛ (t)| 2 represents the probability of observing the system in state Φ NΛ (r 1 , r 2 ) at time t or alternatively the state population. Solving the TDSE by substitution into Equation (1) and multiplying from the left by Φ * N Λ (r 1 , r 2 ) and integrating over the entire coordinate space transforms the TDSE into a set of coupled ordinary differential equations, where D NΛN Λ are the dipole matrix elements between the states Φ NΛ and Φ N Λ . Standard angular momentum algebra allows for the expression of the two-electron dipole matrix elements in terms of the ν NΛ n 1 l 1 n 2 l 2 CI coefficients and the one-electron dipole matrix elements; the latter are calculated numerically given the P nl (r) radial functions [22]. Thus, solving the TDSE amounts to calculating only the two-electron dipole matrix elements and integrating Equation (11) to find the time dependent coefficients, C NΛ (t).
For the reader who is interested, it is advised to look at the works in [20,21] for a more detailed exposition of the method for QDs. We are now ready to proceed with the main subject of the present work which are the radial and angular distributions of the doubly ionized QD system.
The time-dependent coefficients, C NΛ (t), along with the CI coefficients, ν NΛ n 1 l 1 n 2 l 2 , allow the calculation of various observables of the system, via the expansions (10) and (7); ultimately, the radial and angular distributions are expressed in terms of these coefficients and radial wavefunctions P nl (r), albeit in a complicated fashion.

Two-Electron Joint Radial Probability Distributions
We start by giving the proper definitions for the probabilities to be calculated. The radial probability density is defined to be, With the density automatically normalized to unity, dr 1 dr 2 P r (r 1 , r 2 ) = 1, since, for the total time-dependent wavefunction, d 3 r 1 d 3 r 2 |Ψ(r 1 , r 2 , t)| 2 = 1.
Calculating the two electron radial distribution amounts to evaluating the square modulus of the time-dependent wavefunction, expanded on the spectral basis (which is itself expanded as a sum of determinants) and integrating over all angular coordinates. Starting with the time-dependent wavefunction expansion, Equation (10) and the two-electron eigenstate expansions, Equations (7) and (5), we have, This is eventually written as, with the partial-wave channels, f Λ l 1 l 2 (r 1 , r 2 , t), defined as, and the two-electron radial functions by, ρ n 1 l 1 ;n 2 l 2 (r 1 , r 2 ) ≡      P n 1 l 1 (r 1 )P n 2 l 2 (r 2 ), n 1 l 1 = n 2 l 2 1 2 P n 1 l 1 (r 1 )P n 2 l 2 (r 2 ) + (−) l 1 +l 2 +L P n 1 l 1 (r 2 )P n 2 l 2 (r 1 ) . n 1 l 1 = n 2 l 2 (16) In the above, the exchange term results from the anti-symmetrization of the two-electron eigenstates. The appearance of the phase factor (−) l 1 +l 2 +L is attributed to the symmetry properties of the Clebsch-Gordon coefficients under particle exchange. For the case where we want to examine only the singlet symmetry, combinations of l 1 , l 2 and L are taken such that the phase factor is always +1. Now, using the definition of the radial probability density and the orthonormality of the bipolar spherical harmonics, we have We then arrive at our main equation for the radial probability distribution as, This form of the radial distribution equation is particularly useful as it allows the decomposition of the full radial distribution into partial waves, (l 1 , l 2 ). The above equation in combination with Equation (15) provides the total radial distribution including both the bound and continuum states of the QD. For our purposes, in the summation of Equation (15), we must take care to include only index combinations corresponding to double ionization (DI) states. The energy of the states in the double continuum is, E N = n 1 + n 2 ; thus, the minimum values of N, n 1 , n 2 are set by the conditions E N > 0, and n 1 > 0, n 2 > 0, which effectively exclude the bound and the single ionization states from the summation.

Two-Electron Joint Angular Probability Distribution
Attempting to calculate the square of the wavefunction expanded on the spectral basis and integrating over radial coordinates leaves an equation that is not practical to evaluate. Instead, to calculate the angular distributions, an asymptotic momentum basis is employed. It is valid to assume asymptotic conditions after the laser pulse has finished and sufficient time has elapsed, as, if t → ∞, then r 1 , r 2 , |r 1 − r 2 | → ∞, i.e the two electrons asymptotically move as free electrons.
A method is described below to express the DI probability distributions following double electron ejection in a photoionization process. For this, we need to define the asymptotic momentum two-electron wavepackets, in the context of a doubly-ionized QD. Asymptotically (at detector's location), the QD's asymptotic Hamiltonian and the momentum operator commute with each other and therefore their eigenstates may serve as a common basis to express the quantum mechanical state of the system of the two-electrons; these are denoted here as φ (3), of the ith electron. Note that the superscript (−) refers to ingoing asymptotic behavior, but for the sake of conciseness φ (−) k 1 k 2 (r 1 , r 2 ) is denoted by φ k 1 k 2 (r 1 , r 2 ) in the following. From the above, it is easily concluded that the two-electron states φ k 1 k 2 (r 1 , r 2 ) are expressed as the product of states of φ k (r): The positive-energy eigenstates of this far-distance Hamiltonian φ k describe states of energy = k 2 /2 with momentum vector, k: The next step is to express φ k (r) in the partial wave basis of φ nlm (r), defined by the common set of commuting observables h, l 2 , l z and denoted by φ nlm = φ n lm , Now, we need to express φ k 1 k 2 (r 1 , r 2 ) on the two-electron 'partial wave basis' Φ (Λ) k 1 l 1 ;k 2 l 2 (r 1 , r 2 ), which in turn is properly interpreted as the CI basis Φ (Λ) n 1 l 1 ;n 2 l 2 (r 1 , r 2 ) of Equation (5). The short-range phase shifts δ kl are here defined from the calculated radial eigenstates of Equation (6) relative the corresponding asymptotic Bessel functions (with the same indices). For the imposed zero boundary conditions, P nl (R) = 0, we have, P nl (r) → 2 k n π sin(k n r − l π 2 + δ k n l ), n > 0.
The wavenumber, k, is forced to take discrete values, k → k n . This is due to the fact that, assuming the QD is enclosed in a spherical box, all states (both negative and positive energy) are actually discrete. Thus, with the convention that, for positive energies, 1 > 0, 2 > 0: In fact, the values of the discretized wavenumber values depend on the l angular quantum number, k n = k n (l). In the following, we suppress this dependence. We make the final results insensitive to this dependence by enlarging the box so that to have a sufficient positive-energy density of states. Thus, the indices n and k are equally valid to denote the same state with n = k 2 n /2, namely φ nlm (r) ≡ φ klm (r). Direct replacement of (21) into φ k 1 k 2 (r 1 , r 2 ) gives, Now, we can replace 1 → k 1 and 2 → k 2 and then express the product of the single-electron states φ k 1 l 1 m 1 (r 1 )φ k 2 l 2 m 2 (r 2 ) on the basis of the two-electron states with E = 1 + 2 , total angular momentum L = l 1 + l 2 , andL z = l z 1 + l z 2 , namely the Φ LM L (r 1 , r 2 ). Following a standard angular momentum algebra, we arrive at the expansion below for the asymptotic solutions: where C l 1 m 1 ;l 2 m 2 LM L is a Clebsch-Gordon angular momentum coupling coefficient. Φ Λ k 1 l 1 ;k 2 l 2 (r 1 , r 2 ) are angular momentum coupled two-electron basis functions.
We are now ready to proceed with calculation of the joint two-electron angular probability. The differential probability of measuring the system in the state φ k 1 k 2 is given by the square of the projection of the time-dependent wavefunction evaluated after the end of the pulse, written formally as The time-dependent wavefunction is written as, |Ψ(r 1 , r 2 , t) = ∑ NΛ C NΛ (t) ∑ l 1 n 1 ;l 2 n 2 ν NΛ n 1 l 1 n 2 l 2 A 12 |φ Λ l 1 n 1 l 2 n 2 , where A 12 is the anti-symmetrization operator. The projection is calculated as, Again, following standard algebraic manipulations, we arrive at the below compact expression, with the CI matrix elements V NL n 1 l 1 n 2 l 2 defined as, , n 1 l 1 = n 2 l 2 , 2ν NL n 1 l 1 n 2 l 2 , n 1 l 1 = n 2 l 2 . (30) Now, by replacing dk i = k 2 i dk iki , i = 1, 2, we end up at, The final step is integration over all possible kinetic energies, and, since these are discretized, integration is equivalent to a summation over the n 1 , n 2 indices, which gives: Equations (29) and (32) constitute the two main equations used to calculate the joint two-electron angular distributions. We extract the distributions corresponding to the DI by imposing the conditions described in the case of the DI radial distribution for Equation (15).

Results and Discussion
The quantum dot structure is solved for using expansion (7) with angular momentum configurations for total angular angular momenta L = 0, 1, 2, and 3, chosen such that we are examining only the singlet states. We have kept all configurations with l up to 3, which are adequate for our purposes. The box radius chosen for the calculation was R b ≈ 160 nm.
As the QD structure is size dependent, the radial and angular distributions are presented for two different sizes. The first quantum dot (Q 1 ) is chosen with radius of 4.6 nm, while the second (Q 2 ) has radius of 3.2 nm. Both quantum dots are built from the same semiconductor crystal with an electron effective mass and dielectric constant of m * e = 0.1 × m e and κ = 5 × 1/4π 0 , respectively, with m e as the vacuum electron mass and 0 as the vacuum permittivity. The TDSE is then propagated for a laser pulse of central carrier frequency ω = 304.7 meV, at an intensity I 0 = 6.397 × 10 5 W/cm 2 . The pulse duration is ≈ 0.163 ps (12 cycles). Given these laser pulse parameters, there are two distinct double ionization mechanisms that can occur. For Q 1 , assuming the quantum dot is initially in its ground state (E 0 =−582.8 meV), then E 0 + 2ω > E 2 and thus double ionization via the simultaneous absorption of two photons is an open ionization channel, known as direct double ionization [see left sketch of Figure 1]. For E 1 = −334.4 meV, it holds that E 1 + ω < E 2 and the probability for further ionization of Q + 1 , via single-photon absorption, is significantly diminished as it is not energetically favorable. Double ionization can then occur, only if the Q + 1 ionizes by absorbing two photons or more. We then expect for the chosen range of intensities the dominant double ionization channel to be the two-photon direct ionization from the ground state (since the sequential channel necessitates a three-photon absorption which becomes important only towards higher intensities).
In contrast, the situation for the Q 2 is very different. For Q 2 with ground state energy E 0 = −365 meV, the direct two-photon double ionization is again energetically possible, since E 0 + 2ω > E 2 . However, now, double ionization can also proceed sequentially from Q 2 to Q + 2 by one photon absorption and then one further ionization of Q + 2 by absorbing one more photon (since for E 1 = −232.9 meV we have E 1 + ω > E 2 ); the important difference is that the sequential ionization via Q + 2 requires two-photon as well [See right sketch of Figure 1]. In the left sketch (R Q = 4.6 nm) is the direct double ionization (DDI) mechanism, where one-photon ionization of the singly-ionized QD is not energetically allowed. The right sketch (R Q = 3.2 nm) is the sequential double ionization mechanism, where one-photon ionization from the ground state of the neutral and the singly-ionized QD is energetically allowed.
Thus, for Q 1 , the dominant ionization regime is the direct double ionization (DDI) channel and for Q 2 , sequential double ionization (SDI). The main difference is that for the former, while the primary interaction is the radiation field, the electron-electron interaction may contribute significantly in the overall double ionization process, whereas for the SDI case the electron-electron interaction is greatly diminished in significance; theoretically, it is not a prerequisite for double ionization to occur.
As a result of the importance of the electron-electron correlations between the SDI and DDI, we want to examine what quantitative differences show up in the radial and angular probability patterns for the ejected electrons following double ionization between Q 1 and Q 2 .

Direct Double Ionization Regime (R Q = 4.6 nm)
The form of the radial distribution Equation (17) is such that the total radial distributions can be decomposed into partial wave contributions (l 1 , l 2 ). While there are 10 different partial wave contributions present in the calculations, only the four most dominant partial waves are presented here, as shown in Figure 2. As discussed above, in the DDI case, the absorption of two photons by the QD directly from its ground state will lead to the simultaneous ejection of the two interacting electrons, and thus inter-electronic interactions could play an important role in the double ionization process.
It is not a trivial task to demonstrate undoubtedly the individual role the laser field and the electron-electron (e-e) interaction in the overall double ionization process but we try below to interpret in a qualitative fashion the radial distribution patterns. We notice that the most dominant configuration channels are those of the pp, and then that of the sd configurations, followed by the ss and dd channels. This order of significance could be attributed to the increased need of electron-electron correlation interactions going from the pp to the other channels.
In the case of the pp channel the angular momentum and energy transfers are somewhat consistent and may be completed without the need of e-e interactions; this is because each photon carries an angular momentum of one and as such one could think of the pp process as each absorbed photon transferring its angular momentum to the individual electrons in QD's ground state and the electrons liberate (with s → p for each ejected electron); the residual e-e interactions do cause changes only on the (primary acquired) electron's energies on their way out the QD region and thus the peak structure in the spectrum, leading to a highly peaked radial pattern. Thinking along similar lines for the other peaks, we may infer that e-e interactions must occur in order to be observed. For example, in the case of the sd peak, it appears that one of the electrons has 'absorbed' the two-photon's angular momenta (s → d); this is in contrast with the energy sharing pattern where the two-electrons appear to have acquired similar energies. Alternatively, the sd configuration in the DI channel may result from p → s and p → d one-photon absorption from the pp component of the ground state, which is discussed in Section 6. The inter-electronic interactions allow for the energy of the two photons to be shared between the two electrons, as shown in [21], where it is found that in the direct regime the kinetic energy distributions tend to have their peaks at equal energies, 1 ≈ 2 . In accordance to this observation, it is not surprising that the structure of the partial radial distributions in Figure 2 show peaks centered around equal distances, r 1 ≈ r 2 . For example, for the dominant peak of the, (p, p)-wave, we see that two-electron wavepacket is about 25 nm away from the QD's center travelling outwards.

Sequential Double Ionization Regime (R Q = 3.2 nm)
As mentioned above, the smaller quantum dot has lower ionization thresholds, thus providing access to the sequential double ionization mechanism. While both sequential and direct double ionization contribute, it is immediately clear in Figure 3 due to the drastic changes in structure of the radial distributions that the sequential mechanism dominates the double ionization process. In addition, comparing, for example, the pp contribution in both DDI and SDI, the scale factor for DDI is 10 −12 , whereas, for SDI, it is 10 −9 . An in depth discussion on the relative contributions of the partial waves and 1 L symmetries to the total kinetic energy distribution is provided in [20,21], and the same considerations can be carried here, specifically the dominance of the pp channel by three orders of magnitude over the others. It is clear from the separation of the two peaks that the pp channel is populated mainly by SDI, while for the other channels these peaks are affected by the electronic correlation, suggesting competition between SDI and DDI. In terms of e-e correlation effects the DDI processes are certainly of less interest, as they are overwhelmed by the primary interaction with the THz laser radiation. Finally, we can draw our attention to the total distributions presented in Figure 4, where, in the SDI regime (right), the dominance of the pp wave can be clearly seen by the fact that, even though the total distributions is represented simply as a sum of the partial distributions, the total distribution looks almost exactly the same as the pp distribution which has mostly washed out the effect of the other distributions (the main difference being a higher probability along equal positions, attributed to e-e correlation). On the other hand, in the DDI regime, the total distribution (left of Figure 4) results from four dominant partial distributions with comparable relative probabilities, giving a distribution that cannot be associated to an overwhelming partial wave, bringing to light again the role of e-e correlation.
We may, however, try and associate the peak positions in the distributions relative to the expected peak positions resulting from a simplified energy calculation. With the excess energy (kinetic energy) of the electrons after photoionization has occurred, E = k 2 1 /2 + k 2 2 /2, and assuming that the peak strength of photoionization process occurs at the electric field maximum of the laser pulse at a time ∼ τ p /2, we can find the classical velocity of the electron from the kinetic energy and calculate how far the electron is expected to travel after the second half of the laser pulse has finished, r i ∼ v i τ p , i = 1, 2.
For the SDI case (right panel of Figure 4) the kinetic energy peaks of the two electrons for the pp-wave are 174.14 meV (1.6 s.a.u) and 65.3 meV (0.6 s.a.u.) [21]. In addition, for reference, 1 s.a.u of time is 0.006047 ps. For the laser pulse given here with a pulse duration of 0.163 ps, the expected peak position in the joint radial distribution is placed at (r 1 , r 2 ) ≈ (61 nm, 39 nm). Comparing this with the pp radial distribution, we can see that the peaks are approximately present at this position. On the other hand, in the case of DDI (left panel of Figure 4), with both electrons having a kinetic energy of 16.33 meV (0.15 s.a.u) for the pp-wave, the resultant expected position for both electrons is ≈20 nm and comparing with the pp radial distribution we can see again that the peak falls approximately on this position. As the expected positions are the result of a rough estimation based on classical kinematics relations, a discrepancy between the radial distributions and the expected positions is to be expected.

Angular Distributions
Equation (32) is evaluated using the asymptotic basis states forming ψ k 1 k 2 (r 1 , r 2 ) for momenta calculated from the relation i = k 2 i /2 where i = 1,2. Since the joint two-electron angular probability distribution has four angle parameters (ϑ 1 , ϕ 1 , ϑ 2 , ϕ 2 ) there are various ways of plotting. The choice here is to fix the ejection angles of one of the electrons (ϑ 1 , ϕ 1 ) = 0 and then to plot the angular probability for the second electron for a fixed ϕ 2 = 0 [see Figure 5]. As can be seen in Figure 6, the direction of ejection of the second electron relative to the first is primarily along the parallel (0 • ) and anti-parallel (180 • or back-to-back) directions. It is clear that the total distribution is dominated by the contribution of the pp channel, which seems to imply that the two electrons escape primarily due to absorption of a single photon each. Attention should now be drawn to the fact that the peak ratio of anti-parallel to parallel ejection in the direct regime is ≈2.5 times greater. Intuitively, this would be expected in a classical picture of the simultaneous ejection of two electrons, where one can expect the electrons to move in opposite directions due to the interelectronic interaction. Corroborated by the discussion of the radial distributions in the DDI case, it should be concluded that the ejection pattern of the two-electrons is mainly determined by channels due to the interaction with the laser field with a lower contribution originating from the e-e interaction channels. These e-e channels contribute destructively for electrons ejected along the same direction θ 1 ∼ θ 2 ∼ 0 as they lead to a reduction of the pp pattern in the overall angular distribution probability along these angles, in contrast to the back-to-back ejection pattern for θ 2 ∼ 180 • .

Sequential Double Ionization Regime (R Q = 3.2 nm)
Again, decreasing the quantum dot size to access the SDI regime, we note an overall peak increase of the distributions in Figure 7 due to the overall higher probability of ionization in the SDI regime, with the multiplier increasing by two orders of magnitude. In view of the opening of the SDI channel, the total angular distribution is naturally dominated by the pp partial wave, as ionization mostly occurs via the independent ejection of the two electrons, having absorbed a single photon each with no angular momentum transfer required. The domination of the pp ionization channel over the others is even more striking in this case. Another clear difference in the SDI regime is the ratio of peak anti-parallel to parallel ejection has drastically decreased, from ≈2.5 in the DDI regime, to ≈1.2 in the SDI regime. This leads to more symmetric pattern as concerns the probability to have two-electrons in the same (θ 2 ∼ 0 • ) and opposite (θ 2 ∼ 180 • ) directions (given that θ 1 = 0 • is fixed). Considering that the electric field of the radiation alternates its direction evenly, this is what should be expected for independent ejection of the two-electrons. In this sense, the observed reduction ratio of the parallel/back-to-back ejection relative to the DDI regime is again a manifestation of the role of the e-e interaction in the QD's double ionization.

Partial Wave DI Probabilities
One question of research interest is the relative role of the interelectronic interactions and the radiation field in the double ionization process; this certainly depends on the target system and the radiation properties. The size variability of the QD offers the possibility to switch between ionization regimes and keeping the external field constant. In our case, For the larger QD, R Q = 4.6 nm, the DDI mechanism prevails in the double ionization, whereas, for the smaller one, R Q = 3.2 nm, the SDI mechanism is the main double ionization channel. Some insight can be shed by calculating the dynamic development of the four dominant partial DI yields during the pulse at times where the external field vanishes; thus, in Figure 8, we plot the partial yields for the (ss), (sd), (pp), and (dd) channels when the field has completed its 3rd, 6th, 9th, and 12th cycle, for both the larger (left plot) and the smaller (right plot) QD. Figure 8. The partial-wave DI yields for the dominant channels for the R Q = 4.6 nm (left) and the R Q = 3.2 nm QDs (right) plotted versus the field cycles. The total duration was 0.136 ps and the peak intensity was 10 5 W/cm 2 .
It facilitates the discussion if we assume the following symbolic equation of a double-ionization event where the QD absorbs two photons and excites at a state where two electrons are liberated with suitable energies and angular momenta: Generally, the |s 2 CI component is stronger than the |p 2 component (C 1 >> C 2 ). On the other hand, the absolute amplitudes of |C i | 2 , i = 1 − 4 are related with the partial DI yields shown in Figure 8. The procedure followed to calculate the partial DI yields is shown in [21].
Given the above considerations, in the case of the larger QD (where direct DI prevails), the dominant production channel (kp, k p) is a direct absorption of 1-photon from each of the 1s orbitals of the |s 2 of the ground state, ψ g . This transition mode, at rise side of the pulse, generates the (p, p) wavepacket, whereas, later on, the e-e interaction 1/r 12 potential takes over and results to a redistribution of the DI probability to the other partial waves, (ss), (sd), (dd); this redistribution corresponds mainly to angular momentum exchange rather than energy exchange so that both electrons carry the same kinetic energy at the end, as suggested by inspection of the radial distributions in Figure 2 (and noticing that this reflects directly to the joint two-electron kinetic energy spectrum) [21].
It is interesting to comment here that this is in contrast with the atomic Helium (a purely two-electron system) where the trend is to get a flat or U-shaped kinetic energy distribution (in contrast to the peaked one, observed here) depending on the photon energy [23]. At this point, we only add the speculative comment that one reason for this difference could be attributed to the very rapid decay of the Gaussian central potential in QD (∼e −βr 2 relative to the ∼1/r Coulombic potential of the He core).
Similarly, at the beginning of the pulse, the (0,0), (0,2), and (2,2) may be produced directly from weaker component |p 2 of the ground state. The latter therefore suggests that the stronger the correlation in the ground state of the QD the larger the relative production of these channels. Along these lines of thinking, the decrease of the (1,1) population after the ninth cycle (where the laser is almost absent) with the simultaneous increase of (0,2) and (0,0) is the manifestation of e-e correlation interactions, ∼ 1/r 12 (this is a common feature in the case of DI of atomic He, as reported in [14]), while the radiation is the primary interaction during the initial stages of the ionization (when the pulse is rising).
Turning now to the case of the smaller QD (where the SDI is the dominant DI channel), by inspection of the right plot of Figure 8, we observe that the decrease of the (1,1) double ionization channel when the pulse is falling (after the ninth cycle) is not observed. The transition mechanisms discussed above (photon absorption from the weaker components of the ground state at the initial stages and the 1/r 12 e-e interactions at later times) may play a role in the production of the two-electron partial wavepackets of (0,0), (0,2), and (2,2). However, the initially emitted electrons have different kinetic energies (which is reflected in different localization in space) in (1,1), resulting in a weaker 1/r 12 potential and as such production of (0,0), (0,2), and (2,2) via the (1,1) partial wave is no longer effective. Thus, our conclusion in this case is that the most probable DI channel for the latter (weaker) partial waves is the result of direct two-photon absorption from the |p 2 component of the ground state. Such a mechanism may explain the constant larger (more than double) yield of the (0,2) channel relative to the (0,0) given the propensity rule for one-photon electric transitions favouring p → d over p → s (about a factor of 2). We finalize this discussion by commenting that the present observations are in good agreement with the analysis of Barna et al. [24] in the atomic case.

Conclusions
An ab-initio method for the calculation of the joint radial-and angular-electron distributions of the double ionization of a two-electron QD in THz fields was developed and applied to two QDs of different size. The computational method assumes a spherical QD with no other parameters introduced. We developed in detail the non-trivial case of calculating the two-electron continuum QD wavefunction and the associated method of extracting experimentally accessible observables. We related the dominant double ionization mechanisms with the contribution of the electron-electron and radiation interaction potentials in the process and found a similarity in the angular distribution patterns but not in the radial distribution patterns. We believe that the conclusions as they stand can be carried forwards and could be used to interpret experimental ejection patterns for more QDs of more complicated electronic structures. Other aspects of the interaction of a QD with a THz field which could be of interest are the properties of harmonic generation, which was left out of the current discussion, as it involves a different formulation to be developed and constitutes a research project by itself.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

SDI
Sequential Double Ionization DDI Direct Double Ionization e-e electron-electron CI Configuration Interaction