Peter M Pinsky
email: pinsky@stanford.edu

Keywords: anisotropic elasticity, biological materials, computer modelling, diffusion, elasticity, ionic transport, porous media porous media, nonlinear, multi-phasic, ionic transport, corneal swelling

A nonlinear macroscopic multi-phasic model for describing interactions between solid, fluid and ionic species in biological tissue materials

Long-Yuan Li, Peter

Multiphysics problems involving fluid-solid, chemical-fluid-solid, and electrochemical-fluid-solid interactions often arise in porous materials. A porous material is a solid (often called the matrix or skeleton) permeated by an interconnected network of pores (voids) filled with a fluid which may be a pure solvent or an ionic solution. Usually both solid matrix and the pore network (also known as the pore space) are assumed to be continuous, so as to form two interpenetrating continua such as in a sponge. Many natural substances such as rocks, soils, biological tissues, and manmade materials such as foams and ceramics can be considered to be porous materials. Porous materials whose solid matrix is elastic with a viscous fluid occupying the pores are referred to as poroelastic materials. A poroelastic material is characterized by its porosity and permeability, as well as the properties of its constituents (solid matrix and fluid).

Poroelasticity theory describes the interaction between the flow and deformation in porous materials. The concept of a porous material originally emerged in soil mechanics, and in particular in the works of Terzaghi (1923), who introduced the concept of effective stress to describe the equilibrium problem of consolidation in soil mechanics. However a more general concept of a poroelastic material, independent of its nature or application, is attributed to [START_REF] Biot | General theory of three-dimensional consolidation[END_REF][START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF][START_REF] Biot | Nonlinear and semi-linear rheology of porous solids[END_REF]. The Biot theory of 2 poroelasticity provides a complete and general description of the mechanical behaviour of poroelastic materials, which has been widely used in civil, mining and petroleum engineering, acoustic wave propagation in saturated media and recently in biology and biomedical engineering.

Biot's equations for the linear theory of poroelasticity are derived from equations of linear elasticity for the solid matrix together with Darcy's law for the flow of fluid through the porous matrix. Since the early 1980s, interest has intensified concerning the role of poroelasticity on the mechanical behaviour of biological tissues. Several models have been proposed to investigate the biomechanics of soft tissues based on Biot's poroelastic theory. They include, among others, successful applications in the study of articular cartilage [START_REF] Mow | Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments[END_REF][START_REF] Mow | Fluid transport and mechanical properties of articular cartilage: a review[END_REF], cardiac muscle [START_REF] Yang | The possible role of poroelasticity in the apparent viscoelastic behavior of passive cardiac muscle[END_REF], arterial walls [START_REF] Simon | Finite element models for arterial wall mechanics[END_REF][START_REF] Simon | Porohyperelastic-transportswelling theory, material properties and finite element models for large arteries[END_REF], skin and bones [START_REF] Mak | A biphasic poroelastic analysis of the flow dependent subcutaneous tissue pressure and compaction due to epidermal loadings: issues in pressure sore[END_REF][START_REF] Mak | Numerical simulation of streaming potentials due to deformation-induced hierarchical flows in cortical bone[END_REF]. All these models have been based on the concept of tissues as structurally biphasic, i.e. composed of a solid matrix saturated by an interstitial fluid. By incorporating material nonlinearity of the solid matrix into the biphasic framework, the model has been able to describe the nonlinear response of porous materials [START_REF] Holmes | Finite deformation of soft tissue: analysis of a mixture model in uniaxial compression[END_REF][START_REF] Cohen | A transversely isotropic biphasic model for unconfined compression of growth plate and chondroepiphysis[END_REF][START_REF] Soulhat | A fibril-network-reinforced biphasic model of cartilage in unconfined compression[END_REF][START_REF] Soltz | A conewise linear elasticity mixture model for the analysis of tension-compression nonlinearity in articular cartilage[END_REF].

For most biological tissues, the fluid phase contains ionic species and the solid phase may also carry fixed electric charges. In some cases the spatial concentration gradient of ionic species in the fluid phase can generate significant osmotic pressure, which in magnitude may be comparable to the fluid pressure. In these cases the transport of ionic species in the fluid phase also needs to be considered in the model as the osmotic pressure can influence the fluid flow and thus the deformation of the porous material. [START_REF] Lai | A triphasic theory for the swelling and deformation behaviours of articular cartilage[END_REF] proposed a triphasic model of articular cartilage as an extension of the biphasic model, where negatively charged proteoglycans are modelled to be fixed to the solid matrix, and monovalent ions in the interstitial fluid are modelled as additional fluid phases. This model was later extended to incorporate multiple polyvalent ions by [START_REF] Gu | A triphasic analysis of negative osmotic flows through charged hydrated soft tissues[END_REF][START_REF] Gu | Transport of multi-electrolytes in charged hydrated biological soft tissues[END_REF]. Similar modification of the biphasic model to include thermo-chemo-electro-mechanical interactions in saturated charged porous solids was made by Huyghe andJanssen (1999), van Meerveld et al. (2003) and [START_REF] Zhang | Transportation of neutral solute in deformable, anisotropic, soft tissue[END_REF]. Finite element analyses based on triphasic models were undertaken by [START_REF] Ferguson | Fluid flow and convective transport of solutes within the intervertebral disc[END_REF] and [START_REF] Snijders | Triphasic finite element model for swelling porous media[END_REF] for the prediction of biomechanical behaviour of intervertebral disc tissues and other swelling porous media.

In this paper, a nonlinear, macroscopic multi-phasic model is proposed for describing the interactions between solid, fluid, and ionic species taking place in porous materials. In contrast to existing multi-phasic models, the present model not only considers the geometrical and material nonlinearities of the porous material but also uses a more general flow theory describing the transport of ionic solution and ionic species in the medium with or without membranes. The model can be regarded as a direct extension of Biot's poroelasticity from the linear interaction between solid and fluid phases to the nonlinear interactions between solid, fluid and ionic species. The model can be used to predict the interaction between mechanical and physiochemical properties in biological tissues. 3 Consider a porous material that is made of a microscopically inhomogeneous solid matrix filled with a homogeneous fluid undergoing large displacements, which deforms from an initial undeformed configuration into the current deformed configuration (see Figure 1). Let Ω be a representative elementary volume of the porous material. The representative elementary volume is defined as being fixed to the continuous solid frame; hence the solid mass is conserved within the representative elementary volume when it deforms from Ω o in the undeformed configuration to Ω in a deformed configuration. The fluid part of the representative elementary volume occupies the pore space enclosed within the solid frame, and its mass is not necessarily conserved.

Governing equations in the biphasic model

Let Ω s and Ω f be the volumes occupied by the solid and fluid phases in Ω, such that

f s Ω ∪ Ω = Ω
, and Ω so and Ω fo be the volumes occupied by the solid and fluid phases in Ω ο , such that

fo so o Ω ∪ Ω = Ω
. The external surfaces Γ of the representative elementary volume Ω are the sum of the solid portion s Γ and the fluid portion f Γ , such that

f s Γ ∪ Γ = Γ
. The porosity φ of the representative elementary volume thus is defined as [START_REF] Lopatnikov | Variational formulation of fluid infiltrated porous material in thermal and mechanical equilibrium[END_REF])

Γ Γ = Ω Ω = f f φ (1)
Let X denote the position vector of a material point of the solid matrix in the porous material in the undeformed configuration and x=ϕ ϕ ϕ ϕ(X) represents the position vector of the material point X in the deformed configuration. The function ϕ ϕ ϕ ϕ(X) is called the deformation mapping from the initial undeformed configuration to the current deformed configuration. By introducing the deformation gradient tensor F, where

X φ F ∂ ∂ = (2)
the right and left Cauchy-Green tensors, C and B, can be expressed in terms of the deformation gradient tensor as follows

F F C T = , T FF B = (3)
The displacement vector u of the material point X is given by the difference between its current position and its initial position, so

