General scheme of the optical device analysis

Here we summarize the models and equations that is used for the optical device analysis in nextnano++.

Table of contents

Determination of carrier densities and current densities

Quantum mechanical calculation of charge carrier densities

Multi-band model (\(\mathbf{k}\cdot\mathbf{p}\) model)

Once the \(\mu\)-th component envelope function of the \(j\)-th eigenstate of electron (\(l=\text{c}\)) or hole (\(l=\text{v}\)) in the \(i\)-th band is obtained as \((F_\mu)^{i}_{l,j}(\mathbf{x})\) from the multi-band Schrödinger equation, the probability distribution of this \(j\)-th eigenstate reads

(2.2.1)\[p_{l,j}^{i}(\mathbf{x}) = \sum_{\mu}\bigg|(F_\mu)^{i}_{l,j}(\mathbf{x})\bigg|^2.\]

where we are assuming 3D structure so far.

Then the quantum mechanical carrier densities for 3D structure are defined from these probability densities, energy eigenvalues \(E_{\text{c},j}\) and \(E_{\text{v},j}\), position-dependent quasi-Fermi levels \(E_{\text{F},n}(\mathbf{x})\) and \(E_{\text{F},p}(\mathbf{x})\) as

(2.2.2)\[n(\mathbf{x}) = \sum_{i\in{\text{CB}}} g_\text{c}^{i}\sum_{j} p_{\text{c},j}^{i}(\mathbf{x})\ f\bigg([E_{\text{c},j}^i-E_{\text{F},n}(\mathbf{x})]/kT\bigg)\]
(2.2.3)\[p(\mathbf{x}) = \sum_{i\in{\text{VB}}} g_\text{v}^{i}\sum_{j} p_{\text{v},j}^{i}(\mathbf{x})\ f\bigg([-E_{\text{v},j}^i+E_{\text{F},n}(\mathbf{x})]/kT\bigg)\]

where \(f(E)\) is the Fermi-Dirac distribution at temperature \(T\), \(g_\text{c}^{i}\) and \(g_\text{v}^{i}\) represent the possible spin and valley degeneracies.

When the simulation is over 1D structure, the wave function can be separated into the plane wave specified with the lattice wave vector \(\mathbf{k}_\parallel\) in the lateral 2D direction and the quantized wave function in the growth direction, which has the \(\mathbf{k}_\parallel\)-dependency. Then the charge carrier densitiy is obtained by the following integral over \(\mathbf{k}_\parallel\):

(2.2.4)\[n(x) = \sum_{i\in{\text{CB}}} g_\text{c}^{i}\sum_{j}\frac{1}{(2\pi)^2}\int_{\Omega_{BZ}}\ d^2\mathbf{k}_\parallel\ p_{\text{c},j}^{i}(x,\mathbf{k}_\parallel)\ f\bigg([E_{\text{c},j}^i(\mathbf{k}_\parallel)-E_{\text{F},n}(x)]/kT\bigg)\]
(2.2.5)\[p(x) = \sum_{i\in{\text{VB}}} g_\text{v}^{i}\sum_{j}\frac{1}{(2\pi)^2}\int_{\Omega_{BZ}}\ d^2\mathbf{k}_\parallel\ p_{\text{v},j}^{i}(x,\mathbf{k}_\parallel)\ f\bigg([-E_{\text{v},j}^i(\mathbf{k}_\parallel)+E_{\text{F},n}(x)]/kT\bigg)\]

Here the integration is over the two-dimensional Brillouin zone \(\Omega_{BZ}\).

Similarly, the charge carrier densities for 2D structure is calculated by the integral over the 1-dimensional Brillouin zone as

(2.2.6)\[n(\mathbf{x}) = \sum_{i\in{\text{CB}}} g_\text{c}^{i}\sum_{j}\frac{1}{2\pi}\int_{\Omega_{BZ}}\ d \text{k}\ p_{\text{c},j}^{i}(\mathbf{x},\text{k})\ f\bigg([E_{\text{c},j}^i(k)-E_{\text{F},n}(\mathbf{x})]/kT\bigg)\]
(2.2.7)\[p(\mathbf{x}) = \sum_{i\in{\text{VB}}} g_\text{v}^{i}\sum_{j}\frac{1}{2\pi}\int_{\Omega_{BZ}}\ d \text{k}\ p_{\text{v},j}^{i}(\mathbf{x},\text{k})\ f\bigg([-E_{\text{v},j}^i(k)+E_{\text{F},p}(\mathbf{x})]/kT\bigg)\]

