Spin Relaxation in GaAs: Importance of Electron-Electron Interactions

We study spin relaxation in n-type bulk GaAs, due to the Dyakonov–Perel mechanism, using ensemble Monte Carlo methods. Our results confirm that spin relaxation time increases with the electronic density in the regime of moderate electronic concentrations and high temperature. We show that the electron-electron scattering in the non-degenerate regime significantly slows down spin relaxation. This result supports predictions by Glazov and Ivchenko. Most importantly, our findings highlight the importance of many-body interactions for spin dynamics: we show that only by properly taking into account electron-electron interactions within the simulations, results for the spin relaxation time—with respect to both electron density and temperature—will reach good quantitative agreement with corresponding experimental data. Our calculations contain no fitting parameters.


Introduction
Recently, spin coherence in semiconductors has been the focus of both theoretically [1][2][3] and experimentally [4,5] research. A key aim is to achieve a clear understanding of spin decoherence phenomena. This is very important for the emerging field of spintronics, whose goal is to exploit the electron spin, in addition to its charge, within electronic devices. Spin-based devices promise new useful applications in electronics and quantum information [6]; therefore, we wish to control, manipulate and detect electronic spins efficiently, provided that their lifetimes are long enough. The long spin decoherence times measured in semiconductors [7] have made them the subject of intense research. In particular, the n-type bulk GaAs semiconductor has been shown to be a suitable material for spintronics, as it provides the easy availability of high-quality samples and the possibility of using time-resolved optical techniques for exciting and detecting spin-polarized electrons [8].
In n-type bulk GaAs, the main sources of spin relaxation are Elliot-Yafet and Dyakonov-Perel (DP) mechanisms [8], which depend on the spin-orbit interaction. In the Elliot-Yafet mechanism, spin-orbit interaction causes spin depolarization via spin-flips during the carrier scattering events. The spin relaxation due to the DP mechanism follows from the energy splitting, for any non-zero value of the wavevector, of the spin-up and spin-down states. This is present in solids that lack bulk inversion symmetry, like GaAs [9]. This energy splitting gives rise to an effective magnetic field, whose Larmor frequency depends on the carrier's momentum. Therefore, each electronic spin will precess at a different, momentum-dependent rate. In the range of the low-to-medium doping concentrations and high temperatures considered in this paper, DP becomes the dominant mechanism for spin decoherence.
Very recently, the ensemble Monte Carlo (EMC) method has been equipped for dealing with spin transport [10][11][12]. EMC is a stochastic method devised to solve numerically the Boltzmann equation for charge transport in semiconductors [13]. Here, we improve the treatment within EMC of many-body interactions; see Sections 3 and 5. We will use EMC to estimate the effect of electron-electron scattering on the spin relaxation time (SRT) and present results for n-type bulk GaAs at relatively high temperatures (280 ≤ T ≤ 400 K) and low-to-moderate doping concentrations (n i = 10 16 to 2.5 × 10 17 cm −3 ).
Our results display good to very good agreement with the experimental results by Oertel et al. [4] and with no adjustable parameters. In particular, we use the same spin-orbit coupling value, 21.9 eVÅ 3 , suggested in the experimental paper. To our knowledge, this is the first time that EMC simulations can quantitatively reproduce spin-relaxation experimental results, and we will discuss the importance, to this aim, of properly taking into account electron-electron interactions.
Our results also confirm that, in the non-degenerate regime, SRT increases with the electron density, both including or excluding electron-electron scattering [3]. The latter is due to the fact that by increasing the doping concentration, the electron-impurity scattering rate increases and, consequently, the related motional narrowing effect.
Finally, our findings suggest that the prediction made for two-dimensional systems by Glazov and Ivchenko [1,2], that electron-electron scattering slows down the SRT via motional narrowing, can be extended to the three-dimensional case.