X x X X φ X u - = - = ) ( ) ( (4)
The Green strain tensor E can be expressed in terms of the deformation gradient tensor or the displacement gradients as follows where the operator o ∇ is the gradient with respect to the material coordinates and I is the second order unit tensor. It is clear from equations (2)-( 5) that if the displacement vector at a material point is known then the deformation gradient tensor, the Cauchy-Green tensors and the Green strain tensor at that material point can be determined.

] ) )( ( ) ( ) [( 2 1 ) ( 2 1 ) ( 2 1 T o o o T o T u u u u I C I F F E ∇ ∇ + ∇ + ∇ = - = - = (5)
The equilibrium of the porous material is expressed in terms of the total Cauchy stress tensor σ σ σ σ by,

0 = + ⋅ ∇ b σ in Ω (6)
where b is the body force field per unit volume. Based on the assumption of Biot's poroelasticity, the constitutive equation for a fluid-filled porous material can be expressed as follows [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF][START_REF] Detournay | Fundamental of poroelasticity[END_REF])

( ) p J W T o α - - - ∂ ∂ = - F F E S S 1 (7)
where S is the second Piola-Kirchhoff stress tensor, S o is the initial second Piola- Kirchhoff stress tensor, W is the strain energy density function of the porous material, J=det(F) is the Jacobian of the deformation gradient tensor F, α is a material constant, and p is the fluid pressure. Equation ( 7) is a generalization of the constitutive equation for linear and isotropic porous materials [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF][START_REF] Detournay | Fundamental of poroelasticity[END_REF] suitable for nonlinear and anisotropic porous materials. The first and second terms in the right-hand-side of equation ( 7) represent the stress contributions from the solid matrix and fluid pressure, respectively. The factor in parenthesis is the pull-back of the fluid pressure through the motion to the reference configuration. For anisotropic porous materials, examples of the strain energy density function W can be found in [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF] and [START_REF] Pinsky | Computational modeling of mechanical anisotropy in the cornea and sclera[END_REF]. The relationship between the total Cauchy stress tensor in equation ( 6) and the second Piola-Kirchhoff stress tensor in equation ( 7) is given by

T J FSF σ 1 - = (8)
Equations ( 2)-( 8) are similar to those describing the nonlinear elastic deformation of hyperelastic materials under externally applied loads. The only difference here is the constitutive equation ( 7) involves the fluid pressure. It is for this reason that additional equations are required in order to determine the stresses, strains and displacements. It should be emphasized that the displacements, strains and stresses used here are based on a macroscopic description and referenced to the porous material. Thus the deformations or strains used in the present formulation may also include the contribution from the change of porosity of the porous material, which is not properly addressed in Biot's theory.

Fluid pressure in a porous material depends not only on the microstructure of the porous material but also on the flow of the fluid through the connective pores. Fluid flow in porous materials has been studied for more than a century. The interaction between solid and liquid phases in porous materials was first quantified by Darcy in 1856. In a study of water percolating through sand, Darcy discovered that the flow rate was proportional to the pressure gradient. This empirical relationship was theoretically derived later by other investigators on the basis of mechanical analysis of fluid flow in porous materials [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF]de Boer, 2000), which demonstrated that Darcy's law is valid for most Newtonian liquids at low flow velocity. Fluids in many biological materials and in most rock and soil materials belong to this category and therefore Darcy's law can be applied in these fields.

Let q be the specific discharge vector, which describes the motion of the fluid relative to the solid in a porous material, which is formally defined as the rate of fluid volume crossing unit area of the porous material, and ζ be the change of fluid content of the porous material with respect to the initial state, which is defined as the variation of fluid volume per unit volume of porous material due to permeable fluid mass transport. A positive ζ corresponds to a "gain" of fluid by the porous solid. The relationship between q and ζ can be obtained by considering the volumetric conservation of the fluid in the representative elementary volume as follows [START_REF] Biot | The elastic coefficients of the theory of consolidation[END_REF][START_REF] Detournay | Fundamental of poroelasticity[END_REF])

q ⋅ -∇ = ∂ ∂ o t ζ (9)
where t is time. The flow of fluid within the porous material can be described by Darcy's law. By ignoring the fluid density variation effect, Darcy's law can be simply expressed as follows

p ∇ - = K q ( 10
)
where K is the permeability tensor which may reduce to diagonal form if the porous material is orthotropic or isotropic. Substituting equation ( 10) into (9) yields, ( )

p t o ∇ ⋅ ∇ = ∂ ∂ K ζ (11)
In addition to equation ( 11), an equation of state that links the fluid pressure to the volume deformations of the fluid phase and the porous material is needed in order that both the variation of fluid content and fluid pressure can be determined. The equation of state for a general porous material can be expressed as follows [START_REF] Detournay | Fundamental of poroelasticity[END_REF],

[ ] ) 1 ( - - = J M p α ζ (12)
where M and α are the material constants. For linear problems of small deformation J = 1+Tr(Ε Ε Ε Ε) where Tr(Ε Ε Ε Ε) is the volume dilatation of the porous material. Equations ( 11) and ( 12) provide two additional equations which together with equations ( 5)-( 7) can be used to determine the stresses, strains, displacements, variation of fluid content and fluid pressure in the porous material.

For most biological materials, both the solid and fluid phases of the porous material can be regarded as incompressible. In this case the change of volume of the porous material is caused only by the "gain" or "loss" of the fluid, that is,

ζ = -1 J
. In this

p t J o ∇ ⋅ ∇ = ∂ ∂ K (13)
Equations ( 5), ( 6), ( 7) and ( 13) are the governing equations for determining stresses, strains, displacements and fluid pressure in fluid-filled porous materials whose solid and fluid phases are incompressible. The porous material can be nonlinear and anisotropic elastic material. The boundary conditions for stresses and displacements are the same as those normally used in solid mechanics. The boundary conditions for the fluid pressure require that the pressure by prescribed if the boundary surface is permeable, otherwise a zero normal derivative of fluid pressure can be assumed. The governing equations involve the equilibrium, geometrical and constitutive equations of solid materials and the moment balance equation (Darcy's law) of fluids, and so the model is usually referred to as the biphasic model. The linear version of the present biphasic model has been proposed by [START_REF] Naili | A method for the determination of mechanical parameters in a porous elastically deformable medium: applications to biological soft tissues[END_REF] and [START_REF] Swartz | Mechanics of interstitial-lymphatic fluid transport: theoretical foundation and experimental validation[END_REF].

However, the present model is more general as not only is it nonlinear but also anisotropic.

Governing equations in the triphasic model

For most biological materials, the fluid is an ionic solution and thus the transport in the porous material involves not only the flow of ionic solution but also the diffusion, convection and migration of ionic species. Dependent on the properties of the porous material, the permeability of a porous material may be different for the solvent and solutes. If so, osmotic pressure may be generated when the ionic species have concentration gradients. This osmotic pressure can affect the flow of solution and thus the overall deformation of the porous material.

The governing equations describing the transport of ionic solution and ionic species in porous materials can be derived in terms of the mass conservation of individual components in the representative elementary volume as follows

) ( ) / ( q f o o f f t ρ ρ ⋅ -∇ = ∂ Ω Ω ∂ (14) k o o f k t C N ⋅ -∇ = ∂ Ω Ω ∂ ) / ( ( 15
)
where ρ f is the density of the ionic solution, C k is the concentration of species k, and N k is the molar flux of species k flowing in unit time through unit area of porous material. For a dilute solution ρ f can be taken as the density of the solvent and thus is a constant. Note that, 16) that equation ( 14) is equivalent to equation ( 9) when ρ f is constant. However, for the porous material of incompressible solid and fluid phases, ζ can be replaced with J-1. Also, for a macroscopic model it is more convenient to use ζ than Ω f . Note that

o o fo o fo f o f φ ζ + = Ω Ω + Ω Ω - Ω = Ω Ω (16)
so s Ω ≡ Ω
for a porous material with an incompressible solid phase. In this case the change of porosity can be calculated directly from the variation of fluid content ζ as follows

φ φ φ φ φ ζ - - = - Ω Ω = Ω Ω - Ω = 1 o o o o fo f (17)
Note that equation ( 17) is slightly different from that proposed by [START_REF] Gu | Transport of fluid and ions through a porouspermeable charged-hydrated tissue, and streaming potential data on normal bovine articular cartilage[END_REF][START_REF] Gu | Transport of multi-electrolytes in charged hydrated biological soft tissues[END_REF]. In [START_REF] Gu | Transport of fluid and ions through a porouspermeable charged-hydrated tissue, and streaming potential data on normal bovine articular cartilage[END_REF][START_REF] Gu | Transport of multi-electrolytes in charged hydrated biological soft tissues[END_REF], the volume dilatation is taken as the solid phase volume dilatation, whereas in the present formulation, it is taken as the porous matrix volume dilatation. It is believed that the present volume dilatation definition is more appropriate because the strains are calculated from the displacements which are the macroscopic average displacements of the porous material rather than the displacements of the solid phase of the porous material in the representative elementary volume.

The fluxes of the solution and solutes which appear in equations ( 14) and ( 15) can be expressed in two different ways. One is to use the common approach in which the flux of the solution is defined by Darcy's law, that is, equation ( 10), and that of the solutes is defined by the Nernst-Planck equation (Newman, 1991), that is,

q D D N k k k k k k k C C RT F z C + Φ ∇       - ∇ - = (18)
where D k and z k are the diffusion tensor and charge number of species k, respectively, F is the Faraday constant, R is the universal gas constant, T is the absolute temperature, and Φ is the electrostatic potential. For an isotropic porous material D k reduces to diagonal form. The fluxes defined in equations ( 10) and ( 18) do not include those generated by the osmotic pressure and thus they can be applied to porous materials that do not have membrane-like feature, which has no osmotic pressure. The other is based on the theory of irreversible thermodynamics [START_REF] Kedem | Thermodynamic analysis of the permeability of biological membranes to non-electrolytes[END_REF][START_REF] Li | Transport of multicomponent ionic solutions in membrane systems[END_REF] in which the fluxes of the solution and solutes across a permeable membrane are expressed as follows C is the average concentration of species k (see Figure 2), and M is the total number of species involved in the solution. For a semipermeable membrane of thickness h m , the gradient of a variable across the membrane can be approximated as

            ∆Φ + ∆ - ∆ - = ∑ = M k k k k k p C RT F z C RT p L 1 σ q (19)       ∆Φ + ∆ - - = k k k k k k k C RT F z C RT C ω σ q N ) 1 ( (20
m h / (.) (.) ∆ = ∇
. Thus equations ( 19) and ( 20) can be re-written as

            Φ ∇ + ∇ - ∇ - = ∑ = M k k k k k p C RT F z C RT p 1 σ L q (21)       Φ ∇ + ∇ - - = k k k k k k k C RT F z C RT C ω q N ) 1 ( σ (22) where n L m p p h L =
is the equivalent hydraulic conductivity coefficient matrix,

n ω m k k h ω =
is the equivalent permeability coefficient matrix of species k, and n is the normal of the membrane. Equations ( 21) and ( 22) can be compared with Darcy's law (equation ( 10)) and the Nernst-Planck equation (equation ( 18)), where approaches is that in the Nernst-Planck case σ k = 0, while in the current approach it is not necessarily zero. Thus, the former can be regarded as a special case of the latter approach. For a porous material which has equal permeability to the solvent and solutes (i.e., there is no membrane-like feature) one can simply take σ k = 0. Equations ( 21) and ( 22) imply that a porous material with membrane-like feature can be regarded as an assembly of a large number of membranes.

Equations ( 21) and ( 22) show that the fluxes of solution and solutes are coupled. The first two terms in the right-hand-side of equation ( 21) represent the fluxes resulting from the gradients of fluid pressure and osmotic pressure, respectively. It is seen from equation ( 22) that the convective velocities of solutes are not necessarily equal to the solution velocity. This is because the permeability of the solvent and solutes in the porous material may be different. The solutes are normally more hindered by the pore surfaces and thus 0 ≤ σ k ≤ 1.

Note that for charged ions the molar fluxes of ionic species should satisfy the conservation equation of current, which is expressed by

0 ) 1 ( 1 1 =             Φ ∇ + ∇ - - = = ∑ ∑ = = M k k k k k k k k M k k k C RT F z C RT C z F z F ω q N i σ ( 23
)
where i is the current density. Equation ( 23) can be used to determine the gradient of the electrostatic potential. Substituting equations ( 21) and ( 22) into ( 14) and ( 15), yields,

            Φ ∇ + ∇ - ∇ ⋅ ∇ = ∂ ∂ ∑ = M k k k k k o C RT F z C RT p t J 1 σ K K (24) F o r P e e r R e v i e w O n l y [ ]                   Φ ∇ + ∇ - ∇ - ⋅ ∇ +             Φ ∇ + ∇ ⋅ ∇ = + - ∂ ∂ ∑ = M k k k k k k k o k k k k o k o C RT F z C RT p C C RT F z C C J t 1 ) 1 ( ) 1 ( σ σ φ K D (25)
where

RT k k ω D = , p L K =
, and 1 -= J ζ have been applied. Equations ( 5), ( 6), ( 7), ( 24) and ( 25) are the governing equations for determining stresses, strains, displacements, fluid pressure, and ionic concentrations in ionic solution-filled porous materials. Unlike the biphasic model where only the interaction between the solid and fluid is considered, the above model involves the interactions between the solid, ionic solution and ionic species, and so is referred to as the triphasic model. If finite element methods are used to solve the coupled system, the primary variables can be taken to be the displacements of the matrix u, the variation of fluid content ζ, and the concentrations of ionic species C k ; all other quantities can be expressed in terms of these primary variables.

Application of the triphasic model in prediction of corneal swelling

As an application of the proposed triphasic model, we consider a case where a human cornea is subject to an osmotic perturbation of NaCl solution applied to its endothelial surface. Assume that the tissue is initially equilibrated with the external bathing solution. Thus an increase in concentration in the bathing solution will cause an outward flow from the tissue and thus a reduction of the corneal thickness. For simplicity, we assume that the deformation is small and the strains parallel to the corneal surface are negligible and thus the problem can be simplified into a onedimensional problem in which all unknowns are only a function of the coordinate along the thickness direction. In this case, the strain-displacement relationship (5), equilibrium equation ( 6) and constitutive equation ( 7) become

dx du xx = ε (26) = xx σ constant (27) p xx o xx - + = - ε µ λ σ σ ) 2 ( ) ( ( 28
)
where u is the displacement, ε xx is the strain, σ xx is the total stress, σ o is the initial stress which is generated by the swelling pressure p sw , λ and µ are the Lame constants.

The swelling pressure can be expressed in terms of the stromal hydration H w which is defined as the mass of water in unit mass of dry tissue material. Hence, the initial stress can be expressed as ) exp( 10 41 . 2 (Li and Tighe, 2006). The normal thickness of the cornea is about 0.5 mm, the maximum increase or decrease of the corneal thickness in a perfusion experiment is only about 20 µm (Li and Tighe,200), which gives a maximum strain of 4%. This justifies the above assumptions of small deformation and linear version of constitutive equations when considering that the matrix is a soft material which has very low elastic modulus. The inclusion of swelling pressure in equation ( 28) is based on the fact that, in the absence of fluid pressure, an external compression stress must be applied in order to 27) can be identified, which in value is equal to the intraocular pressure (IOP = 2.67x10 4 dyn-cm -2 ) but in an opposite direction (Li and Tighe, 2006). Thus, from equations ( 27) and ( 28) the fluid pressure can be expressed as follows

6 w sw o H p - × - = - = σ in dyn-cm -2
sw xx p IOP p - + + = ε µ λ ) 2 ( ( 29
)
Equation ( 29) implies that the strain-induced stress plus IOP will balance the swelling pressure and the fluid pressure. It is obvious that if the elastic stiffness of the tissue is very low, then the first term in the right-hand-side of equation ( 29) may be negligible.

In this case the fluid pressure is only dependent on the IOP and swelling pressure, and therefore the tissue deformation can be determined purely from the transport model of ions and ionic solution which has been studied by several researchers [START_REF] Klyce | Numerical solution of coupled transport equations applied to corneal hydration dynamics[END_REF][START_REF] Li | Mathematical modelling of corneal swelling[END_REF]Li and Tighe, 2006). Assume that the reflection coefficients for chloride and sodium ions are identical, i.e., σ Cl = σ Na = σ kk . By using electro-neutrality condition (i.e., C Cl = C Na = C), equation ( 23) can be simplified as