Single-band model

Things are simpler.

When the single-band Schrödinger equation is set to be solved, the envelope function of the \(j\)-th eigenstate has only one component \(F^{i}_{l,j}(\mathbf{x})\). Also, the k-integration in (2.2.4) to (2.2.7) can be done analytically due to the parabolic dispersion according to the effective mass tensor \(\underline{m}_e^{*i}\) and \(\underline{m}_h^{*i}\).

Thanks to this simpicity the quantum mechanical charge carrier densities for \(d\)-dimensional simulation can be written up by the following expression:

(2.2.8)\[n(\mathbf{x}) = \sum_{i\in{\text{CB}}} g_\text{c}^{i}\ \bigg( \frac{m_{\text{dos,e}}kT}{2\pi\hbar^2} \bigg)^{(3-d)/2}\ \sum_{j} p_{\text{c},j}^{i}(\mathbf{x})\ \mathcal{F}_{(1-d)/2}\bigg([E_{\text{c},j}^i-E_{\text{F},n}(\mathbf{x})]/kT\bigg)\]
(2.2.9)\[p(\mathbf{x}) = \sum_{i\in{\text{VB}}} g_\text{v}^{i}\ \bigg( \frac{m_{\text{dos,h}}kT}{2\pi\hbar^2} \bigg)^{(3-d)/2}\ \sum_{j} p_{\text{v},j}^{i}(\mathbf{x})\ \mathcal{F}_{(1-d)/2}\bigg([-E_{\text{v},j}^i+E_{\text{F},p}(\mathbf{x})]/kT\bigg)\]

TODO: The sign in the fermi-dirac integral might be opposite. check the source code.

Here \(\mathcal{F}_n(E)\) denotes the Fermi-Dirac integral of order \(n\) and \(m_{\text{dos,}\lambda}^i\) is so-called density-of-states mass defined as

(2.2.10)\[m_{\text{dos,}\lambda}^i=\big( \text{det}\ \bar{m}_\lambda^{*i} \big)\ \ \ \ \lambda=\text{e,h}\]

where \(\bar{m}_\lambda^{*i}\) describes the \(2\times 2\) or \(1\times 1\) submatrix of the effective mass tensor \(\underline{m}_\lambda^{*i}\) in the direction of \(\mathbf{k}_\parallel\).

In any cases, the carrier densities are dependent on the electrostatic potential \(\phi(\mathbf{x})\) through the wave function, which is obtained from the \(\phi\)-dependent Hamiltonian \(H(\phi)\). Thus we can also write them as \(n(\mathbf{x},\phi)\) and \(p(\mathbf{x},\phi)\), which enters into the non-linear Poisson equation introduced later.

Moreover, when the current equation is included in the calculation scheme, seeing the carrier densities as \(n(\mathbf{x}, \phi, E_{\text{F},n})\) and \(p(\mathbf{x}, \phi, E_{\text{F},p})\) makes it easy to understand what the self-consistent calculation is actually doing.

Classical calculation of charge carrier densities

Things are much more simpler.

When any kind of Schrödinger equation is not solved, the charge carrier densities are estimated from the position-dependent conduction and valence band edges \(E_\text{c}^i(\mathbf{x})\) and \(E_\text{v}^i(\mathbf{x})\), quasi-Fermi levels, and the electrostatic potential \(\phi(\mathbf{x})\) in the context of Thomas-Fermi approximation.

These classical charge carrier densities are calculated as

(2.2.11)\[n(\mathbf{x}) = \sum_{i\in{\text{CB}}} N_\text{c}^i(T)\ \mathcal{F}_{1/2}\bigg([-E_\text{c}^i(\mathbf{x})+e\phi(\mathbf{x})+E_{\text{F},n}(\mathbf{x})]/kT\bigg)\]
(2.2.12)\[p(\mathbf{x}) = \sum_{i\in{\text{VB}}} N_\text{v}^i(T)\ \mathcal{F}_{1/2}\bigg([E_\text{v}^i(\mathbf{x})-e\phi(\mathbf{x})-E_{\text{F},p}(\mathbf{x})]/kT\bigg).\]

