[ Home ] [ News ] [ Contact ] [ Search ]

 About us
 Useful Links
 Publications
 Copyright notice
 nextnano³ documentation

 -> Download Software
 * password protected

 

 
 

Generate new band structure

We use van de Walle model to generate the new band structure that results from strained heterostructures. A strain causes splitting of degenerate bands and a shift in their energies.

Please also have a look at the FAQ section about strain and deformation potentials and their effects on the band shifts.

Van de Walle (PRB 39, 1871 (1989)) describes the dependence on strain of the valence and conduction band edges (zincblende). In this paper strain is decribed along [001] and [111] directions. In a general case, when we have [a b c] growth direction you will find the equations below (will have to be written down!).

I made some derivations in order to understand what the single band model for holes is, when we have an arbitrary growth direction. The results may be summarized as follows (I follow a paper of Chao and Chuang PRB 46, 4110):

  • 1) In the arbitrary growth direction strain tensor eij can have 6 independent components.
  • 2) Starting point is the Hamiltonian of Bir and Pikus, that takes into account eij.
  • 3) We shall neglect split-off band.
  • 4) In the case of [001] growth direction, when strain tensor has diagonal components only with exx=eyy<>ezz, we find that there are two bands, each of them is twofold degenerated. Effective mass tensor is diagonal in the crystal system and does not depend on strain.
  • 5) But in the case of arbitrary growth direction the situation is more complicated:
    - There are two bands, as it was before, each of them is twofold degenerated.
    - Effective mass tensor is not diagonal in the crystal system. Its elements depend on strain.

Taking into account all above, I'd like to ask you: Can we use this model in nextnano3?
It seems to me, that we need three modifications only:

  • 1) Correct dependence of the band extrema on strain are implemented already.
  • 2) Dependence of the effective mass on strain - not very difficult.
  • 3) Nondiagonal effective mass tensor - this is more difficult, because it requires coordinate transformation of such terms as kxky, kykz and kxkz. But such problems must be solved already in k.p. As for the discretization of nondiagonal terms, I think that this is already done, because even if in the crystal system we don't have them, they appear in the simulation system.

In the MODULE mod_vb_abs_energy band edges are calculated starting from k.p Hamiltonian. So, band edges are correct for any orientation. In order to get the band edge shift due to strain, we diagonalize numerically at each grid point the strain dependent k.p Hamiltonian for k=0.

k.p parameters in the database are used to determine band shifts due to strain (MODULE mod_vb_abs_energy). Program computes the k.p and strain dependent Hamiltonian for each grid point, diagonalizes it, and takes the eigenvalues for the band edges. So, band edges for single band calculations are k.p parameter dependent.

Last days I studied how 1-band Schrödinger is implemented in nextnano3. I tried to understand how band edges are determined. As far as I understood, in order to determine shifts due to strain, you diagonalise k.p Hamiltonian for every grid point. Eigenvalues are taken for the band edges.
Despite of the fact that FUNCTION vb_abs_energy(grid_point,band_number) has detailed comments that have helped me a lot, there are some points that I'd like to clarify.

  • 1) What are the basis vectors that are defined by CALL get_basis_vectors(strain_tensor_cxyz,basis_vec,Spin_orbit_HamM)?
  • 2) Can you please explain how do you choose energy of "heavy" hole and "light" hole among the eigenvalues. I'd be very happy to know just a physical algorithm of yours.

I diagonalize the k.p Hamiltonian including deformation potentials at k=0 at every spatial grid point. The eigenvalues are the local valence band edges. The problem now is to find out which eigenvalue belongs to which band:
What I do is, I project the eigenstates onto the spin-orbit eigenstates named as basis_vec. One problem though is the quantization axis of these states which is the z axis by default. Since the strain is breaking the symmetry the states are getting mixed. I tried to overcome this problem for some extent by trying to find some symmetry axis of the strain tensor. Then, anticipating that this is the best possible quantization axis, I rotate the basis_vec such that the quantization axis coincides with the symmetry axis of the strain. Then I project the eigenstates from the k.p diagonalization onto the rotated basis_vec and just make a majority vote: If a state is 60% heavy and 40 % light hole I call it "heavy hole". I checked this procedure on some structures but mainly in 1D, so maybe you can tell me if it works properly in your case.

