M Zheng

Hongli Cathelinaud

X Ma

Xianghua Fan

Zhang
email: xzhang@univ-rennes1.fr

Shuo Chen

Xvsheng Qi

Zhuanghao Zhe

Michel Cathelinaud

Hongli Ma

Xianping Fan

Xiangh Zhang

Enhanced electrical conductivity and photoconductive properties of Sn-doped Sb 2 Se 3 crystals †

Enhanced electrical conductivity and photoconductive

properties of Sn-doped Sb2Se3 crystals

Introduction

Antimony selenide (Sb 2 Se 3 ) is an important V-VI binary semiconductor with excellent photovoltaic, photoconducting, photocatalytic and thermoelectric properties. [1][2][3][4][5][6] It has a high absorption coefficient higher than 10 5 cm À1 in the visible region and an optimum bandgap width of around 1.1 eV for single junction solar cells. In addition, it is stable and nontoxic and composed of earth abundant elements. These properties make it a promising candidate for various important applications in diverse areas such as solar energy conversion, photodetector technology, thermoelectric power generation, and opto-electronics in the near-infrared region. 5,[7][8][9][10][11] Despite the important research efforts, the progress in application development is still limited. One major issue with Sb 2 Se 3 is its intrinsic low electrical conductivity (10 À6 O À1 m À1 , in the bulk state), which results in low thermoelectric power factor, low spectral response and low external quantum efficiency. 12,13 A few studies to improve the electrical conductivity of Sb 2 Se 3 have been reported. For example, Choi et al. developed Sb 2 Se 3 nanowires decorated with Ag 2 Se nanoparticles, which increased the electrical conductivity of the Sb 2 Se 3 nanowires. 12 Zhang et al. explained an interpenetrating iodine-doped-Sb 2 Se 3 /Cu 2 GeSe 3 heterojunction network fabricated by controlling the crystallization of chalcogenide glass which can improve the photoelectric performance. 14 Our previous work including the synthesis of hybrid nanorods with the Sb 2 Se 3 /AgSbSe 2 heterojunction structure also shows a higher electrical conductivity and photodetector performance. 10 Despite these endeavors, the synthesis of Sb 2 Se 3-based crystals with controllable electrical conductivity has been rarely reported and still needs further exploration.

Doping is an efficient way to tune the electrical and photoelectrical performances of semiconductors. It has been widely utilized in the semiconductor industry. [15][16][17] For example, S-doped In 2 Se 3 nanowires show a much higher conductivity and spectral responsivity to visible light than undoped In 2 Se 3 nanowires. 17 A small amount of Sn leads to a significant increase in the electrical conductivity and photoconductivity of the amorphous As 2 Se 3 films. 18 Herein, we demonstrate a conventional melt-quenching method to controllably synthesize high-quality and variably doped (Sn x Sb 1Àx ) 2 Se 3 semiconducting crystals with the Sn ratio x ranging from 0.00 to 0.10. By incorporating additional Sn atoms into the binary Sb 2 Se 3 , a significant improvement in the electrical conductivity and

Experimental

Synthesis of (Sn x Sb 1Àx ) 2 Se 3 crystals A conventional melt-quenching method using high purity tin (99.99%), antimony (99.999%) and selenium (99.999%) as raw materials was used to synthesize (Sn x Sb 1Àx ) 2 Se 3 crystals. A schematic illustration of the preparation process is shown in Fig. S1 (ESI †). At first, a mixture of 12 g of raw materials according to the chemical composition (Sn x Sb 1Àx ) 2 Se 3 (where x = 0.00, 0.03, 0.05, 0.07 and 0.10, respectively) was loaded into a fully cleaned silica tube of 1 cm internal diameter. A vacuum trap that was placed in liquid nitrogen was used to condense any volatile products and also to improve the vacuum. The silica tube was sealed under a typical vacuum of around 10 À5 mbar and then was placed in a rocking furnace. The mixture was heated to 700 1C at a rate of 1.5 1C min À1 and held for about 8 hours to ensure complete melting and homogenization. Then it was cooled down to 650 1C during 2 hours and subsequent quenching in water at room temperature. Finally, in order to reduce the stress generated from the thermal shock in the quenching process, the as-synthesized crystals were annealed at 350 1C for 3 hours.

Characterization

X-ray diffraction (XRD) analysis was carried out using a PANalytical's X-ray diffractometer with a step width of 0.021 and a counting duration of 0.13 s per step. Scanning electron microscopy (SEM) images and the corresponding Energy Dispersive Spectrometer (EDS) elemental mappings were obtained using a JEOL JSM-7100F thermal field emission electron microscope equipped with an EDAX instrument. TEM, HRTEM and SAED images were obtained using a JEOL JEM-2100 field emission transmission electron microscope. The bulk sample was milled and then dissolved in chloroform before being placed on the surface of a copper grid for further TEM measurements. The quality and composition of the crystals were also analyzed by X-ray photoelectron spectroscopy (XPS, ESCALAB). Before electrical conductivity and photoconductivity measurements, the as-synthesized (Sn x Sb 1Àx ) 2 Se 3 crystals were sliced into small discs with a thickness of about 2 mm also with surface polishing. Electrical resistivity and bulk carrier density were measured using a van der Pauw Ecopia HMS-3000 Hall Measurement System (Ecopia Corporation, Anyang, South Korea). The temperature-dependent resistivity was measured by the conventional four-probe method in the temperature range of 300-400 K. To characterize the photoelectric performance, the photo-electro-chemical (PEC) measurement has been carried out using the classical 3-electrode method, where Ag/AgCl in saturated KCl was used as the reference electrode, Pt-wire as the counter electrode and our sample as the working electrode. Additionally, to prepare the working electrode, the as-synthesized bulk (Sn x Sb 1Àx ) 2 Se 3 crystals were sliced into small discs with a diameter of 10 mm and a thickness of about 2 mm, and then the surface was finely polished and cleaned. After this, silver paint was homogeneously applied on one face of the disc-like sample, acting as the back contact, which was connected by a copper wire and then covered with a layer of epoxy resin to be isolated from the electrolyte during the PEC test. Then these electrodes were immersed into the electrolyte (0.5 mol L À1 LiClO 4 solution), and a white light tungsten halogen lamp with controllable light intensities was used as the light source. A mechanical chopper was used to obtain chopped light and an Autolab Metrohm potentiostat was used to record the current-voltage characteristics.

Results and discussion

Fig. 1a shows the XRD patterns of the as-synthesized (Sn x Sb 1Àx ) 2 Se 3 polycrystalline semiconductors (x = 0.00, 0.03, 0.05, 0.07 and 0.10) acquired on sliced and polished disc-like samples, which all exhibit prominent peaks in agreement with the JCPDS standard card (15-0861) of the orthorhombic phase of Sb 2 Se 3 without a second phase, indicating the absence of any detectable impurity. The absence of multiple sets of peaks also indicates that the (Sn x Sb 1Àx ) 2 Se 3 crystals do not contain a second crystalline phase after Sn doping. Then the XRD patterns of the corresponding powdered samples obtained by grinding the bulk crystals into fine powders can be indexed to the same crystalline phase, as shown in Fig. 1b, suggesting that the crystallization process takes place uniformly inside the bulk sample and is not restricted to the surface. The relative peak intensity is different between the powder and the bulk because the preferential orientation in the bulk is destroyed during the grinding of the powders. It is noted that further were determined by refining the diffraction patterns through a Rietveld method implemented in the FULLPROF program, 19 which show a slight increase upon increasing the Sn doping concentration (Fig. 1c). This can be attributed to the small difference in the atomic radius between Sb (0.160 nm) and Sn (0.158 nm). And this slight difference satisfies the criterion for substitution solubility proposed by Hume-Rothery who pointed out that an extensive solid solubility of one metal in another only occurs if the diameter of the atoms differs by less than 15%. 20 Furthermore, the calculated unit cell volume of the (Sn x Sb 1Àx ) 2 Se 3 crystals gradually increases from 545.23 Å 3 (x = 0.00) to 546.53 Å 3 (x = 0.10), and this increase is known to occur due to doping. 21 Fig. 1d shows the XPS spectrum of Sn-3d, compared to the pure Sb 2 Se 3 . Two obvious peaks located at 494.7 and 486.2 eV from (Sn 0.10 Sb 0.90 ) 2 Se 3 crystals indicate the existence of Sn, further confirming the successful doping.

