L Charleux
email: ludovic.charleux@univ-savoie.fr

V Keryvin

M Nivard

J.-P Guin

J.-C Sangleboeuf

Y Yokoyama

A method for measuring the contact area in instrumented indentation testing by tip scanning probe microscopy imaging

Keywords: Nanoindentation, Atomic force microscopy, Hardness, Elastic behavior, Finite element analysis

The determination of the contact area is a key step to derive mechanical properties such as hardness or an elastic modulus by instrumented indentation testing. Two families of procedures are dedicated to extracting this area: on the one hand, post mortem measurements that require residual imprint imaging, and on the other hand, direct methods that only rely on the load vs. the penetration depth curve. With the development of built-in scanning probe microscopy imaging capabilities such as atomic force microscopy and indentation tip scanning probe microscopy, last generation indentation devices have made systematic residual imprint imaging much faster and more reliable. In this paper, a new post mortem method is introduced and further compared to three existing classical direct methods by means of a numerical and experimental benchmark covering a large range of materials. It is shown that the new method systematically leads to lower error levels regardless of the type of material. Pros and cons of the new method vs. direct methods are 1 also discussed, demonstrating its efficiency in easily extracting mechanical properties with an enhanced confidence.

Introduction

Over the last two decades, Instrumented Indentation Technique (IIT) has become a widespread procedure that is used to probe mechanical properties for samples of nearly any size or nature. However, the intrinsic heterogeneity of the mechanical fields underneath the indenter prevents from establishing straightforward relationships between the measured load vs. displacement curve and any expected mechanical properties as it would be the case for a tensile test. Many models have been published in the literature in order to enable the measurement of properties such as an elastic modulus, hardness or various plastic properties. Despite their diversity, most of these models deeply rely on the accurate measurement of the projected contact area between the indenter and the sample's surface. The existing methods that are dedicated to estimating the true contact area can be classified into two subcategories: the direct methods which rely on the sole load vs. displacement curve [START_REF] Doerner | A method for interpreting the data from depth-sensing indentation instruments[END_REF][START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF][START_REF] Loubet | Vickers indentation curves of magnesium oxide (MgO)[END_REF] and the post mortem methods that use additional data extracted from the residual imprint left on the sample's surface. For example, Vickers, Brinell and Knoop hardness scales rely on post mortem measurements of the geometric size of the residual imprint. However, in the case of Vickers hardness, the contact area is only estimated through the diagonals of the imprint, the possible effect of piling-up or sinking-in is then neglected.

Other post mortem methods use indent cross sections to estimate the projected contact area [START_REF] Beegan | The nanoindentation behaviour of hard and soft films on silicon substrates[END_REF]5]. In the 1990s, the development of nanoindentation led to a growing interest in direct methods because they do not require time consuming post mortem measurement of micrometer or even nanometer scale imprints, typically using Atomic Force Microscopy (AFM) or Scanning Electron Microscopy (SEM). Uncertainty level on direct measurements remains high, mainly because of the difficulty to predict the occurrence of piling-up and sinking-in. Oliver and Pharr have eventually considered this issue as one of the "holy grails" in IIT [START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF]. Recent development in Scanning Probe Microscopy (SPM) using the Indentation Tip (ITSPM) brought new interest in post mortem measurements. Indeed, ITSPM allows systematic imprint imaging without manipulating the sample or facing repositioning issues to find back the imprint to be imaged. Yet ITSPM imaging technique suffers from drawbacks when compared to AFM: it is slower, it uses a blunter tip associated with a much wider pyramidal geometry and a higher force applied to the surface while scanning. While the later may damage delicate material surfaces, the formers will introduce artifacts. Nonetheless, these artifacts will not affect the present method. In addition, ITSPM only allows for contact mode imaging, non contact or intermittent contact modes are not possible.

As a consequence, only the techniques based on altitude images can be used with ITSPM and there is a need for new methods as very recently reviewed by Marteau et al. [6]. This article introduces a new post mortem procedure that relies only on the altitude image and that is therefore valid for most types of SPM images, including ITSPM. In this paper, a benchmark based on both numerical indentation tests as well as experimental indentation tests on properly chosen materials to span all possible behaviors is first introduced.

