Open access

Combination of Acoustic Methods and the Indentation Technique for the Measurement of Film Properties

Written By

Fabio Di Fonzo, Francisco García Ferré, Dario Gastaldi, Pasquale Vena and Marco G. Beghi

Published: 28 August 2013

DOI: 10.5772/56815

From the Edited Volume

Modeling and Measurement Methods for Acoustic Waves and for Acoustic Microdevices

Edited by Marco G. Beghi

Chapter metrics overview

2,472 Chapter Downloads

View Full Metrics

1. Introduction

New techniques are continuously being developed to produce films and thin films, which have a variety of properties and are exploited by an ever increasing number of technologies. The properties of films typically depend on the preparation process, and can be significantly different from those of the material in bulk form, due to the different micro- and nano- structures obtained by the various production processes. Although new techniques are also introduced to measure the film properties, the characterization of thin layers remains an open issue. A precise knowledge of the mechanical properties is crucial in several cases. Namely, whenever the films have structural functions, as it happens in several micro electro-mechanical systems (MEMS), and when the operation of a device is based on acoustic waves, as it happens in various types of micro-devices, for either sensing or signal processing purposes. But since the mechanical properties depend on the microstructural properties, they are also of more general interest, because they can be a good indicator also for properties which are not strictly mechanical, but for instance depend on the degree of compactness of the film.

A full mechanical characterization includes the determination of both the elastic properties, which characterize the reversible deformations, and the properties which characterize the irreversible behaviour. In most cases the elastic behaviour can be completely characterized by the elastic moduli, or equivalently by the components of the elastic tensor. It is well known that also in the simplest case, the homogeneous isotropic continuum, the elastic stiffness cannot be characterized by a single parameter, but needs two independent parameters. These can be taken as, for instance, the bulk modulus B, which characterizes the stiffness with respect to volume changes, and the shear modulus G, which characterizes the stiffness with respect to changes of shape, or can be taken (as it happens more frequently) as Young modulus E and Poisson’s ratioν. In the case of anisotropic solids the number of independent parameters increases. The inelastic behaviour is typically more complex, and only some overall parameters, like the yield stress or the hardness, is usually measured.

A wide variety of methods has been developed to perform the mechanical characterization of bulk specimens. Among them, a specific class exploits vibrations of acoustic nature as a probe of the material behaviour. These methods are non destructive, and involve only elastic strains; therefore, they are intrinsically unable to give indications about any inelastic behaviour. On the other hand, due to the complete absence of inelastic strains, the relationship between the raw measurement results and the stiffness parameters can be more straightforward, and less subjected to uncertainties or to spurious effects, possibly allowing better accuracies.

The mechanical characterization of supported films typically requires specific methods. The most widespread technique is indentation, for which a specific standard exists: instrumented indentation, ISO 14577. Indentation induces both elastic and inelastic strains, and supplies significant information about irreversible deformation, namely the hardness value, but the extraction of the information concerning the elastic behaviour is non trivial. The most frequently adopted analysis technique (Oliver & Pharr, 1992; Bhushan & Li, 2003) gives access, once the compliance of the indenter is taken into account, to a single parameter, usually referred to as “indentation modulus”. If a reasonable assumption about the value of Poisson’s ratio is available, a value of Young modulus can be derived, which obviously depends on the reliability of the adopted assumption. In the case of films, since the nano- and micro-structure can be different from that of bulk samples, a well grounded assumption about the value of Poisson’s ratio might be not available. It is also well known that, when supported films are measured, care must be exercised to avoid the influence of the substrate properties.

Methods which exploit vibrations of acoustic character have been developed also for supported films. Acoustic properties depend on stiffness and inertia; therefore, as it happens for bulky samples, acoustic methods require a value of mass density, independently measured. However in acoustic methods the intrinsic absence of inelastic strains makes the derivation of the stiffness parameters less subjected to spurious effects, and less dependent on specific modelling assumptions.

Among the techniques based on acoustic excitations, the so called laser ultrasonics techniques rely on impulsive, therefore broadband, excitation, requiring an analysis of the response in the frequency domain. Excitation is performed by a laser pulse, and the response is measured by piezoelectric or optical means. Quantitative acoustic microscopy relies instead on monochromatic excitation by a piezoelectric transducer, which typically is also exploited as sensor. In the detection of vibrational excitations, substantial advantages are offered by light, a contact-less and inertia-less probe; such advantages are particularly relevant in the measurement of films and small structures. They are exploited by Brillouin spectroscopy, which relies on Brillouin scattering: the inelastic scattering of light by acoustic excitations. Brillouin spectrometry, like Raman spectrometry, is not performed by specifically exciting vibrations, but relies on thermal excitation (anti-Stokes events) and on the excitation produced by the scattering process itself (Stokes events). In both cases excitation has a small amplitude, but it has the broadest band, allowing access to the GHz and multi GHz band, accessible by piezo-electric transducers only with specific transducer structures.

For all these methods, the outcome is the measurement of the propagation velocity of one or more acoustic modes. If sufficient information is gathered, because several acoustic modes are detected and/or their velocities are measured for a sufficiently wide interval of frequencies, a full elastic characterization can be achieved by purely vibrational means, if an independent value of the mass density is available (Comins et al., 2000; Zhang et al., 2001; Djemia et al., 2004; Beghi et al., 2011).

However, a complete elastic characterization by only acoustic means is not always achievable, as it happens when only one acoustic mode is measurable. The results of acoustic methods and of indentation can therefore be combined, with the purpose to obtain a complete elastic characterization, not achievable by each of the techniques alone. This can be particularly useful in the case of new materials or of films of unconventional structures, for which a reliable assumption about the value of Poisson’s ratio, needed by indentation, is not available. And even when a full elastic characterization by acoustic techniques is achievable, the combination with indentation offers a useful cross-check among techniques based on completely different principles and also supplies the hardness value, intrinsically out of reach of the acoustic techniques.

This chapter is devoted to this combination of indentation with acoustic techniques, namely quantitative acoustic microscopy (Bamber et al., 2001; Goodman & Derby, 2011) and Brillouin spectroscopy (Garcia Ferré et al, 2013).


2. Acoustic modes in elastic solids

