Updated: January 23, 2008   
Elastic Normal Modes in an Anisotropic Elastic
Medium with Cubic Crystal Elastic Coefficients

Calculations of the modes of vibration of a sphere are made taking into account the anisotropicity of the elastic constants. This is applied to silicon. The results are compared with molecular dynamics simulations. The error resulting from approximating the material as isotropic is found to be significant.

   Consider a continuous elastic medium. Assume it is homogeneous (material properties independent of position). But it may be anisotropic (material properties vary depending on direction). The density of the medium is ρ. In continuum elastic mechanics it is necessary to be careful about the distinction between "material coordinates" and "space coordinates." Position in material coordinates is specified by r. If the medium is unstressed and at equilibrium then material coordinates r are the same as space coordinates. Depending on the situation, the Cartesian components of r may be denoted as x, y and z, or else as x1, x2 and x3. In spherical coordinates, the components of r are r, θ and φ, where x = r sin(θ) cos(φ), y = r sin(θ) sin(φ) and z = r cos(θ). The displacement of the point of material whose equilibrium position is r is u(r,t). The Cartesian components of u are u1, u2 and u3. The spherical components of u are ur, uθ and uφ.
   The strain tensor is a measure of material squeezing or twisting. Unfortunately, there are two common definitions of the strain tensor. To avoid confusion here, symbols e and ε will be used, following the approach of Ashcroft and Mermin ("Solid State Physics" page 444, 1976). The first definition (for e) is very commonly used (C. Kittel, etc.)
   ∂ui
-----
∂xj		 ( i = j )
eij =    
  
∂ui
-----
∂xj		 + ∂uj
-----
∂xi
 ( i ≠ j )
   The second definition (for ε) has the advantage that the 3×3 matrix of coefficients transforms under rotation as a second rank tensor.
εij =
1
--
2		( ∂ui
-----
∂xj		 + ∂uj
-----
∂xi		)
   In tensor notation, εij = (1/2)(ui, j + ui, j). The difference between the two definitions is that the off diagonal components of eij are double those of εij.
   The stress tensor σij is related to the strain tensor through
σij = cijkl εkl
where it is implicit that a summation is to be taken over k = 1, 2, 3 and l = 1, 2, 3.
   Since σij, eij and εij are all symmetric 3×3 matrices, it is convenient also to introduce:
σ1 = σ11= σxx   e1 = e11 = exx
σ2 = σ22= σyy   e2 = e22 = eyy
σ3 = σ33= σzz   e3 = e33 = ezz
σ4 = σ23= σyz   e4 = e23 = eyz = 2 ε23
σ5 = σ13= σxz   e5 = e13 = exz = 2 ε13
σ6 = σ12= σxy   e6 = e12 = exy = 2 ε12
   In terms of these,
σα = 6
Σ
β=1		Cαβ eβ
where Cαβ are called the "elastic stiffness constants" or "elastic moduli." If α corresponds to ij or ji and β corresponds to kl or lk then cijkl = Cαβ. [Note that this differs from equation (22.83) on page 446 of Ashcroft and Mermin "Solid State Physics" 1976]
   For clarity, equations that hold only in isotropic elastic media will be coloured in this shade of purple.
   For an isotropic material,
