text:
# Exciton spin relaxation in single semiconductor quantum dots

## Abstract

We study the relaxation of the exciton spin (longitudinal relaxation time $`T_1`$) in single asymmetrical quantum dots due to an interplay of the short–range exchange interaction and acoustic phonon deformation. The calculated relaxation rates are found to depend strongly on the dot size, magnetic field and temperature. For typical quantum dots and temperatures below 100 K, the zero–magnetic field relaxation times are long compared to the exciton lifetime, yet they are strongly reduced in high magnetic fields. We discuss explicitly quantum dots based on (In,Ga)As and (Cd,Zn)Se semiconductor compounds.

The current interest in the manipulation of spins in semiconductors is based on the ability to control and maintain spin coherence over practical length and time scales. Of particular interest for possible applications in quantum computing is the storage of spins in zero–dimensional semiconductor structures. Recent theoretical investigations of the single carrier–spin relaxation in semiconductor quantum dots (QDs) found that the electron spin–flip transitions are very slow with relaxation times of milliseconds and longerKha ; Woods . But since most concepts for quantum information processing involving QDs are based on optical generation, manipulation and read-out of spins, one has to investigate the dynamics of the excitonic spin. In our previous communicationTBK  (motivated by experimental results for InAs QDsPai  at high temperatures $`T40`$ K), we studied the exciton–spin relaxation in QDs related to a second–order process governed by optical phonons. Further detailed theoretical studies of the exciton–spin relaxation in QDs are still missing.

Experimental investigations on the exciton–spin dynamics in QDs refer mostly to the spin–coherence problemGupta ; Woggon ; Hvam  i.e., they determine the transverse relaxation time $`T_2`$. Studies of the longitudinal spin–relaxation time $`T_1`$ are rare since they require strict resonant excitation conditions and/or high magnetic fields. Recent resonant–excitation experiments on InAs QDs indicate no exciton–spin relaxation at low temperatures: the exciton–spin is totally frozen during the radiative lifetime even for high magnetic fields up to 8 TPai . Similar results are reported for CdSe QDs with large lateral dimensions (comparable to twice the exciton Bohr radius of $``$ 5 nm in bulk CdSe crystals)Flis ; Hundt . However, longitudinal spin–relaxation times comparable to the exciton lifetime were deduced from high magnetic field experiments for much smaller CdSe QDsRes . The observed magnetic field dependence (the relaxation rate strongly increases with magnetic fields) points on an acoustic–phonon mediated exciton spin–flip as underlying mechanism. The same conclusion is drawn from the linear temperature dependence of the spin–relaxation rate reported for non–resonant excitation in InGaAs disksGotoh .

We propose in this paper an intrinsic mechanism for excitonic spin–flip transitions at low temperatures resulting from a deformation–induced exchange interaction. Our model calculations are able to qualitatively reproduce the above mentioned trends of the experiments. We study the dependence of the exciton spin relaxation time $`T_1`$ on temperature, magnetic field, and the QD size. Finally, we draw conclusions on the suitability of various materials for the storage of optically generated spins.

Extensive experimental studies have identified the main features of the exciton fine structure in self–organized QDs by means of single–dot spectroscopyBayer ; BK . Such QDs are usually strained and have an asymmetrical shape with a height smaller than the base size. It has been shown that a reduction of the QD symmetry lifts degeneracies among the exciton states and results, in particular, in a splitting of the exciton ground state. Thus, as a consequence of strain and confinement in the growth direction of a QD, the ground states of the heavy–hole (hh) and light–hole (lh) excitons are well-separated and the hh–exciton has the lowest energy. The hh– and lh–exciton quartets are characterized by the projections $`J_z=\pm 1`$, $`\pm 2`$ and $`J_z=\pm 1`$, $`0`$ of the total angular momenta $`J=1,2`$, respectively. The short-range exchange interaction splits the ground states of both hh– and lh–excitons into doublet states (so–called singlet–triplet splitting), as is shown in Fig. 1(a). The lateral anisotropy of a QD leads to a further splitting of the opticaly allowed doublets $`|\pm 1`$ into two levels (labeled $`|X`$ and $`|Y`$) with dipole moments along the two nonequivalent in–plane QD axes. This splitting is due mostly to the long–range exchange interaction and originates from the lateral elongation of the QD’sBayer ; Iv1 ; Ivchenko . Continuous wave single–dot spectroscopy experiments have clearly evidenced the two related, linearly polarized optical transitions e.g., in GaAs interfacial dots Gm1  or self-organized $`\mathrm{InGaAs}`$ QD’sBK . Measured magnitudes of this splitting reach some tens or even hundreds of $`\mu `$eVPai ; BK ; Gm1 . Relaxation processes between the radiative states changes the occupation of the exciton levels and are manifested in a change of the linear degree of polarization of the luminescencePai ; TBK .

In what follows we consider (longitudinal) spin–relaxation processes between the allowed $`|X`$– and $`|Y`$– states of the hh–exciton. In such processes, the hh–exciton stays in the same spatial state and just flips its spin (i.e., the electron and the hole spin flip simultaneously), by emitting or absorbing an acoustic phonon. Since the bare electron/hole–phonon interaction does not contain spin operators, it cannot directly couple the involved exciton states, however, this can occur via the hh and lh exciton mixing due to the interplay of the the short–range exchange interaction and the lattice deformations.

The short–range (isotropic) exchange interaction is given by the Hamiltonian BP

$`_{ex}={\displaystyle \frac{2}{3}}\mathrm{\Delta }_{st}\stackrel{}{\sigma }\stackrel{}{J},`$ (1)
where $`\stackrel{}{\sigma }`$ and $`\stackrel{}{J}`$ are the electron spin and the hole total angular momentum operator, respectively, and $`\mathrm{\Delta }_{st}`$ is the (singlet–triplet) exchange energy. Phonon–induced deformations come into play via the off–diagonal terms in the Bir–Pikus HamiltonianBP

$`(\epsilon )=b{\displaystyle \underset{i}{}}J_i^2\left(\epsilon _{ii}{\displaystyle \frac{1}{3}}\epsilon \right)+{\displaystyle \frac{2}{\sqrt{3}}}d{\displaystyle \underset{i>j}{}}\left[J_iJ_j\right]\epsilon _{ij},`$ (2)
where b and d are the exciton deformation potentials, $`\epsilon _{ij}`$ is the deformation tensor, and $`[J_iJ_j]=J_iJ_j+J_jJ_i`$. (An overall shift of $`a\epsilon `$ of the lh and hh exciton levels is omitted).

In the basis of the hh– and lh–exciton states ($`|1_{hh}`$, $`|+1_{hh}`$ and $`|1_{lh}`$, $`|+1_{lh}`$, respectively), the total Hamiltonian $`=(ϵ)+_{ex}`$ is represented by the following matrix