Here \(N_\text{v}^i(T)\) and \(N_\text{v}^i(T)\) are the equivalent density of states at the conduction and valence band edges, which are given by

(2.2.13)\[ N_l^i(T) = g_{l}^i\bigg(\frac{m_{\text{dos,}\lambda}^ikT}{2\pi\hbar}\bigg)^{2/3}\ \ \ \ \ \ (l,\lambda)=(\text{v,h}),\ \text{or}\ (\text{c,e}).\]

Here \(m_{\text{dos,}\lambda}^i\) is the density-of-mass for \(d=3\) defined in (2.2.10).

This calculation of carrier densities is much faster than the quantm mechanical calculation, but the quantum effect such as energy quantization, carrier leackage into the barrier, etc. cannot be taken into account.

Also in this case, the carrier densities can be written as \(n(\mathbf{x},\phi)\) and \(p(\mathbf{x},\phi)\), which enters into the non-linear Poisson equation introduced next.

Moreover, when the current equation is included in the calculation scheme, seeing the carrier densities as \(n(\mathbf{x}, \phi, E_{\text{F},n})\) and \(p(\mathbf{x}, \phi, E_{\text{F},p})\) makes it easy to understand what the self-consistent calculation is actually doing.

Poisson equation

Poisson equation

This equation governs the relation between the electrostatic potential \(\phi(\mathbf{x})\) and total charge density distribution \(\rho(\mathbf{x},\phi)\) as follows:

(2.2.14)\[-\nabla\cdot[\varepsilon_0\varepsilon_r (\mathbf{x}) \nabla\cdot\phi(\mathbf{x})] = \rho(\mathbf{x},\phi)\]

where \(\varepsilon_0\) is the vacuum permittivity, \(\varepsilon_r\) is the material dependent static dielectric constant. And the total charge density distribution consists of the densities of ionized donors \(N_D^+\), ionized acceptors \(N_D^-\), piezoelectric and pyroelectric charge \(\rho_{pz}\) and \(\rho_{py}\), besides the carrier densities \(n(\mathbf{x},\phi)\) and \(p(\mathbf{x},\phi)\), which are calculated either classically or quantum mechanically:

(2.2.15)\[\rho(\mathbf{x},\phi)=e[-n(\mathbf{x},\phi)+p(\mathbf{x},\phi)+N_D^+(x)-N_A^-(\mathbf{x})+\rho_\text{pz}(\mathbf{x})+\rho_\text{py}(\mathbf{x})]\]

When the Schrödinger-Poisson equation is solved, i.e. quantum_poisson{ } is specified in run{ } section, the carrier densities defined in either multi-band model or single-band model are substituted into this \(\rho(\mathbf{x},\phi)\) and the Poisson equation is solved accordingly. Then the resulting \(\phi(\mathbf{x})\) is returned into the Schrödinger equation and the carrier densities are calculated once again.

This cycle is continued until the carrier densities satisfies the convergence criteria, which can be tuned by the users from run{ poisson{ } }. The final result of \(n(\mathbf{x},\phi)\), \(p(\mathbf{x},\phi)\) and \(\phi(\mathbf{x})\) must satisfy both Schrödinger and Poisson equations, or we can say that the Schrödinger equation and Poisson equation are self-consistent with respect to the resulting carrier densities and electrostatic potential.

On the other hand, when only the Poisson equation is solved, i.e. only poisson{ } is specified run{ } section, the carrier densities are calculated according to (2.2.11) and (2.2.12) instead. We can say in other words that the carrier density calculation in the context of Thomas-Fermi approximation and the Poisson equation are self-consistent with respect to the resulting carrier densities and electrostatic potential.

Ionized donor/acceptor densities

The densities of ionized impurities are calculated in the context of Thomas-Fermi approximation with these formulas:

(2.2.16)\[N_\text{D}^+(\mathbf{x})=\sum_{i\in \text{Donors}}\frac{N_{\text{D},i}(\mathbf{x})}{1+g_{\text{D},i}\exp((E_{\text{F},n}(\mathbf{x})-E_{\text{D},i}(\mathbf{x}))/k_\text{B}T)}\]
(2.2.17)\[N_\text{A}^-(\mathbf{x})=\sum_{i\in \text{Acceptors}}\frac{N_{\text{A},i}(\mathbf{x})}{1+g_{\text{A},i}\exp((E_{\text{A},i}(\mathbf{x})-E_{\text{F},p}(\mathbf{x}))/k_\text{B}T)}\]