Physical Model
We study carrier and spin dynamics in n-type bulk GaAs considering a single parabolic energy band (the central Γ valley), which determines the effective isotropic electron mass m * lab = 0.067m e , m e being the bare electron mass. This approximation is justified, as we do not consider highly energetic electrons excited by a strong electric field, so that inter-valley scattering can be discarded. We include only normal-type scattering events, such as Umklapp processes, that are negligible in direct-gap doped semiconductors. The scattering mechanisms considered are electron-longitudinal acoustic phonon scattering, electron-longitudinal polar optical phonon (POP) scattering, electron-single charged ionized impurity scattering in the Brooks-Herring approach [13,14] and finally electron-electron scattering. Piezoacoustical interaction is not included, because it becomes relevant for GaAs samples only at low temperatures [15]. The scattering rate for the electron-longitudinal acoustic phonon collisions is determined by the acoustic deformation potential [13] in elastic approximation, as inelastic absorption/emission processes are important only at low temperatures [13]. The electron-longitudinal POP scattering rate (Fröhlich interaction) [16] includes absorption and emission processes with a threshold energy of 35 meV, making it the only dissipative process in our model. Phonons are considered at equilibrium at the lattice temperature, T . Following [17,18], the screening of the electron-phonon interactions are not taken into account in the present work.
The scattering rates of electrons from an initial state, |i , to a final state, |f , are calculated to first order, according to the Fermi golden rule: where V is the scattering potential (considered as a perturbation) and E f and E i are the final and initial energy, respectively. By neglecting the exchange and correlation effects, Coulomb interaction between two charges in a homogeneous electron gas is usually estimated using the random phase approximation (RPA) [20], giving rise to an effective screened potential, V sc , whose Fourier components are: Here, (q, ω, T ) is the temperature-dependent dielectric function, v q = e 2 /εq 2 the Fourier components with wavevector q of the bare Coulomb potential and ε is the material dielectric constant, ε = 12.9ε 0 , for GaAs. We approximate (q, ω, T ) with the long-wavelength limit of its static counterpart at finite temperature (q = 0, ω = 0, T ). This is equivalent to the long-wavelength limit of the linearized Thomas-Fermi approximation (LFTA) (q, ω = 0, T ) = 1 + (β 2 T F /q 2 ). We use Dingle's finite temperature LFTA for n-type semiconductors [14,21], which determines the inverse screening length, β T F , from the following equation: Here, n e is the electronic concentration, e the elementary charge, k B the Boltzmann constant, µ is the chemical potential of the electronic ensemble and F j denotes the Fermi-Dirac integral of order j with x ∈ R and Γ Euler's Gamma function. In the non-degenerate regime, Equation (3) reduces to the usual Debye-Hückel inverse screening length. In Figure 1, we plot the values of the screening length λ T F = 1/β T F calculated according to Equation (3) against the electron density. For Equation (3) to hold, the momentum, q, transferred between colliding electron and impurity must remain small [22]. As the electron-impurity process is treated as elastic, q is given by: where v is the magnitude of the electron (group) velocity and θ is the scattering angle. Insofar as the electron-impurity scattering favours small scattering angles, q remains small, and therefore, Equation (3) gives an accurate approximation of the dielectric function in RPA. The electron-impurity scattering angular distribution from our simulations confirms that the LTFA is a good approximation in the regime under investigation, especially at low densities, as the majority of the scattering events happens at small angles; see Figure 2.