( ) ( ) x C C D D D D C D z C D z x C D z x C D z x RT F Na Cl Na Cl Na Na Na Cl Cl Cl Na Na Na Cl Cl Cl ∂ ∂ + - = + ∂ ∂ + ∂ ∂ - = ∂ Φ ∂ 1 2 2 (30)
where z Na =1 and z Cl =-1are the charge numbers of Na and Cl ions, D Na and D Cl are the diffusion coefficients of Na and Cl ions. Using equation ( 30), equations ( 24) and ( 25) can be simplified as

        ∂ ∂ - ∂ ∂ = ∂ ∂ 2 2 2 2 2 x C RT x p k t kk xx xx σ ε (31) [ ]         ∂ ∂ - ∂ ∂ - + ∂ ∂ + = + ∂ ∂ 2 2 2 2 2 2 2 ) 1 ( 2 ) ( x C RT x p C k x C D D D D C t kk xx kk Cl Na Cl Na o xx σ σ φ ε (32)
where k xx is the permeability coefficient in the thickness direction. Using equation (31) to simplify (32), it yields, 29), ( 31) and ( 33) are the three equations used to solve for the fluid pressure, strain, and ionic concentration. Note that the variation of the stromal hydration is related to the strain. The relationship between the two is as follows,

        ∂ ∂ - ∂ ∂ - ∂ ∂ + = ∂ ∂ + 2 2 2 2 2 2 2 2 ) ( x C RT x p C k x C D D D D t C kk xx kk Cl Na Cl Na o xx σ σ φ ε (33) Equations (
d w wo wo w d w wo w wo w xx H H H H H H ρ ρ ρ ρ ρ ε / / 1 / / ) ( + - = + - = (34)
where ρ w = 1 g-cm -3 is the density of water, ρ d = 1.49 g-cm -3 is the density of the dry tissue material, and H wo =3.4 is the initial stromal hydration. For given initial condition (stromal hydration H wo and ionic concentration C o ) and boundary conditions (fluid pressure p or dp/dx and ionic concentration C or dC/dx), one can solve the partial differential equations ( 31) and ( 33) where the fluid pressure and stromal hydration can be expressed in terms of the strain by using equations ( 29) and ( 34). It should be pointed out that although the governing equations ( 29) and ( 31) for the solid and fluid phases are linear, the transport equation ( 33) for the ionic species is nonlinear due to the involvement of the change of porosity. Therefore, the coupled system of equations is not only time-dependent but also nonlinear.

Figure 3 shows a typical example of the corneal response to the shock of NaCl hypertonic solution applied (after the steady state has been reached) to the endothelial surface for a hour duration, while the epithelial surface of the cornea is left intact. The transport parameters employed here are given in Table 1 which are the same as those used in [START_REF] Li | Transport of multicomponent ionic solutions in membrane systems[END_REF]Tighe (2004, 2006a). The only difference is the equation ( 29) for fluid pressure; here it is also dependent on the volume dilatation of the tissue, while in [START_REF] Li | Transport of multicomponent ionic solutions in membrane systems[END_REF]Tighe (2004, 2006a) 0 2 ≡ + µ λ was assumed. The Lame constants for a highly hydrated proteoglycan-containing matrix are estimated as 15 2 ≈ + µ λ kPa [START_REF] Pinsky | A microstructurally based finite element model of the incised human cornea[END_REF][START_REF] Hjortdal | Regional elastic performance of the human cornea[END_REF]Li and Tighe, 2006b). In practice, however, the value of λ+2µ in the thickness direction would be greater than 15 kPa because of the lamellae interlacing. Therefore, it is thought to be more reasonable to use λ+2µ=20 kPa in the example. In order to show the influence of the elastic properties of the tissue material on its swelling in the collagen fibril direction, a large value of λ+2µ=100 kPa is also used. The plotted results show that the thickness changes to the shock are smaller when the strain-induced stress is included in the model. This is to be expected. As the tissue shrinks (or swells) due to the outward (or inward) volumetric flow, the solid matrix of the tissue will provide a "spring-like" mechanical resistance. This mechanical resistance will produce an opposite direction strain, which leads the overall strain be reduced. The results shown in figure 3 indicate that, the larger the value of λ+2µ, the smaller the changes in thickness. It is expected that, when the value of λ+2µ becomes large enough so that the mechanical stress and the swelling pressure are in balance, there will be no change in the thickness.

It should be mentioned that the example shown here only demonstrates the influence of the elastic modulus of the matrix material on the response of the corneal thickness to the perfusion of external solutions. The interaction feature involving the nonlinear material characteristics of the cornea has to be investigated through a full threedimensional model which is currently being developed. Also, it should be pointed out that, because of the lack of experimental data all material constants used here in the example are the same as those used in the previous transport models (see, for example, [START_REF] Klyce | Numerical solution of coupled transport equations applied to corneal hydration dynamics[END_REF]Li and Tighe, 2006a). Since most of these material constants were determined by using a best fit between the transport model, which did not consider the influence of solid phase, and the perfusion experiments, this makes it difficult to validate the present model. Nevertheless, the present results have clearly shown that if the solid phase influence is considered the corneal thickness response will be reduced remarkably. This suggests that the previously determined material constants need be re-assessed. It is noted than even for the one-dimensional model considered here, there are about 20 material constants that are required to determine. A nonlinear, macroscopic multi-phasic model has been proposed for describing the interactions between solid, fluid, and ionic species taking place in porous materials. Governing equations have been derived based on nonlinear theories of solid mechanics, linear flow theory of Newtonian fluids, and the theory of irreversible thermodynamics for the transport of ions and ionic solution. The model is built up on Biot's poroelasticity but considers the nonlinearity of the porous material and the influence of transport of ionic species on the flow of ionic solution and the deformation of the porous material.

Conclusions

The present model has shown that the coupling between solid, fluid and ionic species exists when the porous material involves membranes or has membrane-like features, which could be inside or just on the boundaries of the material. Otherwise, the coupling occurs only between the solid and fluid phases in which case the problem can be solved using the biphasic model.

As the application of the present multi-phasic model, a numerical example of the human cornea under the shock of NaCl hypertonic solution applied to its endothelial surface has been investigated. The results have demonstrated that the mechanical properties of the tissue have important influence on the swelling of the cornea. Without taking into account this influence, the predicted swelling could be exaggerated. Thank you very much for sending me the valuable comments and suggestions made by the three reviewers. We have revised our paper in the light of these comments and suggestions and the following are the responses (text with blue colour) to the reviewers' comments and the new texts added in the revised paper are highlighted with yellow colour.

Referee: 1

"The some derived equations in Section 2 and Section 3 have nonlinear expressions, but the discussion in Section 4 basically made them to be linear expressions. Hence Figure 3 is similar to that of [START_REF] Li | Transport of multicomponent ionic solutions in membrane systems[END_REF]Tighe (2004, 2006a), which the nonlinear characteristics focused on was not displayed completely."

There are three sets of governing equations (each phase has one set). Because of the coupling, governing equations in liquid and ionic species phases are also nonlinear.

Only simplification used in the example section is the governing equations in solid phase which are linearised, while in sections 2 and 3 they are nonlinear. This simplification is purely because the use of one-dimensional model (normal to the thickness), in which case the elastic modulus is only due to the matrix material (which is almost linear) and the collagen fibrils have no contribution to the modulus of elasticity. The main purpose of Fig. 3 is not to show the nonlinear characteristics but to demonstrate how the elastic modulus can affect the thickness change during a numerical perfusion experiment. These two points have been re-addressed in the text of the example section (see added text in pages 9 and 11). It should be pointed out that [START_REF] Li | Transport of multicomponent ionic solutions in membrane systems[END_REF]Tighe (2004, 2006a) model did not have the solid phase while the present model including the example has although in the example only linear equations of the solid phase are considered (see also our response 3 to Referee 2).

"Journal readers may not have more profound physical or mechanical basis, and are not good at understanding complex formula. If it is allowed, the manuscript may give some clearer example of verification, then the readers can easily understand."

See our response to Referee 3.

Referee: 2

1) Some of the equations in section 2 are textbook material and can be taken out by referring to a textbook.