Some basic elements of the acoustics of homogeneous and of layered solids are briefly recalled here for sake of completeness. For an elastic continuum undergoing deformation the mechanical state can be represented by the displacement vector fieldu(r,t), where r=(x1,x2,x3) is the position vector and t is time. The local state of the solid being represented by the strain and stress tensors, if the behaviour is linear elastic the stiffness is represented, for any type of anisotropy, by their proportionality: the tensor of the elastic constantsCijmn, which can be conveniently cast in the form of the matrix of the elastic constants Cij; inertia is represented by the mass density ρ. In a homogeneous linear elastic continuum with no body forces the equations of motion for the displacement vector field are homogeneous, and have the form (Auld, 1990; Every, 2001; Kundu, 2012)

ρ2uit2=j,m,nCijmn2umxjxn,  i=1,2,3.E1

These equations admit solutions in the form of travelling plane acoustic waves, or modes (Every, 2001; Kundu, 2012):


where k is the wavevector, ω = 2π f the circular frequency, f the frequency, A˜an arbitrary complex amplitude, and e the polarization vector, which is normalized. The elastic continuum description, underlying Eqs. (1) and (2), is appropriate whenever the wavelength λ=2π/|k| is much larger than the interatomic distances. The three translational degrees of freedom of each infinitesimal volume element correspond, for each wavevector k, to three independent modes, having different polarization vectors e1, e2, and e3, and frequencies. In general the phase velocity v=ω/|k|=λf depends on the directions of both k and e. In an infinite homogeneous medium, waves of the form of Eq. (2) exist for any frequency f compatible with the mentioned lower limit for wavelength.

In the simplest case, the isotropic continuum, which we will mainly consider, the matrix of the elastic constants is fully determined by only two independent quantities; the only non null matrix elements are C11 = C22 = C33, C44 = C55 = C66, C12 = C13 = C23 = C11 - 2C44. In this case the shear modulus G coincides with C44, while Young modulus E, Poisson’s ratio ν and bulk modulus B are respectively given by (Every, 2001; Kundu, 2012)


In an isotropic continuum the phase velocities do not depend on the direction of k, but only on the relative orientation of e with respect to k; one of the three modes is longitudinal (e1||k) and has velocityvl, the other two are transversal (e2k,e3k), are independent (e2e3) and degenerate: they have the same velocity vt (Auld, 1990; Kundu, 2012). The two velocities are

            vl=C11ρ=Eρ(1ν)(1+ν)(12ν)             (a)and     vl=C44ρ=Eρ12(1+ν)                        (b)E7

In non isotropic continua more than two independent quantities are needed to determine the matrix of the elastic constants, and the phase velocities, beside depending on the direction of k, have a more complex dependence on the Cij values. They however remain of the formv=C/ρ, where C is the appropriate combination of elastic constants and, possibly, direction cosines of k.

In a homogeneous semi-infinite solid the free surface breaks the translational symmetry in the direction perpendicular to the surface, causing new phenomena: the reflection of bulk waves, and the existence of surface acoustic waves (SAWs). Namely, at a stress free surface, Eqs. (1) admit, for an isotropic continuum, a further solution: the Rayleigh wave, the paradigm of SAW. Such waves have peculiar characters (Farnell & Adler, 1972): a displacement field confined in the neighborhood of the surface, with the amplitude which declines with depth, a wavevector k parallel to the surface, and a velocity lower than that of any bulk wave, such that the surface wave cannot couple to bulk waves, and does not lose its energy irradiating it towards the bulk. Pseudo-surface acoustic waves can also exist, which violate this last condition. The velocity vR of the Rayleigh wave cannot be given in closed form; a good approximation is (Farnell & Adler, 1972)


The continuum model of a homogeneous solid does not contain any intrinsic length scale. Accordingly, all the solutions for this model are non dispersive, meaning that the velocities (Eqs. (7) and (8)) are independent from wavelength (or frequency).

More complex modes occur in non homogeneous media. Layered media, like a supported film, are a particularly relevant case, in which new types of acoustic modes can occur (Farnell & Adler, 1972); namely, modes confined around an interface (Stoneley waves) and modes which are essentially guided by one layer (Sezawa waves). All these modes have a wavevector k parallel to the surface. In this case the medium does have an intrinsic length scale, given by the layer thicknesses, and the modes become dispersive: their velocities depend on wavelength, or, more precisely, on the wavelength to thickness ratio(s). When the wavelength is comparable to the layer thicknesses, the displacement fields of the acoustic modes extend over several layers: the modes are modes of the whole structure. Their velocities therefore depend on the properties of at least two adjacent layers, and can be computed only numerically, as non trivial functions of the properties of the substrate and the layer(s), and of the wavelength to thickness ratio. The dispersion relations ω(|k|) orv(f) can thus be obtained.

In particular, in the design and analysis of measurement methods for a supported film, two basic cases can be distinguished: ‘high’ and ‘low’ thicknesses. Thickness is ‘high’ or low’ in comparison with the acoustic wavelength involved in the measurement; therefore a same thickness can be ‘low’ in an experiment which exploits relatively low frequencies, and ‘high’ in an experiment performed at higher frequencies, i.e. smaller wavelengths.

The thickness is ‘high’ when it is significantly larger than the involved acoustic wavelength. In this case the film behaves acoustically as a semi-infinite medium: it supports acoustic waves of bulk type, both longitudinal and transverse, as if it was infinite, and, at the free surface, supports the Rayleigh surface wave, as if it was semi-infinite. In this case all these acoustic modes are non dispersive, their properties depend on the material properties of the film alone, and in the analysis of experimental results the thickness value is irrelevant. It turns out, experimentally, that this regime is achieved for thicknesses of the order of a few wavelengths, without requiring thicknesses which exceed the wavelengths by orders of magnitude.

Thicknesses are ‘low’ when they are comparable to or smaller than the acoustic wavelength. The bulk waves cannot be fully developed in the film, the Rayleigh wave becomes a so called ‘modified Rayleigh wave’, whose displacement field penetrates in the substrate, and Sezawa waves can be supported, whose displacement fields also penetrate in the substrate. Therefore the properties of all these acoustic modes depend on the material properties of the film, but also on the film thickness, whose value becomes a crucial parameter, and on the properties of the substrate. In particular, the velocities of both the modified Rayleigh wave and the Sezawa waves depend on the wavelength to thickness ratio: these modes are dispersive. This must be taken into account in the analysis of experimental results.

