Electron-phonon relaxation and excited electron distribution in gallium nitride
We develop a theory of energy relaxation in semiconductors and insulators highly excited by the long-acting external irradiation. We derive the equation for the non-equilibrium distribution function of excited electrons. The solution for this function breaks up into the sum of two contributions. The low-energy contribution is concentrated in a narrow range near the bottom of the conduction band. It has the typical form of a Fermi distribution with an effective temperature and chemical potential. The effective temperature and chemical potential in this low-energy term are determined by the intensity of carriers’ generation, the speed of electron-phonon relaxation, rates of inter-band recombination and electron capture on the defects. In addition, there is a substantial high-energy correction. This high-energy ’tail’ covers largely the conduction band. The shape of the high-energy ’tail’ strongly depends on the rate of electron-phonon relaxation but does not depend on the rates of recombination and trapping. We apply the theory to the calculation of a non-equilibrium distribution of electrons in irradiated GaN. Probabilities of optical excitations from the valence to conduction band and electron-phonon coupling probabilities in GaN were calculated by the density functional perturbation theory. Our calculation of both parts of distribution function in gallium nitride shows that when the speed of electron-phonon scattering is comparable with the rate of recombination and trapping then the contribution of the non-Fermi ’tail’ is comparable with that of the low-energy Fermi-like component. So the high-energy contribution can affect essentially the charge transport in the irradiated and highly doped semiconductors.
The study of relaxation processes in crystals is indispensable for the understanding of a number of physical phenomena induced by external sources at high excitation levels. First of all this concerns to the experiments with the impact by electron beams and powerful lasers [1, 2]. The knowledge of the energy and momentum relaxation is also indispensable for study of the high-field transport of charge carriers [3, 4, 5], including the case of electric breakdown [6, 7]. One should also mention such phenomena as laser ablation , abrupt thermal impact [9, 10], ultra-fast phase transitions , photocatalysis  and the solar energy conversion [13, 14]. In the investigation of the energy relaxation processes in [1, 15], [14, 16, 17, 18], , [19, 20] and  the non-equilibrium distribution of electrons usually is simulated as a quasi-Fermi function which is localized close to the bottom the conduction band.
In our works [12, 22, 23, 24] on the example of the compounds and it was established that in a wide-bandgap material the evolution of energy relaxation has some peculiarities. It has been shown that the form of quasi-equilibrium distribution differs markedly from the case of narrow-bandgap compounds, and covers a wide range of energy within the conduction and valence bands over band edges at a distance comparable to the width of the band gap.
In the current article we will consider the quasi-equilibrium distribution in more general terms that allow us to incorporate into the theory  the processes of radiative recombination and trapping of carriers. This makes it possible to specify the energy dependence of the distribution function either directly near the band edges as well as away from them. We apply the proposed theory to the case of gallium nitride in the wurtzite structure.
The paper is organized as follows. In Sec.2, we present the analytical methods we use. First, we derive the equation for the quasi-equilibrium distribution function. Next, its formal solution is given in Subsec.2.1. On this basis in Subsec.2.2 we develop an effective phonon model that allows an exact solution which is presented in Subsec.2.3. Section 3 describes the technics of numerical calculation for . The necessary material dependent quantities, namely, optical excitation probabilities and electron-phonon coupling characteristics are calculated ab initio within the pseudopotential approach and density functional theory. Section 4 presents the application of Sec.2 and Sec.3 to the non-equilibrium distribution of carriers in the conducting band of GaN. Finally, the results are discussed and summarized in Sec.5.
2 Calculation methods
The kinetic equation for the time-dependent distribution function of a non-equilibrium electron state created by the powerful external source of high-energy irradiation was investigated in Ref.. Here is the density of electron states in the conduction band, is the occupation number of non-equilibrium electrons, the energy is taken with respect of the bottom of conduction band (hereafter ’excess energy’). The excitation of electrons from the valence to conduction band by external source is followed by the fast momentum randomization, within several femtoseconds . That is why the occupation number of electron’s band state can be regarded just as a function of energy , the phonon occupation number is considered to be energy-dependent as well: .
We restrict ourselves to the case when the maximum excess energy of external excitation is located lower than the threshold of impact ionization . Numerically it roughly equals to the band-gap value above the bottom of the conduction band. Once the energy of an electron is lower than the production of secondary electrons and holes is prohibited since it is impossible to simultaneously satisfy the energy and momentum conservation laws. Therefore inelastic electron-electron processes are absent in this range, and energy relaxation can only be associated with electron-phonon scattering. For this threshold corresponds to the excess energy equal to 3.25 eV.
Retaining terms linear and quadratic by phonon energy which is a small quantity relative to the band energy one comes to the following equation for a non-equilibrium distribution of electrons in the conduction band :
A spectral function of electron-phonon interaction is defined as:
is the matrix element of electron-phonon interaction operator. is the energy of a phonon of branch with the wave vector . The term describes the distribution of excited carriers created by external sources of generation. The expression takes into consideration the flow of electrons from the conduction band to the valence bands and to the impurity levels. The range of integration in Eqs.(2) extends from zero up to the maximal phonon energy . The function
specifies the probability per unit time for an electron to leave a given excess energy level , so it is the electron-phonon relaxation rate.
The well-known Eliashberg spectral function of the electron-phonon interaction, which relates to the Fermi-level in superconductors , can be generalized for an arbitrary excess energy level in semiconductor as
This function, in its turn, determines in semiconductor the electron-phonon coupling constant , relating to a given energy :
2.1 Stationary solution for the distribution function
If the action of external excitation continues for a time long enough then the quasi-stationary distribution become established. Then a search of a quasi-stationary distribution function reduces itself to the solution of the equation for the time-independent occupation number . It is convenient to seek for a stationary solution of (1) in the form . Here is the solution for the stationary equation (1) but without sources and sinks. One can check by a direct substitution that there exists a partial solution for in the following form
Here is an integration constant.
Under conditions of low intensity of external exposure one has to regard as a small quantity. In this case keeping the leading terms we come from Eq.(1) to the linearized equation with a stationary solution
Both distributions (7) and (9) should be localized in the energy region where the electron-phonon relaxation is inefficient, i.e. at , where is the maximal energy of the phonon spectrum. An equation for the correction comes out by substitution of into Eq.(1). We restrict ourselves by the examination of excitation regimes when can be treated as a small and smooth quantity. Retaining the leading terms the equation for the correction term in the stationary regime acquires the following form:
The solution of Eq.(10) which satisfies the condition = 0 at can be written as
Here is the maximal level of electrons excitation.
In a stationary regime the total number of electrons generated per unit time by the external source should be equal to the number of electrons that leave the conduction band during the same time due to the inter-band recombination and trapping by the impurity levels:
As is mostly localized in a small region of order of phonon spectrum width above the conduction band bottom then one can neglect the value at .
The expression similar to Eq.(13) for the high energy ’tail’ of the distribution function was formerly discussed in Refs.[22, 23, 12] where the electron-hole recombination has not been formally taken into consideration. One can see now that the electron-hole recombination actually does not affect the shape of the high-energy ’tail’ of distribution.
2.2 The ’effective phonon’ model
Let us define an averaged energy loss per one electron in the process of a single transition act from an excess level as
In calculations of Refs.[22, 23] it was shown that for realistic band structures in and this quantity manifests a weak energy dependence. Below we will show that this is also valid for . Hence we neglect this dependence and introduce the Einstein-like lattice vibrational model considering the electron scattering by a single ’effective’ phonon with energy , and identify with . In this model we assume that
Here we introduce the averaged electron-phonon coupling factor
Then the definition (4) implies that
Expressions for , , follow from definitions (2):
Taking these approximations into account we find also the expression for as defined by the Eq.(8)
Then Eq.(7) for acquires the form of quasi-Fermi distribution
where effective temperature of the excited electrons near the bottom of the conduction band is . The effective chemical potential is defined implicitly by equation (12).
2.3 The case of weak degeneracy
We discuss hereafter the case of week degeneracy, when . This is a realistic case for the exposure by the sunlight or a mercury lamp, when the typical number of photons per unit cell of a crystal does not exceed . In this case Eq.(9) reduces to the quasi-Boltzmann distribution:
For further simplification let us introduce the separable approximation for the emission term in Eq.(13)
Here the function , in general time-dependent, is defined to be equal to the total number of excited electrons per unit of time, that depends on the intensity of exposure. In the quasi-stationary regime changes slowly, , and can be considered as a constant . The function describes the distribution of the excited electrons over the energy scale of the conduction band. It has to be normalized to unity: . Employing this definition and effective phonon approximation one can write the expression for the ’tail’ distribution as
The first-principle calculations demonstrate that in near the bottom of the conduction band has almost linear energy dependence. We approximate this dependence as .
In general, is also a non-equilibrium distribution of phonons, and it would be necessary to investigate its evolution in conjunction with the evolution of the electron distribution. However, if the duration of action of external excitation exceeds the characteristic relaxation time of the phonon subsystem, it can be regarded as Boltzmann function.
In the low-temperature limit when we can neglect a phonon occupation number. We also neglect the energy dependence of the recombination and trapping rate in the small energy range near the bottom of the conduction band, hence . In the high excitation regime the value and in Eq. (12) it can be replaced by infinity. This helps us, employing Eq.(12), to find the normalization coefficient in expression (21) for :
Now we see that the depends on the intensity of exposure, via the function , as well as on the rate of trapping and recombination processes.
Hence the total electron distribution function
is written as
3 Technical details
We apply the described approach to the gallium nitride compound in the wurtzite structure. Our calculations have been performed by employing the ab initio pseudo-potential computer code Quantum Espresso (QE) based on the density functional theory for the electron band structure and density-functional perturbation theory for the phonon band states . The pseudo-wave functions were expanded in plane waves with energy cutoff 820 eV. In the calculations of the electron states we employed a set of 50 wave vectors per irreducible part of the Brillouin zone (IPBZ). The calculations of the phonon states and characteristics of the electron-phonon coupling were performed for a set of 12 wave vector in the IPBZ. Gallium norm-conserving atomic pseudo-potential was calculated using the Bachelet-Hamann-Schlüter method  with Perdew-Zunger exchange-correlation potential . The pseudo-potential for nitrogen was constructed within the approach of Von Barth and Car . We will show later that such a way of calculations provide very good results both for electron and phonon band structures.
We compute the distribution function also basing on the band structure calculations by means of the QE computer code. Namely, if the energy of the quantum of optical excitation is then for the excess energy we sum the probabilities of all direct excitations from the electronic states at the energy to the states at the energy . Hence, the non-normalized function is
where is the probability of the transition between the states and . For the calculations of the transition probabilities we apply dipole approximation. The details of such an approach can be found in Ref..
4 Results and discussions
The calculated dispersion curves and the corresponding total density of electron states are given for in Fig.1. The direct band gap in -point is equal to 3.25 eV, in good agreement with experimental value of 3.4 eV (Ref.). The top of the valence band is formed by -like states of nitrogen, the bottom of the conduction band corresponds to the -like states of . Also the energy differences between the conduction band states at symmetry points of the Brilloiune zone agree with the results of previous calculations . E.g., our calculated , , differences are equal to 1.8, 1.55, 2.06 eV, whereas Bulutay et al. obtain 1.6, 1.75, 1.87 eV, respectively. The conduction band at displays the non-parabolic dispersion, so the density of states versus the excess electron energy deviates from the free-electron-like law practically immediately above the bottom of the conduction band. The non-parabolicity of the conduction band states has been revealed experimentally in Ref.. Up to the excess energy of about 1.25 eV, that is up to the energy in the central valley lying between the side and states, the total DOS displays almost linear dependence (see for inset in Fig.1), so near the bottom of the conduction band we approximate it with where =0.025.
In Fig.2 our calculated phonon dispersion curves are compared with the experimental data of Ref.. Similarly to the theoretical results of this paper, our computed frequencies are systematically about 5 % higher than the experimental ones. We demonstrate in this figure the calculated data scaled with the factor 0.95. In Fig.2 we also show the corresponding density of phonon states. In the following calculations the highest energy of the phonon spectrum is taken to be equal to 785 = 0.097 eV.
Hereafter we discuss the electron-phonon relaxation at excess energies which do not exceed the impact ionization level 3.25 eV in . The inter-valley electron-phonon scattering events happen at higher energies and are not a subject of the current study.
In Fig.3 we show the Eliashberg function calculated for three different values of excess energy in interval from 0.03 to 1.83 eV. We see that with increasing excess energy the electron-phonon scattering becomes stronger. For the states near the bottom of the conduction band the scattering of excess electrons energy occurs mainly via emission of optical phonons, both transversal () and longitudinal () branches, being in correspondence with the common opinion on the main role of such phonons in the processes of electron energy relaxation . We see, however, that with rising electron excess energy the contribution of acoustic branches essentially increases. In contrast to Ref. where the dominating scattering via emission of the phonons was supposed to occur, we find that at any value of the excess energy the contribution of the -branches is almost equal to that of longitudinal ones.
In Fig.4 we demonstrate the energy dependence of the , , values. The dependencies of and also demonstrate the strengthening of electron-phonon scattering with the rise of the excess energy. However, the effective phonon energy only slightly changes, so in the estimation of the relaxation rate and distribution curves we use the energy-average value = 0.085 eV. We see in the Eliashberg curves that the value lies just between the bands of the and branches which at any value of the excess energy have almost equal contributions to the electron relaxation rate.
According to Eq.(17) the change of average lifetime of a single electron with energy is determined by the changes of the density of states and the averaged probability of phonon emission . In Fig.5 we show , and versus the excess energy.
We see very rapid change of at the energy less than 0.25 eV. But at a higher energy rather slowly decreases with energy while rapidly increases that leads to the reduction of . In other words, the reduction of the electron relaxation time is determined by the expansion of the space of electronic states available for the electron that looses its energy via phonon emission. Note that the changes of in Fig.5 are very similar to that observed in our previous study of the electron dynamics.
The available literature data on the electron relaxation time in are contradictory. In Refs.[37, 38] the electron-phonon relaxation time has been evaluated using Fröhlich theory of electron-phonon scattering which employs the experimental values of dielectric constants and phonon frequencies. This estimation of relaxation time, 10 fs, corresponds well to our energy-averaged data. However, the energy dependence of the relaxation time has not been investigated.
The energy dependence of a relaxation time was investigated in Refs.[33, 39]. Bulutay et al , also relying on the Fröhlich’s electron-phonon matrix elements, performed the calculations of phonon-assisted scattering rate in taking into account the coupling of excited electrons with long-wave and phonon modes. According to these calculations, the relaxation rate increases with the excess energy, although this growth is not as rapid as in our first-principle calculations. Dasgupta et al. in their experimental work came to the conclusion that the relaxation rate decreases with excess energy. However, their estimation of the relaxation rate and time was not direct since they were based on the measurements of the inelastic mean free path (IMFP) of electrons in a hot electron transistor with as the base layer. The increase of the IMFP with the excess energy can be associated both with the increase of the electron velocity and the reduction of the electron-defect scattering. The discrepancy between our ab initio calculated data and experimental data of Ref. deserves a special study.
Let us note in this connection that Fröhlich’s interaction is limited to the interaction with the long-wavelength optical phonon. So the phenomenological electron-phonon Fröhlich’s Hamiltonian involves only comparatively small values of the wave vector. That could be enough in the case of a weak excitation because the scattering events are concentrated at the bottom of the conduction band. In our study of large deviation from equilibrium, we should not be limited to a single long-wavelength optical phonon, instead we should take into account all phonon states, including all nine optical branches over the entire Brillouin zone.
In Fig. 6 the results of our calculations for three distribution curves, , and are shown. Our aim is to compare values of and , so we omit the factor . Two parameters in our evaluations are determined by the experimental conditions, namely the excitation energy and the rate of electron-hole recombination and trapping . We analyse the variation of the distribution functions with respect to by calculating the and functions for three values of : 3.48, 4.08 and 5.08 eV, that correspond to the lower edge, the center and the top of our energy interval of interest. The value of is a ’technological’ parameter that strongly depends on the concentration of impurities and intrinsic defects in a crystal. We study its influence on the distribution functions by varying it from 0.04 eV, the value typical for the rate of electron-phonon scattering, up to the values of 10 and 100 times less, that is going to the pico-second times of recombination and trapping processes.
We see that the function has a simple gaussian shape with a peak energy increasing almost linearly with a rise of , . Respectively, with increasing the upper boundary of the distribution shifts to a higher energy, but the values of the function do not change markedly, keeping at the level about 40 . The functions and of appear to be similar to those calculated before for , in spite of a different way of calculating, based on the LMTO approach. We also see that when is taken to be close to the rate of electron-phonon scattering , then the values of and appear to be comparable. When is 10 times smaller than , then is also smaller than , but still not negligible. When is of two orders lower than , then the value of becomes relatively small, it does not exceeds 1.5% of .
On the basis on the kinetic theory we derived the equation for the quasi-equilibrium distribution function of electrons in semiconductors excited by strong long-lasting external irradiation. The leading role in the formation of the carriers’ distribution is associated with inelastic scattering of electrons by phonons. The solution essentially describe a stable flow of electrons downward in an energy scale. In the case of wide-bandgap material we can allocate two parts. One of them describes the accumulation of electrons above the bottom of the conduction band and has a conventional Fermi-like shape with a width of order of phonon spectrum. Its size depends both on the level of excitation and on the rate of ’leakage’ of electrons through the process of electron-hole recombination and electron trapping by the defects. The effective temperature and chemical potential generally depend on an occupation number of excited optical phonon.
When the rate of recombination of electrons and holes and capture electrons on defects is much smaller than the characteristic time of electron-phonon relaxation then a kind of ’bottleneck’ arises in that flow. The high-energy ’tail’ appears as a result of a balance between the supply from an external source and leakage. This second contribution is much smaller in size but extends much further over the conduction band.The form of the ’tail’ depends on the band structure and on the rate of electron-phonon energy relaxation. Interestingly, its shape and scale do not depend on the rate of recombination and trapping. The depth of their penetration into the conduction band depends on the rate of excitation of electrons from the valence band to the conduction one and can be comparable with the magnitude of the band gap above the bottom of the conduction band.
Based on the theory we have carried out the first-principle calculations for the characteristics of electron-phonon scattering and the distribution functions for in the wurtzite structure. A good quality of the calculated electron and phonon band structures is confirmed by the comparison with available experimental and theoretical data. By calculating the spectral distribution of the electron-phonon scattering (Eliashberg function) we have revealed that the contribution of the LO-phonon branch to the electron relaxation rate is comparable with that of the TO-phonon branch; the contribution of the acoustic branch, although markedly smaller, is not negligible. We also find that the electron-phonon relaxation time decreases with the rise of the electron excess energy, which is determined by the increase of the electron density of states.
Calculating the shape of electron distributions we find that when the rate of electron-phonon scattering is comparable with the rate of recombination and electron trapping by lattice defects then the values of the low-energy component are comparable with those of the high-energy ’tail’ of distribution. However, when the rate of recombination decreases then the value of the low-energy term becomes bigger. When the rate of electron recombination/trapping is two orders smaller, i.e. these processes occur in the picosecond energy range, the tail distribution does not exceed 1.5 % of the low-energy distribution.
Considering as a typical example we conclude that the high-energy correction has to be taken into account in the study of the behavior of highly excited charge carrier in a wide-bandgap semiconductor compounds owing a high rate of electron-phonon scattering. The length of the high-energy ’tail’ may cover a substantial part of the conduction band. Its contribution to the total concentration of non-equilibrium carriers may be comparable to that of the conventional Fermi part and thus can have a significant influence on the transport properties of the irradiated crystal.
-  A.Othonos. J.Appl.Phys, 83:1789, 1998.
-  J.Shah. Ultrafast Spectroscopy of Semiconductors and Semiconductor Nanostructures. Springer,New York, 1999.
-  E.Bringuier. Phys.Rev.B, 57:2280, 1998.
-  W.Tisdale, K.Williams, B.Timp, D.Norris, E.Aydil, and X.Y. Zhu. Science, 328:1543, 2010.
-  M.Dur, S.Goodnick, S.Pennathur, J.Wager, M.Reigrotzki, and R.Redmer. J.Appl.Phys., 83:3176, 1998.
-  Y.Sun, S.Boggs, and R.Ramprasad. Appl.Phys.Lett., 101:132906, 2012.
-  F.Bertazzi, M.Moresco, and E.Bellotti. J.Appl.Phys., 106:063718, 2009.
-  C.Wang, L.Jiang, F.Wang, , X.Li, Y.Yuan, and H.Tsai. Phys.Lett.A, 375:3200, 2011.
-  A.A.Grinberg and S.Luryi. Phys.Rev.Lett., 65:1251, 1990.
-  Y.Gurevich, G.Logvinov, Y.Drogobitskiy, and A.Carballo-Sanchez J. Salazar. J.Appl.Phys., 91:183, 2002.
-  D.von der Linde, K.Sokolowski-Tinten, and J.Bialkowski. Appl.Surf.Sci., 109-110:1, 1997.
-  V.P.Zhukov, V.G.Tyuterev, E.V.Chulkov, and P.M.Echenique. International Journal of Photoenergy, 2014:738921, 2014.
-  M.Neges, K.Schwarzburg, and F.Willig. Solar Energy Materials and Solar Cells, 90:2107, 2006.
-  E.Tea, H.Hamzeh, and F.Aniel. J.Appl.Phys., 110:113108, 2011.
-  M.Harb, R.Ernstorfer, T.Dartigalongue, C.Hebeisen, R.Jordan, and R.Miller. J.Phys.Chem.B, 110:25308, 2006.
-  J.Callan, A.-T.Kim, L.Huang, and E.Mazur. Chem.Phys., 251:167, 2000.
-  G.Fasol, W.Hackenberg, H.P.Hughes, K.Ploog, E.Bauser, and H.Kano. Phys.Rev.B, 41:1461, 1990.
-  P.Supancic, U.Hohenester, P.Kocevar, D.Snoke, R.Hannak, and W.Ruhle. Phys.Rev., 53:7785, 1996.
-  C.Collier, B.Born, X.Jin, and J.Holzman. Appl.Phys.Lett., 103:072106, 2013.
-  U.Hohenester, P.Supancic, P.Kocevar, X.Q.Zhou, W.Kutt, and H.Kurz. Phys.Rev.B, 47:13233, 1993.
-  S.S.Prabhu, A.S.Vengurlekar, S.K.Roy, and J.Shah. Phys.Rev.B, 51:14233, 1995.
-  V.Zhukov and V.Tyuterev. Phys.Sol.State, 54:2173, 2012.
-  V.Zhukov, V.Tyuterev, and E.Chulkov. J.Phys.:Cond.Matt., 24:405802, 2012.
-  V.G.Tyuterev, V.P.Zhukov, P.M.Echenique, and E.V.Chulkov. J.Phys.: Cond.Matt., 27:025801, 2015.
-  G.M.Eliashberg. Interactions between electrons and lattice vibrations in a superconductor. Sov. Phys. JETP, 11:696–702, 1960.
-  R.Loudon. The quantum theory of light. Oxford University Press Inc., New York, 1983.
-  S.Baroni, S.de Gironcoli, and A.Dal Corso. Rev.Mod.Phys., 73:515, 2001.
-  G.B.Bachelet, D.R.Hamann, and M.Schlüter. Phys.Rev.B, 26:4199, 1982.
-  J.P.Perdew and A.Zunger. Phys.Rev.B, 23:5048, 1981.
-  A.DalCorso, S.Baroni, R.Resta, and S.deGironcoli. Phys.Rev.B, 47:3588, 1993.
-  A.Benassi. PW SCF s epsilon.x users manual. Tech. rep., Physics Department, Universita degli Studi di Modena e Reggio Emilia, and INFM/S3 (Nanostructure and Biosystem at Surfaces), Viale A. Allegri 9, 42121, Reggio Emilia, Italy, 2008.
-  B.G.Streetman and S.Banerjee. Solid State electronic Devices, page 524. Prentice Hall, New Jersey, 5th edition, 2000.
-  C.Bulutay, B.K.Ridley, and N.A.Zakhleniuk. Phys.Rev.B, 62:15754, 2000.
-  S.Syed, J.B.Heroux, Y.J.Wang, M.J.Manfra, R.J.Molnar, and H.L.Stormer. Appl.Phys.Lett., 83:4553, 2003.
-  T.Ruf, J.Serrano, M.Cardona, P.Pavone, M.Pabst, M.Krisch, M.DAstuto, T.Suski, I.Grzegory, and M.Leszczynski. Phys.Rev.Lett., 86:906, 2001.
-  B.C.Lee, K.W.Kim, M.Dutta, and M.A.Stroscio. Phys.Rev.B, 56:997, 1997.
-  N.Stanton, A.Kent, A.Akimov, P.Hawker, T.Cheng, and C.T.Foxon. J.Appl.Phys., 89:973, 2001.
-  S.J.Sheih, K.Tsen, D.Ferry, A.Botchkarev, B.Sverdlov, A.Salvador, and H.Morko. Appl.Phys.Lett., 71:1852, 1995.
-  S.S.Dasgupta, J.L.Nidhi, A.Raman, C.Hurni, G.Gupta, J.Speck, and U.K.Mishra. Appl.Phys.Express, 6:034002, 2013.