The (Sn x Sb 1Àx ) 2 Se 3 crystals' (x = 0.00, 0.03, 0.05, 0.07 and 0.10) morphology of the freshly fractured surface observed by SEM is shown in Fig. S3 (ESI †). All the samples consisted of large crystals with a layered-like structure and clear strips, indicating a preferential orientation. According to the theoretical and experimental studies, bulk Sb 2 Se 3 has an orthorhombic crystal structure and is composed of 1D ribbon-like (Sb 4 Se 6 ) n units stacking together via van der Waals interactions. 22 Its preferential growth is indeed determined intrinsically by the anisotropic Sb-Se atom chain or the layer structure. 7,10 To further investigate the doping nature and also the grain size of the (Sn x Sb 1Àx ) 2 Se 3 crystals, the finely polished disk-like samples were etched chemically by 1 mol L À1 NaOH solution for 20 min at room temperature. As shown in Fig. 2a, the Sb 2 Se 3 crystal consists of some big micron-sized grains. Fig. 2b shows a representative etched microstructure of the (Sn 0.05 Sb 0.95 ) 2 Se 3 crystal; the grains are clearly shown with an average size of approximately 15 mm. In addition, compositional spatial distribution was determined by EDS elemental mappings marked with Sn, Sb and Se, which reveals a uniform distribution of Sb and Se within the crystal. And a few Sn-rich zones between the grain boundaries can also be observed after chemical etching. Then for the (Sn 0.10 Sb 0.90 ) 2 Se 3 crystal with a higher Sn doping concentration, more obvious Sn-rich regions can be observed, as shown in Fig. 2c. This could be explained by phase separation, which is reasonable to occur, the minor Sn-rich phase solidifies first, and then is pushed to the grain boundaries upon freezing of the major phase. To further reveal the elemental composition of these Sn-rich regions, EDS analysis has been performed. Fig. 2d shows the EDS spectra recorded from two different positions (Point 1 in Fig. 2b and Point 2 in Fig. 2c); the atomic ratio of Sn : Sb : Se is 14.54 : 23.74 : 61.72 and 15.39 : 22.22 : 62.39 for Point 1 and Point 2, respectively. It is noted that the composition analysis of these Sn-rich regions performed on different (Sn x Sb 1Àx ) 2 Se 3 crystals show quite similar results, close to the ratio of 2 : 4 : 8 for Sn 2 Sb 4 Se 8 , which is also consistent with the XRD results. Fig. 3a shows the TEM image of an isolated particle milled from the (Sn 0.10 Sb 0.90 ) 2 Se 3 crystal. The selected area electron diffraction (SAED) pattern (Fig. 3b) from the red marked area (in Fig. 3a) is a periodic spot pattern, suggesting its highly crystalline nature. Fig. 3c shows the corresponding HRTEM image with obvious lattice fringes associated with the (020) planes (d-spacing of 0.591 nm) of the orthorhombic phase Sb 2 Se 3 . Then the EDS spectrum (Fig. 3d) exhibits strong Sb and Se peaks, and also an observable Sn peak located at 3.44 eV, indicating the successful Sn doping in this crystal.

Fig. 4a shows the electrical resistivity of the as-synthesized (Sn x Sb 1Àx ) 2 Se 3 crystals with different x values. It can be seen that the pure Sb 2 Se 3 exhibits extremely high resistivity. With Sn doping, a significant decrease of resistivity is observed, and the decrease is almost exponential with the x value. The temperature dependent electrical conductivity (s) in the temperature (T) range of 300-400 K has also been studied. As shown in Fig. 4b, the plots of ln s versus 1/T indicate that the electrical conductivity increases with increasing temperature for all samples, confirming semiconductor characteristics of these samples. Moreover, the conductivity exhibits single activation energy in the temperature range of measurement, according to the Arrhenius equation:

s = s 0 exp(ÀDE/kT) (1)
where s 0 is the pre-exponential factor, DE is the activation energy and k is the Boltzmann constant. 23 Thus the thermal activation energy DE can be calculated from the slope, while the intercept on the y-axis gives the s 0 value. Then the bulk carrier density and the doping type were determined by Hall measurements. All the electrical parameters of (Sn x Sb 1Àx ) 2 Se 3 crystals are summarized in Table 1. It is obvious that the room temperature electrical conductivity shows an increase with several orders of magnitude after Sn doping, accompanied by a synchronous increase of bulk carrier density, which change from 1.44 Â 10 11 cm À3 (x = 0.00) to 1.94 Â 10 16 cm À3 (x = 0.10). Importantly, such a great improvement will provide a bright future for widening applications of Sb 2 Se 3 . Since the activation energy varies inversely with the charge carrier density, Sn incorporation in this system should theoretically lead to a decrease of activation energy, which is consistent with the calculated DE values. The pre-exponential factor has a close relationship with the defect or impurities in the crystals, Sn doping into Sb 2 Se 3 also creates a compositional and configuration disorder, leading to an increase in the localized state density, which is an important factor for the change of s 0 values. Table 1 also shows the conduction type of (Sn x Sb 1Àx ) 2 Se 3 crystals; Sn incorporation has been found to result in P-type behavior, suggesting Sn substitution at the Sb site. Sn possesses one fewer valence electron than Sb, thus replacing Sb with Sn constitutes hole doping. 24 Photoconductivity is a well-known property of semiconductors which describes the electrical conductivity changes with the incident radiation. The photoconductive characteristics of the semiconducting (Sn x Sb 1Àx ) 2 Se 3 crystals are investigated here by using the classical 3-electrode photo-electro-chemical (PEC) measurements. A schematic illustration of this PEC measurement is shown in Fig. 5a. It is noted that with this 3electrode technique which measures the generated and separated electron-hole pairs on the surface, it can eliminate the influence of the bulk resistance of the sample, facilitating therefore the comparison between different samples. Fig. 5b shows the current density versus voltage under chopped light using (Sn x Sb 1Àx ) 2 Se 3 (x = 0.00, 0.03, 0.05, 0.07 and 0.10) as the working electrode, respectively. It is obvious that both dark current and photocurrent density of pure Sb 2 Se 3 are very weak, which can be explained by its intrinsic low conductivity and low carrier concentration. For Sn-doped crystals, a significant increase of current density is observed, and this increase is proportional to the Sn doping concentration. Up to x = 0.10, the dark current density increased approximately 10 times (from À16 mA cm À2 to À160 mA cm À2 at a bias voltage of À0.6 V) and the photocurrent density increased approximately 14 times (from À18 mA cm À2 to À250 mA cm À2 at À0.6 V). This improvement can be considered as an important step to promote Sb 2 Se 3 -based semiconductors for application in high performance photoelectric devices. Fig. 5b also shows that the Sn-doped Sb 2 Se 3 crystals are identified as p-type semiconductors, due to their photocurrent density increase with a negative bias, which is confirmed by the Hall measurement results. Then the light intensity dependence of current density is also investigated. As shown in Fig. 5c, using (Sn 0.10 Sb 0.90 ) 2 Se 3 as an example, upon illumination, the current density is increased, and then it increases gradually with increasing light intensity. Additionally, current-voltage curves of (Sn x Sb 1Àx ) 2 Se 3 at other x values (0.00, 0.03, 0.05 and 0.07) are shown in Fig. S4 (ESI †), where similar phenomena were observed. All the samples show that the photocurrent strongly depends on the light intensity, which can be attributed to the change in photo-generated carriers. The corresponding photocurrent density as a function of light intensity is shown in Fig. 5d, and dependence can be described by the power law:

J ph = AP y (2)
where J ph represents the photocurrent density; P is the power density of the incident light; A is a constant at a specific wavelength; and y is an exponent. 25,26 The fitting curves showed that y = 0.27, 0.52, 0.65, 0.60 and 0.54 for x values of 0.00, 0.03, 0.05, 0.07 and 0.10, respectively. All these y values are non-unity, indicating that some surface defects and carrier trapping at the localized states could result in a complex electron-hole non-radiative recombination process. 4 To evaluate the stability and the response speed of the photocurrent, chronoamperometry measurements have been performed. Fig. 6a andb show two plots of the time-resolved PEC response using Sb 2 Se 3 and (Sn 0.10 Sb 0.90 ) 2 Se 3 as the working electrode, respectively. As shown in Fig. 6a, when the light is on, the current density rapidly increases from À10.75 mA cm À2 (dark condition) to À12.25 mA cm À2 (with 22.51 mW cm À2 light illumination), and then sharply returns to its initial value as the light is turned off. This light response remains the same after multiple cycles, indicating the excellent reversibility and stability of this Sb 2 Se 3 PEC electrode. Then a similar photoresponse can be seen for (Sn 0.10 Sb 0.90 ) 2 Se 3 , except for a larger photocurrent density in the latter case (Fig. 6b). When using the (Sn 0.10 Sb 0.90 ) 2 Se 3 crystal as the working electrode, a huge photocurrent density spike before arriving at a steady value can also be observed, which was caused by recombination in the presence of surface states. 27 To extract the exact response time, one of the response periods is magnified, as shown in Fig. 6c andd, respectively. Under closer examination, the response time (t on ) and the recovery time (t off ) of the Sb 2 Se 3 working electrode are 0.05 s and 0.07 s, respectively. For (Sn 0.10 Sb 0.90 ) 2 Se 3 , both values show a decrease in value to 0.03 s, which can be attributed to much more surface dangling bonds, high crystal quality and preferential orientation. 28,29 It is worth noting that these short response times are comparable to the very good results of some reported high performance photodetectors, 30,31 further conforming the high potential application of Sn-doped Sb 2 Se 3 crystals as excellent photoelectric devices.

The optical absorption spectrum is one of the most effective tools for understanding and/or engineering the band structures of semiconductor materials. It is known that an appropriate band gap of the semiconductor plays a crucial role in fabricating high performance photoelectric devices. In this work, the optical absorption spectrum has been employed to characterize the as-synthesized (Sn x Sb 1Àx ) 2 Se 3 crystals. And the absorption data were calculated from diffuse reflectance data using Kubelka-Munk equations: F(R) = a/L = (1 À R) 2 /(2R), where R is the reflectance, a and L are the absorption and scattering coefficients, respectively. 32 As shown in Fig. 7a, the onset of the absorption of Sb 2 Se 3 crystals starts near 1200 nm, and a shift toward a longer wavelength is observed after Sn doping. The higher absorption efficiency for the doped crystals in the NIR range is interesting for broad spectrum applications. Fig. 7b shows a plot of [F(R)hn] 2 versus energy (hn); the direct band gap of Sb 2 Se 3 is 1.17 eV, and then for the Sn-doped Sb 2 Se 3 crystals, the band gap is similar and close to 1.10 eV. The dwindling band gaps of the doped (Sn x Sb 1Àx ) 2 Se 3 crystals can be explained by the cooperation of the Burstein-Moss shift and the renormalization effect. This mechanism exists in narrow band gap semiconductors with a specific doping level. 33,34 Finally, it is worth noting that the band gap of the (Sn x Sb 1Àx ) 2 Se 3 crystals is close to the best value for absorbing efficiently the solar spectrum, making them attractive candidates for photovoltaic application.

Conclusions

An effective melt-quenching method has been used to synthesize high-quality and variably doped (Sn x Sb 1Àx ) 2 Se 3 crystals with the dopant Sn ratio x ranging from 0.00 to 0.10. The successful doping is verified by XRD, XPS and EDS results. With increasing Sn doping concentration, (Sn x Sb 1Àx ) 2 Se 3 crystals exhibit an improvement of electrical conductivity with several orders of magnitude due to a much higher charge carrier concentration. Temperature dependent electrical conductivity displays a typical semiconductor characteristic and a decrease of activation energy is also confirmed with Sn incorporation. Sn-Doped Sb 2 Se 3 crystals are identified as p-type semiconductors according to the photoconductive characteristics. Notably, compared to pure Sb 2 Se 3 , the dark current density of a representative (Sn 0.10 Sb 0.90 ) 2 Se 3 increased approximately 10 times and the photocurrent density increased approximately 14 times. In addition, it shows a short response/recovery time (0.03/0.03 s) and a reversible and stable photoresponse. Therefore, the Sn-doped Sb 2 Se 3 crystals have an optimal band gap, (1.10 eV) which overlaps well with the solar spectrum, and excellent electrical conductivity and photoconductive properties, combined with their easy and relatively large-scale synthesis features. These properties confirm the high application potential of Sn-doped Sb 2 Se 3 as highly efficient photoelectric devices.

a

  State Key Laboratory of Silicon Materials, School of Materials Science and

A c c e

  photoconductivity is observed. Moreover, the as-synthesized bulk samples with easy and relatively large-scale synthesis features show very interesting properties for fabricating Sb 2 Se 3 -based photoelectric devices.

Fig. 1

 1 Fig. 1 (a) XRD patterns of the bulk (Sn x Sb 1Àx ) 2 Se 3 crystals (x = 0.00, 0.03, 0.05, 0.07 and 0.10). (b) XRD patterns of the powdered (Sn x Sb 1Àx ) 2 Se 3 samples. (c) Lattice constants a, b and c derived from powder XRD peaks, plotted as a function of Sn doping concentration (x value). (d) XPS spectrum of Sn-3d for the (Sn x Sb 1Àx ) 2 Se 3 crystals (x = 0.00 and 0.10).

Fig. 2

 2 Fig. 2 SEM images of (a) the Sb 2 Se 3 crystal, (b) the (Sn 0.05 Sb 0.95 ) 2 Se 3 crystal and (c) the (Sn 0.10 Sb 0.90 ) 2 Se 3 crystal surface with chemical etching by NaOH solution and following the corresponding EDS elemental mappings of Sn, Sb and Se. (d) EDS spectra recorded from two different positions (Point 1 in (b) and Point 2 in (c)).

Fig. 3

 3 Fig. 3 (a) TEM image of the representative (Sn 0.10 Sb 0.90 ) 2 Se 3 crystal. SAED pattern (b) and the corresponding HRTEM image (c) from the rectangle marked area in (a). (d) EDS spectrum.

Fig. 4

 4 Fig. 4 (a) Electrical resistivity of (Sn x Sb 1Àx ) 2 Se 3 crystals with different x values. (b) Temperature dependent electrical conductivity (s) of (Sn x Sb 1Àx ) 2 Se 3 crystals in the temperature (T) range 300-400 K.

Fig. 5

 5 Fig. 5 (a) Schematic diagram of the 3-electrode photo-electro-chemical (PEC) measurement. (b) Current-voltage characteristics using (Sn x Sb 1Àx ) 2 Se 3 crystals as working electrodes. (c) Dark current density and photocurrent densities at different power densities of the representative (Sn 0.10 Sb 0.90 ) 2 Se 3 crystal. (d) Photocurrent densities measured as a function of incident light density at a bias voltage of À0.7 V.

Fig. 6

 6 Fig. 6 Time-resolved PEC photoresponse by using Sb 2 Se 3 (a) and (Sn 0.10 Sb 0.90 ) 2 Se 3 (b) as the working electrode, respectively. (c) Magnified plot of one response cycle in (a), and (d) one response cycle in (b), used to obtain the response time and the recovery time.

Fig. 7

 7 Fig. 7 (a) Optical absorption spectra (calculated from diffuse reflectance data) of the (Sn x Sb 1Àx ) 2 Se 3 crystals (x = 0.00, 0.03, 0.05, 0.07 and 0.10) (b) a plot of [F(R)hn] 2 vs. energy, from which direct band gap was obtained.

Table 1

 1 A summary of electrical conductivity (s), bulk carrier density (n), pre-exponential factor (s 0 ), activation energy (DE) and conduction type of (Sn x Sb 1Àx ) 2 Se 3 crystals with varied xx value s (S cm À1 ) n (cm À3 ) s 0 (S cm À1 ) DE (eV) Conduction type

	0.00 0.03 0.05 0.07 0.10	2.71 Â 10 À6 1.44 Â 10 11 1.31 Â 10 3 0.52 1.30 Â 10 À3 1.52 Â 10 15 2.65 Â 10 3 0.38 4.12 Â 10 À3 3.59 Â 10 15 9.96 Â 10 2 0.32 9.65 Â 10 À3 6.90 Â 10 15 1.49 Â 10 3 0.30 7.50 Â 10 À2 1.94 Â 10 16 3.48 Â 10 2 0.22	N P P P P