I'm studying your subroutine vb_abs_energy. Let me ask a "technical" question.
In the subroutine get_basis_vectors you rotate eigen vectors of the spin-orbit coupling Hamiltonian. Right?
In principle, eigenvector is a spinor. You independently rotate the first 3 components of it and the last 3 components of it. Am I right?
It looks like (maybe I'm wrong) that you don't rotate spin part of the wave function.
As for me this is a mistake. As far as I remember from Landau & Lifshitz, when you rotate axis you also mix different spins.

 

I made a routine that can calculate effective mass and band shifts of holes for zinc-blende. My idea is very similar to what Matthias has done, but there are some changes. There are two possible situations:

  • a) There is no strain, or strain ellipsoid is a sphere.
  • b) otherwise

The problem: Diagonalize strain dependent k.p Hamiltonian and distinguish heavy, light and split-off hole states.
Solution: Distinguish them by projecting them on the states |j,m>
Quantization axis is chosen in such a way:
case a) z-axis coincides with k vector;
case b) quantization axis coinsides with the principal axis of the strain tensor that corresponds to the maximum eigenvalue.
In situation a) we cannot talk about effective-mass tensor because of warping. However, it is possible to define effective mass along any direction.
I hope, of course, that k.p calculation will work fine.
But now I'd like to plug this routine into the code and see what will be the result. For that I need routine that can discretise 1/mxy terms.

 

Let me inform you about a test that I performed in order to understand if band shifts due to strain are calculated correctly.
The structure that I consider is a 1D quantum well GaAs-4K/InAs-4K/GaAs-4K.
In order to check the program, I wrote my own subroutine and compared the band shifts.
I considered 2 growth direction: [0 0 1] and [3 1 1].
Let us compare band edges in InAs (inside the well, where we have nonzero strain):
In the [001] direction our results are the same:
Matthias, hh,lh,s/o 0.286451116426106 8.643207818113535E-002 -0.568763735470029
Michael,  hh,lh,s/o 0.286451116426106 8.643207818113536E-002 -0.568763735470029
But for [311] growth direction they are different:
Matthias, hh,lh,s/o 0.304729495783759 9.995444713679591E-002 9.995444713679591E-002
Michael,  hh,lh,s/o 9.995444713679575E-002 0.304729495783759 -0.561455237788657
Namely, your hh is equal to my lh, your lh and s/o are the same (!) and equal to my hh. I'd like to say that your results are not accurate enough, on my opinion. If you don't agree, please, let me know. Basically, version of Matthias produced correct band edge for heavy and light hole, but noncorrect for s/o band (it was the same as for the light hole).

 

Module that I sent you executes subroutine of Matthias, then executes my subroutine, prints two values to compare and returns values of Matthias for the function value. So, it only compares. I left it so because I didn't want to change code too much. You know that my works for zincblende only yet. When you have time, you change it in such a way that for wurzite it takes Matthias value, for zb it takes mine. As for me, I have my version where my values are taken instead of the values of Matthias. After your make a new release you should update tutorial 1D #5.

--

I send you my routine. Actually, I wanted to spend some time and made also the modification for wurtzite. So, you will get:

  • 1) mod_kp2sb.f90
    It contains subroutines that transfer parameters from your code to my routine, subroutine that calculates band shifts and effective masses.
  • 2) I send you module vb_abs_energy.f90.
    I added lines in order to compare your and my results:
    ! my routine
    call init_strain(strain_tensor_cxyz) ! passes strain to my code
    as = (2*ms1+ls1)/3                   ! conversions to "traditional" parameters
    bs = (ls1-ms1)/3
    ds = ns1/sqrt(3d0)
    print *,'d',d
    call init_LMN_param( L,M,N,as,bs,ds,d) !
    passes k.p - parameters to my code
    call get_shifts(hh_energy1,lh_energy1,so_energy1) ! calculates band edges
    print *,'-----------------------------------------------------'
    print *,'Matthias, hh,lh,s/o',hh_energy, lh_energy, so_energy
    print *,'Michael,  hh,lh,s/o',hh_energy1,lh_energy1,so_energy1
    print *,'-----------------------------------------------------'
    !
    end of my routine