where the summation is over all different donor or acceptors, \(N_\text{D},N_\text{A}\) are the doping concentrations, \(g_\text{D},g_\text{A}\) are the degeneracy factors (\(g_\text{D}=2\) and \(g_\text{A}=4\) for shallow impurities), and \(E_{D},E_{A}\) are the energies of the neutral donor and acceptor impurities, respectively.

These energies of neutral impurities \(E_{\text{D},i},E_{\text{A},i}\) are determined by the ionization energies \(E_{\text{D},i}^\text{ion},E_{\text{A},i}^\text{ion}\) , the bulk conduction and valence band edges (including shifts due to strain) and the electrostatic potential.

(2.2.18)\[E_{\text{D},i}(\mathbf{x})=E_\text{c}(\mathbf{x})-e\phi(\mathbf{x})-E_{\text{D},i}^\text{ion}(\mathbf{x})\]
(2.2.19)\[E_{\text{A},i}(\mathbf{x})=E_\text{v}(\mathbf{x})-e\phi(\mathbf{x})+E_{\text{A}}^\text{ion}(\mathbf{x})\]

Piezoelectric and pyroelectric charge densities

\(\rho_{pz}\) and \(\rho_{py}\) are calculated according to the result of strain equation. (TO be updated)

Current equation

Current equation

The continuity equations in the presence of creation (generation, \(G\) ) or annihilation (recombination, \(R\) ) of electron-hole pairs read

(2.2.20)\[\begin{split}\begin{aligned} -e\frac{\partial n}{\partial t} + \nabla\cdot \big(-e\mathbf{j}_n(\mathbf{x})\big) &= -e\big(G(\mathbf{x})-R(\mathbf{x})\big),\\ e\frac{\partial p}{\partial t} + \nabla\cdot e\mathbf{j}_p(\mathbf{x}) &= e\big(G(\mathbf{x})- R(\mathbf{x})\big), \end{aligned}\end{split}\]

where the current is proportional to the gradient of quasi Fermi levels \(E_{\text{F},n/p}(\mathbf{x})\)

(2.2.21)\[\begin{split}\begin{aligned} \mathbf{j}_n(\mathbf{x}) &= -\mu_n(\mathbf{x})n(\mathbf{x})\nabla E_{\text{F},n}(\mathbf{x}),\\ \mathbf{j}_p(\mathbf{x}) &= \mu_p(\mathbf{x})p(\mathbf{x})\nabla E_{\text{F},p}(\mathbf{x}). \end{aligned}\end{split}\]

Here the charge current has the unit of (area)\(^{-1}\)(time)\(^{-1}\). \(\mu_{n/p}\) are the mobilities of each carrier. In nextnano++, \(\mu_{n/p}\) are determined using the mobility model specified in the input file under currents{ }.

Hereafter we consider stationary solutions and set \(\dot{n}=\dot{p}=0\). The governing equations then reduce to

(2.2.22)\[\begin{split}\begin{aligned} \nabla\cdot\mu_n(\mathbf{x})n(\mathbf{x})\nabla E_{\text{F},n}(\mathbf{x})&=-(G(\mathbf{x})-R(\mathbf{x})),\\ \nabla\cdot\mu_p(\mathbf{x})p(\mathbf{x})\nabla E_{\text{F},p}(\mathbf{x})&=G(\mathbf{x})-R(\mathbf{x}), \end{aligned}\end{split}\]

which we call current equation.

We can also say that the current equation governs the relationship between the carrier densities \(n(\mathbf{x})\), \(p(\mathbf{x})\) and quasi Fermi levels \(E_{\text{F},n/p}(\mathbf{x})\).

The nextnano++ tool solves this equation and Poisson equation (and also Schrödinger equation) self-consistently.

In their solution, the corresponding calculation of the carrier densities \(\big(n(\mathbf{x}, \phi, E_{\text{F},n}),p(\mathbf{x}, \phi, E_{\text{F},p})\big)\) and Poisson equation are firstly iterated for a given quasi-Fermi levels until the carreir densities converge. Then the resulting carrier densities are substituted into the current equation and the quasi-Fermi levels are updated. This whole cycle is iterated until the quasi-Fermi levels satifies the convergence criteria, which can be tuned by the users from run{ current_poisson{ } } or run{ quantum_current_poisson{ } }.