Acknowledgements

The authors gratefully acknowledge support for this research from the Program for Innovative Research Team in the University of Ministry of Education of China (IRT13R54), the Fundamental Research Funds for the Central Universities (2016FZA4007) and the Zhejiang Provincial Natural Science Foundation of China (LY16E020003). We also thank the Chinese Scholarship Council and the Eiffel Excellence Scholarship Programme for financial support.

Notes and references

Conflicts of interest

There are no conflicts to declare.


A Bel-Brunon

S Kehl

C Martin

S Uhlig

W A Wall
email: wall@lnm.mw.tum.de

Numerical identification method for the non-linear viscoelastic compressible behaviour of soft tissue using uniaxial tensile tests and image registration -Application to rat lung parenchyma

Keywords: lung parenchyma, non-linear viscoelasticity, numerical identification, image registration, uniaxial tensile test, compressibility

Numerical identification method for the non-linear viscoelastic compressible behavior of soft tissue using uniaxial tensile tests and image registration -Application to rat lung parenchyma

Introduction

Characterizing the mechanical behavior of lung tissues has been the subject of several scientific studies as it is required to model the global and local behavior of the lung. A computational model of the whole lung could help understanding several phenomena related to lung diseases. For instance, acute respiratory distress syndrome (ARDS) is a disease that alters the mechanical behavior of the lung parenchyma [START_REF] Kallet | Respiratory system mechanics in acute respiratory distress syndrome[END_REF]. When a patient suffering from ARDS is mechanically ventilated to help him breathe, healthy and diseased areas are deformed but the compliance difference of these two tissues can lead to a local overstretching in the parenchyma. A lung model including mechanical properties of both healthy and diseased parenchyma could help locating and quantifying the possible overstretching occurring within the parenchyma and adapting the curing protocol to avoid it. A comprehensive model of the lung is currently developed in our group including among others, structure and fluid phenomena [START_REF] Wall | Towards a comprehensive computational model for the respiratory system[END_REF].

Lung parenchyma consists of alveolar tissue, including respiratory bronchioles, alveolar ducts and terminal bronchioles. At the macroscale (a few cm), it can be considered as homogeneous but at the microscale (alveolar scale, i.e. a mm), it exhibits a foam-like structure due to the alveoli, which are inflated with air during breathing. Therefore, unlike most of the soft tissues commonly studied (e.g. abdominal organs, skin, aorta), it cannot be considered as incompressible at the macroscale as its main function is to be inflated. Most of the previous studies to characterize the lung parenchyma have investigated its non-linear elasticity and viscoelasticity but the compressible properties have not been characterized sufficently so far.

To characterize lung parenchyma, the few available previous studies describe mostly uniaxial tensile tests on tissue slices. Briefly, in terms of elasticity, either an exponential potential [START_REF] Navajas | Dynamic viscoelastic nonlinearity of lung parenchymal tissue[END_REF] or a polynomial potential [START_REF] Rausch | Material model of lung parenchyma based on living precision-cut lung slice testing[END_REF] describe rather accurately the uniaxial behavior of the lung parenchyma. Gao et al. performed biaxial tensile tests and concluded that an exponential law was the best fit for the tissue behavior [START_REF] Gao | Mechanical properties of human lung parenchyma[END_REF]. A volumetric component of the constitutive law is only assessed in [START_REF] Rausch | Material model of lung parenchyma based on living precision-cut lung slice testing[END_REF]; however the method to identify the volumetric parameters of the constitutive law only relied on the displacement of a single location in the sample and the measurement of this displacement was done by hand.

In terms of viscoelasticity, uniaxial tensile tests were used to characterize the dynamic response of the tissue. Different frameworks helped quantifying the tissue viscosity, the most common [START_REF] Mijailovich | Dynamic moduli of rabbit lung tissue and pigeon ligamentum propatagiale undergoing uniaxial cyclic loading[END_REF][START_REF] Navajas | Dynamic viscoelastic nonlinearity of lung parenchymal tissue[END_REF][START_REF] Romero | Dynamic nonlinearity of lung tissue: frequency dependence and harmonic distortion[END_REF][START_REF] Pinart | Dynamic nonlinearity of lung tissue: effects of strain amplitude and stress level[END_REF][START_REF] Yuan | Dynamic properties of lung parenchyma: mechanical contributions of fiber network and interstitial cells[END_REF] being quasilinear viscoelasticity proposed by Fung that decouples the elastic and viscous phenomena and allows accounting for elastic non-linearities [START_REF] Fung | Biomechanics: Mechanical Properties of Living Tissues[END_REF]. In these studies, harmonic distorsion, loss and storage moduli as well as elastance, resistance and hysteresivity of the tissue were computed with the help of a Fourier transform. These parameters allow the investigation of the influence of loading factors -like stress amplitude, strain amplitude, frequency -or pathologies -like fibrosis (Dolhnikoff et al., 1999), acute lung injury [START_REF] Ingenito | Effects of acute lung injury on dynamic tissue properties[END_REF][START_REF] Rocco | Lung tissue mechanics and extracellular matrix remodeling in acute lung injury[END_REF] -on the material parameters. A review of the elastic and viscous characterization of lung tissue, in particular using tests on excised samples, is given in [START_REF] Suki | Lung tissue mechanics as an emergent phenomenon[END_REF].

Another framework is proposed by [START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF] to account for the three-dimensional non-linear viscoelasticity of fiber-reinforced composites: the hyperelastic potential associated with the elastic behavior of the tissue is modified by a viscous contribution. This framework is particularly suited for a direct use in nonlinear continuum mechanics based simulations, this is why this approach is chosen in the present paper.

The aim of this paper is to determine a non-linear viscoelastic constitutive employing novel experimental and identification methods. The paper is divided into five parts: the first part covers the description of the experimental methods, particularly the optical method to measure the displacement field on the sample. In the second part, the numerical model and the constitutive law proposed for the lung parenchyma are detailed. The third, fourth and fifth parts are focused on the identification method, the results and the discussion, respectively.

Uniaxial tensile tests on living PCLS

Since the preparation of the samples and the experimental set-up are described in a previous paper [START_REF] Rausch | Material model of lung parenchyma based on living precision-cut lung slice testing[END_REF], only the main features are recalled here.

Preparation of the samples

The Precision-Cut Lung Slices (PCLS) are prepared from three isolated rat lungs as previously described in [START_REF] Martin | Videomicroscopy of methacholineinduced contraction of individual airways in precision-cut lung slices[END_REF]. The lungs are dissected from the animal, filled with an agarose solution (1.5%) via the trachea and put on ice to allow the agarose to cool down and solidify. The lung lobes are separated and cut into tissue cylinders using a coring tool.

The obtained cylinders have a diameter of 14 mm; they are then cut into slices using a Krumdieck tissue slicer (Alabama Research and Development, Munford, AL). The thickness is set to 400 µm. The two sides of the round slices are trimmed with two parallel razor blades giving the tissue strip a width of 7.0 mm. The strips are incubated in minimal essential medium (MEM) and the agarose is washed out by frequently changing the medium within the first 4 hours and completing an overnight incubation. As shown in [START_REF] Martin | Videomicroscopy of methacholineinduced contraction of individual airways in precision-cut lung slices[END_REF], the obtained PCLS are viable for three days. In this protocol, the rat samples are tested within 48 hours after death, which guarantees that the tissue is still living while being tested [START_REF] Martin | Videomicroscopy of methacholineinduced contraction of individual airways in precision-cut lung slices[END_REF] ; preliminary studies showed that the mechanical properties are not changing within 3 days (data not shown). The samples are kept in saline solution in the fridge (4 • C) until they are tested. Tab. 1 sums up the specimen numbers.

Experimental set up