$``$ $`=`$ $`\left|\begin{array}{cccc}0& 0& \frac{\mathrm{\Delta }_{st}}{\sqrt{3}}& j^{}\\ 0& 0& j& \frac{\mathrm{\Delta }_{st}}{\sqrt{3}}\\ \frac{\mathrm{\Delta }_{st}}{\sqrt{3}}& j^{}& E_{lh}& 0\\ j& \frac{\mathrm{\Delta }_{st}}{\sqrt{3}}& 0& E_{lh}\end{array}\right|.`$ (7)
The origin of energy axis is fixed at the hh–exciton states which are separated from the lh– exciton by $`E_{lh}`$, see Fig. 1(a). The deformation–dependent off–diagonal term is $`j=\sqrt{3}b(\epsilon _{xx}\epsilon _{yy})/2id\epsilon _{xy}`$. Here, we neglect the anisotropic part of the splitting of the allowed doublets, since the energy difference $`E_{lh}`$ is, typically, much larger (of the order of several tens of a meV) than the fine structure energiesBayer . For the same reason, the hh–lh–exciton mixing can safely be regarded as a perturbation. It follows from Eq. (7) that perturbations which mix hh– and lh–exciton states result in a coupling of the $`|1_{hh}`$– and $`|+1_{hh}`$–states and, consequently, lead to the mixing of the $`|X`$– and $`|Y`$– states of the hh–exciton. For QDs which are elongated in the  direction, this coupling is due to vibrations causing $`(\epsilon _{xx}\epsilon _{yy})`$ deformations

$`X||Y=ı{\displaystyle \frac{(j+j^{})\mathrm{\Delta }_{st}}{\sqrt{3}E_{lh}}}ı{\displaystyle \frac{\mathrm{\Delta }_{st}}{E_{lh}}}b(\epsilon _{xx}\epsilon _{yy}).`$ (8)
Thus, a coupling of the $`|X`$– and $`|Y`$–states is a result of an interplay between the short-range exchange interaction and deformations: neither the exchange interaction nor the deformation perturbation $`j`$ alone can couple these states. According to Eq. (8), the matrix element for the relaxation processes among the $`|X`$– and $`|Y`$– states of the hh–exciton is

$`M_{XY}=ı{\displaystyle \frac{\mathrm{\Delta }_{st}}{E_{lh}}}b(\epsilon _{xx}\epsilon _{yy}),`$ (9)
where the expectation value is taken with respect to the spatial part of the exciton wave function. The elastic strain is given by

$`\epsilon _{ij}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{u_i}{r_j}}+{\displaystyle \frac{u_j}{r_i}}\right),`$ (10)
where the phonon displacement field is

$`\stackrel{}{u}={\displaystyle \frac{1}{\sqrt{V}}}\stackrel{}{e}(\stackrel{}{q})\sqrt{{\displaystyle \frac{\mathrm{}}{2\rho \omega _q}}}\left(e^{i\stackrel{}{q}\stackrel{}{r}}a_\stackrel{}{q}+e^{i\stackrel{}{q}\stackrel{}{r}}a_\stackrel{}{q}^+\right).`$ (11)
Here $`\stackrel{}{q}`$, $`\stackrel{}{e}`$ and $`\omega _q=sq`$ are the phonon momentum ($`q=|\stackrel{}{q}|`$), polarization vector and frequency, respectively, V is the normalization volume, $`\rho `$ is the mass density, and $`s`$ is the sound velocity. Upon substitutions given by Eqs. (10) and (11) we obtain for the matrix element in Eq. (9)

$`|M_{XY}|^2`$ $`=`$ $`{\displaystyle \frac{1}{V}}\left({\displaystyle \frac{\mathrm{\Delta }_{st}}{E_{lh}}}\right)^2{\displaystyle \frac{\mathrm{}b^2}{2\rho \omega _q}}(q_xe_xq_ye_y)^2|M_{or}|^2,`$ (12)
$`M_{or}(\stackrel{}{q})`$ $`=`$ $`\mathrm{\Phi }_0|e^{i\stackrel{}{q}\stackrel{}{r_e}}+e^{i\stackrel{}{q}\stackrel{}{r_h}}|\mathrm{\Phi }_0,`$ (13)
where indices $`e`$ and $`h`$ stand for electrons and holes, $`M_{or}`$ is the orbital part of the matrix element, and $`\mathrm{\Phi }_0(\stackrel{}{r_e},\stackrel{}{r_h})`$ is the exciton ground state envelope wave function.

The spin–relaxation rates accompanied by phonon emission and absorption are given by Fermi’s golden rule

$`{\displaystyle \frac{1}{\tau _{}}}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{\stackrel{}{q}}{}}{\displaystyle \underset{\stackrel{}{e}}{}}|M_{XY}|^2\delta (\mathrm{}\mathrm{\Omega }\mathrm{}\omega _q)(N_{\omega _q}+1),`$
$`{\displaystyle \frac{1}{\tau _{}}}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{\stackrel{}{q}}{}}{\displaystyle \underset{\stackrel{}{e}}{}}|M_{XY}|^2\delta (\mathrm{}\mathrm{\Omega }\mathrm{}\omega _q)N_{\omega _q},`$ (14)
where $`\mathrm{}\mathrm{\Omega }`$ is the energy splitting between the $`|X`$– and $`|Y`$–states (see Fig. 1(b)), and $`N_{\omega _q}=1/(e^{\mathrm{}\omega _q/k_BT}1)`$ is the thermal phonon distribution function.

We discuss now what determines the value of the orbital part of the matrix element Eq. (13). Because of energy conservation, $`|M_{or}||M_{or}^e+M_{or}^h|`$ has to be calculated for the phonon momentum $`q=\omega _q/s\mathrm{\Omega }/s`$, while the characteristic extent of the orbital wave function $`\mathrm{\Phi }_0`$ is, roughly speaking, the lateral size $`L`$ of the QD. (For flat QDs, the QD hight $`L_zL`$.) Due to the oscillatory behavior of the exponential functions in $`M_{or}`$, we have $`|M_{or}^{e,h}|1`$, if $`\mathrm{\Omega }L/s1`$. To illustrate this approximation, we consider the model of a strongly confined QD with the shape of an parallelepiped and infinite barrier. Then the exciton states are reasonably well approximated by noninteracting electron–hole pair states. In addition, the single carrier–phonon matrix element $`|M_{or}^e(\stackrel{}{q})|=|M_{or}^h(\stackrel{}{q})||M_{or}(\stackrel{}{q})|`$ separates in the $`x,y,z`$ coordinates. For a confined direction (say, the $`x`$–direction), the following expression holds

$`|M_{or}(q_x)|=\left|{\displaystyle \frac{\mathrm{sin}Q_x}{Q_x}}+{\displaystyle \frac{Q_x\mathrm{sin}Q_x}{Q_x^2\pi ^2}}\right|,`$ (15)
where $`Q_x=q_xL_x/2`$ and $`L_x`$ is the QD size in the $`x`$–direction. As seen from Eq. (15), the matrix element $`|M_{or}(q_x)|`$ is smaller than unity for any $`q_x0`$, but it approaches unity for small $`q_xL_x1`$ ( $`|M_{or}(q_x)||1Q_x^2/\pi ^2|`$). With increasing $`q_xL_x1`$, $`|M_{or}(q_x)|`$ decreases rapidly ( $`|M_{or}(q_x)||\mathrm{sin}Q_x/Q_x|`$). Thus, if we replace the orbital matrix element in Eq. (12) by unity, $`|M_{or}(\stackrel{}{q})|^21`$, the exciton–spin relaxation rates Eq. (14) will be overestimated. Note that for typical values of $`\mathrm{\Omega }`$ and $`s`$ we have $`(s/\mathrm{\Omega })40`$ nm, $`L`$ is usually about $`10`$$`20`$ nmtv , so that the above mentioned assumption $`\mathrm{\Omega }L/s1`$ is reasonable. As a result, a lower limit for the relaxation times is given by