A gradual transition from the ‘low’ to ‘high’ thickness behaviour occurs at intermediate values. For instance, when the film thickness is very small the Rayleigh wave has a displacement field which mainly extends in the substrate, and a velocity which approaches the Rayleigh wave velocity of the substrate. For increasing thickness the fraction of the acoustic energy which is in the substrate decreases, and the velocity gradually approaches the Rayleigh wave velocity of the film material; the latter is reached when the displacement field amplitude becomes negligible before reaching the substrate.


3. Measurement techniques

Some observations are recalled here concerning two methods to measure the acoustic velocities, namely quantitative acoustic microscopy (AM) and Brillouin spectroscopy (BS), and micro-indentation. It should be recalled that, unlike polymers which exhibit a frequency dependent mechanical response, the viscous effects in ceramic-like materials and in metals is generally negligible, such that the high and low frequency properties are essentially indistinguishable. It can also be remembered that the acoustic methods measure the adiabatic stiffness, which does not coincide with the isothermal one, measured in monotonic tests (if strain rate is not too high); however, in elastic solids the difference between adiabatic and isothermal moduli seldom exceeds the measurement uncertainties (Every, 2001). It is also worth remembering that in most metals and ceramics the acoustic velocities are of the order of a few km/s = μm×GHz.

3.1. Acoustic microscopy

Acoustic microscopy (Zinin, 2001; Zinin et al., 2012) is one of the techniques which adopt harmonic excitation. It exploits a piezoelectric actuator coupled to sapphire acoustic lens, which operates also as a transducer, in a pulse-echo operation mode. The acoustic lens is mechanically coupled to the sample by a liquid drop, often distilled water. Acoustic microscopy can be operated with imaging purposes; in the quantitative acoustic microscopy version (Zinin, 2001) it allows to measure the acoustic properties of the sample.

Experiments are performed by first adjusting the lens to sample distance z to give the maximum intensity of the reflected signal: at this distance the lens is said to be at focus. The distance is then decreased (the lens is defocused), measuring the intensity of the reflected signal. Since the direct outcome of the transducer is a voltage V, the so called V(z) curve is obtained. Such a curve has characteristic oscillations, of period Δz, arising from the interference between two acoustic paths. One is the direct back reflection at the liquid/sample interface, due to the mismatch of acoustic impedances, the second involves a refraction into the solid sample: acoustic wave are excited in the solid, which propagate irradiating part of their energy towards the coupling liquid, back to the lens. In homogeneous solids, and in “high thickness” layers, the refracted bulk acoustic waves propagate towards the depth of the solid, and the energy irradiated back into the fluid mainly comes from the Rayleigh wave, propagating along the specimen surface. In specimens with thin coatings also other acoustic modes can play a significant role.

When the Rayleigh wave has the dominant role, the measurement of the period Δz allows to derive its velocityvR, since

Δz=v02f(1cosθR)   with   sinθR=v0vRE9

where v0 is the acoustic velocity in the fluid and f is the acoustic frequency. Piezoelectric transducers coupled to acoustic lenses typically operate in the frequency range from tens of MHz up to 1 GHz, resulting in acoustic wavelengths from a few to many micrometres. Therefore only films having a thickness of several micrometres can be in the ‘high thickness’ regime mentioned above.

3.2. Brillouin spectroscopy

Brillouin scattering is the inelastic scattering of light by vibrational excitations of acoustic nature, or by long wavelength acoustic phonons; the scattering mechanisms are discussed elsewhere. In BS measurements a focused laser beam, of angular frequency Ωi, wavelength λ0 and wavevectorqi, impinges on the sample, and the light scattered with wavevector qs is collected and spectrally analyzed. The strongest feature in the spectrum is due to the elastically scattered light, at frequency Ωi, but also Stokes/anti-Stokes doublets can be present, due to inelastic scattering by thermally excited vibrations (Sandercock, 1982; Comins, 2001; Grimsditch, 2001; Every, 2002; Beghi et al., 2012). The wavevectors qi and qsselect the probed wavevectors k=qsqi andk=(qsqi). Inelastic scattering by a bulk wave (see Eq. (2)) or a SAW of such wavectors gives one or more doublets, at frequencies Ωs = Ωi ± ω. Detection of a doublet thus allows to measure the velocity of an acoustic mode, either v=ω/|k| orv=ω/|k|. In the evaluation of k the effects of refraction must be taken into account, as recalled below.

BS thus measures the velocities of bulk and surface acoustic waves in a fully optical, and therefore contact-less, way. Obviously, bulk waves are detectable only in sufficiently transparent media. The probed acoustic wavelengths λ=2π/|k| and λ=2π/|k| are determined by the optical wavelength and the scattering geometry: with visible light they are typically sub-micrometric, resulting in acoustic frequencies in the GHz to tens of GHz range. Wavelengths in this range are three orders of magnitude larger than interatomic distances, making the continuum description (Eq. (1)) fully appropriate.

Although BS can be performed in various geometries (Grimsditch, 2001; Beghi et al., 2011; 2012), the most frequent one is backscattering (qs=qi), which maximizes|k|. This geometry selects, for scattering by bulk waves in isotropic media, the specific probed wavevectork=±2qin, where n is the refractive index, while, for scattering by SAWs, it can explore a range of wavevectors of modulus |k|=±2|qi|sinθwhich depend on the incidence angle θ (the angle between the incident beam and the surface normal) but do not depend on n. Therefore the analysis of experimental results for bulk scattering requires an independent value of the refractive index. Since the probed acoustic wavelengths are typically sub-micrometric, films having a thickness of a couple of micrometres already are in the ‘high thickness’ regime mentioned above.

3.3. Nanoindentation

During nano-indentation tests, an indenter is pressed onto the surface of the material that has to be characterized. The tip of the indenter, that can have different shapes, is made of a hard material, usually diamond. During the test, the load P and the tip penetration h are continuously recorded during a loading-unloading cycle, building up the P-h curve, whose analysis allows the computation of indentation modulus and indentation hardness, attributable to the elastic and post-elastic response of the material, respectively.

The key parameter when analyzing indentation is the reduced modulus Er. Generally Er depends on the elastic properties of the material as well as on the shape of the indenter and, therefore, the shape of the contact area. With the assumption of axisymmetric indenters and homogeneous, elastic, isotropic materials, Er can be computed with the simple relation:


Where E and ν are Young modulus and Poisson’s ratio of the material, respectively.

From the theory of contact mechanics between an isotropic elastic half space and an elastic sphere (Hertz, 1881), the relationship between the contact force P and the relative displacement h between the bodies dependes on the elastic properties of both bodies, through one single parameter E* according to the relationship:


in which R is the indenter radius and E*, also referred to as indentation modulus, is such that:


where Ei and νiare the Young modulus and the Poisson’s ratio of the indenter. In case of diamond, it is Ei=1141GPa andνi=0.07.

The above solution, holding for a spherical indenter, has been extended to conical indenters by Sneddon (Sneddon, 1965). By using simple geometrical arguments on the geometry of the tip, he obtained the following force-displacement relationship:


in which αis the opening half angle of the conical tip and a is the contact radius.

A simple relationship between the elastic parameter E* and the the unloading contact stiffness S=dP/dh is:


in which A is the projected contact area having radius a. The above solution holds for an elastic contact between a perfect conical surface and the elastic isotropic semi-infinite solid.

Oliver and Pharr (Oliver and Pharr 1992) introduced a generalized method to determine the reduced modulus Er for pyramidal indenters and non linear materials with irreversible response. When elastic-plastic response is obtained, the unloading branch of the P-h curve will be lower with respect to the loading branch, and will exhibit non zero displacement at zero load, denoting a residual inprint on the material surface. The method of Oliver and Pharr is based on the basic hypothesis that the initial portion of the unloading P-h curve is essentialy driven by elastic response of the material. Thus the contact stiffness S calculated at the top portion of the unloading curve can be used in a generalized version of equation (14) to estimate the indentation modulus:


in wich an adjusting empirical parameter β has been introduced with the purpose to account for specific pyramidal shapes. For the three sided pyramid (Berkovich tip) β = 1.034 is usually assumed.

The contact area A is the projected area of the material in contact with the indenter surface and it is a function of hc which is the height of the contact surface in the vertical direction (the indentation direction). The function A(hc) is generally calibrated for low depth indentations by running experiments on reference materials (typically fused silica); calibrated functions A(hc) accounts for the blunt geometry of real used tips and therefore it must be re-calibrated periodically as the number of idents made with the tip increases.

The Oliver and Pharr method provides a simple relationship for the experimental determination of the contact depth hc


which is used to calculate the contact area and the reduced modulus Er using (15) and (12).

The ability of the nanoindentation technique to probe small amount of material by shallow penetration depths enables one to probe the mechanical properties of thin coatings without removing the film from the substrate.

The main difficulty in nanoindentation of thin films is to ensure that the properties of the substrate are not affecting the characterization of the coatings. To achieve this aim, a simple rule to restrict the maximum penetration depth to no more than 1/10 of the coating thickness is usually applied.

However, for measurement of the elastic properties, the influence of the substrate compliance is unavoidable. Despite this intrinsic difficulties, methods to determine film thickness out of nanoindentation on film/substrate systems are available. As an example, Doener and Nix (Doener and Nix 1986) determined a relationship between an effective reduced modulus (which is strictly related to E*) and the properties of film, substrate and indenter, by using finite element simulations (Fischer-Cripps, 2011):


in which the subscripts f and s denote film and substrate, respectively. In (17) α is a parameter which governs the decay of the effect of the substrate for increasing penetration depths of indentations and t/Ais the geometric scaling parameter. For ideal conical indentation, A is proportional to h2, so the scaling parameter is proportional to t/h.

For deep indentations, i.e. for small t/h, equation (17) provides the properties of the substrate; whereas, for shallow indentation (i.e. at the limitt/h0)equation (17) provides the properties of the film.


4. Data analysis

Data analysis is first discussed for the case of an isotropic semi-infinite medium. This analysis applies to the ‘high thickness’ films as well, either supported or free standing; with AM this regime can be achieved only by thicknesses of many micrometres, while with BS it is already achieved by thicknesses of a couple of micrometres.

The assumption of a semi-infinite medium, or of a film of ‘high thickness’, implies that the velocities of all the acoustic modes do not depend on the wavevector (Farnell & Adler, 1972). The crude outcome of measurements is a set of NE values for the reduced modulus {Er,m} m = 1,2,...,NE, and a set of NR values for the Rayleigh velocity{vR,i}, i = 1,2,...,NR; if the bulk acoustic wave is measurable, a third set of Nl values for the bulk longitudinal velocity{vl,j}, j = 1,2,...,Nl, is also available. For each value an uncertainty can be evaluated, like σvR,i forvR,i; these uncertainties can be statistically evaluated a posteriori, and therefore be the same for all the values belonging to the same set, or can be evaluated for each single value, if a sufficiently detailed analysis of the measurement process is available. With AM only the vR values are accessible, the bulk waves being not measurable; with BS the bulk waves are measurable only for sufficiently transparent media (silicon, for example), and in this case also the bulk transversal wave can be measurable, beside the bulk longitudinal one, giving a fourth set of Nt values for the bulk transversal velocity{vt,j}. This does not always happen, because the spectral doublet due to scattering by the transversal wave is typically weaker than that due to the longitudinal wave. From the set of NE values for the reduced modulus {Er,m} an average value E¯r can be derived, with its uncertaintyσE¯, and the same can be said forv¯R, v¯landv¯t, when they are measurable.

The assumption of an isotropic medium implies that the stiffness is fully identified by two independent parameters: therefore it is intrinsically a two-dimensional quantity, which can be represented by a point in a two dimensional ‘stiffness space’. The ‘stiffness space’, in turn, can be represented by the (E,ν) couple, i.e. the (E,ν) plane, or, equivalently, by other planes, like the (B,G) plane or the (C11, C44) plane, or the plane defined by the two Lamé constants. Any of these planes can be mapped into any other (see Eqs. (3)-(6)); the representation by the (E,ν) plane is first considered here. For any point of this plane the value Er(E,ν) is immediately computed, and, if the mass density is known, the value vR(E,ν)is also computed, as well as vl(E,ν) andvt(E,ν).

4.1. Graphical method

The condition Er(E,ν)=E¯r identifies a line in the ‘stiffness space’, and the condition E¯rσE¯Er(E,ν)E¯r+σE¯ identifies a confidence band, or a ‘stripe’. The same happens forv¯R, v¯landv¯t, when available. A first analysis procedure is graphical, and basically consists of drawing the available curves, or bands, in the stiffness space. If two of them are available, typically E¯r and v¯R from AM (Bamber et al., 2001; Goodman & Derby, 2011), their intersection fully identifies the stiffness, and allows a semi-quantitative estimate of its uncertainty. If three or more lines are available, as it can happen with BS from sufficiently transparent samples (Beghi et al., 2011; Garcia Ferré et al, 2013), all their intersections typically do not coincide exactly, but the amplitude of the region which contains them allows a qualitative estimation of the consistency of the different measurement techniques. In particular, a confidence band which lies beside the intersection of the others is a hint to a possible systematic error in one of the measurement techniques.