The samples are tested using a Bose ElectroForce 3100 uniaxial device (Bose Cooperation, Eden Prairie, USA). The samples are fixed onto the machine using specifically designed clamps to prevent sliding and damage (Fig. 1). The moving clamp displacement is recorded using the internal displacement sensor of the tensile machine (range of ± 2.5 mm with a resolution of 12.5 µm); the other clamp is connected to the machine through a force sensor (range of ± 0.5 N with a resolution of 2.5 mN). Finally, a digital microscope is placed perpendicular to the sample plane to record the test. the gray-levels are affected by a modification of the lightening or by light reflections at the surface of the sample, artifacts will occur in the computation of the displacement field.

In our case, the registration procedure relies on the microstructure of the lung parenchyma. By using a lightening that goes partially through the sample, this microstructure can be visualized; hence, the application of an additional random pattern using paint -as it is usually necessary when using an optical method to compute the experimental displacement -is superfluous.

However, due to the low density of the tissue, the light rays intensity is hardly lessened when going through the tissue. Putting the light source and the microscope directly opposite to each other on both sides of the sample leads to a general overexposure, as well as a modification in the overall light intensity in the image during the sample stretching (since the sample gets thinner while being stretched, the penetration of light gets stronger). The solution to have a good insight of the microstructure and a satisfactory light intensity for the images was to put the light source on one side and slightly above the sample (Fig. 1); therefore the light is slightly scattered which reduces reflecting signals. The drawback of this method is that it does not completely prevent the occurrence of shining points, where the light is reflected on the wet surface of the sample. Even if not numerous, these parasite pixels can create artifacts during the registration procedure. However, the image registration algorithm helps smoothing these possible artifacts.

Testing protocol

Thickness measurement and placement on the machine

The first step of the experimental protocol is to measure the actual thickness of the sample. Although the cutting tool is very accurate, the softness of the lung tissue and its compressibility inevitably lead to a variability in the thickness. To accurately measure the global sample thickness, a light microscope and a micrometer are used as shown in Fig. 2. Two microscopic slides are placed on each other and under the microscope; a pattern is drawn on the upper one. The microscope is focused on the pattern and the vertical position of the microscope plate is measured by a micrometer. The sample is then placed between the two microscopic slides and the water is wiped up so as not to create a water film between the two microscopic slides. At this stage the upper slide is in a higher position due to the thickness of the sample and the microscope does not focus on the pattern anymore. Then, the focus is adjusted on the upper slide pattern again; the resulting displacement of the microscope plate is the thickness of the sample, measured by the micrometer.

The thickness is measured within an accuracy of 10 µm.

Once the thickness has been measured, the sample is hydrated again The initial length of the sample, i.e. the initial distance between the clamps, is set to 2.0 mm.

Preconditioning and tensile test

In this work, the aim is to investigate the tissue behavior in the range of large strains. Preliminary tensile tests conducted on lung samples showed that when stretched to 50% strain and released, the sample exhibits a residual strain. Besides, during cycling tensile loading, the viscoelastic nature of the tissue leads to a transitory phase where the force-displacement cycles cannot be superimposed. One can assume that during breathing, the tissue exhibits a rather repeatable behavior from one cycle to the other. Therefore, a preconditioning phase is applied to the tissue to first, stabilize its residual strain and second, to stabilize the hysteresis loop obtained while cycling the sample.

In the present case, the preconditioning includes two main steps: first, the sample is stretched from 0 to 0.2 mm at a constant velocity of 0.5 mm/s and maintained 10 s at 0.2 mm; then this phase is repeated for 0.4, 0.6 and 0.8 mm. The tissue is then released to a stretch of 0.2 mm for 30 s. Second, the tissue is cycled between 0.2 and 1.2 mm for 10 cycles. The nine first cycles are used for preconditioning while the tenth is the actual tensile test.

We found that nine cycles were enough to stabilize the cycle (see Sec.5.1).

Reference state

The preconditioning is performed to stabilize the residual strain; it leads to the definition of the stress-free reference state which is different from the initial state. We choose to perform the tensile test in a range that makes the sample slightly buckle at the minimum stretch. Therefore, when cycling, the force is slightly negative (compressive force) and becomes positive. The state corresponding to a force equal to zero is defined as the reference state for the image registration procedure. In practice, the zero-force state does not necessarily match one of the registration images for sampling frequency reasons; in that case, we choose the image that corresponds to the smallest force to define the reference geometry of the sample. Over all the samples, the initial force is always below 5% of the maximal force of the test. We also assume that the difference between reference and initial configurations is small enough to consider that the initial thickness is still applicable to the reference configuration. Indeed, the reference state corresponds experimentally to a clamp displacement of around 0.1 mm.

Postprocessing by image registration

Using a full-field measurement method when testing soft material has several advantages. First, cutting samples of soft tissues of accurate dimensions is a difficult task; even with precision-slicing apparatus and procedure as described in Section 2.1, one cannot always obtain perfectly rectangular samples. Then, such thin samples are difficult to handle and to place in the clamps. The method described previously eases this procedure but we can for instance not perfectly control the orientation of the sample with respect to the loading direction. Finally, screwing the clamps to properly prevent sliding without damaging the sample requires some care; it can happen that the sample is not squeezed enough and experiences sliding during the test.

For these reasons, it makes sense to use optical measurements additional to the force and global displacement information to know more accurately the actual displacement of the sample during the experiment.

Besides, the present work is focused on a compressible tissue. The identification procedure therefore requires information about the evolution of the width of the sample, which is directly related to the compressibility of the tissue. Global longitudinal measurement would not be sufficient to accurately identify the compressibility material constant.

As mentioned before, image registration consists in mapping two images of the sample that represent two different deformation states. These two images are called source image and target image; they both contain a scalar field that is a gray-level distribution over the pixels of the image. Image registration consists in finding the displacement field that minimizes the gap between the target image and the source image once it is deformed by this displacement field. We describe here the main features of this minimization procedure; more detailed information can be found in the cited papers.

The source and target images are denoted by I s (X) and I t (x) respectively where x and X are the spatial and material coordinates respectively.

For a so-called forward transformation, registering the two images consists in minimizing the error D with respect to the transformation ϕ [START_REF] Haber | A framework for imagebased constrained registration with an application to local rigidity[END_REF]:

D [ϕ(X)] = Ωx I s ϕ -1 (x) -I t (x) 2 dx (1)
where ϕ : I s → I t (Fig. 3) so that x = ϕ(X) = X + u with u a spatially varying displacement field. The solution u * of the minimization is written: u * = arg min

φ I s (X) I t (x) Ω 2mm
u∈L 2 D (2)
An elastic regularization R is added to ensure the well-posedness of the optimization problem [START_REF] Modersitzki | Numerical Methods for Image Registration[END_REF][START_REF] Fischer | Ill-posed medicine -an introduction to image registration[END_REF].

The minimization is therefore:

u * = arg min u∈L 2 [D [I s , I t ] + αR(u)] (3)
where the weighting factor α is set to 0.1 in this work.

We minimize D using a Gauss-Newton optimization scheme [START_REF] Haber | A multilevel method for image registration[END_REF] as well as a multiresolution approach (Lester and Arridge, 1999) (Fig. 4). First, from the pixel-size grid, i.e. the image recorded during the experiment, several coarser and coarser grids are successively created by filtering and sampling the gray-level field of the image at the one-level finer scale. We obtain n levels of resolution. The minimization starts at the coarser level n; the result u n of the minimization of D is linearly interpolated over the grid at level n -1; u n-1 n is used as the initial guess for the minimization at level n -1. This procedure is repeated until u 1 = u * -at the pixel levelis determined. In our case, five levels of resolution are used. The image registration algorithm is implemented in Matlab, as well as the post-processing of the registration results (creation of the input file for the simulation).

The overall registration procedure can be computed in two different ways:

sequentially or incrementally. The sequential method consist in defining a reference image (corresponding to the reference state of the sample) and register all the other images (while the sample is stretched) with respect to this reference image. The incremental method consists in updating the reference image: image n + 1 is registered with respect to image n. In this case, the Finite Element mesh (see next section) needs to be updated at each new reference image: the position of the nodes of the FE-mesh at step n is updated using the result of the registration between images n -1 and n. Sec.

5.2 gives a comparison of the two methods.

The images provided by the microscope are 640 × 480 px 2 , corresponding to 15.2 × 11.4 mm 2 . The registration procedure guarantees an accuracy of a pixel size of 0.023 mm.