$`\tau _{}`$ $`=`$ $`{\displaystyle \frac{9}{4}}{\displaystyle \frac{\mathrm{}\rho s^2}{b^2}}\left({\displaystyle \frac{E_{lh}}{\mathrm{\Delta }_{st}}}\right)^2\left({\displaystyle \frac{s}{\mathrm{\Omega }}}\right)^3\left(1e^{\mathrm{}\mathrm{\Omega }/k_BT}\right),`$ (16)
$`\tau _{}`$ $`=`$ $`\tau _{}e^{\mathrm{}\mathrm{\Omega }/k_BT}.`$ (17)
In order to obtain a numerical estimate, we use typical parameterstv , together with an estimated hh–lh–exciton splitting of $`E_{lh}=10`$ meV, i.e., $`E_{lh}/\mathrm{\Delta }_{st}=50`$. As a result, we get for the spin relaxation time scale $`T_1(\tau _{}^1+\tau _{}^1)^11500`$ ns (at $`T=10`$ K) and $`T_1150`$ ns (at $`T=100`$ K). Thus, for typical quantum dots and low temperatures, the zero–magnetic field relaxation times are very long compared to the exciton lifetimes of $`1`$ ns. This is in qualitative agreement with resonant–excitation experiments at low temperaturesPai ; Flis . Please note, that because of the small interlevel splitting $`\mathrm{}\mathrm{\Omega }0.1`$ meV, the relaxation times $`\tau _{}`$ and $`\tau _{}`$ are already comparable at $`T2`$ K. For higher temperatures $`T_1`$ is inversely proportional to the temperature. But for $`T100`$ K, the interaction with LO phonons via a second–order process will be the dominant relaxation mechanismPai ; TBK .

The calculated spin relaxation times according to Eqs. (16,17) are large partly because of the small phonon density of the states (which is $`(\mathrm{\Omega }/s)^2`$) at the scale of the interlevel splitting $`\mathrm{}\mathrm{\Omega }`$. But, considerably large interlevel splittings of $`\mathrm{}\mathrm{\Omega }_B=g\mu _BB1.5`$ meV (where $`g`$ is the exciton effective $`g`$–factor) arise in high–field magneto–optical experiments. Here an enhancement of spin–relaxation rate is expected. The spin flip in Faraday geometry ($`\stackrel{}{B}z`$) occurs between the $`|1_{hh}`$– and $`|+1_{hh}`$– states of the hh–exciton (see Fig. 1(c))Bayer . Consequently, exciton–spin relaxation processes between the Zeeman sublevels are reflected in a change of the circular degree of polarization of the luminescence.

For magnetic fields and zero temperature, a lower limit of the spin relaxation time is given by

$`T_1(B)`$ $``$ $`\tau _{}(B)={\displaystyle \frac{9}{4}}{\displaystyle \frac{\mathrm{}\rho s^2}{\stackrel{~}{b}^2}}\left({\displaystyle \frac{E_{lh}}{\mathrm{\Delta }_{st}}}\right)^2\left({\displaystyle \frac{s}{\stackrel{~}{\mathrm{\Omega }}_B}}\right)^3,`$ (18)
$`\stackrel{~}{b}^2`$ $`=`$ $`b^2+{\displaystyle \frac{\mathrm{\Omega }_B}{3\sqrt{\mathrm{\Omega }_B^2+4\mathrm{\Omega }^2}}}d^2,`$ (19)
$`\stackrel{~}{\mathrm{\Omega }}_B`$ $`=`$ $`\mathrm{\Omega }+\mathrm{\Omega }_B.`$ (20)
The relaxation time in Eq. (20) varies as the square of the ratio $`E_{lh}/\mathrm{\Delta }_{st}`$ and depends strongly on the QD size because both the splitting $`E_{lh}`$ and the exchange energy $`\mathrm{\Delta }_{st}`$ are size–dependentTag ; rem . If the splitting $`E_{lh}`$ is due mostly to the strain (large lattice mismatch), $`T_1`$ changes with the lateral and vertical size as $`L^4`$ and $`L_z^2`$, respectively. For small lattice mismatch, $`E_{lh}`$ is due mostly to the vertical confinement and, therefore, $`T_1L^4L_z^2`$. In any case, the relaxation time $`T_1`$ strongly decreases when decreasing the lateral size $`L`$ of a QD. This result can explain qualitatively the observations for the CdSe QDs of different lateral sizes in high magnetic fields. As noted above, no spin relaxation was found for large CdSe QDs (with $`L`$ about 10 nm)Hundt , while for small CdSe QDs (with $`L`$ about 3 nm - 5 nm) an efficient exciton–spin relaxation was observed Res .

To obtain numerical estimates for the spin relaxation in magnetic fields, we use the experimental values for $`E_{lh},\mathrm{\Delta }_{st}`$ and $`g`$, and we take the typical valuestv  for other factors in Eq. (18). For InAs/GaAs QDs, we use $`\mathrm{\Delta }_{st}=0.2`$ meV, $`g=3`$BK  and $`E_{lh}=30`$ meV ($`E_{lh}`$ of several tens of a meV is reported in Ref.Bayer ). As a result we get a lower limit for the relaxation time $`T_113/B^3`$ ($`\mu `$s T3) which is much larger than the exciton lifetime scale $`1`$ ns even at high $`B`$ (e.g., $`T_125`$ ns at $`B=8`$ T). Indeed, no spin relaxation is found for low temperature experiments in high magnetic fields Pai .

This case is different for small CdSe/ZnSe QDs investigated in Ref.Res ; Verbin . Here, values of $`\mathrm{\Delta }_{st}1`$ meV, $`g2`$ and $`E_{lh}40`$ meV are reported. As a result we calculate $`T_15`$ ns at $`B=6`$ T which is close to values of 2.5 ns - 3.4 ns deduced from the experimentsRes . For large CdSe/ZnSe QDs investigated in Ref.Hundt  no data on the exchange energy $`\mathrm{\Delta }_{st}`$ are reported but we still can estimate the relaxation rate. The lateral size of the QDs in this case is about a factor of 2.5 larger than in Res . The consequence is a much smaller exchange splitting $`\mathrm{\Delta }_{st}0.2`$ meV rem  which results in much longer relaxation times of $`100`$ ns (at $`B=6`$ T). As already mentioned, no spin relaxation is observed in the experiments of Ref.Hundt .

Let us now consider the magnetic field dependence of the spin–relaxation time which was measured in small CdSe/ZnSe QDs Res . The experimental results indicate a roughly anti–proportional decrease $`T_1B^1`$ up to 8 T. As seen from Eqs. (18,20), $`T_1`$ scales as $`B^1`$ for small $`B`$ ($`\mathrm{\Omega }_B<\mathrm{\Omega }`$), in agreement with the experiment. But it would depend much stronger on the magnetic field (as $`B^3`$) for large $`B`$. This difference to the experiment can originate from the treatment of the orbital matrix element, which up to now has been fixed as 1. But a factor of $`1/|M_{or}|^2`$ has to be accounted for explicitly in Eq. (18) for high magnetic fields. The reason is, that the relevant phonon momentum $`q(B)\mathrm{\Omega }_B/s`$ is linearly proportional to $`B`$. Thus, $`1/|M_{or}|^2`$ increases with magnetic field and becomes important at $`q(B)L1`$. This is the case e.g., for $`B>5`$ T for CdSe QDs with a lateral size of $`L4`$ nm and is the regime of the experiments in Res . So we can confirm the general trend of the experimentRes , namely that $`T_1`$ decreases with magnetic field. An analytical dependence of $`T_1`$ on $`B`$ (at high $`B`$) would require an explicite calculation of $`|M_{or}|^2`$ with the QD exact parameters, which for the QDs used in Ref.Res , unfortunately, are not known.