Validations of the data analysis procedure have been performed exploiting fused silica samples supplied for calibration of the indentation equipment (Bamber et al., 2001; Garcia Ferré et al, 2013). In the derivation of stiffness values from the measured acoustic velocity the widely accepted value of 2200 kg/m3 was adopted for the mass density of fused silica. Fig. 1 shows a recent result obtained by BS, which allows to measure also the doublets due

Figure 1.

adapted from Garcia Ferré et al. 2013). Representation in the (E,ν) plane of the measurements on a bulk silica platelet of nominal properties E = 72 GPa and ν = 0.17. Results from indentation (reduced modulus Er, green lines) and from BS (bulk longitudinal wave velocity vl, blue lines; Rayleigh wave velocity vR, red lines; ratio of bulk longitudinal to bulk transversal wave velocities vl/vt, violet lines). For each quantity: central values and uncertainty bands. The conversion from raw spectral data to vl is performed exploiting the measured value of refractive index nsi = 1,464; the conversion from velocities vl and vR to the elastic moduli is performed by the widely accepted value of 2200 kg/m3. The ratio vl/vt is directly obtained from raw spectral data, and is a function of ν which can be plotted without needing any of the above data.

to bulk waves (Garcia Ferré et al, 2013), In the derivation of the corresponding velocity the value of the refractive index n is needed. The latter was measured by an ellipsometer, obtaining nsi = 1,464 ± 0.003. An analogous outcome was obtained by AM (Bamber et al., 2001), measuring only the Rayleigh velocity but not needing the value of the refractive index.

4.2. Least squares estimation method

A more detailed statistical analysis can be performed by a least squares minimization procedure. Following standard non linear estimation theory (Seber & Wild, 2003), a least square estimator S(E,ν) is built as the weighted sum over each single result:


the weights being determined by the uncertainties. Obviously, the sum over the vl values is present only if the measurement of the bulk longitudinal wave is available. If the bulk transversal wave is also measurable, a further summation over the vt values is also present in Eq. (18). The minimum Smin of S(E,ν) identifies the most probable value of the (E,ν) couple, and the isolevel curves of the normalized estimator


identify the confidence region at any predetermined confidence level (Lefeuvre et al., 1999; Beghi et al., 2001, 2011, 2013). The value of S’(E,ν) that identifies a confidence region is determined by the confidence level (typically the 68%, 90% and 95% confidence levels are considered) and by the available number of measurements (velocities of acoustic modes and indentation data) (Seber & Wild, 2003).

The estimation of stiffness is performed computing the normalized estimator S’(E,ν) at the nodes of a discrete mesh, which is refined until discretization effects become negligible, and then identifying the minimum ofS'(E,ν), as well as the confidence regions. It is often found that the precise position of the minimum, which by definition falls is in region in which the gradient of the estimator is small, is relatively more sensitive to numerical noise. The confidence regions typically have a regular shape; their boundaries, which fall in regions in which the gradient of the estimator is significant, turn out to be a more robust result.

The robustness of this outcome was exploited to define the following procedure to pick the final results of the evaluation: for each parameter, the minimum and the maximum of the values falling within the confidence region are taken. Their average value is taken as the estimate of the parameter, and the semi-amplitude of the interval is taken as the estimate of its uncertainty. Graphically, this means taking, instead of the point at which the estimator is at its minimum and of the precise shape of the confidence region, the rectangle, with sides parallel to the axes, circumscribed to the confidence region, and its centre (see Fig. 2a). This gives a simple and clearly defined algorithm to identify the measured value and its uncertainty, for each of the parameters, namely for E and ν in the representation considered until here. However, it must be remembered that these values are only the one-dimensional projections of a quantity, the stiffness, which is intrinsically two-dimensional. A statistically more rigorous procedure to take into account its two-dimensional nature would require considering it as a bivariate normal distribution, and implementing an algorithm to estimate its whole covariance matrix, which includes the expectation values for each parameter and all the correlations. This implementation is left to future developments.

The two-dimensional nature of stiffness must however be taken into account when seeking the values of the other parameters (such as G, B …). Once E and ν are identified with their uncertainties σE and σν, as indicated above, the functional dependences on E and ν, like Eqs. (5) and (6) can be exploited to derive the values of other parameters, like G, B, C11 and C44. However, care must be exercised in the evaluation of their uncertainties. The standard form of the error propagation formula holds for parameters which are uncorrelated. Making this assumption to estimate the uncertainties of, e.g., G and B, from σE and σν, would treat E and ν as uncorrelated, disregarding the information about their correlation which is intrinsically given by the shape of the confidence regions. This might lead to a substantial overestimation of the uncertainties of all the other parameters.

The most rigorous way to estimate the values and the uncertainties of all the parameters would require, as mentioned above, the estimation of the whole covariance matrix. A simpler procedure, which however fully takes into account the two-dimensional nature of stiffness, replicates the ‘circumscribed rectangle’ algorithm indicated above. Considering for instance G and B, it can be implemented in two different, but equivalent, ways. The first one substitutes the functional dependences of E and ν on G and B (Eqs. (5) and (6)) into the definition of the estimator S(E,ν) (Eq. (18)) obtaining its expression as S(B,G), and then repeats the whole procedure outlined above: computation of the normalized estimator at the nodes of a rectangular mesh in the (B,G) plane, determination of the confidence regions in this plane, identification of the parameters and of their uncertainties by the ‘circumscribed rectangle’ algorithm. In a completely equivalent way, once the confidence regions are determined in the (E,ν) plane, their frontiers can be mapped into the (B,G) plane by the functional dependences of G and B on E and ν, and the parameters and their uncertainties are then identified by the same algorithm. The two implementations only differ by the discretization errors, because a mesh which is rectangular and regular in the (E,ν) plane is neither rectangular nor regular in the (B,G) plane, and vice versa.

4.3. Anisotropic and non homogeneous cases