Numerical model of the lung parenchyma samples

Finite Element model

The numerical identification procedure is based on the simulation of the experimental test and the minimization of the error between the simulated and the experimental displacements. In order to simulate the test, a zone of interest (ZOI) is defined on the reference image and a Finite Element mesh is created on this ZOI. The shape of the ZOI is quadrangular but not necessarily rectangular. The registration procedure is computed on the entire image but only the results in the ZOI are further used for the identification. The Matlab function interp2 is used to interpolate the computed displacements on the nodes of the FE-mesh. The thickness that was experimentally measured is defined as the initial thickness of the simulated sample.

The simulation is performed using our in-house multi-purpose research software platform BACI. The input data for the simulation of the test are the FE-mesh, the material law with the initial guess of parameters and the boundary conditions.

Constitutive law

The aim of this study is to identify a suitable visco-hyperelastic constitutive model for lung parenchyma and determine the corresponding model parameters. The model that we use in the present case is a generalized

Maxwell model with different hyperelastic contributions.

The hyperelastic contribution is described in details in [START_REF] Rausch | Material model of lung parenchyma based on living precision-cut lung slice testing[END_REF]. Briefly, a hyperelastic toolbox is implemented in BACI which allows the combination of several different Strain Energy Functions (SEF). We only mention here the isotropic SEF as the lung parenchyma is isotropic. The SEF are of two types: coupled and decoupled. A coupled SEF is based on the overall strain tensor, while a decoupled SEF is split into volumetric and isochoric contributions. The decoupled contributions are all polyconvex and stress-free in the reference state, which allows to sum them up without further condition.

In our case, only decoupled potentials are used due to the higher variety of available potentials of this form than for coupled potentials.

The viscous contribution is based on the work of [START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF]. A viscous part is added to the hyperelastic contribution based on a generalized Maxwell model. In this paper, we consider the basic form of the generalized Maxwell model, i.e. a dashpot in parallel with a spring. The general form of the potential is the following:

ψ = ψ ∞ (I i ) + γ (C, Γ) (4)
where ψ ∞ is the hyperelastic potential and γ is the dissipative potential.

This theoretical framework is detailed for the isotropic case in Appendix A; details for anisotropic cases can be found in [START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF].

Similar to the right Cauchy-Green strain tensor C being the kinematic conjugate of the elastic second Piola-Kirchhoff stress tensor S, Γ is the kinematic conjugate of the fictitious stress tensor Q. These stress tensors are expressed by:

S ∞ = 2 ∂ψ ∞ ∂C (5) Q = -2 ∂γ ∂Γ (6)
Q is a non-equilibrium stress; therefore, Q = 0 when t → ∞. The evolution of Q is governed by an equation of evolution:

Q + 1 τ Q = β Ṡ (7) Q| t=0 = 0 (8)
where τ defines a relaxation time and β is a non-dimensional free-energy factor. The global stress tensor is therefore S = S ∞ + Q. The elasticity tensor is also expressed as the sum of an elastic summand and a viscous summand (Eq. 9).

C = C ∞ + C vis (9)
Appendix A and Appendix B sum up the formulas used to compute stress and elasticity tensors for the isotropic case and with the assumption that the viscous effect is the same for isochoric and volumetric contributions, that is to say that a single (τ , β) couple is used.

The hyperelastic potentials tested in the present work are summed up in Tab. 2. Only isochoric potentials involving the first invariant are identified as it was shown in [START_REF] Rausch | Material model of lung parenchyma based on living precision-cut lung slice testing[END_REF] that potentials with both the first and the second invariants did not lead to a significant decrease of the error in the identification. We compared several combinations of hyperelastic potentials to find out which one gives the best fit. Tab. 3 displays the tested combinations.

Identification method

The identification procedure consists in minimizing the error between computed and experimental displacements for the last cycle, i.e. both the loading and unloading curves are taken into account. The input data for this inverse analysis are the simulation data mentioned before and the experimental displacements at the nodes of the FE-mesh that are taken into account for the cost function calculation -typically all the degrees of freedom of the FE-mesh nodes that are not kinematically constrained (Fig. 5).

Vol.

Sussman-Bathe

ψ vsb = κ 2 (J -1) 2 Ogden ψ vo = κ 4 (-2 ln J + J 2 -1) Penalty ψ vp = ǫ J γ + 1 J γ -2 Iso.
Polynomial order n

ψ ipn = c n ( Ī1 -3) n Exponential ψ ie = k 1 2k 2 e k 2 Ī1
3 -1 -1

Table 2: Volumetric and isochoric hyperelastic potentials tested in this work. Parameters in red are to be identified. Sussman-Bathe from [START_REF] Doll | On the development of volumetric strain energy functions[END_REF]. Ogden from [START_REF] Ogden | On the overall moduli of non-linear elastic composite materials[END_REF]. Penalty from (Balzani et al., 2006). Polynomial from a split and extension of Yeoh's potential [START_REF] Yeoh | Some forms of the strain energy function for rubber[END_REF]. Exponential from (Delfino et al., 1997). The minimization follows a pseudo-transient continuation type scheme.

A ψ = ψ vo + ψ ip1 + ψ ip2 + ψ ip3 B ψ = ψ vo + ψ ip1 + ψ ip3 C ψ = ψ vo + ψ ip1 + ψ ip2 K ψ = ψ vo + ψ ip1 + ψ ip4 P ψ = ψ vo + ψ ie
The constitutive parameters of each summand of the material law and the associated viscous contribution are stored in the parameter vector p; its size is n p . For each loop of the inverse analysis, n p + 1 simulations are performed: n p simulations with only one modified element of p (p(i) = p(i)+ a + b × p(i), where a and b are user-defined perturbation parameters) and one simulation with unperturbed p. These simulations allow the determination of the Jacobian matrix J of the problem (Eq. 10) and the computation of the new set of parameters p = p + ∆p (Eq. 11,12)

J =      ∂u 1 ∂p 1 • • • ∂u 1 ∂pn p . . . . . . ∂un ∂p 1 • • • ∂un ∂pn p      (10) ∆p = J T J + µdiag J T J -1 J T R (11) R =      r 1 . . . r n      =      u 1,sim -u 1,exp . . . u n,sim -u n,exp      (12)
where u i is the displacement at node i. The criterion of convergence of the minimization is based on the gradient of the error ∇(e) where e is the cost function, i.e. the sum of the error at all nodes and all time steps.

The parameters of the inverse analysis are therefore the convergence criterion, the maximum number of optimization loops, the perturbation parameters a and b and an initial regularization factor µ.

A nested parallelization implemented in BACI allows computing these simulations in parallel, each simulation being itself also computed in parallel; regardless of n p , one loop of the inverse analysis lasts as long as the longest of the n p + 1 simulations, instead of n p + 1 times the time of one simulation.

For each combination of one volumetric potential and possibly several isochoric potentials of Tab. 2, a purely hyperelastic identification is done first, then the viscous contribution is added, which means that an additional parameter τ is to be identified -in this work we set β to 1 and do not identify it. The results of the first identification are used as initial guesses for the viscoelastic identification. During this second phase, both the hyperelastic and the viscous parameters are optimized.

Results

The identification results are statistically analyzed using the non-parametric Mann-Whitney test [START_REF] Mann | On a test of whether one of two random variables is stochastically larger than the other[END_REF]; the p-factor is indicated when a tendency is significant.

Force versus displacement curves

Fig. 6 shows the evolution of the hysteresis loop during the ten cycles of the tensile loading. The hysteresis loop decreases and stabilizes; only the last cycle is taken into account for the identification.

Comparison between sequential and incremental image registration

The displacements obtained with the incremental and sequential methods are compared in Fig. 7 and8. Fig. 7 shows that when the stretch increases (from image 4 to image 14 in that case), the gap between the two methods increases but remains small (less than 2% of the maximum displacement in x-direction). Fig. 8 shows that the incremental method gives a slightly better result as the new position of the FE-mesh nodes corresponds better to the visually expected one (the displacement of the moving clamp nodes is better). That is why we choose the incremental method.

Convergence of the identification procedure

Parametric studies have been conducted to ensure that the present identification method is robust and valid. In the following, the criterion of convergence is ∇(e) ≤ 10 -6 .