It is true that Eqs.( 2)-( 6) can be found in some high level textbooks. However, they are used a few times during the derivation of other equations in sections 3 and 4 and therefore it is better to keep them. If they were removed it would be difficult for readers to follow.

2) It would be useful to point out what are the basic field variables (such as displacement in an elasticity problem) that need to be solved for after Eqn. (25).

These have been added in text (see page 9). The normal thickness of the cornea is about 0.5mm. The increase or decrease of the thickness during a perfusion experiment is about 20 micro-meters (see Fig. 3), which gives the maximum strain of about 4%. Thus small deformation assumption is justified. This is particularly true for the matrix material which has small elastic modulus. Some text has been added in section 4 (page 9) to justify this.

4) References should be provided for the source of data in Table 1 This has been added in text (see page 16). This has been added in text (see page 14).

Referee: 3 . "The governing equations of the model seem correct to me, but it would be helpful if the authors included a plot, table, or at least a sentence about validation. The previous model (Li and Tighe, 2006a) was able to reproduce experimental data (Fischbarg, J. 1973, Exp. Eye Res. 15, 615-638.) very well and there did not seem to be much room for further improvement (in terms of agreement with experiment). It would be interesting to see how the present model compares with the same or other experiment. This might help readers decide whether to adopt the presented model.

