- Email: [email protected]
Investigations of the recovery of the phonon density of states from specific heat: numerical study
Physics Letters A (1992) 201—209 North-Holland
PHYSICS LETTERS A
Investigations of the recovery of the phonon density of states from specific heat: numerical study Ben-yuan Gu, Bi-zhen Dong and Guo-zhen Yang China Center ofAdvanced Science and Technology (World Laboratory), P.O. Box 8730, Beijing 100080, China and Inst itute ofPhysics, Academia Sinica, P.O. Box 603, Beijing 100080, China
Received 10 February 1992; revised manuscript received 13 April 1992; accepted for publication 24 August 1992 Communicatedby A.R. Bishop
An iterative algorithm for the determination of the phonon density of states (DOS) from the specific heat is presented based on the general theory ofthe amplitude—phase retrieval problem. From the model calculation, it is found that a convergent solution can always be obtained. However, the profile of the recovered DOS depends on the initial frequency distribution used in the algorithm since the inverse problem is inherently ill-posed. In order to obtain the ideal reconstruction of the DOS, it is necessary to introduce additional information about the character of the DOS to guide the choice ofthe initial distribution. The algorithm can also give a robust reconstruction of the DOS from specific heat data contaminated with random noise. The recovery of the DOS based on specific heat data in a confined temperature range has also been examined in detail.
1.
Introduction
The inverse problem for determining the phonon density of states (DOS) from specific heat has been extensively investigated for a long time. Various approaches have been presented to solve this inverse source problem in the past [1—81,but the inversion for the DOS has not been successfully solved up to date. At a very early time, the Fourier deconvolution method has been proposed for solving this inversion [1—3], but the numerical implementation of this method has not been carried out until 1990. Chen et al. implemented a detailed numerical calculation in terms of a FFT algorithm and showed the existence and uniqueness of this inverse solution [6]. However, they found that the Fourier deconvolution procedure is extremely unstable when dealing with simulation data contaminated with random noise, which is inevitable in measurements. The reason leading to the instability of the solution is that this inverse problem is inherently ill-posed. Any inversion method always suffers from this knotty problem. The Tinldionov regularization approach has been
applied to handle this kind of inverse problems [7— 101. This approach demands the construction of an auxiliary nonnegative function incorporated into the relevant equation, which serves to introduce additional information about the unknown solution. The approach also needs the determination of the value of the adjustable regularization parameter a. There exists no general criterion for the choice of an optimum a. It is determined from experience. The inverse solution obtained by this approach is unstable with the parameter a. Recently, based on the Möbius inverse formula, Chen has derived the exact closed-form solution formulated in terms of the inverse Laplace transform [6,111. The expression for this exact solution suffers from a defect with regard to the implementation with realistic data because it requires a numerical cornputation of the inverse Laplace transform. However, the Laplace inverse is numerically unstable and not every function is Laplace invertible [121. For typical experimental data this leads to intractable problems from both an analytical and a numerical point of view. The purpose of the present paper is to offer a new
0375-9601/92/$ 05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved.
201
Volume 170, number 3
PHYSICS LETTERS A
algorithm for the determination of the DOS from specific heat based on the general theory of the amplitude—phase retrieval problem [13—1 5]. From the model calculation it shows that a convergent solution can always be obtained, but the solution depends on the initial frequency distribution used in the algorithm as this inverse problem is inherently ill-posed. In order to get an ideal reconstruction it is necessary to have additional information about the characteristic of the DOS to guide the choice of the initial frequency distribution. Our algorithm can offer a robust construction of the DOS from specific heat data corrupted with random noise. The recovcry of the DOS based on the specific heat data in a confined temperature range has also been examined. The outline of this paper is as follows. In section 2 we describe the relevant formulas and an iterative algorithm. The numerical results of the simulation calculation for the reconstruction of the DOS are presented in section 3. Finally, we summarize the important findings and give concluding remarks.
2 November 1992
of the first kind (eq. (1)), which is known to be inherently ill-posed due to the instability of its solution [2,31. Therefore, this problem has to be handled with great care. No systematical method can be available. Introducing the new variables x and v as ~ k ~— T e = e — ( eq. (1) can be rewritten [1—3,61: .
I
C(.v) =
J
~i(y—x)G(x)
dx=~(v)*G(v),
(4)
where C( v) =
-~
exp ( — i~y) C
1. ( e ~)
k~
(5a)
~ —ln(llw 1k) —
and G(x) =exp[ (1 —~j)x]g((k/h)ex)
,
(Sb)
and the integral kernel is the known Planck—Boltzmann one:
2. Description of the formulas and iterative algorithm The specific heat of a lattice vibration can be expressed as [16]
7 (hw/kT)2exp(hw/kT)
C~(F)=rkJ
[exp(hw/kT)—l]2
g(w)dw, (1)
where h is the Planck constant, k is the Boltzmann constant, T is the absolute temperature, w is the circular frequency and Wm the maximum phonon circular frequency. In addition, the phonon density of states is normalized to 3N:
J
g(w)dw=3N.
(2)
Equation (1) is valid for a crystalline lattice with r atoms per unit cell within the harmonic approximation even for high temperature. The inverse prob1cm can be addressed as follows: How to determine the DOS g(w) from the measured lattice specific heat data at constant volume C~(T)? This problem corresponds to solving the Fredholm integral equation 202
—
2
~
(6) [exp(e”)—11 The asterisk stands for the Fourier convolution between the functions 1 and G. The adjustable parameter ij is chosen as 0
or finite.
(7b)
For convenience of the description below, we reform eq. (4) to the operator form (8) where w expresses the linear integral operator. In general Z may be a nonunitary operator. We then have (9) where the plus stands for the Hermitian conjugation operator and I represents the identity transform. A is a Hermitian operator. As the operator cui belongs to the so-called “ill-posed” ones, in general its in-
Volume 170, number 3
PHYSICS LETTERS A
verse operator is nonexistent. Consequently, the task to seek the inverse solution becomes difficult, From the viewpoint of the general amplitude— phase retrieval problem, this inverse source problem can be described as: If we know the nonunitary linear transform ui and have the information on C, how can we recover the unknown information on G? Moreover, the recovered DOS, G, must satisfy eq. (8) to high accuracy. In other words, we need to reconstruct a G to make d~Gapproach closer to C to high accuracy. To describe the degree of approach of hG to C we introduce the quantity called “approach distance”, defined by [13—15] 1/2
D(G)=~IC—hGI~= (Jdy C(y)...hG~2~ / (10) If D=0 then C=hG. Consequently the retrieval problem for the DOS is to search the extrema of the functional D with respect to the function argument G, which requires that
Then we can derive the following equation satisfied by G (see refs. [13—15] for details), [p(x)+AD]G(x)=cui~C(y)~ANG(x),
(12)
where A~= (h ~ contains only diagonal elements; this is in general a diagonal, not an identity operator, and AN = (cui ~ represents an operator that contains only all the nondiagonal elements of ~ + h, all the diagonal elements being zero. p(x) is any arbitrary nonnegative real function. When G(x) approaches the actual solution, one has p(x)—*0. We assume that p(x) =0. The final equation can be written \
From model calculations we find that the convergent solution always exists to different degrees of accuracy. Although the recovered DOS deviates from the model DOS and takes different shapes, depending on the initial frequency distribution, all the solutions indeed satisfy eq. (8) to an acceptable accuracy, substantially exceeding the precision of the measurement for the specific heat. Consequently, we believe that the recovered DOS has in part brought reasonable information about the real DOS, even losingthe sharp spike structure arising from the van Hove singularities in the frequency spectra when using an inappropriate initial distribution. As thermodynamical quantities, for example, specific heat, electrical resistivity of metals and carrier concentration in the conduction band or valence band of semiconductors, are often insensitive to detailed structures in g(w) and the DOS always contains van Hove singularities in the high frequency regime [161, it is not expected that the DOS can be exactly reconstructed . from specific heat data without any additional information about the DOS. .
.
(11)
öGDO.
i’i
2 November 1992
2—i ±r’( \ ~ \ 1 D r~?~ ~ ~~~kY)—AN~.J’,X)J -
~,l3,
This equation can be solved by an iterative algonthm, as described by refs. [13—151,which starts with some initial distribution function for G (x). The iteration is terminated when the predetermined accuracy is reached. By choosing an appropriate initial DOS distribution for G and performing the iteration procedures, we can obtain the inverse solution for the DOS from specific heat to a certain accuracy.
3. Numerical results based on the YIG algorithm To demonstrate the application of the Yang—Gu (Y/G) algorithm to the inversion of the phonon density of states, we need to perform a particular calculation. We first reform eq. (8) to an approximate form to easily handle the discrete specific heat data. As described by ref. [6], the specific heat of a real system is finite, h(u) should be an integrable function for the whole range of u, all the 1 (u), C(y) and G (x) are Fourier band-limited functions, so one can carry out a Fourier transform. According to the sampling theorem, the original continuous functions, C(y) and G(x), can be represented by a set of discrete sampling values. If the numbers of sampling points are N1 and N2 for x and y, respectively, then G(x) and C(y) can be expressed as single column matrices, and the linear operator h corresponds to an N~x N2 matrix. Thus we have Ni
C~=
~
~
G, j
1, 2, 3,
..., N2,
(14)
and the normalized condition for the DOS g(w) can be expressed as 203
Volume 170, number 3
PHYSICS LETTERS A
v
~ g, = 3N.
(15)
i=1
In the finite-dimensional algebra, the relevant iterative equation can be reformed to
=( ~
m+
(5J
ti5~
C,
—
/
~ A,
1G~m)) A,,
I +G
(16)
~
and the error is estimated by the sum squared error (SSE) defined by
r
SSE=[
.v2 / ~
~C,
~ ~v1
—
]/
\21 cI,1Gimo))
.~, ~
C~-
(17)
where G,~’”°~ is the final convergent solution. To guarantee obtaining the convergent solution when using the Yang—Gu algorithm, the property of the transform kernel h plays a critical role. Although the transform can be nonunitary, one demands the mean value of the diagonal elements of A = h should be larger than that of the nondiagonal elements. If the kernel h(x, ~j’) is a space-confined function in the x-space for any p the above condition
values of the parameter ij. All the shown curves are space-confined functions. The larger the value of ij, the narrower the width is. The application of the Y/G algorithm to the recovery of the DOS is demonstrated with materials with a face-centered-cubic structure. A very similar DOS distribution g(w) to the experimental ones for Cu, Al and Pb can be seen [6]. In this calculation. we produce an artificial simulation DOS as a representative for these materials, as shown in fig. 2. Curve (a) represents a completed DOS coveringthe finite frequency regime [0.0, 4.53] in units of 10” rad/s. Hereafter we always measure the circular frequency in this unit. Curve (b) corresponds to a lowfrequency-cut DOS case with the spectral regime [1.0. 4.53 1. The two curves are vertically offset by an amount of 1.0 for clarity. The modified specific heat C(v) as a function of ~vis depicted in fig. 3, taking i~= 1.2. The temperature is scanned from 10 to 1000 K. Curve (a) is related to the completed DOS case but curve (b) corresponds to the low-frequency-cut DOS case. It is evident that curve (a) has a longer
is satisfied.
I
The behavior of the transform kernel cI(u) as a function of u=y—x is shown in fig. I for different
3
2 November 1992
I
I
6
I ~
3riil/.1) ..,~..- -~-.--
0 -4
0
4
-
.
S
Fig. 1. Behavior of the transform kernel ti(u) for several values of the parameter ~j. 204
(x t()’ Fig. 2. The model frequency distribution g(w) for the phonon density of states (DOS) of materials with a face-centeredcubic structure. Curve (a) corresponds to a completed DOS distribution, curve (b) to a low-frequency-cut DOS with spectral regime [1.0, 4.53] in units of l0’~rad/s. The curves are veriically offset by an amount of 1 for clarity. .
.
Volume 170, number 3
PHYSICS LETTERS A
,.“
temperature from 10 to 1000 K. The solid line corresponds to the completed DOS case and the dots to the low-frequency-cut DOS case. We find that they are nearly coincident. Ofcourse, the behaviorof C(y) depends on ~, but C~(T) is independent of it. Increasing i~the tail in C(y) becomes longer and its magnitude increases in the small y region. To reduce the number of sampling points, one demands that all functions G(x), C(y) and 1~(u)are confined in finite regions as narrow as possible. Generally speaking, the parameter ‘i can take any value within the range [0,3]. As shown in fig. 1, the larger the value
I
I
I
I
of ,j, the narrower the spread width of cb. From this consideration, a larger ~jseems to be favorable. However, the extended region of C(y) becomes broader
0.0
2.5
5.0
7.5
y
(In K
I
I
I
I
I
~
2 November 1992
0.1
a 0.0
A )
increasing the value of ‘i, Although the calculation results are not affected by the choice of the value of the number of sampling points needed in the numerical computation should be changed and can be determined by the sampling theorem to guarantee the accuracy of the solution. We prefer to choose a compromise value for i~, i.e. ,~= 1.2. ~
Fig. 3. Reformed specific heat C(y) against y (=ln T), calculatedbased on the model g(w). ~= 1.2 and the temperature range is [10, 1000] K. Curve (a) corresponds to the completed DOS and curve (b) to the low-frequency-cut DOS. 30
E
I
I
20
~
0 0
I
I
250
500
750
1000
T(K) Fig. 4. Calculated specific heat C,~(T) corresponding to the model DOS. The solid line is for the complete DOS case and the dots are related to the low-frequency-cut DOS.
tail in the small value region of y compared with curve (b). We display the calculated specific heat C~(T) as a function of the temperature Tin fig. 4, scanning the
3.1. Dependence of the convergent solution on the initial frequency distribution ofthe DOS low-frequency-cut DOS in the spectral regime [1.0, COrn], with Wrn=4.53. The temperature is scanned over the range [10, 1000] K. The numbers of sampling points are N~= 128 and N2 = 64, according to In the calculation we set j= 1.2, r= 1, and use the the sampling theorem. The number of iterations is usually 300. The CPU time is only 1—2 mm for 300 iterations on a typical personal computer, for example, an IBM PC 386/33 MHz. We display the convergent solution obtained in figs. 5—7 for different initial distributions. In all these figures curves (a) correspond to the initial distribution, curves (b) to the recovered function of G(x). For comparison, the model function of G(x) is also shown, by curves (c). The initial frequency distribution is chosen as follows: A tilted straight line (fig. 5), a single Gaussian function (fig. 6) and a dual Gaussian function (fig. 7). It is clearly seen that the shape of the recovered G(x) depends on the initial frequency distribution. However, all the recovered G(x) show a space confinement function which reveals the substantial characteristic of the DOS, and 205
PHYSICS LETTERS A
Volume 170, number 3
I
2 November 1992
0.3
I
I
I
0.2
1
~0.2
Fig. 5. Profile of the recovered G(x) from the known C(y) starting with a tilted line for the initial distribution. Curve (a) is the initial distribution and curve (b) the recovered G(x). Curve (c) is the model function of G(x) for comparison. ~= 1.2 and SSE=0.234x l0~. _____________________________________
0.1
.,~‘
,/
I
a
I’
~‘,
:
I~
•‘“. 0.0
b
4
5
______
~-.---.—
6
(ln K)
Fig. 6. The recovery ofG(x) starting with a single Gaussian func. Curve (a) is the initial distribution curve (b) is the reconstructed function of G(x) and curve (c) is the model function of 0(x). ~= 1.2 and SSE=0.903x ~
tion.
keep the total area below the curve constant. This means that the normalization condition forg(w) remains unchanged. All the sum squared errors are less than l0~.On the other hand, the sharp peak structure in the DOS which corresponds to the van Hove singularities is always rounded. The dependence of the shape ofthe recovered G(x) on the choice of the 206
(In K
Fig. 7. Behavior of the recovered function forG(x) starting with a dual Gaussian function. Curve (a) is the initial distribution, curve (b) the recovered G(x) and curve (c) the model function ofG(x). ,~= 1.2 and SSE=0.588x i0~.
To seek the appropriate initial distribution one can take some measures. For instance, one first carries out the recovering procedure for the DOS in the finite low-frequency regime without the van Hove sin7.
‘,.
I
x
exist two peaks in the DOS, one can then use a dual Gaussian function as amay favorable initial distribution, a satisfactory solution be found as shown in fig.
“ ,
6
initial frequency distribution reflects the fact that this inverse problem is inherently ill-posed. In order to obtain the ideal reconstruction of the DOS, one has to introduce additional information about the DOS structure to serve as a guide on the choice of the mitial distribution. For instance, if one knows that there
0.2
I
4
gularities, then uses this result to serve as the initial . . . distribution in the low-frequency regime complementing the other part of the function in the highfrequency regime. Starting with this initial distribution, it is expected that a nearly ideal solution can be obtained. In order to visualize the difference between the model function of C(y) and the calculated C(y) by substituting the recovered G(x) into eq. (14), we present a comparison in fig. 8. The recovered G(x) used in the calculation is taken from curve (b) in fig. 5. The SSE is 0.234x l0~. In fig. 8 the solid line
Volume 170, number 3
PHYSICS LETTERSA
0.2
2 November 1992
0.2
I
I
0.1 I
00.1 0
0.0 2
4
6 y(InK)
I I
Fig. 8. Comparison between the model C(y) (solid line) and the calculated C(y) (circles) by substituting the recovered 0(x) into eq. (14). The corresponding recovered G(x) is shown by curve (b) in fig. 5. The SSEO.234x 10g.
corresponds to the model C(y) and the circles represent the calculated one. It is evident that both data sets are in good agreement to high accuracy, with errors much smaller than the practical measurement errors of specific heat, about 1~Solo. 3.2. Recovery of the DOS from specific heat data contaminated with random noise The realistic specific heat data inevitably suffer from the contamination ofrandom noise arising from the measured error in experiment or the difference of the quantity of the samples. To simulate it, we introduce noise into C(y) at random. The shape of the contaminated C(y) with a noise/signal ratio of 5% is depicted in fig. 9, shown by circles. The original C(y) without noise is shown by the solid line. We chose the dual Gaussian function as the initial distribution, as indicated by curve (a) in fig. 10. The convergent solution is shown by curve (b). For comparison, the model G (x) is also shown by curve (c). The sum squared error reaches 0.875 X 10~.It is clearly seen that the recovered G(x) exhibits a smooth behavior without any observable oscillations and its profile is very similar to that of curve (b) in fig. 7. This example sufficiently reveals the fact that the Y/G algorithm offers a relatively robust recon-
2
I
4
6
y (In K) Fig. 9. Modified function C(y) contaminated with random noise. The solid curve is the ideal model function of C(y) without any noise and the circles represent the noise corrupted data. 0.3
I
I
~I I
0.2
a
~ o.i
0.0
6 x (In K
)
Fig. 10. The recovery of the function G(x) from the noise corrupted data C’(y). The noise/signal ratio is 5%. The initial distribution is taken as a dual Gaussian function as shown in curve (a). The recovered G(x) is indicated by the solidline (b ).Curve 4. (c)isthemodelfunctionofG(x)., 7=1.2andSSEO.875xl0-
struction for the DOS from realistic data contaminated with random noise. Consequently, it can be concluded that the Y/G algorithm is very effective 207
Volume 170, number 3
PHYSICS LETTERS A
and vital in solving the inverse problem in practice.
~
3.3. Reconstruction of the DOS from specific heat data in a finite temperature range An interesting problem is: Can we recover a part of the phonon density of states from specific heat data in a finite temperature range? Due to the transform function cb involved, it exhibits the character of the space confinement with a narrow width. This means that, roughly speaking, the value of specific heat at low temperatures basically depends on the profile of g(w) in the low-frequency regime, and specific heat at high temperatures essentially depends on the shape of g( w) in the high-frequency regime [16]. In other words, the specific heat at low temperatures is dommated by the behavior of g(w) for small w. Troublesome singularities in g(w) often appear at high w. Consequently, it is expected that the reconstruction of the DOS at small w may be easy to perform. We now investigate this problem. We assume that the involving temperature range is [1, 300] K. We use the relation KTm = hWrn to determine the upper cut-off frequency COrn if the upper limit of the temperature ~5 Tm because the vibration modes with frequency higher than COrn do not seem to be excited at this point. For instance, when the temperature range covers [1, 300] K the active spectral regime is approximately [0.1, 4] in units of 1013 rad/s. In this case, where both the temperature and the DOS are confined in some regimes, the recovered functions of G(x) are shown in figs. 11—13 for several initial distributions. In all these figures curves (a) correspond to the initial distribution, curves (b) to the recovered function, curves (c) to the original model of G(x). Both the recovered functions for G(x) in figs. II and 12 show a similar shape even starting with quite different initial distributions. Although their behaviors deviate from the ideal one, the important peak structure in the DOS has been basically recovered. The spike-like structure in the DOS is rounded. The area under the recovered function is kept constant, which means that the normalization condition for the DOS holds well. When choosing a single Gaussian function as the initial distribution the recovered G(x) exhibits a single high peak with a small bump in the low-frequency regime, as shown in fig. 13. 208
2 November 1992
~
0.2
~
~
0.1
0.0 3.0
4.5
x
6.0
(In K)
Fig. 11. Reconstruction of the DOS from a finite segment of the specific heat curve. The confined temperature range is [1. 3001 K. The corresponding frequency is confined in the regime [0.1, 4.01 in units of 103 rad/s. Curve (a) is the initial distribution, curve (b) the recovered function and curve (c) the original model function of G(x). ~= 1.2 and SSE = 0.445 x I 0~.
I
- I
0.2
0.1
0.0
1
I
I
I
I
2
3
4
5
6
x (In K
Fig. 12. Same as fig. 11, except that the initial distribution function is a parabolic function. SSE=0.457x l0~.
Volume 170, number 3
I
PHYSICS LETTERS A
I
inverse problem difficult. In this sense, our algorithm has offered much useful information about the DOS from specific heat data. To gain the ideal reconstruction of the DOS it is necessary to have additional information about the DOS to serve as a guide to choose a favorable initial frequency distribution. Our algorithmcan also deal with realistic data
.....
able reconstruction can This be obtained in the of noise corrupted data. reflects the factpresence that our algorithm has a large capability to anti-interference contaminated with random noise. A relatively relinoises. Consequently, we deem that our algorithm is valuable and vital to carry out the recovery of the DOS from the really measured specific heat data in practice.
I
I
~
0.2
.‘.
,“~
‘I I
0.1
~I
.‘
‘I
/
/
I
I~
~
0.0
1.5
‘. I
3.0
2 November 1992
I
4.5
I
6.0
(In K)
Fig. 13. Same as fig. 11, except that the initial distribution is a single Gaussian function. SSE=0.234x i0~.
Acknowledgement This work was supported by the Chinese National Natural Science Foundation.
4. Discussions and concluding remarks In this work we have presented a numerical investigation of the reconstruction of the DOS from specific heat. This iterative algorithm is based on the general theory of amplitude—phase retrieval and an iterative procedure. We used a model phonon density of states to demonstrate the new algorithm. The recovered function of G(x) depends on the initial distribution since this inverse problem is inherently ill-posed. All the recovered DOS show a bandlimited characteristic and the normalization condition holds. The inverse problem for the DOS is known to be ill-posed with inherent instabilities of the solution. The general conclusion is that only a limited amount of useful information about the phonon density of states can be extracted from the lattice specific heat. On the other hand, thermodynamical quantities, such as specific heat, electrical resistivity of metals and carrier concentrations in the conduction band or valence band of semiconductors, are usually insensitive to the detailed structures in the DOS and the van Hove singularities always appear in the DOS which leads to the appearance of a spike-like structure in the DOS. Al] these reasons make the solution of the
References [1] E.W. Montroll, J. Chem. Phys. 10 (1942) 218. [2]I.M. Lifshitz, Zh. Eksp. Teor. Fiz. 26 (1954) 551. [3] R.G. Chambers, Proc. Phys. Soc. 78 (1961) 941. [4] 0. Weiss, Prog. Theor. Phys. Japan 22 (1959) 526. [5] J.W. Loram, J. Phys. C 19 (1986) 6113. [6] N.X. Chen, Y. Chen and G.Y. Lin, Phys. Lett. A 149 (1990) [7] C.W. Groetsch, The theory ofTinkhonov regularization for Fredhoim equation of the first kind (Pitman, London, 1984).
[81A.J.
Pindor, in: Lecture notes in physics, Vol 115. Modern trends in the theory of condensed matter, eds. A. Pekalski andJ. Przystawa (Springer, Berlin, 1980) pp. 563—570. [9] X.G. Sun and DL. Jaggard, J. AppI. Phys. 62 (1987) 4382. [10] J. Igalson, A.J. Pindor and L. Sinadower, J. Phys. F 11 (1981)995. [11] N.X. Chen, Phys. Rev. Lett. 64 (1990) 1193. [12] R.V. Churchill, Operational mathematics (McGraw-Hill, New York, 1958). [13] BY. Gu and G.Z. Yang, Acta Opt. Sin. 1 (1981) 517; G.Z. Yang and BY. Gu, Acta Phys. Sin. 30 (1981) 410. [14]G.Z. Yang, L. Wang, B.Z. Dong and BY. Gu, Optik 75 (1987) 68. [15] B.Z. Dong, (1991) 1048.BY. Gu and G.Z. Yang, J. Opt. Soc. Am. A 8 [16] W. Jones and N.H. March, Theoretical solid state physics, Vol. 1 (Wiley, New York, 1973) pp. 238—247.
209
PHYSICS LETTERS A
Investigations of the recovery of the phonon density of states from specific heat: numerical study Ben-yuan Gu, Bi-zhen Dong and Guo-zhen Yang China Center ofAdvanced Science and Technology (World Laboratory), P.O. Box 8730, Beijing 100080, China and Inst itute ofPhysics, Academia Sinica, P.O. Box 603, Beijing 100080, China
Received 10 February 1992; revised manuscript received 13 April 1992; accepted for publication 24 August 1992 Communicatedby A.R. Bishop
An iterative algorithm for the determination of the phonon density of states (DOS) from the specific heat is presented based on the general theory ofthe amplitude—phase retrieval problem. From the model calculation, it is found that a convergent solution can always be obtained. However, the profile of the recovered DOS depends on the initial frequency distribution used in the algorithm since the inverse problem is inherently ill-posed. In order to obtain the ideal reconstruction of the DOS, it is necessary to introduce additional information about the character of the DOS to guide the choice ofthe initial distribution. The algorithm can also give a robust reconstruction of the DOS from specific heat data contaminated with random noise. The recovery of the DOS based on specific heat data in a confined temperature range has also been examined in detail.
1.
Introduction
The inverse problem for determining the phonon density of states (DOS) from specific heat has been extensively investigated for a long time. Various approaches have been presented to solve this inverse source problem in the past [1—81,but the inversion for the DOS has not been successfully solved up to date. At a very early time, the Fourier deconvolution method has been proposed for solving this inversion [1—3], but the numerical implementation of this method has not been carried out until 1990. Chen et al. implemented a detailed numerical calculation in terms of a FFT algorithm and showed the existence and uniqueness of this inverse solution [6]. However, they found that the Fourier deconvolution procedure is extremely unstable when dealing with simulation data contaminated with random noise, which is inevitable in measurements. The reason leading to the instability of the solution is that this inverse problem is inherently ill-posed. Any inversion method always suffers from this knotty problem. The Tinldionov regularization approach has been
applied to handle this kind of inverse problems [7— 101. This approach demands the construction of an auxiliary nonnegative function incorporated into the relevant equation, which serves to introduce additional information about the unknown solution. The approach also needs the determination of the value of the adjustable regularization parameter a. There exists no general criterion for the choice of an optimum a. It is determined from experience. The inverse solution obtained by this approach is unstable with the parameter a. Recently, based on the Möbius inverse formula, Chen has derived the exact closed-form solution formulated in terms of the inverse Laplace transform [6,111. The expression for this exact solution suffers from a defect with regard to the implementation with realistic data because it requires a numerical cornputation of the inverse Laplace transform. However, the Laplace inverse is numerically unstable and not every function is Laplace invertible [121. For typical experimental data this leads to intractable problems from both an analytical and a numerical point of view. The purpose of the present paper is to offer a new
0375-9601/92/$ 05.00 © 1992 Elsevier Science Publishers B.V. All rights reserved.
201
Volume 170, number 3
PHYSICS LETTERS A
algorithm for the determination of the DOS from specific heat based on the general theory of the amplitude—phase retrieval problem [13—1 5]. From the model calculation it shows that a convergent solution can always be obtained, but the solution depends on the initial frequency distribution used in the algorithm as this inverse problem is inherently ill-posed. In order to get an ideal reconstruction it is necessary to have additional information about the characteristic of the DOS to guide the choice of the initial frequency distribution. Our algorithm can offer a robust construction of the DOS from specific heat data corrupted with random noise. The recovcry of the DOS based on the specific heat data in a confined temperature range has also been examined. The outline of this paper is as follows. In section 2 we describe the relevant formulas and an iterative algorithm. The numerical results of the simulation calculation for the reconstruction of the DOS are presented in section 3. Finally, we summarize the important findings and give concluding remarks.
2 November 1992
of the first kind (eq. (1)), which is known to be inherently ill-posed due to the instability of its solution [2,31. Therefore, this problem has to be handled with great care. No systematical method can be available. Introducing the new variables x and v as ~ k ~— T e = e — ( eq. (1) can be rewritten [1—3,61: .
I
C(.v) =
J
~i(y—x)G(x)
dx=~(v)*G(v),
(4)
where C( v) =
-~
exp ( — i~y) C
1. ( e ~)
k~
(5a)
~ —ln(llw 1k) —
and G(x) =exp[ (1 —~j)x]g((k/h)ex)
,
(Sb)
and the integral kernel is the known Planck—Boltzmann one:
2. Description of the formulas and iterative algorithm The specific heat of a lattice vibration can be expressed as [16]
7 (hw/kT)2exp(hw/kT)
C~(F)=rkJ
[exp(hw/kT)—l]2
g(w)dw, (1)
where h is the Planck constant, k is the Boltzmann constant, T is the absolute temperature, w is the circular frequency and Wm the maximum phonon circular frequency. In addition, the phonon density of states is normalized to 3N:
J
g(w)dw=3N.
(2)
Equation (1) is valid for a crystalline lattice with r atoms per unit cell within the harmonic approximation even for high temperature. The inverse prob1cm can be addressed as follows: How to determine the DOS g(w) from the measured lattice specific heat data at constant volume C~(T)? This problem corresponds to solving the Fredholm integral equation 202
—
2
~
(6) [exp(e”)—11 The asterisk stands for the Fourier convolution between the functions 1 and G. The adjustable parameter ij is chosen as 0
or finite.
(7b)
For convenience of the description below, we reform eq. (4) to the operator form (8) where w expresses the linear integral operator. In general Z may be a nonunitary operator. We then have (9) where the plus stands for the Hermitian conjugation operator and I represents the identity transform. A is a Hermitian operator. As the operator cui belongs to the so-called “ill-posed” ones, in general its in-
Volume 170, number 3
PHYSICS LETTERS A
verse operator is nonexistent. Consequently, the task to seek the inverse solution becomes difficult, From the viewpoint of the general amplitude— phase retrieval problem, this inverse source problem can be described as: If we know the nonunitary linear transform ui and have the information on C, how can we recover the unknown information on G? Moreover, the recovered DOS, G, must satisfy eq. (8) to high accuracy. In other words, we need to reconstruct a G to make d~Gapproach closer to C to high accuracy. To describe the degree of approach of hG to C we introduce the quantity called “approach distance”, defined by [13—15] 1/2
D(G)=~IC—hGI~= (Jdy C(y)...hG~2~ / (10) If D=0 then C=hG. Consequently the retrieval problem for the DOS is to search the extrema of the functional D with respect to the function argument G, which requires that
Then we can derive the following equation satisfied by G (see refs. [13—15] for details), [p(x)+AD]G(x)=cui~C(y)~ANG(x),
(12)
where A~= (h ~ contains only diagonal elements; this is in general a diagonal, not an identity operator, and AN = (cui ~ represents an operator that contains only all the nondiagonal elements of ~ + h, all the diagonal elements being zero. p(x) is any arbitrary nonnegative real function. When G(x) approaches the actual solution, one has p(x)—*0. We assume that p(x) =0. The final equation can be written \
From model calculations we find that the convergent solution always exists to different degrees of accuracy. Although the recovered DOS deviates from the model DOS and takes different shapes, depending on the initial frequency distribution, all the solutions indeed satisfy eq. (8) to an acceptable accuracy, substantially exceeding the precision of the measurement for the specific heat. Consequently, we believe that the recovered DOS has in part brought reasonable information about the real DOS, even losingthe sharp spike structure arising from the van Hove singularities in the frequency spectra when using an inappropriate initial distribution. As thermodynamical quantities, for example, specific heat, electrical resistivity of metals and carrier concentration in the conduction band or valence band of semiconductors, are often insensitive to detailed structures in g(w) and the DOS always contains van Hove singularities in the high frequency regime [161, it is not expected that the DOS can be exactly reconstructed . from specific heat data without any additional information about the DOS. .
.
(11)
öGDO.
i’i
2 November 1992
2—i ±r’( \ ~ \ 1 D r~?~ ~ ~~~kY)—AN~.J’,X)J -
~,l3,
This equation can be solved by an iterative algonthm, as described by refs. [13—151,which starts with some initial distribution function for G (x). The iteration is terminated when the predetermined accuracy is reached. By choosing an appropriate initial DOS distribution for G and performing the iteration procedures, we can obtain the inverse solution for the DOS from specific heat to a certain accuracy.
3. Numerical results based on the YIG algorithm To demonstrate the application of the Yang—Gu (Y/G) algorithm to the inversion of the phonon density of states, we need to perform a particular calculation. We first reform eq. (8) to an approximate form to easily handle the discrete specific heat data. As described by ref. [6], the specific heat of a real system is finite, h(u) should be an integrable function for the whole range of u, all the 1 (u), C(y) and G (x) are Fourier band-limited functions, so one can carry out a Fourier transform. According to the sampling theorem, the original continuous functions, C(y) and G(x), can be represented by a set of discrete sampling values. If the numbers of sampling points are N1 and N2 for x and y, respectively, then G(x) and C(y) can be expressed as single column matrices, and the linear operator h corresponds to an N~x N2 matrix. Thus we have Ni
C~=
~
~
G, j
1, 2, 3,
..., N2,
(14)
and the normalized condition for the DOS g(w) can be expressed as 203
Volume 170, number 3
PHYSICS LETTERS A
v
~ g, = 3N.
(15)
i=1
In the finite-dimensional algebra, the relevant iterative equation can be reformed to
=( ~
m+
(5J
ti5~
C,
—
/
~ A,
1G~m)) A,,
I +G
(16)
~
and the error is estimated by the sum squared error (SSE) defined by
r
SSE=[
.v2 / ~
~C,
~ ~v1
—
]/
\21 cI,1Gimo))
.~, ~
C~-
(17)
where G,~’”°~ is the final convergent solution. To guarantee obtaining the convergent solution when using the Yang—Gu algorithm, the property of the transform kernel h plays a critical role. Although the transform can be nonunitary, one demands the mean value of the diagonal elements of A = h should be larger than that of the nondiagonal elements. If the kernel h(x, ~j’) is a space-confined function in the x-space for any p the above condition
values of the parameter ij. All the shown curves are space-confined functions. The larger the value of ij, the narrower the width is. The application of the Y/G algorithm to the recovery of the DOS is demonstrated with materials with a face-centered-cubic structure. A very similar DOS distribution g(w) to the experimental ones for Cu, Al and Pb can be seen [6]. In this calculation. we produce an artificial simulation DOS as a representative for these materials, as shown in fig. 2. Curve (a) represents a completed DOS coveringthe finite frequency regime [0.0, 4.53] in units of 10” rad/s. Hereafter we always measure the circular frequency in this unit. Curve (b) corresponds to a lowfrequency-cut DOS case with the spectral regime [1.0. 4.53 1. The two curves are vertically offset by an amount of 1.0 for clarity. The modified specific heat C(v) as a function of ~vis depicted in fig. 3, taking i~= 1.2. The temperature is scanned from 10 to 1000 K. Curve (a) is related to the completed DOS case but curve (b) corresponds to the low-frequency-cut DOS case. It is evident that curve (a) has a longer
is satisfied.
I
The behavior of the transform kernel cI(u) as a function of u=y—x is shown in fig. I for different
3
2 November 1992
I
I
6
I ~
3riil/.1) ..,~..- -~-.--
0 -4
0
4
-
.
S
Fig. 1. Behavior of the transform kernel ti(u) for several values of the parameter ~j. 204
(x t()’ Fig. 2. The model frequency distribution g(w) for the phonon density of states (DOS) of materials with a face-centeredcubic structure. Curve (a) corresponds to a completed DOS distribution, curve (b) to a low-frequency-cut DOS with spectral regime [1.0, 4.53] in units of l0’~rad/s. The curves are veriically offset by an amount of 1 for clarity. .
.
Volume 170, number 3
PHYSICS LETTERS A
,.“
temperature from 10 to 1000 K. The solid line corresponds to the completed DOS case and the dots to the low-frequency-cut DOS case. We find that they are nearly coincident. Ofcourse, the behaviorof C(y) depends on ~, but C~(T) is independent of it. Increasing i~the tail in C(y) becomes longer and its magnitude increases in the small y region. To reduce the number of sampling points, one demands that all functions G(x), C(y) and 1~(u)are confined in finite regions as narrow as possible. Generally speaking, the parameter ‘i can take any value within the range [0,3]. As shown in fig. 1, the larger the value
I
I
I
I
of ,j, the narrower the spread width of cb. From this consideration, a larger ~jseems to be favorable. However, the extended region of C(y) becomes broader
0.0
2.5
5.0
7.5
y
(In K
I
I
I
I
I
~
2 November 1992
0.1
a 0.0
A )
increasing the value of ‘i, Although the calculation results are not affected by the choice of the value of the number of sampling points needed in the numerical computation should be changed and can be determined by the sampling theorem to guarantee the accuracy of the solution. We prefer to choose a compromise value for i~, i.e. ,~= 1.2. ~
Fig. 3. Reformed specific heat C(y) against y (=ln T), calculatedbased on the model g(w). ~= 1.2 and the temperature range is [10, 1000] K. Curve (a) corresponds to the completed DOS and curve (b) to the low-frequency-cut DOS. 30
E
I
I
20
~
0 0
I
I
250
500
750
1000
T(K) Fig. 4. Calculated specific heat C,~(T) corresponding to the model DOS. The solid line is for the complete DOS case and the dots are related to the low-frequency-cut DOS.
tail in the small value region of y compared with curve (b). We display the calculated specific heat C~(T) as a function of the temperature Tin fig. 4, scanning the
3.1. Dependence of the convergent solution on the initial frequency distribution ofthe DOS low-frequency-cut DOS in the spectral regime [1.0, COrn], with Wrn=4.53. The temperature is scanned over the range [10, 1000] K. The numbers of sampling points are N~= 128 and N2 = 64, according to In the calculation we set j= 1.2, r= 1, and use the the sampling theorem. The number of iterations is usually 300. The CPU time is only 1—2 mm for 300 iterations on a typical personal computer, for example, an IBM PC 386/33 MHz. We display the convergent solution obtained in figs. 5—7 for different initial distributions. In all these figures curves (a) correspond to the initial distribution, curves (b) to the recovered function of G(x). For comparison, the model function of G(x) is also shown, by curves (c). The initial frequency distribution is chosen as follows: A tilted straight line (fig. 5), a single Gaussian function (fig. 6) and a dual Gaussian function (fig. 7). It is clearly seen that the shape of the recovered G(x) depends on the initial frequency distribution. However, all the recovered G(x) show a space confinement function which reveals the substantial characteristic of the DOS, and 205
PHYSICS LETTERS A
Volume 170, number 3
I
2 November 1992
0.3
I
I
I
0.2
1
~0.2
Fig. 5. Profile of the recovered G(x) from the known C(y) starting with a tilted line for the initial distribution. Curve (a) is the initial distribution and curve (b) the recovered G(x). Curve (c) is the model function of G(x) for comparison. ~= 1.2 and SSE=0.234x l0~. _____________________________________
0.1
.,~‘
,/
I
a
I’
~‘,
:
I~
•‘“. 0.0
b
4
5
______
~-.---.—
6
(ln K)
Fig. 6. The recovery ofG(x) starting with a single Gaussian func. Curve (a) is the initial distribution curve (b) is the reconstructed function of G(x) and curve (c) is the model function of 0(x). ~= 1.2 and SSE=0.903x ~
tion.
keep the total area below the curve constant. This means that the normalization condition forg(w) remains unchanged. All the sum squared errors are less than l0~.On the other hand, the sharp peak structure in the DOS which corresponds to the van Hove singularities is always rounded. The dependence of the shape ofthe recovered G(x) on the choice of the 206
(In K
Fig. 7. Behavior of the recovered function forG(x) starting with a dual Gaussian function. Curve (a) is the initial distribution, curve (b) the recovered G(x) and curve (c) the model function ofG(x). ,~= 1.2 and SSE=0.588x i0~.
To seek the appropriate initial distribution one can take some measures. For instance, one first carries out the recovering procedure for the DOS in the finite low-frequency regime without the van Hove sin7.
‘,.
I
x
exist two peaks in the DOS, one can then use a dual Gaussian function as amay favorable initial distribution, a satisfactory solution be found as shown in fig.
“ ,
6
initial frequency distribution reflects the fact that this inverse problem is inherently ill-posed. In order to obtain the ideal reconstruction of the DOS, one has to introduce additional information about the DOS structure to serve as a guide on the choice of the mitial distribution. For instance, if one knows that there
0.2
I
4
gularities, then uses this result to serve as the initial . . . distribution in the low-frequency regime complementing the other part of the function in the highfrequency regime. Starting with this initial distribution, it is expected that a nearly ideal solution can be obtained. In order to visualize the difference between the model function of C(y) and the calculated C(y) by substituting the recovered G(x) into eq. (14), we present a comparison in fig. 8. The recovered G(x) used in the calculation is taken from curve (b) in fig. 5. The SSE is 0.234x l0~. In fig. 8 the solid line
Volume 170, number 3
PHYSICS LETTERSA
0.2
2 November 1992
0.2
I
I
0.1 I
00.1 0
0.0 2
4
6 y(InK)
I I
Fig. 8. Comparison between the model C(y) (solid line) and the calculated C(y) (circles) by substituting the recovered 0(x) into eq. (14). The corresponding recovered G(x) is shown by curve (b) in fig. 5. The SSEO.234x 10g.
corresponds to the model C(y) and the circles represent the calculated one. It is evident that both data sets are in good agreement to high accuracy, with errors much smaller than the practical measurement errors of specific heat, about 1~Solo. 3.2. Recovery of the DOS from specific heat data contaminated with random noise The realistic specific heat data inevitably suffer from the contamination ofrandom noise arising from the measured error in experiment or the difference of the quantity of the samples. To simulate it, we introduce noise into C(y) at random. The shape of the contaminated C(y) with a noise/signal ratio of 5% is depicted in fig. 9, shown by circles. The original C(y) without noise is shown by the solid line. We chose the dual Gaussian function as the initial distribution, as indicated by curve (a) in fig. 10. The convergent solution is shown by curve (b). For comparison, the model G (x) is also shown by curve (c). The sum squared error reaches 0.875 X 10~.It is clearly seen that the recovered G(x) exhibits a smooth behavior without any observable oscillations and its profile is very similar to that of curve (b) in fig. 7. This example sufficiently reveals the fact that the Y/G algorithm offers a relatively robust recon-
2
I
4
6
y (In K) Fig. 9. Modified function C(y) contaminated with random noise. The solid curve is the ideal model function of C(y) without any noise and the circles represent the noise corrupted data. 0.3
I
I
~I I
0.2
a
~ o.i
0.0
6 x (In K
)
Fig. 10. The recovery of the function G(x) from the noise corrupted data C’(y). The noise/signal ratio is 5%. The initial distribution is taken as a dual Gaussian function as shown in curve (a). The recovered G(x) is indicated by the solidline (b ).Curve 4. (c)isthemodelfunctionofG(x)., 7=1.2andSSEO.875xl0-
struction for the DOS from realistic data contaminated with random noise. Consequently, it can be concluded that the Y/G algorithm is very effective 207
Volume 170, number 3
PHYSICS LETTERS A
and vital in solving the inverse problem in practice.
~
3.3. Reconstruction of the DOS from specific heat data in a finite temperature range An interesting problem is: Can we recover a part of the phonon density of states from specific heat data in a finite temperature range? Due to the transform function cb involved, it exhibits the character of the space confinement with a narrow width. This means that, roughly speaking, the value of specific heat at low temperatures basically depends on the profile of g(w) in the low-frequency regime, and specific heat at high temperatures essentially depends on the shape of g( w) in the high-frequency regime [16]. In other words, the specific heat at low temperatures is dommated by the behavior of g(w) for small w. Troublesome singularities in g(w) often appear at high w. Consequently, it is expected that the reconstruction of the DOS at small w may be easy to perform. We now investigate this problem. We assume that the involving temperature range is [1, 300] K. We use the relation KTm = hWrn to determine the upper cut-off frequency COrn if the upper limit of the temperature ~5 Tm because the vibration modes with frequency higher than COrn do not seem to be excited at this point. For instance, when the temperature range covers [1, 300] K the active spectral regime is approximately [0.1, 4] in units of 1013 rad/s. In this case, where both the temperature and the DOS are confined in some regimes, the recovered functions of G(x) are shown in figs. 11—13 for several initial distributions. In all these figures curves (a) correspond to the initial distribution, curves (b) to the recovered function, curves (c) to the original model of G(x). Both the recovered functions for G(x) in figs. II and 12 show a similar shape even starting with quite different initial distributions. Although their behaviors deviate from the ideal one, the important peak structure in the DOS has been basically recovered. The spike-like structure in the DOS is rounded. The area under the recovered function is kept constant, which means that the normalization condition for the DOS holds well. When choosing a single Gaussian function as the initial distribution the recovered G(x) exhibits a single high peak with a small bump in the low-frequency regime, as shown in fig. 13. 208
2 November 1992
~
0.2
~
~
0.1
0.0 3.0
4.5
x
6.0
(In K)
Fig. 11. Reconstruction of the DOS from a finite segment of the specific heat curve. The confined temperature range is [1. 3001 K. The corresponding frequency is confined in the regime [0.1, 4.01 in units of 103 rad/s. Curve (a) is the initial distribution, curve (b) the recovered function and curve (c) the original model function of G(x). ~= 1.2 and SSE = 0.445 x I 0~.
I
- I
0.2
0.1
0.0
1
I
I
I
I
2
3
4
5
6
x (In K
Fig. 12. Same as fig. 11, except that the initial distribution function is a parabolic function. SSE=0.457x l0~.
Volume 170, number 3
I
PHYSICS LETTERS A
I
inverse problem difficult. In this sense, our algorithm has offered much useful information about the DOS from specific heat data. To gain the ideal reconstruction of the DOS it is necessary to have additional information about the DOS to serve as a guide to choose a favorable initial frequency distribution. Our algorithmcan also deal with realistic data
.....
able reconstruction can This be obtained in the of noise corrupted data. reflects the factpresence that our algorithm has a large capability to anti-interference contaminated with random noise. A relatively relinoises. Consequently, we deem that our algorithm is valuable and vital to carry out the recovery of the DOS from the really measured specific heat data in practice.
I
I
~
0.2
.‘.
,“~
‘I I
0.1
~I
.‘
‘I
/
/
I
I~
~
0.0
1.5
‘. I
3.0
2 November 1992
I
4.5
I
6.0
(In K)
Fig. 13. Same as fig. 11, except that the initial distribution is a single Gaussian function. SSE=0.234x i0~.
Acknowledgement This work was supported by the Chinese National Natural Science Foundation.
4. Discussions and concluding remarks In this work we have presented a numerical investigation of the reconstruction of the DOS from specific heat. This iterative algorithm is based on the general theory of amplitude—phase retrieval and an iterative procedure. We used a model phonon density of states to demonstrate the new algorithm. The recovered function of G(x) depends on the initial distribution since this inverse problem is inherently ill-posed. All the recovered DOS show a bandlimited characteristic and the normalization condition holds. The inverse problem for the DOS is known to be ill-posed with inherent instabilities of the solution. The general conclusion is that only a limited amount of useful information about the phonon density of states can be extracted from the lattice specific heat. On the other hand, thermodynamical quantities, such as specific heat, electrical resistivity of metals and carrier concentrations in the conduction band or valence band of semiconductors, are usually insensitive to the detailed structures in the DOS and the van Hove singularities always appear in the DOS which leads to the appearance of a spike-like structure in the DOS. Al] these reasons make the solution of the
References [1] E.W. Montroll, J. Chem. Phys. 10 (1942) 218. [2]I.M. Lifshitz, Zh. Eksp. Teor. Fiz. 26 (1954) 551. [3] R.G. Chambers, Proc. Phys. Soc. 78 (1961) 941. [4] 0. Weiss, Prog. Theor. Phys. Japan 22 (1959) 526. [5] J.W. Loram, J. Phys. C 19 (1986) 6113. [6] N.X. Chen, Y. Chen and G.Y. Lin, Phys. Lett. A 149 (1990) [7] C.W. Groetsch, The theory ofTinkhonov regularization for Fredhoim equation of the first kind (Pitman, London, 1984).
[81A.J.
Pindor, in: Lecture notes in physics, Vol 115. Modern trends in the theory of condensed matter, eds. A. Pekalski andJ. Przystawa (Springer, Berlin, 1980) pp. 563—570. [9] X.G. Sun and DL. Jaggard, J. AppI. Phys. 62 (1987) 4382. [10] J. Igalson, A.J. Pindor and L. Sinadower, J. Phys. F 11 (1981)995. [11] N.X. Chen, Phys. Rev. Lett. 64 (1990) 1193. [12] R.V. Churchill, Operational mathematics (McGraw-Hill, New York, 1958). [13] BY. Gu and G.Z. Yang, Acta Opt. Sin. 1 (1981) 517; G.Z. Yang and BY. Gu, Acta Phys. Sin. 30 (1981) 410. [14]G.Z. Yang, L. Wang, B.Z. Dong and BY. Gu, Optik 75 (1987) 68. [15] B.Z. Dong, (1991) 1048.BY. Gu and G.Z. Yang, J. Opt. Soc. Am. A 8 [16] W. Jones and N.H. March, Theoretical solid state physics, Vol. 1 (Wiley, New York, 1973) pp. 238—247.
209