Frédéric Hardy
email: frederic.hardy@kit.edu

Robert Eder

Martin Jackson

Dai Aoki

Carley Paulsen

Thomas Wolf

Philipp Burger

Anna Böhmer

Peter Schweiss

Peter Adelmann

Robert A Fisher

Christoph Meingast

Multiband Superconductivity in KFe 2 As 2 : Evidence for one Isotropic and several Liliputian Energy Gaps

Keywords: KFe 2 As 2, iron pnictide, multiband superconductivity, heat capacity, magnetization

We report a detailed low-temperature thermodynamic investigation (heat capacity and magnetization) of the superconducting state of KFe 2 As 2 for H || c axis. Our measurements reveal that the properties of KFe 2 As 2 are dominated by a relatively large nodeless energy gap (∆ 0 = 1.9 k B T c ) which excludes d x 2 -y 2 symmetry. We prove the existence of several additional extremely small gaps (∆ 0 < 1.0 k B T c ) that have a profound impact on the low-temperature and low-field behavior, similar to MgB 2 , CeCoIn 5 and PrOs 4 Sb 12 . The zero-field heat capacity is analyzed in a realistic self-consistent 4-band BCS model which qualitatively reproduces the recent laser ARPES results of Okazaki et al. (Science 337 (2012) 1314). Our results show that extremely low-temperature measurements, i.e. T < 0.1 K,will be required in order to resolve the question of the existence of line nodes in this compound.

Introduction

The pairing mechanism in iron-pnictide superconductors is still a subject of intense debate. Similarly to heavy fermions, cuprates and ruthenates, the proximity of these materials to a magnetic instability naturally suggests that spin fluctuations can mediate the formation of Cooper pairs, although other scenarios involving orbital fluctuations are possible. [1][2][3][4][5] In this context, the symmetry of the superconducting-state order parameter can have either an s± or a d-wave symmetry. Unfortunately, these states are almost degenerate and the realization of one of these two states is material-specific, depending on the number and position of Fermi-surface sheets in the Brillouin zone and their mutual interactions. In this context, the interpretation of experimental data is very complicated since the existence (absence) of nodal behavior does not permit the ruling out of s-wave (d-wave) symmetry. The Ba 1-x K x Fe 2 As 2 series is a prominent example. Indeed, at the optimal concentration (x ≈ 0.4), heat-capacity 6 and ARPES 7,8 measurements give strong evidence of an s-wave state while in the strongly correlated end-member KFe 2 As 2 , that has only hole pockets (see Fig. 1(a)). The situation remains highly controversial. Recently, thermal-conductivity measurements of Reid et al. 9 were found to extrapolate at T → 0 to a finite residual term κ(0)/T, independent of sample purity. This was interpreted as a signature of universal heat transport, a property of superconductors with symmetry-imposed line nodes. This hypothetical change from s-to d-wave symmetry as a function of doping is allowed theoretically via an intermediate s + id state that breaks time-reversal symmetry. [10][11][12][13] On the other hand, laser ARPES 14 have revealed the existence of accidental line nodes on only one of the zone-centered pockets in KFe 2 As 2 , which is only compatible with a nodal s± state.

However, none of these methods are bulk probes of the superconducting state. In this Article, we report a detailed lowtemperature thermodynamic investigation (heat capacity and magnetization) of the superconducting state of KFe 2 As 2 . We show quantitatively that the properties of KFe 2 As 2 , including the upper critical field (H c2 ), are dominated by a relatively large nodeless energy gap of amplitude 1.9 k B T c which excludes de facto d x 2 -y 2 symmetry (see Fig. 1(b)). We prove the existence of several additional extremely small gaps (∆ 0 < 1.0 k B T c ) and show that they have a profound impact on the low-temperature and low-field behavior, as previously shown experimentally [17][18][19][20][21] and theoretically [22][23][24][25][26] for MgB 2 , CeCoIn 5 and PrOs 4 Sb 12 . The zero-field heat capacity is analyzed in a realistic self-consistent 4-band BCS model which qualitatively reproduces the recent laser ARPES results of Okazaki et al. 14 We also find that extremely low-temperature measurements, i.e. T < 0.1 K, are required to observe the signature of possible line nodes in KFe 2 As 2 . In accord with recent angleresolved heat-capacity experiments, [START_REF] Kittaka | SCES conference[END_REF] our results are compatible with either a d xy or a nodal s± state.

Experimental details

Single crystals of KFe 2 As 2 were grown in alumina crucibles using a self-flux method with a molar ratio K:Fe:As=0.3:0.1:0.6. The crucibles were put and sealed into an iron cylinder filled with argon gas. After heating up to 700 • C and then to 980 • C, the furnace was cooled down slowly at a rate of about 0.5 • C/h. The composition of the samples was checked by energy-dispersive x-ray analysis and four-circle diffractometry. The specific heat was measured with a commercial Quantum Design Physical Property Measurement System (PPMS) for T > 0.4 K and with a homemade calorimeter for T < 0.4 K. For T > 2 K, we used a vibrating sample magnetometer to measure the magnetization. At lower temperature, magnetization measurements were per- terference device (SQUID magnetometer) equipped with a miniature dilution refrigerator developed at the Institut Néel-CNRS Grenoble. [START_REF] Burger | [END_REF] 3. Zero-field electron specific heat, C e (T, 0) Figure 2(a) shows the low-temperature heat capacity of KFe 2 As 2 . We find a large Sommerfeld coefficient γ n = 103 mJ mol -1 K -2 and T c = 3.4 K, in agreement with our previous studies. 29 Below 0.2K, the high-temperature tail of a Schottky anomaly, probably due to paramagnetic impurities, is observed (see inset Figure 2(a)). The electronic contribution C e , shown in 2(b), is obtained by subtracting, from the measured data, a Debye term (inferred from the 5 T data), and the Schottky contribution.

The overall curve bears a strong similarity with that of MgB 2 . 30 In particular, we observe that: (i) the jump at T c , ∆C/γ n T c ≈ 0.54, is substantially smaller than the BCS value (∆C/C n = 1.43) for a weakly coupled single band s-wave superconductor and (ii) there is a steep quasi-linear decrease of C e /T with decreasing temperature for T/T c ≤ 0.1. A similar linear dependence of the penetration depth was reported by Hashimoto et al. 31 and was interpreted as evidence of line nodes. In this Article, we argue that this steep feature is instead related to the existence of small energy gaps (∆ S /k B T c ≈ T/T c ≈ 0.2), as inferred from small-angle neutron scattering (SANS) experiments. 32 Assuming that all 3 sheets around the Γ point exhibit this tiny gap and using the expression of the heat-capacity jump of a two-band s-wave superconductor in the weak coupling limit, 33,34 ∆C

k B T c = 1.43 • N S ∆ 2 S + N L ∆ 2 L 2 (N S + N L ) N S ∆ 2 S + N L ∆ 2 L , (1)
(where the subscripts S and L refer to the small and large gaps, respectively). We estimate ∆ L /k B T c ≈ 1.8 on the band using the individual density of states inferred from dHvA and ARPES measurements 15,16 (see Table I). Interestingly, we obtain from this simple estimation a remarkably large gap anisotropy ∆ L /∆ S ≈ 9, in comparison with MgB 2 where it is about 4. 30 Thus, KFe 2 As 2 represents a somewhat extreme case of multiband superconductivity similar to the heavy-

i = πN A k 2 B a 2
3 2 m * i (with a = 3.84 Å). 15,16 The last column contains the densities of states used in the 4-band BCS model. m e is the bare electron mass. 4. Mixed-state reversible magnetization, M rev (H)

dHvA & ARPES C(T) m

* i γ i γ i /γ n N i (0)/N(0) (m e ) (mJ mol -1 K -2 ) α 6.
µ 0 H ( T ) M ( 1 0 3 A / m ) H | | c ( a ) ( c ) ( d ) ( b ) 0 . 6 K 1 . 5 K 2 . 0 K 2 . 6 K 4 . 0 K M r e v ( 1 0 3 A / m ) µ 0 H ( T ) 0 . 6 K 1 . 5 K 2 . 0 K 2 . 6 K 3 . 0 K M r e v -M [ 4 K ] ( 1 0 3 A / m ) µ 0 H ( T )
Figure 3(a) shows magnetization curves of KFe 2 As 2 for H || c down to 0.6 K. In the normal state (T = 4 K), a sizeable paramagnetic signal is observed with a susceptibility of about 3.3 × 10 -4 in agreement with our previous report. 29 At lower temperatures, the magnetization curves are reversible over a wide field interval (e.g. H c2 /2.5 < H < H c2 = 1.4 T at 0.6 K), indicating a small concentration of pinning centers as confirmed by the observation of a well defined hexagonal vortex lattice by SANS. 32 Together with the observation of quantum oscillations and the absence of significant residual density of states in the limit T → 0 (see Fig. 2(b)) this indicates that our KFe 2 As 2 single crystals are weakly disordered and are in the clean limit. Note that this is at odds with Co-doped BaFe 2 As 2 35, 36 As a result, accurate reversible magnetization curves can be obtained for our sample by averaging the increasing and decreasing branches of the magnetization loop, as illustrated in Fig. 3(b).

In single-band type II superconductors, M rev (H) is entirely defined by H/H c2 and the Ginzburg-Landau parameter κ = λ/ξ, 37 with