The main purpose of this paper is to demonstrate how important the interaction between solid, liquid and ionic species phases is. As it can be seen from the present tri-phase model, there are two sets of material constants, one for the transport of solution and solutes and the other for the mechanical properties. Currently, the constants for the transport were determined based on the best fit between the transport model which did not consider the solid phase influence and the perfusion experiments (see, for example, [START_REF] Klyce | Numerical solution of coupled transport equations applied to corneal hydration dynamics[END_REF]. This is why "The previous model (Li and Tighe, 2006a) was able to reproduce experimental data (Fischbarg, J. 1973, Exp. Eye Res. 15, 615-638.) very well". The present results show that if the solid phase influence is considered the corneal thickness response will be reduced remarkably. This suggests that the previously determined constants need be re-assessed. Since even for the one-dimensional model there are still about 20 material constants required to determine, this requires more experimental data and is obviously out of the scope of this paper. More explanation has been added in text (see pages 11/12)

The manuscript is organized in a logical way and is written in good English. I have noticed some minor typographical errors: 1)

In abstract: "typical example of …" instead of "typical example about …" 2) Section 4, last paragraph: "Lame constants" instead of "lame constants". Better yet "Lamé constants" These errors have been corrected.

  of the material constants M and α can be taken as ∞ Cheng, 1993). Eliminating ζ from equations (12), equation (11) yields ( )

  porosity of the porous material. It is clear from the relationship between ζ and Ω f given by equation (

  ) where L p is the hydraulic conductivity coefficient, σ k and ω k are the reflection coefficient and permeability coefficient of species k, respectively, , p ∆ , k C ∆ and ∆Φ are the difference of pressures, concentrations and electrostatic potentials in the side of the membrane wall, k

  to D k . The only difference between the two

  swelling. Based on the stress boundary condition at the endothelial surface the stress constant in equation (

Figure 2 .

 2 Figure 2. Sign definition of fluxes in membrane wall (arrow direction represents q and N being positive).

Figure 3 .

 3 Figure 3. Stromal thickness response to the shock of 15 mOsm NaCl hypertonic solution applied (after a steady state) to the endothelial surface for an hour duration while the epithelial surface is left intact.

  section 4, a justification for the small deformation assumption would be helpful, what is the actual level of deformation during the corneal swelling?

5)

  Fig 2, sign definition -sign convention.

Table 1 .

 1 Parametric values employed in the numerical studies for figure 3[START_REF] Li | Transport of multicomponent ionic solutions in membrane systems[END_REF] Tighe, 2004, 2006a).

	Parameter		Na +			Cl -		P
	Charge number			1				-1		-
	Molar mass		11 g/mole		17 g/mole	-
	Initial									2.67x10 4
	conditions	145x10 -6 mole/cm 3	145x10 -6 mole/cm 3	dyne/cm 2
	BCs at epithelium		dC/dx=0		dC/dx=0		dp/dx=0
	BCs at	(145+15)x10 -6	(145+15)x10 -6	2.67x10 4
	endothelium		mole/cm 3		mole/cm 3		dyne/cm 2
	Permeability coefficient in epithelium						2.90x10 -15 cm 4 /(s-dyne)
	Permeability coefficient in stroma						4 8.63x10 -15 H w cm 4 /(s-dyne)
	Permeability coefficient in endothelium						2.65x10 -15 cm 4 /(s-dyne)
	Reflection
	coefficient in		0.800			0.580		-
	epithelium,σ kk
	Reflection
	coefficient in			0				0		-
	stroma,σ kk
	Reflection
	coefficient in		0.685			0.450		-
	endothelium,σ kk
	Diffusion
	coefficient in		4.95x10 -10		7.5x10 -10		-
	epithelium, D		cm 2 /s			cm 2 /s
	Diffusion coefficient in stroma, D	w 91 H . 5	w ρ / 10 6 × -+	H ρ	d w	cm 2 /s	w 0 . H 9	d w ρ H ρ / w 10 6 + × -	cm 2 /s	-
	Diffusion coefficient in endothelium, D		2.63x10 -8 cm 2 /s		4.00x10 -8 cm 2 /s		-
	Active pump rate at endothelium, J			0				-4.0x10 -10 mole/(s-cm 2 )	-

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgement -This work was supported by the UK EPSRC research grant (EP/G056501/1).


Pawel Kuczera

Janusz Wolny

Frank Fleischer

Walter Steurer

Structure refinement of decagonal Al-Ni-Co, superstructure type I

Keywords: quasicrystals, structure analysis user supplied): Edagawa phase, Superstructure type I, decagonal quasicrystals Edagawaphase, superstructure type I, decagonal quasicrystal

come

Introduction

