Materials Science and Engineering A 448 (2007) 221–228
Phenomenological synthesis of nanometric patterning in ballistically formed MBE silicon between 450 and 713 K J.S. Kirkaldy a , A. Perovic a,∗ , D.D. Perovic b , A. Duft a , G.C. Weatherly a a
b
Brockhouse Institute for Materials Research, McMaster University, Hamilton, Ontario, Canada L8S 4M1 Department of Metallurgy and Materials Science, University of Toronto, Toronto, Ontario, Canada M5S 1A4 Received 12 July 2006; accepted 4 October 2006
Abstract The comprehensive 673 K steady-state data set of Perovic et al. (D.D. Perovic, G.C. Weatherly, J.-P. Noel, D.C. Houghton, J. Vac. Tech. B9 (1991) 2034 and D.D. Perovic, G.C. Weatherly, P.J. Simpson, P.J. Schulz, T.E. Jackman, G.C. Aers, J.-P. Noel, D.C. Houghton, Phys. Rev. B43 (1991) 14257) was later expanded to substrate temperature variations from 450 to 687 K with transients up to the steady state. Our particular experimental emphasis on temperature dependence and atomic force microscopy (AFM) at the nanostructure scale is combined and organized so as to exhibit complete systematics as a function of temperature. We show that the kinetics and instability patterns map with consistency to the steady-state 3D microstructures observed and predicted in binary cellular forced velocity solidification of metallic alloys. Adapting the solidification phenomenological algorithm to the MBE case, we are able to extract a plausible and useful Arrhenius expression for the diffusion of occluded vacant sites in ballistically formed Si. It is further illustrated that the MBE instability margin maps in some detail to the interface Rayleigh droplet-forming morphology of 3D binary solidification. We offer this example of universality in kinetics as a contribution to the non-linear science of pattern formation, as well as complementing the often preferred analyses within a molecular dynamics paradigm. © 2006 Elsevier B.V. All rights reserved. Keywords: Molecular beam epitaxy; Silicon; Microstructure; Rayleigh instability
1. Introduction With the earlier observation of Fig. 1a and b of crystalline MBE Si growth on a substrate temperature of 673 K the investigators first conceived of an analogy with a forced velocity eutectic [1,2]. It later became apparent that a better microstructural comparison pertained to the less constrained, 3D forced velocity binary alloy cells in both elevation and plan which, except for absolute scale and rates of formation, are in qualitative correspondence (Figs. 1 and 2). The observations of Fig. 1a and b in exhibiting spontaneous deflection of ballistically implanted free volume by bulk or phase boundary diffusion offer a venue for establishing a better understanding of diffusion in ballistically formed Si and the pattern-forming process generally. Fig. 1b is unique in the literature for it provides a direct observation of ballistically created free volume production condensing as pores.
∗
Corresponding author. Tel.: +1 905 525 9140x24980; fax: +1 905 521 2773. E-mail address:
[email protected] (A. Perovic).
0921-5093/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.msea.2006.10.017
Noteworthy to the comparison of these forced velocity irreversible processes [1–12] were the precise {1 0 0} habits (Fig. 1a and b), the indications for actual and possible pinching-off of phase boundary intrusions (Figs. 1 and 2) by Rayleigh instabilities on boundary cylinders, the creation of aligned 3D strings of voids as a pure “lattice of vacancies” (Fig. 1a and b) and the 3D, close to hexagonal arrays in plan (Figs. 1b, 2 and 3). The well-known decanted liquid pattern of Figs. 2 and 3 is known as “pox” [4,7] and is characteristic of the near-margin in 3D solidification [7]. The high temperature trend towards (1 1 1) frontal facets (Fig. 1) is also seen in the steady-state solidification of both alloyed Ge and Si [7]. Subsequent to Refs. [1,2], the observations were extended to the substrate temperature range 450–713 K incorporating a transition from steady-state crystalline to amorphous Si at approximately 573 K which obviously complicates both an atomistic and phenomenological representation. The systematics defined by Figs. 4–12 will occupy our attention later. In the isothermal analysis, we assume that following frontal relaxation of a supersaturation induced instability the Si deposited ballistically on a substrate injects a steady state of free
222
J.S. Kirkaldy et al. / Materials Science and Engineering A 448 (2007) 221–228
Fig. 1. (a) Cross-section of MBE generated silicon on (1 0 0) Si substrate at 673 K where free volume precipitation produces pores in a patterned array. During production the pores are continuous with cell intrusions or pox [1,2], frontal projections have near (1 1 1) facets. (b) Pox array of (a) in plan view; spacing ∼25 nm. Note a tendency for pox alignment which anticipates morphology changes (cf. Fig. 3).
Fig. 2. Interfacial breakdown of the planar surface of 3D Al–Zn alloy. (a) Transverse section at arrow in (b) longitudinal section. Note the important innovation over earlier experiments that the liquid in the untransformed pattern at the margin has been transformed out by a global quench [3].
volume or vacant sites, which subject to capillarity, thermokinetic optimization and symmetry must migrate by volume or surface diffusion towards the vacuum of the pores. On this basis, we seek a theoretical mapping to the similar structures that appear in well characterized 2D binary forced velocity
Fig. 3. Decanted interface in marginal forced solidification of a dilute Pb–Sn alloy [4]. Dimple or pox spacing ∼100 nm. Although not to our knowledge reported this 3D manifestation is certainly subject to immediate sub-surface pinching off of the cylindrical liquid intrusion (compare Figs. 1 and 2).
cellular solidification due to constitutional supercooling (CS) (Figs. 2 and 3). Two very early molecular dynamics models of MBE Si exhibit conserved well defined occluded vacant sites [11] and non-ordered pores [12] which follow from the directed bonding of the structure. We will verify the mapping in Section 2 by testing the scaling formula in approximation for cellular solidification against that for cellular MBE Si and its observations. Restricting the ensuing scaling comparison to the margin of the latter instability on the circumstantial grounds of the microstructure comparisons we ultimately obtain an expression for the vacancy diffusion length (DS /vS ) in terms of recorded parameters. With such an unstable excess of over-saturated surface or volume vacancies considered as the solute, anticipating diffusive segregation to the pores, this configuration is the analogue of dispersed constitutional supercooling which enters the alloy solidification problem, necessarily recognizing the antisymmetry of this MBE problem since the main diffusion effect is in the near-crystalline rather than the viscous “liquid” state. While this change in symmetry does not alter the mutual dissipative character of the processes and the deflected diffusion trajectories there is a change in the symmetry of the thermal systematics. This semi-quantitative phenomenological modelling serves to focus the popular molecular dynamics approaches towards the same problem [11,12] within the conception of ballistic formation of a pseudo-alloy of vacant sites in or on Si, including the all-important inclusion of surface tension effects.
J.S. Kirkaldy et al. / Materials Science and Engineering A 448 (2007) 221–228
223
Fig. 4. MBE Si at 489 K: (a) bright-field electron micrograph demonstrating the presence of a thick amorphous layer with aligned pores and thin highly defective crystalline layer next to the substrate, (b) the free surface of the amorphous layer indicating rounded fronts, (c and d) SAD patterns taken from the substrate displaying defective crystalline layer and amorphous layer, respectively.
Fig. 5. MBE Si at 608 K: (a) the pores and (b) a highly defective layer.
2. MBE diagnostics based upon a scaling formula for alloy solidification cells The observations of Figs. 4–8 have revealed that the Si deposited is primarily amorphous changing to crystalline above about 573 K. Notwithstanding, a substantial amorphous film persists at the ballistically impacted front extending to the higher temperatures. Such an amorphous front exhibits a welldefined surface tension and a supersaturation of vacancies to
some depth which is analogous to constitutional supercooling in the solidification case, thus abetting an instability proceeding by lateral isothermal diffusion and informing us that the empirical diffusion-initiated parameter is to be identified as the pore volume fraction. Then, provided the partition coefficient k in solidification is sufficiently less than unity (typically <0.2), the experimentally tested two-dimensional low solubility, phase diagram independent scaling formula for 2D (approximating to 3D) spacing based upon a minimum undercooling
Fig. 6. MBE Si at 647 K demonstrating: (a) the layer with much reduced defects and (b) the free surface develops the faceted, predominantly {1 1 1} planes.
224
J.S. Kirkaldy et al. / Materials Science and Engineering A 448 (2007) 221–228
Fig. 7. MBE Si at 687 K: defect-free layer with {1 1 1} and {1 0 0} facets formed at the free surface with aligned pores.
Fig. 8. MBE Si at 713 K demonstrating the presence of an almost defect-free layer.
Fig. 9. AFM images of Si grown at: (a) 489, (b) 647 and (c) 687 K, respectively.
J.S. Kirkaldy et al. / Materials Science and Engineering A 448 (2007) 221–228
225
[7,10]. We record the CS criterion β∼ =
Fig. 10. Graph of pore spacing vs. temperature extracted from both TEM and AFM data. Squares represent AFM data.
Fig. 11. Graph showing the thickness of the MBE Si vs. temperature.
criterion [9] 3 1/2 λL
XL
γ T m DL ∼ = 24 L GL vL
(1)
gives an estimate of the solidification wavelength λL which is within 20% of the observed value at the margin. In this, XL is the dilute alloy mole fraction, vL the front velocity, Tm the solvent melting temperature, GL the temperature gradient, γ the surface tension, L the latent heat per unit volume and DL is the liquid diffusion coefficient. Also, Tm /GL is a characteristic thermal length, γ/L the capillary length and DL /vL is the diffusion length. The spacing magnitude is close to the geometric mean of these three characteristic lengths, and may be taken as universal provided we can enter the appropriate lengths for the MBE Si case of Fig. 1. To complete a semi-quantitative mapping between the two processes, we need one more constraint which we now identify as the experimentally tested marginal instability criterion for solidification cells within an accuracy of ±20% at the margin
(XL (1 − k)|m|)/kGL ≥1 DL /vL
(2)
Here m is the liquidus slope and k ∼ 0.1 is the partition coefficient. For dilute XL ∼ k the usual eutectic phase diagram implies |m| ∼ Tm , so near the margin the trailing product in (1) approximates to (DL /vL )2 at the margin. We further approximate that the capillarity lengths in both cases are of equal magnitude since these quantities generally tend to approach the lattice parameter in magnitude (∼0.3 × 10−9 m). Hence, presuming that the process of Si spatial rearrangement can be represented in terms of supersaturated vacancy diffusion from the analogue of constitutional supercooling, we form the unit invariant ratio from (1) eliminating the common constant 24 to yield the scaling relation between diffusion lengths DS XS 1/4 λS 3/2 DL (3) = vS XL λL vL approximately valid for mapping at the margin. It remains to identify the corresponding onset condition for the MBE cellular instability which at fixed fluence and time turns out to be a critical high temperature where the pox of Fig. 1 vanishes with the increasing temperature (Fig. 8). This means that increasing GL in solidification maps to increasing substrate temperature in the MBE process. Then, for the widely studied dilute organic alloy Succinonitrile–Salol (S–S) [5], typical near-marginal λL = 2 × 10−4 m, XL = 0.0025, k = 0.2, DL = 0.4 × 10−9 m2 /s, vL = 2 × 10−6 m/s and GL = 10,000 K/m. We obtain from (3) the near-marginal MBE thermal or diffusion length of 0.4 × 10−9 m and near-marginal DS ∼ 0.3 × 10−18 m2 /s corresponding to observed near-marginal XS ∼ 0.01 at 673 K and vS = 0.6 nm/s from the original data set [1,2]. It must be noted at this point that due to marginal faceting the instability seems to converge to an infinite rather than finite wavelength as in solidification. While it has not been directly proven theoretically that there exists an MBE Si analogue of β = 1 we here have experimental proof that near 713 K, the inherently non-steady initial condition spontaneously resolves into a frontal steady state of unique planar morphology at the margin (Fig. 8). A further careful examination of Fig. 1a indicates that the frontal intrusions and spacings are nearly equal, which is also the case for the marginal finite amplitude and unit aspect ratios representing the 2D analogue of the 3D pox of Fig. 3 [3]. The near-hexagonal distribution of the intrusions or pox in the plan of Figs. 1b and 2 compared with non-marginal solidification hexagonals according to Tiller’s Fig. 5.5 [7] confirm that in terms of relative magnitudes and shapes, a mapped marginality is a valid constraint. 3. Temperature dependence of the transient between 450 and 713 K
Fig. 12. Graph of 2 ln S (cell spacing) vs. temperature demonstrating the Arrhenius behaviour with activation energy of 9500 cal/mol.
At the lower temperatures, 450–529 K, the MBE Si consists of a transition layer comprised of stacking faults, nano-twins
226
J.S. Kirkaldy et al. / Materials Science and Engineering A 448 (2007) 221–228
and amorphous patches beyond which a final homogeneous amorphous state continues to grow as shown by the brightfield electron micrograph in Fig. 4a and the associated selected area diffraction patterns, Fig. 4b–d, taken from the substrate, defective and amorphous region, respectively. The selected area diffraction pattern from the faulted layer (Fig. 4c) shows the extra spots due to micro-twinning on {1 1 1} planes whereas the diffuse rings in Fig. 4d indicate the amorphous state. The interface of crystalline–amorphous transition of Fig. 4a is very uneven. At higher temperatures (608 K), the amorphous state gradually disappears leaving defective Si only (Fig. 5a). The defects become less dense as the temperature is raised further to 647 K (Fig. 6a). At the highest temperatures, 687 and 713 K, defect-free Si is formed (Figs. 7 and 8). Our experimental observation of the general microstructure of the films deposited in the temperature range between 450 and 713 K agrees well with the findings of other studies reported in the literature [13]. Besides the high temperature features described above, we must turn our attention to the pores observed in the whole temperature range studied, 450–713 K, and extending throughout both the amorphous stage and crystalline MBE Si growth. Apart from the early study by Perovic et al. [1,2], this feature seems to have attracted little attention in the literature. In the following, we will focus on the size and shape of these pores clearly seen in Fig. 4a where the frontal amorphous surface is a smoothly curved cellular interface (Fig. 4e). The cylindrical intrusions are broken to elongated pores due to the Rayleigh instability (Fig. 5a) and the mainly curved frontal interface persists at higher temperatures where except for a thin amorphous surface layer, almost complete crystallinity near the front prevails (Fig. 5b). At still higher temperatures, the pores are broken into rows of nearspherical voids (Fig. 6a and b) while the frontal surface (Fig. 6b) assumes a profile with predominantly imperfect {1 1 1} facets attending the pore entries. At 687 K, the spheroidized pores are more regularly spaced (Fig. 7a) whereas the frontal interface (Fig. 7b) shows a tendency to change from the predominantly {1 1 1} facets persisting at lower temperatures to the development of broad imperfect {1 0 0} facets with shorter {1 1 1} facets entering the substrate pores. Finally, as noted above, the 713 K grown Si shown in Fig. 8 has become an almost perfect single crystal with pores aligned at a much larger spacing than in the sample grown at 687 K, broken into regularly spaced spherical droplets and evidently about to disappear. We believe this signals the free volume margin where supersaturation and associated instability vanishes. In other words, this sudden change heralds the onset of a planar {1 0 0} perfect crystal interface at approximately 723 K allowing perfect {1 0 0} faceting. This is when the incoming kinetic energy no longer sufficiently perturbs the thermal spectrum of surface atoms, thus suppressing the vacancy supersaturation-induced instability. This secondary element of mapping of the margin is a strong validation of our mapping exercise in illustrating the universality in pattern formation. We initially attempted to quantify the cell spacing in order to establish the temperature dependence. However, this spacing of low temperature samples could only be roughly estimated from the electron micrographs due to the overlapping of the pores
through the thickness of the foil. In order to get a better quantification, we turned to the atomic force microscopy technique (AFM), which rendered clear images of the surfaces. For example, Fig. 9a–c shows the AFM images of Si grown at 489, 648 and 687 K. The cell spacing (S) extracted from AFM data at given temperature up to 687 K is plotted in Fig. 10. It is seen that the spacing decreases as the temperature is decreased. This behaviour is to be attributed entirely to the Arrhenius thermal effect, since amorphization has a negligible volumetric effect. From our micrographic measurements which were in agreement with the sample preparation data provided by Dr. J.-M. Baribeau from National Research Council of Canada, we have constructed the systematics of total MBE Si thickness which is almost constant (Fig. 11) and an Arrhenius plot of cell spacings S via the dimensional DS ∼ S2 /t at fixed fluence and time (Fig. 12). Fig. 12 is a plot of 2 ln S versus 1/T resulting in an Arrhenius behaviour with activation energy ∼9500 cal/mol. This value is not inconsistent with a liquid-like volume diffusion mode for vacancies in the amorphous state on lowering the temperature. The Arrhenius expression is then obtainable employing DS at 673 K, viz., −9500 DS (m2 /s) = 3.44 × 10−15 exp (4) RT assuring us that we have identified he correct number of constraints. At temperatures higher than about 687 K, near-surface lateral diffusion must prevail in a very thin amorphous layer, being ultimately blocked by faceting. From (4) we can derive the ∗ ∼X D ∼ Si self-diffusion coefficient as approximately DSi = S S −20 2 10 m /s. This is of comparable magnitude to chemical diffusion coefficients in amorphous Si at 673 K (Ref. [14], Figs. 3 and 5). We will next explore the significance of the Rayleigh instability which incorporates an activation energy for the pore wall free volume surface diffusion. 4. Further diagnostics based upon a Rayleigh instability in cellular MBE silicon The instability provides a powerful lesson on the delicate role of surfaces in phase transformations. It will be noted that the funnel-shaped, impact-generated free volume or vacancy rich poxed surface displays a sharp increase in curvature as it penetrates so is poised to act against surface diffusion delivery of vacant sites to form pores. Accordingly, the curved surface must also act to change the surface concentration distribution sufficiently to define a countering downhill diffusion profile emerging at the pore wall. But the walls themselves are unstable to spontaneous pinching by reduction of radius at a sufficient distance from the deepening stable tip of the cylinder initiating the periodic process of Fig. 8. The instability associated with pox evidently occurs in solidification near the 3D margin (Fig. 2 in comparison with present Fig. 1a and b) but originally this was not quantified. Unlike the stable MBE pores, evidence always vanished with completion of solidification. Fortunately, Fig. 2 provides the exception which proves the point. The instability also occurs as an artifact for
J.S. Kirkaldy et al. / Materials Science and Engineering A 448 (2007) 221–228
high aspect ratio cells in 2D thin film organic alloy solidification where it is associated with weak wetting of the confining glass plates by the solid so as to create the necessary cylindrical roots. We juxtapose the foregoing MBE Si bulk diffusion model to the proposition that vacancies depositing into upper cylindrical pore walls define a high density of dangling bonds which initiate a high free volume surface diffusion rate which in turn activates the Rayleigh instability. The following analysis explores the quantitative elements of this phenomenon providing a comparison insight into the invariable Rayleigh root instabilities in 3D solidification. Nichols and Mullins [15] (see also Ref. [16]) have offered a quantitative model of the surface self-diffusive and capillary controlled ripening of a cylinder. The time to formation of a sphere near the end of a half-cylinder can be approximated as τ=
6.56r04 kT DB γΩ4/3
(5)
where for our conditions the solid-vacuum surface energy (γ) is estimated as 1.5 J/m2 , the atomic volume (Ω) is 10−29 m3 and kT is 9.28 × 10−21 J for a growth temperature of 673 K. From Fig. 1 with a front velocity of 6 × 10−10 m/s and a radius (r0 ) of the cylinder of 3 × 10−9 m, τ can be estimated from the delay distance for bridging of the voids as 100 s. Substituting these values in Eq. (5) gives surface DB ∼ 1.6 × 10−14 m2 /s (see below). The value is plausible relative to inferred surface structure and the much smaller bulk DS value of ∼10−18 m2 /s deduced above. Furthermore, since pertaining to an observably smooth pore surface at a temperature near 673 K it is, as expected, in turn significantly larger than pre-occlusion liquid-like surface values for ballistically perturbed planar frontal atoms (10−12 m2 /s) estimated by Gilmer and Roland [11]. The latter large value undoubtedly played a part in the transient instability response to a 3D sinusoidal perturbation [7] which led to the 3D stable pox array. Further to the formation of the microstructure, the vacancies are not an equilibrium solute so the pores cannot homogenize out as the Rayleigh droplets associated with Figs. 2b and 3 must do during uninterrupted solidification. We must also emphasize that the foregoing calculations are valid only as to magnitudes. Notwithstanding, the alternative molecular dynamics numerics have much further to go before reaching a quantitative science. 5. Discussion It is worth repeating that all of the free volume created ballistically and absorbed at the front in a perfect steady state below the temperatures of extreme faceting must be deposited in the trailing pores. This is simply because, with the substrate lattice match preserved, steady-state free volume created cannot, beyond a negligible amount, be given up laterally from the external surfaces. The concept of vacancy diffusion and clustering in such a globally near-conserved process is therefore a sound one. The ∼1% free volume implied by the pores in Fig. 1a and b at 673 K fixes the occluded frontal vacancy concentration at about the same average bulk value in the diffusion zone and
227
is comparable to that found in the molecular dynamics experiments of Gilmer and Roland [11]. It is furthermore interesting that Smith and Strolowitz [12] have generated voids using the molecular dynamics methodology. However, their simulation, which exhibits a fractal character, refers to a process which is more like direct precipitation than the presently proposed continuous and ordered two-stage capillary pinching process. The enumeration of elements of universality at equilibrium or in kinetics is generally regarded as a desirable outcome of physical theory. The extended theme of this contribution has been that the observational equivalence between two processes that share a similar pattern can be sufficiently strong as to imply a numerical mapping between selected operating parameters. This appears to be the consequence of a patterning similarity principle in Newtonian mechanics of particle manifolds [17] which does not require dimensional analysis by power law expansions. Such modelling by mapping yields evaluations and phenomenological implications which in one or the other case are not directly accessible to experiment nor easily modelled atomistically. It furthermore provides for the first time an estimate of the bulk self-diffusion coefficient of Si in impacted MBE in the range of 10−19 m2 /s. This is not to discount the intelligence that molecular dynamics programs could be formulated to yield the vacancy diffusion effects generated herein. To recapitulate the ascendant role of universality in spontaneous pattern formation, we refer to the precise marginal microstructure of Fig. 8 in the near approach to steady-state stability. In this, the pressure due to root curvature exactly balances a thermally constructed driving force which in 3D designates the onset of the Rayleigh instability, precisely represented in full detail, not only by the MBE experiment (Fig. 7), but also by a dripping water tap near its marginal state. Acknowledgements We dedicate this paper to the memory of our colleague and mentor, George Weatherly, who passed away during a lengthy and much-argued theoretical evolution. It was George who first pointed out the solidification analogy together with the water tap analogy. We are indebted to Dr. J.-M. Baribeau of the National Research Council of Canada who provided us with the samples and the sample preparation data of Fig. 11. References [1] D.D. Perovic, G.C. Weatherly, J.-P. Noel, D.C. Houghton, J. Vac. Tech. B9 (1991) 2034. [2] D.D. Perovic, G.C. Weatherly, P.J. Simpson, P.J. Schulz, T.E. Jackman, G.C. Aers, J.-P. Noel, D.C. Houghton, Phys. Rev. B43 (1991) 14257. [3] T. Sato, K. Ito, G. Ohira, Trans. JIM 21 (1980) 441. [4] G.S. Cole, Ph.D. Dissertation, University of Toronto, 1963, in press. [5] L.X. Liu, J.S. Kirkaldy, J. Cryst. Growth 144 (1994) 335. [6] J.S. Langer, Rev. Mod. Phys. 52 (1980) 1. [7] W.A. Tiller, The Science of Crystallization: Macroscopic Phenomena, Cambridge University Press, Cambridge, 1992, p. 236. [8] P. Kurowski, S. de Cheveign´e, G. Faivre, C. Guthmann, J. Phys. Fr. 50 (1989) 3007. [9] L.X. Liu, J.S. Kirkaldy, A. Kroupa, Acta Metall. Mater. 43 (1995) 2905.
228
J.S. Kirkaldy et al. / Materials Science and Engineering A 448 (2007) 221–228
[10] W.A. Tiller, K.A. Jackson, J.W. Rutter, B. Chalmers, Acta Metall. 1 (1953) 428. [11] G.H. Gilmer, C. Roland, Appl. Phys. Lett. 65 (1994) 824. [12] R.W. Smith, D.J. Strolowitz, J. Appl. Phys. 79 (1996) 1448. [13] D.J. Eaglesham, H.-J. Gossman, M. Cerullo, Phys. Rev. Lett. 65 (1990) 1227.
[14] P.A. Stolk, S. Coffa, J.M. Poate, in: H. Pain, D. Gupta (Eds.), Diffusion in Amorphous Materials, TMS, 1994, p. 177. [15] A. Nichols, W.W. Mullins, J. Appl. Phys. 36 (1965) 1826. [16] K. Brattkus, J. Phys. Fr. 5 (1989) 2999. [17] E.T. Whittaker, A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, Cambridge University Press, Cambridge, 1952, p. 47.