M rev = H -H c2 (2κ 2 -1)β A + 1 , (2)
at high field (Abrikosov regime, with β A the Abrikosov coefficient). In the intermediate field range (London regime), [START_REF] De Gennes | Superconductivity of metals and alloys[END_REF][START_REF] Kogan | [END_REF] the reversible magnetization is linear in the logarithm of the applied field with

µ 0 M rev = - φ 0 8πλ ln b H c2 H , (3)
where λ = λ ab is the in-plane penetration depth for H || c and b a constant. 40 Thus, in close analogy to the case of MgB 2 , 41 Fig. 3(c) shows that the linear evolution of M rev (H) expected near H c2 is not observed and the London dependence dominates up to H c2 ,which therefore allows the determination of λ ab (T) using Eq.( 3). The derived superfluid density, defined as ρ ab (T) = [λ ab (0.6 K)/λ ab (T )] 2 , is shown in Fig. 3(d) together with the calculation for an s-wave gap of amplitude ∆ 0 = 1.9 k B T c , which accurately reproduces the data. As a consequence, our analysis firmly establishes the existence of a relatively large nodeless gap in KFe 2 As 2 , in agreement with our rough above-mentioned heat-capacity analysis. Contrary to direct penetration-depth measurements, 31 our estimate of ρ ab (T) was inferred from high-field data where the large vortex cores related to the smaller gaps have already overlapped, as observed in MgB 2 , 42 and discussed hereafter. This explains why this indirect derivation of ρ ab (T) is only sensitive to the larger gap as found for MgB 2 . 43,44 5. Mixed-state specific heat, γ (H) Evidence for the existence of tiny energy gaps can be found using heat-capacity measurements in the mixed state as previously shown for pure, Al-and C-doped MgB 2 . 17,19,45,46 Figure 4(a) shows the field dependence of the electron heat capacity γ(H) at 0.12 K (i.e. T/T c = 0.035) for H parallel to c. Similar to MgB 2 , we find that γ(H) is very non-linear with applied magnetic field. In very low fields, γ(H) increases abruptly and reaches γ(H)/γ n ≈ 0.6 at only H/H c2 ≈ 0.1. In larger fields, γ(H) closely follows the behavior expected for an individual band (magenta line), 47 indicating that the interband couplings between the sheet with the largest gap and the other bands is rather small. 23 Using the densities of states derived from dHvA and ARPES (see Table I), we can unambiguously ascribe the largest gap to the band and subtract its contribution from γ(H) to obtain the mixed-state heat capacity of the remaining α, β and ζ bands (green curve in Figs. 4(a) and 4(b)). We find that these bands have almost recovered their normal-state value in a 'crossover' field H S c2 ≈ 0.1 × H c2 which is defined here as,