Until now, the data analysis procedure was discussed under the assumptions of isotropic and homogeneous medium, remembering that ‘homogeneous’ also includes layers in the ‘high thickness’ regime, in which the presence and the properties of the substrate are essentially irrelevant. The extension of the procedure beyond these assumptions is briefly indicated here.

In the case of an anisotropic, but still homogeneous, medium, the number of independent parameters needed to fully identify the tensor of the elastic constants is larger than two, meaning that stiffness is not a two-dimensional quantity, but a multi-dimensional one. Considering here for example the cubic symmetry, the independent parameters are three, which can be taken as C11, C12 and C44; for the hexagonal symmetry the number of independent parameters is five.

From the acoustic point of view, in a homogeneous anisotropic medium the acoustic waves can be not exactly longitudinal and transverse, and, more importantly, their velocity depends on the direction of the wavevector k, although not on its modulus. Therefore, any acoustic velocity, like e.g.vR, is a function of not only E and ν, but of more variables:vR=vR(C11,C12,C44,k/|k|). In the same way, each experimental value, like vR,i in Eq.(18), must be associated to the propagation direction (k/|k|)i along which it was measured. From the indentation point of view, anisotropy means that each measurement depends on the direction of the normal to the surface on which indentation is performed, and, if the indenter is not axisymmetric, it might also depend on the angular position of the indenter around its axis.

It is thus evident that the anisotropic case requires a data analysis which is substantially more complex, and, perhaps even more importantly, it requires measurements performed on surfaces of various orientations. This typically requires crystals with different cuts, which, in the case of films, might be not achievable. In the case of supported films, and remembering that we are now dealing with the ‘high thickness’ regime, it seems that the only type of anisotropy which can be characterized in detail is the transverse isotropy, meaning a film which is isotropic in its plane and has different properties in the perpendicular direction.

Remaining instead in the isotropic case, but relaxing the assumption of homogeneous medium, a couple of cases seem to be tractable in some detail. The first one is that of a homogeneous substrate with a homogeneous supported film in the ‘low thickness’ regime, in which the substrate becomes relevant. In this case ‘non homogeneous’ means a sharp discontinuity between two homogeneous media. The second case is that of a mild gradient of properties, distributed over a superficial layer.

In a homogeneous supported film, beside the Rayleigh wave of velocityvR, other surface waves can exist (Sezawa waves). For all of them, for given substrate properties (E,ν)subs and film thickness h, the velocity depends on the film properties and on the modulus k=|k| of the wavevector:vR=vR(E,ν,k|h,(E,ν)subs). By the same token, each measured value, like vR,i in Eq.(18), must be associated to the modulus ki at which it was measured, obtaining the dispersion relationsv(k). By AM the dispersion relation is obtained by measurements at different frequencies, while with BS it is obtained by measurements at different incidence angles.

From the indentation point of view, the critical point for a supported film in the ‘low thickness’ regime is the avoidance of the influence of the substrate. With BS the ‘low thickness’ regime occurs at thicknesses of the order of one micrometre or less, for which the avoidance of substrate influences is a crucial point, but it can be dealt with using sharp tips with accurate diamond area function calibrations (Barone et al. 2010) and shallow penetration depths. With AM, as discussed above, the ‘low thickness’ regime can occur at thicknesses of several micrometres, for which indentation can be easily performed avoiding the substrate influences. Therefore, for films of this thickness, the characterization by AM can be substantially improved by the combination with indentation.

The second case, that of a mild gradient of properties distributed over a finite depth, can be treated in an approximate way considering the average properties over the explored depth. In this case it is crucial that the indentation depth be comparable to the exploited acoustic wavelength(s) (Goodman & Derby, 2011). Although the two techniques do not average in the same way the properties at different depths, the above condition ensures that averaging is performed over the same depth interval. This procedure approximates the actual superficial gradient by a homogeneous equivalent film of average properties, over a perfectly homogeneous substrate of slightly different properties.


5. Experimental results

The combination of AM and indentation, after successful validation by a fused silica sample, was exploited to characterize TiN/NbN multilayers deposited by closed field magnetron sputtering (Bamber et al., 2001), assuming for the mass density of the multilayer the same value of TiN, 5210 kg/m3. Since the thicknesses of the single layers are much smaller than the acoustic wavelengths at which they are probed, the whole multilayer coating acoustically behaves as an equivalent homogeneous medium. In this case, the same authors state that the results were not conclusive, and they present a detailed analysis, focusing on the interpretation of the indentation results: they show that a main source of uncertainty lies in the analysis of indentation results from a multilayer.

Other two possible sources of uncertainty can be considered. Firstly, since the thickness of the whole multilayer is smaller than the acoustic wavelength, the layer is not in the ‘high thickness regime’, meaning that also for the acoustic behaviour the influence of the substrate is not negligible, and that the speed of the Rayleigh wave depends on frequency. The authors mention that by performing AM measurements at various frequencies the dispersion relation can be measured, but they do not give details on how they exploited this possibility. Secondly, by its same geometry a multilayer is not isotropic. It has in-plane isotropy, but different properties in the normal direction. This type of anisotropy corresponds to hexagonal symmetry (transverse isotropy), which requires five independent parameters to fully specify the elastic tensor. An isotropic model is the best approximation currently available with the experimental information, but it obviously introduces approximations.

A similar investigation (Goodman & Derby, 2011) was aimed at comparing the properties of the air side and the tin side of float glass, which is produced by floating molten glass on molten tin. The mass density values were taken from previous studies on the properties of soda-lime glass, and on the influence of the tin content. In this case, despite the need of considering a possible gradient of properties in the superficial layer, due to the occurrence of microdefects and to diffusion of tin, higher precision results were obtained, allowing to detect a measurable difference between the two sides. Due to the possible gradient of properties close to the surface, it was ensured that the depths probed by indentation and by AM be strictly comparable.

In a recent investigation (Garcia Ferré et al, 2013), the combination of indentation and BS was exploited to achieve a complete mechanical characterization of Al2O3 coatings deposited by pulsed laser deposition (PLD). Coatings were deposited on silicon wafer and on stainless steel substrates, ablating a 99.99% pure polycrystalline alumina target, and obtaining ultra-smooth samples (RMS <0,1 nm). SEM observations show a uniform, compact and fully dense microstructure; Raman spectroscopy and X-ray diffraction supply evidence of the amorphous structure of the coatings. Due to the absence of a crystalline structure and of a structure at mesoscopic level, like e.g. a columnar structure, the coating properties are assumed to be isotropic. The thickness of the measured coatings were between 4 and 8 μm. This means that by BS the layer is acoustically in the ‘high thickness’ regime, and that indentation can be performed with penetration depths, ranging from 300 to 500 nm, below one tenth of the layer thickness, thus avoiding the effect of the substrate.

