A nonlinear Bloch model for Coulomb interaction in quantum dots
AA nonlinear Bloch model for Coulomb interactionin quantum dots
Brigitte Bidegaray-Fesquet ∗ and Kole Keita Univ. Grenoble Alpes, LJK, BP 53, 38041 Grenoble Cedex, FranceCNRS, LJK, BP 53, 38041 Grenoble Cedex, France
[email protected] , [email protected]
November 19, 2018
Abstract
In this paper, we first derive a Coulomb Hamiltonian for electron–electron interaction in quantum dots in the Heisenberg picture. Then weuse this Hamiltonian to enhance a Bloch model, which happens to benonlinear in the density matrix. The coupling with Maxwell equations incase of interaction with an electromagnetic field is also considered fromthe Cauchy problem point of view. The study is completed by numericalresults and a discussion about the advisability of neglecting intra-bandcoherences, as is done in part of the literature.
Keywords: Maxwell–Bloch model, quantum dot, Coulomb interaction, Cauchyproblem, Liouville model, positiveness properties.
Bloch model is a very common model to describe the time evolution of a systemof electrons in different contexts such as gases of electrons, glasses, or crys-tals. The very classical case of gases and glasses involves isotropic media. Theelectrons are supposed to be localized and non interacting. Their behavior isaveraged at the mesoscopic scale. This leads to relatively simple models wherematter energy levels are quantized and labelled by integers. The case of crystals[BBFB+04] also involves integer indexed levels, but symmetries and directionsin matter have to be taken into account.Bloch model has also been extended to the description of quantum wells[HK96, KR92, HK09], and quantum dots [GH02, BF10]. In these models, matteris described by the state of two species of particles (electrons and holes, orequivalently conduction and valence electrons). In quantum wells, energy levels ∗ Corresponding author.
[email protected] , tel: +33 4 76 51 49 94, fax.: +334 76 63 12 63 a r X i v : . [ m a t h - ph ] F e b re indexed by vectors, which correspond to displacements over the underlyinglattice. In contrast, the confinement of electrons in quantum dots leads tointeger indexed levels like in gases, which often leads to consider quantum dotsas pseudo-atoms, but this is a very raw vision. In particular, among otherdifferences, electrons at the same mesoscopic location do interact directly via Coulomb interaction.[BF10] is a preliminary paper that derives basic Bloch equations for twospecies of electrons (conduction and valence) only taking into account the freeelectron Hamiltonian and the interaction with a laser electric field. The aimof the present paper is to include properly Coulomb interaction in this model.Beyond the sole derivation of the model, we want in particular to study its math-ematical properties. In the continuation of [Bid01, BBR01], we want to showthat a certain number of properties are preserved through the time evolution,such as Hermicity and positiveness of the density matrix.
The outline of this paper is as follows. We devote Section 1.2 to the descriptionof the basic Bloch model which does not include Coulomb interaction but onlythe free energies of the electrons and the action of an electromagnetic field. InSection 2, we derive the Coulomb Hamiltonian in terms of the conduction andvalence operators. The associated Heisenberg equation is derived in Section3, but it ends up with an open system of equations. The system is closedusing the Wick theorem, and the final Bloch equation has a Liouville form,but is nonlinear. This nonlinearity does not allow to use previous literaturedirectly and we prove anew the Hermicity, positiveness and boundedness resultsin Section 4. In Section 5, we show the impact of the Coulomb contributionin numerical results and also compare our model with a vanishing intra-bandcoherence model defined in [GH02].
Let us first recall the main results obtained in [BF10] and fix the notations.
Let A and B be two operators, we define their commutator by [ A, B ] = AB − BA and their skew-commutator by { A, B } = AB + BA . For an operator A , we definethe associated observable (cid:104) A (cid:105) = Tr( S A ) by averaging with respect to the initialstate density S of the system. If the system is described by a Hamiltonian H ,the time-evolution for this observable is given by the Heisenberg equationi (cid:126) ∂ t (cid:104) A (cid:105) = (cid:104) [ A, H ] (cid:105) , (1)where (cid:126) is the reduced Planck constant. When the observable is the densitymatrix, the Heisenberg equation of motion is called the Bloch equation.2 .2.2 Operators for quantum dots A quantum dot is defined as a collection of conduction and valence electrons.There is of course no conduction in quantum dots since the electrons are confinedin every direction, but this terminology is useful to distinguish between thevalence electrons — which are in fact the absence of holes in the valence band— and the free, but confined, electrons. For each species, energy levels arequantized and indexed by a set of integers, I c and I v , for conduction and valenceelectrons, respectively. For i ∈ I c , we define the creation and annihilationoperators c † i and c i . Likewise, for valence electrons, we define the creation andannihilation operators v † i and v i . Property 1. { c i , c † j } = δ i,j , { c i , c j } = { c † i , c † j } = 0 , { v i , v † j } = δ i,j , { v i , v j } = { v † i , v † j } = 0 , where δ i,j denotes the Kronecker symbol.This implies in particular that c i c i = c † i c † i = 0, which means that it isimpossible to create twice or annihilate twice the same electron. This is thePauli exclusion rule: electrons are fermions. Of course any conduction operatorcommutes with any valence operator. The observable we are interested in is the density matrix . It includes a con-duction density matrix , which elements are the ρ c ij = (cid:104) c † j c i (cid:105) . This matrix isHermitian and positive semi-definite. Its diagonal terms ρ c ii = (cid:104) c † i c i (cid:105) are alsocalled populations and give the probability to find an electron in state i . Theoff-diagonal terms, ρ c ij , i (cid:54) = j are called (intra-band) coherences . Of course wealso define a valence density matrix , which elements are the ρ v ij = (cid:104) v † j v i (cid:105) . Be-sides we are interested in inter-band coherences defined by ρ cv ij = (cid:104) v † j c i (cid:105) . Theentries of these matrices are the variables of the Bloch equation. They are castin a single density matrix ρ = (cid:18) ρ c ρ cv ρ vc ρ v (cid:19) where ρ vc = ρ cv ∗ , which ensures that ρ is Hermitian and positive semi-definite. The basic Bloch equation for quantum dot is derived in [BF10]. It takes intoaccount two types of Hamiltonians in the Heisenberg equation for the densitymatrix, namely free electron Hamiltonians and interaction Hamiltonians with alaser field. The free electron Hamiltonians read H c0 = (cid:88) k ∈ I c (cid:15) c k c † k c k , H v0 = (cid:88) k ∈ I v (cid:15) v k v † k v k , (cid:15) c k and (cid:15) v k arethe free electron energies associated with each electron level. The electron levelsare described by the wave functions ψ c k and ψ v k , solutions to the free electronSchr¨odinger equation subject to the boundary conditions of the quantum dot(see [HK09]). The interaction with a laser characterized by its time-dependentelectric field E ( t ) is described by the Hamiltonians H Lc = 12 (cid:88) ( k,l ) ∈ ( I c ) ( E ( t ) · M c kl c † k c l + E ∗ ( t ) · M c ∗ kl c † l c k ) ,H Lv = 12 (cid:88) ( k,l ) ∈ ( I v ) ( E ( t ) · M v kl v † k v l + E ∗ ( t ) · M v ∗ kl v † l v k ) ,H Lcv = (cid:88) ( k,l ) ∈ I c × I v ( E ( t ) · M cv kl c † k v l + E ∗ ( t ) · M cv ∗ kl v † l c k ) , where the dipolar moment matrices are matrices with entries in C and may beexpressed in terms of the wave functions associated to each level: M c kl = (cid:90) d r ψ c ∗ l ( r ) er ψ c k ( r ) , M v kl = (cid:90) d r ψ v ∗ l ( r ) er ψ v k ( r ) , M cv kl = (cid:90) d r ψ v ∗ l ( r ) er ψ c k ( r ) , where e is the unsigned charge of the electron. Injecting these Hamiltonians inthe Heisenberg equation, the basic Bloch equations can be cast in Liouville formi (cid:126) ∂ t ρ = [ V ( t ) , ρ ] . (2)In equation (2), V ( t ) = V F + V E ( E ( t )) is a sum of a constant term V F stemmingfrom the free energies collected in diagonal matrices E c0 = diag( { (cid:15) c i } i ∈ I c ) and E v0 = diag( { (cid:15) v i } i ∈ I v ), and a time dependent term due to the interaction withthe electric field: V F = (cid:18) E c0 E v0 (cid:19) and V E ( E ( t )) = (cid:18) (cid:60) E ( t ) · M c E ( t ) · M cv E ∗ ( t ) · M cv ∗ (cid:60) E ( t ) · M v (cid:19) . The scalar product of the electrical field and the dipolar moment matrix, is ascalar product in C which yields a matrix with entries in C , with the samedimension as ρ , as is necessary to give a meaning to the right-hand side of (2). Equation (2) clearly preserves the Hermitian structure of ρ . Its exact solutionis ρ ( t ) = exp (cid:18) − i (cid:126) (cid:90) t V ( τ ) d τ (cid:19) ρ (0) exp (cid:18) i (cid:126) (cid:90) t V ( τ ) d τ (cid:19) . (3)4his expression allows to prove a certain number of properties (see [BF10]).Let d = card( I c ) + card( I v ) be the total number of levels, given a positivesemi-definite initial data ρ (0) ∈ M d ( C ) and a continuous electric field E ( t ), • equation (2) is globally well-posed, i.e. there exists a unique solution ρ ∈ C ( R + ; M d ( C )) that exists for all times t ≥ • for all time ρ ( t ) is a positive semi-definite matrix; • its trace is conserved through the time evolution. Coulomb interaction can be introduced using field operators. We denote byˆ ψ † c ( r ) and ˆ ψ † v ( r ) the creation field-operators of respectively a conduction elec-tron and a valence electron at the space location r , and ˆ ψ c ( r ) and ˆ ψ v ( r ) thecorresponding annihilation field-operators. We consider that there are N rele-vant electrons in the quantum dot, and can write the Coulomb Hamiltonians as H c − c = 12 N (cid:88) i,j =1 (cid:90) (cid:90) d r i d r j ˆ ψ † c ( r i ) ˆ ψ † c ( r j ) V c ( r i , r j ) ˆ ψ c ( r j ) ˆ ψ c ( r i ) , (4a) H v − v = 12 N (cid:88) i,j =1 (cid:90) (cid:90) d r i d r j ˆ ψ v ( r i ) ˆ ψ v ( r j ) V v ( r i , r j ) ˆ ψ † v ( r j ) ˆ ψ † v ( r i ) , (4b) H c − v = N (cid:88) i,j =1 (cid:90) (cid:90) d r i d r j ˆ ψ † c ( r i ) ˆ ψ v ( r j ) V c − v ( r i , r j ) ˆ ψ † v ( r j ) ˆ ψ c ( r i ) , (4c)where V c , V v and V c − v are the conduction–conduction, valence–valence andconduction–valence Coulomb potentials. The Coulomb potentials have the form V c ( r, r (cid:48) ) = V v ( r, r (cid:48) ) = k C | r − r (cid:48) | and V c − v ( r, r (cid:48) ) = − k C | r − r (cid:48) | , where k C is Coulomb’s constant. The difference of treatment of conduction andvalence electrons stems from the fact that Coulomb interaction describes theinteraction of electrons and holes (see e.g. [HK09]) and that the presence of anelectron in the valence band is indeed the absence of the corresponding hole,and vice-versa, which inverts the role of creation and annihilation operators.The total Coulomb Hamiltonian is H C = H c − c + H v − v + H c − v .We want to derive Bloch-type equations including the Coulomb interaction.Bloch equations have the advantage not to depend explicitly on the exact formof the field-operators. To this aim, we writeˆ ψ c ( r ) = (cid:88) α ∈ I c ψ c α ( r ) c α , ˆ ψ † c ( r ) = (cid:88) α ∈ I c ψ c ∗ α ( r ) c † α , (5a)5 ψ v ( r ) = (cid:88) α ∈ I v ψ v α ( r ) v α , ˆ ψ † v ( r ) = (cid:88) α ∈ I v ψ v ∗ α ( r ) v † α . (5b)They are weighted by the conduction and valence electron wave functions ψ c α ( r )and ψ v α ( r ), which are the same as those who occurred in the expression of thedipolar moment matrices. In the sequel, to avoid unnecessary written complex-ity, we will often omit to specify which set the indices belong to.Inserting decompositions (5) in Hamiltonians (4) we obtain H c − c = (cid:88) α ,α ,α (cid:48) ,α (cid:48) R c α α α (cid:48) α (cid:48) c † α c † α c α (cid:48) c α (cid:48) , (6a) H v − v = (cid:88) α ,α ,α (cid:48) ,α (cid:48) R v α α α (cid:48) α (cid:48) v α (cid:48) v α (cid:48) v † α v † α , (6b) H c − v = − (cid:88) α ,α ,α (cid:48) ,α (cid:48) R c − v α α α (cid:48) α (cid:48) c † α v α (cid:48) v † α c α (cid:48) , (6c)where R c α α α (cid:48) α (cid:48) = N (cid:90) (cid:90) d r d r (cid:48) ψ c ∗ α ( r ) ψ c ∗ α ( r (cid:48) ) V c ( r, r (cid:48) ) ψ c α (cid:48) ( r (cid:48) ) ψ c α (cid:48) ( r ) , (7a) R v α (cid:48) α (cid:48) α α = N (cid:90) (cid:90) d r d r (cid:48) ψ v α (cid:48) ( r ) ψ v α (cid:48) ( r (cid:48) ) V v ( r, r (cid:48) ) ψ v ∗ α ( r (cid:48) ) ψ v ∗ α ( r ) , (7b) R c − v α α α (cid:48) α (cid:48) = − N (cid:90) (cid:90) d r d r (cid:48) ψ c ∗ α ( r ) ψ v α (cid:48) ( r (cid:48) ) V c − v ( r, r (cid:48) ) ψ v ∗ α ( r (cid:48) ) ψ c α (cid:48) ( r ) . (7c)The symmetries in the integrands of (7) induce the following properties. Property 2.
Since V ( r, r (cid:48) ) is an even function of r − r (cid:48) , variables r and r (cid:48) playthe same role and R c α α α (cid:48) α (cid:48) = R c α α α (cid:48) α (cid:48) , R v α α α (cid:48) α (cid:48) = R v α α α (cid:48) α (cid:48) . Since V ( r, r (cid:48) ) is a real valued function R c α α α (cid:48) α (cid:48) = R c ∗ α (cid:48) α (cid:48) α α , R v α α α (cid:48) α (cid:48) = R v ∗ α (cid:48) α (cid:48) α α ,R c − v α α α (cid:48) α (cid:48) = (cid:16) R c − v α (cid:48) α (cid:48) α α (cid:17) ∗ . The Pauli exclusion principle (skew-commutation, Property 1) also inducesthat some terms in H c − c and H v − v are necessarily zero. Property 3. If α = α or α (cid:48) = α (cid:48) , c † α c † α c α (cid:48) c α (cid:48) = 0 , v α (cid:48) v α (cid:48) v † α v † α = 0 . r i and r j .Then they are annihilated while interacting and recreated at the same locations. Definition 1.
A product of operators will be said to be in the normal order ,if the annihilation operators are on the right and the creation operators on theleft. H c − c already follows a normal ordered form and we can keep it untouched: H c − c = (cid:88) α ,α ,α (cid:48) ,α (cid:48) R c α α α (cid:48) α (cid:48) c † α c † α c α (cid:48) c α (cid:48) . (8a)Hamiltonians H v − v and H c − v given by equations (6b) and (6c) do not follownormal ordered forms.To express H v − v we need to compute a normal ordered form of v α (cid:48) v α (cid:48) v † α v † α : v α (cid:48) v α (cid:48) v † α v † α = δ α α (cid:48) δ α α (cid:48) − δ α α (cid:48) v † α v α (cid:48) − δ α (cid:48) α δ α (cid:48) α + δ α (cid:48) α v † α v α (cid:48) + δ α (cid:48) α v † α v α (cid:48) − δ α α (cid:48) v † α v α (cid:48) + v † α v † α v α (cid:48) v α (cid:48) . Thanks to Property 2, R v α α α (cid:48) α (cid:48) = R v α α α (cid:48) α (cid:48) , therefore − δ α α (cid:48) v † α v α (cid:48) and − δ α α (cid:48) v † α v α (cid:48) lead to the same contribution. The same argument can be appliedto δ α (cid:48) α v † α v α (cid:48) and δ α (cid:48) α v † α v α (cid:48) . Hence H v − v = 2 (cid:88) α,α (cid:48) ,β ( R v βαα (cid:48) β − R v βαβα (cid:48) ) v † α v α (cid:48) + (cid:88) α ,α ,α (cid:48) ,α (cid:48) R v α α α (cid:48) α (cid:48) v † α (cid:48) v † α (cid:48) v α v α . (8b)In the definition of H v − v we have dropped the δ α α (cid:48) δ α α (cid:48) and − δ α (cid:48) α δ α (cid:48) α terms which would lead to zero contributions in the Heisenberg equation.In the same way c † α v α (cid:48) v † α c α (cid:48) = δ α α (cid:48) c † α c α (cid:48) − c † α v † α v α (cid:48) c α (cid:48) , hence H c − v = − (cid:88) α,α (cid:48) ,β R c − v αβα (cid:48) β c † α c α (cid:48) + (cid:88) α ,α ,α (cid:48) ,α (cid:48) R c − v α α α (cid:48) α (cid:48) c † α v † α v α (cid:48) c α (cid:48) . (8c) We now write Heisenberg equation (1) where A are operators c † j c i , v † j v i or v † j c i and H are the Hamiltonians defined by Equation (8). For the free Hamiltonians and the laser interactions, the computation of thecommutators led to expressions in terms of the two-operator densities, and7herefore to a closed set of equations [BF10]. This will not be the case anymore here, since the commutators stemming from Coulomb Hamiltonians willgive rise to four-operator densities. To go further we should a priori have evo-lution equations for these observables via the Heisenberg equation, computingcommutators with the already defined Coulomb Hamiltonians. This would leadinevitably to six-operator densities, and so on. To avoid this endless procedure,we have to close the system at some point. This is the goal of the Wick theorem[Wic50] which amounts in our case to write the four-operator densities as sumsof products of two-operator densities following e.g. the rule (cid:10) c † α c † α c α (cid:48) c α (cid:48) (cid:11) (WT) = (cid:10) c † α c α (cid:48) (cid:11) (cid:10) c † α c α (cid:48) (cid:11) − (cid:10) c † α c α (cid:48) (cid:11) (cid:10) c † α c α (cid:48) (cid:11) , where the symbol (WT) = means ”is approximated through Wick theorem”. c † j c i In order to derive the Heisenberg equation, we have to compute [ c † j c i , H c − c ] and[ c † j c i , H c − v ] (since [ c † j c i , H v − v ] is clearly zero). H c − c According to Equation (8a)[ c † j c i , H c − c ] = (cid:88) α ,α ,α (cid:48) ,α (cid:48) R c α α α (cid:48) α (cid:48) [ c † j c i , c † α c † α c α (cid:48) c α (cid:48) ] . Remark . We already know many situations where [ c † j c i , c † α c † α c α (cid:48) c α (cid:48) ] is nec-essarily zero: • if none of the indices α , α , α (cid:48) , α (cid:48) is equal either to i or j , • if α = α or α (cid:48) = α (cid:48) (see Property 3).We compute separately each commutator:[ c † j c i , c † α c † α c α (cid:48) c α (cid:48) ] = δ iα c † j c † α c α (cid:48) c α (cid:48) + δ iα c † j c † α c α (cid:48) c α (cid:48) − δ jα (cid:48) c † α c † α c α (cid:48) c i − δ jα (cid:48) c † α c † α c α (cid:48) c i . Using Property 2, we see that[ c † j c i , H c − c ] = (cid:88) α ,α ,α (cid:48) ,α (cid:48) R c α α α (cid:48) α (cid:48) (cid:16) δ iα c † j c † α c α (cid:48) c α (cid:48) − δ jα (cid:48) c † α c † α c α (cid:48) c i (cid:17) = 2 (cid:88) α,α (cid:48) ,α (cid:48) R c iαα (cid:48) α (cid:48) c † j c † α c α (cid:48) c α (cid:48) − (cid:88) α ,α ,α (cid:48) R c α α jα (cid:48) c † α c † α c α (cid:48) c i .
8e now apply Wick theorem which leads to (cid:68) [ c † j c i , H c − c ] (cid:69) (WT) = 2 (cid:88) α,α (cid:48) ,α (cid:48) R c iαα (cid:48) α (cid:48) (cid:16)(cid:68) c † j c α (cid:48) (cid:69) (cid:10) c † α c α (cid:48) (cid:11) − (cid:68) c † j c α (cid:48) (cid:69) (cid:10) c † α c α (cid:48) (cid:11)(cid:17) − (cid:88) α ,α ,α (cid:48) R c α α jα (cid:48) (cid:0)(cid:10) c † α c i (cid:11) (cid:10) c † α c α (cid:48) (cid:11) − (cid:10) c † α c α (cid:48) (cid:11) (cid:10) c † α c i (cid:11)(cid:1) = 2 (cid:88) k,α,α (cid:48) ( R c iαkα (cid:48) − R c iαα (cid:48) k ) (cid:10) c † α c α (cid:48) (cid:11) (cid:68) c † j c k (cid:69) − (cid:88) k,α,α (cid:48) (cid:0) R c kαjα (cid:48) − R c αkjα (cid:48) (cid:1) (cid:10) c † α c α (cid:48) (cid:11) (cid:68) c † k c i (cid:69) . Defining matrix Λ c ( ρ ) asΛ c ik ( ρ ) = 2 (cid:88) ( α,α (cid:48) ) ∈ ( I c ) ( R c iαkα (cid:48) − R c iαα (cid:48) k ) ρ c α (cid:48) α , (9)we can cast the result as (cid:68) [ c † j c i , H c − c ] (cid:69) (WT) = [Λ c ( ρ ) , ρ c ] ij . H c − v The same sort of computation as in Section 3.2.1 is performed on Equation (8c)to compute [ c † j c i , H c − v ].Using the matrices ζ v ( ρ ), γ c − v ( ρ ), and η c − v ik where ζ v ik ( ρ ) = (cid:88) ( α,α (cid:48) ) ∈ ( I v ) R c − v iαkα (cid:48) ρ v α (cid:48) α , (10) γ c − v ik ( ρ ) = − (cid:88) ( α,α (cid:48) ) ∈ I v × I c R c − v iαα (cid:48) k ρ cv α (cid:48) α , (11) η c − v ik = − (cid:88) β ∈ I v R c − v iβkβ , (12)we obtain that (cid:68) [ c † j c i , H c − v ] (cid:69) (WT) = [ ζ v ( ρ ) + η c − v , ρ c ] ij + (cid:88) k γ c − v ik ( ρ ) ρ vc kj − (cid:88) k ρ cv ik γ c − v ∗ kj ( ρ ) . v † j v i H v − v According to Equation (8b), we have to evaluate two types of commutators tocompute [ v † j v i , H v − v ]. The first commutator is clearly computed in the same9ay as [ c † j c i , H c − c ] replacing conduction electron operators by valence ones. Forthe second part, we have to compute [ v † j v i , v † α v α (cid:48) ]. We obtain (cid:68) [ v † j v i , H v − v ] (cid:69) (WT) = [Λ v ( ρ ) + κ v , ρ v ] ij , where Λ v ( ρ ) and κ v are defined asΛ v ik ( ρ ) = 2 (cid:88) ( α,α (cid:48) ) ∈ ( I v ) ( R v iαkα (cid:48) − R v iαα (cid:48) k ) ρ v α (cid:48) α , (13) κ v ik = 2 (cid:88) β ∈ I v ( R v βikβ − R v βiβk ) . (14) H c − v According to Equation (8c), the commutator [ v † j v i , H c − v ] only involves the firstterm of H c − v , which we can write as[ v † j v i , H c − v ] = (cid:88) α ,α ,α (cid:48) ,α (cid:48) R c − v α α α (cid:48) α (cid:48) c † α c α (cid:48) [ v † j v i , v † α v α (cid:48) ] . The simplifaction of the commutators and the Wick theorem leads to (cid:68) [ v † j v i , H c − v ] (cid:69) (WT) = [ ζ c ( ρ ) , ρ v ] ij + (cid:88) k γ c − v ∗ ik ( ρ ) ρ cv kj − (cid:88) k ρ vc ik γ c − v kj ( ρ ) , using the previously defined γ c − v ( ρ ) and the new matrix ζ c ( ρ ) defined by ζ c ik ( ρ ) = (cid:88) ( α,α (cid:48) ) ∈ ( I c ) R c − v αiα (cid:48) k ρ c α (cid:48) α . (15) v † j c i The commutators involving v † j c i are all possible to write with already definedmatrices, indeed we first compute [ v † j c i , H c − c ] according to Equation (8a) andrecognize (cid:68) [ v † j c i , H c − c ] (cid:69) (WT) = (cid:88) k Λ c ik ( ρ ) ρ cv kj . Then Equation (8b) leads to (cid:68) [ v † j c i , H v − v ] (cid:69) (WT) = − (cid:88) k ρ cv ik Λ v kj ( ρ ) − (cid:88) k ρ cv ik κ v kj , and Equation (8c) to (cid:68) [ v † j c i , H c − v ] (cid:69) (WT) = (cid:88) k ζ v ik ( ρ ) ρ cv kj + (cid:88) k γ c − v ik ( ρ ) ρ v kj − (cid:88) k ρ c ik γ c − v kj ( ρ ) − (cid:88) k ρ cv ik ζ c kj ( ρ ) + (cid:88) k η c − v ik ρ cv kj . .5 Matrix formulation We would like to cast the former results as (cid:42)(cid:34)(cid:32) c † j c i v † j c i c † j v i v † j v i (cid:33) , H C (cid:35)(cid:43) (WT) = [ V C ( ρ ) , ρ ] , where V C ( ρ ) = (cid:18) V c ( ρ ) V c − v ( ρ ) V v − c ( ρ ) V v ( ρ ) (cid:19) , which implies (cid:68) [ c † j c i , H C ] (cid:69) (WT) = V c ( ρ ) ρ c + V c − v ( ρ ) ρ vc − ρ c V c ( ρ ) − ρ cv V v − c ( ρ ) , (cid:68) [ v † j v i , H C ] (cid:69) (WT) = V v − c ( ρ ) ρ cv + V v ( ρ ) ρ v − ρ vc V c − v ( ρ ) − ρ v V v ( ρ ) , (cid:68) [ v † j c i , H C ] (cid:69) (WT) = V c ( ρ ) ρ cv + V c − v ( ρ ) ρ v − ρ c V c − v ( ρ ) − ρ cv V v ( ρ ) . Identifying the coefficients computed in Sections 3.2, 3.3, and 3.4, we obtainthat the result can indeed be cast as [ V C ( ρ ) , ρ ] where V c ( ρ ) = Λ c ( ρ ) + ζ v ( ρ ) + η c − v ,V c − v ( ρ ) = γ c − v ( ρ ) ,V v − c ( ρ ) = γ c − v ∗ ( ρ ) = V c − v ∗ ( ρ ) ,V v ( ρ ) = Λ v ( ρ ) + ζ c ( ρ ) + κ v , and the various matrices have been defined by Equations (9) to (15).The evolution equation for the density matrix including also the free electronHamiltonians, the interaction with a laser field, and Coulomb interaction cantherefore be cast in a Liouville formi (cid:126) ∂ t ρ = [ V ( t, ρ ( t )) , ρ ] , (16)where V ( t, ρ ( t )) = V ( t ) + V C ( ρ ( t )) and V ( t ) has been introduced in Equation(2). Recall (see Equation (2)) the basic Bloch equationi (cid:126) ∂ t ρ = [ V ( t ) , ρ ] , where V ( t ) = (cid:18) E c0 + (cid:60) E ( t ) · M c E ( t ) · M cv E ∗ ( t ) · M cv ∗ E v0 + (cid:60) E ( t ) · M v (cid:19) . The symmetries in the definition of M c and M v imply that their diagonal entriesare zero. On the contrary for E c0 and E v0 there are only diagonal entries. In the11volution equation of ρ c ij , the ρ c ij term only involves E c0 , and the other entriesof ρ play a role through M c and M cv .We can therefore easily analyze the Coulomb contributions in terms of shiftson the free electron energies and off-diagonal terms. Hence we compute δ(cid:15) c i ( ρ ) = 2 (cid:88) ( α,α (cid:48) ) ∈ ( I c ) ( R c iαiα (cid:48) − R c iαα (cid:48) i ) ρ c α (cid:48) α + (cid:88) ( α,α (cid:48) ) ∈ ( I v ) R c − v iαiα (cid:48) ρ v α (cid:48) α − (cid:88) β ∈ I v R c − v iβiβ ,δ(cid:15) v i ( ρ ) = 2 (cid:88) ( α,α (cid:48) ) ∈ ( I v ) ( R v iαiα (cid:48) − R v iαα (cid:48) i ) ρ v α (cid:48) α + (cid:88) ( α,α (cid:48) ) ∈ ( I c ) R c − v αiα (cid:48) i ρ c α (cid:48) α + 2 (cid:88) β ∈ I v ( R v βiiβ − R v βiβi ) . We can define the energy shift matrices as δE c ( ρ ) = diag( { δ(cid:15) c i ( ρ ) } i ∈ I c ) and δE v ( ρ ) = diag( { δ(cid:15) v i ( ρ ) } i ∈ I v ). Hence V ( t, ρ ( t )) = E ( ρ ( t )) + R ( t, ρ ( t )), where E ( t, ρ ( t )) = diag( E c ( t, ρ ( t )) , E v ( t, ρ ( t ))) and E c ( ρ ( t )) = E c0 + δE c ( ρ ( t )) and E v ( ρ ( t )) = E v0 + δE v ( ρ ( t )) . The model for quantum dots given in [GH02] has clearly been derived mimickingquantum well models as in [KR92]. In quantum wells, electrons and holes caninteract only if they ”see each other” long enough, i.e. if they have (and areindexed by) the same wave vector. This leads morally to weakly coupled two-level systems. Therefore, in [GH02], the variables are the level populations andthe inter-band coherences. In our model, this means that intra-band coherences ρ c ij and ρ v ij for i (cid:54) = j are not considered.We can therefore wonder what becomes of our model if we set intra-bandcoherences to zero. First we notice that in the evolution equation of, e.g., ρ c ij , there is, e.g., a contribution of V c ij ( ρ c jj − ρ c ii ). Hence, even if intra-bandcoherences are initially zero, they are not always zero through the evolutionwith Equation (16). Setting artificially intra-band coherences to zero thereforedestroys the Liouville structure of the system.In our model, this hypothesis also changes the definition of Λ c ( ρ ), ζ v ( ρ ) andΛ v ( ρ ), which becomeΛ c ik ( ρ ) = 2 (cid:88) α ∈ I c ( R c iαkα − R c iααk ) ρ c αα ,ζ v ik ( ρ ) = (cid:88) α ∈ I v R c − v iαkα ρ v αα , Λ v ik ( ρ ) = 2 (cid:88) α ∈ I v ( R v iαkα − R v iααk ) ρ v αα . In the same way as in the usual Bloch case [BBR01], the Liouville structure al-lows to state that the density matrix remains Hermitian through time evolution.12ts trace is conserved and it remains a positive semi-definite matrix. Hence fora given electric field E ( t ), the elements of the density matrix are bounded, moreprecisely populations are bounded | ρ c ii ( t ) | , | ρ v jj ( t ) | ≤ Tr( ρ ( t )) = Tr( ρ (0)) , for all i ∈ I c and j ∈ I v , as well as coherences | ρ c ij ( t ) | ≤ (cid:113) ρ c ii ( t ) ρ c jj ( t ) ≤ Tr( ρ (0))2 , for all i, j ∈ I c , and likewise | ρ v jk ( t ) | , | ρ c − v ij ( t ) | ≤ Tr( ρ (0))2 , for all i ∈ I c and j, k ∈ I v . The density matrix governed by Bloch equations in Section 3 was only dependingon time. We can now consider a collection of quantum dots which are scatteredin space and interacting, not directly but through the interaction with an elec-tromagnetic wave that propagates through the medium. This can be modeledby a density matrix, that now depends on time and space, which is coupled withMaxwell equations for the laser field through the expression of polarization. Wetherefore address the system µ∂ t H = − curl E ,ε∂ t E = curl H − ∂ t P , P = N b Tr( M ρ ) ,∂ t ρ = − i (cid:126) [ V ( ρ ) , ρ ] , (17)where all the variables depend on time and space in 3 dimensions: the electricand magnetic fields E and H , the polarization P , and the density matrix ρ . InEquation (17), ε and µ denote the electromagnetic permittivity and permeabilityof the underlying medium. They both can depend on the space variable. Thedensity of quantum boxes is given by N b . The Bloch and Maxwell equationsare coupled via the polarization that involve the dipolar moment matrix M .Such models have already been written and studied mathematically andnumerically in a few physical contexts. Here the specificity is the fact that theBloch equation is nonlinear in ρ . Note that even in the case when V does notdepend on ρ , the full coupled model is already nonlinear since V is affine in E .This system can be cast in the abstract setting of [DS12]. In this paper,they introduce a general abstract setting able to treat both Maxwell–Landau–Lifschitz and classical Maxwell–Bloch equations. In this setting, the electromag-netic field is supposed to exist in all space R . Matter described by the densitymatrix is only occupying a subdomain Ω of R . The variables are gathered inone variable U = ( u, v ), where u = ( u , u ) = ( H , E ) and v = ρ . The variable13 can be viewed as 6 real variables and variable v as d real variables. Thisvariable U is supposed to be in L = L ( R ; R ) × L (Ω; R d ). The abstractsystem reads (cid:40) ( ∂ t + B ) u = ( κ − · l ) F (¯ v, u ) , for x ∈ R ,∂ t v = F ( v, u ) , for x ∈ Ω . (18)In this formulation κ ( x ) = ( κ ( x ) , κ ( x )) = ( µ, ε ). It is uniformly positiveas needed in [DS12]. Let H curl be the space of functions f ∈ L ( R ; R ) withcurl f ∈ L ( R ; R ). The linear differential operator B is defined on H curl × H curl by B ( u , u ) = ( κ − curl u , − κ − curl u ). The variable v is extended by ¯ v onthe whole R and is zero outside Ω. We can identify l = 0, l = −N b Tr( M · ),and F ( v, u ) = − i (cid:126) [ V ( ρ ) , ρ ]. The system (17) verifies the hypotheses given in[DS12], namely, • F is affine in u : F ( v, u ) = F ( v ) + F ( v ) u , • for j = 0 , F j (0) = 0, • for all R > C F ( R ) such that for all v ∈ B R (ball of radius R in R d ), | F j ( v ) | + | ∂ v F j ( v ) | ≤ C F ( R ), • there exists K ≥ u, v ) ∈ R × R d , F ( v, u ) · v ≤ K | u | .We have in particular used the L ∞ bounds of Section 4.1, more precisely, welook for v ∈ L ∞ ((0 , ∞ ); L ∞ (Ω; R d )) . (19)Besides we suppose to have at time t = 0 the conditionsdiv( κ j u j − l j ¯ v ) = 0 , for j = 1 , , (20)which are indeed the physical conditions div( µ H ) = 0 and div( ε E + P ) = 0[Dum05]. The structure of Equation (18) ensures that this condition holds forall time if it is valid at the initial time. In this section, we state without proof the results obtained in [DS12] and thatwe can apply to our context. The first result addresses the existence of globalfinite energy solutions.
Theorem 1 (Theorem 3, [DS12]) . For any initial data U = ( u , v ) ∈ L ( R ; R ) × ( L (Ω; R d ) ∩ L ∞ (Ω; R d )) satisfying (20) , there exists U ∈ C ([0 , ∞ ); L ) whichis a solution to (18) – (20) , and satisfies the finite energy condition (19) . More-over, for all T > , there exists a constant C that only depends on T , F , l and (cid:107) v (cid:107) L ∞ , such that for all t ∈ [0 , T ] , (cid:107) U ( t ) (cid:107) L ≤ C (cid:107) U (cid:107) L .
14o have a uniqueness result, we need some regularity on ε and µ . Usuallyin physical contexts, ε and µ may be discontinuous across the boundary of Ω,but we do not know how to tackle with this problem. We hence assume thatΩ is bounded and κ i − ∈ C ∞ ¯Ω ( R ) , for i = 1 , , (21)the space of C ∞ functions on R with compact support included in ¯Ω, whichmeans in particular that κ i is 1 outside Ω and the transition across the boundaryis smooth. We also assume that the initial data for the electromagnetic wave issmooth enough. Theorem 2 (Theorem 5, [DS12]) . Under the assumptions of Theorem 1 and (21) , and assuming that curl u i ∈ L ( R ) for i = 1 , , there exists only onesolution to (18) – (20) with initial data U , given by Theorem 1. It satisfies curl u i ∈ C ([0 , ∞ ); L ( R )) for i = 1 , . This stems from the fact that l = 0 and F does not depend on u but only u . Self-Induced Transparency (SIT) is a typical two-level phenomenon: using alight pulse which is resonant with the transition, absorption and stimulatedemission are combined to obtain exact population inversion and an unchangedelectric field. This phenomenon has been predicted theoretically and confirmedexperimentally [AE87, MH67, GS70].The propagating field is a pulse given by E ( t, z ) = E ( t, z ) sin( ω ( t − z/v )) , where v is the velocity of the pulse, ω is both the center frequency of the pulseand the transition frequency of the medium, and E ( t, z ) is the pulse envelope.It is shown that the envelope is not reshaped by the medium, only if it is asymmetric hyperbolic secant E ( t, z ) = E sech (cid:18) t − z/vτ (cid:19) , where E = 2 mτ . In this expression, τ is the pulse duration and m = M/ (cid:126) , where M is the dipolarmoment associated to the transition. According to the Area Theorem [MH67],the medium undergoes k exact inversions if A = m (cid:90) ∞−∞ E ( t, z ) dt = kπ. The corresponding pulse is called a kπ -pulse. It is easy to compute that A = mτ E [arctan(sinh t )] ∞−∞ = mτ E π, kπ -pulse is obtained for the amplitude E = k (cid:126) /M τ . For ourtest-cases, we will use a 2 π -pulse, for which the medium is inverted and goeseventually to its original state. The return to the initial state is an easy-to-check criterion to validate numerical approaches, as has been already done in[BF06, ZAG95]. In this paper we want in particular to investigate the validity of the vanishingintra-band coherence assumption. To this aim we need a minimum of three levelsand we therefore adapt the SIT experiment to our framework. We absolutely donot claim that SIT has any practical application in the quantum dot context,but only choose this test-case because of the easiness to interpret the results. ω c1 ω v1 ω ω c1 ω v2 ω ω v1 ω ω v1 ω (a) (b)Figure 1: Adaption of the SIT test case to the quantum dot context. (a) Originaltwo-level case; (b) 2 three-level test cases.In Figure 1(a), the original two-level test case is represented, for which thereis a single conduction level and a single valence level, separated by the energycorresponding to the field frequency. The upper plot represents the (normalized)time-evolution of the electric field. The time-evolution of the population of theinitially empty conduction level is given by the lower curve. We observe thatthe medium undergoes two complete population inversions.In Figure 1(b), we have two three-level test cases with a single conductionlevel and two valence levels. In the first place we do not take Coulomb interac-tion into account. In the first test case (represented by solid lines both on theplot and on the scheme) the transition between the two valence levels is also16esonant with the field and this destroys the SIT phenomenon. It suffices to getboth valence levels far apart enough (e.g. 2 ω as in the second case, representedby dashed lines) to recover SIT. We use this last configuration as basis test casefor the following numerical experiments. The simulations are performed using a code based on a finite difference Yeescheme and a relevant choice for the time discretization of the Bloch equation(see [Bid03]). It allows to keep the good properties of the original Yee scheme:second order and explicitness. A splitting scheme, first described in [BBR01],allows to preserve positiveness at the discrete level. It is strongly based onthe exact solution given by Equation (3). It has been adapted to include alsoCoulomb interaction, still preserving positiveness, but at the cost of a loss ofapproximation order, which becomes one.Integrating the zero intra-band coherences assumption is a priori a problemsince it destroys the Liouville structure and an exact solution is no more avail-able. It is however possible to solve a Liouville-like equation and set artificiallyintra-band coherences to zero. This adds a step at each time iteration but allowsto preserve the general structure of the numerical code.To determine the right envelope amplitude for numerics, we use the argumentof [ZAG95]: in practice the input pulse is cut off on an interval t ∈ [ − τ, τ ],therefore the numerical area is A n = mτ E [arctan(sinh t )] − = mτ E (0 . π ) , which slightly changes the value of E . To include Coulomb interaction in the SIT test case, we have to give valuesto the coefficients given by Equation (7). Their exact computation is not inthe scope of the present paper. We choose to take them of the same order R , taking into account the symmetries described by Property 2, but not equal(which would lead to vanishing Λ c , Λ v and κ v ). The test is performed using thefull Coulomb terms (no vanishing intra-band coherence assumption).For small values of R (see the evolution of ρ c11 described for R = 10 − in Figure 2(a), solid plot) SIT is only slightly affected. In this figure, thedashed plot corresponds to the reference case ( R = 0) and is the same as thedashed plot of Figure 1(b). The effect is clearer for stronger values of R (e.g., R = 3 × − in Figure 2(b)). The total inversion is prevented by Coulombinteraction. Although inversion is not complete, the return to zero of ρ c11 isobserved in this test case. We notice that the first valence level (which is notsupposed to take part in the SIT experiment) is slightly populated, which makesthis test case not really a two-level experiment.17a) R = 10 − (b) R = 3 × − Figure 2: Impact of Coulomb terms for two interaction strengths. (a) R =10 − – comparison with the Coulomb-free case; (b) R = 3 × − – excitationof the conduction level and of the first valence level. We first test the impact of the vanishing intra-band coherence assumption onthe Coulomb-free model. We always use the same experimental setting (seeFigure 3(a)) and this assumption amounts to taking ρ v12 = 0.conduction ρ c11 ρ v22 ρ c − v12 ρ v11 ρ v12 ρ c − v11 valence (a) (b)Figure 3: Impact of vanishing intra-band coherences on the Coulomb-free model.The result is displayed in Figure 3(b). The final equilibrium state for matteris slightly changed and ρ c , which is given by the solid curve, does not eventuallyreturn to zero. Inversion is not total.Now we combine both Coulomb interaction and the vanishing intra-bandcoherence assumption. If R is low (e.g., 10 − ), the result is not much affectedby this assumption and is essentially the same as that plotted on Figure 2(a).In the case when R = 3 × − , we see in Figure 4 that the final equilibrium18igure 4: Impact of vanishing intra-band coherences on the full model for R =3 × − .state is not physical ( ρ v11 < R is low. Whenwe are close to the SIT experiment we have a typical two-level phenomenon: ρ v11 (cid:39) • the positiveness of each diagonal term (populations), • the estimation of coherence by populations, here: | ρ c − v12 | ≤ ρ c11 ρ v22 , if thesecond valence level would be the only relevant one (see [BBR01]).Setting intra-band coherences to zero within the numerical process does notaffect these properties, and the iteration used in the proof of the positiveness ofthe density matrix applies.But for a three-level system (and the case R = 3 × − is a true three-level case), the positiveness of the matrix involves some more properties, whichare affected by setting ρ v12 to zero. Although trace is still conserved, the posi-tiveness of the population is affected. We would of course have the same resultwith a dedicated code where intra-band coherences would simply not be com-puted. Besides the effect is clear enough in Figure 4 not to be attributed tosimple round-off errors (there are only 600 time-steps in this computation).The conclusion is that even if intra-band coherences seem not to be very rele-vant for some physical applications, it is very important to include them in the19athematical description and in the numerical computation to keep the naturalmathematical structure of the density matrix. Remark . In absence of electromagnetic field the evolution equation for theconduction electrons is reduced toi (cid:126) ∂ t ρ c11 = 0 . This is only due to the fact that there is only one conduction level in our testcase, and does not depend on the intra-band coherence vanishing assumptionor on specific values of the Coulomb coefficients. Hence, when the populationof the conduction level has been set into a non-physical state (and this is dueto the intra-band coherence vanishing assumption), it remains in this state forever.
In this paper, Bloch-type equations have been derived considering Coulombeffects in quantum dots. We have shown analytically and numerically thatCoulomb effects are not negligible in some quantum dot structures, and wehave given the link between mathematical properties and physical relevancy ofthe Bloch model and more specifically in the treatment of intra-band coherences.Then this model has been coupled with the description of laser propagation inthe quantum dot structures, leading to a Maxwell-Bloch system for which wehave studied the global Cauchy problem. This system has been implementednumerically and simulations have been performed on a self-induced transparencytest-case. In particular, we have tested the impact of Coulomb parameters andintra-band coherences. We have illustrated numerically that the modificationof the equation structure when intra-band coherences are neglected can lead tonon-physical solutions. Further work will include additional effects in the sameBloch-type framework.
Acknowledgements
LJK is partner of the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) fundedby the French program Investissement davenir. The authors wish to thank EricDumas for fruitful discussions.
References [AE87] Allen, L. and Eberly, J.H., ”Optical resonance and two-level atoms”,(Dover, 1987).[BBFB+04] Besse, C., Bid´egaray-Fesquet, B., Bourgeade, A., Degond, P., andSaut, O., ”A Maxwell–Bloch model with discrete symmetries for wavepropagation in nonlinear crystals: an application to KDP,” MathematicalModelling and Numerical Analysis, , 321 (2004).20Bid01] Bid´egaray, B., Contributions `a l’´electromagn´etisme dans le domainetemporel. Mod´elisation classique et quantique en optique non lin´eaire,
Ha-bilitation thesis, Universit´e Paul Sabatier, Toulouse, France (2001).[Bid03] Bid´egaray, B., ”Time discretizations for Maxwell–Bloch equations,” Nu-merical Methods for Partial Differential Equations, , 284 (2003).[BF06] Bid´egaray-Fesquet, B., ”Hi´erarchie de mod`eles en optique quantique. DeMaxwell–Bloch `a Schr¨odinger non-lin´eaire,” volume 49 of Math´ematiqueset Applications , (Springer, 2006).[BF10] Bid´egaray-Fesquet, B., ”Positiveness and Pauli exception principle inraw Bloch equations for quantum boxes,” Annals of Physics, , 2090(2010).[BBR01] Bid´egaray, B., Bourgeade, A., and Reignier, D., ”Introducing physicalrelaxation terms in Bloch equations,” Journal of Computational Physics, , 603 (2001).[Dum05] Dumas, E., ”Global existence for Maxwell–Bloch systems,” Journal ofDifferential Equations, , 484 (2005).[DS12] Dumas, E. and Sueur, F., ”Cauchy problem and quasi-stationary limitfor the Maxwell–Landau–Lifschitz and Maxwell–Bloch equations,” Annalidella Scuola Normale Superiore de Pisa, Classe di Scienze, XI , 503 (2010).[GH02] Gehrig, E. and Hess, O., ”Mesoscopic spatiotemporal theory forquantum-dot lasers,” Physical Review A, , 033804 (2002).[GS70] Gibbs, H.M. and Slusher, R.E., ”Peak amplification and breakup ofa coherent optical pulse in a simple atomic absorber,” Physical ReviewLetters, , 638 (1970).[HK09] Haug, H. and Koch, S.W., ”Quantum Theory of the Optical and Elec-tronic Properties of Semiconductors,” Fifth edition, (World Scientific,2009).[HK96] Hess, O. and Kuhn, T., ”Maxwell–Bloch equations for spatially inhomo-geneous semiconductor lasers I. Theoretical formulation,” Physical ReviewA, , 3347 (1996).[KR92] Kuhn, T. and Rossi, F., ”Monte Carlo simulation of ultrafast pro-cesses in photoexcited semiconductors: Coherent and incoherent dynam-ics,” Physical Review B, , 7496 (1992).[MH67] McCall, S.L. and Hahn, E.L., ”Self-induced transparency by pulsedcoherent light,” Physical Review Letters, , 908 (1967).[Wic50] Wick, G.-C., ”The evaluation of the collision matrix,” Physical Review, , 268 (1950). 21ZAG95] Ziolkowski, R.W., Arnold, J.M., and Gogny, D.M., ”Ultrafast pulseinteraction with two-level atoms,” Physical Review A,52