H S c2 ≈ H c2 ξ S ξ L 2 , (4)
and which would correspond to the upper critical field of the small gaps in the absence of interband couplings. 23 Thus, the disappearance of the small gaps associated with the α, β and ζ sheets in an applied field is more rapid than that of the band which shows a conventional individual dependence. Following Klein et al., 41 we assume that all the excitations are localized in the vortex cores (i.e. we neglect the small-gap Doppler shift 48 ) and that the system can be described by only one field dependent quantity, ξ c (H), which is a measure of the vortexcore size. In this context, γ(H) ∝ γ n • (ξ c (H)/d) 2 (where d ∝ 1/ √ H is the intervortex distance), and we obtain directly ξ c (H) as shown in Fig. 4(c), with ξ c (H=H c2 )= Φ 0 /2πµ 0 H c2 . We find that the vortex-core size smoothly decreases from 50 to 15 nm in high fields, explaining the smooth evolution of the contributions of the α, β and ζ bands to γ(H) near H S c2 . Thus, the small gaps on these sheets remain finite due to nonzero interband coupling even for H >> H S c2 where their vortex cores overlap. In the opposite limit, i.e. where the Doppler shift (Volovik effect) dominates, we note that a similar dependence of γ(H) for the small gaps (see Fig. As 2 (blue symbols) at T = 0.12 K for H || c. The magenta line is the theoretical curve of the mixed-state heat capacity of an s-wave superconductor 47 normalized by the density of states of the band (see Table I). The green curve is the resulting contribution of the 3 small energy gaps α, β and ζ obtained by subtracting the heat capacity of the band from the data. (b) Close-up of the low-field region. (c) Field dependent vortex core size derived from the mixed-state heat capacity.

Comparison with thermal-conductivity measurements, κ (T, H)

In light of our results, we comment here on the interpretation of recent heat-transport experiments in KFe 2 As 2 . In Refs. 9 and, 49 thermal-conductivity measurements κ(T)/T, performed for T > 0.1 K, were found to extrapolate at T → 0 to a finite residual term κ(0)/T independent of sample purity. This was interpreted as a signature of universal heat transport, a property of superconductors with symmetry-imposed line nodes such as d-wave states. Experimentally, these measurements were not strictly realized in zero magnetic field because it was necessary to apply a small field of 0.05 T (i.e. H/H c2 ≈ 0.03) to suppress superconductivity of the soldered contacts. However, as shown in Fig. 4(b), this small field is large enough to produce an enhancement of the density of states, reaching 40% of the normal-state value at 0.12 K, which inexorably leads to a finite value of κ 0 /T. In addition, our specific-heat measurements show that a significant increase of κ(H)/T is also to be expected for H/H c2 < 0.1 due to these small gaps. This feature was not observed in Refs., 949 or in the more recent data of Watanabe et al. 50 while it clearly appears in many other multiband superconductors including MgB 2 , CeCoIn 5 and PrOs 4 Sb 12 . 18,20,21 Thus, the origin of the finite κ 0 /T cannot be attributed to d-wave superconductivity in KFe 2 As 2 in these experimental conditions. We note that the use of superconducting solder was already pointed out to produce spurious results in Ref.. 21 On the other hand, our observation of a relatively large isotropic gap does not rule out definitively d-wave superconductivity in KFe 2 As 2 . Actually, only the d x 2 -y 2 order parameter, with nodes located on the diagonals of the Brillouin zone, is excluded while d xy symmetry remains possible if the large gap effectively occurs on the sheet. This conclusion is corroborated by recent angleresolved heat-capacity measurements of Kittaka et al. [START_REF] Kittaka | SCES conference[END_REF]

Four-band BCS analysis of C e (T, 0)

Recently, laser ARPES measurements, 14 performed at 1.5 K, revealed angle-dependent gaps of the 3 bands around the Γ point with accidental line nodes only on the ζ sheet. These results are compatible with a nodal s-state and exclude a possible change of symmetry of the superconducting order parameter as a function of K doping. To check whether these results can be confirmed by bulk measurements, we analyze our zero-field heat capacity, taking into account the observed Fermi surface. In the absence of a sizeable residual density of states for T → 0, the modest specific-heat jump clearly shows that KFe 2 As 2 is close to the weak coupling limit. Therefore, we can model the temperature dependence of C e in a pure 4band BCS model. Assuming 2D Fermi-surface pockets, we obtain the following system of gap equations:

∆ i (φ, T ) = -4 j=1 N j (0) 2π 2π 0 dφ • c 0 d V i j (φ,φ )∆ j (φ ,T ) √ 2 +|∆ j (φ ,T )| 2 tanh β 2 2 + |∆ j (φ , T )| 2 , (5)
where N i (0) is the density of states of the i-th band with i ∈ {α, ζ, β, }, V i j (φ, φ ) are the intraband (i = j) and interband (i j) pairing potentials, β=1/k B T, and φ and φ the azimuthal angles on the sheets i and j, respectively. In the swave channel, 51,52 we write:

V i j (φ, φ ) = V (0) i j + V (1) i j • (cos 4φ + cos 4φ ). (6)
Such interactions lead to anisotropic gaps of the form:

∆ i (φ, T ) = ∆ (0) i (T ) + ∆ (1) i (T ) • cos(4φ), (7)
which are calculated self-consistently from Eqs. ( 5) and used to compute the superconducting-state heat capacity. We constrain all the N i (0) to match as closely as possible the values inferred from dHvA and ARPES measurements. In this form, the model is parametrized with 20 interaction constants and this number is reduced to 5 by assuming that:

V (0) i j = δ i j • V 1 + (1 -δ i j ) • V 2 , (8)
V

i j = 0, (1)
V β (φ, φ ) = V (1) = V (1) α = 0, (9)
with i, j ∈ {α, ζ, β}. Here, Eqs. ( 8) and ( 9) impose that the intra-and interband interactions of the zone-centered bands are angle-independent and equal to V 1 and V 2 , respectively because these bands have a quasi-2D morphology and are centered around the same point. 52 On the other hand, inelastic neutron scattering experiments have revealed the persistence of resonant spin excitations in heavily overdoped Ba 1-x K 1-x Fe 2 As 2 (x ≈ 0.9) 53 and incommensurate spin fluctuations in KFe 2 As 2 that approximately connect the Γ and X bands. 54 These observations convincingly indicate that the Γ-X interband interactions remain significant in KFe 2 As 2 , even in the absence of electron pockets. This is particularly true for the α and ζ pockets which are strongly involved in nesting in Ba 0.6 K 0.4 Fe 2 As 2 . In KFe 2 As 2 , these bands share a dominant xy/yz orbital character with the band while the β pocket exhibits mainly x 2y 2 component. This implies that the β sheet plays no decisive role in pairing, as illustrated by the small gap observed on this pocket in both compounds. 8,14 Furthermore, the ζ band shows an additional finite z 2 component which is absent from the α and sheets. Consequently, the z 2 component in the ζ pocket has no counterpart for the sign change of the superconducting gap in the sheet and this can lead to a sign change of the gap in the parts of the ζ band where the z 2 contribution dominates. 14 Thus, the only prominent angle-dependent interaction is V (1) ζ . All these assumptions are summarized in Eq. (10) and Fig. 5.

The heat capacity, as well as the angular and the temperature dependence of the gaps calculated with this model are shown in Fig. 6, using the density of states and the remaining tively. These parameters are all given in units of V 1 because it is not their absolute values that matters, but rather their relative weights. At the Γ point, our results are in good agreement with the laser ARPES experiments. As shown in Fig.

Table II. Fit parameters of the 4-band BCS model. The pairing potentials are given in units of V 1 (with V 1 =V (0) αα =V (0) ζζ =V (0) ββ and V 2 =V (0) αζ =V (0) αβ =V (0) ζβ ). Comparison of the energy gaps derived from laser ARPES (at 1.5 K), SANS (at 0.1 T), with those obtained from the 4-band BCS analysis of the zero-field heat capacity. The gaps are given in units of k B T C . For the ζ and bands, the mean values are given. Assignments for the SANS data is arbitrary. 6(b), the larger (smaller) of the 3 energy gaps is found on the α (β) band while the ζ gap exhibits accidental nodes which arise from the angle-dependent interband interaction with the pocket. These results are at odds with all theoretical calculations [55][56][57][58] which predict the largest gap on the β band. Moreover, at the X point, we recover the large gap ∆ =1.9 k B T c inferred from our magnetization and field-dependent heatcapacity data. This gap was not observed in any other experiments and is due to a significantly larger intraband constant on the band (see Table II). We stress that our analysis is not unique and other sets of parameters could fit C e (T) equally well. However, they would result necessarily in gap amplitudes close to the values we obtain because we constrain the individual densities of states to match approximately the values inferred from dHvA and ARPES. Moreover, as illustrated in Fig. 6(a), each gap has its own role in the temperature dependence of C e (T) in the superconducting state. Indeed, the gap alone is responsible for the jump at T c and has a vanishing contribution at low temperatures, while the β gap is predominantly responsible for the hump observed around T/T c ≈ 0.2. However, its decrease below this temperature is too steep to reproduce the experimental data. It is smoothened by the nodal contribution of the ζ gap. The slight maximum in the contributions from the α and ζ bands moreover lessen the dip of the high-temperature side of the shoulder due to the β band.

V 2 V (0) V (0) α V (0) ζ V (1
Although our results do not bring direct evidence of the existence of line nodes, they firmly establish the existence of tiny energy gaps with ∆/k B T c < 1.0. Their small amplitude imposes the requirement of cooling the sample below 80 mK to be able to observe the linear nodal behavior, as shown in the inset of Fig. 6(a). To our knowledge, no measurements in this temperature range were ever reported. Thus, our results do not exclude d xy symmetry. On the other hand, the agreement with laser ARPES is only qualitative. As shown in Table II temperature of the latter is 10 times larger. 8, 59

Conclusions

We have shown the existence of a relatively large nodeless energy gap of amplitude 1.9 k B T c that excludes the possibility of d x 2 -y 2 symmetry for the superconducting-state order parameter in KFe 2 As 2 . Our results do not bring direct evidence for line nodes, they clearly prove the existence of tiny energy gaps (∆/k B T c < 1.0) which strongly govern the low-field and low-temperature heat capacity, much like MgB 2 . Furthermore, the small amplitudes of the gaps indicate that very lowtemperature measurements (T < 80 mK) will be required in order to observe the possible signatures of line nodes in this compound; a restriction that also applies to other probes like penetration depth and heat transport. Our results shows qualitative agreement with recent laser ARPES measurements, and strongly suggest the superconducting-state symmetry to be swave.

Fig. 1 .

 1 Fig. 1. (Color online) (a) Schematic Fermi surface of KFe 2 As 2 inferred from dHvA and ARPES measurements. 15, 16 (b) Possible symmetry of the superconducting-state order parameter of KFe 2 As 2 (only one band is shown at the Γ point).

Fig. 2 .

 2 Fig. 2. (Color online) (a) Heat capacity of KFe 2 As 2 . The inset is a close-up of the low-temperature region showing the high-temperature tail of a Schottky anomaly. (b) Zero-field electronic heat capacity C e of KFe 2 As 2 . The green line is the weak-coupling BCS heat capacity for an s-wave superconductor (∆ 0 =1.764 k B T c ).

Fig. 3 .

 3 Fig. 3. (Color online) (a) Raw magnetization curves of KFe 2 As 2 measured at several temperatures for H || c. (b) Reversible magnetization curves. (c) Difference of the superconducting-and normal-states reversible magnetizations. (d) In-plane superfluid density derived from the reversible magnetization curves. The solid line is the superfluid density of an s-wave gap with ∆ 0 = 1.9 k B T c .

2 Fig. 4 .

 24 Fig.4. (Color online) (a) Field dependence of the heat capacity of KFe 2 As 2 (blue symbols) at T = 0.12 K for H || c. The magenta line is the theoretical curve of the mixed-state heat capacity of an s-wave superconductor47 normalized by the density of states of the band (see TableI). The green curve is the resulting contribution of the 3 small energy gaps α, β and ζ obtained by subtracting the heat capacity of the band from the data. (b) Close-up of the low-field region. (c) Field dependent vortex core size derived from the mixed-state heat capacity.

Fig. 5 .

 5 Fig. 5. (Color online) Schematics of the 4-band BCS model used to analyze C e (T, 0).

Fig. 6 .

 6 Fig. 6. (Color online) (a) Temperature dependence of the heat capacity of KFe 2 As 2 derived in the 4-band BCS model (black lines). Contributions from the individual bands are also given. The inset shows T/T c < 0.1 K. (b) Angular dependence of the individual gaps at T = 0 K. (c) Temperature dependence of the indivdual gaps for θ = 0.

Table I .

 I Parameters derived from the dHvA and ARPES measurements assuming 2D Fermi-surface sheets, with γ

  , Okazaki et al. reported overestimated gap values in comparison to heat-capacity and SANS measurements. Particularly, they find ∆ α = 3.8 k B T c , which is conspicuously comparable in amplitude to the largest gap observed close to

We thank J.-P. Brison, M. Lang, M. Ichioka, K. Machida, A. Chubukov, T. Shibauchi and S. Kittaka for stimulating and enlightening discussions. This work was supported by the Deutsche Forschungsgemeinschaft through DFG-SPP 1458 "Hochtemperatursupraleitung in Eisenpniktiden". The work performed in Grenoble was supported by the French ANR Projects (SINUS and CHIRnMAG) and the ERC starting grant NewHeavyFermion.


K Bougot-Robin

J.-P Hugonin

M Besbes

And H Benisty
email: henri.benisty@institutoptique.fr

Broad working bandwidth and "endlessly" singlemode guidance within hybrid silicon photonics

The successes of nonlinear photonics and hybrid silicon photonics with a growing variety of functional materials entail ever-enlarging bandwidths. It is best exemplified by parametric comb frequency generation. Such operation challenges the dielectric channel waveguide as the basis for guidance, because of the adverse advent of higher order modes at short wavelengths. Surprisingly, the popular mechanism of endlessly single-mode guidance [Opt. Lett. 22, 961 (1997).] operating in photonic crystal fibers has not been transposed within silicon photonics yet. We outline here the strategy and potential of this approach within planar and hybrid silicon photonics, whereby in-plane and vertical confinement are shown to be amenable to nearsingle-mode behavior in the typical silicon band, i.e., λ 1.1 μm to ∼5 μm. © 2015 Optical Society of America http://dx.doi.org/10.1364/OL.40.003512

The leveraging of nonlinear photonics and hybrid silicon photonics is well exemplified by comb-frequency generation through parametric [1,2] or cascaded interactions [3]. These advances entail large frequency bandwidths, requiring clever dispersion management [4][5][6][7] on about two octaves (λ 1.1 μm to ∼5 μm), and can benefit from functional materials hybridized onto a silicon or silicon-on-insulator (SOI) platform, such as III-V semiconductors (see [8] and references therein). Such broadband regimes challenge the dielectric channel waveguide as the basis for guidance: depending on the guide dimension, the risk arises of perturbing higher order modes at shorter wavelengths or deconfinement at longer ones.

We seek inspiration from the popular mechanism of "endlessly single-mode guidance" [9] of photonic crystal fiber (PCF), Fig. 1(a), to address this issue. The crucial point is that the cladding, carved in the same material as the core, is designed with periodic holes. Consequently, high-frequency modes see it as a higher effective index medium than low-frequency ones: the ability of these modes to squeeze their lobes repulsively off the low index inclusions (for PCF, air holes), raises the highest effective longitudinal-propagation constant of the cladding. For higher modes in a large core, Fig. 1(a), this suffices to provide leakage channels and maintain single-mode behavior "endlessly." This simple scheme did not find a genuine one-dimensional equivalent, in spite of clever early attempts [10]. There is a rich index map for integrated optics, whereby three refractive indices exist in any realistic channel/rib guidance (substrate, core, and cover). This blurs attempts to transpose the teaching of PCF whose cladding exploits only silica and air.

In this Letter, we focus on the SOI whose "index engineering" is mature [11]: we first show that the application to the SOI geometry of the "endlessly" single-mode guidance concept works with "endless" guidance on a broad range of about two octaves. We then argue that the addition of a hybrid stack based on III-V semiconductors can play most of the desired lightmatter interaction role, by relaxing the influence of vertical higher order modes if the layers are thick enough to accommodate the largest wavelengths of interest (∼5 μm).

We analyze guidance in structures such as that of Fig. are not directly related to photonic-crystal waveguides but are nevertheless the logical transposition of the fiber concept to an integrated-optics concept. That of Fig. 1(b) has an index map nx; y, and classical Si thicknesses for operation around 1.5 μm. We use a variant of the rigorous coupled-wave analysis (RCWA) [12][13][14][15][16][17][18][19], whereby the finite vertical dimension of the system, x, is padded with perfectly matched layer (PML) mimicked by a complex-coordinate transform [12], being hence equivalent to infinite claddings. Technically, the use of this transform along x with its reduction in computational domain size limits the need of higher Fourier terms (typically 35 terms here). We can also get the field lateral profile without an extended ydomain, thanks to Bloch boundary conditions [13]: the field Bloch parameters are known at the start of the periodic domain, and the profile can be obtained up to arbitrary y in the periodic domain with the solution, even in the case of slow or absent decay. Numerically, the aperiodic Fourier modal method was benchmarked as very accurate in [18], see also [14]. Material dispersion is taken into account. We take the tabulated data of Si, n Si λ and silica n silica λ as in [4]. We study the resonant response to a point source S to have a picture of guidance and leakage, following [16][17][18]. The system is pictured in Fig. 2(a). The periodic cladding of period h determines leakage through its band structure. It can be studied as ω ωK y ; γ, Fig. 2(b), with K y the transverse vector, γ the longitudinal wavevector along z, and ω the frequency (vacuum wavelength λ 0 ). At a given frequency and for a given mode branch, each interval γ min ; γ max reached at said frequency defines an allowed cladding band. Our goal is to maximize such bands to ensure leakage for all but a fundamental mode of the whole guide. The source S of Fig. 2(a) is either of three electric and three magnetic dipoles along the three axes, addressing any mode symmetry. Compared to [19], we rather impose the frequency ω and scan γ to find the scattering matrix poles. This matrix is derived from the source through the core part to each side of the periodic cladding. Precisely, by calculating the field's imaginary part jImfieldSj at the source location S, we get profiles of resonances versus γ such as Fig. 2(c) for the six source terms (in log scale), and by the same token, information on the bands as discussed later.

We start with a 50/50 duty cycle and a cladding period h 0.6 μm, fairly larger than a typical half-wavelength at λ 0 1.5 μm in Si or also in planar SOI with typical n eff 2.8: λ 0 ∕2n eff 0.268 μm. The silicon thickness is here d 500 nm, etched by e 250 nm groove. The defect size is b 1.2 μm, large enough to otherwise accommodate a quasi-TE mode up to large wavelengths (λ 4.5 μm) without too much deconfinement. In Fig. 2(a), we fill the grooves with a composite air-silicon material of index n EMTd λ, EMT referring to "effective material theory" [8,[20][21][22], and d to "dispersive": a broad range of indices can be made by nano-pattering silicon at modestly sub-wavelengths scales. We refer for simplicity to these dispersive media by the index n EMT for n Si 3.5 exactly 3.5 n Si 1.343 μm, i.e., we use n EMTd λ 1 n EMT -1n Si λ -1∕3.5 -1 with fixed n EMT . Note that a dispersionless EMT material gives a very similar behavior and a fully dispersionless stack still quite so. Conversely, an accurate modeling is nevertheless relevant for nonlinear comb generation [1][2][3].

Figure 2(c) holds all we need: we checked that the claddingallowed bands tightly correspond to gray areas of Fig. 2(c), with higher response to the source. In forbidden bands, except discrete mode poles, the response is low as the energy leaving the system vanishes. Thus, truly guided modes occur at peaks within the white areas signaling forbidden bands. Concerning our quest of "endlessly" single mode guidance, while the TE0 mode shown in Fig. 2(e) fits in the defect, we see in Fig. 2(c) how the limited extent of allowed bands leaves enough room for our large defect to support a second TE1 guided mode, shown in Fig. 2(d). We now want to extend allowed bands till such short-wavelength TE1 modes merge in them.

In Fig. 3(a), we study at λ 0 1.5 μm a 25/75 duty cycle with still n EMT 1.5. We, hence, extend the allowed band up to γ 13.1 μm -1 instead of 12.75 μm -1 in Fig. 2(c), thus, up to a higher effective index n eff γλ 0 ∕2π. Correspondingly, only a TE0 mode lies in the upper forbidden area, the TE1-like peak being repelled into the continuum at γ < 13.1 μm -1 . The TE0 mode profile in Fig. 3(c) fully resembles that of Fig. 2(e): as expected, a relatively large mode area provides this insensitivity.

An unbalanced duty cycle is not the only way to weaken periodicity: a lower index contrast also does. We thus try a rather high EMT index n EMT 2.9 keeping a 50/50 duty cycle, in Fig. 3(b). The allowed band now reaches nearly γ 13.25 μm -1 . Besides the single mode in this area, the TE2 mode now also has a very narrow window between allowed bands. Together with the large offset in γ, this makes unlikely the coupling to TE2, so we effectively approach a single-mode operation.

Checking how much room is left for suppressing TE1 is the issue addressed by Fig. 3(d). For the n EMT 2.9 [b] cladding of Fig. 3(b), we are at the vertical dashed line. By pushing up the frequency to λ 0 1.42 μm, we facilitate squeezing the TE1 mode in the large defect. As shown in Figs. 3(e)-3(f ), we skirt the single-mode limit there: the TE1 mode is extremely deconfined, while the TE0 mode is as confined as in Fig. 2(e). By performing a similar study at several n EMT values from 2.8 to 3.4, we outline the domain of TE single-mode operation in Fig. 3(d).

We now visualize what the TE0 mode looks like at large wavelengths where a single-mode behavior is natural, exploring mode profiles up to λ 0 4.5-5 μm in a hope to cover 1.5 octaves. Beyond ∼5 μm, the silica absorption onset of SOI is detrimental. In Fig. 4, we depict the TE0 mode for the same configurations (a), (b) as seen in Fig. 3. The mode has swollen similarly in both cases (we normalized the color map to silicon's strongest field to ease comparison). There is more confinement in (a) (cladding[a]) as could be expected by the larger forbidden band. We are thus reaching a vertical deconfinement regime, not so surprisingly since the silicon slab with d 500 nm λ 0 ∕9 is below the typical λ 0 ∕2n 650 nm limit.

By a spectral quest of the fundamental mode, we can draw the phase dispersion of the TE0 mode, n eff γλ∕2π, Fig. 4(c), as a function of wavelength (upper axis) and frequency ν c∕λ in THz (lower axis). Of great interest comb and parametric generation is the group index n eff νdn eff ∕d ν, it is relatively stable around 3.5-3.9 over the two octaves.

We can now wonder whether we can compensate also for the vertical swelling (as well as manage dispersion) by getting closer to the original fiber endlessly single-mode concept: to this end we cannot afford a thicker silicon (a vertical TE1 mode would appear), but we can attempt to surround the central core by "allowed modes" in the z-direction as well, not only in the SOI xy plane. This is similar to the low-index "horizontal slot" of [5][6][7], albeit with less index contrast: this option is made possible in the moderate index contrast regime by the advent of hybrid oxide-free InP-on-Si silicon photonics, whereby bonding of a III-V heterostructure (InP-based) to nanostructured silicon has been demonstrated [20,21].

The stack we explore is that of Fig. 5(b), with a thicker silicon slab (700 nm) with 50/50 duty cycle grooves identical in shape to those of [b] above. It consists of an extra 400-nm-thick InP layer (index n ≈ 3.17 [20]), with a 300-nm top-quaternary layer (GaInAsP alloy with bandgap wavelength < 1.42 μm) having the same index as Si for simplicity. No patterning is made in these layers. To compensate the expected TE0 profile "swelling" toward the III-V and maintain in-plane confinement, a lower n EMT 2.6 instead of 2.9 is adopted. , but it is quite deconfined versus TE0. As for the apparently localized "vertical TE2 mode," Fig. 5(h), we checked that it leaks as well. Finally, we did a scan for this structure and got the effective indices of Fig. 4(d) (thin lines). We note a weaker phase-index dispersion and more regular group-index dispersion, with very weak second-order terms.

In this Letter, we do not explore the full parametric study that would map, for instance, the exact onset loci of higher order modes in the rich parameter space, as this obviously depends on application and technology-related choices. We believe that the more important point is to expose an embodiment strategy for the "endlessly" (broad working range) singlemode-guidance principle, first in a single transverse direction, but also in the "vertical" direction, similarly to ongoing efforts [5][6][7]. Even if our strategy a mitigation rather than a total cut-off of the higher-order modes, it would have a clear interest given the large sensitivity of nonlinear processes to overlap integrals underlying light-matter interactions.

To conclude, we have successfully exported to the "silicon realm" the strategy of "endlessly" single-mode guidance developed in PCF. We adapted the clever principle of these PCF studies concerning the shaping of large "escape windows" in the cladding. In terms of spectral width, a broad working range was obtained, from ∼5 μm the onset of silica absorption, to ∼1.15 μm the silicon bandgap. With the rest of the photonic toolbox, planar-based channel guide geometries can challenge current fiber functions. Our strategy holds for nonlinear hybrid as well as standard silicon photonics.

Funding. French Agence Nationale de la Recherche (L'Agence Nationale de la Recherche) (COHEDIO).

Fig. 1 .

 1 Fig. 1. (a) Photonic crystal fiber (PCF) and how a higher order resonance leaks into a high-effective-index cladding mode whose lobes are squeezed off air holes. (b) SOI structure with central ridge and laterally, periodic ridges. (c) SOI+III-V structure.

Fig. 2 .

 2 Fig. 2. System studied, with (a) source S and periodic cladding; (b) generic cladding band structure of a given mode ω ωK y ; γ in first Brillouin zone K y , parametrized by γ: at a given frequency, an allowed γ band arises in an interval γ min ; γ max . (c)-(e) λ 1.5 μm study of the basic structure with n EMT 1.5 and 50/50 duty-cycle. (c) Modal behavior, seen through the imaginary field response jImfieldSj at point S to the six kinds of source terms. The gray areas with always high response are photonic bands. The area of low response are forbidden bands with localized modes. (d) TE1 mode profile, E y component; (e) TE0 mode profile, same. Note in (d) that S stands at a zero of E y , hence the distinct source responses versus TE0.

6 E 3 .

 63 EM MT T i in nd de ex x cut-off wavelength (µm) (a) Modal study with only electrical response shown for 1.5 μm and claddings with 75/25 duty cycle and n EMT 1.5; (b) same for 50/50 duty cycle and smaller contrast, n EMT 2.9; insets [a,b] on the right side show sketches of these cladding; (c) mode TE0 in the [a] cladding case looks identical with TE0 of more confined structures; (d) TE1 cut-off wavelength of structures with variable n EMT , the value n EMT 2.9 of [b] cladding is indicated; (e), (f ) cut-off situation at 1.42 μm, the second-mode TE1 is just above cut-off (profile with slow decay). This TE1 mode would thus leak in a cladding mode at 1.5 μm, as desired, but would remain guided at λ < 1.42 μm. The TE0 mode remains well confined.

Fig. 4 .

 4 Fig. 4. (a), (b) Modes at 4.5 μm for the [a,b] cladding of Fig. 3 are given in (a), (b), respectively. Note the pronounced field extrema within the EMT in (a). (c) Phase and group index dispersions of SOI (thick line) and hybrid (thin line) structures.

Figure 5 (

 5 Figure5(a) reminds the mode profile of the stand-alone SOI stack, whereas Figs. 5(c)-5(e) show the TE0 mode profile at 4.5, 2.5, and 1.5 μm. The avoidance of deconfinement is most visible on the silica side below silicon with a weaker tail there (for these maps, the coordinate-transform x boundaries are set beyond the map limits, while they are closer for extensive scans). And overall, there is a better overlap of the TE0 profiles at different frequencies. Both factors are positive toward enhanced nonlinear mid-infrared operation. As for the higher-order mode patterns, Figs. 5(f )-5(i) show a selection at the three same wavelengths as above, hence the rows [(c) and (f ); (d), (g), and (h); (e) and (i)]. The occurrence of a "TE1 vertical mode profile" is clear [Figs. 5(f ), 5(g), 5(i)]. But (1) its profile nests in the top 300-nm layer even up to 4.5 μm [Fig. 5(f )], and (2) it propagates laterally (allowed bands) except at short wavelength [Fig.5(i)], but it is quite deconfined versus TE0. As for the apparently localized "vertical TE2 mode," Fig.5(h), we checked that it leaks as well. Finally, we did a scan for this structure and got the effective indices of Fig.4(d) (thin lines). We note a weaker phase-index dispersion and more regular group-index dispersion, with very weak second-order terms.In this Letter, we do not explore the full parametric study that would map, for instance, the exact onset loci of higher order modes in the rich parameter space, as this obviously depends on application and technology-related choices. We believe that the more important point is to expose an embodiment strategy for the "endlessly" (broad working range) singlemode-guidance principle, first in a single transverse direction,

Fig. 5 .

 5 Fig. 5. Various modes and their y component. (a) Reference mode at 4.5 μm in SOI stack; (b) sketch of III-V on SOI; (c)-(e) mode in SOI+III-V stack at 4.5, 2.5 and 1.5 μm; (f), (g) the second "vertical 1" mode at 4.5 or 2.5 μm tends to peak in the high-index layer at shorter wavelengths, but has no in-plane confinement; (h) third vertical mode; (i) top-split lobe of the TE high-order mode at 1.5 μm.


M L Ju

H Jmal

R Dupuis

E Aubry

Visco-hyperelastic constitutive model for modeling the quasi-static behavior of polyurethane foam in large deformation

Keywords: Polyurethane foam, Quasi-static behavior, Hyperelastic behavior, Viscoelastic behavior, Identification method

Flexible polyurethane foam is widely used in numerous applications such as seats and mattresses, due to its low stiffness and its ability to absorb deformation energy. The main objective of this paper is to model the quasi-static mechanical behavior of three types of polyurethane foam in large deformation and to compare these three foams with three proposed models. The uniaxial compression/ decompression tests at three different strain rates were performed. The test results show that the three foams present different plateau stresses, maximum stresses and abilities to absorb energy. Moreover, polyurethane foam also presents a nonlinear hyperelastic behavior and a viscoelastic behavior in large deformation. Three visco-hyperelastic models which include a hyperelastic component and a memory component are proposed to model these behaviors. Model parameters were identified using the experimental data and a proper identification method. These models were validated on these three types of foam with the aim to present comparison results. The comparison results show that Ogden's viscoelastic model best agrees with the experimental results.

INTRODUCTION

Flexible polyurethane foams are widely used in a variety of engineering applications ranging from medical fields to sport, cushioning applications, apparel products, packaging, footwear and sound and vibration control [1,2]. Because polyurethane foam is a very important material in polymer engineering, it is essential to study its mechanical behavior to improve the product quality, which has originally led to the investigation presented in this paper.

The properties of polyurethane foams include light weight, low stiffness, excellent insulation abilities, a low Poisson ratio, low density (less than 80kg• m -3 for flexible foam), high compressibility, a slow recovery rate and the ability to absorb strain energy. Polyurethane foams are cellular materials and can be categorized as open or closed cell materials depending on the shape and the connectivity of the cells. If foams allow fluid to flow throughout the cellular structure, these foams are called open cell foams [3]. The variety of applications of open cell foams is enormous and the mechanical behavior of open cell foams requires better understanding in order to increase the number of applications. In this study, three types of open cell polyurethane foams are investigated.

The mechanical behavior of polyurethane foams under different strain rates and large strain levels is highly complicated [4,5]. The description of the mechanical response of foams required considerable time and effort. Konstantinidis et al. used two 3D analytical models to determine the mechanical behavior of closed and open cell foams [6]. Amaral-Labat et al.

presented the mechanical properties of heat-treated organic foams [7]. Zaretsky et al. studied the impact response of high density flexible polyurethane foams [8]. Saint-Michael et al. studied the mechanical properties of high density polyurethane foams [9]. Numerous other researchers worked on the mechanical properties of low density foams [10,11]. More details related to mechanical properties of foams are given by Gibson and Ashby [START_REF] Gibson | Cellular Solids:Structure and Properties[END_REF]. The stress-strain response of polyurethane foams shows a nonlinear hyperelastic behavior which can be described by a number of hyperelastic models [START_REF] Zhang | [END_REF]14]. In the literature, a variety of models have been proposed to fit experimental results for hyperelastic materials. Mooney proposed a model with two parameters [15]. The Neo-Hookean model has only one material parameter [16]. Rivlin made a modification of the Mooney model and proposed a general expression which is named Mooney-Rivlin model [START_REF] Rivlin | Philosophical Transitions of the Royal Society of London[END_REF]18]. The Blatz-Ko model [19] was also used to describe the hyperelastic properties of compressible materials and this model was revised and generalized by Willson and Myers [20,21]. In 1972, Ogden [22,23] proposed a strain energy function expressed in terms of principal stretches, which is a very general expression for describing hyperelastic materials. Other models include those of Yeoh [24], Beatty [25], Arruda-Boyce [26], Bischoff et al. [27], and Attard [28] also helped to predict the hyperelastic behavior of polyurethane foams. Besides hyperelasticity, polyurethane foams show strong viscoelasticity which is an important aspect of polyurethane foam behavior [29,30]. A material is called memory viscoelastic if the present state of response depends, not only on the present state of loading, but also on the previous states [START_REF] Ferry | Viscoelastic Properties of Polymers[END_REF][START_REF] Singh | [END_REF]. Polyurethane foams exhibit a significant hysteretic response during unloading after prior loading in uniaxial compression. The stress during unloading is significantly lower than the one during loading under the same strain. To describe this phenomenon, Yang et Shim [33] proposed a visco-hyperelastic model for foams under strain rates to capture the large three-dimensional compression behavior. In the work presented here, a memory model is proposed to study the viscoelastic behavior of polyurethane foams. This paper describes the mechanical behavior of three different polyurethane foams. The purpose is to combine a hyperelastic equation and a viscoelastic constitutive model to create a visco-hyperelastic model and to estimate the model parameters so as to allow a good correlation between the model and the experimental results. A non-linear identification method has been used out to obtain the parameters. The determination of model parameters is based on the use of test samples under a simplified strain state. The unknown model parameters are obtained using curve fittings from experimental data. The corresponding identification errors are also taken into account to analyze the results. Numerous compression tests were performed to obtain a sufficient number of experimental results for parameter identification. The determination of these parameters through foam testing is important for the provision of data for the simulation of their service performance, product design and quality control.

EXPERIMENTAL

Materials

Three types of polyurethane foam, which are usually used as core material in a sandwich structure, were considered and are referred to here as foam A, foam B and foam C. The properties of the three types of polyurethane foam are summarized in Table 1.

Cell type open open open

In order to study the stress-strain relation of these foams, a series of loading-unloading uniaxial compression tests were performed at a constant temperature and humidity. The polyurethane foam test specimens were cut from a foam block 2000mm×1200mm×75mm (length, breadth, thickness). The band saw method was used to cut the foam specimens in the present work. Following specific track, specimens were cut with desired dimension. All specimens were put on table for the rest of 24 hours to eliminate the biased effects which occur in the cutting process. The specimens were not presented exactly the same dimension in height. For this reason, the samples which have the difference lower than 2 mm from the desired height were chosen to perform the tests. All specimens were provided with the same mechanical and environmental histories. They were original specimens and each specimen was compressed and decompressed only one time.

Quasi-static experiment and device

All the tests were carried out with a standard "INSTRON 3300R4204" compression-tension testing device driven with BLUEHILL software (Figure 1).

Figure 1 experimental device

The device includes a basis frame and an upper block which moves vertically. Two compression plates were installed on this device: one on the base of the machine and the other on the force sensor of the crosshead. The two plates were carefully checked to be strictly parallel. Before the beginning of the experiments, the top plate had to be moved down slightly for full contact with the material because the top and bottom of each foam sample were not exactly parallel. In order to establish a homogeneous deformation field, the shear stresses between the two plates on the top and the bottom sides and the test specimen at the uniaxial compression test had to be ignored [34]. There are numerous norms for the tests of polyurethane foam and the norm 'ASTM D 3574' norm was used here as reference to perform the tests.

For the sake of minimization of the noise contribution, the maximum experimental response force of foam had to be slightly less than the load cell maximum capacity. The test conditions including the strain rate, the maximum compression level, the number of cycles, the sampling period and the test mechanical parameters exported were controlled using the BLUEHILL software configuration window.

At the beginning of the experiment, the samples were set into position on the bottom plate.

The test started when the upper plate began to affect the foams and then the upper block moved down to compress the samples to the final level. The final compression ratio was 80% of the original thickness. At the end of the loading process, the upper block took a converse compression direction and returned to its initial position. The software BLUEHILL recorded all test data such as time, stress, strain, force and displacement. Each specimen had been quasi-statically compressed and then decompressed at a constant speed during all the test process. At the beginning of the test, 15 test samples for each foam at three strain rates were taken for the preliminary test. Then, equation (22) helped to determine the minimum number of test samples which were to ensure the statistical quality of the parameters. The quasi-static compression test conditions and the minimum number of tests are given in Table 2; more details will be given in the discussion section.

Experimental results

Figure 1 shows the experimental stress-strain results of foam A, foam B and Foam C in test 1, test 2 and test 3, respectively. It can clearly be seen in Fig. 1a that the deformation of polyurethane foam has three regions which can be distinguished during uniaxial compression: initial elastic deformation (Region 1), collapse deformation (Region 2) and compaction deformation (Region 3). In the first stage, the polyurethane foam deforms in a linear elastic manner due to cell wall bending, which accounts for less than 5% of the entire deformation.

In the course of the second stage, there is a plateau of deformation at almost constant stress.

The cell walls which are like thin tubes or plates, lose their stability and cause large deformations. During this stress plateau phase, the polyurethane foam undergoes large compressive strains and absorbs a considerable amount of specific energy. In the third stage, a region of densification occurs where the cell walls crush together, resulting in a rapid increase in compressive stress. Figure 2 shows the experimental stress-strain results of three types of foam (foam A, foam B and foam C) in test 1. As can be seen, there is a significant stress difference between the loading and unloading processes corresponding to the same strain level for three foams. This phenomenon can be explained by the viscoelastic behavior of polyurethane foams. It also indicates that the area between the loading curve and unloading curve for the three foams, which represents the energy absorbed by foams, is different. This shows that the three foams have different abilities to absorb the strain energy.

FIG. 2. Three types of foam compression stress-strain results in test 1

Figure 2 also demonstrates the comparison of experimental results of three types of foam in test 1. It can be seen from this figure that foam A has both the highest plateau stress and the maximum stress of the three foams and it also has the biggest energy absorption ability. Foam B has a slightly higher plateau stress but has a lower maximum stress than foam C.

It can be seen in Fig. 1 and Fig. 2 that the strain rates affect the experimental curves and that the experimental stress-strain for the loading and unloading phases is quite different due to the viscoelastic properties. This means that the model parameters for the loading and unloading phases should be calculated separately.

From the above experimental results, the conclusion corresponds to the empirical theories: polyurethane foams exhibit both a nonlinear hyperelastic behavior and a viscoelastic behavior in the quasi-static compression test. It will be of help in the modeling of the mechanical behavior of polyurethane foams. The proposed model should consist of hyperelastic and viscoelastic behaviors. The modeling method is discussed in the following section, based on the experimental conclusion of the present paper.

MODELING

Polyurethane foams are assumed to be homogeneous with a constant-section in small deformation [START_REF] Dupuis | SEM XI International Congress & Exposition on Experimental & Applied Mechanics[END_REF]. This assumption was extended to large deformation by researchers who worked on the macroscopic memory models [START_REF] Deng | [END_REF]37]. According to the experimental results of this paper and the empirical results, polyurethane foams exhibit nonlinear hyperelastic and viscoelastic properties. In order to predict the mechanical behavior of polyurethane foam, a visco-hyperelastic model which can be decomposed into the sum of a viscoelastic component and a nonlinear hyperelastic component, is proposed here. The stress response of the viscohyperelastic model in uniaxial compression is represented as:

(

where is the entire stress response, and are the hyperelastic stress and viscoelastic stress, respectively. The nonlinear hyperelastic constitutive models which are commonly used in the study of rubber-like materials are considered to represent the hyperelastic stress here.

For the viscoelastic part, a hereditary integral model is used to describe the viscoelastic stress.

Hyperelastic component

Hyperelastic constitutive laws are usually used to model materials that respond elastically when subjected to very large strains. In the theory of hyperelasticity, there exists a strain energy function, denoted or , where F is the deformation tensor, where are the principal stretches and are given by , where are strains (for compression ), to describe the hyperelastic phenomenon [22,23,38]. In the present paper, the strain energy functions with principal stretches are used and the stress can be derived from the function by [39]: ( 2 ) where is the determinant of the deformation gradient. It can be expressed as: (3) where are the densities in the reference and deformed configurations respectively.

In this study, the deformation is uniaxial compression and the compression direction is also a principal direction of the deformation. For the development of the models, suppose that the principal stretch corresponding to (Y-coordinate) corresponds to the compression direction; so, the stress state can be deduced:

(4)

Ogden's model

The hyperelastic model developed by Ogden [22] is widely used for compressible materials and the stored energy, W, is expressed in the form of a series:

(

where are material parameters, is the volumetric function. For practical purposes, the sum in ( 5) is restricted to a finite number of terms, while, for consistency with the classic theory, the shear modulus of the material can be deduced with:

(6) Based on equations ( 2) and ( 5), the stress for the Ogden model can be obtained: (7) According to equation ( 4), the solution of and substitution for equation (7) gives

(8)

Polynomial model/ Mooney-Rivlin model

The polynomial hyperelastic material model is a phenomenological model of rubber elasticity. This model can also describe the elastic properties of materials. In this model, the strain energy density function is of the form of a polynomial in the two invariants of the left Cauchy-Green deformation tensor. For compressible materials, a dependence of volume is added and the strain energy density function for the polynomial model is: (9) where (10) F is the deformation gradient tensor. When N=1, the model is reduced to the compressible Mooney-Rivlin model: (11) For N=2, the energy potential is as follows: [START_REF] Gibson | Cellular Solids:Structure and Properties[END_REF] To allow a comparison with other models, the higher order model (N=3), although rarely used, is also considered here and has the following form: [START_REF] Zhang | [END_REF] The stress-strain relations can be derived for the Mooney-Rivlin model (N=1), the polynomial form with N=2 and N=3 by considering equations (2) (4) (11) (12) [START_REF] Zhang | [END_REF].

Reduced polynomial model / Neo-Hookean model

Particular forms of equation ( 9) can be obtained by setting specific coefficients to zero. If all with are set to zero, the reduced polynomial form is obtained:

By setting N=1, the Neo-Hookean form is obtained:

This form is the simplest hyperelastic model and is often used as a prototype for elastomeric materials.

With N=2, the reduced polynomial strain energy potential is represented as: (16) and with N=3, it can be written as:

The stress strain relation is easily obtained with equations (2) (4) (15) (16) [START_REF] Rivlin | Philosophical Transitions of the Royal Society of London[END_REF].

Viscoelastic component

In the literature, a variety of macroscopic models are used to describe the viscoelastic behavior, such as the Maxwell model, the Kelvin Voigt model and the Zener model [40]. In this paper, the viscoelastic stress is represented by a hereditary integral model with different loading and unloading parameters: (18) where T is the test period and are the viscoelastic parameters in the loading and unloading phases, respectively. is a complex number indicating the viscoelastic residue. It is expressed in . represents the viscoelastic mode.

It is expressed in . This model takes account of the memory criterion of polyurethane foams and indicates the viscoelastic properties during the compression test process.

Entire stress response

The hyperelastic models describe the nonlinear elastic behavior and the hereditary integral model indicates the viscoelastic behavior of polyurethane foams. Combining the elastic component and the viscoelastic component, the global model of polyurethane foams undergoing large compressive deformation is given by equation (1).

Identification method

In the quasi-static compression test, the strain rates are a constant in the loading and unloading process. So, the relation between strain and time is linear and the curve obtained is as an isosceles triangle. The function of in the loading and unloading phases in equation ( 18) can be defined as: (19) where is the maximum strain.

Based on the symmetry of the elastic stress during the loading and unloading phases, the difference-stress method is used to identify the viscoelastic parameters [40]. According to equation (18), the difference stress between the loading and unloading phases is given by: (20) where and are the loading time and unloading time corresponding to the same strain level, respectively. Equation (20) helps to identify the viscoelastic parameters. After obtaining the viscoelastic stress in the loading and unloading phases, the elastic stress is determined using the hyperelastic models mentioned above. The identification method is summarized in Fig 3.

FIG. 3. Identification method

Parameter optimization approach

Optimization methods are used as basic tools to identify both viscoelastic parameters and elastic parameters. The stress-strain experimental data can be obtained from the uniaxial compression tests. The mean square error between the experimental data and model data can be calculated using a function : (21) where is the model stress and is the experimental stress. The objective of the optimization is to find the best combination of model parameters which minimize function . MATLAB's optimization toolbox is considered here for minimization. There are deterministic and random methods for optimization The trust region reflective, Levenberg-Marquardt and Gradient methods are three examples of deterministic methods which are effective when the objective function (function to optimize) changes rapidly. In the present paper, the trust region method is used as the option for the descent algorithm. The best parameters must be chosen to minimize the least mean square error between experimental and analytical data and satisfy the boundary condition, using the LSQNONLIN solver with tight stopping criteria.

MODEL RESULTS

Viscoelastic results

The difference-stress method is used here to identify the viscoelastic parameters, and then the elastic parameters are determined following the results of viscoelastic stress. In the study by Ippili R.K. [37], the viscoelastic parameters depend on the test conditions. But the present paper considers the dimensionless viscoelastic parameters to predict the foam behaviors under any test conditions (Table 3)[40]. TABLE 3. Relationship between dimensional and dimensionless viscoelastic parameters

Dimensional viscoelastic parameters Dimensionless viscoelastic parameters

With equation ( 20) and the identification method mentioned above, the dimensionless viscoelastic parameters values are summarized in Table 4. With the viscoelastic parameters in Table 4, the difference-stress results of the three foams in The viscoelastic stress results for the three foams in test 1 are given in Fig. 5. It can seen from these comparison results that foam A has the highest viscoelastic stress of the three foams and the viscoelastic stress of foam B is higher than that of foam C. It shows that the viscoelastic stress follows a reverse orientation of the global stress. This describes the relaxation effect on the global behavior of foams.

FIG. 5. Viscoelastic stress results of the three foams in test 1

Ogden-viscoelastic results

As viscoelastic stress has already been obtained with the model (Table 4, Fig. 5), the elastic behavior described by the Ogden model is now considered. The simulation of three foams at three strain rates was performed with different N values (N=1, N=2 and N=3). The comparison results (Fig. 6) help to determine the Ogden model order in equation ( 5). As shown in Fig. 6, the model prediction results become better and better with the increase of N values. For a precise analysis, the errors of the Ogden viscoelastic model for the three foams, at three strain rates and for different N values, are considered and summarized in Table 6. According to Fig. 6 and Table 6, it is clear that the three-term expression of the Ogden hyperelastic model is required. The results of the Ogden elastic parameter for the three foams in three tests are given in Table 7.

Polynomial-viscoelastic results

The polynomial model (equation 9) is also used here to predict the elastic properties of the foams. Using the same viscoelastic stress determined by the hereditary integral model, the results of the different orders of N of the polynomial-viscoelastic model are given in Fig. 7. As shown in Fig. 7, the Mooney-Rivlin model (N=1) cannot reach the experimental test results, which indicates that this model cannot predict the elastic behavior for the three polyurethane foams. On the other hand, it can be seen from Fig. 7 that the polynomial model with N=2 and N=3 well agrees with the experimental results. The relative errors of the polynomial model with N=1, N=2 and N=3 for three foams in three tests are given in Table 8. According to Table 8 and Fig. 7, the polynomial viscoelastic model of the order N=2 and N=3 better agrees with the experimental data than that of the order N=1(Mooney Rivlin model). Table 8 also shows that the errors between form N=2 and form N=3 do not present a great difference. Moreover, form N=3 has more parameters (12 parameters) than form N=2 (7 parameters). Considering the number of parameters and the accuracy of the simulation results, the polynomial model of the order N=2 is used and the values of the model parameters estimated from the compression test are presented in Table 9.

Reduced polynomial viscoelastic results

The reduced polynomial model is also considered here to simulate the elastic properties of three polyurethane foams. It can be seen in Fig. 8 that the Neo-Hookean model (N=1) does not agree with the experimental curves. With the order N=2, the model results cannot give an acceptable comparison either. On the contrary, the model with the order N=3 is in good agreement with the experimental results. Table 10 shows the errors of the reduced polynomial-viscoelastic model for the three foams in three tests with different values of N.

DISCUSSION

Validation of the parameters

As mentioned in the experimental section, it is necessary to determine the minimum of N test samples (Table 2) to ensure the statistical quality of all identified parameters. This quality is reviewed through the set at a 95% confidence level and the statistical limit error "SLE" shall not exceed 10%. In order to obtain the minimum number of test samples for each strain rate in the three different polyurethane foams, the following equation is required [40]:

where and standard deviation values and the estimated average corresponding to th parameters. These two parameters are obtained from the identification results of the preliminary test with 15 test samples. is a coefficient which can be determined with the probability table suitable for the estimated probability law of th parameters. Ceil(x) is a function to round to the smallest integer not less than x in MATLAB. The Student law is used for all parameters in this paper. In the identification method section, a difference-stress method is used to identify the viscoelastic parameters based on the symmetry of the elastic stress. So, it is necessary to check this assumption. Figure 9 indicates that the model elastic stress verifies the elastic symmetry. Except for the standard deviation analysis, the model results also show a good agreement with the experimental data. According to Fig. 10, the Ogden-hyperelastic model curves lie between the minimum and maximum experimental envelopes. This is another criterion for the validation of the parameters.

Comparison of three models

In the present paper, three visco-hyperelastic models have been proposed. According to the results of Fig. 6 and Table 6, three-term expression is required for Ogden's viscoelastic model. Figure 7 and Table 8 indicate that the value of N=2 of the polynomial viscoelastic model is considered here to predict the mechanical behavior of polyurethane foams. The order of N=3 of the reduced polynomial model is used here according to the comparison results in Fig. 8 and Table 10. Considering the number of model parameters and the errors between the model results and experimental results (see Table 12), Ogden's viscoelastic model with form N=3 has the least number of parameters (6 parameters) and the minimum mean errors for the three foams in three tests. Therefore, the Ogden viscoelastic model is more efficient than the other two models.

Comparison of three polyurethane foams

Finally, Ogden's viscoelastic model has been chosen to perform the comparison results of three different types of foam (Foam A, Foam B and Foam C). The comparison of the results of the viscoelastic stress (for these three foams) which were determined with the hereditary integral model, is given in Fig. 5.The comparison response of the elastic stress of the three foams obtained with Ogden's model in test 1 is given in Fig. 9. It can be seen in Fig. 9 that Foam A has the highest elastic stress of the three foams. Foam B has a slightly higher elastic response than Foam C. From the comparison results, it appears that Ogden's viscoelastic model proposed in this paper can be used to predict the mechanical behavior of different types of polyurethane foam with different densities.

CONCLUSION

This paper has presented numerous experiments with three different polyurethane foams at three different strain rates for loading-unloading uniaxial compression tests. Three viscohyperelastic models have been proposed for the description of the elastic behavior of polyurethane foam and the hereditary integral model is used here to predict the viscoelastic behavior of foams. A difference-stress method based on the elastic symmetry between the loading and unloading phases helps to identify the viscoelastic dimensionless parameters which do not depend on the test conditions. The results of the hereditary integral model agree with the experimental results particularly well, which shows that the model can predict the viscoelastic behavior of foam. Three different strain energy density functions for compressible materials have been applied to predict the hyperelastic behavior of the foams. These models consist of an incompressibility component and a compressibility component and the stress is derived from the model in terms of principal stretches. Several experimental data for compressible polyurethane foams have been used to identify the model parameters. The predictions based on the proposed visco-hyperelastic constitutive model compare well with the experimental data. The errors between the model and experimental data have been calculated. The results show that the Ogden-viscoelastic model with a three-term expression, the polynomial viscoelastic form with N=2 and the reduced polynomial viscoelastic model with N=3 can predict the quasi-static mechanical behavior of polyurethane foams even under large strain compression. Among the three models, Ogden's viscoelastic model provides the best comparison between the errors and the number of parameters. The first advantage of this model is that the dimensionless viscoelastic parameters do not depend on the test conditions. These parameters can be used to predict the viscoelastic behavior of polyurethane foam in other different test conditions, such as different strain rates.

The second advantage is that the model cannot only be used to predict the mechanical behavior of polyurethane foam, but also to predict the hyperelastic behavior of materials, such as rubber-like materials, or viscoelastic behavior materials or both of these two behavior materials. Based on these two main advantages, this model will be used for further research work.

FIG. 1 .

 1 FIG. 1. Experimental stress-strain results: (a) Foam A in test 1, (b) Foam B in test 2, (c) Foam C in test 3

three tests are shown in Fig. 4 .FIG. 4 .

 44 FIG. 4. Model results of difference-stress between the loading and unloading phases of the three foams in test 1 The viscoelastic relative errors of the three foams in three different tests are given in Table5. TABLE5. Difference-stress errors for three foams at three strain rates

FIG. 6 .

 6 FIG. 6. Comparison of the Ogden-viscoelastic model results with N=1, N=2 and N=3.

FIG. 7 .

 7 FIG. 7. Comparison of Polynomial viscoelastic model results with N=1, N=2 and N=3.

  FIG. 8. Comparison of the reduced polynomial viscoelastic model results with N=1, N=2 and N=3

FIG. 9 .

 9 FIG. 9. Comparison of the elastic results of Ogden's model with the three foams in test 1

Figure 10

 10 Figure 10 shows the global results of Ogden's viscoelastic model. From the figure, the model results agree with the experimental results for the three different foams in three different tests. The model response lies between the maximum and minimum envelopes.

TABLE 1 .

 1 Chemical and morphological properties of the three types of polyurethane foam

		Foam A	Foam B	Foam C
	Foam type	Flexible polyurethane	Flexible polyurethane	Flexible polyurethane
		foam	foam	foam
	Isocynate	Toluene diisocynate Toluene diisocynate TDIToluene diisocynate
		TDI		TDI
	Polyol	Polyether	Polyether	Polyether
	Expansion gas	CO2	CO2	CO2
	F a b r i c a t i o n	Free rise	Free rise	Free rise
	process
	Density	26 kg/m 3	39 kg/m 3	47 kg/m 3
	Cell population	2.94/mm 2	3.01/mm 2	3.27/mm 2
	density
	Density standard 0.67 kg/m 3	0.64 kg/m 3	0.84 kg/m 3
	deviation
	Average cell size 882 µm	531 µm	457 µm

TABLE 2

 2

	. Quasi-static compression test conditions and minimum number of test samples (N)
					minimum number of tests N minimum number of tests N minimum number of tests N
					Foam A Foam B Foam C
	Test 1 1	1.06 10 -2 0	80	150	17	17	15
	Test 2 1	5.33 10 -3 0	80	300	37	37	26
	Test 3 1	6.66 10 -4 0	80	2400	15	17	17

TABLE 4 .

 4 Results of dimensionless viscoelastic parameters

		Load parameters Load parameters Load parameters Load parameters			Unload parameters Unload parameters Unload parameters Unload parameters
	Test 1 41.72	0.174	-219.81	22023.94	42.64	0.178	-223.62	22027.07
	Foam A Foam A Test 2 40.59 Foam A	0.173	-215.79	22023.94	41.42	0.174	-217.36	22027.07
	Test 3 43.38	0.173	-214.46	22023.94	44.38	0.178	-218.35	22027.07
	Test 1 29.26	0.082	-147.49	22238.66	30.59	0.084	-150.01	22238.66
	Foam B Foam B Test 2 29.72 Foam B	0.081	-147.14	22238.66	30.99	0.082	-149.44	22238.66
	Test 3 31.84	0.096	-155.69	22238.66	32.89	0.097	-157.10	22238.66
	Test 1 37.72	0.094	-118.87	22024.80	38.51	0.099	-124.67	22027.34
	Foam C Foam C Test 2 41.35 Foam C	0.133	-153.02	22024.79	41.93	0.141	-159.06	22027.34
	Test 3 44.73	0.128	-136.02	22024.79	45.05	0.136	-141.78	22027.34

Table 5 .TABLE 5 .

 55 Difference-stress errors for three foams at three strain rates

		Foam A	Foam B	Foam C
	Test 1	5.77%	6.06%	6.92%
	Test 2	6.25%	5.32%	7.65%
	Test 3	6.79%	4.66%	7.97%

TABLE 6 .

 6 Ogden-viscoelastic model errors for three foams in three tests with different N

		Ogden viscoelastic model Ogden viscoelastic model Ogden viscoelastic model
	N=1	N=2	N=3

TABLE 7 .

 7 Ogden model elastic parameters for the three foams in three tests

		Test 1	2.81	1.66	-2.80	1.61	3.1E-3	38.28
	Foam A Foam A Foam A	Test 2	9.22	1.64	-9.21	1.63	2.8E-3	38.23
		Test 3	12.06	1.64	-12.05	1.63	3.0E-3	38.21
		Test 1	9.27	1.63	-9.26	1.62	1.1E-3	23.64
	Foam B Foam B Foam B	Test 2	11.66	1.61	-11.65	1.61	1.1E-3	17.74
		Test 3	7.27	1.63	-7.26	1.63	8.5E-4	45.75
		Test 1	11.13	1.76	-11.12	7.15	9.5E-4	21.73
	Foam C Foam C Foam C	Test 2	9.86	1.76	-9.85	1.76	6.9E-4	21.53
		Test 3	10.54	1.88	-10.53	1.87	2.4E-4	27.62

TABLE 8 .

 8 Errors of the Polynomial-viscoelastic model for the three foams in three tests with

	different N
				Polynomial-viscoelastic model Polynomial-viscoelastic model Polynomial-viscoelastic model
			N=1	N=2	N=3
		Test 1	15.3%	4.25%	2.85%
	Foam A Foam A Foam A	Test 2	14.3%		2.81%
		Test 3	15.5%	4.34%	3.07%
		Test 1	12.8%	5.27%	5.22%
	Foam B Foam B Foam B	Test 2	11.8%	4.89%	4.84%
		Test 3	10.6%	3.40%	3.41%
		Test 1	9.8%	2.57%	2.54%
	Foam C Foam C Foam C	Test 2	8.6%	3.01%	3.03%
		Test 3	6.6%	2.62%	2.63%

TABLE 9 .

 9 Polynomial model elastic parameters for the three foams in three tests

	C10	C01	C11	C20	C02	D1	D2
	Test 1 -2.02	1.22	-3.8E-2 -0.221	5.6E-3	1.182	0.251
	Foam A Foam A Foam A Foam A Test 2 -1.44 Foam A Foam A	0.86	-2.9E-2 -0.151	4.2E-3	0.879	0.205
	Test 3 -1.83	1.11	-3.4E-2 -0.199	5.1E-3	1.069	0.232
	Test 1 -0.29	0.17	-6.9E-3 -0.026	10 E-4	0.203	0.054
	Foam B Foam B Foam B Foam B Test 2 -0.23 Foam B Foam B	0.13	-5.7E-3 -0.021	8.3E-4	0.168	0.045
	Test 3 -0.05	0.03	-2.6E-3 -0.001	3.6E-4	0.069	0.023
	Test 1 -0.27	0.15	-6.4E-3 -0.023	9.2E-4	0.184	0.052
	Foam C Foam C Foam C Foam C Test 2 -0.15 Foam C Foam C	0.08	-4.6E-3 -0.011	6.5E-4	0.131	0.041
	Test 3 -0.07	0.04	-2.6E-3 -0.028	3.7E-4	0.076	0.026

TABLE 10 .

 10 Errors of the reduced polynomial-viscoelastic model for the three foams in three tests .8 and table 10, the form of order N=3 is required here for the prediction of the elastic behavior of polyurethane foams with the reduced polynomial viscoelastic model. The elastic parameters of the reduced polynomial model with N=3 are given in Table11.TABLE 11. Reduced polynomial model elastic parameters for the three foams in three tests

	with different N
			Reduced polynomial viscoelastic model Reduced polynomial viscoelastic model Reduced polynomial viscoelastic model
			N=1	N=2	N=3
		Test 1	23.0%	12.9%	5.33%
	Foam A Foam A Foam A	Test 2	22.5%	12.4%	4.64%
		Test 3	23.9%	13.4%	5.38%

TABLE 12 .

 12 Error results of comparison three visco-hyperelastic models

		Ogden's	Polynomial	Reduced polynomial
		viscoelastic model	viscoelastic model	viscoelastic model
	Number of parameters Number of parameters	6	7	6
	Test 1	2.57%	4.25%	5.33%
	Foam A Foam A Test 2 Foam A	2.96%	3.64%	4.64%
	Test 3	3.14%	4.34%	5.38%
	Test 1	4.74%	5.27%	5.69%
	Foam B Foam B Test 2 Foam B	4.55%	4.89%	5.26%
	Test 3	3.18%	3.40%	3.62%
	Test 1	2.34%	2.57%	3.01%
	Foam C Foam C Test 2 Foam C	2.92%	3.01%	3.26%
	Test 3	3.33%	2.62%	2.76%
	From

Table 12 ,

 12 it can be seen that the errors of the polynomial viscoelastic model (N=2) are smaller than those of the reduced polynomial viscoelastic model (N=3), but it has more model parameters. So, it is difficult to decide which one is better -the polynomial viscoelastic model of order N=2 or the reduced polynomial viscoelastic model with form N=3. According to the comparison of the three models and three-term expression, Ogden's viscoelastic model is the best model to predict the mechanical behavior of polyurethane foam.

40. Jmal H, Dupuis R, Aubry E. J Cell Plast., 47, 447-465 (2011).