Screened Electron-Electron Interaction
Within the RPA, Bohm and Pines [23] have shown that it is possible to split the Coulomb interaction between electrons in two contributions: the first from the collective long-range behaviour (the electron-plasmon interaction) and the second equivalent to a screened Coulomb interaction between individual electrons. In the present work, we consider only the latter, as the electron-plasmon scattering becomes important in GaAs for higher electronic concentrations than considered here [13]. Using LTFA, the electron-electron interaction may then be approximated by the following screened (Yukawa-type) Coulomb potential: where r 1 , r 2 are the spatial coordinates of the colliding electrons. Only binary electron-electron collisions are considered here, as they are the most likely and effective scattering events. The quantum states of mobile electrons should be localized wave packets, but from the perspective of scattering theory [26], the results are equivalent to those obtained using plane waves. Using this property, electron-electron scattering rates in the non-degenerate regime could be calculated using [13]: where V cr is the volume of the crystal, f k is the electron occupation probability (or distribution function), in general, unknown, k 0 is the wave-vector of the colliding electron and the sum runs over all the other electrons in the ensemble.
Within the EMC method, for any given scattering event, once the electron partner of wavevector k, involved in the collision is chosen, the final states, k 0 , k , of the colliding electrons can be determined from the conservation of total energy and momentum and from the scattering angular distribution, P (θ) [13]: Here, g denotes the magnitude of the vector g = k − k 0 , θ is the angle between g and its final state g = k − k 0 and C is a normalization constant: The expression for the scattering rate in Equation (7) arises from our ignorance about the scattering partner in electron-electron collisional events. This explains the presence of the distribution function in Equation (7). However, in our simulations, after having determined the scattering type, we explicitly determine the electron partner from the ensemble. We do so choosing the second electron via a flat distribution within a sphere of radius λ T F centred on the colliding electron; see Section 5. This procedure removes our ignorance about the scattering partner involved in the collisional event and, at the same time, allows us to retain the Bohm and Pines physical picture of individual particles involved in collisions. Then, it follows that we can compute the e-e (electron-electron) scattering rate in a simpler way, using two other ingredients: the Born Approximation and the non-degenerate nature of the system at hand.
First of all, we note that the Fermi golden rule entails first order Born approximation (BA) (usually simply referred to as "Born approximation "). From Equation (1), we can conclude that the scattering rate in BA must not be sensitive to the sign of the potential, i.e., there is no difference between an attractive and repulsive potential, as long as the magnitude of the charges involved is the same. Secondly, if the antisymmetry of the colliding electrons (non-degenerate regime) and the internal structure of the single-ionized impurities may be ignored, we may use the electron-impurity (e-i) scattering rate in the Brooks-Herring approach also for the case of the e-e scattering rate. What now differentiates between e-i and e-e scattering rates are the different effective masses involved in the collision and, consequently, the different reduced masses and energies associated with the particles' relative motion.
All of this considered, assuming a parabolic band and a Yukawa screened potential, the formula for e-e scattering we use in our simulation is the following [13]: where E = E lab /2 and m * = m * lab /2 are the energy and the effective mass of the colliding electron associated with the relative motion, E lab is the energy in the laboratory frame and E β T F is defined by: In the following, we shall use w ee (k) or w ee (E) or w ee (v) when referring to e-e scattering rates written in terms of wavevector, energy or velocity variables, respectively. The wavevectors, energies and velocities of the electrons involved in e-e scattering should be thought of as wavevectors, energies and velocities associated with the relative motion.

The Born Approximation
There are some important consequences about using the BA that we wish to recall. The BA is well satisfied for sufficiently fast carriers assuming a weak scattering potential. It is indeed a high-energy approximation. At low energy (ka 0 1, where k is the magnitude of the colliding electron wavevector and a 0 is the range of the scattering potential), a sufficient condition for the validity of the BA for a central potential (square well) is [24]: where V 0 is the typical strength of a short-range central scattering potential, V . For an attractive potential, the inequality Equation (12) means that the potential, V , is not strong enough to form bound states.
In the case of electron-electron scattering, assuming a screened Coulomb potential, the Equation (12) becomes [25]: with a * B = (4π 2 ε)/(e 2 m * lab ) the effective Bohr radius. The inequality in Equation (13) is not satisfied for the range of densities considered here, where R ∼ > 1 (see Figure 1). Here, the Born series, which solves the Lippmann-Schwinger equation by iteration, may need more terms to converge, and for that reason, BA might give values for the differential cross-section that are not entirely reliable. Indeed, a comparison of differential cross-sections obtained from a Yukawa potential with BA and with exact results obtained by the partial wave method shows that they may be significantly different, depending on the energy, scattering angle, screening and strength of the potential [26]. Unfortunately, the same study shows that results are not improved by including the second term of the Born series for a Yukawa potential, as the differential cross-sections worsen [26].
There are indications that, when BA is not valid, it tends to overestimate the electron-electron total cross-section and, hence, the e-e scattering. Kukkonen and Smith [27], using the method of phase shifts, have found that the electron-electron total cross-section in a metal, like Na (whose average inter-electronic distance, r s , is 3.96 in Bohr radius units), is overestimated by a factor of two, when assuming a scattering potential, like Equation (6), and including the antisymmetry of the wavefunction of the colliding carriers. This improves over previous results, which did not include the antisymmetry and gave an overestimation of a factor of five [27]. The system we are considering is at high temperatures and in a non-degenerate regime; so, the antisymmetry of the wavefunction may be neglected. However, the value of its electron gas parameter, r s , in effective Bohr radius units is similar, with r s > ∼ 1; see Figure 3. We might then expect BA to overestimate e-e scattering also in our case. Clearly, how much the scattering is overestimated is a complicated issue, which strongly relies on the knowledge of the true interelectronic potential.