Then, the existing direct methods are reviewed and a complete description of the proposed method is given. These methods are then confronted using the above mentioned benchmark and the results are finally discussed.

Numerical and experimental benchmark

A typical instrumented indentation test features a loading step where the load P is increased up to a maximum value P max , then held constant in order to detect creep and finally decreased during the unloading step until contact is lost between the indenter and the sample. A residual imprint is left on the initially flat surface of the sample. During the test, the load P as well as the penetration of the indenter into the surface of the sample h is continuously recorded and can be plotted as shown in Figure 1. For most materials, the unloading step can be cycled with only minor hysteresis, it is then assumed that only elastic strains develop in the sample. As a consequence, the initial slope S of the unloading step is called the elastic contact stiffness. Useful data can potentially be extracted from both the load vs. displacement curve and the residual imprint. The contact area A c is defined as the projection of the contact zone between the indenter and the sample at maximum load on the plane of the initially flat surface of the sample. The displacement of the indenter h is controlled and the force P is recorded.

Numerical approach

The dimensions of the mesh are chosen to minimize the effect of the farfield boundary conditions. The typical ratio of the maximum contact radius and the sample size is about 2 × 10 3 . The problem is solved using the commercial software ABAQUS (version 6.11, 3ds.com). The numerical model is compared to the elastic solution from [START_REF] Love | Boussinesq's problem for a rigid cone[END_REF] (see [9, 10]) using a blunt conical indenter (Ψ = 89.5˚) to respect the purely axial contact pressure hypothesis used in the elastic solution. The relative error is computed from the load vs. penetration curve and is below 0.1%. Pre-processing, post-processing and data storage tasks are performed using a dedicated framework based on the open source programming language Python 2.7 [START_REF] Van Rossum | Interactively testing remote servers using the Python programming language[END_REF]12,13] and the database engine SQLite 3.7 [14]. The indented material is assumed to be isotropic, linearly elastic. The Poisson's ratio ν has a fixed value of 0.3 and the Young's modulus is referred to as E. The contact between the indenter and the sample's surface is taken as frictionless. Two sets of constitutive equations (CE1 and CE2) are investigated in order to cover a very wide range of contact geometries and materials: CE1 This first constitutive equation used in this benchmark is commonly used in industrial studies and in research papers on metallic alloys [START_REF] Bucaille | Determination of plastic properties of metals by instrumented indentation using different sharp indenters[END_REF][START_REF] Dao | Computational modeling of the forward and reverse problems in instrumented sharp indentation[END_REF][START_REF] Ma | New relationship between Young's modulus and nonideally sharp indentation parameters[END_REF]18]. It uses J 2 -type associated plasticity and an isotropic Hollomon power law strain hardening driven by the tensile behavior (stress σ T , strain ε T ) given by Eq. 1:

σ T =    Eε T , for: σ T ≤ σ Y T σ Y T Eε T σ Y T n , for: σ T > σ Y T (1)
Plastic parameters are the tensile yield stress σ Y T and the strain hardening exponent n.

CE2

The second constitutive equation is the Drucker-Prager law [19] which was originally dedicated to soil mechanics but was also found to be relevant on Bulk Metallic Glasses (BMGs) [20,[START_REF] Patnaik | Spherical indentation response of metallic glasses[END_REF]22,23] and some polymers [24]. The yield surface is given by Eq. 2 where q is the von Mises equivalent stress in tension and p the hydrostatic pressure.

Perfect plasticity is used in conjunction with an associated plastic flow.

The plastic behavior is controlled by the compressive yield stress σ Y C and the friction angle β that tunes the pressure sensitivity.

