Modeling the indentation size effects of polymers, based on couple stress elasticity and shear transformation plasticity

This work attempts to construct models which can describe the size effects in modulus and hardness of polymers measured by indentation tests. Firstly, the origin of size effects in the procedure of nanoindentation is analyzed. The finding of the literature that the indentation size effects of polymers are of significant elastic nature is further discussed and confirmed. The elastic size effects are described by a model, through introducing a model of elastic unloading load with consideration of couple stress elasticity into the Oliver-Pharr indentation approach. The accordingly proposed modulus model and hardness model agree excellently with a large amount of experimental data obtained from literatures. The models show that the elastic size effects of polymers and their experimental observations are mainly determined by the molecular structures. The fitting results verify that the size effects in indentation hardness of polymers with complex molecular structures are significantly elastic. It is postulated that the plastic size effects in indentation hardness of polymers are only derived from their glassy components. A shear transformation plasticity formula of glassy polymers recently proposed by literature is slightly extended to characterize the size dependent deformations in indentation tests. A hardness model with consideration of size effects is accordingly obtained and agrees well with related experimental data.


Introduction
The nano-indentation approach of Oliver and Pharr [1] has been widely used to measure modulus and hardness of crystals, polymeric materials or thin films. In such a test, the phenomenon that the measured modulus and hardness increase with decreasing indentation depths is called indentation size effects (simplified as ISEs). For crystals, size effects almost do not emerge in indentation modulus (simplified as modulus ISEs), but mainly in indentation hardness (simplified as hardness ISEs), and have been well studied and mostly described by the model proposed by Nix and Gao [2]. Nix and Gao found ISEs of crystals are induced by the different distribution of geometrically necessary dislocations along the indentation depths. For polymers, size effects are however extensively observed in both indentation modulus and indentation hardness [3][4][5][6][7][8][9][10]. But the notion of dislocations cannot be applied to polymers due to the lack of long-range order. It is necessary to investigate the essential mechanism of the size effects of polymers, from which we may get some favorable results conducive to the development of MEMS, sensors, precision instruments.
The present work attempts to construct models which can describe the ISEs of polymers. The first study in this respect is the hardness model of Lam and Chong [11] which adopts the notion of the geometrically necessary kinks of molecular chains in polymer plasticity. This is similar to the innovation of Nix and Gao C. Peng · F. Zeng (B) Department of Astronautic Science and Mechanics, Harbin Institute of Technology, Harbin 150001, People's Republic of China e-mail: zengfanlin@hit.edu.cn [2] which adopts the notion of the geometrically necessary dislocations in crystal plasticity. A large number of researches have revealed that the elastic size effects in deformation of solids can be interpreted by couple stress elasticity theory [12][13][14][15][16][17][18][19]. In order to employ the couple stress elasticity to construct the ISEs model of polymers, Han and Nikolov [20] simplified the indentation problem as the Boussinesq problem. The couple stress elasticity is introduced into the simplified problem and a hardness model is accordingly proposed. Both the model of Lam and Chong [11] and that of Han and Nikolov [20] involve rotation gradients in deformations of microstructures. Voyiadjis et al. [10] suggested that the size effects of polymers may be at least partly independent of rotation gradients, since size effects were also found in rotation-free situations such as uniaxial tensile experiments of polymeric nanofibers. They accordingly proposed a hardness model, by extending a shear transformation plasticity theory of homogenous materials developed by Voyiadjis and Samadi-Dooki [21], Spaepen [22] and Argon [23].
The present work is an extension of the studies of Han and Nikolov [20] and Voyiadjis and Samadi-Dooki [21]. In Sect. 2, the Oliver-Pharr indentation method is briefly recalled and the origin of the ISEs is discussed. In addition, large amount of experimental data is comparatively analyzed, which supported the finding of Alisafaei et al. [24] that the ISEs of polymers are of significant elastic nature. In Sect. 3, the elastic ISEs of glassy polymers are accordingly modeled. The contact problem of a rigid indenter on polymers is complicated because one need to consider the effects derived from the viscoelasticity. Zhang et al. [25,26] reported some models to solve this problem on different layered polymeric materials, but the focus was not on ISEs. Before we attempted to solve the spherical contact problem in the context of couple stress elasticity a number of solutions regarding similar problems have been reported [27][28][29][30][31][32]. Based on these works, similar to the work of Han and Nikolov [20], the couple stress elasticity was introduced into Hertz contact problem [33], and a theoretical model of contact load which can characterize the size effects in contact tests of elastomers was derived in our previous work [34]. Then the proposed load model is introduced into Oliver-Pharr approach. Accordingly, two models, respectively, characterizing the modulus ISEs and the elastic hardness ISEs are derived. These two models are applied to a large number of experimental data to illustrate their validity. In Sect. 4, the shear transformation plasticity theory proposed by Voyiadjis and Samadi-Dooki [21] is briefly recalled. In order to employ their plasticity formula to characterize the size dependent deformation in indentation, a slight extension of the plasticity formula is conducted. The accordingly proposed hardness model contains the indentation depths. A good agreement between the hardness model and related experimental data is obtained. For simplicity, we suppose the contact process is static or quasi static, thus the viscoelasticity effect was not considered here.