Ensemble Monte Carlo Method
To model electronic and spin dynamics in GaAs and to estimate the spin relaxation time, we employ the ensemble Monte Carlo method [3,11,12]. This is a semiclassical method in that the simulation has both classical and quantum features. Such a semiclassical approach is valid in the case that the built-in and applied electromagnetic fields are spatially slowly varying.
EMC consists of particles' classical "free-flights", during which the particles may be accelerated by classical forces, interrupted by scattering events, which alter the particles' momentum. The probability of such scatterings and the momentum for each particle following a collision is determined computationally using stochastically generated random numbers.
Among the scattering mechanisms we consider (see Section 2), the scattering of carriers with the longitudinal polar optical phonons is the only source of thermal contact with the lattice. For convenience, we also introduce a fictitious scattering, known in the literature as "self-scattering", which does not affect the particle, but simply ensures that the total scattering rate remains constant [13].
The free flight time, τ , for a particle is calculated as: with r 1 a random number generated from a flat distribution between (0, 1) and Γ = Σ i w max i the total scattering rate, a constant. Here, i enumerates possible scattering types, and w max i is the maximum scattering rate possible for process i. A particle undergoes classical motion for a time, τ ; upon free flight termination, a scattering process is identified for that particle by generating another random number 0 ≤ r 2 ≤ Γ, and the scattering mechanism, i, is chosen when w max i ≤ r 2 < w max i+1 . We then calculate w i (E), the scattering rate associated with the energy, E, of the incoming particle. Given another random number r 3 , if r 3 ≤ w i (E), the particle indeed undergoes a scattering of type i, and the scattered particle is then assigned an outgoing momentum, according to the conservation of energy and momentum for the selected scattering process; else, a "self-scattering" is assumed, and no momentum update is necessary [13].

Electron-Electron Scattering
The electron-electron scattering has to be handled somewhat differently, as it involves two particles. Traditionally, a number of approaches have been used, including treating the electron as scattering with a fictitious partner chosen from a Boltzmann distribution or being allowed to scatter with an actual simulated particle, whose momentum, though, was not updated. This second particle has been usually chosen irrespectively of its distance from the first particle.
In this work, we improve over previous EMC schemes and allow e-e scattering only between electrons that are within one screening length of each other. In our scheme, both electrons scatter, and their momenta are both updated. This approach prevents the unphysical accumulation of energy or momentum prevalent in other methods, as well as the scattering of electrons at opposite ends of the device.
To implement this, we effectively discretized the space into a grid of cubes of one screening length side. We keep track of the number of potential scattering events, which include each particle scattering off any in the same cube or in any of the neighboring cubes. Each time an electron-electron event is required, we choose randomly from each of these potential pairings and check that they are within one screening length of each other, and if they are, we carry out the scattering; if they are not within one screening length of each other, we choose a different electron as the second particle in the scattering.

