Phase-field modeling of multivariant martensitic microstructures and size effects in nano-indentation

A finite-strain phase-field model is developed for the analysis of multivariant martensitic transformation during nano-indentation. Variational formulation of the complete evolution problem is developed within the incremental energy minimization framework. Computer implementation is performed based on the finite-element method which allows a natural treatment of the finite-strain formulation and of the contact interactions. A detailed computational study of nano-indentation reveals several interesting effects including the pop-in effect associated with nucleation of martensite and the energy-lowering breakdown of the symmetry of microstructure. The effect of the indenter radius is also examined revealing significant size effects governed by the interfacial


Introduction
Instrumented micro/nano-indentation is a powerful and highly popular experimental technique for characterization of material behavior at small scales at which other techniques are not applicable or are more difficult (Oliver and Pharr, 2004;Schuh, 2006;Fischer-Cripps, 2011). Instrumented indentation has been applied to virtually any material system. This includes shape-memory alloys (SMAs) which are the subject of the present work (e.g. Gall et al., 2001;Zhang and Komvopoulos, 2006;Frick et al., 2006;Muir Wood and Clyne, 2006;Crone et al., 2007;Amini et al., 2011). The interest in the SMAs is due to their spectacular features, notably pseudoelasticity and shapememory effect, that result from the crystallographically reversible martensitic phase transformation (Bhattacharya, 2003).
In the pseudoelastic regime, which is of great importance due to numerous engineering and biomedical applications, the inelastic deformation resulting from the stress-induced martensitic transformation vanishes (completely or nearly completely) upon unloading as a result of the reverse transformation. This concerns also the deformation during indentation. Accordingly, unlike in other materials, in a pseudoelastic SMA, the indentation does not leave any residual imprint. Thus the measured load-indentation depth curve is the only response that is available for characterization of the material. Note that the topography of the residual imprint is an important material characteristic that can be used, for instance, to support the indentation-based identification of mechanical properties, such as the hardening curve in plasticity (e.g. Bolzon et al., 2004;Kucharski and Mróz, 2007) and crystal plasticity (e.g. Petryk et al., 2017). Note also that the martensitic microstructure developed in a pseudoelastic SMA during indentation disappears upon the reverse transformation during unloading, see e.g. Pfetzing-Micklich et al. (2013), and is not available for experimental examination. Modeling is thus apparently the only way to examine the microstructure, clearly with all the related limitations.
Computational modeling and simulation of indentation in SMAs have been the subject of a number of studies reported in the literature. Macroscopic models of polycrystalline SMAs are relevant at higher scales at which the indenter radius (or the deformed volume in sharp indentation) is large compared to the grain size. The respective simulations can be found, e.g., in Muir Wood and Clyne (2006); Yan et al. (2007); Zhang et al. (2007). Simulations of indentation of SMAs at lower scales are much more scarce. Micromechanical models, such as the crystal-plasticity-like models, are applicable to single crystals or to individual grains in polycrystalline aggregates, see Dhala et al. (2019). Models of this class assume that individual phases or variants occupy the same (macroscopic) volume and are represented by the respective volume fractions. They are thus relevant for intermediate length scales, still considered macroscopic from the point of view of spatial arrangement of individual variants and phase boundaries. On the other hand, atomistic simulations employing molecular dynamics (MD) can be applied to simulate nano-indentation taking into account the discrete atomic structure, however, such simulations are limited to very small spatial and temporal scales (Pfetzing-Micklich https://doi.org/10.1016/j.mechmat.2019.103267 Received 2 July 2019; Accepted 27 November 2019 computational study is then carried out aimed at the analysis of the indentation-induced martensitic microstructures. The study reveals several interesting effects including the pop-in effect associated with nucleation of martensite and formation of energy-lowering non-symmetric microstructures in otherwise symmetric problems. Finally, the effect of the indenter radius on the microstructure and on the load-indentation depth response is examined revealing significant size effects. Results of such scope have not been reported so far.

The multiphase-field model
The finite-strain multiphase-field model is presented in this section. The present model is essentially an extension of the model of  to the case of multiple martensite variants, which has been achieved by employing the multiphase double-obstacle potential (Steinbach, 2009). In the computational treatment, the inequality constraints on the order parameters, which are essential in the double-obstacle potential, are enforced using the penalty method, which proves to perform very well. Additional enhancements of the model include the anisotropic elastic-strain energy that is quadratic in the elastic Hencky (logarithmic) strain, which makes the model more robust than in the case of the popular St. Venant-Kirchhoff model. The model assumes + N 1 phases, i.e. the parent phase (austenite) and N variants of the product phase (martensite). Each phase is characterized by the respective order parameter η i , = … i N 0, , . The order parameters represent the phase volume fractions and are subject to the following inequality and summation constraints, Note that the summation constraint (1) 2 implies also the fulfillment of the upper-bound constraints η i ≤ 1.
A finite-deformation framework is adopted, and the configuration of the undeformed austenite is taken as the reference configuration. Denoting by φ the mapping between the reference placement of the material point X and its current placement x at time t, thus = φ t x X ( , ), the deformation gradient F is defined as where u is the displacement, I is the second-order identity tensor, and ∇ denotes the spatial gradient with respect to the reference configuration. The deformation gradient F is multiplicatively decomposed into the elastic part F e and the transformation part F , t i.e. = F F F e t . The transformation part F t is obtained by applying the logarithmic mixing rule (cf.,  to the symmetric transformation (Bain) stretch tensors of individual phases U , i t viz.
is the vector of all order parameters η i . In view of the assumption concerning the reference configuration, we have = U I, 0 t and thus = U 0 log 0 t . In Eq. (3), the transformation stretch tensors are symmetric, In some cases, see e.g. the case with = N 4 martensite variants in Section 4.2, it is convenient to adopt a non-symmetric deformation tensor to describe the deformation associated with the transformation from the parent phase to the pure phase i, thus for = η 1 In such cases, the logarithmic mixing rule (3) is replaced by the linear mixing rule t . The Helmholtz free energy function F is additively split into the bulk contribution F B and the interfacial contribution F , Γ thus The bulk energy F B is adopted in the following form, where the first term represents the weighted sum of the chemical energies F i 0 and the second term represents the elastic strain energy. Here, H e is the elastic Hencky (logarithmic) strain, eT e is the elastic right Cauchy-Green tensor, with = − F F F ( ) , e t 1 and L is the elastic stiffness tensor. Applying the Voigt scheme, the elastic stiffness tensor is given by where L i represents the elastic stiffness tensor of the pure phase i, which can be different for each phase, in particular when elastic anisotropy is accounted for. Note that since the martensite variants are symmetry related, the anisotropic elastic stiffness tensor of one variant can be obtained by applying an adequate rotation to the elastic stiffness tensor of another variant.
The elastic strain energy in Eq. (5) has been assumed quadratic in the elastic Hencky strain H e . When elastic strains are sufficiently small, the simple and popular St. Venant-Kirchhoff model, in which the energy is quadratic in the elastic Green strain proves to perform satisfactorily (e.g. Kružík et al., 2005;Maciejewski et al., 2005;. However, according to our preliminary studies, the condition of small elastic strains is not satisfied in the range of applications studied in this work, in particular, in view of high compressive stresses beneath the indenter (leading to non-physical behavior and severe convergence problems). On the contrary, the energy quadratic in H e has proven to perform well, even if such behavior is not guaranteed for arbitrarily large elastic strains. The reader is referred to Neff et al. (2015) for an overview of the properties of the quadratic Hencky model in the isotropic case.
The interfacial contribution F Γ represents the energy of diffuse interfaces and is expressed in terms of the order parameters and their gradients. Adopting the standard double-obstacle potential, F Γ can be expressed in the following form (cf. Steinbach, 2009;2013), where γ ij is the interfacial energy density (per unit area) of the interface separating the phases i and j, and λ ij is the respective interface thickness. In the notation adopted, for each pair (i, j) of phases, the respective interface is represented by parameters γ ij and λ ij defined such that i < j. The interfacial energy F Γ involves the gradient of the order parameters η, and is accompanied by a homogeneous Neumann boundary condition ∇ = η 0 n on the entire boundary ∂B. The dissipation potential D of the viscous type is assumed in the following quadratic form, where m i are the mobility parameters. In order to examine how the mobility parameters m i , defined individually for each phase, are related to the mobility of the actual interfaces, consider the interface between two phases k and l, where The complete evolution problem is now formulated by following the variational framework developed by Hildebrand and Miehe (2012), see also Tůma et al. ( , 2018 for the treatment of the inequality constraints. The global potential energy functional is first introduced in the form, where is the total free energy and Ω is the potential energy of the external load. For instance, in the particular case when the nominal surface traction t* is prescribed over ∂B t , we have However, in this work, we have = Ω 0, and the load is applied through a contact interaction with the indenter. The details are discussed in Section 3.1. In the time-discrete (finite-step) setting, the following global (unconstrained) incremental potential Π τ is introduced, where the fields with subscript n denote the known solution at the previous time instant t n . For the sake of brevity, the subscript + n 1 indicating the fields at current time instant + t n 1 is omitted. In Eq. (10), τ is the global incremental dissipation potential that, upon applying the backward Euler method to the rate-potential (8), takes the following form, where = − + τ t t n n 1 denotes the time increment. The constraints to be enforced on the order parameters η i , cf. Eq. (1), are not included in the incremental potential Π τ . In order to introduce these constraints into the formulation, a suitable indicator function can be defined. Assume that is an arbitrary convex set in n . The indicator function → = ∪ +∞ I :¯{ } n of the set is defined as, Accordingly, the indicator function I corresponds to the admissible set of the order parameters, cf. Eq. (1), in the form of the standard simplex, is its global counterpart, Finally, the time-discrete constrained evolution problem is formulated as the minimization of the global constrained incremental potential Π* τ with respect to the fields of φ and η, The above compact formulation implies mechanical equilibrium (by minimization with respect to φ) and defines the time-discrete evolution of the order parameters η (by minimization with respect to η). The actual governing equations are discussed in Section 3.3, after the penalty regularization of the inequality constraints is introduced in Section 3.2. Note that the finite-step incremental formulation (15) can be derived from the corresponding rate formulation (Hildebrand and Miehe, 2012;, which is omitted here for brevity.

Contact problem
In the present simulations of nano-indentation, the external load is applied through frictionless contact interaction with the indenter. The indenter is assumed rigid, thus it is represented by a rigid surface denoted by Γ . Formulation of the corresponding contact problem is standard, and only the most important details are provided here, see e.g. Wriggers (2006) for an overview.
A part of the boundary of B, denoted by Γ , c is considered as the potential contact surface. For each point = φ x X ( ), ∈ X Γ , c the associated point x on Γ is found by the closest-point (orthogonal) projection. The normal gap g N is then defined according to where n is the unit normal to Γ at x . Frictionless contact interaction is then introduced into the formulation by enforcing the non-penetration constraint ≥ g 0 N on Γ c . The complete evolution problem including the contact interaction is then formulated as the following constrained incremental minimization problem and + denotes the set of all non-negative real numbers so that the indicator function

Penalty regularization
In the present computational scheme, the inequality constraints on the order parameters (η i ≥ 0) and on the contact normal gap ( ≥ g 0 N ) are enforced using the classical penalty regularization technique.
Recall that in the present model, austenite and N variants of martensite are considered. In the actual computational implementation, the volume fraction of austenite is not deemed an independent variable and, in view of the summation constraint (1) 2 , is defined as a function of the other order parameters = ⋯ η i N ( 1, , ) i . The system with + N 1 phases is therefore modeled using N independent order parameters, assembled in the condensed vector η, The obvious advantage of this treatment is that the total number of degrees of freedom is reduced, which is computationally beneficial. Upon introducing the penalty regularization, the minimization problem (17) is transformed to an unconstrained minimization problem where the penalty-regularized incremental potential Π τ pen takes the form where ϵ η > 0 and > ϵ 0 N are the penalty parameters associated with the order parameter constraints and with the contact constraint, respectively, and the following notation has been adopted, As it is well-known, the penalty parameters must be sufficiently large in order to avoid excessive violation of the constraints. At the same time, too large penalty parameters may lead to numerical problems, such as ill-conditioning of the tangent matrix and poor convergence of the Newton method. The related effects are studied in Section 4.4.

Governing equations in the weak form
Stationarity of the incremental potential Π , τ pen which is the necessary condition for the minimum of Π , τ pen yields the governing equations of the incremental evolution problem. The stationarity of Π τ pen with respect to the field of φ gives the standard weak form of the mechanical equilibrium, i.e. the virtual work principle, where φ δ is the test function that vanishes on the Dirichlet boundary ∂B φ on which φ is prescribed, and Here, P is the first Piola-Kirchhoff stress and t N is the normal contact traction. The formula (24) 3 for the variation δg N of the normal gap is a standard result in contact mechanics (e.g. Wriggers, 2006). On the other hand, the condition of stationarity of Π τ pen with respect to the field of η yields the evolution equation for η in weak form, viz.
where η δ^is the test function. The partial derivatives of the interfacial energy F Γ and the dissipation potential D τ are easily obtained in the explicit form, where for i > j the notation is adopted such that = γ γ ij ji and = λ λ ij ji . The term ∂F B /∂η i can also be readily obtained. However, for the sake of brevity, its explicit expression is not presented here. The reader may refer to Appendix A in , where the corresponding derivatives are provided in explicit form for the model with two hierarchical order parameters.

Remark 1.
The local form of the evolution Eq. (25) can be obtained in a standard manner by applying the Gauss theorem and by exploiting the homogeneous Neumann boundary conditions imposed on η. As a result, the local evolution equation is obtained in the form of the following Ginzburg-Landau equation, where η δ δ /^ is the functional derivative of the free energy functional , and the rate η˙is approximated by the backward-Euler formula (27) 2 , consistent with the incremental dissipation potential (11). In the notation employed above, the total free energy functional^ is expressed in terms of η and includes the penalty regularization of inequality constraints η i ≥ 0,

M. Rezaee-Hajidehi and S. Stupkiewicz
Mechanics of Materials 141 (2020) 103267 and M is a symmetric positive-definite mobility matrix such that the dissipation potential can be written as a quadratic function of η˙,

Finite-element implementation
The finite-element implementation of the model is briefly discussed here. The global unknowns of the problem are the fields of the displacement = − φ u X and the order parameters η. Four-node quadrilateral elements are used with piecewise-bilinear interpolation functions for all unknowns. Standard Gaussian quadrature is used for the numerical integration of the governing Eqs. (23) and (25).
In order to ensure a consistent approximation of the deformation gradient = ∇φ F and its transformation part and thus to avoid spurious stresses within the diffuse interfaces, the transformation part F t is considered constant within each element and is calculated at the element center, see . The matrix exponential, cf. the logarithmic mixing rule (3), and the matrix logarithm, cf. the Hencky strain in Eq. (5), along with their first and second derivatives are computed using the respective closed form representations, cf. Korelc and Stupkiewicz (2014), Hudobivnik and Korelc (2016).
The global coupled nonlinear equations that result from the finiteelement discretization are solved simultaneously with respect to all unknowns using the Newton method. The required tangent matrix is derived by linearizing the coupled equations using the automatic differentiation (AD) technique. The resulting exact linearization guarantees that the quadratic convergence rate of the Newton method is achieved. To this end, the AceGen system is employed, which provides an efficient implementation of the AD technique (Korelc, 2009;Korelc and Wriggers, 2016). For the computations, AceFEM has been used, which is a flexible finite-element environment that is fully integrated with AceGen. A direct linear solver (Intel MKL PARDISO) has been used in the computations. For the 2D problems considered in this study, with the largest problem of about 7.2 million unknowns, the direct linear solver proved to be more efficient than the iterative solvers available in the MKL library.

Preliminaries
In this section, the multiphase-field model presented above is used to study the microstructure evolution during nano-indentation. The present study is restricted to 2D plane-strain problems, and two corresponding transformations are considered, namely the square-to-rectangle transformation with = N 2 martensite variants and the squareto-parallelogram transformation with = N 4 martensite variants. The respective microstructures are examined in Sections 4.2 and 4.3. Next, in Section 4.4, a parametric study is performed with the aim to investigate the effect of different model parameters on the simulation results. Finally, in Section 4.5, predictions of the indentation size effect (ISE) are presented. Specifically, the influence of the indenter radius on hardness and on microstructure evolution is investigated.
According to the classical crystallographic theory of martensite, a compatible interface between stress-free austenite and a single martensite variant does not exist for the majority of materials undergoing martensitic transformation, and compatibility is usually achieved through twinning (Bhattacharya, 2003). However, in 2D, the kinematic compatibility condition is automatically satisfied for isochoric or nearly isochoric transformations, and thus a single martensite variant may form a stress-free interface with austenite (the corresponding condition for compatibility is that one eigenvalue of the transformation stretch is less than unity and the other one is greater than unity).
In view of the essential differences between the martensitic microstructures in 2D and in 3D, as discussed above, the analysis of 2D transformations is necessarily associated with some simplifications. Nevertheless, such an analysis constitutes an important intermediate step in developing a full 3D model (which is the subject of our ongoing work). It is also believed that the analysis of 2D problems may provide valuable results of general interest, for instance, concerning the size effects, as studied in Section 4.5.
One way of interpreting the results obtained for 2D transformations is to treat a single martensite variant in 2D as a so-called habit-plane variant of a 3D transformation. The austenite-martensite interface would then correspond to the austenite-twinned martensite interface and the martensite-martensite interface would correspond to the interface between two twinned martensites. In the light of this interpretation, the corresponding interfacial energies are here adopted higher than those of the atomic scale interfaces, as they are assumed to include the energy of elastic micro-strains (cf. Maciejewski et al., 2005;Petryk et al., 2010).

Microstructure evolution: square-to-rectangle transformation
The first numerical example concerns a preliminary study of the microstructure evolution under nano-indentation for the square-torectangle transformation, i.e. for the system with austenite and = N 2 martensite variants, which is the simplest transformation in 2D. The two martensite variants are characterized by the transformation stretch tensors where the values of the stretch parameters = α 0.95 and = β 1.05 have been adopted such that U 1 t and U 2 t correspond to a nearly isochoric transformation. Recall that for the austenite we have = U I 0 t . According to the crystallographic theory of martensite, the austenite-martensite interfaces and the martensite-martensite interface in this transformation are oriented at approximately ± 45 ∘ , with respect to the coordinate system aligned with the square lattice of austenite.
The material parameters are taken as follows. The interfacial energy density is adopted equal to = = γ γ 0.4 i (MPa s) −1 has been adopted, which provides reasonable results for a realistic indentation speed of = v 5 nm/s, see below. The interface thickness parameters are assumed equal for all interfaces, = = λ λ 12 ij nm. Note that the finiteelement size must be sufficiently smaller than the interface thickness λ so that the mesh can correctly resolve the diffuse interfaces. Here, the element size is taken approximately equal to λ/7.
Finally, the penalty parameters = ϵ 1000 η GPa and = ϵ 1000 N GPa/ nm are chosen large enough to enforce the respective inequality constraints with an adequate accuracy and, at the same time, not to affect the computational efficiency, see the related discussion in Section 4.4. x y nm 2 is considered. The vertical displacement of the bottom boundary and the horizontal displacement of the lateral boundaries are constrained. A rigid circular indenter with radius = R 50 nm is considered, and indentation is applied at the center of the top surface at a constant speed = v 5 nm/s up to the maximum indentation depth of = h 50 max nm. Recall that contact is frictionless. In this preliminary study, unsymmetric microstructures are excluded by simulating one half of the rectangular domain with the adequate symmetry conditions imposed on the unknowns along the symmetry axis. The computational domain is discretized with 168 500 quadrilateral elements, leading to the total number of degrees of freedom of approximately 674 000. Fig. 1c shows the detail of the deformed finite-element mesh in the vicinity of the indenter.
The snapshots of the microstructure evolution at selected indentation depths are shown in Fig. 2b. The red marks superimposed on the load-indentation depth (P-h) curve in Fig. 2a indicate the instants at which the snapshots are taken. The martensite variants 1 and 2, denoted by V1 and V2, respectively, are represented by the red and blue domains that correspond to η 1 ≥ 0.5 and η 2 ≥ 0.5, respectively. The austenite, which occupies the remaining part of the computational domain, is not shown explicitly. Grey lines representing the interfaces have been laid over the microstructures for a clear distinction of the phase boundaries, but they do not represent the actual diffuseness of the interfaces (the actual diffuseness of the interfaces is visible in Fig. 1b). Note that in this example only one half of the domain has been simulated, but for a better visualization of the microstructures in Fig. 2, the other half has been replicated by exploiting the symmetry.
The transformation initiates with the nucleation of martensite variant V2 at the indentation depth of about = h 7 nm. The preference towards the formation of variant V2 is due to the fact that a compressive stress is applied beneath the indenter, which complies with the transformation stretch tensor U 2 t . As the indentation proceeds, the domain of variant V2 grows and eventually, at the indentation depth of about = h 22 nm, variant V1 appears. Further increase of the indentation depth leads to the growth of both variants, while the shape of the respective domains changes. Note that the orientation of the actual interfaces is close to ± 45 ∘ , as predicted by the crystallographic theory.
During unloading, for which the corresponding snapshots are not reported here, the two martensite variants shrink simultaneously, while preserving the final microstructure, i.e. the one that corresponds to = h 50 nm. Only at the final stage, starting from = h 10 nm, the microstructure changes and the evolution follows that of the loading stage in the reverse order. At the indentation depth of about = h 7 nm, variant V1 vanishes completely followed by the disappearance of variant V2 at the indentation depth of about = h 3 nm. A movie showing the complete microstructure evolution is provided as a supplementary material 1 accompanying this paper.

Microstructure evolution: square-to-parallelogram transformation
As the next example, the microstructure evolution under nano-indentation is studied for the square-to-parallelogram transformation, i.e. for the system with austenite and = N 4 martensite variants. As discussed in Section 2, in order to describe the deformation associated with this transformation, it is convenient to adopt non-symmetric deformation tensors F i t . Then, the transformation deformation gradient F t within diffuse interfaces is obtained by applying the linear mixing rule, t . In the present case, the following deformation tensors are used to describe the four martensite variants, where the stretch parameters = α 0.95 and = β 1.05 are taken as those in the previous example and the shear parameter = γ 0.05 is assumed. The remaining model parameters are the same as those in the previous example, c.f. Section 4.2.
The geometry of the problem is the same as in the previous example, see Fig. 1a, except that the computational domain of the size × = × L L 1000 1000 x y nm 2 is considered. In this case, unlike in the previous example, the symmetry is not exploited and the simulation is carried out for the full domain. Accordingly, non-symmetric microstructures may develop, which indeed is observed in some cases, as illustrated in Sections 4.4 and 4.5. The radius of the rigid circular indenter = R 50 nm and the indentation speed = v 5 nm/s are set equal to those in the previous example. The indentation continues up to the maximum indentation depth of = h 30 max nm. The computational domain is discretized by keeping the element size equal to λ/7, as in the previous example, giving the total number of about 300 000 elements and 1 800 000 degrees of freedom.
Four representative snapshots of the microstructure evolution are shown in Fig. 3 along with the P-h curve. Each martensite variant i, denoted by Vi, is represented by the respective colored domain that corresponds to η i ≥ 0.5. Again, the domain of the austenite is not shown explicitly.
All martensite variants nucleate nearly at the same indentation depth of approximately = h 4 nm and subsequently grow together as the indentation proceeds. The same is also true during unloading, where all the variants shrink and annihilate simultaneously. A movie showing the complete microstructure evolution is provided as a supplementary material 2 accompanying this paper.
In Fig. 4, a comparison has been made between the orientations of the interfaces in the simulated microstructure and those predicted by the crystallographic theory (see the orientations shown in the circles included in Fig. 4). It can be seen that the orientations predicted by the phase-field model are in a good agreement with the theoretical ones. Even a better agreement can be observed when the microstructure occupies a larger domain, see the microstructure for = R 100 nm in Fig. 9. Parts of the fully developed microstructure at = h 30 nm, see Fig. 3b, apparently resemble the wedge-like microstructures considered by Bhattacharya (1991). Note, however, that these are not stress-free wedge-like microstructures. It can be checked that they do not satisfy the corresponding compatibility conditions (Bhattacharya, 1991), even if the orientations of the individual interfaces are close to the theoretical orientations of stress-free interfaces, as illustrated in Fig. 4.
The P-h curve displays a sudden load drop at the instant of nucleation of martensite. The degree of abruptness as well as the change of the slope of the P-h curve in the post-nucleation stage are much higher than those in the previous example, see Fig. 2a. It has been observed in a preliminary study that, as the value of the shear parameter γ increases, a sharper nucleation occurs, and the slope of the P-h curve in the post-nucleation stage deviates more from the elastic one. The square-to-parallelogram transformation involves four variants of martensite, as compared to two variants in the square-to-rectangle transformation, so that more complex microstructures may develop with higher flexibility for accommodation of possible incompatibilities. 3 Clearly, the corresponding effects are more pronounced as γ is increased. The studies carried out in the subsequent sections are limited 1 http://bluebox.ippt.pan.pl/~sstupkie/files/SuppSq2Rec.gif 2 http://bluebox.ippt.pan.pl/~sstupkie/files/SuppSq2Par.gif 3 For similar reasons, the NiTi shape-memory alloy, which undergoes the cubic-to-monoclinic transformation with 12 martensite variants, is superior to other alloys, in particular, to those undergoing the cubic-to-tetragonal transformation with only three variants.
to the square-to-parallelogram transformation.

Parametric studies
In this section, we examine whether and how the microstructure and the P-h curve are sensitive to the choice of selected material and process parameters. In particular, the effect of the chemical energy and the effect of the indentation speed are investigated for the system with = N 4 martensite variants. In addition, to assess the accuracy of the finite-element solution and the efficiency of the numerical scheme, the effect of the penalty regularization parameter ϵ η on these characteristics is examined. The effect of the contact penalty parameter ϵ N is not included in this study, since the corresponding effect is very well known in computational contact mechanics.
First, the effect of the chemical energy is studied. To this end, the computation has been additionally carried out for two values of the chemical energy, namely does not reach the boundary of the computational domain. It is immediate to see that the higher the chemical energy, the higher the load at the nucleation event and the higher the slope of the P-h curve in the transformation stage.
The other major effect of the chemical energy is related to the change in the microstructure evolution, see   load plateau, as indicated by the arrow in Fig. 5a. Furthermore, at the indentation depth of = h 10 nm, the separation of the indenter from the surface occurs. Hereafter, the reverse transformation proceeds at zero load with no considerable change in the microstructure, and ultimately the martensite domains annihilate.
Next, the effect of the indentation speed is examined. The simulation has been performed for two indentation speeds, namely = v 0.5 nm/s and = v 50 nm/s (in the sequel, referred to as 'slow' and 'fast', respectively), and the results are compared with the reference case of = v 5 nm/s. Note that exactly the same results can be obtained by varying the mobility parameter m, while keeping the indentation speed v constant. For instance, the slow case with = v 0.5 nm/s and = m 1 (MPa s) −1 is equivalent to the case with = v 5 nm/s (as in the reference case) and = m 10 (MPa s) −1 . Fig. 6 shows the snapshots of the microstructure evolution for the three indentation speeds. The general pattern is not much influenced, although the details are different. During loading, for a fixed indentation depth, the size of the transformed domain decreases with increasing indentation speed. The effect is opposite during unloading. In the slow case, the driving force for interface propagation is close to zero and the microstructure is thus close to the equilibrium microstructure. As the loading rate is increased, the driving force increases in accordance with the viscous evolution law, and the microstructure grows and shrinks (during loading and unloading, respectively) with a delay with respect to the equilibrium microstructure, which explains the observation concerning the size of the transformed domain.
Moreover, in the fast case, the separation of the indenter during unloading occurs at the indentation depth of approximately = h 10 nm. Henceforth, due to the release of the external load, a relaxing microstructure develops in the system, such that additional martensite domains appear, similar to the case of = F 0 m 0 in Fig. 5b. Another interesting observation is that the microstructure is not symmetric in the slow case. At the initial stage of the transformation, the domains grow in a symmetric manner, however, later, the symmetry is broken and further evolution during loading and unloading proceeds in an unsymmetric manner, see Fig. 6. This effect is not observed in the reference and fast cases, and the respective microstructures evolve in a symmetric manner during the complete forward and reverse transformation. It has been checked that the non-symmetric microstructure is energetically more favorable than the symmetric one, i.e. at the instant of the symmetry breakdown the incremental energy is lowered thanks to non-symmetric evolution, and this is captured by the present incremental energy minimization framework. It has also been confirmed that the symmetry breakdown is not a numerical artifact. The simulations of the slow case have been repeated several times with different time stepping, and the symmetry breakdown has been observed in all those simulations.
The effect of the indentation speed on the P-h curve is depicted in Fig. 7. As expected, the decrease of the indentation speed leads to the decrease of the maximum indentation load and also to the decrese in the hysteresis width. The resulting rate-dependent response is evidently the outcome of the assumed rate-dependent viscous-type dissipation, c.f. Eq. (8). Note that with further reduction of the indentation speed the dissipated energy does not decrease to zero. As discussed by Tůma et al. (2018), this is because the nucleation event, associated with a sudden load drop, proceeds dynamically with a non-zero local strain rate that is independent of the indentation speed and is associated with non-zero dissipation. Fig. 7 illustrates also the effect of the stiffness of the loading device. So far, the simulations have been performed by assuming an infinite stiffness of the loading device, i.e. the position of the indenter was controlled directly by prescribing a constant indentation speed. However, in practice, the indentation machine possesses a finite stiffness that may affect the load-indentation depth response, due to the release of the accumulated elastic energy at the instant of nucleation. In order to study the impact of the loading-device stiffness, additional simulations have been performed in which the position of the indenter is controlled through an elastic spring with the stiffness K, see the inset in  The case of the infinitely stiff loading device is denoted by = ∞ K . For a finite, but very high stiffness = K 1 GPa, the predicted response is essentially identical to that corresponding to = ∞ K and is not included in Fig. 7. Note that the indentation load and the spring stiffness are referred to a unit thickness in the out of plane direction, and hence the unit of the stiffness is 1 Pa. Upon reducing the stiffness K, the sudden load drop at the instant of nucleation is replaced by a gradual transition zone in the P-h curve. The corresponding response resembles the pop-in effect, which is commonly observed in nano-indentation of materials that deform through plastic slip. Pop-in has also been observed in materials undergoing martensitic transformation, but the corresponding experimental results are much more scarce (e.g. Caër et al., 2013;Dar and Chen, 2017).
Finally, the effect of the penalty regularization parameter ϵ η , c.f. Eq. (21), on the accuracy and efficiency of the computational scheme is examined. To this end, the simulations have been performed for several values of the penalty parameter, = ⋯ ϵ 10, 10 , , 10 η 2 6 GPa. Fig. 8a shows the profile of the order parameters η i along a horizontal line at the distance of 10 nm from the contact surface (taken in the reference configuration), where the stresses are the highest and thus the violation of the bound constraints 0 ≤ η i ≤ 1 is more severe. It follows that for = ϵ 10 η GPa (dashed lines in Fig. 8a) the constraints are significantly violated. For = ϵ 10 η 2 GPa (dotted lines in Fig. 8a), the violation is visible, but small (of the order of 0.01). With further increase of ϵ η , the violation becomes insignificant. The results corresponding to = ϵ 10 , 10 η 4 5 and 10 6 GPa are barely distinguished from those for = ϵ 10 η 3 GPa and thus are not included in the figure. The impact of the penalty parameter ϵ η on the P-h curve is qualitatively similar to that revealed in Fig. 8a, and the corresponding results are thus not shown here.
The bar chart in Fig. 8b shows the simulation wall-clock time and the number of time steps needed to complete the simulation. Note that an adaptive time-stepping algorithm has been employed in the simulations. Quite surprisingly, the computational efficiency is not significantly affected when the penalty parameter is varied between 10 2 and 10 6 GPa. Considering that the penalty parameter = ϵ 10 η 3 GPa guarantees satisfactory accuracy, as discussed above, this value has been used for all simulations.

Size effects
In this section, the effect of the indenter radius R on the microstructure and on the P-h curve is studied for the system with = N 4 martensite variants. In order to quantitatively interpret the related size effects, the indentation hardness H is calculated at the maximum indentation depth according to where P max is the maximum indentation load and A is the corresponding projected contact area, here both quantities are referred to unit thickness in the out-of-plane direction. In the present computations, the finite-element mesh has been refined in the vicinity of the contact surface. This does not visibly influence the load-indentation depth response, but it improves the accuracy with which the area A and thus the hardness H are determined.
The simulation is performed for the indenter radius R varied between 25 nm and 800 nm. To make the comparisons meaningful, the ratio between the maximum indentation depth h max and the indenter radius R is kept constant. Furthermore, in order to preserve the geometrical similarity, the ratio between the indenter radius and the dimensions of the computational domain, i.e. R/L x and R/L y , is kept constant for all simulations. This, in particular, ensures that upon proper normalization the elastic response is identical in all cases. In addition, to avoid mesh bias, the element size is also kept equal to λ/7. Thus, for a fixed interface thickness parameter λ, upon increasing the indenter radius, and accordingly increasing the size of the computational domain, the number of elements increases, and the computational cost increases (the largest computation in this section involves about 1.2 million elements and 7.2 million degrees of    9. The effect of the indenter radius R on the microstructure evolution for the system with = N 4 martensite variants with = λ 12 nm. The spatial dimensions have been normalized with respect to the indenter radius R (note the scale makers at the bottom). freedom). This restrains the computations to a limited range of the indenter radius. To evade this limitation and to study the size effects for a wider range of the indenter radius, several values of the interface thickness, i.e. = λ 24, 48 and 96 nm, are used in addition to the reference case of = λ 12 nm. The element size (kept equal to λ/7) is then increased accordingly, thus reducing the number of elements for a given size of the computational domain. Upon increasing the interface thickness parameter λ by the factors of 2, 4 and 8, the mobility parameter m must be reduced by the respective factors, so that the effective mobility of the interfaces is not influenced, see the related discussion in Tůma et al. (2018). Note that λ is bounded from above, since too large values of λ result in microstructures with excessively diffuse interfaces.
The effect of the indenter radius R on the microstructure evolution for the case with = λ 12 nm is depicted in Fig. 9. Similar to the effect of the indentation speed in Fig. 6, the general pattern of the microstructure is not affected much, however, some small differences are noticeable. For the cases with = R 25 nm and = R 50 nm, the microstructure evolves in a symmetric manner. On the other hand, for the case with = R 100 nm, the symmetry of the microstructure breaks at the final stage of loading, and then the microstructure evolution proceeds in an unsymmetric manner, see the related discussion in Section 4.4.
The microstructure evolution is also influenced by the interface thickness parameter λ. For a fixed indenter radius R, increasing λ results in microstructures with more diffuse interfaces. Moreover, for some cases, it has been observed that increasing λ induces some minor changes in the microstructure pattern, more specifically, small martensite domains appear in the vicinity of the indenter (see e.g. the martensite domains that appear in the vicinity of the indenter in the case of zero chemical energy, Fig. 5b). The microstructures obtained for = λ 24, 48 and 96 nm are not presented here for the sake of brevity. The reader is referred to , where a comprehensive study has been performed to investigate the effect of the interface thickness parameter on the microstructure of the austenite-twinned martensite interface. Fig. 10a shows the effect of the indenter radius R on the P-h curve. It can be seen that, as the indenter radius R decreases, the transformation initiates at a lower indentation load P. This is because, for the same indentation load, a smaller indenter radius leads to higher stresses beneath the indenter. A slight difference can also be noticed in the slope of the curves in the transformation stage.
In order to capture the size effect in the P-h curves, in Fig. 10b the load P and the indentation depth h are normalized by the indenter radius R. In the elastic stage, the normalized P-h curves coincide. In the transformation stage, however, a size-dependent response is observed.
As the indenter radius R decreases, the normalized load, P/R, required to initiate the transformation increases (contrary to the actual nucleation load P in Fig. 10a), and subsequently the transformation proceeds at a higher load. Fig. 11 shows the dependence of the hardness H on the indenter radius R for a fixed normalized indentation depth = h R / 0.5. The general trend is that hardness increases with decreasing indenter radius R. This size effect is consistent with the well-known indentation size effect observed in metals that deform through plastic slip. However, the mechanism is here different. In plasticity, the indentation size effect is attributed to the geometrically necessary dislocations (GNDs), the density of which increases as the physical dimension of the deformed zone decreases with decreasing indentation depth (Nix and Gao, 1998;Pharr et al., 2010). In martensitic transformation, the size effect is governed by the interfacial energy. The inelastic deformation proceeds by formation and evolution of microstructure. As the indenter radius R decreases (and so does the indentation depth h for a fixed ratio of h/R), the size of the transformed domain decreases, and the contribution of the energy of the interfaces to the total free energy increases. As a result, at a smaller scale, a relatively higher load, i.e. a higher hardness, is needed to induce the interfaces. The related effects are correctly represented by the phase-field model.
Due to the reasons discussed earlier, the hardness H cannot be determined in the whole range of R for a fixed interface thickness Fig. 10. Indentation size effect: the effect of the indenter radius R on (a) the P-h curve and (b) the normalized P-h curve for the case with = λ 12 nm. Fig. 11. Indentation size effect: the dependence of the hardness H on the indenter radius R for the fixed ratio of = h R / 0.5.
parameter λ. As shown in Fig. 11, parameter λ influences the results. However, the related effect is moderate and the general dependence of H on R can be deduced from Fig. 11. The effect of the interfacial energy manifests itself also in the shape of the microstructure, see Fig. 9. When the transformed domain is relatively large, the interfaces are mostly planar with only small portions of high curvature. The orientations of the planar interfaces are then close to those predicted by the crystallographic theory, i.e. the orientations are governed by the kinematic compatibility and by the related elastic strain energy. On the other hand, when the transformed domain is relatively small, the interfacial energy becomes dominant, thus resulting in curved interfaces. The related effects are visible in Fig. 9 when comparing the microstructures for a fixed indenter radius R and varying the normalized indentation depth h/R, and for a fixed normalized indentation depth h/R and varying the indenter radius R.

Conclusions
A finite-strain phase-field model of multivariant martensitic transformation has been developed, implemented and successfully applied to simulate spatially-resolved microstructure evolution during nano-indentation. In view of considerable number of global unknowns, which comprise displacements and order parameters, the analysis is limited to 2D problems. Extension to 3D problems, which leads to demanding large-scale computations, is in progress.
The computational study focused on pseudoelastic response of SMAs under nano-indentation has revealed several interesting effects. We are not aware of any published results of similar scope. Firstly, a significant indentation size effect has been predicted which is governed by the interfacial energy. The increase of hardness with decreasing indenter radius (this effect is often referred to as 'smaller is stronger') results from the related increase of the contribution of the interfacial energy to the total free energy. The mechanism is thus different from that in plasticity, but the overall effect is similar. Secondly, nucleation of martensite may be associated with a load drop (for a stiff loading device) or with a sudden increase of the indentation depth (for a compliant loading device). The latter behavior resembles the pop-in effect that is often observed in nano-indentation testing.
It has also been observed that in some cases (e.g., for slower indentation or for larger indenter radius) the symmetry of the microstructure is broken. It has been checked that the development of the non-symmetric microstructure is energetically preferable, which is captured by the present computational scheme. Note that the computational scheme is based on the incremental energy minimization approach, however, the actual governing equations are derived from the condition of stationarity of the incremental potential, which is only a necessary condition for the minimum.
One of the features of the present computational model is that the inequality constraints imposed on the order parameters are enforced using the penalty method. Explicit treatment of the inequality constraints is necessary because the model employs the double-obstacle potential and the computational treatment is based on an implicit monolithic scheme. A parametric study has shown that the penalty parameter can be varied within a very wide range of values without visible impact on the overall efficiency of the computational scheme, while the constraints are enforced with a satisfactory accuracy already for a moderate penalty parameter. The penalty regularization proves thus to perform very well for the problem at hand.