The elastoplastic nature of the ISEs of polymers
We start by recalling the indentation approach proposed by Oliver and Pharr. Their approach adopts a diamond Berkovich indenter. A single indentation cycle contains three periods: loading to a peak load P max , as shown in Fig. 1a; holding for a period of time at the peak load; unloading from the peak load to zero with a plastically residual depth h f being left, as shown in Fig. 1b. The holding period is conducted in order to minimize any non-elastic effects upon unloading period, assuring a completely elastic deformation with respect toh − h f . Figure 1c shows a typical load-curve of a single indentation cycle. Three important quantities are recorded: the maximum deformation depth h max , the peak load P max and the plastically residual depth h f .
Oliver-Pharr approach starts with the elastic contact analysis at the initial unloading stages. They found that the elastic unloading curve can be best fitted by a power law relation: where P is the elastic contact load at the unloading stage, α and m are fitted parameters and h f is exactly the plastic residue. The critical property estimated from the unloading curve is the elastic contact stiffness at initial unloading stages, which is defined as The contact depth h c shown in Fig. 1a is then estimated by The projected area of elastic contact shown in Fig. 1a is determined by substituting h c into the area function of the cross section of indenter, i.e., A c 24.5h 2 c (4) and the elastic modulus is related to A c by where β 0 is a dimensionless parameter accounting for indenter geometry. Finally, the hardness is determined by The critical advance in Oliver-Pharr approach is their understanding of the effective behavior of Berkovich indenter at the initial unloading stages. Firstly, on the basis of a large number of experimental observations, they found the elastic unloading curves at the initial stages can be exactly represented by Eq. (1), in which the values of m are around 1.5. The contact load of linear elastic half space by a rigid paraboloid of revolution has been analytically expressed by Sneddon [35] as a power law with the exponent being exact 1.5. Therefore, Oliver and Pharr suggested the effective behavior of the Berkovich indenter in contact with the elastic recovery of deformed surface can be approximated by the behavior of a paraboloid of revolution in contact with the flat surface of half space, as schematically shown in Fig. 2. Further, on the basis of several finite element simulations, Pharr and Bolshakov [36] verified this inference and constructed the concept of "effective indenter shape".
In Sneddon's work, the contact load of linear elastic half space by a rigid paraboloid of revolution is expressed as where P 0 is the elastic contact load, μ is the shear modulus of half space, R is the curvature radius of indenter tip, v is Poisson ratio, h e is the normal displacement of the contact center. In the limit of small displacement, this contact problem is the same as the contact of half space by a rigid sphere, which is also the famous Hertz contact problem [33]. As indicated in Eq. (7), one can obtain the modulus of the material by employing a rigid sphere to contact with it in the limit of small displacement. During this kind of tests, the contact load P 0 and the elastic contact   [37] versus h e , in which obvious size effects is presented within h e < 10 μm. The size effect model Eq. (9) is correspondingly applied depth h e are recorded. The modulus data can be readily estimated by Eq. (7) with known R and v. Recently, Han et al. [37] have conducted the contact tests on PDMS, and the results are plotted here in Fig. 3 (the scattering points). Figure 3 shows that the modulus of PDMS estimated by Eq. (7) increase with the decreasing contact depths h e , i.e., obvious size effects are presented. This kind of size effects were also observed in the contact tests of various polymers conducted by other authors [38][39][40]. The radius R and Poisson ratio v are constant in a specific test. Therefore, according to Eq. (7), the size effects are probably transferred by the contact load P 0 and contact depth h e from the elastic resistance. Specifically, in fact a contact load larger than the P 0 decided by Eq. (7) is recorded at each contact depth, which indicates that the contact load itself contains size effects.
Recently, in order to characterize the size effects observed in contact tests of polymers, we investigated the Hertz contact problem in the context of couple stress elasticity [34]. A model of contact load with consideration of couple stress effects is analytically derived and expressed as where P is the couple stress-based contact load, l is a characteristic length of deformation body employed in couple stress elasticity to characterize size-dependent deformations, similar to the two Lamé constants in classical linear elasticity. It can be seen from Eq. (8) that the couple stress-based contact load P contains an augmenting term relative to the classical spheric contact load P 0 . Besides, the augmenting effects increase with the decreasing deformation depth h e , which denotes the size effects existing in contact load. Incorporating Eq. (7) into Eq. (8), we can obtain a modulus model: where E 0 denotes the constant modulus theoretically estimated by Eq. (7). This model has been successfully used to characterize the size effects found in spherical contact test of PDMS [37]. As shown in Fig. 3 (the solid line), a good agreement between the model and the experimental data is obtained (see Ref. [34] for details).
According to the Oliver-Pharr approach, during the initial unloading stages, the effective contact between the Berkovich indenter and the elastic recovery of deformed surface is equivalent to that between a paraboloid of revolution and the flat surface of half space. The above discussion has indicated that the modulus measured by spheric contact tests also contains size effects and these size effects are found to derive from the contact load. Note that the indentation modulus in Oliver-Pharr approach is directly decided by elastic contact stiffness S and contact depth A c which are ultimately decided by the contact load P and contact depth h c . Therefore, it is reasonable to consider that the modulus ISEs probably also derive from the elastic contact load during the initial unloading stages. It can be also inferred that the hardness ISEs derive from the peak load P max . Specifically, in indentation tests, the recorded elastic contact load and peak load at each contact depth are, respectively, larger than their counterparts theoretically decided by the fundaments of Oliver-Pharr approach.
Recently, Alisafaei et al. [24] discussed the elasto-plastic nature of the size effects of polymers. The modulus ISEs are of only elastic nature, since the measurement of modulus in Oliver-Pharr approach is only associated with elastic deformation. Through comparing the indentation modulus data of metals and polymers, they found that the modulus ISEs are extensively observed in polymers whereas almost no modulus ISEs is expected in metals since size effects in metals are associated with plastic deformation. Alisafaei et al. accordingly suggested that the size effects of polymers are dominated by elastic deformation, specifically, the hardness ISEs of polymers are of significant elastic nature. In addition, they strengthened the inference since it is found that for a polymer the increase of indentation modulus is approximately proportional to the increase of indentation hardness with decreasing depths.
It can be seen from these two figures that for PS, PMMA, Epoxy (Alisafaei et al. [8]), Nylon66 and LDPE the hardness ISEs are always accompanied by the modulus ISEs, while for UHMWPE and PTFE the ISEs emerge in neither modulus nor hardness. In other words, for these polymers, ISEs would emerge in hardness as long as they emerge in modulus, or vice versa. As discussed above, the modulus ISEs probably derive from the elastic contact load during the initial unloading stages. Additionally, the reversible elastic contact load P mediating modulus ISEs is a part of the peak load P max utilized to estimate hardness. Therefore, for a polymer if the contact load mediates modulus ISEs, the same elastic size effects would emerge in the hardness. Among the polymers displayed in Figs. 4 and 5, none has only hardness ISEs but no modulus ISEs. Thus, it can be considered that for these polymers the hardness ISEs is dominated by the elastic part. It can even be observed that for those polymers with size effects their modulus and hardness have similar shapes of increase tendency. This kind of similarity is consistent with the approximately proportional relation suggested by Alisafaei et al. [24] between their modulus ISEs and hardness ISEs.