Thermalization
In order to start collecting data, we need to wait for the electronic system to thermalize to the chosen lattice temperature. To do so, we initially allow the system to evolve for an appropriate time (thermalization run), during which the only source of thermal contact with the environment is provided by polar optical phonons. We note that the thermalization and the corresponding data collection runs always include the same type of scatterings; in particular, they will both include (or not include) electron-electron interactions.
The initial particle configuration (positions and momenta) for the thermalization run is chosen in the following way. The electron positions are generated randomly inside the bulk semiconductor using uniform pseudorandom numbers. Their velocity distribution is generated by choosing the x, y and z components independently from a random Gaussian distribution to reproduce a Boltzmann distribution with an arbitrary temperature of 130 K, which allows us to check that the system correctly relaxes to the lattice temperature.
In order to ascertain that thermalization is reached, we have checked when the energy distribution of the carriers becomes a Boltzmann distribution function corresponding to the lattice temperature. Our simulations show that for the range of parameters of interest in this work and when electron-electron interactions are included, discarding the first 30 picoseconds from the simulation is sufficient to ensure thermalization: in particular, close to room temperature, the thermalization for the runs, including electron-electron interactions, appears to be completed after less than five picoseconds. This confirms the crucial role of electron-electron interactions in the thermalization process [18,19]. We note that when electron-electron interactions are not included, proper thermalization is not reached, as energy gets hardly redistributed within the electron ensemble.
For the results shown in this work, after the thermalization run, we have reset the electronic spins to be fully polarized along one direction, namely, the z-axis, and then started data collections.

Spin Dephasing: The Dresselhaus Term
In bulk n-GaAs at room temperature and for the range of doping densities here considered, the main source of spin relaxation is the Dyakonov-Perel mechanism, a type of spin-orbit interaction. It is due to bulk inversion asymmetry, and it acts as an effective, momentum-dependent magnetic field, via the so-called Dresselhaus Hamiltonian [9,28]: where σ are the Pauli matrices and the Larmor precession frequency vector, Ω(k), is: Here, k i are the wavevector components along the cubic crystal axes, i = x, y, z, and γ so is known as the Dresselhaus coefficient, whose values are determined using different methods. In GaAs, γ so values have been suggested that range from 8.5 to 34.5 eVÅ 3 [29]. The Dresselhaus Hamiltonian causes the electron spins to dephase with respect to each other, as each electron spin in the conduction band precesses with a different Larmor frequency Ω(k), which depends on the specific electron's momentum.

Spin Evolution
In the following, we neglect dipole-dipole interaction between spins. In this way, during free-flight, the spin of each electron undergoes an individual coherent evolution according to the Schrödinger equation.
Initially, each electron spin is assumed to be polarized in the z direction, after which, the spin relaxes via the Dyakonov-Perel mechanism, whereby each spinor wavefunction is acted upon by the time evolution operator generated by the Dresselhaus Hamiltonian, H D , in Equation (15).
The time-evolution operator, U , in spin space for a single particle spinor wavefunction, Ψ, over the time step, δt, is: so that the spinor wavefunction, Ψ (t), at time δt is related to its value at the initial time, t = 0, by: In order to integrate numerically Equation (18), we resort to the Crank-Nicolson (C-N) method [30]. This numerical method integrates by interpolating between two consecutive time steps; hence: where Ψ n = Ψ (nδt) denotes the spinor wavefunction at the n-th-time step. Then, the C-N method leads to the solution: which is correct up to O(δt) 4 . This method is particularly convenient, as the inverse of the spin Hamiltonian can be written analytically, allowing for a significant improvement in computational efficiency compared to the exact solution, with insignificant loss of precision. The C-N method is particularly good for the problem of spin evolution, as it gives a unitary evolution of the spinor wavefunction in time; hence, it conserves its norm. In contrast to the commonly used Heun scheme, the C-N method has the advantage that we do not need to renormalise the spinor wavefunction after each time step. The explicit numerical scheme is: where: h i are the i-components of the effective field, given by the Hamiltonian in Equation (15), and h 2 = 3 i=1 h 2 i . At any given time, we can extract the expectation values of the S x , S y and S z components of the individual electron spin operator, S, to get the probability that the spin is aligned along each direction. Finally, this can be averaged over all spins to give the net spin in any direction. As in this work, we are starting from an electronic ensemble fully polarized in the z direction, we will be interested in looking at the time evolution of the z-component of the total spin, S z,tot . At the n-th time step, this is given by: where N is the number of electrons in the simulation and σ z the z-Pauli matrix.