Recombination/Generation

The recombination mechanisms that nextnano++ takes into account for the right-hand-side of (2.2.20) are

  • Shockley-Read-Hall (SRH) recombination

  • Auger recombination

  • Radiative recombination

  • “fixed (applied)”

The equations and parameters used for the three recombination mechanisms on the top are explained here: recombination_model{ }.

The last one “fixed (applied)” is the contribution defined from structure{region{generation{}}} and optics{ photogeneration{ } }. These typically represent generation instead of recombination and used for the simulation of the devices under irradiation such as solar cells or CCDs. (For example, see nextnano++ tutorial GaAs solar cell.)

Optoelectronic characteristics based on the semi-classical model

According to the specification in the section classical{ }, nextnano++ can calculate optoelectronic characteristics of the arbitrary structure by means of the so-called semi-classical model.

In this model, various quantities are calculated from the spontaneous emission rate, which is calculated at each position \(\mathbf{x}\) for the photons with each energy \(E\) based on the energy-resolved carrier densities \(n(\mathbf{x},E)\) and \(p(\mathbf{x},E)\) obtained in the forgoing simulation.

  • Spontaneous emission rate

(2.2.23)\[R^{spon}_{\mathrm{rad}}(\mathbf{x}, E)=C(\mathbf{x})\int dE_\text{h}\int dE_\text{e}\ n(\mathbf{x},E_\text{e})p(\mathbf{x},E_\text{h})\delta(E_\text{e}-E_\text{h}-E).\]

Here \(C(x)\) [\(\mathrm{cm}^3\mathrm{s}^{-1}\)] is the (material-dependent) radiative recombination parameter which is proportional to the one specified in the database (Radiative recombination)

Then the other optical characteristics like stimulated emission rate, absorption/gain spectrum, and the imaginary part of the dielectric constant are calculated according to this \(R^{spon}_{\mathrm{rad}}(\mathbf{x}, E)\).

  • Stimulated emission rate

Stimulated emission rate is calculated here as the net emission rate containing both the generation by the stimulated absorption and the recombination by the spontaneous and stimulated emission according to the following equation:

(2.2.24)\[R^\text{stim}_{rad,net}(\mathbf{x},E)=\left(1-e^\frac{{E-(E_{\text{F}n}-E_{\text{F}p})}}{k_\text{B}T} \right)R_\text{rad}^{spon}(\mathbf{x},E)\]

The reference equation is eq.(9.2.39) of [ChuangOpto1995].

\(R^{spon}_{\mathrm{rad}}(E)\) and \(R^{stim}_\text{rad,net}(E)\) are output on Optical/semiclassical_spectra_photons_~.dat and stim_emission_photons_~.dat as the integral of the above two quantities over \(\mathbf{x}\), i.e.

(2.2.25)\[R^{spon}_\text{rad}(E)=\int d\mathbf{x}\ R^\text{spon}_{rad}(\mathbf{x},E),\ \ \ \ \ R^{stim}_\text{rad,net}(E)=\int d\mathbf{x}\ R^{stim}_\text{rad,net}(\mathbf{x},E)\]

On the other hand, \(R^{stim}_\text{rad,net}(\mathbf{x})\) is obtained as the integral of \(R^{stim}_\text{rad,net}(\mathbf{x},E)\) over the photon energy, which is written as “radiative” in the output file recombination.dat, i.e.

(2.2.26)\[R^{stim}_\text{rad,net}(\mathbf{x})=\int dE\ R^{stim}_\text{rad,net}(\mathbf{x},E).\]

Note

Precisely speaking, \(R^{stim}_\text{rad,net}(\mathbf{x})\) is not directly integrated from \(R^{stim}_\text{rad,net}(\mathbf{x},E)\) but calculated according to the equation

\[R^{stim}_\text{rad,net}(\mathbf{x})=C(x)n(x)p(x)\left(1-e^\frac{E_{\text{F}p}-E_{\text{F}n}}{k_\text{B}T} \right).\]

This is meanwhile equivalent to (2.2.26).

  • Generation by the irradiation (fixed(applied))