Model development
The above discussions indicate the ISEs of polymers are of significant elastic nature. If the plastic hardness ISEs are temporarily neglected, the modulus ISEs and the elastic hardness ISEs may be theoretically modeled by introducing the couple stress-based elastic contact formula Eq. (8) into the Oliver-Pharr indentation approach. Specifically, the elastic contact load at the initial unloading stages, i.e., Eq. (1), may be decomposed into two parts: an elastic load and a corresponding augment mediated by couple stress effects, similar to Eq. (8). For simplicity, the elastic part in Eq. (1) is assumed to be in the same analytical form as Eq. (7), with the exponent m in Eq. (1) being sustained. Besides, the couple stress-mediated load augment is assumed to be the same as that in Eq.
Based on these considerations, the couple stress-based unloading curve can be expressed as Incorporating Eq. (10) into Eq. (2), the couple stress-based contact stiffness can be expressed as where S 0 is the elastic contact stiffness, which is obtained by incorporating the leading term in Eq. (10) into Eq. (2), and expressed as According to Eq. (5), the size effects in modulus can only come from S and A c , because other parameters in Eq. (5) are prescribed constants. Equation (11) has introduced the size effects versus h max − h f into S. Similarly, the size effects could be introduced into A c by incorporating Eqs. (10) and (11) into Eq. (3). However, the respective size effects in P and S may be largely "faded" after they are incorporated into the term P/S, due to the quotient form of this term. Additionally, this incorporation would make the expression of the resultant model of A c much complex. Based on these considerations, for simplicity the size effects in A c are neglected here. Therefore, the size effects in modulus are modeled by directly incorporating Eq. (11) into Eq. (5), and expressed as where E 0 is the elastic indentation modulus, which is obtained by incorporating Eq. (12) into Eq. (5). Equation (13) describes the modulus ISEs as that the modulus varies with the completely elastic depth h max − h f . However, in an indentation test, the measured modulus and hardness are recorded as a function of h max . Therefore, in order to apply Eq. (13) to the experimentally obtained moduli, it should be also expressed as a function of h max . In the indentation tests of UHMWPE, PS and PMMA [3] and PAI [7], the fraction of elastic deformation work (the area encompassed by CDE in Fig. 1c) in total deformation work (the area encompassed by ABCE in Fig. 1c) is found to approximately remain constant as h max increases, as can be seen in Fig. 3. Therefore, for simplicity we assume that the elastic recovery where the constant η e essentially denotes the fraction of elastic depth, ranging from 0 to 1. Based on these approximations, we can express Eq. (13) in the form of where the relevant parameters are combined into a higher-order parameter l R by As the constant fraction η e is dimensionless, the higher-order parameter l R has a dimension of length. So far, the modulus ISEs have been characterized by the exponent m and the higher-order parameter l R . According to the understanding of effective behavior of indenter in the work of Oliver and Pharr, the theoretical range of m is considered to be around 1.5, at most ranging from 1 to 2. In a difference to m, to determine the theoretical range of l R is of difficulties, because l R contains five pieces of information, either from the material itself or from the indenter.
As discussed in Sect. 2, the hardness ISEs of polymers can be decomposed into elastic part and plastic part, and the elastic part may even be dominant. If the plastic hardness ISEs are temporarily neglected here, the elastic part can be similarly modeled by incorporating the couple stress-based elastic contact load, i.e., Eq. (10), into the hardness formula Eq. (6). Therefore, we can obtain where H 0 represents the macro hardness or constant hardness, which is obtained by substituting elastic contact load [the leading term in Eq. (10)] into Eq. (6) with h being h max . The parameters l R and m are exactly those in Eq. (14). However, it should be noted that if this hardness model is directly applied to the experimental data of indentation hardness of a polymer, the estimated values of l R and m may, respectively, deviate from those estimated by Eq. (14), because the plastic hardness ISEs in experimental data may be non-negligible and even significant. The same values as those estimated by Eq. (14) can be obtained only if the plastic hardness ISEs of this polymer can be extremely negligible. So far, the modulus ISEs and the elastic hardness ISEs are, respectively, modeled by Eqs. (14) and (16). It can be found from Eqs. (14) and (16) that a relation exists between the two augmenting effect terms in the two models, which can be expressed as In this equation, the term (H − H 0 )/H 0 and (E − E 0 )/E 0 , respectively, represent the increasing rate of hardness and modulus. Thus, Eq. (17) may exactly characterize the proportional relation (with the proportionality constant being 2m) suggested in Sect. 2. If the plastic hardness ISEs of a polymer can be extremely negligible, the values of l R and m estimated by Eq. (14) will be the same as those estimated by Eq. (16). Under Table 1 Values of E 0 , l R and m obtained by fitting Eq. (14) to the size-dependent experimental data in Fig. 7 Experimental data  Table 2 Values of H 0 , l R and m obtained by fitting Eq. (16) to the size-dependent experimental data in Fig. 8 Polymers this condition, if the experimentally obtained hardness and moduli are, respectively, processed in the form of Eq. (17) and plotted in a form versus 1/ h max , two power law curves with uniform exponent m − 0.5 will be obtained. In particular, these two curves become two straight lines if the values of m are uniformly estimated as 1.5, which requires the "effective indenter shape" in tests is equivalent to a paraboloid of revolution. Additionally, there is a proportional relationship of 3 times between the slopes of the two lines.