There are many experiments, which determine the polarization degree of the luminescence from QDs under nonresonant excitation conditions. We just want to discuss the results of Ref.Gotoh  as a typical representative. Rather short polarization decay times on the order of 1 ns have been found here. But definite conclusions on the underlying spin–relaxation mechanisms are not possible in a scenario where the electron–hole pairs relax through a multitude of barrier states and/or excited quantum–dot states. To identify spin–relaxation processes, the splitting of the exciton ground state has to be large (as for large magnetic fields) or strictly resonant excitation has to be used (see Pai ).

The above discussions show that a direct quantitative comparison of our proposed relaxation model to experimental data is at the moment difficult at best. One obvious reason is that the relaxation times, once they are much larger than the exciton lifetime, cannot be quantified experimentally. The second reason is, that only very few experiments under strict resonant excitation and/or in high magnetic fields have been performed. But still, there are a couple of experimental trends which are explained consistently by our model. The value and magnetic–field dependence of $`T_1`$ in small CdSe QDs can be reproduced. We correctly predict the experimental findings that long relaxation times are expected for InAs QDs as well as large CdSe QDs at low temperatures even in strong magnetic fields.

The experimental differences found for InAs vs CdSe QDs and for large vs small CdSe QDs strongly support our proposed model. Let’s consider alternative spin–relaxation mechanisms within the radiative doublet. Processes relying on the exciton motion like for quantum-well excitons Maialle  are suppressed in QDs. Assume, the electron and hole would individually flip their spins. In contrast to our proposed mechanism, this is a higher–order process involving intermediate dark exciton states (the excitons finally have to return to radiative states to contribute to a depolarized luminescence). The relaxation time here is determined by the longest spin–flip time for a single carrier (electron or hole). The spin flip of the electron has to rely on mechanisms based on the spin–orbit coupling. These relativistic effects have been shown by Khaetskii and Nazarov Kha  to be inefficient in QDs. Finally, the spin–orbit interaction and carrier–phonon coupling are not different enough in InAs and CdSe or for large and small QDs to explain that the spin relaxation times change by orders of magnitude. This dramatic change found in the experiments is, however, readily explained by the strongly differing exchange splitting in InAs and CdSe and its significant enhancement in small dots. Since the exchange splitting enters quadratically in our proposed relaxation scheme, we are quite convinced to have identified the most prominent spin–relaxation process for QDs.