Estimating the Spin Relaxation Time
Using the above methodology, we are capable of simulating the time evolution of the total electronic spin and of its components in the sample. The quantity of interest to us is the characteristic spin relaxation time of the material. This can be extracted from the time evolution of S z,tot .
We assume that, after a transient period, the spin relaxation behavior in the bulk semiconductor takes the form: It is then possible to fit the data from the simulation of the spin time evolution to such a curve (an example is plotted in Figure 4) and produce values for the parameters, A and B, in the exponential fit. In particular, the parameter, B, has units of s −1 and is identified as the characteristic spin relaxation time of the sample, B = 1/τ s . The spin relaxation curve has a behaviour different from an exponential during the first picoseconds; for example, it starts from a maximum at t = 0, where it then displays a quadratic behavior. We then fit the simulation data exponentially only after this initial transient time. From the analysis of the data in the parameter range we are interested in, we see that neglecting the first 10 ps of the spin depolarization curve is sufficient for this scope.

Results and Comparison with Experiments
In this section, we present and discuss our results for the spin relaxation time and compare them to experimental data.
Apart from assuming an exponential decay of the total spin polarization in the z-direction, we note that our simulations have no fitting parameters. In particular, the spin orbit coupling value used is not fitted, but we use the value suggested by Oertel et al. [4] for their experimental data: γ so = 21.9 eVÅ 3 .
In Figure 5, we plot results from simulations with (τ ee s ) and without (τ no ee s ) electron-electron scattering to examine the effect that the inclusion of electron-electron scattering has on τ s at room temperature. In the same figure, we also plot the experimental results obtained by Oertel et al. in [4] (empty square symbols). When we plot τ s against the range of densities n e = 1 × 10 16 to 2.5 × 10 17 cm −3 , we see that the inclusion of electron-electron scattering causes a net increases of τ s at all densities. Glazov and Ivchenko [2] predicted a similar result in the case of a two-dimensional non-degenerate electron gas in GaAs, explaining it with additional motional narrowing caused by the e-e scattering. Our result suggests that this effect is present also in the three-dimensional case. Figure 5. Results for τ s vs. electronic density calculated with and without electron-electron interaction. Here: N = 25 × 10 3 , T = 300 K and γ so = 21.9 eVÅ 3 . The experimental data from [4] are plotted, as well (empty square symbols). ] no e-e e-e exp Additionally, we notice that the percentage increase of τ s with respect to its non-interacting approximation decreases with increasing density, from about ∼90% to about ∼70%, remaining though always very substantial, even for n e = 2.5 × 10 17 cm −3 . Its absolute increment τ ee s − τ no ee s instead increases with the electronic density.
We observe that, when including e-e interaction, our results for densities 10 16 cm −3 ≤ n e ∼ < 10 17 cm −3 are in very good agreement with the experimental data reproduced in Figure 5.
However, at higher densities, our results for τ ee s start to overestimate the experimental data for τ s , reaching a ∼20% overestimate when n e = 2.5 × 10 17 cm −3 .
We suggest that the overestimate of τ s for n e ∼ > 10 17 cm −3 is due to the BA overestimating the e-e scattering rate, as discussed in Section 4.
We focus now on the effect of temperature on the spin relaxation time.
In order to compare our calculations with other experimental data from [4], we consider the temperature range 280 K ≤ T ≤ 400 K and two (fixed) densities, n e = 2.7 × 10 16 cm −3 and 3.8 × 10 16 cm −3 . In both cases, we will consider interacting carriers.
In Figures 6 and 7, we present our results for n e = 2.7 × 10 16 cm −3 and n e = 3.8 × 10 16 cm −3 alongside the corresponding experimental data (empty square symbols). We find good agreement over the entire temperature range between τ ee s and the experimental data. Figure 6. The spin relaxation time, τ s , versus temperature for the carrier density n e = 2.7×10 16 cm −3 from simulations, including electron-electron interaction, and from the experimental results, as obtained in [4] (empty square symbols). The simulations are done with N = 25 × 10 3 and γ so = 21.9 eVÅ 3 and include e-e scattering. Following [4], data are plotted on a log-log scale. As noted before, the values of the spin orbit coupling for GaAs found in the literature vary greatly [3]; one of the main points in our work is that we do not treat γ so as an adjustable parameter, but simply use the value provided by experimentalists.
In order to let the reader appreciate how valuable this is, and in this respect, how relevant is the good agreement between our data and the experimental ones, in this section, we wish to show how sensible our simulations are with respect to the value of γ so .
In Figure 8, we plot τ ee s for three different values of γ so , all within the range suggested in the literature. It can be seen that by varying γ so , the results for the spin relaxation time would vary within one order of magnitude, and this for the whole density range here considered. 11.0 We think that this is a convincing proof that the very good agreement between our results and the experiments is not accidental, but derives from the improvements we have devised in treating the e-e interaction within the EMC method. These improvements allow us to account properly for the electron-electron interaction within the simulations.