Model verification
In order to illustrate the validity of Eqs. (14) and (16) Tables 1 and 2, respectively. It can be seen that all m values in Table 1 are in the theoretical range of m. Specially, the m values of PS and PMMA (Briscoe et al. [3]) are significantly close to 1.5, which implies that the "effective indenter shape" used in these experiments can be well identified as a paraboloid of revolution. Additionally, in Table 1, the l R value of LDPE is larger than those of other polymers. Comparing Table 2 with Table 1, it can be found that for polymers whose modulus and hardness are measured, the values of l R and m in Table 2 are, respectively, greater than those in Table 1. Especially, in Table 2, the l R of PAI reaches about 500 nm and that of LDPE can even reach near 90,000 nm. The m value of LDPE in Table 2 can even reach 2.9 which significantly deviates from the theoretical range of m.

Discussion
Before proceeding to the discussion of results in Tables 1 and 2, let us discuss the influences of characteristic length l and elastoplastic properties of polymers on l R and m, because these features of different polymers are significantly discrepant. Regarding to the characteristic length l of polymers, Nikolov et al. [42] and Han [43] have studied its physical mechanisms and related it to the bending stiffness of polymeric chains. The phenomenological rotation gradient energy in couple stress elasticity is related by Nikolov et al. [42] to the Frank elasticity energy, which is always present as long as the chains possess finite bending stiffness and neighboring interactions. They worked out a relation: K 3 μl 2 linking the phenomenological characteristic length l to the effective Frank elasticity constant K . Han [43] subsequently illustrated in detail the dependence of characteristic length l on the molecular bending stiffness, although the exact l value of each polymer is not estimated. They found that polymers containing complex molecular structures in chains always have larger characteristic length l than those having flexible chains. For example, PAI, Epoxy, PC and PS have the longest Fig. 6 The fraction of elastic deformation work in total deformation work of UHMWPE, PS and PMMA [3] and PAI [7] l due to the stiffening of aromatic rings in the backbone or side group; PMMA and Nylon66 have shorter l, since PMMA has complex side groups and the Nylon66 has stiff areas of chains due to the presence of the O and N atoms [20], although they both lack aromatic rings; UHMWPE and PTFE have shortest l which can even be close to 0 nm due to their highly flexible and linear chains. LDPE was not considered in the work of Han. The l of LDPE is presently considered to be longer than those of UHMWPE and PTFE due to its easily branched chains, but shorter than those of other polymers due to the flexibility of chains.
The elasto-plastic properties of polymers mainly influence the parameter η e which denotes the fraction of elastic recovery in total elastoplastic deformation depth. In general, polymers with higher crosslink/entanglement density can provide higher elastic resistance than those containing flexible chains, because the flexible chains can easily adjust themselves to plastically dissipate the indentation work [43]. This kind of dependence is approximately consistent with the dependence of characteristic length l on the molecular bending stiffness. As can be seen in Fig. 6, the fraction of elastic work of PAI is larger than those of other polymers, and the corresponding fraction of UHMWPE is the lowest. According to the definition of η e in Sect. 3.1, it can be postulated that the parameter η e also approximately conforms the size sequence of molecular stiffness. Additionally, the results of finite element simulations for a variety of elastic to plastic materials [36] show that the exponent m obviously increases with decreasing η e . Furthermore, according to the analytical solutions of Sneddon [35], the increasing m implies that the effective shape of indenter is close to a cone, resulting in a smaller effective curvature radius R.
The l R and m values in Table 1 are consistent with the above considerations. According to the mathematic form of Eq. (15), the l R increases with increasing l, however decreases with increasing η e . This is the reason why the l R values in Table 1 do not strictly conform the size sequence of l. For example, the PS and Epoxy containing rigid aromatic rings in chains should theoretically have the largest l R than other polymers, but their η e value may also be the largest, thus the resulted l R value is not the largest one. For LDPE, if its η e value is considered to be close to that of UHMWPE, it also has the lowest η e value. According to the finite element simulations of Pharr and Bolshakov [36] mentioned above, its lowest η e causes the largest m value among these polymers. More importantly, the lowest η e and the largest m of LDPE, along with the correspondingly smaller effective curvature radius R, make the resulted l R value the largest among these polymers, although the LDPE has a shorter l than them. Conclusively, with these considerations, the results shown in Fig. 7 and Table 1 demonstrate that the modulus model Eq. (14) can be successfully applied to characterize the modulus ISEs of polymers.
The corresponding deviations between the values of l R and m in Table 2 and those in Table 1 may be attributed to the discrepancy that the experimentally obtained hardness ISEs contain both elastic part and plastic part, while Eq. (16) considers only the elastic part. As discussed above, the plastic dissipation of LDPE in indentation is significant due to its relatively flexible chains, thus the plastic hardness ISEs of LDPE may even be dominant. This is exactly the reason why the values of l R and m of LDPE significantly deviate from those predicted by Eq. (14), and from the theoretical range. For other polymers in Table 1, the fraction of plastic deformation is smaller than that of LDPE, thus the plastic hardness ISEs are slighter than those of LDPE. Therefore, the corresponding deviations are not as significant as that of LDPE. It can even be postulated that for all polymers shown in Fig. 8, if the elastic hardness ISEs can be separated from the experimental data and fitted by Eq. (16), the same values of l R and m as those predicted by Eq. (14) will be obtained. That is to say, despite these deviations, the hardness model Eq. (16) can still be applied to characterize the elastic hardness ISEs. In particular, the PAI has the lowest fraction of plastic deformation among all polymers. Therefore, the plastic hardness ISEs of PAI may be trivial. If its modulus data are obtained by experiment and fitted by Eq. (14), the same values of l R and m as their counterparts [estimated by Eq. (16)] in Table 2 will be obtained.
The influences of the plastic hardness ISEs can also be explicitly reflected by relation Eq. (17). For each polymer depicted in Fig. 7, the (E − E 0 )/E 0 values where E 0 comes from Table 1 Table 2 are also plotted. According to the discussion of Eq. (17) in Sect. 3.1, for each polymer whose m values in both Tables 1 and 2 are uniformly close to 1.5, the data points of either (H − H 0 )/H 0 or (E − E 0 )/E 0 will line up more straightly. Additionally, a triple relation will exist between the slopes of the two linear trends. For a more explicit comparison, a straight line passing through the point (0, 0) is, respectively, fitted to these data points. As can be seen in Fig. 9, this kind of tendency is approximately reflected by the data points of PS, PMMA, and Nylon66. As deviant E 0 and H 0 values are predicted for Epoxy (Alisafaei et al. [8]) in Figs. 7 and 8, its data points in Fig. 9 slightly deviate from the linear tendency. However, it can be found that the data points of LDPE significantly deviate the linear trend. The corresponding slopes of these linearly fittings are shown in Table 3. As can be seen therein, for PS, PMMA and Nylon66, the triple relation is obvious. However, for LDPE the relation in slopes significantly deviates from the triple relation, which can be exactly attributed to the significant fraction of plastic hardness ISEs. Finally, it is noted that the linear fits in Fig. 9 are adopted to illustrate not the validities of Eqs. (14) All results present above indicate that for polymers shown in Fig. 8 other than LDPE, the hardness ISEs are significantly elastic, while the hardness ISEs of LDPE are significantly plastic. Additionally, these results yield the presentiment that if the plastic hardness ISEs are to be characterized by a model, this model may be in a similar mathematical form to Eq. (16), but with an exponent larger than m − 0.5. Firstly, as can be seen in Fig. 8, although the experimentally obtained hardness of all polymers have already contained plastic ISEs, the mathematical form of Eq. (16) can still give an exact graphical description of the whole hardness ISEs. This implies the plastic hardness ISEs can also be described by a formula in the similar form to Eq. (16), i.e.,  Fig. 7, respectively, plotted with respect to 1/ h max . The E 0 value and H 0 value of each polymer, respectively, come from Tables 1 and 2. A straight line passing through (0, 0) is, respectively, fitted to these data points the augmenting effect term is proportional to the power of 1/h max . Secondly, for polymers whose moduli are also measured, the m value estimated by Eq. (16) is, respectively, larger than those estimated by the modulus model Eq. (14). This implies that the corresponding exponent in plastic hardness ISEs will be larger than m − 0.5. According to the fitted results of LDPE, this exponent may be larger than 2.4. Thirdly, this inference is compatible with a categorization of the hardness ISEs suggested by Han [43], which says that for a polymer the depth range presenting elastic hardness ISEs is always larger than that presenting the ISEs mediated by other factors. It can be easily seen from Eq. (15) that as the exponent m − 0.5 becomes larger, the augmenting term will become smaller at the depths beyond m−0.5 √ l R , and larger at the depths within m−0.5 √ l R . That is to say, at the depths beyond m−0.5 √ l R the elastic hardness ISEs will be more obvious than the plastic hardness ISEs, and consequently more easily observed by experiments.