q -p tan β -1 -1/3 tan β σ Y C = 0 (2)
Dimensional analysis [25,26,[START_REF] Cheng | Scaling approach to conical indentation in elastic-plastic solids with work hardening[END_REF] is used to determine the influence of elastic and plastic parameters on the contact area A c :

A c =    h 2 max Π C E1 (ν, σ Y T /E, n) h 2 max Π C E2 (ν, σ Y C /E, β) (3)
In this equation, h max is the maximum value of penetration of the indenter into the sample's surface. In both cases, the dimensionless functions show that, since the Poisson's ratio has a fixed value (ν = 0. 1 and2. After each simulation, a load vs. displacement into surface curve and an altitude SPM like image using the Gwyddion (http://gwyddion.net/) GSF format are extracted. The use of such a procedure allows one to consider both numerical and experimental tests in the benchmark and to derive mechanical properties in the same way.

Since the simulations are two-dimensional axisymmetric, the contact area is computed as A c = πr 2 c where r c stands for the contact radius of the contact zone (see Fig. 2).

Experimental testing

The tested materials (see Table 3) are chosen in order to cover a very wide range of contact geometries, from sinking-in (FQ), intermediate behavior (WG), and to piling-up (BMG). Glasses are chosen over metallic alloys because they exhibit negligible creep for temperatures well below the glass transition temperature, no visible size effect and are very homogeneous and isotropic in the test conditions. The FQ and WG samples are tested as received (please note the WG sample was test on it's "air" side) whereas the BMG sample is polished. Nanoindentation testing is performed using a commercial Hysitron TI950 triboindenter. During each test, the load is increased up to P max = 10 mN with a constant loading rate d P/d t = 5 × 10 -5 N/s. The load is then held for 10 s and relieved with a constant unloading rate d P/d t = -1 × 10 -4 N/s. Four tests are performed on each sample and each residual imprint is scanned with the built-in ITSPM device with an applied normal force of 2 µN as summarized in Figure 3. Tests are load controlled and the maximum load is set to P max = 10 mN. The true contact area A c is not known as in the case of the numerical simulations. It is then estimated through the Sneddon's Eq. 4 [START_REF] Sneddon | The relation between load and penetration in the axisymmetric Boussinesq problem for a punch of arbitrary profile[END_REF][START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF] and is called A c,SN . Young's moduli E, and the Poisson's ratios ν of each sample are known prior to indentation testing from the literature or from ultrasonic echography measurements (cf.

Table 3). Recalling that S is the initial unloading contact stiffness (cf. Figure 1), we have:

A c,SN = π 4 S 2 β 2 E 2 eq
where:

   E eq = E 1 -ν 2 β = 1.05 (4)

Methodology review

Direct methods

Direct methods rely on the sole load vs. penetration curve (P, h) to determine the contact height h c using equations given in Table 4. Let us recall that h c < h in the case of sinking-in (as seen in Fig. 2) and h c > h for piling-up.

Three direct methods are investigated in this paper : DN The Doerner and Nix method [START_REF] Doerner | A method for interpreting the data from depth-sensing indentation instruments[END_REF] was one of the first to be published (along with similar work done by Bulyshev et al. [28]) and it provided the basic relationships later improved by the two other methods.

OP The Oliver and Pharr method [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation[END_REF][START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF] is an all-purpose method that is widely used in the literature, commercial software and standards. The main drawback of this method is that it cannot take piling-up into account.

LO The Loubet method [START_REF] Loubet | Vickers indentation curves of magnesium oxide (MgO)[END_REF] is an alternative to the OP method, especially for materials exhibiting piling-up.

Regardless of the chosen method, the value of h c is used to compute the value of the contact area A c thanks to the Indenter Area Function A c (h c ) (IAF). The IAF depends on the theoretical shape of the indenter as well and on its actual defects which are measured during a calibration procedure. Different tip calibration methods are used in the literature:

• Measurement either of the indenter geometry or the imprint geometry made on soft materials for multiple loads using AFM or other microscopy techniques [START_REF] Pethica | Hardness measurement at penetration depths as small as 20 nm[END_REF].

• The IAF introduced by Oliver and Pharr [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation[END_REF][START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF] requires a calibration procedure on a reference material using only the (P, h) curve:

A c,OP (h c ) = C 0 h 2 c + C 1 h c + C 2 h 1/2 c + C 3 h 1/4 c + . . . + C 8 h 1/128 c (5)
Where the (C i ) 0≤i≤8 factors are fitting coefficients obtained from a calibration procedure on fused quartz. For a given indenter, the value of the C i coefficients depend on the penetration depth range used for the calibration procedure. In the case of a perfect modified Berkovich tip, C 0 = 24.5 and C i>0 = 0.

• The method introduced by Loubet (see [START_REF] Loubet | Vickers indentation curves of magnesium oxide (MgO)[END_REF]) :

A c,LO (h) = k h c + ∆h 2 (6)
It is assumed that the only origin of the defects is tip blunting. Then, k comes from the indenter's theoretical shape (k = 24.5 here) and ∆h is the offset caused by the tip defect and is calibrated using a linear fit made on the upper portion of the ( P, h) curve. This procedure can be performed on any material exhibiting neither significant creep nor size effect, typically fused quartz. This method is intrinsically very efficient when the penetration is high compared to ∆h.

In the experimental benchmark, all tests are performed at h max ≥ 250 nm using a diamond modified Berkovich tip that exhibits a truncated length ∆h = 17.8 ± 1.74 nm 1 . Theses values were calibrated on the FQ sample. As a consequence, the IAF introduced by Loubet is used on every direct method.

By contrast, numerical simulations use a perfect tip so that the IAF is A c (h) = 24.5h 2 .

Proposed Method (PM)

SPM imaging grants access to a mapping of the altitude of the residual imprint. It is assumed that the surface of the sample is initially plane and remains unaffected far from the residual imprint. This plane is extracted from the raw image using a disk shaped mask centered on the imprint and a scan by scan linear fit on the remaining zone. It is considered as the reference surface and is subtracted to the raw image to remove the tilt of the initial surface. Under maximum load, the contact contour can exhibit either sinking-in and piling-up. In the first case, its altitude is decreased vis-à-vis the reference surface. This behavior is typical of high yield strength materials such as fused quartz. On the opposite, piling-up occurs when the the increased and is typically triggered by unconfined plastic flow around the indenter as usually observed on low strain hardening metallic alloys. When the indenter is not axisymmetric (as it is the case for pyramids), both sinking-in and piling-up may occur simultaneously. For example, pyramidal indenters can produce piling-up on their faces and sinking on their edges (or no piling-up to the least). During the unloading step, the whole contact contour is pushed upward with only minor radial displacement. A residual piling-up may form even if sinking-in initially occurred under maximum load. We now focus on a half cross-section of the imprint starting at the bottom of the imprint and formulate two assumptions:

1. The highest point of any half cross-section is the summit of the residual piling-up.

2. The summit of the residual piling-up indicates the position of the contact contour.

As a consequence, the highest point of the cross-section gives the radial position of the contact zone boundary. However, from an experimental point of view, the roughness of the sample's surface may make the true position of the residual piling-up's summit unclear. This issue is particularly true for materials exhibiting high levels of sinking-in such as fused quartz but also for most materials along the edges of pyramidal indenters. In order to limit the effect of surface roughness on the radial localization of the contact zone boundary, each profile is slightly rotated by a small angle value α along an axis per-pendicular to the cross-section plane and running through the bottom of the imprint (i. e. the lowest point of the profile). The whole contact contour is then determined by repeating the process in all directions (θ = 0˚to 360˚)

then the contact area A c,P M is calculated. A graphical representation of key steps of the method is made in Figure 4.

[PLEASE INSERT FIGURE 4 HERE]

The optimum tilt value of the tilt angle α is chosen in order to deal with three potential artifacts detailed below and emphasized using three experimental residual imprints in Fig. 5: 1. Materials exhibiting low or no residual piling-up such as FQ cannot be treated with the proposed method when α = 0˚because the positions the highest points of each the extracted half cross sections are driven by the roughness. In this case, positive values of the tilt angle (α ≥ 2.5˚) fully address the issue.

2. The surface roughness of the sample affects the accuracy of the method, especially on materials that show no residual piling-up because there is a competition between the is a competition between piling-up and roughness in the highest point identification (see FQ and BMG on Fig.

5.

Again, positive values of the tilt angle (α ≥ 2.5˚) solve the problem by lowering the roughness peaks located outside the imprint as visible.

3. A particular artifact is visible in the directions of the indenter's edges.

Indeed, the residual piling-up is always very low in these directions because the displacement field tends to push the material towards the faces of the indenter. The attack angle between the edges and the sample's surface is also very low 2 As a consequence, a high value of the tilt angle such as α ≥ 5˚creates artifacts in this zone as clearly visible on FQ and WG.

Given the three above points, the optimum value of the tilt angle is found to be α = 2.5˚. The three values of the tilt angle have also been tested in on the numerical simulations. Even if the overall effect of the tilt angle is lower in this case because the simulations are two dimensional and also because the numerical samples have no roughness, only the first issue is observed and the optimum value is the same the one found experimentally.

Results and discussion

FEM benchmark results

The four methods are confronted on the numerical benchmark and their ability to accurately compute the contact area A c is challenged. The results focus on the relative error e between A c and the contact area predicted by each method. Please note that a 10% relative error on the contact area means roughly a ≈ 10% relative error on hardness but only ≈ 5% on the elastic modulus (see Eq. 4). The results are plotted in Figures 6 and7 for constitutive equations CE1 and CE2 respectively and a summary of the key statistics is given in Table 5. For Fig. 6 and 7, (a) and (b) represent the relative error and its absolute value respectively. The later was chose to emphasize the magnitude of the error while (c) indicates whether piling-up or sinking-in occurs. It is also chosen to measure the success rate of the methods through their ability to match A c within a ±10% error. This value was chosen since even if challenging, it is still realistic from an experimental point of view.

These data can be discussed individually for each method:

• The DN method systematically tends to underestimate A c regardless of the type of constitutive equation and of the occurrence of piling-up or sinking-in. The magnitude of the relative error is the highest among the four tested methods. This lack of accuracy can be put into perspective by recalling that the DN method states that the contact between the indenter and the sample behaves as if the indenter was a flat punch during the first stages of the unloading process. This approach was later proved to be too restrictive by Oliver and Pharr [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation[END_REF] who improved it by taking into account the actual shape of the indenter through the ε coefficient. As ε < 1 in the case of the modified Berkovich tip, the value of the contact depth h c is systematically increased (see Table 4).

• As stated above, the OP method drastically improves the overall performances of the DN method. However, its error level depends strongly on the type of contact behavior (i. e. piling-up or sinking-in) and the mechanical properties of the tested material. Typically, it performs well for the CE1 law when the strain hardening exponent n verifies n > 0.2. It also performs well (relative errors below 10%) on materials exhibiting very high yield strains (higher than 4%) in the case of CE2. The main drawback of the method is its intrinsic inability to cope with piling-up since h c /h can never be higher than 1. This is particularly visible for low values of the strain hardening exponent (CE1: n ≤ 0.1) and with CE2 when the compressive yield strain σ C /E is lower than 3%. The OP method has a low success rate (see Table 5) but this tendency has to be mitigated by the fact that it is very efficient for a large number of metallic alloys, which can be described by CE1 and exhibit moderate values of hardening exponents.

• The LO method allows h c /h > 1 values and is then recommended for piling-up materials; it is overall very efficient with CE2 type materials.

The drawback is that it tends to overestimate the contact area when sinking-in occurs; this is particularly true in the case of CE1 with moderate to high hardening exponents (n ≥ 0.1). These observations are in agreement with the results of Cheng and Cheng [START_REF] Cheng | Effects of 'sinking in' and 'piling up' on estimating the contact area under load in indentation[END_REF] regarding the influence of piling-up and sinking-in on the direct estimation of the contact depth. The success rate of this method is the highest among the direct methods and it is clearly the best available direct method for CE2 type materials and for low hardening CE1 type materials.

• The proposed method exhibits a 100% success rate (with the ±10 % relative error target) and an average absolute relative error of 2.5%.

The error level remains stable regardless of both the type of constitutive equation and its parameters. This result highlights the fact that when experimentally possible, the use of such a post mortem method will improve drastically the overall error level of the contact area measurement.

Experimental benchmark results

The results of the experimental benchmark are represented in Fig. 8.

The tendencies observed in the numerical benchmark are confirmed. The DN method systematically underestimates the contact area. The OP and LO methods perform well only on a given spectrum of contact behaviors: both methods give accurate results on the FQ sample; this is consistent with the fact that both of them were optimized to use this material as a reference.

While the LO method also exhibits a low error level on the WG sample, the OP method leads to an unexpected high error level. It is supposed that even if the WG sample has a very high yield strain, it has no strain hardening mechanism and it is then out of the scope of the OP method. The BMG which exhibits a large residual piling-up obviously leads the OP method to underestimate drastically the contact area. The LO method performs better although it also underestimates the contact area. This later method systematically exhibits relative errors of ±10% while the method proposed in this paper is even more reliable with errors lower than 5%. We observe that the direct methods overall performances are better than in the case of the numerical benchmark.

The contact friction, which is neglected in the numerical benchmark, may improve the accuracy of the direct methods without affecting the proposed method.

Pros and cons

Both benchmarks highlight the precision gap between the new method and the existing direct methods. However, the proposed method differs by nature from the three direct methods it is compared to. This section emphasizes the pros and cons of this method:

Disadvantages:

• The proposed method relies on SPM imaging of the residual imprint while direct methods do not. However, indentation devices tend to be equipped with ITSPM capability that can be used automatically in conjunction with the indentation testing itself with only a small increase in test duration.

• The method relies on the assumption that the imprint is unchanged between the end of the test and the imaging procedure. This means that the proposed method should not be applied to materials exhibiting time-dependent mechanical behavior such as polymers (see [START_REF] Chatel | Original in situ observations of creep during indentation and recovery of the residual imprint on amorphous polymer[END_REF]) and even pure metals. In order to demonstrate this limitation, nanoindentation tests have been run an electropolished pure single crystal aluminum sample (Al) and a mechanically polished pure tungsten sample (W). Those two materials were chosen because of their low elastic anisotropy [START_REF] Hirth | Theory of dislocations[END_REF]. A loading function similar to the one used in [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation[END_REF] was used to limit the effect of creep on the unloading stiffness S. The method appears to work flawlessly on both samples as none of the three artifacts mentioned in part 3.2 are observed. However, the value of the contact area is systematically overestimated by a factor ≈ 1.6 and ≈ 2.7 respectively on the Al and W samples. The evolution of the residual imprint after the indentation test is made possible by creep caused by residual stress.

Advantages:

• The proposed method can be run automatically, it requires no adjustable parameters and is user independent.

• Sample holder and machine stiffness affect the measurement of the penetration into surface h as well as the measured contact stiffness S and, as a consequence, they also affect all direct methods. The value of the machine stiffness can be measured once and for all while the sample holder's stiffness may change between two samples and requires systematic calibration. This concern is particularly true in the case of small samples such as fibers as well as very hard materials (such as carbides). The contact area measurement provided by the proposed method does not rely on h and is then insensitive to the effect of those spurious stiffness issues.

Yet, let us note that while the value of the contact area A c,P M is unaffected by stiffness issues, the value of the contact stiffness S is of course affected. As a consequence, the value of the hardness probed with the proposed method is free from any stiffness concern (as H = P max /A c ) while the value of the reduced modulus E * still requires stiffness calibration (as E * ∝ S/ A c ).

• The method does not require any tip calibration procedure and is compatible with all tip shapes.

• The method is unaffected by erroneous surface detection also because it does not rely on h.

Conclusion

We have proposed a new procedure to estimate the indentation contact area based on the residual imprint observation using altitude images produced by SPM. This area is the key component of instrumented indentation testing for extracting mechanical properties such as hardness or elastic modulus. For the estimation of this contact area, the method has been confronted with three widely used direct methods. We have showed, by means of an experimental and numerical benchmark covering a large range of contact geometries and materials, that the proposed method is far more accurate than its direct counterparts regardless of the type of material. We have also discussed the fact that such post mortem procedures are indeed more time consuming than direct methods; yet they are the future alternative to direct methods with the development of indentation tip scanning probe microscopy techniques. We have also emphasized the fact that this new method has numerous advantages: it can be automated, it is user independent, it is unaffected by stiffness issues and does not require any indenter calibration.

Parameter

Load on sample, P [N]

Loading data Unloading data Unloading power law fit:

P = P max (h -h f )/(h max -h f ) m
Unloading tangent: The experimental curve generally also includes an holding step which is not represented in this case. The contact stiffness S is the unloading step's initial slope. However, the direct determination of S via the upper part of the step is unreliable as it uses only a small part of the curve. For increased accuracy, the whole step is systematically fitted by a power law function which is used to compute back the contact stiffness S as initially recommended in [START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF]. Von Mises eq. stress to tensile yield stress ratio, q/σ Y T [%]

P = P max + S h -h max
Figure 2: Representation of the 2D axisymmetric FEM model including a rigid conical indenter and a deformable elastic-plastic sample. In this particular case, the sample's material is based on CE1 (σ Y T /E = 0.07, n = 0.3). The model is represented at maximum penetration into surface h ma x . (a) The deformed mesh is plotted showing the finely meshed zone containing the contact zone and most of the plastic zone. This zone is initially filled with square shaped elements. The mesh is gradually coarsened away from the contact zone. The whole mesh is not represented since its total size is about 10 3 times the size of the finely mesh zone. (b) The gradient represents the ratio of the von Mises equivalent stress in tension q to the tensile yield stress σ Y T with the corresponding scale given by the bottom color bar. (Right) ITSPM of one of the four tests is plotted for each material. Altitudes below nm are masked in order to emphasize the shape and size of the residual piling-up. There is no visible residual piling-up on FQ, a low altitude circular one on WG and a higher altitude one with summits on the faces on BMG. These observations justify the choice of the three selected materials as they cover all possible contact geometries. Figure 5: Contact contours produced by the proposed method on the three residual imprints measured on each samples composing the benchmark. The imprints were produced using the experimental protocol described in the section 2.2 and the Fig. 3. Three values of the tilt angle α are investigated on each imprint: no tilt (α = 0˚), the proposed value (α = 2.5˚) and last higher one (α = 5˚).

Finite

  element modeling (FEM) simulations are performed using a twodimensional axisymmetrical model represented in Figure 2. The sample is meshed with 3316 four-noded quadrilateral elements. The indenter is considered as a rigid cone exhibiting an half-angle Ψ = 70.29˚to match the theoretical area function of the Vickers and modified Berkovich indenters [7].

  3), only the yield strains (ε Y T = σ Y T /E, in the case of CE1 and ε Y C = σ Y C /E, in the case of CE2) and the dimensionless plastic parameters (n in the case of CE1 and β in the case of CE2) have an influence on the contact area A c . As a consequence, the value of the Young's modulus E has a fixed arbitrary value E = 1 Pa and only the values of the yield stresses σ Y T and σ Y C , the hardening exponent n and the friction angle β are modified. The simulated range of these parameters are given in Tables

Figure 1 :

 1 Figure1: Typical sharp indentation load on sample vs. displacement into surface curve. The test is split into a loading step and an unloading step. The experimental curve generally also includes an holding step which is not represented in this case. The contact stiffness S is the unloading step's initial slope. However, the direct determination of S via the upper part of the step is unreliable as it uses only a small part of the curve. For increased accuracy, the whole step is systematically fitted by a power law function which is used to compute back the contact stiffness S as initially recommended in[START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF].

Figure 3 :

 3 Figure3: Representation of the experimental tests carried out on the three selected materials FQ, WG and BMG (Four tests per sample). (Left) The load vs. displacement into surface (P, h) curves are plotted. The three loading steps are rather similar while the unloading steps are very different. The FQ exhibits very high elastic recovery while the BMG has a low one. (Right) ITSPM of one of the four tests is plotted for each material. Altitudes below nm are masked in order to emphasize the shape and size of the residual piling-up. There is no visible residual piling-up on FQ, a low altitude circular one on WG and a higher altitude one with summits on the faces on BMG. These observations justify the choice of the three selected materials as they cover all possible contact geometries.

  , θ i ) z(r, θ i )r tan α z(x, y) [nm]

Figure 4 :

 4 Figure 4: Key steps of the proposed method exemplified on a Berkovich indent performed on window glass with P ma x = 10 mN . (a) Raw ITSPM image on which the tilt has to be corrected using a linear least squares fit on each individual scan with a circular mask around the indented region. (b) An half cross section at a given angle θ i is extracted. (c) The half cross section z(r, θ i ) is rotated by a tilt angle α = 2.5˚in order to isolate more efficiently the edge of the contact zone. (d) The process is repeated for all required values of θ i and the whole contact zone is extracted (in black).

β

  = 30 o A = A c,DN A = A c,OP A = A c,LO A = A c,P M Compressive yield strain, σ Y C /E [%]

Figure 7 :

 7 Figure 7: Benchmark results of 3 direct methods and of the Proposed Method (PM) in the case of CE2. See Figure 6 for the complete description of the structure of the figure.

Table 1 :

 1 Simulated range of the dimensionless ratios for the numerical simulations using the constitutive equation CE1. The "Number" column stands for the number of values chosen as simulation inputs in the given range.

		Description	Range	Number
	σ Y T /E n	Tensile yield strain Strain hardening exponent	0.001 → 0.01 0.0 → 0.3	10 4
	Parameter	Description	Range	Number
	σ Y C /E β	9 0.0˚→ 30˚4 Compressive yield strain 0.01 → 0.05 Friction angle

Table 2 :

 2 Simulated range of the dimensionless ratios for using the constitutive equation CE2.The "Number" column stands for the number of values chosen as simulation inputs in the given range.

	Label E (GPa)	ν	Provider/Description	Ref.
	FQ	73	0.16	Hysitron Inc.
	WG	71.5	0.23	Planilux, Saint Gobain Inc.	[34]
	BMG	89.3	0.363 Zr 55 Cu 35 Al 10 , 385˚C annealed [35]

Table 3 :

 3 Description of the three tested samples: a fused quartz sample (FQ), a window glass sample (WG) and a zirconium based BMG. Materials properties used in this article are detailed: Young's modulus E, Poisson's ratio ν. Reference to the literature are provided when available. In the case of fused quartz, the properties are given by the provider and are in agreement with the literature.

	Method	Contact depth	Coeff.	Ref.
	DN OP LO	h c,DN = h max -P max /S h c,OP = h max -εP max /S h c,LO = α h max -P max /S	ø ε = 0.72 [2] [1] α = 1.2 [3]

Table 4 :

 4 Equations used to compute the indentation contact depth h c for the three direct methods used in this article.

	Method min(e) [%] max(e) [%] |e| [%] |e| ≤ 10% [%]
	DN	-46.1	-10.4	26.3	0.0
	OP	-39.4	-4.9	16.8	30.2
	LO	-22.3	29.1	10.5	60.5
	PM	-5.7	8.9	2.5	100.0

Table 5 :

 5 Numerical benchmark statistics of the three direct methods and the Proposed Method (PM) on both CE1 and CE2. Please note that 40 simulations were run for CE1 and 36 for CE2. As a consequence, the weight of CE1 is slightly higher than the weight of CE2 in the statistics. Mathematical notations are: e = (A -A c )/A c is the relative error on the contact area A computed by each method relatively to the the true contact area A c , |e| is its absolute value and |e| the arithmetic mean value of its absolute value. The last row displays the success rate of each method which is the proportion of the simulations on which the relative error is in the ±10% range.

	0

This value is the average value of ∆h on all the twelve tests performed on the three samples used in the experimental benchmark and the error is the represented ± one standard deviation.

The angle is 13˚in the case a of modified Berkovich tip.

Acknowledgements

The authors would like to thank the Brittany Region for its support through the CPER PRIN2TAN project and the European University of Brittany (UEB) through the BRESMAT RTR project. Vincent Keryvin also acknowledges the support of the UEB (EPT COMPDYNVER).

 1). On each simulation, the true projected contact area (A c ) computed by FEM, the contact areas estimated from the three direct methods (A c,DN , A c,OP and A c,LO ) and the contact area from the proposed method (A c,P M ) are calculated. (a) The relative error e between the true projected contact area and its four estimations is plotted. Each plot represents a different value of the hardening exponent n. (b) The absolute value of the the relative error |e| is represented in order to emphasize the accuracy of each method (see Table 5). (c) The contact depth h c stands for the axial distance between the edge of the contact zone and the the summit of the indenter. The relative difference between the contact depth and the penetration h is plotted to indicate where the piling-up occurs ((h ch)/h > 0) and when sinking-in occurs ((h ch)/h < 0). This subplot helps in understanding the relationships between the occurrence of piling-up and the accuracy of a given method.