Statistical Analysis of the Spin Relaxation Time Using Coulomb Differential Cross-Sections
We wish to understand better the results relative to the e-e curve in Figure 5. To do so, we focus only on the e-e scattering mechanism, assuming that the other scattering mechanisms give a correct collisional probability. By comparing our calculations with the experimental data, we see that the e-e scattering overestimates τ s at higher concentrations. This may be due to the fact that the e-e scattering itself is overestimated, being that BA is not such a good approximation for low energy carriers; see Section 4. Surprisingly, though, we find very good agreement with the experimental data for densities lower than 10 17 cm −3 , while, as BA worsens at lower densities, we would expect that the SRT curve we obtained from our calculations lies above the experimental curve for the entire range of densities.
To explain this good agreement in the low density limit, we make some general considerations about Coulomb scattering, RPA and screening. Going towards low densities, the RPA starts to break down, which means that in our model, we are no longer allowed to split the e-e interaction into two parts. This can be also understood by looking at r s , as a criterion for the validity of RPA is [32]: From Figure 3, we see that in our system RPA criterion starts to break down for n e ∼ < 1.5×10 17 cm −3 , interestingly a range comparable to the one in which we find agreement between our results for τ ee s and the experimental data. This shows that, at low densities, the potential energy starts to dominate over the kinetic energy. In other words, the long-range component of the Coulomb interaction becomes relevant, and a Yukawa-type potential may be no longer sufficient to realistically describe the inter-electronic potential.
The breakdown of RPA in low electronic densities may affect also the screening length, whose calculation strongly relies on this approximation in our model and, consequently, making the scattering probability less reliable.
The RPA breaking down means that e-e scattering should, in the real system, be more effective. However we still use a Yukawa potential in our calculations; so, as the density decreases, we should be underestimating the e-e scattering and, so, should obtain a τ ee s smaller than the real τ s . However the lower range of density we consider corresponds to the regime where RPA starts to break down (which is compatible with the system r s values), so that the e-e scattering, which results from our simulations, is accidentally correct. We can think of three regimes. In the first, with n e ∼ > 1.5 × 10 17 cm −3 , RPA is appropriate as r s ∼ < 1. BA works well enough as R ∼ 1, and as a result, our simulations overestimate the e-e scattering, i.e., τ ee s > τ s . In the opposite limit (r s 1), RPA is completely inadequate: here, the dominant part of the e-e scattering comes from the long range component of the Coulomb interaction, and if a Yukawa potential would still be used, the simulations would underestimate the e-e scattering; and as a result, τ ee s < τ s . From the trend of r s (see Figure 3), this should happen for densities n e < 10 16 cm −3 , which we do not simulate and which are not realistic, because the system becomes an insulator.
The third regime is intermediate and corresponds to r s of the order of one, with r s > 1. In this regime, RPA has not completely broken down, but the long-range part of the Coulomb interaction starts to become relevant. Using a Yukawa potential then underestimates the e-e interaction, but at the same time, the use of BA (which overestimates the e-e interaction) compensates for this; and we get as a result that τ ee s ∼ τ s . By looking at the values of r s versus density (Figure 3), r s ∼ > 1 for the density range 1 × 10 16 cm −3 ∼ < n e ∼ < 1.5 × 10 17 cm −3 . We indeed find that τ ee s ∼ τ s for the density range 1 × 10 16 cm −3 ∼ < n e ∼ < 1.2 × 10 17 cm −3 (see Figure 5). Another way of looking at the previous considerations is that, for low electronic densities, the system differential cross-section, as described by our simulations, is in some way mimicking a bare Coulomb potential one. Because the later is the exact differential cross-section of the system [31], if our simulations are mimicking it, the related scattering probability would not be overestimated and the quantitative agreement with the experimental result explained.
We will now demonstrate that indeed in our simulations and for the low density range, σ bare ≈ σ Y , with σ bare the bare Coulomb differential cross-section and the σ Y Yukawa differential cross-section in BA. To do so, we will then consider the ratio σ ratio = σ Y /σ bare and determine the conditions, such that σ ratio ≈ 1.
The Yukawa differential cross-section in BA is given by [24]: where θ is the scattering angle associated with the relative motion. The bare Coulomb differential crosssection is obtained from Equation (28) in the limit of β T F → 0, i.e., no screening. Then, σ ratio is: where we have defined the dimensionless quantities ξ = 4E/E β T F = 2E lab /E β T F and s = ξ sin 2 (θ/2). From Equation (29), we note that σ ratio ≤ 1 always. For a given set of energies, E, of the collisional electrons and a threshold value σ * of σ ratio close to unity, there may exist the set of scattering angles I θ * = {θ ∈ [θ * , π] : σ * ≤ σ ratio ≤ 1}. From the angular probability distribution P (θ, E) (see Equation (8)), we can determine the probability, F θ * , that, for a given E, σ * ≤ σ ratio ≤ 1: This integral can be solved analytically, and we get: F θ * (ξ) = cos 2 θ * 2 1 + ξ 4 sin 2 θ * 2 (32) = 4 4 + s(ξ; θ * ) ξ − s(ξ; θ * ) ξ (33) Because the system is at equilibrium, we can use the Boltzmann distribution, f B (E), to weight the function, F θ * , over the whole energy spectrum, giving the probability that an e-e collisional event gives σ * ≤ σ ratio ≤ 1. By using that for fixed σ * , s(ξ; θ * ) becomes a constant, s * , which can be determined from Equation (30); this probability is given by the integral: where α = E β T F /(2k B T ) and the lower integral limit, E * lab , must be determined. This is the the smallest energy for which it is still possible to obtain σ ratio (θ, ξ) as large as σ * . In the limit E lab → E * lab , we have that θ * → π. Imposing the condition σ ratio (θ * , ξ) ≥ σ * , we obtain s(ξ; θ * ) ≥ σ * + √ σ * 1−σ * from which, in the limit θ * → π, we get: The integral, I, is a function of the electronic density through α, so it is possible to compare the probabilities for different electronic densities at T = 300 K. We calculate I for σ * = 0.7 and 0.9; see Figure 9. The results indicate that the e-e collisions with a differential cross-section close to the bare one are more favored at lower densities, which proves our point. The curves in Figure 9 show that, for a density as high as n e = 2.5 × 10 17 cm −3 , only 8% of the total number of carriers would scatter with a differential cross-section, such that 0.7 ≤ σ ratio ≤ 1.

Conclusions
We have improved the treatment of e-e scattering in ensemble Monte Carlo and shown that our method allows one to reproduce, with no fitting parameters, the experimental results for spin relaxation by Oertel et al. We obtain good agreement over the whole range of electron densities and temperatures considered experimentally. Our results show that, in order to achieve quantitative agreement with the experiment, it is crucial to properly include e-e interactions within the simulations. Failure to include many-body interactions leads to greatly underestimating the spin relaxation time.
For the highest electron densities considered, the Born approximation slightly overestimates the e-e scattering rate and, hence, the corresponding scattering cross-section. This implies a higher probability of having a third electron within the scattering cross-section. As future work, we wish to study the importance of this spurious "third-body" effect on spin dynamics in semiconductors and evaluate if an appropriate treatment of it can further improve the agreement with the experimental results.