Models of the plastic ISEs of polymers
This section attempts to model the plastic hardness ISEs of polymers. The indentation test begins with the elastic deformations of the polymer solid. As the indenter gradually penetrates into the sample, the initial reversible elastic resistances of the region beneath the sharp indenter eventually give way to plastic behaviors. Then, the elastic and plastic resistance against the indenter exist meanwhile in the sample in which a large elastic field encompasses the plastic deformation zone. The deformation constitutive of all isotropic solids can be simply formulated. As illustrated in Sect. 3, the size effects in the elastic resistance can also be characterizing by introducing size characteristics into the constitutive formula (see the couple stress theory for more details). The plastic resistance can also yield size effects. A common example is the size effects of metals. As discussed above, especially for LDPE the plastic dissipation in indentation is significant due to its relatively flexible chains, thus the plastic hardness ISEs of LDPE may even be dominant. However, the mechanisms of plastic deformation or the carriers of plastic flow of different materials are completely different. Therefore, the characterization of plastic size effects of polymers is highly dependent on the development of plastic deformation theory of polymers.
Maybe we just need to pay attention to the plastic deformation of the glassy polymers The non-existence of hardness ISEs of UHMWPE and PTFE displayed in Figs. 4 and 5 yields the presentiment that the plastic hardness ISEs of semi-crystalline polymers may be derived from only their glassy components, but not their crystal components. In semi-crystalline polymers, the deformations begin with the elastic and subsequently plastic deformations of glassy components. In a local region, after glassy components are exhausted, the deformations are transferred to the crystal components by various crystallographic processes up to large plastic strains [44]. UHMWPE and PTFE have the highest degree of crystallinity due to their highly linear chains, thus their plastic deformations beneath the indenter mainly occur in crystal components. However, as can be seen from Fig. 8, both UHMWPE and PTFE exhibit no hardness ISEs. This indicates that the plastic deformation of their crystal components mediates no size effect. It can be postulated that for semi-crystalline Nylon66 and LDPE, the plastic deformation of crystal components may also mediate no size effect, and thus their plastic hardness ISEs are only related to glassy components.