Nanoindentation measurements were carried out by a nanoindenter (Micro-Materials, Ltd., Wrexham, UK) equipped with a diamond Berkovich tip. The reduced Young’s modulus and hardness were assessed from the indentation curves following the Oliver and Pharr approach discussed above. BS measurements were performed in the backscattering geometry, with incidence angles from 30° to 70°, using an Ar+ laser at the wavelength of 514,5 nm, and analyzing the scattered light by a tandem multipass Fabry-Perot interferometer of the Sandercock type.

Since the coatings are transparent, BS can measure bulk acoustic modes, but only the longitudinal one could be reliably measured. In the analysis of the corresponding data the value of the refractive index is needed: ellipsometric measurements supplied na = 1.647 ± 0.003. The analysis of acoustic data also requires the value of the mass density. For Al2O3 data are available data for α–alumina, i.e. sapphire: mass density ρs between 3970 and 3980 kg/m3, with index ns, at a wavelength of 514,5 nm, between 1,770 and 1,773 (Malitson, 1962; Defranzo & Pazol, 1993). The significant difference between na and ns, for the same nominal composition, indicates an appreciably different structure, such that the above value of ρs cannot be considered a sufficiently accurate estimate for the amorphous alumina: the coating mass density ρa was therefore derived from ρs, ns and na by the Clausius-Mossotti formula, also called Lorentz-Lorenz formula when written in terms of the refractive indexes (Beghi et al., 2011)


obtaining ρa = 3471 ± 17 kg/m3.

This combination of ellipsometry, BS and nanoindentation was called ‘EBN procedure’ (Garcia Ferré et al., 2013). Its validation was achieved by a bulk silica sample supplied by the producer of the indenter for calibration purposes, of nominal properties E = 72 GPa and ν = 0,17 (see Fig. 1). On this sample BS measurements could detect, beside the Rayleigh wave, both the bulk longitudinal and the bulk transversal waves. When both types of bulk waves are measurable in the same spectrum, the ratio of their velocities vl/vt can be directly obtained from the raw spectral data, irrespective of their calibration and of the values of the mass density and refractive index, and without being affected by possible imperfections of the scattering geometry. The ratio of velocities is thus measured with better precision than the single velocities.

For the silica sample the results from BS and from indentation are shown in Fig.1: the confidence bands resulting from the various measurements all intersect in a small region, showing a high degree of consistency among measurement techniques of completely different nature.

For the alumina coatings indentation measurements gave for the reduced modulus Er = 212.3 ± 5.4 GPa, and for hardness H = 10,32 ± 1,03 GPa. BS measurements could reliably measure the Rayleigh wave and the bulk longitudinal waves, of velocities vR = 4328 ± 42 m/s and vl = 8655 ± 36 m/s. The conversion from raw spectral data to vl is performed exploiting the measured value of refractive index na = 1.647; the conversion from velocities vl and vR to the elastic moduli is performed by the value of ρa = 3471 kg/m3 obtained, as indicated above, by the Lorentz-Lorenz formula. The procedures outlined above are performed in the (E,ν) plane, as well as in the (B,G) and the (C11,C44) planes, as shown in Fig. 2. The stiffness is thus characterized. The results obtained on the basis of the 95% confidence regions are summarized in Table 1. A detailed analysis of the uncertainties of statistical and of systematic nature was performed (Garcia Ferré et al., 2013). Basically, the errors in the measurement of the reduced modulus and the spectral frequencies are of statistical nature, and lead to a broadening of the confidence regions. An uncertainty in the value of mass density or of the refractive index acts instead like a calibration error, leading to a shift of the confidence regions, without appreciable broadening.

The results of Table 1 show that the stiffness of these PLD deposited alumina coating is significantly different from that of most ceramic coatings: namely, Young modulus is lower and Poisson’s ratio is higher, resulting in a stiffness which is close to that of steel. This finding could be rationalized in the light of TEM analysis, which shows that the coatings consist of a homogeneous dispersion of ultra-fine (2-5 nm) Al2O3 nanoparticles in an amorphous alumina matrix (Garcia Ferré et al., 2013). This metal-like stiffness is interesting, because at the substrate/coating interface a stress concentration occurs, due to the discontinuity of mechanical properties. The stress concentration can promote delamination, therefore a smaller stiffness discontinuity, reducing the stress concentration, can favour the durability of adhesion.

Figure 2.

is adapted from Garcia Ferré et al. 2013). Stiffness of the alumina coating, from results of indentation (reduced modulus Er, green lines) and BS (bulk longitudinal wave velocity vl, blue lines; Rayleigh wave velocity vR, red lines). For each quantity: central values and uncertainty bands. The same stiffness data are represented in the (E,ν) plane (Fig. 2a), the (B,G) plane (Fig. 2b) and the (C11, C44) plane (Fig. 2c). Thin black lines: the 68%, 90% and 99% confidence regions; thick black line: the 95% confidence region. In Fig. 2a the ‘circumscribed rectangle algorithm’ (see text) is also shown.

E (GPa) 195.3 ± 7.6
ν0.294 ± 0.020
G = C44 (GPa)75.5 ± 3.8
B (GPa)159.2 ± 11.8
C11 (GPa)259.8 ± 9.7

Table 1

adapted from Garcia Ferré et al. 2013). Stiffness parameters for an alumina coating, obtained by the combination of ellipsometry, BS and indentation (see Fig. 2).


6. Conclusions