During the last two decades the phase diagram of Al-Ni-Co has been thoroughly investigated. There are eight different modifications of the decagonal phase reported in this system [1 and references therein]. One of them is the so called superstructure type I [2]). The diffraction pattern of the superstructure type I phase consists of strong main reflections and rather weak satellite reflections. The satellite reflections cannot be indexed in the standard decagonal reciprocal lattice basis, i.e. the one used for indexing the so called basic Ni-rich phase. Instead, one has to use a rotoscaled basis [3] (rotated by an angle of π/10 and rescaled by the factor of .Insuch a basis all reflections can be divided into three classes according to:

. Reflections with are called main reflections, with first order satellites (S1) and with second order satellites (S2).No quantitative description of this phase has been published yet.

We present the first refinement results of the Edagawa phase conducted in the physical space only. The structure factor used for the modelling process was calculated on the basis of a statistical approach described in [4,5], which allows a purely 3-dimensional, real space optimization of a quasicrystalline structure. This approach has already been successful in the structure refinement of the basic Ni-rich phase [6,7]. Single crystal X-ray diffraction experiments were performed at the Swiss-Norwegian beam line (SNBL) at the European Synchrotron Radiation Facility (ESRF) in Grenoble, France, using a mar345 image plate detector system (wavelength 0.69393Å, oscillation angle 0.7 o ). Two data sets were collected, one data set (280 frames) with a shortexposure time to prevent saturation of strong reflections and one data set with a longexposure time (244 frames) in order to collect sufficient data for rather weak superstructure reflections. The strongest reflections in the synchrotron diffraction data were still in saturation. Therefore, an in-house data set (974 frames, oscillation angle 0.5 o ) was collected with a four-circle diffractometer equipped with a charge-couple detector system (Xcalibur Oxford Diffraction, MoK α radiation, graphite monochromator).

Diffraction data set and sample details

All three data sets were processed and integrated with the CrysAlisRED software (Oxford Diffraction) and scaled afterwards. The number of reflections present in all three datasets was sufficiently large for obtaining a reliable scale factor.A total number of over 97,000 reflections was merged in the Laue group 10/m with R int = 0.084. As a result, we obtained a set of 6843unique reflections(1387 main reflections and 5456 satellite reflections). There were 702 reflections with zero intensity, which were also used in the refinement process. The optimization was performed against structure factor amplitudes, the standard deviations for strong reflections were calculated as and for weak ones as [8]. For 2284 reflections F was greater than 3σ F .

The structure model

The initial model was built based on the electron density maps obtained from the charge flipping procedure. These maps correspond well to the model proposed by Yamamoto and Weber [9]. This model consists of decagonal clusters decorating the vertices of a rhombic Penrose tiling. The edge length of rhombuses is around 19 Å andthe proposed space group is centrosymmetric P10 5 /mcm.

Based on two independent structure units, i.e. a thick and a thin rhombus, reasonable results can be obtained only for a refinement against the main reflections. However for main reflections only the systematic absences rule was verified.Refinements based on models with non-centrosymmetric space group and four independent structure units, i.e. two pairs of thick and thin rhombuses, gave much better results. A modification of Penrose tiling, that is based on two different pairs of rhombuses is called a Two Colour Penrose Tiling [10]. Regarding the systematic absences, a possible space group for main reflections could be . Nevertheless, for the S1 satellite reflections, the rule is not verified (it holds for the S2 satellites). Interestingly however S1 reflections fulfil another systematic absences rule, namely .This could beexplained by colour symmetry [11]. In this paper, however, the space group P is chosen (which explains 10/m Laue group choice). To keep the reflections-toparameters ratio reasonable, we decided to deflate the rhombuses (withrespect to those proposed by [9]) according to the inflation rules for Penrose tiling. The edge length of our structure units equals12.1Å ( ), the period c along the quasi-tenfold axis equals 4.085 Å. The initial decoration of both pairs of rhombuses was the same. Figure 1 shows the initial positions of atoms and illustrates the refined TM concentration. Squares and circles distinguish atoms at z = 0.25c and

Structure factor

The structure factor was calculated inphysical space on the basis of our statistical approach. The concept of average unit cell for Penrose tiling and a detailed derivation of the structure factor for any decoration is extensively discussed in [5]. Here we shall only show the resulting formula:

F t α (k)
is the Fourier transform of a triangular probability distribution associated with a given structure unit at a given orientation. There are four structure units -two thick and two thin rhombuses (sum over t) and five possible orientations of each structure unit (sum over α).The position of an atomj in a given structureunit t at orientation α is represented by r j α,t , p j t is the weight of a given atom in a structure unit (calculated as a fraction of the atom, which is inside a given rhombus, e.g for group1 in the Figure 1a this weight is equal to , etc.), n t is the number of atoms decorating a structure unit, k is a scattering vector and f j t is the average atomic form factor. A formula for f j t reads: The space group chosen for the refinement is P . This requires that all rhombuses within a given pair have the same decoration in all 5 possible orientations (i.e. all parameters describing atoms are the same in all 5 possible orientations). Therefore it is enough to describe a decoration of a given rhombus in just one chosen orientation.

p calc is

Optimized function

We used a conjugated gradient minimization algorithm. The optimized function reads: where: F obs and F calc are the amplitudes of observed and calculated structure factor respectively, σ is the standard deviation of a structure factor amplitude, N atoms is the number of atoms decorating our structure units, d r andd t are the calculated and theoretical concentration of TM atoms respectively, finally c r and c t are the calculated and theoretical point density.

Optimization results

Due to length restrictions of this paper we will not list all 691 parameters. Only the most important parameters i.e. the occupation probabilities and TM atoms concentration are listed in Table 1. All other parameters together with more detailed structure analysis will be given in a separate paper. In Table 1 only few standard deviations are marked. If a standard deviation is not given, then this means that either a parameter was fixed or thestandard deviation was lower than 0.01. All parameters which equal exactly 0 or 1 were fixed during the refinement process. This is due to the fact that after first stages of refinement, if a parameter value was greater than 0.95, then it was fixed to 1, and if it was lower than 0.05, then it was fixed to 0.In Table 1 some positions are marked with -or + sign. These are vertex split positions. Rows marked with -give parameters values layer at z = 0.25c, and rows marked with + for layer at z = 0.75c

In Figure 1 and Table 1 it is seen, that there are 73 atoms decorating pair A of our structure units and 71 atoms decorating pair B. Out of 143 atomic positions there are only 48 not fully occupied (and only 16 with occupation probability less than 0.5). What is more, there are only 22 positions, for which a mixed (Al and TM) occupation is possible. For each of 143 atomic positions an anisotropic Debye-Waller factor as well as a shift form ideal position (in the quasiperiodic plane) were also refined.

The Debye-Waller coefficientb xy i.e. the component in the quasiperiodic plane, is lower than 3 Å 2 for all atomic positions (for most of them it is lower than 2 Å 2 ), which is an acceptable value for inorganic materials and can be interpreted as arising solely from the thermal motion of atoms.However for some atomic positions theb z coefficient, i.e. the Debye-Waller factor component associated with z direction, exceeds the value of 3 Å 2 (for some atoms even exceeds 5 Å 2 ). This could mean that the atomic layers of the structure are puckered rather than flat.Such an assumption is reasonable because of presence of diffuse maxima in the half-integer h 5 layers of the diffraction pattern [1]. In [12] even an 8 Å superstructure for the basic Ni-rich is 1.00 0.00 +28 0.2(2) 0.00 54 1.00 0.00 3 1.00 0.00 27 0.81 0.47 51 1.00 0.24 5 0.00 29 1.00 0.28 55 1.00 0.00 4 1.00 0.00 28 0.74 0.20 52 1.00 0.00 6 1.00 0.29 0.00 30 1.00 0.00 -56 0.2(5) 0.00 5 1.00 1.00 29 1.00 0.00 53 1.00 0.00 7 0.46 0.00 31 1.00 0.00 +56 0.2(9) 0.00 6 1.00 0.00 30 0.48 0.00 54 1.00 0.39 8 1.00 0.00 32 1.00 0.69 57 1.00 1.00 7 1.00 0.00 31 0.50 0.00 55 0.52 0.00 9

1.00 1.00 33 0.73 0.00 58 0.52 0.00 8 1.00 0.00 32 1.00 0.00 56 1.00 0.00 10 1.00 0.00 34 1.00 0.00 59 1.00 1.00 9 1.00 0.00 33 1.00 0.00 57 1.00 0.00 11 1.00 0.00 35 1.00 0.00 60 1.00 0.00 10 1.00 0.00 34 1.00 1.00 58 0.76 0.00 12 1.00 1.00 36 0.90 0.00 61 1.00 0.00 11 0.53 0.00 35 0.76 0.00 59 1.00 1.00 13 0.69 0.00 37 1.00 1.00 62 1.00 1.00 12 1.00 1.00 36 1.00 1.00 60 1.00 0.00 14 1.00 0.00 38 1.00 0.00 63 0.50 0.00 13 0.77 0.00 37 1.00 0.00 61 0.63 0.00 15 1.00 0.00 39 0.43 0.00 64 1.00 1.00 14 1.00 1.00 38 1.00 0.00 62 1.00 assumed. This means that the period along z axis is doubled. However, there are no half-integerh 5 reflections in the diffraction pattern, thus in experiment only an average 4 Å structure is observed. In this approach larger b z values result from overlapping layers of the superstructure having opposite z displacements. In our paper only the average 4 Å structure is investigated. A Fourier map of the refined structure is shown in the Figure 2.Layers at z = 0.25c and z = 0.75c are superimposed.Grey lines show rhombuses chosen as structure units in the optimization process. Letters A and B correspond to pair A and pair B of our structure units respectively. Black lines show rhombuses chosen as structure units in [9], which are inflated (i.e. rescaled by τ) with respect to the grey ones. It is clearly seen that the vertices of these rhombuses (black-lined) forming a Penrose tiling are decorated withquasitenfold clusters. TM atoms are likely to arrange around circles of about 12 Å and 19 Å in diameter, centred at the vertices of the Penrose tiling.

The shortest interatomic distance based on Patterson map analysis is about 2.43 Å. In the refined structure there are some distances shorter than this value. However distances shorter than 2.15 Å exist only between pairs of atoms, for which the sum of occupation probabilities is lower than 1. The distances in the range between 2.15 Å and 2.43 Å could be explained by the puckering of the atomic layers, which is not taken into account in the model and thus in the calculation of the interatomic distances. In the presented approach to structure refinement solely parameters related to physical space were refined, i.e. phasonic displacement parameters were not used neither individually for each atomic position nor globally. Extinction was also negligible. In the end we obtained R = 34.8 % and wR = 6.0 % for all 6843 reflections (R = 29.5 %, wR = 5.9 % for 6141 non-zero structure factor amplitudes), R = 26.8 % and wR = 5.4% for all 4132 main+S1 reflections (R = 23.2%, wR = 5.4 % for 3790 non-zero structure factor amplitudes) and R = 13.3%,wR = 3.6% for all 1387 main reflections (R = 12.1 %, wR = 3.6 % for 1317 non-zero structure factor amplitudes). A detailed R and wR statistics is given in the Table 2.The values of R and wR for satellite reflections only are rather high. One should notice however, that F/σ(F) ratio for these reflections is generally low e.g. there are 2711 S2 reflections but only for 386 of them F is larger than3σ(F).

F calc /F obs plot is shown in the [13] and are equal to 3.94 g/cm 3 ,0.0661 Å -3 respectively.

  An alloy with nominal composition Al 72.0 Co 12.4 Ni 15.6 was prepared by arc melting of the pure elements and heat treated afterwards. According to the phase diagram, the sample was quenched at 730 º C. The chemical composition was determined to Al 70.6(3) Co 13.3(3) Ni 16.1(2) using a JEOL JXA-8200 WD/ED combined microprobe. Backscattered electron images identified the alloy as single phase.

  .75c respectively. At some vertex positions both a circle and a square is visible. These are split positions -partial occupation in both layers is possible for these positions. Grey symbols denote TM atoms, light grey Al/TM atoms and white Al atoms. For simplicity atoms, whose occupation probability decreased to zero during the refinement process are not shown in the figure. The numbers indicated in Figure1are used in the section 'Optimization results' to list the occupation probabilities and TM concentration at each position. For the purposes of this paper the pair of rhombuses in Figure1awill be called pair A and in the Figure1bpair B.

Figure

  Figure 1a.-'pair A' of structure units Figure 1b.-'pair B' of structure units.

  the occupation probability of a given position by an atom, p TM is the concentration of TM atoms at a given position, f TM is the atomic form factor for TM Co and Ni), f Al atomic form factor for Al atoms. D l(s) is the Debye-Waller factor: b xy is the average displacement parameter in the quasiperiodic plane, b z is the displacement parameter in the z-direction (periodic direction). A relation between these parameters and he square mean shift of the atom from the equilibrium position reads:

PF 1 ,

 1 PF 2 , PF 3 , PF 4 keep the concentration of TM atoms, Debye-Waller coefficients, chemical composition and point density respectively in the physical boundaries. Numbers C 1 , C 2 , C 3 , C 4 , are the weights of a given component in the optimized function.

Figure 2 .

 2 Figure 2. A Fourier map of the refined structure. Layers at z = 0.25c and z = 0.75c are superimposed. Letters A and B correspond to pair A and pair B of rhombuses from Figure 1.

Figure 4 .

 4 The chemical composition of the refined structure is Al 70.4 TM 29.6 , point density 0.0668 Å -3 and overall density 4.04 g/cm 3 . Chemical composition is almost equal to the measured one. Accurate density and point density values are known for the basic Ni-rich sample of the composition Al 72 Ni 20 Co 8

Figure 4 .

 4 Figure 4. F calc /F obs plot.

Table 1 .

 1 Refined occupation probabilities and TM concentrations.

				Pair A			Pair B
	No	pcalc pTM No	pcalc pTM No	pcalc pTM No pcalc pTM No pcalc pTM No pcalc pTM
	1	0.52 0.00	26	0.87 0.00	51	1.00 1.00 -1 0.05 0.00 24 1.00 0.00	48 1.00 1.00
	2	0.55 0.58	27	1.00 0.00	52	1.00 0.00 +1 0.11 0.00 25 1.00 0.86	49 1.00 0.00
	3	0.87 0.58 -28	0.10 0.00	53	1.00 0.00	2 1.00 1.00 26 1.00 1.00	50 1.00 0.15
	4

Table 2 .

 2 Structure factor statisticsNo of refl.R [%] wR[%]

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgements

PawelKuczera's PhD fellowshipat the Laboratory of Crystallography, ETH Zurich is financed by SciexNMS ch .

This work has been supported by Ministry of Science and Higher Education of Poland.

Structure optimization calculations was performed in Academic Computer Center CYFRONET AGHon SGI Altix 3700 computerand financed under computational grant MNiSW/SGI3700/AGH/019/2010.

We thank the staff of the Swiss-Norwegian beam line (SNBL) at the European Synchrotron Radiation Facility (ESRF) in Grenoble, France, for their support.

Conclusions

In this paper we have described the first results of the structure analysis of decagonal Al-Ni-Co,superstructure type I (Edagawa phase). The refinement was conducted with only real space parameters,based on a decorated rhombic Penrose tiling model. Reasonable R and wR factors were obtained, density and point density are close do experimental values.In the resulting structure there is a certain degree of disorder within the structure units-thin and thick rhombuses (partial and/or mixed occupation). On the other hand, these rhombuses form the Penrose tiling, thus, in our approach, the structure is also perfectly long-range ordered.


P J Nugent
email: p.nugent@liv.ac.uk

J A Smerdon

R Mcgrath

M Shimoda

C Cui

A P Tsai

H R Sharma

Step-Terrace Morphology and Reactivity to C 60 of the Fivefold i -Ag-In-Yb Surface

Keywords: quasicrystals, surfaces (user supplied): Ag-In-Yb, C60, Scanning Tunneling Microscopy Quasicrystal, Ag-In-Yb, C60, Scanning Tunneling Microscopy, Autocorrelation, Fibonacci Sequence, Surface Morphology, Bias Dependence

The surface of the icosahedral i-Ag-In-Yb quasicrystal provides one of the first non Al-based aperiodic surfaces that is suitable for study under ultra high vacuum conditions. We present a scanning tunneling microscopy (STM) study of the fivefold surface of this new quasicrystal demonstrating detailed structure of the terraces and steps. The analysis of the autocorrelation functions of STM images at opposite bias polarities and of the in-plane structure of the bulk model of i-Cd-Yb, which is isostructural to i-Ag-In-Yb, reveals that the surface terminations occur at the centres of the rhombic triacontrahedral (RTH) clusters, that are the basic building blocks of this material. The study further confirms that the occupied electronic states are enhanced at Ag/In sites and unoccupied states are located on Yb sites. Step edges display a Fibonacci sequence of truncated clusters, which can be also explained in terms of the model structure. Occasionally, a single terrace is found to display different structures at negative bias, whereas the same terrace show a uniform structure at positive bias. Depositing C 60 creates a disordered overlayer on the surface with no resulting FFT or LEED patterns.

Introduction

The discovery of the thermally stable binary icosahedral (i ) quasicrystal, Cd-Yb [1], opened a new area of research in the field of aperiodic materials. Unlike previously discovered thermally stable quasicrystals, all of which are ternary, i -Cd-Yb is composed of only two materials, which has allowed for a full structural solution of an icosahedral quasicrystal [2]. The i -Cd-Yb quasicrystal, whose structure can be considered to be an aperiodic lattice of rhombic triacontrahedral (RTH) clusters, is structurally different from previously studied Al-based quasicrystals, which are formed from Mackay clusters [3]. Al-based quasicrystals display many interesting physical and electronic properties such as hardness, low surface friction and high electrical conductivity [START_REF] Dubois | Useful Quasicrystals[END_REF], and the i -Cd-Yb quasicrystal also displays un- usual physical properties such as large electronic specific heat coefficients and high magnetoresistance [START_REF] Muro | [END_REF][6][7].

The i -Cd-Yb quasicrystal is unsuitable for ultra high vacuum (UHV) study due to the high vapour pressure of Cd, which leads to evaporation on heat treatment, and as such it is difficult to prepare surfaces using the usual surface science methods of ion bombardment and annealing. This issue has been resolved by replacing Cd in this system with equal amounts of Ag and In. This new material, i -Ag-In-Yb, maintains the same electron to atom ratio of 2.0 and is isostructural to Cd-Yb, [8,9] and Ag and In are stable in UHV, which allows the study of this quasicrystal in UHV environments.

In previous work, Sharma et al. [10,11] succeeded in characterising the structure of the fivefold surface of i -Ag-In-Yb using scanning tunneling microscopy (STM), low energy electron diffraction (LEED) and reflection high energy electron diffraction (RHEED). In more recent work, the core levels and valence band structure of the same material have been characterised using x-ray photoelectron spectroscopy (XPS) [12] and ultraviolet photoelectron spectroscopy (UPS) [13]. With a combination of STM data of the fivefold i -Ag-In-Yb surface and the model structure of i -Cd-Yb, it has been possible to determine which bulk planes are responsible for surface terminations. Results indicate that surfaces terminate on the bulk planes which intersect the centre of the RTH clusters.

With the clean surface characterised and its structure understood, it is now possible to use this material as a substrate for thin film growth and adsorption of molecules. There is great interest in growing thin films using quasicrystals as templates and thus exploring the properties of two-dimensional single-element or molecular quasicrystalline structures [14,15].

In this paper we present and characterise several structural features of the fivefold i -Ag-In-Yb surface that were not described in the previous publication [11]. In particular, we will focus on the step-terrace morphology studied by using new high resolution STM images. In addition, we will present a study on the adsorption of C 60 on the same surface.

Experimental Details

The Bridgman method was used to produce a single grain sample of composition Ag 42 In 42 Yb 16 [8,9,16], which was cut perpendicular to the fivefold axis determined by Laue back-scattering. It was then mechanically polished using 6, 1 and 0.25 µm diamond paste. Atomically flat terraces were obtained through cycles of 2.5 KeV Ar + sputtering (30 -60 minutes) at room temperature followed by annealing at 715 K (2 -4 hrs) at a pressure of 1.4 × 10 -10 mbar. Temperatures were monitored using an infra-red pyrometer set at an emissivity of 0.35.

STM experiments were performed in Ultra High Vacuum (UHV) using an Omicron variable temperature STM operating at room temperature. The stated bias refers to to the potential applied to the sample relative to the tip. C 60 was deposited using a home built evaporator with a tungsten filament wrapped tightly around a Pyrex crucible containing C 60 powder. The sample remained at room temperature during deposition. The evaporator was operated at temperatures of 485 to 495 K.

