Square Well Single Electron quantum Dot (SWSED)

The package contained in the file "SWSED.m" calculates to high  accuracy the conduction band structure for a single-electron in a  three-dimensional spherical quantum dot within the model of finite  square-well potential of radius R and energy barrier V. Details description of calculations is given in accompanying  paper "Exact calculation of single-electron states in Si-nanocrystal  embedded in SiO2" translated to Mathematica notebook format at the end of this help file. The following line  starts the package.

<< SolidState`SWSED`

A Square Well Single Electron quantum Dot (SWSED), as drawn in  Figure 1, can be specied by 6 real numbers:

The effective mass approximation is proven to be justified for  the determination of the first levels even for dot radius as small as 15  Å (see reference [2], in the end of this manual). The package  is developped by Tran The Trung and Nguyen Hong Quang at Institute of Physics, 46 Nguyen  Van Ngoc, Cong Vi, Ba Ðình, Hà Noi, Vietnam in  2003.

Energy Level

The main function of the package is EnergyLevel. We can get  information about this function by the standard Mathematica command ? placed before the name of the function.

?EnergyLevel

For input to EnergyLevel, we have 2 list. The first list is the 2  real positive numbers that specifies the SWSED. The second list contains 2 integer l and n, being the quantum angular state and quantum energy state  respectively. Note that l must be equal or greater than 0 and n must be equal or greater than 1. The following line give energy  at l = 1, n = 2 of a SWSED with radius 10 Angstrom, barrier 3.1 eV, inside  material is Si and outside material of SiO2 in eV.

EnergyLevel[{20, 3.1}, {1, 2}]

To generate a demonstration of dependence of energy  level on radius and barrier of the SWSED we can have:

$TextStyle = {FontFamily -> "Arial", FontSize -> 12, 
          FontWeight -> "Bold"};