Convergence using artificial data. For this study, artificial data are generated using one of the experiments. A test is simulated with experimental load conditions and an arbitrary set p of material parameters. The computed displacements are used to create a new input file. We then proceed to identify the material parameters of the same law using this file from an arbitrary initial guess, and check that the identified parameters coincide with p. This procedure has been conducted using a hyperelastic potential of the following form :

ψ = ψ vo + ψ ip1 + ψ ip3 (13)
For a total number of elements from 250 to 4410, the maximal value of the average error on the three identified parameters is 0.00731 % (250 elements) which shows convergence to the right values (target values: κ = 30 kPa, C 1 = 5 kPa, C 3 = 10 kPa; initial guess: κ = 50 kPa, C 1 = 0.5 kPa, C 3 = 5 kPa).

These values were obtained after 9 to 11 identification loops.

To study the convergence of the viscoelastic identification, one viscous constant τ is added. In that case, four parameters are to be identified (target values: κ = 70 kPa, C 1 = 3.5 kPa, C 3 = 33 kPa, τ = 0.1; initial guesses: κ = 30 kPa, C 1 = 5 kPa, C 3 = 50 kPa, τ = 0.05). A single FE-mesh size of 1210 elements is used. The identification leads to an error on the target value in the order of 1 % for τ and 0.05 % for the other parameters. The convergence is therefore very satisfactory.

Influence of the FE-mesh size. Fig. 9 shows the result of the identification of the law of Eq. ( 13) for different mesh sizes. All the identifications started from the same initial guess; they are conducted on a real test (not on generated data from known parameters as for the previous study). When the number of elements increases, the identified parameters converge to their optimized value, as the evolution of the error shows. The compromise between computational time and accuracy was set to around 1200 elements. Sensitivity analysis. We evaluated the sensitivity of the error to small variations of the constitutive parameters. Fig. 10 and 11 display the evolution of the error for the hyperelastic law given in Eq. ( 13) with a variation of ±10% around the identified value of the parameters κ, C 1 and C 3 , and for a viscoelastic law (potential of Eq. ( 13) + viscous contribution) with a variation of ±20% around the identified value of the time parameter τ . In both cases, when we look at the scales, the error map is rather flat, but it is however slightly influenced by small variations of the parameters. In particular, the first order parameter (C 1 ) has a strong influence, while C 3 and κ are rather coupled (the error displays a valley). However, we could check that a minimum exists by testing different initial guesses for all the samples and by checking that the identification always converges to the same set of parameters. For the viscoelastic identification, there is a clear minimum of the error with respect to τ .

Hyperelastic identification

For both the hyperelastic and viscoelastic identifications, information on the thickness is required. The measured thickness for all the tested samples was 442 ± 34 µm, 503 ± 57 µm and 519 ± 70 µm for rats G, I and J respectively. For different reasons (rupture, poor quality of the image registration) some of these tests could not be fully post-processed.

Volumetric potential. As mentioned in Sec.3.2, three volumetric potentials are tested. Fig. 12 displays the error obtained with these three volumetric potentials combined with the same isochoric potential of type B (which parameters are not identified), for all the samples of a single rat. We can see that Sussman-Bathes and Ogden potentials give equal results while the penalty potential is either equivalent or worse. Therefore, the volumetric potential used in the following parts will always be the Ogden form.

Identification of volumetric and isochoric parameters. Fig. 13 shows the results obtained for each tested frequency among all the tested samples (three rats). We can see the comparison between the combinations of potentials given in Tab. 3. This shows that the performances of these forms of potential give rather similar results, but the exponential form (law P) gives the best fit for each frequency. A Mann-Whitney statistical test computed on all the samples shows that law C gives significantly poorer results than laws A, B and P.

Tab. 4 gives the value of the identified parameters for these hyperelastic potentials. As there is no statistical difference between the identified parameters for the varying tested frequencies, the values are averaged over all the samples and all the frequencies. The identified value of κ is averaged over all the potentials as this value is independent from the form of the isochoric potential. Top: tested frequencies.

Viscoelastic identification