Model development
On the basis of the above discussions, we can accordingly focus on the plastic deformation of glassy polymers and expect to incorporate size characteristic into the formula of plasticity. The understanding of the basic mechanism of plastic deformation of glassy polymers has always been a subject of intense interest. The early famous studies in this respect include the Eyring model [45], the Roberson model [46] and the Argon's kink pair model [47]. Argon's model successfully characterizes not only the pressure, strain rate, and temperature sensitivity of the plasticity but also the softening and secondary hardening of glassy polymers. These studies interpreted the micro mechanism of plasticity as the motions of molecular chains like rearrangements, straightening and localized twist.
Recently, Voyiadjis and Samadi-Dooki [21] developed a shear transformations (STs) plasticity formula for glassy polymers. The notion of STs was originally proposed by Spaepen [22] and Argon [23] for glassy metals, which represents the atomic glide, slip, or shear rotation inside the lumped localized deformation sites. Oleinik et al. [48] have revealed that within the plastically deformed samples there exists significant amount of stored energy, which can never be carried through chain rearrangements. This is one of the motivations for Voyiadjis and Samadi-Dooki to develop their plasticity formula. Their constitutive formula shows its success in not only capturing the primary softening behavior, but also justifying the secondary hardening of PMMA. Additionally, the nucleation of STs is associated with the effects of the strain rate, temperature, and thermal history of the sample on its post-yield behavior.
According to the shear transformation plasticity theory of Voyiadjis and Samadi-Dooki [21], at unstressed state, a large amount of free volume sites with different sizes are distributed in the microstructures due to the amorphous nature of glassy polymers. Among these free volume sites, these with excessive free volume are more likely to foster atomic STs under the enforcing of shear stress, as shown in Fig. 10. These sites with excessive free volume are only a fraction of all free volume sites. Similar to the notion of potential jump sites in glassy metals, these sites may be called potential transformation sites. It is noted that the STs are not any existing defects but dynamic events only instantaneously nucleating under enforcing of shear stress. Those potential transformation sites where the ST events are occurring are called shear transformation zones (STZs) (see Fig. 10). The STZs have a characteristic average volume and a characteristic shear strain γ T .
The kinetics relation between the global shear strain rateγ P and global shear yield stress τ is expressed aṡ where k B is Boltzmann constant, T is absolute temperature,γ P 0 is a pre-exponential which depends on the shear strain γ T and the frequency of atomic vibration (~Debye frequency) in the range of 10 10 s −1 , F is the average nucleating energy of STs, and c f is a factor depending on the volume fraction of potential transformation sites.
The average nucleating energy of STs with shear strain of γ T and average volume of , is considered equal to the deformation energy in the inclusion model of Eshelby [49] and expressed as where μ is the shear modulus, and f (ν, β) accounts for shear and dilatational components of the transformation strain tensor. For polymers with the Poisson's ratio of 0.35-0.4, f (ν, β) has a value close to 0.5 [10,21]. Voyiadjis et al. [10] subsequently suggested a statistical interpretation to the ISEs of glassy polymers based on the shear transformation theory. They suggested that the nucleation of STs in the deformation zone beneath the indenter is a probabilistic phenomenon. On one hand, the size of a unit potential transformation site could be manifested as a spherical region with a diameter of about 10 nm; on the other hand, the potential transformation sites are only a fraction of all free volume sites and discretely distributed in the material. Therefore, at very Fig. 11 The diagram for the distribution of shear transformation zones along the indentation depth small indentation depths, the deformation zone is not big enough to "see" a potential transformation site, and thus atoms in this zone can hardly adjust themselves (shear transformation) to dissipate the indentation work, as schematically in Fig. 11a, consequently the corresponding deformation zone is highly stressed. As the indenter continues to move down, the region around the highly stressed zone starts to "see" potential transformation sites and form ST events (Fig. 11b). At very deep indentation depths, the ST events could massively nucleate and control the kinetics of deformation (Fig. 11c).
The above statistical interpretation of size effects provides the novel and intuitive insight to understand the plastic ISEs of polymers. In their size effects model, the total yield stress of the region beneath the indenter is formulated as τ total χτ ST + (1 − χ )τ local , where τ ST is the STs mediated stress, τ local is the non-STs mediated local high stress and χ denotes the total probability of finding a fertile site that can undergo shear transformation. The size characteristics are incorporated into the probability of the nucleation of STs. Thus, their proposed hardness model can capture size effects. Following the development process of plasticity formula of Voyiadjis and Samadi-Dooki [21], perhaps we can directly incorporate the size characteristics into Eq. (18), since the volume fraction factor c f has already incorporated the information of deformed volume due to their extension on this factor. The factor c f in Eq. (18) derives from the plasticity theory of glassy metals proposed by Spaepen [22] and Argon [23], which denotes the total probability that an atom is on a potential jump site, or the fraction of potential jump sites. Voyiadjis and Samadi-Dooki extended this factor to denote the volume fraction of potential transformation sites of glassy polymers and expressed it as where λ is a geometrical constant, N is the total number of the STZs in glassy polymers. ν * is the minimum excessive free volume which can act as a potential transformation site, f V p represents the total free volume with f being the fraction of the free volume and V p being the total plastic volume. Equation (20) implies that the nucleation of STs depends on the size characteristics of material required to plastically deform whereas the STs themselves are the carriers of plastic deformations. In the ISEs model of Voyiadjis et al. [10], in order to express the total probability of finding a fertile site that can undergo shear transformation, they assumed that in indentation tests the plastic volume V p and minimum excessive free volume ν * are, respectively, proportional to the volume of the hemisphere under the indenter and the average volume of STZs. The same assumption is employed here. Thus, after other parameters in Eq. (20) are incorporated into a fitting parameter k, it can be rewritten as It should be noted that the c f , which derives from Eq. (18) and depends on the volume fraction of potential transformation sites, is different from the probability factor χ in the ISEs model of Voyiadjis et al. [10], even though the same expression and assumptions are employed.
Incorporating Eq. (21) into Eq. (18) and considering that the parameter γ T τ is far larger than k B T, Eq. (18) can be converted to It can be seen that the indentation size has been successfully introduced into the constitutive formula. When the indentation is implemented at very small depths, the volume fraction factor in Eq. (21) would rapidly decrease due to the exponent of h max being 3, which denotes that the deformation zone is too trivial to foster a ST event. The shear stress in Eq. (22) would consequently be given a rapid increase up to a state which may accordingly represent the non-STs mediated local high stress.
With the relation H 6.6τ [50] a hardness model can be expressed as where the last term denotes the STs plasticity mediated hardness ISEs. When the indentation is implemented at far deeper depths, the STs are massively nucleated and the hardness augment tends to be 0. The hardness model Eq. (23) can be rewritten as where the parameter V is expressed as and the macro hardness H 0 is expressed as It can be found that Eq. (24) has a similar form to Eq. (16) which characterizes the elastic hardness ISEs. More importantly, the exponent of 1/ h max in Eq. (24) is 3, being larger than the exponent m − 0.5 in Eq. (16). These two features are well consistent with the inference in the end of Sect. 3.3. As discussed therein, if the plastic hardness ISEs are to be characterized, the corresponding hardness model should have a similar form to Eq. (16), but with an exponent larger than m −0.5. Additionally, a dimensional analysis to Eq. (25) suggests that the parameter V has a dimension of volume. This parameter can be seen as a higher-order volume parameter similar to the higher-order length parameter l R in Eq. (16).