The second thing that I'd like to do is to implement effective mass calculation.

Strain also enters into the single-band calculation. There are two effects: Band shifts and effective mass changes.
As for the band shift I used my routine that is, presumably, correct. As for the effective mass, I calculated their dependence on strain and put the values into the input file. For that I needed to put them in the input file.
I believe that there should be a flag, that will allow the user to choose between constant effective masses (taken from database) and effective masses calculated from k.p parameters. Of course, one needs nondiagonal masses.
> So for 311 the strain is "more important" than for 001. Maybe this could be a reason for the difference.
> > I solved Schrödinger equation for a rather large well with 1D k.p. After that I calculated effective mass from k.p parameters and solved the same in the effective mass approximation. Since it was 1D, I could use isotropic masses.
> > The interesting thing was that for [001] direction results for the 1st state were almost the same, but for [311] they were very different.
> > May be it should be so.
> > I have an opinion, that when I touch something with [311] direction
> > it can be wrong. This was with strain and with band edges. We have to check k.p.

 

main calls SUBROUTINE gen_new_band_struc.

  • MODULE mod_gen_new_material_band_struc
    • SUBROUTINE gen_new_band_struc
      uses SUBROUTINE cb_new_bands to generate new conduction band structure for each material number to have the same number of bands within one material number.
      Initializes structures: new_mat_band_struc (mod_new_band_struc)
                                     old_mat_band_struc (mod_old_band_struc)
      Note: This subroutine only affects the conduction bands.
      SUBROUTINE has to be called at the beginning of the programm. It also uses routines to free memory (free_mem, free_mem2, free_mem3).
    • SUBROUTINE gen_band_struc
      generates new_mat_band_struc out of old_mat_band_struc
      Note: This subroutine only affects the conduction bands.
    • SUBROUTINE check_consistency
      Checks if the new bands contain the same minima in both structures.
      If necessary it generates additional bands in order to get a consistent band structure.
  • MODULE mod_new_band_struc and MODULE mod_old_band_struc
    New conduction band structure which results from the lifting of degeneracies by applied strain.

    new_mat_band_struc sorted by new band numbers and new minima.
    old_mat_band_struc represents same new band structure as new_mat_band_struc but sorted by old bands and old minima.

    The structure is initialized by subroutine gen_new_band_struc in mod_gen_new_material_band_struc.
  • MODULE mod_cb_new_bands
    SUBROUTINE cb_new_bands
    gets the strain tensor and determines the new number of conduction bands by using the unaxial and absolute (=hydrostatic) deformation potentials.
    If absolute energies of unstrained bands are equal, some modifications are needed to keep at least the same number of new bands as old bands.
    It sorts into new bands and generates new band structure.
    It also uses FUNCTION cb_c_kxyz which provides the k coordinates of the minimum in the crystal system.
    This module also contains the following:
    SUBROUTINE free_mem     : Free allocated memory for new_band_struc.
    SUBROUTINE free_mem2    :
    Free allocated memory for old_band_struc.
    SUBROUTINE free_mem3    :
    Free allocated memory for old_band_struc.
    SUBROUTINE error
    SUBROUTINE energy_sort  :
    Sorts array elements in order of energy.
    FUNCTION   get_max_array:
    Returns the index of the largest element of this array.

    This subroutine is called only once in SUBROUTINE gen_new_band_struc in MODULE mod_gen_new_material_band_struc.
    Compare with FUNCTION cb_new_energy.
  • MODULE mod_cb_c_kxyz
    FUNCTION cb_c_kxyz
    provides the k coordinates of the minimum in the crystal system.
    This function is used by SUBROUTINE cb_new_bands and FUNCTION cb_new_energy.

 


MODULE input_block_two:

The following module is used by SUBROUTINE valence_band in MODULE input_block_two which returns for material grid point the value of the valence band edge for all bands:

  • MODULE mod_vb_abs_energy
    FUNCTION vb_abs_energy
    returns band edge in eV.
    It first gets the strain tensor in the crystal system, the deformation potentials for zincblende and wurtzite and its k.p parameters. The diagonalization of the Hamiltonian is done one step away from the Gamma point to avoid degeneracy. We take two energy values for kx=0 and kx=100 to fit parabola to valence bands.
    SUBROUTINE get_basis_vectors
    diagonalizes the strain matrix, normalizes the eigenvectors and gets the high symmetry direction.