Fig. 16 displays the identification results for the viscoelastic laws for the tested frequencies. Altogether, the performance of law F is significantly lower than all the other laws (the p-factor is ≤ 0.05 with respect to law D, ≤ 0.001

Law Parameter Value Unit κ 7.06e4 ± 2.93e4 Pa D C 1 3.22e3 ± 9.53e2 Pa C 2 2.04e3 ± 2.86e3 Pa C 3 2.92e4 ± 1.52e4 Pa E C 1 3.54e3 ± 9.48e2 Pa C 3 3.30e4 ± 1.58e4 Pa L C 1 4.04e3 ± 1.07e3 Pa C 4 9.62e4 ± 5.09e4 Pa M k 1 5.58e3 ± 1.

Discussion

This paper presents a method to determine material parameters for the lung parenchyma using advanced methods from both the experimental and the theoretical point of view. Thanks to this method, we obtain a reliable viscoelastic constitutive law that is suitable for Finite Element calculations.

The experimental protocol has been improved since [START_REF] Rausch | Material model of lung parenchyma based on living precision-cut lung slice testing[END_REF].

Measuring the thickness of each sample individually revealed its non-negligible variability throughout the set of samples; this distribution is now taken into account in the simulation. The other main improvement is the use of an optical method to determine the displacement field at the sample surface. It is worth noticing that as the samples are really thin, we can consider that this surface information is sufficient to describe the behavior of the tissue. In our case, knowing the displacement field is of high importance as we are particularly interested in the lateral contraction of the sample. When we optically check the results of the image registration by superimposing the image and the corresponding mesh (Fig. 8 for instance), we can see that this optical method applied to the natural pattern of the parenchyma gives fully satisfactory results, despite some light reflections at the sample surface. These artifacts are indeed smoothed by the Finite-Element kinematics inherent to the image registration method.

The theoretical constitutive framework used in this paper is adapted from an existing model proposed in [START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF]. It allows the description of the non-linear elasticity as well as the viscosity for both the isochoric and the volumetric components of the deformation, not only the isochoric one.

In the present work, a quantitative determination of elastic as well as viscous properties of rat lung parenchyma is provided. The order of magnitude of the strain in our protocol is 0 to 45 %. In this range, the non-linear elastic behaviour of lung parenchyma is equally well described by a polynomial potential of degree ≥ 3 associated to a linear component (degree 1) or an exponential potential. No significant influence of the testing frequency was found in the identified elastic parameters. For the volumetric aspect of the law, we could determine that both Ogden and Sussman-Bathe give the same result and that the volumetric parameter κ is rather accurately identified (c.f. low standard deviation). In [START_REF] Rausch | Material model of lung parenchyma based on living precision-cut lung slice testing[END_REF], values for the parameters for a law of type B are given. C 1 is similar to the one of this paper, while κ and C 3 are smaller. The explanation can come from a better description of the lateral contraction in the present paper, that highly influences the determination of κ, as well as the consideration of the actual thickness of the sample.

Adding a viscous component to the identified non-linear elastic constitutive law significantly increases the quality of the identification, especially at low frequency where the error is cut by half compared to the purely elastic case. It allows the description of the experimentally observed hysteresis loop. The elastic contribution of the law is rather insensitive to the frequency, but the time constant τ significantly decreases with increasing frequency. It

shows that the lower the frequency, the higher the viscous effects.

The main limitation of this work is that the uniaxial tension loading mode does not correspond to the mainly volumetric deformations the tissue would be subjected to in reality. It is not guaranteed that the constitutive law determined in this paper for an isochoric-dominated loading would not be influenced by a volumetric-dominated loading. Identifying the volumetric contribution using the width evolution is however a first step towards an accurate estimation of the compressible behavior of the homogeneous lung parenchyma.

A possible improvement for the presented method would be to replace the two optimization processes (image registration and inverse analysis) by a single one. It would consist in improving the mechanical regularization in the image registration, currently purely elastic, by taking into account the actual material law that is to be identified. This type of mechanical/optical identification method is described for instance in [START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF] and could be implemented in our in-house code in a straightforward way.

The present experimental protocol and identification method are now to be applied for more biological tissues. Although they are particularly adapted to compressible material, they could also be used for any soft biological tissue, possibly requiring the addition of a random pattern on the surface of the sample. Future work is to investigate the viscoelastic properties of human lung parenchyma as well as parenchyma samples altered with ARDS. These data could then be included into the computational model of the respiratory system that is developed in our group.

Figure 1 :

 1 Figure 1: Experimental setup. The double arrow shows the tensile direction.

Figure 2 :

 2 Figure 2: Method to measure the thickness of the sample: first, we focus the microscope on the upper slide without sample. Second, once the sample is put between the two slides, the microscope is focused again on the upper slide and the displacement of the microscope plate is measured by a micrometer. The measured displacement is the thickness of the sample.

Figure 3 :

 3 Figure 3: Source image (left) and target image (right).

Figure 4 :

 4 Figure 4: Multiresolution approach to solve the minimization problem of the image registration. u i is the solution of the minimization at level i. u j i is the linear interpolation of u i over level j and is the initial guess for the resolution at level j. u 1 = u * is the solution of the problem. The finer grid (lower left) is the image recorded during the test.

Figure 5 :

 5 Figure 5: Degrees of freedom used for the inverse analysis. The nodes at the fixed clamp are not used. Only the x-displacement of the moving clamp nodes is used; for the other nodes, the two degrees of freedom are used.

Figure 6 :

 6 Figure 6: Example of measured curve over ten cycles of tension loading (sample from rat G, frequency: 0.4 Hz). The displacement displayed here is measured by the sensor of the tensile machine.

Figure 7 :

 7 Figure 7: Difference of the displacements computed with the incremental and sequential image registration. At each image, average values and maximum values (over all the FEmesh nodes) of the gap in the x-and y-directions are normalized by the maximum value of the x-and y-displacements respectively.

Figure 8 :

 8 Figure 8: Position of the FE-mesh nodes superimposed on the test picture: position computed by the incremental (blue) and sequential (yellow) image registration methods, for a deformed image. The moving clamp is the right one.

Figure 9 :

 9 Figure 9: Evolution of the identified parameters for law B (see Tab. 3) with respect to the total number of elements, for one sample of rat G. The evolution of the error e is displayed in red.

Figure 11 :

 11 Figure 11: Sensitivity of the error to small variations of parameter τ of a viscoelastic law. The sensitivity study is conducted on a sample of rat I tested at 0.2 Hz, the constitutive law is of type E. Parameters of the elastic contribution are kept constant.

Figure 12 :

 12 Figure 12: Comparison of the performance of Ogden, Sussman-Bathe and penalty volumetric potentials for parenchyma samples from the same rat. The volumetric potential is always associated to the same isochoric potential of type B (see Tab. 3).

Fig. 14

 14 Fig. 14 displays an example of hyper-versus viscohyperelastic identification. For force versus displacement curves in the direction of the applied load, we can see that a hyperelastic potential represents the non-linearity of the material but only adding a viscous contribution can lead to the correct hysteresis.

Figure 13 :

 13 Figure 13: Comparison of the performance of different forms of hyperelastic potentials for various frequencies, among all the tested samples (12 at 0.8 Hz, 12 at 0.4 Hz, 14 at 0.2 Hz, all from three different rats). The average ± SD of the error is displayed. * corresponds to p ≤ 0.05.

Figure 14 :

 14 Figure 14: Comparison of the force versus displacement curve along axis x for one sample of rat G and for several nodes: experimental, computed with a hyperelastic constitutive law B (κ = 57.6 kPa, C 1 = 1.36 kPa and C 3 = 71.9 kPa) and with a visco-hyperelastic law E (κ = 62.4 kPa, C 1 = 2.78 kPa, C 3 = 41.8 kPa, τ = 0.293). The hyperelastic curve has been smoothed not to represent some artifacts created by vibrations. The experimental displacement is computed by image registration.

Fig. 15

 15 Fig. 15 displays all the results obtained by identification of the parameters of the hyperelastic potentials and of the corresponding viscoelastic laws. We can see the significant increase in accuracy obtained with the viscoelastic laws, especially when the frequency decreases. Mann-Whitney post-processing of the error values reveals that adding a viscous contribution significantly increases the quality of the identification; the p-value is below 0.0001 for all couples of laws (A, B, C, K, L vs. D, E, F, L, M respectively).

Figure 15 :

 15 Figure 15: Distribution of the error obtained for all the samples for the identified hyperelastic potentials, with and without a viscous contribution. Laws D, E, F, L and M (triangles) are viscoelastic laws that correspond to the hyperelastic potentials A, B, C, K, P (circles) modified with a viscous contribution. X-axis: G, I, J are the three tested rats.

Table 1 :

 1 Samples description.

	Name Number of samples per frequency
		0.2 Hz 0.4 Hz	0.8 Hz
	G	7	5	4
	I	4	4	4
	J	3	3	4

Table 3 :

 3 Identified hyperelastic potentials.

Table 4 :

 4 Identified parameters (average ± SD) for purely elastic potentials.

	Law Parameter	Value	Unit
		κ	7.25e4 ± 2.74e4	Pa
		C 1	2.45e3 ± 8.95e2	Pa
	A	C 2	3.69e3 ± 1.96e3	Pa
		C 3	3.22e4 ± 1.81e4	Pa
	B	C 1	2.87e3 ± 9.05e2	Pa
		C 3	3.83e4 ± 1.63e4	Pa
	K	C 1	3.43e3 ± 1.05e3	Pa
		C 4	1.34e5 ± 8.30e4	Pa
	P	k 1	4.70e3 ± 1.60e3	Pa
		k 2	5.90 ± 1.51	-

Table 5 :

 5 Identified parameters for viscoelastic potentials.

			69e3	Pa
		k 2	4.86 ± 1.42	-
	Frequency Parameter	Value	Unit
	0.8 Hz		5.06e-2 ± 2.53e-2
	0.4 Hz	τ	1.23e-1 ± 3.87e-2	s
	0.2 Hz		2.57e-1 ± 4.22e-2

Acknowledgement 1. Support by the German Science Foundation/Deutsche

Forschungsgemeinschaft DFG (Wa 1521/9 and Uh 88/7) is gratefully acknowledged.

with respect to law M and ≤ 0.0001 with respect to laws E and L) . Tab.

5 displays the identified values for the parameters of the viscoelastic laws.

The parameters related to the elastic potentials (C i , κ, k i ) are averaged over all the samples and all the frequencies as the influence of the frequency is The decoupled form of the potential requires expressing modified strain tensors. We define F as the deformation gradient and J = det F is the volume change. The Cauchy-Green strain tensor is therefore C = F T F. In the case of a decoupled hyperelastic potential, we apply a volumetric/isochoric multiplicative split of the deformation.

where F and C are the modified deformation gradient and Green-Lagrange tensor respectively, representing the isochoric part of the deformation. The strain invariants are defined as:

The modified invariants associated with the isochoric contribution are defined as :

. Hyperelastic potentials

Decoupling an isotropic hyperelastic potential leads to the following form:

Adding a viscous contribution to both the isochoric and the volumetric part of the hyperelastic potential corresponds to:

where Γ vol and Γ are kinematic internal variables that characterize the non-equilibrium state of the deformation process. It can be further split the following way:

For the next section, we define the following quantities:

. Stress tensors

According to the work of [START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF], the elastic second Piola-Kirchhoff stress tensors are defined by: .13) where P = I-1 3 C -1 ⊗C is a fourth-order projector. As mentioned before, the fictitious stress tensors Q vol and Q iso associated with these elastic stress tensors follow an evolution equation given by Eq. ( 7).

Appendix A.4. Elasticity tensors

According to the work of [START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF], the elastic contribution of the elasticity tensors are defined by: A.15) where:

The viscous contribution can be calculated from the elastic contribution,

. Fictitious stress tensor

To solve the differential equation ( 7), we apply the following method:

Once the elastic stress tensor S is computed at time step n + 1, we use a Newton scheme to compute the fictitious stress tensor at step n + 1, using the history of these tensors:

with

Viscous elasticity tensor

As described in [START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF], the overall elasticity tensor at time step n + 1 is calculated as the sum of the elastic tensor and the viscous tensor, which leads to:

where δ = β exp(-∆t/2τ ) as we consider the viscous constants τ and β to be the same for isochoric and volumetric contributions.