Model verification and discussion
In principle, the validity of model Eq. (24) should be illustrated by hardness data containing only plastic ISEs. However, to separate the plastic ISEs from hardness data of polymers in Fig. 5 is impossible. A compromise is applying it to the hardness data of LDPE [5] in which the plastic ISEs are dominant. As shown in Fig. 12a, a good agreement between the model and the experimental data is found. The macro hardness H 0 and volume parameter V are, respectively, estimated to be about 0.025 GPa and 2 × 10 6 nm 3 . Then the characteristic volume of STZ is estimated to be about 2718 nm 3 (for choice γ T 0.04) and parameter k is estimated to be about 5 × 10 4 , respectively, from Eqs. (26) and (25). Taking the obtained and k into Eq. (21), the volume fraction factor with respect to h max is depicted in Fig. 12b. It can be seen that the curve is compatible with the hardness curve. When the indentation is implemented at small depths, the deformation zone is too small to foster a ST event. As the indenter continues to move down, more and more potential transformation sites are found and STs are massively nucleated. Accordingly, the hardness data tends to be stable.
The presently obtained of LDPE is close to the STZ volumes of several other polymers given by Mott et al. [51] through molecular simulations, which range from hundreds to thousands of cubic nanometers. The presently obtained k also has the same order of magnitude as that of PMMA estimated by Voyiadjis et al. [10]. Note that the determination of is significantly sensitive to the chosen value of γ T , for example, the will be 1183 nm 3 for choice γ T 0.05. The γ T of PMMA is believed to be in the range of 0.03 to 0.05 by the interpolation of experimental data (Malekmotiei et al. [52]; Voyiadjis et al. [10]), and thus the middle of this range is adopted in the present estimation. The temperature T is set at 295 K for room temperature experiments, as Tavares et al. [5] did not specify that the experiments were carried out under heating or cooling conditions. The shear strain rateγ P has been found to alter the macro hardness [10], as illustrated in Eq. (26). However, its contribution is not able to justify the observed indentation size effects in hardness [53], although it emerges into the higher-order volume parameter V . In above estimation of , the value ofγ P is obtained froṁ γ P √ 3Cε i where C is a constant value equal to 0.09 andε i the indentation strain rate [52,54]. Tavares et al. [5] did not provide the indentation strain rateε i , but stated that a high rate was adopted to minimize the influence of viscoelastic properties in the measurements. The above estimation assumesε i to be 1 s −1 . Additionally, the shear modulus μ is obtained from 2μ E/(1 + ν) where ν 0.4 [55] and E is the macro modulus in Table 1 as Eq. (26) denotes the macro hardness.
It is finally noted that the presently proposed ISEs model is a slight extension of the STs plasticity theory proposed by Voyiadjis and Samadi-Dooki. Specifically, some assumptions also considered in the hardness model proposed by Voyiadjis et al. are incorporated into the STs plasticity formula to make it capable for capturing the size dependent deformations in indentations. Therefore, the proposed model can essentially be treated as an application tentatively demonstrating the validity of the STs plasticity theory of glassy polymers, considering the agreement displayed in Fig. 12 and the obtained acceptable fitting results. However, the shear transformations of glassy polymers require further experimental and simulated evidences. The presently proposed model also needs further experimental verification implemented on various glassy polymers, which would be the focus of our future work.