There is another radiative recombination rate output on recombination.dat called “fixed(applied)”, which should be always negative. This is the contribution of the generation specified from structure{region{generation{}}} and optics{ photogeneration{ } }. When we do not specify either of them, this recombination rate is always 0.

(2.2.27)\[\begin{split}\begin{aligned} R_\text{fixed}(\mathbf{x})=&-\big(G(\mathbf{x})&\text{ specified from structure{ }}\big)\\ &-\big(\int dE\ G(E,\mathbf{x})\text{ calculated according to the configuration in classical{ }}\big). \end{aligned}\end{split}\]

This is mostly used for the analysis of the absorbing devices such as solar cells or CCDs.

  • photocurrent

Then the photocurrent \(I_\text{photo}\) is calculated as the summation of the integration of these “radiative” and “fixed”:

(2.2.28)\[I_\text{photo} =e\cdot\bigg( \int d\mathbf{x}\ R^{stim}_\text{rad,net}(\mathbf{x}) + \int d\mathbf{x}\ R_\text{fixed}(\mathbf{x}) \bigg)\]

  • internal quantum efficiency

is calculated as

(2.2.29)\[\eta_{IQE} = \frac{I_\text{photo}}{I_\text{total}}\]

where \(I_\text{total}\) is the total injected current consisted of both electron and hole currents.

  • volume quantum efficiency

, which is also called as radiative quantum efficiency, is calculated as

(2.2.30)\[\eta_{VQE} = \frac{R^{stim}_\text{rad,net} + R_\text{fixed}}{R_\text{total}}\]

where \(R_\text{total}=R^{stim}_\text{rad,net} + R_\text{fixed} + R_\text{Auger} + R_\text{SRH}\) is the total recombination rate including both radiative and non-radiative recombination.

Both \(\eta_{IQE}\) and \(\eta_{VQE}\) agree if the electrons and holes injected into the active region are fully consumed up by the recombination there. However, if they are not consumed up, \(e\cdot R_\text{total}<I_\text{charge}\) and this results in \(\eta_\text{IQE1}>\eta_\text{IQE2}\)

Note

If you have any comments on the terminologies and definitions of these quantities, please send to support [at] nextnano.com..

Moreover, the electrical power and optical power are calculated and output in power.dat:

  • Power

    (2.2.31)\[\sum_{i}V_\text{i-th contact}\cdot I_\text{i-th contact}\]

  • Absorbed-power

    (2.2.32)\[\int dEd\mathbf{x}\ E\cdot G(E,\mathbf{x})\]

    where \(G(E,x)\) is the generation rate calculated according to the configuration in classical{ }.

  • Emitted-power

    (2.2.33)\[\int dEdx\ E\cdot R^{spon}_\text{rad}(E,x)\]

Optoelectronic characteristics based on the quantum model

The nextnano++ tool has another important calculation scheme of optical properties, which is specified in the section optics{ }. Here nextnano++ calculates them using the Fermi’s golden rule (time-dependent perturbation theory) with 8-band k.p model.

These quantities are now supported

  • Optical absorption coefficient

  • Real/imaginary part of the dielectric constant

  • Refractive index

  • Optical gain as a negative part of optical absorption coefficient

  • Spontaneous emission rate

  • Transition intensity (optical matrix element)

For further detail about this section, please see Optical absorption for interband and intersubband transitions.

References

  • [ZiboldPhD2007]
    Semiconductor based quantum information devices: Theory and simulations
    T. Zibold
    Selected Topics of Semiconductor Physics and Technology (G. Abstreiter, M.-C. Amann, M. Stutzmann, and P. Vogl, eds.), Vol. 87, Verein zur Förderung des Walter Schottky Instituts der Technischen Universität München e.V., München, 151 pp. (2007)
  • [BirnerPhD2011]
    Modeling of semiconductor nanostructures and semiconductor-electrolyte interfaces
    S. Birner
    Selected Topics of Semiconductor Physics and Technology (G. Abstreiter, M.-C. Amann, M. Stutzmann, and P. Vogl, eds.), Vol. 135, Verein zur Förderung des Walter Schottky Instituts der Technischen Universität München e.V., München, 239 pp. (2011)
    ISBN 978-3-941650-35-0
  • [ChuangOpto1995]
    Physics of Optoelectronic Devices
    S. L. Chuang
    John Wiley & Sons, Inc., New York (1995)

Last update: nn/nn/nnnn