L = {
Table[{R, 
              EnergyLevel[{R, 1}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 5}]
, 
          Table[{R, 
              EnergyLevel[{R, 3.1}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 5}]
, 
          Table[{R, 
              EnergyLevel[{R, 30}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 5}]
, 
          Table[{R, 
              EnergyLevel[{R, 130}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 5}]
};
<< Graphics`MultipleListPlot`;
MultipleListPlot[L
, 
        PlotJoined -> True
, 
        PlotRange -> All
, 
        Frame -> True
, 
        PlotLegend -> {"Vo = 1 eV", "Vo = 3.1 eV", "Vo = 30 eV
", "Vo = 130 eV"}
, 
        ImageSize -> 430
, 
        LegendPosition -> {0.3, -0.3}
, 
        FrameLabel -> {"R in Å", "Energy l=0 n =1 in \

eV"}
];

When the barrier becomes too large, the solution of energy level approaches 
that of infinite square well quantum dot. For example:

EnergyLevel[{20, 10000000}, {1, 2}]

One get the same result by setting V = ∞:

EnergyLevel[{20, ∞}, {1, 2}]

Energy Structure

Bound state of the single electron in SWSED is the state that has energy  level lower than barrier. If barrier is sufficiently low while quantum energy  state number n is high, the energy level may exceeds barrier and the state is unstable:

EnergyLevel[{20, 1.2}, {1, 2}]

One can find all bound states of the single electron for a given SWSED by  trying all different pairs of quantum number l and n which gives bound energy levels. However, the SWSED Package provides a  built-in function that give exactly this, the function EnergyStructure.

The output of EnergyStructure is a list of 2 elements. The first element  is the list of 6 real positive numbers that specifies the SWSED (R, V, m inside, m outside, ϵ inside, ϵ outside). The second element has depth of 2. The first depth of this list  correponds to quantum angular number l = 0, 1, 2.... The second depth are list of energy levels that has the  same quantum angular number l but different quantum state number n = 0, 1, ... Options of EnergyStructure also include EMass and  Permeability as for EnergyLevel.

EnergyStructure[20, 3.1]

It is much easier to view the whole energy structure graphically, by using 
the function StructurePlot of SWSED package.

StructurePlot[EnergyStructure[20, 3.1]];

The dashed line shows the energy barrier of the SWSED. No bound energy levels 
can exceed this barrier. StructurePlot has the following options

These options are standard Mathematica options for Graphics, Show and Plot functions. Beside, we add 3 options

The option LevelLabel specifies the text or expression label at each level. 
Default value is Automatic, which mean level are labeled by standard quantum 
physics labeling for electron state in atom, as shown above. One can specify 
different LevelLabel by a text/expression or a list of text/expression.

StructurePlot[EnergyStructure[20, 3.1], 
        PlotLabel -> "Different style",  
        PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.01]}, 
        LevelLabel -> Range[13]];

StructurePlot[EnergyStructure[20, 3.1], ImageSize -> 450, 
        PlotLabel -> "Different style",  
        LevelLabel -> "SiO2"];

If no label is needed, one can set LevelLabel to None. Example of the same 
plot as above with different style and having no label is:

StructurePlot[EnergyStructure[20, 3.1], 
        PlotLabel -> "No Label",  
        LevelLabel -> None];

The option LevelShape can take either Automatic or None value and effective 
only when PlotPDF is False. The option PlotPDF will be discussed in the next 
section.

From the graphics of the energy structure, one can find the order of energy 
level, for example above, as 1s < 2p <3d <2s < 4f < 3d ... The package SWSED 
comes with a function that help automatically give the order of the energy 
level structure: StructureOrder.

StructureOrder[EnergyStructure[20, 3.1]]

Note that StructureOrder has no option:

StructureOrder[EnergyStructure[30, 3.1 ], 
      EMass -> {0.27, 0.27}]

The option EMass is reserved to EnergyStructure function:

?EMass

So the energy order for SWSED with R = 30, V = 3.1eV and electron mass both 
0.27 inside and outside can be found by:

StructureOrder[
      EnergyStructure[30, 3.1 , EMass -> {0.27, 0.27}]]

As input to StructurePlot and StructureOrder are output of EnergyStructure, 
we also provide a function EnergyStructureQ, which tests if a data can be 
possibly an output of EnergyStructure

?EnergyStructureQ

EnergyStructureQ[
      EnergyStructure[30, 3.1 , EMass -> {0.27, 0.27}]]

You can use this function in case you are not sure if input to StructurePlot 
or StructureOrder are corrected.

Wave Function & Probability Density Function

Apart from energy level information, one can also find the wave function  of the single electron in a given SWSED at a bound state. This wave function  is a product of spherical harmonic (Subscript[Y, l, m]), with orbital quantum number l (l = 0, 1, 2, …) and azimuthal quantum number m (-l ≤ m ≤ l), and radial part F(r).
	Ψ((Overscript[r,  ->]) ) = (Subscript[Y, l, m])(φ, θ) F(r)		
F(r) depend on the SWSED parameters and  the state, which can be calculated by the function RadialWaveFunction. See  the appended article in the end of this help file for the analytic form of  F(r).

RadialWaveFunction also has the options EMass and Permeability for  specification of material of the SWSED other than Si-SiO2. Example of  wavefunction for the electron inside SWSED1 (in section Energy Level above)  at state l = 0, n = 1 is

WF1 = RadialWaveFunction[{20, 3.1}, {0, 1}]

Plot[WF1[x], {x, 0, 20}, Frame -> True, 
        FrameLabel -> {"Radius in Angstrom", "Radial Wave Function
"}];

A physically interesting derivation from wavefunction is probability density 
function (PDF) of finding random position of the single electron in an unit 
volume at given radial distance to center of SWSED, which is just the square 
of the complete wavefunction normalized, caculable by DensityFunction

DensityFunction also has the options EMass and Permeability for  specification of material of the SWSED other than Si-SiO2. Here is the  example of PDF for the electron inside SWSED1 at state l = 0, n = 1 is

PDF1 = DensityFunction[{20, 3.1}, {1, 1}]

Plot[PDF1[x], {x, 0, 20}, Frame -> True, 
        FrameLabel -> {"Radius in Angstrom", "PDF"}, 
        PlotRange -> All];

The probability density function can also be used for visualisation in 
EnergyStructurePlot. By setting the option PlotPDF as True, 
EnergyStructurePlot will shows the PDF associated with the energy level. 
Because of visual evaluation purpose, currently EnergyStructurePlot is 
implemented to show only the PDF inside the barrier radius of SWSED.

StructurePlot[EnergyStructure[20, 3.1], PlotPDF -> True, 
        ImageSize -> 600, PlotLabel -> "With PDF", 
        PlotStyle -> {Hue[0.6]}];

If one wants to see more clrearly the PDF at, for example, state {l = 3, n = 2}, one can always use DensityFunction:

PDF32 = DensityFunction[{20, 3.1}, {3, 2}];

Plot[PDF32[x], {x, 0, 20}, Frame -> True, 
        FrameLabel -> {"Radius in Angstrom", "PDF"}, 
        PlotRange -> All];

References

[1] Accompanying paper "Exact calculation of single-electron states in Si-nanocrystal  embedded in SiO2" below
[2] J. Sée, P. Dollfus, and S. Galdin, Phys. Rev. B 66  (2002) 193307
[3] S. Wolfram, The Mathematica Book, Cambridge University  Press, (1999)

Authors:

Nguyen Hong Quang (a,b), Tran The Trung (a,c), Johann Sée (b) and 
Philippe Dollfus (b),

Abstract
We present an exact calculation of the single-electron energies  and wave-functions for any bound state in a realistic Si-SiO2spherical  quantum dot,including the material dependence of the electron effective mass.  The influence of dot radius,confinement barrier potential and barrier-to-dot  electron mass ratio on the electronic structure is investigated in detail.  The results show that the energy structure shifts down from some tens to some  hundreds meV compared to that obtained in the simplified model where the  change in effective mass is neglected. Our exact single-electron calculation  is finally used to verify the accuracy of the results obtained from a  numerical approximate method developed to treat many-electron systems.
PACS  numbers:71.15.La,73.20.DxKeywords:A.Semiconductors; A.Nanostructures;  D.Electronic states
Corresponding author. 	E-mail :  philippe.dollfus@ief.u-psud.fr
				Fax : +33 (0)1 69  15 40 20

1. Introduction

Efforts to integrate sub-10 nm-scaled CMOS devices will face physical 
limitations in the near future. To continue further increases in density with 
low power consumption, the trend is now to develop new device concepts 
turning the size reduction and quantum effects to advantage. In this field, 
one of the most promising idea is to use the Coulomb blockade mechanism in 
conducting island/tunnel junction systems to accurately control the current 
in single electron devices (SEDs).Self-assembled silicon nanocrystals (NCs) 
embedded in SiO2 appear as a good basic structure for the realization of SEDs 
likely to operate at room temperature[1]. The compatibility of such NCs with 
Si technology makes them good candidates for actual applications in 
combination with conventional CMOS circuits. The concept of multi-dot memory 
using Si NCs as a floating gate has been demonstrated[2-4]. It could be 
developed to design nonvolatile multibit memory cells of very high 
scalability and reliability with low power consumption. The Coulomb blockade 
and quantization effects have been used in memory device with self-aligned 
doubly-stacked dots to improve the retention time[5] and demonstrated in 
single electron transistors (SETs) at room temperature[6].

The realisation and optimisation of SEDs requires accurate computational 
simulation of single electron charging effects to improve our understanding 
of the underlying physics. One important preliminary task is to calculate 
exactly the energy band structure inside the semiconductor NC where quantum 
effects become important. Several numerical techniques have been developed 
using either a tight-binding description[7] or an effective mass 
approximation[8-10]. The self-consistent Poisson-Schrödinger 
calculation based on the Hartree approximation or the density functional 
theory is convenient to determine the ground state of a many-electron system 
in the dot as in Ref.10. Simple analytical models have been reported in the 
case of a single electron in a spherical quantum dot[11,12]. However,in these 
works the authors assumed the electron effective mass to be the same in the 
dot as in the surrounding barrier,which can make the results inappropriate to 
realistic applications. The present work deals with the exact calculation of 
single-electron energy and wave-function for all bound states in a realistic 
Si NC embedded in SiO2by including the change in effective mass from one 
material to the other. We give a detailed investigation of the energy 
structure for a wide range of dot radius,barrier height,and barrier-to-dot 
electron mass ratio. Our results are finally used to verify the accuracy of 
the results obtained using the numerical method developed in Ref.10. It 
should be noted that the calculation is exact, fast and applicable for a 
spherical square-well quantum dot system of any size and potential barrier in 
any material system.

2. Model for Exact Calculation

We study the conduction band structure for a single-electron in a  spherical QD within the model of finite square-well potential of radiusRand  energy barrier. Unless otherwise indicated, the effective atomic units are  used in the whole article, i.e.the unit of mass is the free electron mass m0, the unit of length is the effective Bohr radius, \(a\_B\^*\)=(Subscript[a, B])/(Subscript[m, Si]) (with (Subscript[a, B]) = 0.529 Å) and the unit of energy is the effective Hartree  \(TraditionalForm\`\(E\_Hart\^*\) = m\_Si\ E\_Hart\) (with  (Subscript[E, Hart]) = 27.212 eV).The electron effective mass in silicon (Subscript[m, Si]) is approximated as the harmonic mean of the longitudinal and transverse  masses in Δ valleys, i.e. (Subscript[m, Si]) = 3 (((1/(Subscript[m, l])  +  2/(Subscript[m, t])))^(-1))= 0.27. On the basis of Weinberg experimental work [13] the barrier height  (Subscript[V, 0]) and the electron effective mass in SiO2 (Subscript[m, Si02]) are chosen equal to 3.1 eV and 0.5, respectively. These values are  frequently used in tunnelling calculations in the Si/SiO2 system (see[14] and  references therein). 
In the effective mass approximation, which proves to  be justified for the determination of the first levels even for dot radius as  small as 15 Å[15], the single electron Schrödinger  equation is given by
	[-1/2 (Subscript[∇, (Overscript[r,  ->]) ])\((\(1\/\(m\^*\)\) ∇\_\(r\& -> \))\) + V((Overscript[r,  ->]) ) ]Ψ((Overscript[r,  ->]) ) = E Ψ((Overscript[r,  ->]) )	 (1)
where
	V((Overscript[r,  ->]) ) = {({{0 ,  r < R}, {(Subscript[V, 0]) ,  r ≥  R}})		and 		m* = {({{(Subscript[m, Si ,   r  <  R])}, {(Subscript[m, SiO2]) ,  r ≥ R}})		(2)
Due to spherical symmetry of the system, the wavefunction may be  written in the form
	Ψ((Overscript[r,  ->]) ) = (Subscript[Y, l, m])(φ, θ) F(r)		(3)
where (Subscript[Y, l, m]) is a spherical harmonic with orbital quantum number l (l = 0, 1, 2, …) and azimuthal quantum number m (-l ≤ m ≤ l). It can be easily shown from Eq. (1) and (2) that the differential  equations for the radial part F(r) inside and outside the dot are just the  equations of spherical ordinary and modified Bessel functions, (Subscript[j, l])(r) and (Subscript[k, l])(r), respectively. The solutions, which are finite at the center of the  dot (r =0) and tend to zero whenrincreases indefinitely are  given by
-inside the dot:	
	F(r) = A (Subscript[j, l])(α r) = A (Sqrt[(π/(2  αr))]) (Subscript[J, l + 
                1/2])(αr);   with α = (Sqrt[E])		(4)
-outside the dot:
	F(r) = B (Subscript[k, l])(βr)	= B(Sqrt[(π/(2  βr))]) (Subscript[K, l + 
                1/2])(βr) ;   with β = (Sqrt[γ((Subscript[V, 0]) - E)]), γ = ((Subscript[m, SiO2])/(Subscript[m, Si]))	(5)
where (Subscript[J, l + 1/2])(αr) is the Bessel function of the first kind and (Subscript[K, l + 1/2])(βr) is the modified Bessel function of the second kind. The  normalization coefficient B is defined from the continuity condition imposed  on the wavefunction at the dot surface: B = A (Subscript[j, l])(αR)/(Subscript[k, l])(βR). From the normalization of the wavefunction, A is given by
	A  = (Sqrt[T1 + T2]), where T1 = \(TraditionalForm\`∫\_0\%R\(\(|\)\(\(j\_l\)(αr)\)\( | \^2\)\(\(r\^2\) dr\)\)\), T2 = \(TraditionalForm\`\(\(|\)\(\(\(j\_l\)(αR)\)\/\(\(k\_l\)(β 
                 R)\)\)\( | \^2\)\(∫\_R\%∞\(\(|\)\(\(k\_l\)(βr)\)\( | \^2\)\(\(r\^2\) dr\)\)\)\)\)	(6)
To determine the energy levels in the dot, we use the  transcendental equation resulting from the boundary condition imposed on the  derivative of the wave-function at the dot surface, i.e., by taking into  account the change in effective mass between Si and SiO2,
	(Subscript[j, l])'(αR) = \(TraditionalForm\`\(1\/γ\) \(\(j\_l\)(αR)\)\/\(\(k\_l\)(βR)\)\)(Subscript[k, l])'(βR)	(7)
Using the following formulas for derivatives
	 (Subscript[j, l])'(r) = (Subscript[j, l - 1])(r) - ((l + 1)/r)(Subscript[j, l])(r);   and (Subscript[k, l])'(r) = -(Subscript[k, l - 1])(r) - ((l + 1)/R)(Subscript[k, l])(r)		(8)
we get the following equation from which the discrete values  of electron energy E = (α^2) = (Subscript[V, 0]) - (β^2)/γ can be obtained:
	γ\(TraditionalForm\`\([\(α\ \(\(j\_\(l - 1\)\)(αR)\)\)\/\(\(j\_l\)(αR)\) - \(l + 1\)\/R]\)\  = \ \(-\(\(β\ \(\(k\_\(l - 
                          1\)\)(\(\@\(γ(V\_0 - α\^2)\)\) 
                      R)\)\)\/\(\(k\_l\)(\(\@\(γ(
                          V\_0 - α\^2)\)\) 
                    R)\)\)\) - \(l + 1\)\/R\)		(9)
Note that, if we ignore the difference of electron effective mass  between the different materials, the equation (9) simplifies to
	 \(TraditionalForm\`\(\(αj\_\(l - 1\)\)(αR)\)\/\(\(j\_l\)(αR)\)\)=\(TraditionalForm\`\(\(βk\_\(l - 1\)\)(βR)\)\/\(\(k\_l\)(βR)\)\), with (α^2)  +  (β^2)= (Subscript[V, 0])		(10)
which has been used by various authors [11-12]. 
Once the  energies (Subscript[E, l, n]) = ((Subscript[α, l, n])^2) which satisfy Eq. (9) are determined, the exact wave-functions of bound  states are given by	(Subscript[Ψ, l, m, 
            n])((Overscript[r,  ->]) ) = {	({{(Subscript[Aj, l]) ((Subscript[α, l, n]) r) (Subscript[Y, l, 
                    m]) (φ, θ) ,  r < R}, {\(A \(\(\( j\_l\) \((\(α\_\(l, n\)\) 
                          R)\)\)\/\(\(k\_l\) \((\(\@\(γ \((V\_0 - α\_\(l, n\)\^2)\)\)\) 
                          R)\)\)\) \(k\_l\) \((\(\@\(γ \((V\_0 - α\_\(l, n\)\^2)\)\)\) r)\) \(Y\_\(l, 
                    m\)\) \((φ, θ)\)\ , \ \ r ≥ R\)}})   	(11)
with A defined by Eq.(6). Given the radius R, the confinement  barrier height (Subscript[V, 0]) and the mass ratio γ, Eq.(9) is solved to determine all allowed  bound statesin the dot. The solution has been carried out by computation code  written in Mathematica environment [16], yielding high accuracy within manageable calculation  time.

3. Result and Discussion

We are interested in spherical Si-SiO2 quantum dots with typical radius R in 
the range 10Å-60Å. The barrier energy V0 is 3.1 
eV, and the mass ratio is γ =1.85. However, for 
purposes of comparison and to analyse the influence of these parameters, V0 
and γ have been considered to vary in the range 1 eV-30 eV and 0.1 - 
2.5, respectively. We plot in Fig.1 the R dependence of the 1s-level for 
γ =1 and for V0 =
 1 eV, 3.1 eV and 30 eV.

L = {
Table[{R, 
              EnergyLevel[{R, 1}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 2}]
, 
          Table[{R, 
              EnergyLevel[{R, 3.1}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 2}]
, 
          Table[{R, 
              EnergyLevel[{R, 30}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 2}]
};

<< Graphics`MultipleListPlot`;
MultipleListPlot[L
, 
        PlotJoined -> True
, 
        PlotRange -> All
, 
        Frame -> True
, 
        PlotLegend -> {"Vo = 1 eV", "Vo = 3.1 eV", "Vo = 30 eV
"}
, ImageSize -> 400
, 
        LegendPosition -> {0.2, 0}
, 
        LegendSize -> {0.6, 0.4}
, 
        FrameLabel -> {"Dot radius R in Å", "Energy 
in eV"}
, SymbolShape -> None
, 
        PlotLabel -> "Figure 1"
];

L = {
Table[{V, 
              EnergyLevel[{8, V}, {0, 1}, EMass -> {0.27, 0.27}]}, {V, 
              1, 20, 1}]
, 
          Table[{V, 
              EnergyLevel[{10, V}, {0, 1}, EMass -> {0.27, 0.27}]}, {V, 
              1, 20, 1}]
, 
          Table[{V, 
              EnergyLevel[{15, V}, {0, 1}, EMass -> {0.27, 0.27}]}, {V, 
              1, 20, 1}]
, 
          Table[{V, 
              EnergyLevel[{30, V}, {0, 1}, EMass -> {0.27, 0.27}]}, {V, 
              1, 20, 1}]
};

<< Graphics`MultipleListPlot`;
\(\(Show[{
Graphics[{Text["R = 8 Å", {5, 1.6}]
, 
              Text["R = 10 Å", {7, 
                  1.2}]
, 
              Text["R = 15 Å", {10, 
                  0.6}]
, 
              Text["R = 30 Å", {14, 
                  0.2}]}]
, 
          MultipleListPlot[L
, 
            PlotJoined -> True
, 
            SymbolShape -> None
, 
            DisplayFunction -> 
              Identity
]}
, 
        DisplayFunction -> $DisplayFunction
, 
        PlotRange -> {0, 1.8}
, 
        Frame -> True
, 
        ImageSize -> 400
, 
        FrameLabel -> {"Barriere height V0 in eV", "Energy in 
eV"}
, PlotLabel -> "Figure 2"];\)\)

In Fig.2 theV0 dependence of the same level is shown for γ
 = 1 and R =
 8 Å, 1 
Å, 15 Å, 30 
Å. These figures verify that the energy increases when the 
radius decreases and when the barrier potential increases, which is a common 
feature of quantum confinement. Note that the dependence on this parameters 
is all the stronger that the radius and the barrier energy become smaller.

To see the effect of the mass difference between the dot and the surrounding 
barrier we depict in Fig.3 the radius dependence of 1s-level for V0
 = 3.1 eV γ
 = 0.5,1,2.5. It shows that the energy 
decreases when the mass ratio is larger.

L = {
Table[{R, 
              EnergyLevel[{R, 3.1}, {0, 1}, 
                EMass -> {0.27, 0.135}]}, {R, 10, 60, 
              2}]
, 
          Table[{R, 
              EnergyLevel[{R, 3.1}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 2}]
, 
          Table[{R, 
              EnergyLevel[{R, 3.1}, {0, 1}, 
                EMass -> {0.27, 0.675}]}, {R, 10, 60, 
              2}]
};

<< Graphics`MultipleListPlot`;
MultipleListPlot[L
, 
        PlotJoined -> True
, 
        PlotRange -> All
, 
        Frame -> True
, 
        PlotLegend -> {"γ = 0.5", "γ = 1", "
γ = 2.5"}
, 
        ImageSize -> 400
, 
        LegendPosition -> {0.2, 0}
, 
        LegendSize -> {0.6, 0.4}
, 
        FrameLabel -> {"Dot radius R in Å", "Energy 
in eV"}
, SymbolShape -> None
, 
        PlotLabel -> "Figure 3"
];

L = {
Table[{V, 
              EnergyLevel[{30, V}, {0, 1}, EMass -> {0.27, 0.135}]}, {V, 
              1, 20, 1}]
, 
          Table[{V, 
              EnergyLevel[{30, V}, {0, 1}, EMass -> {0.27, 0.27}]}, {V, 
              1, 20, 1}]
, 
          Table[{V, 
              EnergyLevel[{30, V}, {0, 1}, EMass -> {0.27, 0.27}]}, {V, 
              1, 20, 1}]
};

<< Graphics`MultipleListPlot`;
MultipleListPlot[L
, 
        PlotJoined -> True
, 
        PlotRange -> All
, 
        Frame -> True
, 
        PlotLegend -> {"γ = 0.5", "γ = 1", "
γ = 1.85"}
, 
        ImageSize -> 400
, 
        LegendPosition -> {0.2, -0.3}
, 
        LegendSize -> {0.6, 0.4}
, 
        FrameLabel -> {"Barrier height V0 in eV", "Energy in 
eV"}
, SymbolShape -> None
, 
        PlotLabel -> "Figure 4"
, 
        Axes -> False
];

A similar feature is observed in Fig. 4 for the barrier  potential dependence of the 1s-level (R =  30Å, and γ =0.5, 1, 1.85). We still  see that at given R and V0 the energy is lower when γ increases;  moreover, the change of energy as a function of γ is all the higher  that V0 is smaller.
To estimate quantitatively the change of the electron  energy depending on whether the difference of effective mass between  materials is taken into account or not, we plot in Figs.5 and  6 the quantity (Subscript[E, γ = 1])-(Subscript[E, γ = 1.85]) as a function of dot radius and barrier potential, respectively.

L = {
Table[{R, 
              EnergyLevel[{R, 3.1}, {0, 1}, EMass -> {0.27, 0.27}] - 
                EnergyLevel[{R, 3.1}, {0, 1}]}, {R, 10, 60, 
              2}]
, 
          Table[{R, 
              EnergyLevel[{R, 3.1}, {1, 1}, EMass -> {0.27, 0.27}] - 
                EnergyLevel[{R, 3.1}, {1, 1}]}, {R, 10, 60, 
              2}]
, 
          Table[{R, 
              EnergyLevel[{R, 3.1}, {2, 1}, EMass -> {0.27, 0.27}] - 
                EnergyLevel[{R, 3.1}, {2, 1}]}, {R, 10, 60, 
              2}]
};

<< Graphics`MultipleListPlot`;
MultipleListPlot[L
, 
        PlotJoined -> True
, 
        PlotRange -> All
, 
        Frame -> True
, 
        PlotLegend -> {"1s", "2p", "3d"}
, 
        ImageSize -> 400
, 
        LegendPosition -> {0.2, 0}
, 
        LegendSize -> {0.6, 0.4}
, 
        FrameLabel -> {"Dot radius R in Å", "Energy 
in eV"}
, SymbolShape -> None
, 
        PlotLabel -> "Figure 5"
];

L = {
Table[{V, 
              EnergyLevel[{30, V}, {0, 1}, EMass -> {0.27, 0.27}] - 
                EnergyLevel[{30, V}, {0, 1}]}, {V, 0.5, 15, 
              0.5}]
, 
          Table[{V, 
              EnergyLevel[{30, V}, {1, 1}, EMass -> {0.27, 0.27}] - 
                EnergyLevel[{30, V}, {1, 1}]}, {V, 0.5, 15, 
              0.5}]
, 
          Table[{V, 
              EnergyLevel[{30, V}, {2, 1}, EMass -> {0.27, 0.27}] - 
                EnergyLevel[{30, V}, {2, 1}]}, {V, 0.5, 15, 
              0.5}]
};

<< Graphics`MultipleListPlot`;
MultipleListPlot[L
, 
        PlotJoined -> True
, 
        PlotRange -> All
, 
        Frame -> True
, 
        PlotLegend -> {"1s", "2p", "3d"}
, 
        ImageSize -> 400
, 
        LegendPosition -> {0.2, 0}
, 
        LegendSize -> {0.6, 0.4}
, 
        FrameLabel -> {"Barrier height V0 in eV", "Energy in 
eV"}
, SymbolShape -> None
, 
        PlotLabel -> "Figure 6"
, 
        AxesOrigin -> {0, 0}
];

The three first levels 1s, 2p and 3d are considered in these graphs. We can 
see clearly that the magnitude of the energy change is larger for small dots 
and small barrier energies. The typical energy change is about some tens meV 
for large dots and reaches more than 100 meV for small 
dots, which cannot be neglected for most practical applications. Despite the 
effective mass approximation should fail at high energy in Si, we plot in 
Fig.7 the full picture of the energy structure of a single electron quantum 
dot Si-SiO2  as a function of the principal quantum numbern(with γ
 =1, Vo=3.1 eV and 
R = 30Å).

ES = 
        EnergyStructure[30, 3.1, EMass -> {0.27, 0.27}];

<< Graphics`MultipleListPlot`
g1 = 
        MultipleListPlot @@ (Join[
              Table[
Table[{j - 0.58, 
                    ES[[2, i, j]]}, {j, 
                    Length[ES[[2, i]]]}]
, {i, 
                  Length[ES[[2]]]}], {PlotJoined -> True, 
                SymbolShape -> PlotSymbol[Box], 
                PlotStyle -> Thickness[0.01], 
                DisplayFunction -> Identity}]);

Filling[l_List, L_Integer] := 
        If[Length[l] >= L, l, 
          Join[l, Table[,{L - Length[l]}]]];
Choping[l_List] := Select[l, NumericQ];
ES[[2]] = 
        Map[Choping, 
          Transpose[Map[Filling[#, 4] &, ES[[2]]]]];
\(\(g2 = 
        StructurePlot[ES
, 
          LevelLabel -> 
            Join[Table["", {6}], {"i ", "k ", "l ", "m  
", "n "}, Table["", {5}], {"h  "}, 
              Table["", {3}], {"d  ", "f ", "g  ", "s  ", "p  "}]
, LevelShape -> None, 
          DisplayFunction -> Identity];\)\)

Show[{g1, g2}, DisplayFunction -> $DisplayFunction, 
        FrameLabel -> {"Principal quantum number n", "Energy in eV
"}, ImageSize -> 400, PlotLabel -> "Figure 7", 
        Frame -> True];

\(\(Show[{StructurePlot[
            EnergyStructure[30, 3.1, EMass -> {0.27, 0.27}], 
            LevelLabel -> None, DisplayFunction -> Identity], 

          StructurePlot[EnergyStructure[30, 3.1], 
            PlotStyle -> {Dashing[{0.01}]}, 
            DisplayFunction -> Identity]}, 
        DisplayFunction -> $DisplayFunction, ImageSize -> 400, 
        PlotLabel -> "Figure 8"];\)\)

It is also plotted in Fig.8 as a function of the orbital numberl (with V0  = 3.1 eV and R  = 30 Å,  and both γ = 1 and γ  = 1.85). If we use the conventional  spectroscopic notation for the electron energy levels, we obtain an order of  energy levels quite different from that of the hydrogen atom. The sequence is  as follows
	(Subscript[E, 1  s]) < (Subscript[E, 2  p]) < (Subscript[E, 3  d]) < (Subscript[E, 2  s]) < (Subscript[E, 4  f]) < (Subscript[E, 3  p]) < (Subscript[E, 5  g]) < (Subscript[E, 4  d]) < (Subscript[E, 6  h]) < (Subscript[E, 3  s]) ...
The difference in the order of levels between the hydrogen atom and  the semiconductor QD is explained by the fact that in hydrogen atom the  Coulomb potential is deeply singular at r =0, while there  is no singularity in the dot. This sequence agrees with the results reported  in Ref.11. We can also note that the separation of the successive energy  levels increases when the numbernincreases, i.e.when we go up to the top of  the dot well. This feature is quite contrary to the case of a normal atom, in  which the separation between energy levels of highly excited states becomes  smaller and smaller whennincreases. 
In Fig.8 we can compare the energy  structure in two cases: by including the material-dependence of mass in the  calculation (dashed lines) or not (solid line). In accordance with results  presented above, the energy levels resulting from a realistic model (γ  =1.85) are lower than using the simplified model (γ  =1). This situation may result in the emergence of new  bound states which have been omitted in the simplified calculation. For  example, the highly excited level 6d and 9k in Fig.8 are only observed for  γ =1.85. Additionally, the energy shift being  significantly dependent on l value, one may observe a change in the sequence of levels (see e.g. 11n  and 7g). By choosing the appropriate radius, barrier energy and mass ratio,  one can expect this situation to happen for lower levels. 
Finally, we would  like to compare our exact calculation with the results obtained from an  approximate numerical method developed to treat many electron systems. Fig.9  presents the radius dependence of the first energy level calculated either by  our exact treatment (solid line) or by self-consistently solving Poisson-Schr ödinger equations using the Hartree approximation [10] (dashed  line). In Fig.10 the corresponding wave functions resulting from both methods  are shown. The excellent agreement obtained for the ground state is observed  for higher levels as well (not shown), confirming the validity of both  methods.

L = {
Table[{R, 
              EnergyLevel[{R, 3.1}, {0, 1}, EMass -> {0.27, 0.27}]}, {R, 
              10, 60, 2}]
, {{10, 0.911872}, {15, 
              0.462641}, {20, 0.277998}, {25, 0.185126}, {30, 0.132016}, {35, 
              0.098850}, {40, 0.076767}, {45, 0.061331}, {50, 0.050119}, {55, 
              0.041720}, {60, 0.035267}}
};

<< Graphics`MultipleListPlot`;
\(\(MultipleListPlot[L
, 
        PlotJoined -> True
, 
        PlotRange -> All
, 
        Frame -> True
, 
        PlotLegend -> {"Exact", "Poison-Schrödinger"}
, 
        ImageSize -> 400
, 
        LegendPosition -> {\(-0.2\), 0}
, 
        LegendSize -> {1. , 0.4}
, 
        FrameLabel -> {"Dot radius R in Å", "Energy 
in eV"}
, SymbolShape -> None
, 
        PlotLabel -> "Figure 9"
, 
        AxesOrigin -> {0, 0}
];\)\)

L1 = 
        Interpolation[{{0.`, 0.04636525739999999`}, {1.2182629`, 
              0.0461782934`}, {2.5380516000000006`, 
              0.0454304374`}, {3.7563145000000007`, 
              0.044308739399999984`}, {5.1776290000000005`, 
              0.04243909939999999`}, {6.801995099999999`, 
              0.040195617399999994`}, {7.8172064`, 
              0.037952135399999996`}, {9.5430516`, 
              0.034586955399999994`}, {11.471995100000001`, 
              0.030286955399999996`}, {13.400891900000001`, 
              0.025613027399999995`}, {14.923732200000002`, 
              0.0218739194`}, {16.649577400000002`, 
              0.017947847399999998`}, {19.187628999999998`, 
              0.0121521354`}, {21.725680600000004`, 
              0.0074782934`}, {24.6697887`, 
              0.0037390993999999998`}, {27.4108919`, 
              0.0013086964000000004`}, {29.8474177`, 
              0.00018695540000000053`}, {33.197629`, 0.00000001}, {40, 
              0}}];

L2 =  
        DensityFunction[{30, 3.1}, {0, 1}, 
          EMass -> {0.27, 0.27}];

<< Graphics`Legend`

\(\(Plot[{\ 10^3\ L2[x], L1[x]}, {x, 0, 40}
, 
        Frame -> True
, 
        FrameLabel -> {"Radial position in the dot r in 
Å", \*""Square wave function [(nm^(-3))]""}
, PlotLabel -> "Figure 10"
, 
        ImageSize -> 400
, 
        PlotLegend -> {"Exact", "Poisson-Schrödinger"}
, 
        LegendPosition -> {\(-0.15\), 0.2}
, 
        LegendSize -> {1, 0.25}
, 
        LegendShadow -> None
, 
        PlotStyle -> {{}, {Dashing[{0.01}]}}];\)\)

Conclusion

We have presented the exact calculation of energy structure of a single 
electron in Si spherical nanocrystal surrounded by SiO2 forming a square well 
with finite barrier potential.Contrary to the case of simplified models used 
so far,the present model takes into account the material dependence of the 
electron effective mass.Fast and accurate calculations can be carried out to 
generate the energy structure for any combination of three inputs:radius of 
the dot,barrier potential of the dot and the ratio of barrier to dot electron 
effective mass.It is shown that the energy resulting from the simplified 
approach (withγ =1) is overestimated for all 
states.The effect of material-dependence of electron effective mass becomes 
all the more significant that the quantum dot size and the barrier energy are 
smaller:the difference between realistic and simplified calculations of 
energy levels may reach more than 100 meV in small dots.Additionally,our 
calculation may result in new energy levels omitted by the simplified method 
and in a modified order of energy levels.As a future extension of this 
work,the basis of Bessel-like wavefunctions determined for a single electron 
could be generalised to treat the effect of Coulomb interaction in 
many-electron quantum dots.Such an approach should be faster to compute than 
the usual self-consistent Poisson-Schrödinger solver.

Acknowledgement

One of the authors (N.H.Quang) wishes to acknowledge financial support of the 
CNRS (PICS program).He is most appreciative for the hospitality shown to him 
while a visitor at the IEF,Orsay.

References

[1] T.Baron,P.Gentile,N.Magnea and P.Mur,Appl.Phys.Lett.79 (2001) 1175

[2] S.Tiwari,F.Rana,H.Hanafi,A.Hartstein and 
E.F.Crabbé,Appl.Phys.Lett.68 (1996) 1377

[3] E.Kapetanakis,P.Normand,D.Tsoulakas,K.Beltsios,Appl.Phys.Lett.80 (2002) 
2794

[4] S.Banerjee,Solid State Commun.125 (2003) 567

[5] R.Ohba,N.Sugiyama,K.Uchida,J.Koga and A.Toriumu,in Tech.Dig.IEDM'00 
(2000) 313

[6] B.Choi,S.Hwang,I.Kim,H.Shin,Y.Kim and E.Kim,Appl.Phys.Lett.73 (1998) 3129

[7] Y.Niquet,C.Delerue,G.Allan and M.Lannoo,Phys.Rev.B,62 (2000) 5109

[8] A.Thean and J.Leburton,IEEE Trans.Electron Dev.,22 (2001) 148

[9] G.Iannacone and P.Coli,Appl.Phys.Lett.,78 (2001) 2046

[10] J.Sée,P.Dollfus and S.Galdin,J.Appl.Phys.,92 (2002) 3141

[11] Y.F.Zhao et.al.,J.Vac.Sci.Technol.B,15 (1997) 1452

[12] X.Leyronas,M.Combescot,Solid State Commun.119 (2001) 631; 
cond-mat/0201553

[13] Z.A.Weinberg,J.Appl.Phys.53 (1982) 5052

[14] E.Cassan,S.Galdin,P.Dollfus and P.Hesto,J.Appl.Phys.86 (1999) 3804

[15] J.Sée,P.Dollfus,and S.Galdin,Phys.Rev.B 66 (2002) 193307

[16] S.Wolfram,The Mathematica Book,Cambridge University Press,(1999)