Conclusions
The size effects in indentation modulus (modulus ISEs) and indentation hardness (hardness ISEs) of polymers are theoretically modeled. By recalling the Oliver-Pharr indentation method, and comparatively analyzing the experimental from several literatures, the finding of the literature that for polymers the ISEs are of significant elastic nature is confirmed. Specifically, the hardness ISEs of polymers are found to be significantly associated with elastic deformation. The modulus ISEs and the elastic hardness ISEs are modeled, by introducing a couple stress elasticity-based unloading model proposed in our previous work into Oliver-Pharr indentation approach. The resultant modulus model and hardness model, along with the corresponding results of application to experimental data, show that these elastic size effects and their experimental showing are mainly determined by the molecular structures. Specifically, polymers with complex molecules would exhibit significant modulus ISEs. Besides, the hardness ISEs of this kind of polymers are proved to be significantly elastic.
The couple stress elasticity-based hardness model shows a shortage in characterizing the hardness ISEs of LDPE. We postulate that the hardness ISEs of LDPE are significantly plastic, due to its highly flexible molecular structures. As both the UHMWPE and PTFE with high crystallinity exhibit no indentation size effects, we postulate that the plastic hardness ISEs of polymers are only determined by their glassy components and may be described by a model in a similar form to that of the couple stress elasticity-based hardness model. Based on these considerations, a recently proposed shear transformation plasticity theory of the literature is extended to characterize the size dependent deformation in indentation. The accordingly proposed hardness model is consistent with previously mentioned presentiment and successfully applied to the related experimental data. The plastic ISEs of glassy polymers are manifested as that at small indentation depths the sample can hardly nucleate the shear transformations to dissipate the enforcing of the indenter.