Finally, we want to discuss different material systems with respect to their suitability for spin storage. The relevant parameters for several semiconductor compounds forming QD’s and their matrices are presented in Table 1. Among the different material parameters in Eq. (20), the characteristic energies $`E_{lh}`$ and $`\mathrm{\Delta }_{st}`$ are strongly affected by the structural composition. Thus we can only give general trends. We will consider QDs with typical sizes of about $`L_z2`$$`5`$ nm and $`L10`$$`20`$ nm which are in the strong vertical confinement regime. Moreover, $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}`$As QDs are also strongly confined in the lateral plane, and the exchange splitting $`\mathrm{\Delta }_{st}`$ in these QDs significantly exceeds its bulk value of $`\mathrm{\Delta }_{st}^b`$rem . In CdSe QDs the bulk value of the exchange splitting is already quite large. In CdSe QDs with lateral sizes larger than the exciton radius of $`a_{ex}5`$ nm the exchange splitting is additionally affected by the vertical confinement. Since the exchange interaction in (In,Ga)As crystals is much weaker than in CdSe crystals, $`\mathrm{\Delta }_{st}`$ is of the same order of magnituderem1  in strongly confined $`\mathrm{In}_\mathrm{x}\mathrm{Ga}_{1\mathrm{x}}`$As QDs and large, weakly confined CdSe QDs. Considering additionally the strain–induced lh–hh splitting $`E_{lh}`$, which has to be large to suppress spin relaxation, InAs QDs have the longest relaxation time $`T_1`$.

Less suitable for spin storage is the GaAs/Al0.3Ga0.7As system since its lattice mismatch is very small. Thus, among the QD structures listed in Table I, the shortest $`T_1`$ relaxation times can be expected for strongly confined GaAs/Al0.3Ga0.7As QDs. Likewise, small $`T_1`$ times are predicted for small CdSe QDs since in the strong lateral confinement regime, $`\mathrm{\Delta }_{st}`$ is considerably enhanced.

In conclusion, the exciton–spin relaxation within the radiative doublet of the exciton ground state in single asymmetrical QDs is studied. As a possible intrinsic mechanism for such a process, the exciton spin–acoustic phonon coupling via the strain–dependent short range exchange interaction is proposed. For zero-magnetic fields and low temperatures, the calculated $`T_1`$–relaxation times for typical QDs are long compared to the exciton lifetime. A strong reduction of the relaxation times occurs in QDs in high magnetic fields. Nevertheless, numerical estimates for InAs QD’s and large CdSe QDs give large values for the $`T_1`$–relaxation times even for high $`B`$ 8 T (up to a few of tens of nanoseconds). In addition, the relaxation time $`T_1`$ strongly decreases in strongly confined (in the lateral plane) QDs. The $`T_1`$–relaxation times estimated for small CdSe QDs in high magnetic fields reduce to the nanosecond scales. These predictions of our model are in qualitative agreement with experimental findings. Within the considered mechanism, we conclude that InAs QDs and CdSe QDs with typical lateral sizes of $`L10`$$`20`$ nm display a weak exciton spin relaxation, whereas for typical GaAs/Al0.3Ga0.7As QDs and small CdSe QDs the $`T_1`$ times are expected to be rather short.

This work was supported by the Center for Functional Nanostructures (CFN) of the Deutsche Forschungsgemeinschaft (DFG) within project A2.


format: text

text:
# Impurity and spin effects on the magneto-spectroscopy of a THz-modulated nanostructure

## I Introduction

The interest in the properties of radiation induced excited states in bulk semiconductors,Planken et al. (1995); Nie et al. (2003) or a two-dimensional electron gas (2DEG),Mani et al. (2002) and semiconductor nanostructuresLorke et al. (2000) has grown in the last years parallel to advances in high precision experimental techniques and sample making skills. Increased brilliance of coherent far-infrared (THz) synchrotron radiationAbo-Bakr et al. (2003) and sensitive photoconductivity measurements techniquesPlanken et al. (1995); Nie et al. (2003) will certainly intensify research activities on electronic systems driven far from equilibrium and time-dependent phenomena in nanostructures.

We model and investigate here the time evolution of a short nanoscale quantum wire or an elongated quantum dot with one embedded impurity after the system with several electrons has been irradiated by a short intense THz pulse. Corresponding atomic systems have been investigated,Chu and Chu (2001a) but the present model of an artificial two-dimensional (2D) molecule in a GaAs structure allows us to consider the effects of an external magnetic field $`B<2`$ T, that is commonly available in the labs and has large effects on the 2DEG in a AlGaAs/GaAs heterostructure, but very small effects on usual atomic systems.

Having seen that an intensive THz pulse can cause nonlinear effects in quantum dotsPuente et al. (2001) and ringsGudmundsson et al. (2003) and possibly change the binding energy of impurities in bulk materialNie et al. (2003) we study its effects in the present system on the binding energy of the impurity in presence of other interacting electrons, and ask ourselves if it is possible to excite collective spin modes in the system in an external magnetic field. THz radiation with the equivalent energy of few meV is simply on an appropriate scale with the energy spectrum of our system determined by its geometry (finite size) and the host material of GaAs.

Modeling time-dependent phenomena in nanosystems in a magnetic field driven out of equilibrium does not leave many options to a researcher. The Coulomb interaction between the electrons is usually of paramount interest together with the spin structure of the system. This can be emphasized by its finite size even though the effective g-factor may be small as is the case in GaAs. The well established time-dependent Hartree-Fock approximation can be adequate in many cases,Puente et al. (2001) but one would hope for more accurate spin structure delivered by a time-dependent local spin-density approximation (LSDA).Puente and Serra (1999) It has to be kept in mind that both the inclusion of the magnetic field in a LSDAQian et al. (2003) and its extensions to dynamic systemsAryasetiawan and Gunnarsson (2002) are nontrivial tasks and are not fully completed to date. In our approach we shall use the ground state parameterization of the exchange and correlation functionals and ignore the small corrections due to current terms in the functionals in a magnetic field.Vignale and Rasolt (1987)

Commonly, both in calculations for excited atoms and moleculesBandrauk and Lu (2000); Chu and Chu (2001a, b) or nanostructures,Puente and Serra (1999); Puente et al. (2001) a spatial grid has been employed. In this work we resort to recently developed grid-free implementation of LSDAZheng and Almlöf (1993); Glaesemann and Gordon (1998, 1999); Berghold et al. (1998) in the hope to be able to describe complex dynamical systems for a longer period of time. The grid-free methods, as will become clear below, result in the transformation of the relevant equations to compact matrix equations of not too large size. On the other hand the matrix elements are usually much more complicated than in grid methods, even though most of them are known analytically. This compact matrix structure of the problem lends itself very well to the matrix and vector handling of modern computer languages and facilitates parallelization that does neither require lare bandwidth nor high communication speed between the nodes of a CPU cluster.

## II Model

### II.1 The ground state

To calculate the time evolution of the system requires the knowledge of the initial state, its ground state, at a particular time $`t=t_0`$, which within the LSDA is described by the Kohn-Sham equations

$$H|\alpha )=(H_0+H_\sigma +V_\varphi +H_{\mathrm{int}})|\alpha )=\epsilon _\alpha |\alpha ),$$
(1)
where $`H_0`$ is the Hamiltonian for the 2DEG in the perpendicular magnetic field $`𝐁=B\widehat{𝐳}`$, that is confined by the parabolic potential $`V_{\mathrm{conf}}(r)=m^{}\omega _0^2r^2/2`$. A more general shape of the system is attained through the additional confinement potential $`V_\varphi `$ breaking the rotational symmetry of $`H_0`$ and modifying the radial confinement

$$V_\varphi (𝐫)=\frac{1}{2}m^{}\omega _0r^2\underset{p=1}{\overset{p_{\mathrm{max}}}{}}\alpha _p\mathrm{cos}(p\phi )+V_0\mathrm{exp}(\gamma r^2),$$
(2)
where $`p_{\mathrm{max}}`$ is an even integer to ensure confinement. The eigenstates $`|\alpha `$ of the Fock-DarwinFock (1928) Hamiltonian $`H_0`$ are used as a mathematical basis in which the eigenstates $`|\alpha )`$ of $`H`$ are expanded, as a result Eq. (1) can be transformed into a matrix eigenvalue equation. A suitable truncation of the basis $`\{|\alpha \}`$, that can be varied to improve the accuracy, then leads to a numerically tractable eigenvalue problem.

In this LSDA mean-field type of an approximation the effective interaction $`H_{\mathrm{int}}`$ is composed of the direct Coulomb repulsion potential, or the Hartree term

$$V_\mathrm{H}(𝐫)=\frac{e^2}{\kappa }\frac{n(𝐫^{})}{|𝐫𝐫^{}|}𝑑𝐫^{},$$
(3)
and the potential describing the effects of the exchange and correlation

$$V_{xc,\sigma }(r,B)=\frac{}{n_\sigma }(nϵ_{xc}[n_{},n_{},B])|_{n_\sigma =n_\sigma (r)},$$
(4)
where $`n_\sigma (r)`$ is one of the two spin ($`\sigma =,`$) components of the total electron density $`n=n_{}+n_{}`$. The Zeeman energy of the electrons is included in $`H_\sigma `$ with the appropriate effective g-factor $`g^{}`$ and the Bohr magneton, $`\mu _B`$, $`H_\sigma =\pm (1/2)g^{}\mu _BB`$.

In a large 2DEG the natural length and energy scales are the magnetic length $`l=\sqrt{\mathrm{}c/(eB)}`$ and the cyclotron energy $`\mathrm{}\omega _c=\mathrm{}eB/(m^{}c)`$, and instead of the electron density it is convenient to define the local filling factor $`\nu (𝐫)=2\pi l^2n(𝐫)`$. In the confined system described by $`H_0`$ the appropriate length scale is $`a=l\sqrt{\omega _c/\mathrm{\Omega }}`$, an “effective magnetic length”, where $`\mathrm{\Omega }=\sqrt{\omega _c^2+4\omega _0^2}`$ is the new characteristic frequency of the system. We shall thus use a modified filling factor appropriate for a finite system

$$\stackrel{~}{\nu }(𝐫)=2\pi a^2n(𝐫),$$
(5)
and the polarization $`\zeta (𝐫)=[n_{}(𝐫)n_{}(𝐫)]/n(𝐫)`$ as new variables to express the exchange and correlation potentials (4) in. This will allow us to have a unified description of the system in a finite or vanishing magnetic field.

In these new variables the exchange and correlation potentials areLubin et al. (1997)

$`V_{xc,}={\displaystyle \frac{}{\stackrel{~}{\nu }}}(\stackrel{~}{\nu }ϵ_{xc})+(1\zeta ){\displaystyle \frac{}{\zeta }}ϵ_{xc}`$
$`V_{xc,}={\displaystyle \frac{}{\stackrel{~}{\nu }}}(\stackrel{~}{\nu }ϵ_{xc})(1+\zeta ){\displaystyle \frac{}{\zeta }}ϵ_{xc}.`$ (6)
In the magnetic field the exchange and correlation energy density is interpolated between the infinite and the no-field limits asKoskinen et al. (1997)

$`ϵ_{xc}^B(\stackrel{~}{\nu },\zeta )=ϵ_{xc}^{\mathrm{}}(\stackrel{~}{\nu })e^{f(\stackrel{~}{\nu })}+ϵ_{xc}^0(\stackrel{~}{\nu },\zeta )(1e^{f(\stackrel{~}{\nu })}),`$ (7)
with $`f(\stackrel{~}{\nu })=(3\stackrel{~}{\nu }/2)+7\stackrel{~}{\nu }^4`$ and $`ϵ_{xc}^{\mathrm{}}(\stackrel{~}{\nu })=0.782\sqrt{\stackrel{~}{\nu }}e^2/(kl)`$. The low-field energy functional is interpolated between the spin-polarized and the unpolarized limitsvon Barth and Holm (1996)

$`ϵ_{xc}^0(\stackrel{~}{\nu },\zeta )=ϵ_{xc}(\stackrel{~}{\nu },0)+f^i(\zeta )\left[ϵ_{xc}(\stackrel{~}{\nu },1)ϵ_{xc}(\stackrel{~}{\nu },0)\right]`$
$`\text{with}f^i(\zeta )={\displaystyle \frac{(1+\zeta )^{3/2}+(1\zeta )^{3/2}2}{2^{3/2}2}}.`$ (8)
The exchange and correlation contributions are then separated $`ϵ_{xc}(\stackrel{~}{\nu },\zeta )=ϵ_x(\stackrel{~}{\nu },\zeta )+ϵ_c(\stackrel{~}{\nu },\zeta )`$, and the exchange contribution for the unpolarized system is expressed as $`ϵ_x(\stackrel{~}{\nu },0)=[4/(3\pi )]\sqrt{\stackrel{~}{\nu }}e^2/(kl)`$, while for the polarized one it is $`ϵ_x(\stackrel{~}{\nu },1)=[4/(3\pi )]\sqrt{2\stackrel{~}{\nu }}e^2/(kl)`$. We use the parameterization of Ceperley and Tanatar for the correlation contributionTanatar and Ceperley (1989)

$`ϵ_c(\stackrel{~}{\nu },\zeta )=a_0{\displaystyle \frac{1+a_1x}{1+a_1x+a_2x^2+a_3x^3}}E_{Ryd}^{},`$ (9)
where $`x=\sqrt{r_s}=(2/\stackrel{~}{\nu })^{1/4}(a/a_B^{})^{1/2}`$, and $`a_B^{}`$ is the effective Bohr radius. The optimized values of the correlation parameters $`a_i`$ have been found by Ceperley and Tanatar using a Monte Carlo model of a 2DEG.Tanatar and Ceperley (1989)

So stated, the $`V_{xc,\sigma }`$’s are complicated functionals of the the spin densities, or $`\stackrel{~}{\nu }`$ and $`\zeta `$, and thus a common method is to solve the resulting Kohn-Sham equations on a spatial grid. We instead follow the alternative example of a number of researchers to use the mathematical basis $`\{|\alpha \}`$ to cast the LSDA functionals into matrix expressions in order to implement a grid-free LSDA.Zheng and Almlöf (1993); Glaesemann and Gordon (1998, 1999); Berghold et al. (1998) Briefly, the method relies on the following steps: In the end we need the matrix elements of the exchange and correlation potentials, $`\alpha |V_{\mathrm{xc},\sigma }|\beta `$, so we start by constructing $`\alpha |\stackrel{~}{\nu }|\beta `$. The polarization $`\zeta `$ has to be treated as a combination of two functions, $`f=\stackrel{~}{\nu }_{}\stackrel{~}{\nu }_{}`$ and $`g=\stackrel{~}{\nu }^1`$. Due to the completeness of the basis we then have $`\alpha |\zeta |\beta =_\gamma \alpha |f|\gamma \gamma |g|\beta `$, and to evaluate the latter matrix elements we need to find a unitary matrix $`𝐔`$ that diagonalizes the matrix $`\stackrel{~}{\nu }`$, i.e. $`\stackrel{~}{\nu }=𝐔\mathrm{𝐝𝐢𝐚𝐠}(\lambda _1,\mathrm{},\lambda _n)𝐔^+`$, where $`\lambda _i`$ are the eigenvalues of $`\stackrel{~}{\nu }`$. All functions of the matrix $`\stackrel{~}{\nu }`$ are now calculated according to

$$f[\stackrel{~}{\nu }]=𝐔\mathrm{𝐝𝐢𝐚𝐠}(f(\lambda _1),\mathrm{},f(\lambda _n))𝐔^+.$$
(10)
Corresponding methods are used for more complex functionals of $`\stackrel{~}{\nu }`$ and $`\zeta `$.

The matrix elements of the electron density $`n`$ or the effective filling factor $`\stackrel{~}{\nu }`$ can be written in terms of a generalized one-electron overlap integral and the matrix elements of the density operator $`\rho `$

$$\alpha |\stackrel{~}{\nu }|\beta =\underset{p,q}{}\rho _{qp}𝑑𝐫\varphi _\alpha ^{}(𝐫)\varphi _p^{}(𝐫)\varphi _q(𝐫)\varphi _\beta (𝐫).$$
(11)
In the ground state the density matrix is constructed from the expansion coefficients of electronic states in terms of the basis states, the occupation of the states expressed by the Fermi distribution, and the chemical potential $`\mu `$ ensuring the conservation of the number of electrons, $`\rho _{qp}=_\gamma f(\epsilon _\gamma \mu )c_{\gamma p}^{}c_{\gamma q}`$. We use here the density matrix $`\rho _{qp}`$ in anticipation of the situation when the system will be subjected to a strong external time-dependent perturbation and quantities like the single-electron energy spectrum, and the occupation described by the Fermi distribution will have no meaning. The analytic expression for the generalized one-electron overlap integral can be found in Appendix. The matrix elements of the direct Coulomb interaction and the confinement potential have been published elsewhere.Magnúsdóttir (1999); Vasile (2002)

### II.2 Time-dependent excitation

At $`t=t_0`$ the Hamiltonian of the system acquires time dependence, $`H(t)=H+W(t)`$, caused by an external perturbation of a finite duration

$`W(t)`$ $`=`$ $`V_tr^{|N_p|}\mathrm{cos}(N_p\varphi )\mathrm{exp}(sr^2\mathrm{\Gamma }t)`$ (12)
$`\mathrm{sin}(\omega _1t)\mathrm{sin}(\omega t)\theta (\pi \omega _1t).`$
Since the perturbation is of arbitrary strength we have to assume the task of solving the equation of motion for the density operator

$$i\mathrm{}d_t\rho (t)=[H+W(t),\rho (t)].$$
(13)
The structure of this equation is inconvenient for numerical evaluation so we resort instead to the time-evolution operator $`T`$, defined by $`\rho (t)=T(t)\rho _0T^+(t)`$, which has simpler equation of motion

$`i\mathrm{}\dot{T}(t)`$ $`=`$ $`H(t)T(t)`$ (14)
$`i\mathrm{}\dot{T}^+(t)`$ $`=`$ $`T^+(t)H(t).`$
In each time step the density operator $`\rho `$ changes and thus also the DFT Hamiltonian which is a functional of $`\rho `$. Similar nonlinear problem is encountered when the time-dependent Kohn-Sham equations are solved directly on a spatial grid.Puente and Serra (1999) We apply a similar procedure here discretizing time and using the Crank-Nicholson algorithm for the time integration resulting in the coupled equations

$`\rho (t+\mathrm{\Delta }t)=T(\mathrm{\Delta }t)\rho (t)T^+(\mathrm{\Delta }t)`$
$`\left\{1+{\displaystyle \frac{i\mathrm{\Delta }t}{2\mathrm{}}}H[\rho ;t+\mathrm{\Delta }t]\right\}T(\mathrm{\Delta }t)\left\{1{\displaystyle \frac{i\mathrm{\Delta }t}{2\mathrm{}}}H[\rho ;t]\right\},`$
with the initial condition, $`T(0)=1`$, already used on the right hand side of (II.2). In each time-step the equations have to be iterated until self-consistency is attained. This is performed in the truncated Fock-Darwin basis $`\{|\alpha \}`$ for the corresponding matrix version of the equations (II.2).

The accuracy and the fine-tuning of the model has been tested for the ground state with comparison to the exact results calculated for a circular parabolically confined quantum dot,Pfannkuche et al. (1993) and in the time-evolution by convincing us that the model can reproduce the Kohn results for the absorption of a Far-infrared radiation.Maksym and Chakraborty (1990); Gudmundsson and Gerhardts (1991) In addition, the normalization condition for the density matrix is tested regularly during the integration of the time-dependent problem.

### II.3 Specification of model parameters

In order to model a short GaAs quantum wire with or without an embedded impurity and several electrons we select the parameters, $`\mathrm{}\omega _0=3.37`$ meV, and assume $`g^{}=0.44`$, and $`m^{}=0.067m`$. The impurity potential is selected to be Gaussian (2) with $`V_0=10`$ meV, and $`a^2\gamma =1.0`$. The Gaussian shape is used here in order to limit the size of the basis during the CPU-intensive time integration of the problem. The wire shape of the system is attained by setting $`\alpha _2=0.7`$ and $`\alpha _i=0`$ for all $`i2`$. The charge density for 10 electrons in the system is shown in Fig. 1 with and without an embedded impurity. At $`B=0`$ T the length of the system is approximately 200nm, and the difference between the charge densities with and without an impurity in the system turns up in Fig. 1b as a fairly well localized impurity state.

This can be confirmed by the Kohn-Sham energy spectra displayed in Fig. 2, where the main change due to the addition of the impurity potential into the system is the lowering of the energy of the highest almost degenerate doubly occupied state by approximately 1 meV.

The calculations reported here are for $`B=0`$ T, and $`B=1`$ T. In order to avoid numerical problems with the determination of the eigenvectors for exactly degenerate eigenvalues in the diagonalization of the Hamiltonian we actually use $`B=10^8`$ T instead of $`B=0`$ T.

For the radiation pulse (12) of duration just over 3 ps,Planken et al. (1995) we select $`\mathrm{\Gamma }=2`$ ps<sup>-1</sup>. The envelope frequency corresponds to $`\mathrm{}\omega _1=0.658`$ meV, and the base frequency is equivalent to $`\mathrm{}\omega =2.63`$ meV. Meanwhile, we select a $`|N_p|=1`$ and $`s=0`$ in (12) to represent a dipole radiation. But, due to the finite cut-off in the basis $`|\alpha `$ the matrix elements of the excitation can never be purely of the dipole type. We choose the orientation of the initial radiation pulse (12) to be either parallel, $``$, to the long axis of the system or perpendicular to it, $``$. The finite asymmetric cut-off of the basis can introduce a tiny mixing of the perpendicular and the parallel excitation at vanishing magnetic field. In fact, the basis cut-off is performed by specifying the highest radial quantum number $`n_{\text{max}}`$ used in the calculation, by doing so an asymmetry is introduced between the number of negative and positive angular quantum numbers $`M`$. This asymmetry is made unimportant with respect to the energy spectrum and the wave functions of the system by selecting a high enough $`n_{\text{max}}`$, but the representation of a dipole potential is not quite so clean. A simple alternative method is not easy since the effective asymmetry does depend on the magnetic field.

We integrate the equations for the time-evolution operator (II.2) over a time interval of 50 ps using an increment $`\mathrm{\Delta }t=0.0025`$ ps and evaluate the expectation values for the coordinates $`x`$ and $`y`$ for each spin direction, $`x_{}`$, etc. We can thus identify oscillations perpendicular or parallel to the long axis of the system, talking about perpendicular or parallel detection. Resonances in the collective oscillations of the electron spin densities are sought by Fourier transformation of the time series in the interval 5 $``$ 50 ps, when the external excitation has been turned off and switch-on effects have vanished and the system has reached a “steady state”.

## III Results

### III.1 The $`B=0`$ T case

In Fig. 3 the Fourier power spectra of the time series for the coordinates $`x`$ and $`y`$ are used to represent the excitation spectra of the system without an impurity. The $``$ polarized excitation is seen to introduce a tiny $``$ oscillation ($`1/10`$ the amplitude) as has been explained earlier, but at the same frequency. This simple excitation spectrum does not change much with increasing excitation strength $`V_t`$. A second harmonic shows up in the parallel oscillation, but more interestingly the height of the main peak does not grow linearly with the excitation strength when $`V_ta>1`$ meV. The excitation then becomes nonlinear and as $`V_ta=5`$ meV we see that it shifts slightly to lower energy.

Before commenting on the red shift of the peak with increasing $`V_t`$ we show what happens in the system with an embedded impurity in Fig. 4. The main peak in the longitudinal oscillation (Fig. 4b) is slightly blue shifted compared to the system with no impurity. The reason is readily explained by a glance at the Kohn-Sham energy spectra in Fig. 2; the presence of the impurity potential enlarges the energy gap around the chemical potential $`\mu `$. Instead of a spatially fairly extended state just below $`\mu `$ there is now a localized state at lowered energy.

Again, with increasing excitation strength $`V_t`$ we see a red shift of the resonance in the parallel oscillations. At strong excitation we can not refer back to the Kohn-Sham spectra of the ground state, but we can expect that the energy pumped into the system makes the states just below the gap more extended, thus reducing the effective energy gap. The perpendicular oscillation in Fig. 4a is more complex than before. Again, there is a resonance at the same energy as in the parallel mode that is slightly red shifted with increased excitation, but there now appear resonances at higher energies indicating the presence of more localized modes in the system, an extra excitation mode located at the impurity site. This is again supported by a look at the Kohn-Sham energy spectra in Fig. 2a, at the changes just above the gap caused by the impurity.

The case of the initial excitation pulse aligned perpendicular with the long axis is review in Fig. 5. We see that again in the $``$ detection we have the peak structure encountered in Fig. 4a as could be expected, but the slight presence of a $``$ excitation in the initial pulse is not of the dipole kind so mostly higher order modes are excited in the $``$ direction.

Exactly this, is further displayed in Fig. 6a, where, in addition, the relative strength of the response is also visible. In contrast, the strong confinement in the perpendicular direction of the system or the wire ensures the $``$ response to be not strongly dependent on the main excitation polarization, see Fig. 6b.

To clearly demonstrate the collective oscillations in the system we show in Fig. 7 the induced electron density for the system with and without an embedded impurity for the two polarizations of the excitation pulse. The two figures Fig. 7a and b for the $``$ excitation show that the ends of the system oscillate in the $``$ direction while the oscillation around the impurity is at least at this moment parallel to the long axis. In case of no impurity in the system the oscillation of the center part is more perpendicular. When the excitation is applied $``$ to the system (see Figs. 7c and d) the main oscillation mode is also along the system of the fundamental dipole type as can be seen by comparing the induced density to the total density in Fig. 1. In this case the oscillation seems to be more extended to the whole system when there is an impurity present in it, but when there is no impurity present the center of the system displays a tiny bit of a $``$ mode superimposed on its motion. In fact, the induced density can show higher order oscillation modes that we can not discern with our time series of the expectations values for the the coordinates $`x`$ and $`y`$, which are most sensitive to dipole oscillations occuring in the system.Puente and Serra (1999); Puente et al. (2001)

The red shift of the collective resonance for the system with an embedded impurity is further explored in Fig. 8 in comparison with the response of the system without an impurity. Clearly, the resonance is red shifted and broadened with increasing excitation strength $`V_t`$ pointing towards a reduction in the binding energy of the impurity.

Comparison between the sub figures of Fig. 8 reveals that with increasing excitation strength more energy flows into a simple dipole oscillation of the system when an impurity is present in it.

This is confirmed by the induced electron density in this case showing large oscillations of the electron density around the impurity site parallel to the system, and mostly in harmony with the oscillations in the rest of the system, reflected in a single resonance. The broadening of the resonance is indicative of a finer structure that we are unable to resolve in our calculations.

The case is different with the perpendicular oscillation shown in Fig. 9, where in the end a strong excitation will lead to many modes becoming active. So, even though the parallel collective oscillation remains simple the perpendicular one is generally more complex.

### III.2 The $`B=1`$ T case

The collective oscillations in the finite magnetic field $`B=1`$ T are for short times after the initial excitation pulse very similar to the oscillations encountered in the system in no magnetic field as could be guessed by comparing the Kohn-Sham energy spectra for the two cases, see Fig. 2. It has though to be kept in mind that the Lorentz force now ensures a quick communication between the parallel and the perpendicular modes, they are not independent. But more happens: The finite size of the system together with the small Zeeman energy, and the finite temperature $`T=1`$ K establish a tiny spin polarization in the ground state, i.e. the two spin densities do not have exactly the same spatial dependence. The spatially varying exchange force in the LSDA gives a different restoring force to the spin densities. All these factors lead to a slow flow of energy into the individual oscillations of the spin densities as can be seen in Fig. 10. This flow is so slow that after 50 ps the system has still not entered a “steady state” as in the case of $`B=0`$ T.

We are still observing switch-on effects for the collective oscillations of the spin densities after 50 ps. In the spin oscillation we see present both high and low frequency components but the lack of a steady state prevents us from making a Fourier analysis. Just to avoid any misconception we should mention that the time increment of $`\mathrm{\Delta }t=0.0025`$ ps gives us 20,000 points in the range of 50 ps and at a finer level of resolution the spin oscillation looks very smooth with no apparent noise like might be hastily inferred by looking at Fig. 10.

Even more amazing is the fact that to begin with this spin oscillation increases with the strength of the excitation $`V_t`$ but, as is seen in Fig. 10, for a still stronger excitation it again decreases. Stronger excitation of the system prevents the energy in flowing into the collective spin oscillations, but directs it to the center-of-mass motion and the other collective charge oscillations. This behavior is definitely a nonlinear phenomena that can not be observed in a model of linear response.

The presence of the impurity also quite effectively directs more energy into the spin oscillation as is seen in Fig. 10. This can be explained by a stronger spatial variation of the exchange force around the impurity.

The regularity of the collective spin oscillations can be seen in Fig. 11 where simply both spin densities are seen almost overlapping at $`t=50`$ ps. The peak structure of the spin densities closely follows the structure of the ground state charge density displayed in Fig. 1.

## IV Discussion and summary

Without the external magnetic field present the collective charge oscillations excited in the system by the short THz pulse do reach a “steady state” condition within 1 or 2 picoseconds after the excitation is switched off allowing a Fourier analyzes to discern the modes present. At low excitation level we observe a linear regime, but at higher excitation nonlinear effects show up. The excitation gap in the system is reduced indicating a less localized single-electron states and in particular the effective binding energy of the impurity is reduced.

At finite magnetic field ($`B=1`$ T) a slowly increasing collective oscillation of the spin densities has not reached a “steady state” condition after 50 ps. As the spin oscillation is directly coupled with the collective charge oscillations a Fourier analyses is not possible in this case for either component. The excitation of the collective spin oscillations is much more effective when an impurity is present in the system, since then the initial spin densities at $`t=0`$ have a larger spatial variation. The more amazing fact is that with increased intensity of the initial excitation the onset of the spin oscillation can be effectively shut down. We contribute this to the flow of the energy nonadiabatially pumped into the system being more easily directed to the oscillations of the charge density in this case.

For the results shown in this paper we have always used the same type of excitation pulse. What would happen if we use a longer pulse with several oscillations or of a different form. We can not claim that we have thoroughly investigated this question, but our limited experimenting with this seems to indicate that the pulse form is indeed important for the heights of the resonance peaks, but does not affect their location strongly. We are pumping energy into the system by the excitation in varying amounts depending on the peak shape and duration and it can be expected that this also affects the initial distribution of the energy between the modes. By removing the system further from equilibrium than we do we are sure that stronger nonlinear effects can occur with considerable change of peak locations also. We already see the precursors of this in the redshift of the peaks connected to the impurity embedded in the system.

We are aware that our choice of an impurity potential represented by a Gaussian well does not correspond to unscreened Coulomb impurities that might be situated in the plane of the 2DEG in a nanosystems. This form has only been selected in order to keep the truncation of the Fock-Darwin basis within reasonable bounds in order to enable the “long-time” integration of the time evolution operator to extend to 50 ps with 20,000 steps, within each one the self-consistency was attained with further 4 loops to properly update the density operator. With these requirements in mind it should be reiterated that the present model was not designed to describe any known experimental device, but rather it was constructed to give qualitative insight into time-dependent phenomena that can take place in a nanosystem that is irradiated by short and intensive THz pulse, in response to the preparations underway in several experimental groups to study these phenomena. We only study the behavior of the system for a short time after the application of the excitation and neglect all possible damping and dissipation effects that would be caused by the coupling of the system to its immediate environment. These are certainly very important effects that await further investigation.

###### Acknowledgements.

The research was partly funded by the Icelandic Natural Science Foundation, the University of Iceland Research Fund, the National Science Council of Taiwan under Grant No. 91-2119-M-007-004, and the National Center for Theoretical Sciences, Tsing Hua University, Hsinchu Taiwan. VG acknowledges instructive discussions with Llorens Serra and a beneficial exchange of programming ideas with Ingibjörg Magnúsdóttir and Gabriel Vasile.

## The generalized one-electron overlap integral

In order to reduce the memory size needed to store the generalized overlap integral it is convenient to split it into two parts by applying the completeness of the basis once

$$\mathrm{\Phi }_{\alpha pq\beta }=𝑑𝐫\varphi _\alpha ^{}(𝐫)\varphi _p^{}(𝐫)\varphi _q(𝐫)\varphi _\beta (𝐫)=\underset{\delta }{}\mathrm{\Phi }_{\alpha p\delta }\mathrm{\Phi }_{q\beta \delta },$$
(16)
where an analytic expression for the wave functions corresponding to the Fock-Darwin basis states $`|\alpha =|M_\alpha ,n_\alpha `$ gives

$`\mathrm{\Phi }_{\alpha p\delta }=2\pi a^2\delta _{M_\alpha +M_p,M_\delta }\beta _\alpha \beta _p\beta _\delta \mathrm{\hspace{0.33em}2}^{(|M_\alpha |+|M_\beta |+|M_\delta |)/2}`$
$`{\displaystyle \underset{i_\alpha =0}{\overset{n_\alpha }{}}}{\displaystyle \underset{i_p=0}{\overset{n_p}{}}}{\displaystyle \underset{i_\delta =0}{\overset{n_\delta }{}}}{\displaystyle \frac{(1)^{i_\alpha +i_p+i_\delta }}{i_\alpha !i_p!i_\delta !}}\left({\displaystyle \genfrac{}{}{0pt}{}{n_\alpha +|M_\alpha |}{n_\alpha i_\alpha }}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{n_p+|M_p|}{n_pi_p}}\right)`$
$`\left({\displaystyle \genfrac{}{}{0pt}{}{n_\delta +|M_\delta |}{n_\delta i_\delta }}\right)\left({\displaystyle \frac{|M_\alpha |+|M_\beta |+|M_\delta |}{2}}+i_\alpha +i_p+i_\delta \right)!`$
$`\times \left({\displaystyle \frac{2}{3}}\right)^{[(|M_\alpha |+|M_\beta |+|M_\delta |)/2]+i_\alpha +i_p+i_\delta +1},`$
with

$$\beta _\alpha =\frac{1}{2^{(|M_\alpha |+1)/2}}\sqrt{\frac{n_\alpha !}{\pi (|M_\alpha |+n_\alpha )}}.$$
(17)