Results and Discussion

It was shown previously that the fivefold surface of i -Ag-In-Yb quasicrystal yields atomically flat terraces with fivefold symmetry [11]. The surface possesses step heights of 0.28 ± 0.04, 0.58 ± 0.03 and 0.85 ± 0.05 nm, with respective percentage occurrences of 66%, 12% and 22%. These heights occur in the model structure as separations between atomic planes that intersect the centres of the RTH clusters.

STM images exhibit a strong bias dependence. Fig. 1 (a) shows an STM image taken with a negative bias, which is composed of protrusions of size 1.30 ± 0.04 nm located at the vertices of pentagons of edge length 2.40 ± 0.15 nm. STM images can be compared directly to the in-plane structure intersecting the cluster centres shown in Fig. 1 (c). The size of the protrusions corresponds to the diameter of the Cd ring, 1.27 nm. Fig. 1 (b) shows an image obtained by applying a positive bias, which corresponds to unoccupied states, and has been taken from a different part of the surface. This image no longer displays protrusions and instead display rings of individual atoms. The rings have a diameter of 1.80 ± 0.09 nm, which corresponds to the diameter of 1.94 nm of the Yb rings from the model. Fig. 2 shows in greater detail the changes that occur as the bias is changed from negative to positive. These images feature a very common structural motif that exists on the surface. At negative bias, this feature appears as a pentagonal arrangement of protrusions. As the bias is changed, it is possible to resolve inside the protrusions. For instance, at +0.1 V bias, pentagons can be resolved at the vertices of the pentagonal arrangement of protrusions. This does not occur often, and can be related to the unknown distributions of Ag and In in the ring structure that a protrusion is constructed from. From the model structure of the approximant of i -Ag-In-Yb model, the surface layer is expected to be composed mostly of In [17], but the exact distribution of Ag and In in the ring is not known. This is probably due to the density of states becoming localised to specific atoms, allowing for them to be resolved. As the bias is made more positive, these protrusions fade and are replaced by a pentagonal arrangement of rings (Fig. 2 (c)).

Autocorrelation patterns from STM images indicate the presence of short and long range order by displaying the distances between recurring features that occur in the images. Fig. 3 (a,b) demonstrate autocorrelation pattern taken from negative and positive biased STM images respectively. The autocorrelation pattern from the negative image displays many rings of maxima with radii 2.3 nm, 3.8 nm, 6.1 nm and so on. These are related by 1, τ and τ 2 respectively, where τ is the golden mean, 1.618. . . , and confirm aperiodic order. The first maximum corresponds to the protrusion separations observed by STM.

The autocorrelation of a positive-bias STM image consists of smaller and more frequent rings of maxima. The first maximum occur at a radius of 0.6 nm which correspond to the atomic separation of Yb atoms, 0.6 nm. The radius of the second maxima at 1.1 nm originates from the next nearest neighbour separation of the Yb atoms, and is τ -related to the first maximum. The third maximum of 1.9 nm corresponds to the diameter of the Yb rings.

Autocorrelation patterns from the model structure are then compared to each other and the experimental results. Fig 3 (c,d) are autocorrelation patterns taken from the Cd-Yb model. Fig. 3 (c) has been determined using the positions of the cluster centres which are associated with the positions of the protrusions. These have been modeled as a point source as opposed to larger objects that match the size of these protrusions, so that a clearer, sharper autocorrelation pattern is produced. The pattern agrees with that of the negative biased STM image shown in Fig. 3 (a), apart from one feature. The first maximum on the theoretical autocorrelation pattern occurs at 1.4 nm which is smaller than the first maximum from the STM autocorrelation pattern by a factor of τ . This is because the model features a smaller cluster separation distance which is not resolved by STM. To compare the model structure to the positive biased autocorrelation, the Yb atom placements have been used. This produces a complex pattern similar to that of the experimental autocorrelation pattern, and as with the cluster centre position autocorrelation pattern, matches up well to the experimental data with maxima at similar positions. The positive bias autocorrelation patterns are a factor of τ 2 smaller than the negative bias autocorrelation patterns. This is consistent with the autocorrelation patterns from the model, and also consistent with fast Fourier transform (FFT) patterns obtained from the respective STM images (which are not shown). Now we discuss the morphology of steps and terraces. At the edge of 0.28 nm high steps it is possible to observe protruding clusters from STM images taken at a negative bias. These protruding clusters can yield a Fibonacci sequence of L and S segments, extending for several tens of nanometers with a rare defect of a triple L sequence, as shown in Fig. 4 (a). These L and S segments measure 4.24 ± 0.27 nm and 2.53 ± 0.22 nm respectively. These two lengths are τ -related, so that L = τ S. L segments occur more often than S segments, with the ratio of occurrences also being related by τ . The Fibonacci sequence can be explained if the step-edge is parallel to the edge of the pentagon features observed in the substrate, which corresponds to high symmetry directions as shown in Fig. 4 (b). These sequences are only observed on steps of the smallest height. This may be due to the larger step height being the product of two small step heights, causing a reconfiguration of clusters at the step edge to accommodate the height difference.

On several terraces imaged at negative bias, two different types of structure are observed as shown in Fig. 5. One type displays a clearer configuration of protrusions than the other. The same area imaged with the same tip but at a positive bias shows a uniform structure as seen in Fig. 5 (b). This may suggest that the non-uniformity encountered at a negative bias is related to an electronic effect.

We have investigated the reactivity of the surface by depositing C 60 on the surface. Fig. 6 (a) shows a sub-monolayer coverage with the area covered by roughly 10% of C 60 . The protrusions observed on the image are 0.46 ± 0.04 nm or 0.65 suppresses the LEED spots from the substrate as surface order is destroyed.

Summary and Conclusions

Scanning tunneling microscopy has been successfully employed to characterise the fivefold i -Ag-In-Yb surface by identifying surface features and comparing these to the model structure of i -Cd-Yb. These surfaces intersect the centre of the RTH clusters that create the bulk of this material. The surface morphology imaged is dependent on the sample bias. At a negative bias, protrusions relating to Ag/In sites, which corresponds to the truncation of RTH cluster, while at a positive bias the Yb atoms are observed. This suggests that the occupied electronic states are enhanced at Ag/In sites and unoccupied states are located on Yb sites. Truncated rhombic triacontrahedral clusters are visible at step-edges forming a Fibonacci sequence. This is expected from the model structure if the step is parallel to one of the high symmetry directions. Different structures have been observed on the same terrace when the sample is at a negative bias. This lack of uniformity is not observed at a positive bias, and may be an electronic effect.

Depositing C 60 on the surface creates a disordered array of molecules, which provides no FFT pattern. Continued deposition leads to a disordered 3D growth mode.

Figure 1 .

 1 Figure 1. (a, b) High resolution STM images of the clean fivefold i-Ag-In-Yb surface taken at -0.81 V and 0.95 V respectively (30 nm × 30 nm). Image (b) has been enhanced using a fast Fourier filter (30 nm × 30 nm). (c) Surface termination plane taken from Cd-Yb model, with Cd featuring as black dots, and Yb as light grey dots. The plane intersects the middle of the RTH clusters, displayed as larger dark grey dots (20 nm × 20 nm).

Figure 2 .

 2 Figure 2. Highlighted feature of the clean surface of i-Ag-In-Yb at varying bias. (a, b, c) taken at -1.0 V, +0.1 V and +1.2V respectively.

Figure 3 .

 3 Figure 3. (a, b) are auto correlation patterns taken from negative and positive bias STM images respectively. These are accompanied by autocorrelation patterns taken from (c) the positions of clusters and (d) the positions of Yb atoms from the i-Cd-Yb model (20 nm × 20 nm).

Figure 4 . 6 Figure 5 .

 465 Figure 4. (a) Step terraced image from the clean fivefold i-Ag-In-Yb surface, with a step edge containing a 100 nm long sequence of protrusions (250 nm × 250 nm). The highlighted step edge on (a) is magnified on the right-hand-side. The indicated L section demonstrates a break in the sequence. (b) Section of the Cd-Yb model with only the Cd atoms displayed with several potential edge terminations marked by solid lines (20 nm × 20 nm).

Figure 6 .

 6 Figure 6. (a) Image of a 0.1 ML deposition of C 60 on the fivefold i-Ag-In-Yb surface. Molecules have two different heights indicated by arrows (50 nm × 50 nm). (b) Image of > 1 ML of deposited C 60 (100 nm × 100 nm).