Nanoindentation is among the most widespread techniques to perform the mechanical characterization of materials. It is applicable to bulk materials, and it is the standard characterization route for supported films whose thickness is not too low. Beside the value of hardness, instrumented indentation can supply a measurement of stiffness, by the careful analysis of a deformation process which involves both elastic and inelastic strains. However, stiffness is a two-dimensional quantity for isotropic materials, and has a higher number of dimensions in the anisotropic cases. Indentation supplies a one-dimensional estimate, and can achieve a full estimate only by some kind of assumption, typically about the value of Poisson’s ratio. The techniques which exploit propagating acoustic waves involve exclusively elastic strains: they therefore offer the most direct access to the elastic properties, and potentially their most accurate measurements. However, if only one acoustic mode can be detected, as it typically happens with AM and with BS of metallic samples, also acoustic techniques only supply a one-dimensional estimate. The combination of indentation and an acoustic method then allows to obtain the full two-dimensional measurement of stiffness, avoiding any assumption. With sufficiently transparent samples BS can detect further acoustic modes, supplying additional information and allowing to reduce the uncertainties. The consistency of results obtained by the two completely different methods could be checked by reference samples. The combined procedure, which is able to achieve the complete elastic characterization without needing assumptions, is particularly useful in the case of materials with novel properties, either for the novelty of the material or for the novelty of the production process. In these cases any assumption is not sufficiently grounded, and might lead to misleading results; the combined procedure then achieves a characterization of stiffness whose accuracy and reliability would not be reachable by each of the techniques alone. Furthermore, the presence of indentation simultaneously provides the measurement of hardness, intrinsically out of the reach of acoustic methods.


  1. 1. AuldB. A.1990Acoustic fields and Waves in SolidsRobert E. Krieger Publishing Company, Malabar, Florida
  2. 2. BamberM. J.CookeK. E.MannA. B.DerbyB.2001Accurate determination of Young’s modulus and Poisson’s ratio of thin films by a combination of acoustic microscopy and nanoindentationThin Solid Films398-399299305
  3. 3. BaroneA. C.SalernoM.PatraN.GastaldiD.BertarelliE.CarnelliD.VenaP.2010Calibration issues for nanoindentation experiments: direct atomic force microscopy measurements and indirect methods.Micros. Res. Tech. 739961004
  4. 4. BeghiM. G.BottaniC. E.PastorelliR.2001High accuracy measurement of elastic constants of thin films by surface Brillouin scatteringIn Mechanical properties of structural films, C. Muhlstein, C. & Brown, S.B. (Eds.), 109126ASTM STP 1413, American Society for Testing and Materials, Conshohoken, PA
  5. 5. BeghiM. G.EveryA. G.PrakapenkaV.ZininP. V.2012Measurement of the Elastic Properties of Solids by Brillouin Spectroscopy, in Ultrasonic and Electromagnetic NDE for Structure and Material Characterization: Engineering and Biomedical Applications, T. Kundu (Ed.), chapter 10, CRC Press, Boca Raton, FL.
  6. 6. BeghiM. G.Di Fonzo, F., Pietralunga, S., Ubaldi, C. & Bottani, C. E. (2011Precision and accuracy in film stiffness measurement by Brillouin spectroscopy,Review of Scientific Instruments, 102Paper 053107
  7. 7. BhushanB.LiX. D.2003Nanomechanical characterisation of solid surfaces and thin filmsInt. Mater. Rev., 48125164
  8. 8. CominsJ. D.EveryA. G.StoddartP. R.ZhangX.CrowhurstJ. C.HearneG. R.2000Surface Brillouin scattering of opaque solids and thin supported filmsUltrasonics, 38450458
  9. 9. DoernerM. F.NixW. D.1986A method of interpreting the data from depth-sensing indentation instruments,J. Mater. Res. 1n. 4, 601609
  10. 10. EveryA. G.2001The Elastic Properties of Solids: Static and Dynamic Principles, In: Handbook of Elastic Properties of Solids, Liquids, and Gases, M. Levy, H. Bass, R. Stern & V. Keppens (Eds.), Volume I: Dynamic Methods for Measuring the Elastic Properties of Solids, 336Academic Press, New York
  11. 11. FarnellG. W.AdlerE. L.1972Elastic wave propagation in thin layers.In: Physical Acoustics, W.P. Mason & R.N. Thurston (Eds.), 935127Academic, New York.
  12. 12. Fischer-crippsA. C.2011Nanoindentation, Springer (third edition), Chap. 8.
  13. 13. Garcìa Ferré, F., Bertarelli, E., Chiodoni, A., Carnelli, D., Gastaldi, D., Vena, P., Beghi, M.G. & Di Fonzo, F. (2013). The mechanical properties of a nanocrystalline Al2O3/a-Al2O3 composite coating measured by nanoindentation and Brillouin spectroscopy, Acta Materialia, 61 26622670 .
  14. 14. GoodmanO.DerbyB.2011The mechanical properties of float glass surfaces measured by nanoindentation and acoustic microscopyActa Materialia5917901799
  15. 15. HertzH.1881On the contact of elastic solids, J. Reine Angew. Math. 92156171translated and reprinted in English in Hertz’s Miscellaneous Papers, Macmillan & Co., London, 1896, Ch. 5.
  16. 16. KunduT.2012Mechanics of elastic waves and ultrasonics non-destructive evaluation, in Ultrasonic and Electromagnetic NDE for Structure and Material Characterization: Engineering and Biomedical Applications, T. Kundu (Ed.), chapter 1, CRC Press, Boca Raton, FL.
  17. 17. OliverW. C.PharrG. M.1992An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experimentsJ Mater Res, 715641583
  18. 18. SeberG. A. F.WildC. J.2003Nonlinear regressionWiley Interscience, Hoboken, NJ.
  19. 19. SneddonI. N.1965The relation between load and penetration in the axisymmetric Boussinesq problem for a punch of arbitrary profile.Int. J. Engng. Sci., 34757
  20. 20. ZhangX.StoddartP. R.CominsJ. D.EveryA. G.2001High-temperature elastic properties of a nickel-based superalloy studied by surface Brillouin scatteringJournal of Physics: Condensed. Matter, 1322812294
  21. 21. ZininP. V.2001Quantitative acoustic microscopy of solids, in Handbook of Elastic Properties of Solids, Liquids, and Gases, M. Levy, H. E. Bass and R. R. Stern & V. Keppens (Eds.), Volume I: Dynamic Methods for Measuring the Elastic Properties of Solids, 187226Academic Press, New York
  22. 22. ZininP. V.ArnoldW.WeiseW.BerezinaS.2012Theory and Applications of Scanning Acoustic Microscopy and Scanning Near-Field Acoustic Imaging, in Ultrasonic and Electromagnetic NDE for Structure and Material Characterization: Engineering and Biomedical Applications, T. Kundu (Ed.), chapter 11, CRC Press, Boca Raton, FL.

Written By

Fabio Di Fonzo, Francisco García Ferré, Dario Gastaldi, Pasquale Vena and Marco G. Beghi

Published: 28 August 2013