The following modules are used by SUBROUTINE conduction_band in MODULE input_block_two which returns for material grid point the value of the conduction band edge for all bands:

  • MODULE mod_cb_new_energy
    FUNCTION cb_new_energy
    gets the strain tensor and determines the new conduction band energy by using the unaxial and absolute (=hydrostatic) deformation potentials.
    In direct semiconductors (e.g. GaAs) the conduction band minimum lies at the Gamma point and is thus only affected by the absolute (hydrostatic) deformation potential which shifts the band energy accordingly. In indirect semiconductors like Si and Ge the degenerate minima are lying at points along the Delta axis or at the L points respectively. In cubic crystals these symmetry points are equivalent minima in k space and are thus degenerate. These degeneracies are lifted by uniaxial strain. Due to symmetry arguments for both the minima at X and L only two deformation potentials are needed because stress applied along [111] acts on all X minima in the same way and thus no splitting occurs. The same applies to the minima at L and stress along [001]. The perturbation term looks like this:
           He = Xid Tr[eij] + Xiu(k.eij.k).
    Here, k gives the position of the minimum in k space. The shift of the average conduction band energy is given by:
           deltaEc,av = ( Xid + 1/3  Xiu )Tr[eij] = ac Tr[eij]

    av_energy           = cb_energy + abs_defpot_cb  *( exx+eyy+ezz )
    total_energy_of_min = av_energy + uniax_defpot_cb*(  2 kx ky exy + (kx2-1/3)exx +
                                                         2 kx kz exz + (ky2-1/3)eyy +
                                                         2 ky kz eyz + (kz2-1/3)ezz  )
    The splitting of the degenerate minima is proportional to the uniaxial deformation potential Xiu which has different values for the L and for the X valley (XiuDelta or XiuX, XiuL).
    Example: For a uniaxial stress along the [001] direction the X valleys split relative to the average conduction band energy:
           DeltaEc001        =  2/3 XiuX (ezz-exx)
           DeltaEc010,100 = -1/3 XiuX  (ezz-exx)

    This FUNCTION is called only once in SUBROUTINE conduction_band in MODULE input_block_two.
    Compare with SUBROUTINE cb_new_bands.
  • MODULE mod_cb_num_new_bands
    FUNCTION cb_num_new_bands
    [see below]

The following modules are used by SUBROUTINE density_of_state_cond in MODULE input_block_two which returns for material grid point the value of the density of states for each band in [1/m3]:

These modules are used by SUBROUTINE non_parab_cond in MODULE input_block_two:

This function is used by SUBROUTINE get_maxval_bands in MODULE input_block_two which provides maximum number of conduction bands in and maximum number of valence bands (including strain):

  • MODULE mod_get_max_cbands
    FUNCTION get_max_cbands()
    provides maximum number of bands in whole grid.

 


  • MODULE mod_get_alloy_const
    FUNCTION get_alloyc gets alloy concentration.
    This function is never used!!! But it might be useful one day.
  • MODULE mod_cb_new_masses_xyz
    FUNCTION cb_new_masses_xyz
    provides mass tensor in simulation coordinate system for new band structure.
    This function is never used!!! But it might be useful one day.
  • MODULE mod_cb_new_kxyz_min
    FUNCTION cb_new_kxyz_min
    This function is never used!!! But it might be useful one day.

All these functions use the following:

cluster         =  belongs_to_cluster(grid_point)
material        =  cluster_to_material(cluster)
material_number =  material_info(material)%material_number
model           =  material_info(material)%material_model_number
nthm            =  material_info(material)%nth_model
naf             =  material_info(material)%alloy_function_number
ntha            =  material_info(material)%nth_alloy_function
pos             =  get_geo_coordinates(grid_point) !
Get alloy concentration.
alloyc          => alloy(naf,ntha,pos) !
Don't forget to dellocate alloyc !!!

Lots of routines contain these expressions explicitly. By calling one of the "never-called" functions, one should be able to "save" these lines in the code.

 

   
Last modified: 25-Aug-2009   -   optimised for Microsoft Internet Explorer 7®