C11 = C22 = C33
C12 = C21 = C23 = C32 = C31 = C13
C44 = C55 = C66 = (1/2)(C11 - C12)
C14 = C41 = C15 = C51 = C16 = C24 = C42 = C25 = C52 = C26 = C62 = ... = 0
Therefore, there are only two independent elastic constants. Other elastic constants for isotropic materials are related as follows:
λ = C12 (first Lamé constant )
μ = C44 ( shear modulus )
λ + 2μ = C11
B = λ + (2/3) μ ( B = dP/(dV/V) = bulk modulus )
μ = 3 B Y / (9B-Y) ( Y = Young's modulus )
ν = (3 B - Y) / (6B) ( ν = Poisson ratio )
Clong = sqrt(C11/ρ) = sqrt((λ+2μ)/ρ) ( longitudinal speed of sound )
Ctrans = sqrt(C44/ρ) = sqrt(μ/ρ) ( transverse speed of sound )
Y = μ(3λ+2μ)/(λ+μ)
ν = λ / (2(λ+μ))
   For a cubic crystalline material,
C11 = C22 = C33
C12 = C21 = C23 = C32 = C31 = C13
C44 = C55 = C66
C14 = C41 = C15 = C51 = C16 = C24 = C42 = C25 = C52 = C26 = C62 = ... = 0
and the Zener anisotropy factor is defined by
A = 2 C44 / (C11 - C12 )
   A cubic crystal becomes isotropic when A equals 1. Otherwise, there is no way to define the shear modulus, Young's modulus or Poisson ratio of a cubic crystal except in an average sense, such as the Voigt-Reuss-Hill average.
   In an isotropic medium, the differential equation for the displacement field [after Ye 2000 equation (3) ] is:

(λ+2μ) grad(div u) - μ curl(curl u) = ρ (∂2 u/∂t2 )
(1.11)

where ρ is density, λ is the first Lamé constant and μ is the shear modulus.
   I now want to generalize this to an anisotropic material. The left side of the equation is the force per unit volume. Consider some region V with a smooth surface, A. Consider a small element of this surface with area dA and outward pointing normal vector dA, where ||dA|| = dA. The force on the enclosed region exerted by the surrounding material is dF = σ·(dA). Specifically, in Cartesian coordinates,
dFi = σij dAj
   Recall that the familiar Green theorem says that for any vector field B, the surface integral of B·dA is equal to the volume integral of ∇·B. An extension of this is the "generalized Stokes theorem" which says, in particular, that for a tensor field σ, the surface integral of σ·dA is equal to the volume integral of ∇·σ. Specifically, in Cartesian coordinates, the surface integral of σijdAj equals the volume integral of (∂/∂xj)σij which can also be written compactly as σij, j.
   Finally, let the volume V be shrunk down to zero size, in which case the force on the enclosed region per unit volume is σij, j. Since we are using material coordinates, we are free to apply Newton's second law F = ma to a small region, resulting in
σij, j = ρ ∂2ui/∂t2
cijkl εkl, j = ρ ∂2ui/∂t2
(1/2) cijkl( uk,lj+ul,kj) = ρ ∂2ui/∂t2
Now, consider a normal mode u(r) with frequency ω (in rad/s).
(1/2) cijkl( uk,lj+ul,kj) = - ω2 ρ ui
(-1/(2ρ)) cijkl( uk,lj+ul,kj) = ω2 ui
Define H to be the linear operator for which
(H u)i = (-1/(2ρ)) cijkl( uk,lj+ul,kj)
then we can write
H u = ω2 u
and the equation for a normal mode is now in the form of an eigenvalue equation. The eigenvalue is ω2.
   What follows is formally analogous to time independent perturbation theory used in quantum mechanics. For that reason a quantum mechanical style notation is used.
Given some Hamiltonian, H', if we want to find the eigenvalues and eigenvectors of the equation
H' |Ψ> = E |Ψ>
it can be convenient to find Ho and λV such that H' = Ho + λV. Choose Ho so that a complete basis of eigenstates |Ψn> and eigenvalues En are known, for which Ho |Ψn> = En |Ψn>. Then by considering an expansion in the parameter λ the eigenstates and eigenvalues of H' may be found. However, in this case the unperturbed eigenstates are all degenerate, and so perturbation theory as an expansion in λ cannot be carried out.
   We will choose Ho to correspond to some set of elastic constants which have isotropic elasticity. H' corresponds to the actual elastic constants. See reference [2] (Trallero-Giner 1998) for a similar approach applied to polar optical modes.
   We will begin by introducing the eigenstates for Ho. Even though not always indicated, these are simply vector-valued functions of position. Conceivably they could be complex valued, even though that is not done until much later in the discussion. The eigenfunctions of Ho are denoted | p k l m >, where
p = 1, 2 or 3 (polarization)
k = radial wavevector (real values from 0 to ∞)
l = angular momentum quantum number (0, 1, 2, 3, ...)
m = azimuthal quantum number (-l, -l+1, ..., 0, l-1, l)
First, define the function fm(φ) by
fm(φ) = cos(m φ)   ( m ≥ 0 )
fm(φ) = sin(m φ)    ( m < 0 )
For p = 1 ("longitudinal spheroidal" polarization):
| 1 k l m >(r) = ∇( jl (k r) Plm(cos(θ)) fm(φ))
with eigenvalues ω2 = (k Clong)2, where Clong is the longitudinal speed of sound.
For p = 2 ("transverse torsional" polarization):
| 2 k l m >(r) = ∇×(r jl (k r) Plm(cos(θ)) fm(φ))
with eigenvalues ω2 = (k Ctrans)2, where Ctrans is the transverse speed of sound.
For p = 3 ("transverse spheroidal" polarization):
| 3 k l m >(r) = ∇×∇×(r jl (k r) Plm(cos(θ)) fm(φ))
with eigenvalues ω2 = (k Ctrans)2, where Ctrans is the transverse speed of sound.
   jl (·) are spherical Bessel functions of the first kind of order l.
Plm(·) are associated Legendre functions.
   These eigenfunctions are not normalized. It is convenient to leave them this way.
   The next important step is to define an inner product on the space of functions. Consider any two states (not necessarily eigenstates) | α > and | β >. This corresponds to two possibly complex, vector valued functions of position uα(r) and uβ(r). The "truncated inner product", denoted < β | α >Rm is the integral over θ and φ and r from 0 to Rm of uα(r)·(uβ(r))* = uαi (r)(uβi (r))*, where the * denotes taking the complex conjugate. This makes it possible to define the normalized matrix element of H',

< p' k' l' m' | H' | p k l m > =
lim
Rm → ∞		  
< p' k' l' m' | H' | p k l m >Rm
---------------------------------------------------------------
sqrt(< p' k' l' m' | p' k' l' m' >Rm < p k l m | p k l m >Rm)
   The matrix elements of Ho are as follows:
< 1 k l m | Ho | 1 k l m > = (Clong k)2
< 2 k l m | Ho | 2 k l m > = (Ctrans k)2
< 3 k l m | Ho | 3 k l m > = (Ctrans k)2
< p' k' l' m' | Ho | p k l m > = 0   if   p'≠p or k'≠k or l'≠l or m'≠m.
   The matrix elements of H' are calculated numerically using a computer program (in Turbo C++ pert12.cpp).
   The following "selection rules" have been deduced by numerically evaluating many matrix elements.:
(a) p' = p
(b) k' = k
(c) l' = l or l ± 2 or l ± 4 or l ± 6, etc.
(d) m' = m or m ± 2 or m ± 4 or m ± 6, etc.
   These selection rules can at least partly be explained in terms of symmetry. The Hamiltonian H' has cubic symmetry. Therefore matrix elements can only exist between states of related symmetry.
   The ultimate problem of finding normal modes of the elastic medium is now reduced formally to diagonalizing a matrix to obtain eigenvectors and eigenvalues.
   Because of the selection rules, this matrix is block diagonal. Therefore, the problem of diagonalization can be broken down into separate tasks. Unfortunately, these blocks are still infinitely large.
   The following values for elastic constants are assumed for silicon: C11 = 166.0 GPa, C12 = 64.0 GPa, C44 = 80.0 GPa, ρ = 2300 kg/m3. These values are close to those from ioffe.rssi.ru for Si at 300 K.

Table I.  Matrix elements of silicon H' for longitudinal spheroidal states connected to | 1 k 0 0 > for l≤8.
All entries have been divided by 106 k2, so as to have units of (km/s)2.
Calculated and HTML table generated using C++ program pert13.cpp     100000 iterations     Rm = 300
Data file containing numbers in this table: 37table1.txt
  1 k 0 0 1 k 2 0 1 k 4 0 1 k 4 4 1 k 6 0 1 k 6 4 1 k 8 0 1 k 8 4 1 k 8 8
1 k 0 0 82.6 0 -24.7 -20.5 -3.4 7.2 8.7 4.0 7.0
1 k 2 0 0 79.6 13.0 -15.6 -14.0 -4.4 -5.7 1.6 6.6
1 k 4 0 -24.6 13.5 82.0 -1.3 10.4 -16.3 -19.4 -5.5 2.1
1 k 4 4 -20.6 -16.2 -1.9 82.0 3.7 -15.0 1.3 -4.5 -21.2
1 k 6 0 -3.2 -13.9 9.6 3.8 81.6 -13.5 10.0 -13.2 0
1 k 6 4 6.8 -4.4 -16.4 -14.9 -13.4 84.4 0 -5.0 -2.2
1 k 8 0 8.1 -5.6 -19.8 0.8 9.5 0.3 81.9 -11.0 -0.5
1 k 8 4 4.0 1.3 -5.8 -4.4 -12.8 -5.1 -11.0 83.5 1.3
1 k 8 8 7.3 6.7 2.2 -21.7 -0.7 -2.9 -0.2 1.4 80.7

Table II.  Matrix elements of silicon H' for longitudinal spheroidal states for l≤2.
All entries have been divided by 106 k2, so as to have units of (km/s)2.
Calculated and HTML table generated using C++ program pert13b.cpp     300000 iterations     Rm = 300
Data file containing numbers in this table: 37table2.txt
  1 k 0 0 1 k 1 -1 1 k 1 0 1 k 1 1 1 k 2 -2 1 k 2 -1 1 k 2 0 1 k 2 1 1 k 2 2
1 k 0 0 82.6 0 0 0 0 0 0.2 0 0
1 k 1 -1 0 82.6 -0.3 0 0 0 0 0 0
1 k 1 0 0 0 82.6 0 0 0 0 0 0
1 k 1 1 0.2 0 0 82.6 0 0 0 0 0
1 k 2 -2 0 0 0 0 84.0 0 0 0 0
1 k 2 -1 0 0 0 0 0 84.0 -0.4 0 0
1 k 2 0 0 0 0 0 0 0 79.6 0 0
1 k 2 1 0 0 0.2 0 -0.4 0 0 84.0 0
1 k 2 2 0 0 0 0 0 0 0 0.2 79.6

Table III.  Selected diagonal matrix elements of silicon H'
Calculated and HTML table generated using pert14.cpp    
300000 iterations     Rm = 300
Data file containing numbers in this table: 37table3.txt
p k l m <pklm|H'|pklm> sound
speed
1 k 0 0 82.6 (km/s)2 k2 9092 m/s
1 k 2 -2 84.0 (km/s)2 k2 9167 m/s
1 k 2 -1 84.0 (km/s)2 k2 9167 m/s
1 k 2 0 79.6 (km/s)2 k2 8922 m/s
1 k 2 1 84.0 (km/s)2 k2 9167 m/s
1 k 2 2 79.5 (km/s)2 k2 8921 m/s
2 k 1 -1 28.6 (km/s)2 k2 5348 m/s
2 k 1 0 28.6 (km/s)2 k2 5349 m/s
2 k 1 1 28.6 (km/s)2 k2 5349 m/s
2 k 2 -2 27.3 (km/s)2 k2 5227 m/s
2 k 2 -1 27.2 (km/s)2 k2 5224 m/s
2 k 2 0 30.1 (km/s)2 k2 5491 m/s
2 k 2 1 27.2 (km/s)2 k2 5224 m/s
2 k 2 2 30.1 (km/s)2 k2 5492 m/s
3 k 2 -2 32.0 (km/s)2 k2 5659 m/s
3 k 2 -1 31.1 (km/s)2 k2 5576 m/s
3 k 2 0 25.6 (km/s)2 k2 5068 m/s
3 k 2 1 31.4 (km/s)2 k2 5610 m/s
3 k 2 2 25.6 (km/s)2 k2 5067 m/s


Table IV:  Speeds of Sound in Silicon at 300 K
      from cc3bmod.cpp
crystal
direction		ClongCtrans
[1 0 0 ]8440 m/s5859 m/s
[1 1 0 ]9148 m/s5859 m/s, 4678 m/s
[1 1 1 ]9372 m/s5102 m/s

l-Cutoff Approach
   A natural approach to solving this kind of problem is to keep only part of the matrix up to some maximum l, called lcutoff. This results in a finite size matrix which can be diagonalized. The procedure can be repeated for larger and larger values of the l cutoff. If the final results converge, then the procedure is successful.
   The simplest version is to choose an l cutoff of 3. In this case the matrix is already diagonal.
   If the l cutoff is increased to 5 then we have a 4×4 matrix to diagonalize.
   The next interesting task is to apply these new basis functions in order to find normal modes of a vibrating sphere with free boundary conditions. For each mode of p = 1 we need a corresponding mode of p = 3 in order to be able to satisfy the boundary conditions, with the exception of the | 1 k 0 0 > mode. Suppose we chose an l cutoff of 5. In that case we will have four (though approximate) new basis functions. We adjust the k of each function so that all have the same common frequency, ω. The value of ω will later be chosen to satisfy the boundary conditions. Only three of the functions have a θ component in their surface stress σrθ. Therefore σrθ must be set to zero at three unrelated places. There are four different functions contributing to σrr. Therefore, σrr must be set to zero at four unrelated, arbitrarily chosen places. This provides a total of seven equations to satisfy. Four new (approximate) basis functions have been made here, and three new approximate basis functions need to be put together from | 3 k 2 0 >, | 3 k 4 0 > and | 3 k 4 4 >. Thus, there are seven parameters to adjust. Since the linear system is homogeneous, there is always a zero solution. Looking for zeros in the 7×7 determinant will show that for certain frequencies there are also nonzero solutions. These correspond to the normal mode frequencies.
   The results can be compared to frequency analysis of direct simulations of spheres with the same elastic constants. This is a test of correctness as well as a demonstration of the rate of convergence as the l cutoff increases.

Lowest Order Approximation Applied to a Free Silicon Sphere
   Some time ago [md15.htm] I did numerical molecular dynamics simulations of spherical clusters of up to 8000 point particles that allowed me to estimate the vibrational mode frequencies of a free sphere with the anisotropic elastic properties of silicon. So there is an opportunity to test the results of this present calculation against the earlier ones.
    Because the speed of sound depends on the direction in silicon, it no longer makes sense to attempt to use a dimensionless frequency variable, like S = ν d / Clong or η = ω R / Ctrans. So I calculate all frequencies for a standard 10 nm diameter sphere of silicon. The frequencies simply scale as the reciprocal of the diameter.
    For the l = 0 spheroidal "breathing mode" case, if we choose lcutoff = 1 and thus ignore the mixing that actually happens with higher l states, the approximate eigenstate is | 1 k 0 0 >, where k = ω/(9092 m/s), since the longitudinal speed of sound to use (from Table III) is 9092 m/s. The boundary condition is zero stress σrr at the surface of the sphere. Using cubic73.cpp, the frequency is 24.2 cm-1 while the simulations in md15.htm gave 23.8 cm-1. This is a 2% disagreement. It is also interesting to compare this with traditional Lamb calculations made assuming the silicon is isotropic, getting errors of -23% to +11% depending on which crystal axis is selected for the speed of sound. The worst results are obtained if the [100] speeds of sound are assumed.

Table V:  l = 0 spheroidal "breathing mode"
frequency (cm-1)
for d=10 nm Si sphere
methodsound speeds used
23.8molecular
dynamics		correct anisotropic
24.2lcutoff = 19092 m/s (Ct not needed)
19.3isotropic8440 m/s 5859 m/s    [100] axis
23.0isotropic9148 m/s 5859 m/s    [110] axis
26.4isotropic9148 m/s 4678 m/s    [110] axis
26.3isotropic9372 m/s 5102 m/s    [111] axis

    Without anisotropy, the l = 2 spheroidal mode is 5-fold degenerate. The anisotropy of the silicon breaks the degeneracy. The frequencies calculated using cubic74c.cpp are shown in Table VI below:

Table VI:  l = 2 spheroidal "football mode"
mν (cm-1)
-2 16.0
-1 16.0
0 14.0
1 15.9
2 14.1

    The molecular dynamics simulation md15.htm showed that the l = 2 spheroidal mode was split into two frequencies 13.19 cm-1 and 16.30 cm-1. So there is 6% error in one case and 2% error in the other case.
    Initially I thought it was an error that led to m = -2 and m = +2 having different frequencies. But this really does happen.
   The l = 1 spheroidal mode is 3-fold degenerate without anisotropy. The anistropy does not break the degeneracy. There is still a single frequency. (cubic74d.cpp) The present calculation gives 19.3 cm-1 but molecular dynamics simulations got 18.41 cm-1, a 5% difference.
    For the l = 2 torsional mode, there is 5-fold degeneracy without anisotropy. The anisotropy breaks this degeneracy. (cubic77.cpp) The frequencies are

Table VII:  l = 2 torsional "wet dog mode"
mν (cm-1)
-2 14.1
-1 14.2
0 14.7
1 14.2
2 14.7

   This is qualitatively in agreement with md15.htm, but the differences are 7% and 5%.

Lowest Order Approximation Applied to a Silicon Sphere in a Matrix
   For the situation of a sphere surrounded by an infinite elastic "matrix", complex valued mode frequencies may be obtained for the decaying modes [L. Saviot]. The imaginary part Im(ω) of the complex frequency ω corresponds to decay of the mode as e-Im(ω) t or e-t/τ where τ is the damping time of the mode.
   For the lowest l = 2 spheroidal mode of a 10 nm silicon sphere in SiO2 matrix, the frequencies (calculated using cubic74e.cpp) are shown in Table VIII below:

Table VIII:  l = 2 spheroidal mode in matrix
mRe(ν) (cm-1)τ = 1/Im(ω)sound speeds used
lcutoff = 3   (Calculated using cubic74e.cpp)
-2 11.6 0.47 ps 9167 m/s 5659 m/s
-1 11.6 0.47 ps 9167 m/s 5576 m/s
 0 11.6 0.45 ps8922 m/s 5068 m/s
+1 11.6 0.48 ps9167 m/s 5610 m/s
+2 11.6 0.46 ps8921 m/s 5067 m/s
Isotropic method (Also calculated using cubic74e.cpp)
 12.90.42 ps8440 m/s 5859 m/s [100] axis
 12.20.46 ps9148 m/s 5859 m/s [110]a axis
 11.70.46 ps9148 m/s 4678 m/s [110]b axis
 11.70.46 ps9372 m/s 5102 m/s [111] axis

References:

1. Zhen Ye, "On the Low Frequency Elastic Response of a Spherical Particle" Chinese Journal of Physics Vol. 38, pages 103-110, (2000). link to pdf ::

2. A similar approach, applied to polar optical modes, can be found in: "A perturbative treatment of crystal anisotropy" pages 41-44 in "Polar optical modes in heterostructures" chapter three, pages 21-50 in "Long wave polar modes in semiconductor heterostructures" C. Trallero-Giner, R. Pérez-Alvarez and F. García-Moliner Pergamon, New York, 1998

Daniel Murray
Associate Professor
Math, Stats & Physics Unit
University of British Columbia - Okanagan
Kelowna, BC, Canada
daniel "dot" murray "at" ubc "dot" ca
For a list of related articles click here.

"Microscopic calculations on Raman scattering from acoustic phonons confined in Si nanocrystals" Jian Zi, Kaiming Zhang, Xide Xie Phys. Rev. B, Aug. 15, 1998 arXiv:cond-mat/9807293 v1 22 Jul 1998 arxiv.org/pdf/cond-mat/9807293 - c:\cofrest\9807293.pdf [email protected] http://www.fudan.edu.cn/index_ch.php
Ding Haojiang and Chen Weiqiu, "Nonaxisymmetric free vibrations of a spherically isotropic spherical shell embedded in an elastic medium"
Int. J. Solids Structures Vol. 33, No. 18, pp. 2575-2590, 1996 ::
M. Montagna and R. Dusi, "Raman scattering from small spherical particles" Phys. Rev. B 52, 10080 (1995) - about matrix effect, small influences.
"Raman scattering from fractals. Simulation on large structures by the method of moments" G.Viliani, R.Dell'Anna, O.Pilla, M.Montagna, G.Ruoco, G.Signorelli, V.Mazzacurati http://arxiv.org/pdf/cond-mat/9504039 ::
Transverse acoustic nature of the excess of vibrational states in vitreous silica http://arxiv.org/pdf/cond-mat/0209519 ::
Lingjun Wang, Guanghong Wei, Jian Zi "A planar force-constant model for phonons in wurtzite GaN and AlN: Application to hexagonal GaN/AlN superlattices" http://arxiv.org/pdf/cond-mat/9812293 ::
Raman scattering by electron-hole excitations in silver nanocrystals Authors: H. Portales, E. Duval, L. Saviot, M. Fujii, M. Sumitomo, S. Hayashi To be pub.in PRB cond-mat/0101471 http://arxiv.org/pdf/cond-mat/0101471 :: "Phonons in a nanoparticle mechanically connected to a substrate" K. Patton and M. Geller (Dec. 30, 2002 version) ::
V. L. Gurevich and H. R. Schober, Phys. Rev. B 57, 11295 (1998).
Sergey Bastrukov and Pik-Yin Lai, "Oscillatory behavior of a flexible cavity in an elastic medium" J. Moscow Phys. Soc. Vol. 9 (1999) pages 71-76 ::
Hosted by www.Geocities.ws

1