Pierre Lugan

David Clément

Philippe Bouyer

Alain Aspect

Laurent Sanchez-Palencia

Anderson Localization of Bogolyubov Quasiparticles in Interacting Bose-Einstein Condensates

Keywords: numbers: 05.30.Jp, 03.75.Hh, 64.60.Cn, 79.60.Ht

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

An important issue in mesoscopic physics concerns the effects of disorder in systems where both quantum interference and particle-particle interactions play crucial roles. Multiple scattering of non-interacting quantum particles from a random potential leads to strong Anderson localization (AL) [1], characterized by an exponential decay of the quantum states over a typical distance, the localization length. AL occurs for arbitrarily weak disorder in 1D and 2D, and for strong-enough disorder in 3D [2]. The problem is more involved in the presence of interactions. Strong disorder in repulsively interacting Bose gases induces novel insulating quantum states, such as the Bose [3] and Lifshits [4] glasses. For moderate disorder and interactions, the system forms a Bose-Einstein condensate (BEC) [4,5,6], where the disorder induces a depletion of the condensed and superfluid fractions [7] and the shift and damping of sound waves [8].

These studies have direct applications to experiments on liquid 4 He in porous media [9], in particular as regards the understanding of the absence of superfluidity. Moreover, the realization of disordered gaseous BECs [10,11,12,13,14] has renewed the issue due to an unprecedented control of the experimental parameters. Using optical speckle fields [START_REF] Goodman | Statistical Properties of Laser Speckle Patterns in Laser Speckle and Related Phenomena[END_REF], for instance, one can control the amplitude and design the correlation function of the random potential almost at will [14], opening possibilities for experimental studies of AL [START_REF] Kuhn | [END_REF]17]. Earlier studies related to localization in the context of ultracold atoms include dynamical localization in δ-kicked rotors [18] and spatial diffusion of laser-cooled atoms in speckle potentials [19].

Transport processes in repulsively interacting BECs can exhibit AL [17,20]. However, for BECs at equilibrium, interaction-induced delocalizing effects dominate disorderinduced localization, except for very weak interactions [5,6]. The ground state of an interacting BEC at equilibrium is thus extended. Beyond, one may wonder how the many-body (collective) excitations of the BEC behave in weak disorder. In dilute BECs, these excitations correspond to quasiparticles (particle-hole pairs) described by the Bogolyubov theory [21]. In this case, the interplay of interactions and disorder is subtle and strong arguments indicate that the Bogolyubov quasiparticles (BQP) experience a random potential screened by the BEC density [5]. This problem has been addressed in the idealized case of uncorrelated disorder (random potentials with a delta correlation function) in Ref. [22].

In this Letter, we present a general quantitative treatment of the localization of the BQPs in an interacting BEC with healing length ξ in a weak random potential with arbitrary correlation length σ R . For weak disorder, we introduce a transformation that maps rigorously the many-body Bogolyubov equations onto the Schrödinger equation for a non-interacting particle in a screened random potential, which we derive analytically. We calculate the Lyapunov exponent Γ k (inverse localization length) as a function of the BQP wavenumber k for 1D speckle potentials. For a given ratio ξ/σ R , we determine the wavenumber k max for which Γ k is maximum: We find k max ∼ 1/ξ for ξ ≫ σ R , and

k max ∼ 1/σ R for ξ ≪ σ R .
The absolute maximum appears for ξ ∼ σ R , so that the finiterange correlations of the disorder need to be taken into account. Numerical calculations support our analysis. Finally, possibilities to observe the AL of BQPs in BECs placed in speckle potentials are discussed.

We consider a d-dimensional Bose gas in a potential V (r) with weak repulsive short-range atom-atom interactions, characterized by the coupling constant g. Its physics is governed by the many-body Hamiltonian

Ĥ = dr ( 2 /2m)[(∇ θ) 2 n + (∇ √ n) 2 ] +V (r)n + (g/2)n 2 -µn (1)
where m is the atomic mass, µ is the chemical potential, and θ and n are the phase and density operators, which obey the commutation relation [n(r), θ(r ′ )] = iδ(rr ′ ). According to the Bogolyubov-Popov theory [21,[START_REF] Popov | Functional Integrals in Quantum Field Theory and Statistical Physics[END_REF][START_REF] Petrov | [END_REF], for small phase gradients ( 2 |∇θ| 2 /2m ≪ µ) and small density fluctuations (δn ≪ n c , where n c = n and δn = n -n c ), Hamiltonian (1) can be diagonalized up to second order as

Ĥ = E 0 + ν ǫ ν b † ν bν
, where bν is the annihilation operator of the excitation (BQP) of energy ǫ ν . The many-body ground state of the Bose gas is a BEC with a uniform phase and a density governed by the Gross-Pitaevskii equation (GPE):

µ = -2 ∇ 2 ( √ n c )/2m √ n c + V (r) + gn c (r). (2)
Expanding the density and phase in the basis of the excitations,

θ(r) = [-i/2 n c (r)] ν [f + ν (r) bν -H.c.] and δn(r) = n c (r) ν [f - ν (r) bν + H.c.],
the Hamiltonian reduces to the above diagonal form provided that f ± ν obey the Bogolyubov-de Gennes equations (BdGE) [START_REF] De Gennes | Superconductivity of Metals and Alloys[END_REF]:

-( 2 /2m)∇ 2 + V + gn c -µ f + ν = ǫ ν f - ν (3) -( 2 /2m)∇ 2 + V + 3gn c -µ f - ν = ǫ ν f + ν , (4)
with the normalization dr 2)-( 4) form a complete set to calculate the ground state (BEC) and excitations (BQPs) of the Bose gas, from which one can compute all properties of finite temperature or time-dependent BECs.

f + ν f - ν ′ * + f - ν f + ν ′ * = 2δ ν,ν ′ . Equations (
Here, we analyze the properties of the BQPs in the presence of weak disorder. According to Eqs. ( 3) and ( 4), they are determined by the interplay of the disorder V and the BEC density background n c . Let V (r) be a weak random potential ( V ≪ µ, see below) with a vanishing average and a finite-range correlation function,

C(r) = V 2 R c(r/σ R ), where V R =
V 2 is the standard deviation, and σ R the correlation length of V . As shown in Refs. [4,5,6], the BEC density profile is extended for strong-enough repulsive interactions (i.e. for ξ ≪ L, where ξ = / √ 4mµ is the healing length and L the size of the BEC). More precisely, up to first order in V R /µ, the GPE (2) yields

n c (r) = [µ -V (r)]/g (5)
where V (r) = dr ′ G ξ (rr ′ )V (r ′ ) and G ξ , the Green function of the linearized GPE [4,6], reads G ξ (q) = (2π) -d/2 /[1 + (|q|ξ) 2 ] in Fourier space [START_REF] Here | we define G ξ (q) = dr (2π) d/2 G ξ (r) exp(-iq • r)[END_REF]. Then,

V (q) = V (q)/[1 + (|q|ξ) 2 ]. (6)
Thus, ξ is a threshold in the response of the density n c to the potential V , as V (q) ≃ V (q) for |q| ≪ ξ -1 , while V (q) ≪ V (q) for |q| ≫ ξ -1 . The potential V (r) is a smoothed potential [6]. If V is a homogeneous random potential, so is V , and according to Eq. ( 5), the BEC density profile n c is random but extended [4,6].

Solving the BdGEs (3) and ( 4) is difficult in general because they are strongly coupled. Yet, we show that for a weak (possibly random) potential V (r) this hurdle can be overcome by using appropriate linear combinations g ± k of the 2 and k = |k|. For V = 0, the equations for g ± k are uncoupled [see Eqs. ( 7) and ( 8)] and we recover the usual plane-wave solutions of wave vector k and energy [21]. For weak but finite V , inserting Eq. ( 5) into Eqs. ( 3) and (4), we find

f ± k functions, namely g ± k = ±ρ ±1/2 k f + k + ρ ∓1/2 k f - k , with ρ k = 1 + 1/(kξ)
ǫ k = ρ k ( 2 k 2 /2m)
2 k 2 2m g + k = - 2 2m ∇ 2 g + k - 2ρ k V 1 + ρ 2 k g - k + V - 3 + ρ 2 k 1 + ρ 2 k V g + k (7) -ρ 2 k 2 k 2 2m g - k = - 2 2m ∇ 2 g - k - 2ρ k V 1 + ρ 2 k g + k + V - 1 + 3ρ 2 k 1 + ρ 2 k V g - k . (8)
Equations ( 7) and ( 8), which are coupled at most by a term of the order of

V [since | V | ≤ |V | and 2ρ k /(1 + ρ 2 k ) ≤ 1]
, allow for perturbative approaches. Note that the functions g + k and g - k have very different behaviors owing to the signs in the lefthand-side terms in Eqs. (7) and (8). Equation ( 8) can be solved to the lowest order in V R /µ with a Green kernel defined in Ref. [6]: we find g 2 . Equation ( 7) cannot be solved using the same method because the perturbation series diverges. Nevertheless, from the solution for g - k , we find that

- k (r) ≃ 2/ǫ k 1+ρ 2 k dr ′ G ξ k (r-r ′ ) V (r ′ )g + k (r ′ ), where ξ k = ξ/ 1 + (kξ)
|g - k /g + k | 2/ǫ k 1+ρ 2 k | V | < | V |/µ ≪ 1.
The coupling term in Eq. ( 7) can thus be neglected to first order in V R /µ and we are left with the closed equation

-( 2 /2m)∇ 2 g + k + V k (r)g + k ≃ ( 2 k 2 /2m)g + k , (9)
where

V k (r) = V (r) - 1 + 4(kξ) 2 1 + 2(kξ) 2 V (r). (10)
Equation ( 9) is formally equivalent to a Schrödinger equation for non-interacting bare particles with energy 2 k 2 /2m, in a random potential V k (r). This mapping allows us to find the localization properties of the BQPs using standard methods for bare particles in 1D, 2D or 3D [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]. However, since V k (r) depends on the wave vector k itself, the localization of the BQPs is dramatically different from that of bare particles as discussed below.

In the remainder of the Letter, we restrict ourselves to the 1D case, for simplicity, but also because AL is expected to be stronger in lower dimensions [2]. The Lyapunov exponent Γ k is a self-averaging quantity in infinite 1D systems, which can be computed in the Born approximation using the phase formalism [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF] (see also Ref. [17]). We get 6) and (10), we find

Γ k = ( √ 2π/8)(2m/ 2 k) 2 C k (2k), where C k (q) is the Fourier transform of the correlation function of V k (z), provided that Γ k ≪ k [16, 17, 27]. Since C k (q) ∝ |V k (q)| 2 , the compo- nent of V k relevant for the calculation of Γ k is V k (2k). From Eqs. (
V k (2k) = S(kξ)V (2k); S(kξ) = 2(kξ) 2 1 + 2(kξ) 2 , (11)
and the Lyapunov exponent of the BQP reads

Γ k = [S(kξ)] 2 γ k ( 12
)
where

γ k = ( √ 2π/32)(V R /µ) 2 (σ R /k 2 ξ 4 )c(2kσ R )
is the Lyapunov exponent for a bare particle with the same wavenumber k [17,[START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF].

Let us summarize the validity conditions of the perturbative approach presented here. It requires: (i) the smoothing solution (5) to be valid (i.e. V R ≪ µ); (ii) the coupling term proportional to g - k in Eq. ( 7) to be negligible (which is valid if V R ≪ µ); and (iii) the phase formalism to be applicable. The latter requires

Γ k ≪ k, i.e. (V R /µ)(σ R /ξ) 1/2 ≪ (kξ) 3/2 [1 + 1/2(kξ) 2 ], which is valid for any k if (V R /µ)(σ R /ξ) 1/2 ≪ 1.
Applying Eq. ( 12) to uncorrelated potentials

[C(z) = 2Dδ(z) with σ R → 0, V R → ∞ and 2D = V 2 R σ R dxc(x) = cst]
, one recovers the formula for Γ k found in Ref. [22]. Our approach generalizes this result to potentials with finiterange correlations, which proves useful since uncorrelated random potentials are usually crude approximations of realistic disorder, for which σ R can be significantly large. We show below that if ξ σ R , as e.g. in the experiments of Refs. [10,11,12,13,14], the behavior of Γ k versus k is dramatically affected by the finite-range correlations of the disorder.

Let us discuss the physical content of Eqs. ( 9) and ( 10). According to Eqs. ( 3) and ( 4), the properties of the BQPs are determined by both the bare random potential V and the BEC density n c in a non-trivial way. Equation ( 9) makes their roles more transparent. As the occurrence of the smoothed potential V (z) in Eq. ( 10) is reminiscent of the presence of the meanfield interaction gn c in the BdGEs (3) and ( 4), it appears that the random potential V k (z) results from the screening of the random potential V (z) by the BEC density background [5]. More precisely, the expression (11) for the Fourier component V k (2k), relevant for the Lyapunov exponent of a BQP, shows that the screening strength depends on the wavenumber k. In the free-particle regime (k ≫ 1/ξ), we find that the Lyapunov exponent of a BQP equals that of a bare particle with the same wavenumber (Γ k ≃ γ k ), as expected. In the phonon regime (k ≪ 1/ξ), the disorder is strongly screened and we find Γ k ≪ γ k , as in models of elastic media [START_REF] Ishii | [END_REF].

Here, the localization of a BQP is strongly suppressed by the repulsive atom-atom interactions, as compared to a bare particle in the same bare potential. These findings agree with and generalize the results obtained from the transfer matrix method, which applies to potentials made of a 1D random series of δ-scatterers [22].

Our approach applies to any weak random potential with a finite correlation length. We now further examine the case of 1D speckle potentials used in quantum gases [10,11,12,13,14]. Inserting the corresponding reduced correlation function,

c(κ) = π/2(1 -κ/2)Θ(1 -κ/2)
where Θ is the Heaviside function [17], into Eq. ( 12), we find

Γ k = π 8 V R µ 2 σ R k 2 (1 -kσ R ) [1 + 2(kξ) 2 ] 2 Θ(1 -kσ R ) (13)
which is plotted in Fig. 1. To test our general approach on the basis of this example, we have performed numerical calculations using a direct integration of the BdGEs (3) 13) as shown in Fig. 2. These results validate our approach. It should be noted, however, that our numerical calculations return BQP wave functions that can be strongly localized for very small momenta k. This will be discussed in more details in a future publication [29].

Of special interest are the maxima of Γ k , which denote a maximum localization of the BQPs. It is straightforward to show that, for a fixed set of parameters (V R /µ, ξ, σ R ), Γ k is non-monotonic and has a single maximum, k max , in the range [0, 1/σ R ] (see Fig. 2). This contrasts with the case of bare particles, for which the Lyapunov exponent γ k decreases monotonically as a function of k, provided that c(2kσ R ) decreases versus k (which is valid for a broad class of random potentials [START_REF] Lifshits | Introduction to the Theory of Disordered Systems[END_REF]). The existence of a localization maximum with respect to the wavenumber k is thus specific to the BQPs and results from the strong screening of the disorder in the phonon regime. In general, the value of k max , plotted in the inset of Fig. 2 versus the correlation length of the disorder, depends on both ξ and σ R . For σ R ≪ ξ, we find

k max ≃ 1 √ 2ξ 1 -σR/ξ 2 √
2 , so that the localization is maximum near the crossover between the phonon and the free-particle regimes as for uncorrelated potentials [22]. For σ R ≫ ξ however, we find k max ≃ 2/3σ R , so that k max is no longer determined by the healing length but rather by the correlation length of the disorder, and lies deep in the phonon regime. For k > 1/σ R , Γ k vanishes. This defines an effective mobility edge due to longrange correlations in speckle potentials, as for bare particles [17,[START_REF]Beyond the Born approximation, Γ k and γ k may not strictly vanish for k > 1/σR[END_REF].

Finally, let us determine the absolute localization maximum. The Lyapunov exponent Γ k decreases monotonically versus ξ and V R /µ. However, for fixed values of V R /µ and ξ, Γ k has a maximum at σ R = 3/2 ξ and k = ξ -1 / √ 6 (see Fig. 1) and we find the corresponding localization length (L max = 1/Γ max ):

L max (ξ) = (512 √ 6/9π)(µ/V R ) 2 ξ. (14)
At the localization maximum, we have σ R ∼ ξ, so that the disorder cannot be modeled by an uncorrelated potential, and the long-range correlations must be accounted for as in our approach. For σ R = 0.3µm [14] and V R = 0.2µ, we find L max ≃ 280µm, which can be smaller than the system size in disordered, ultracold gases [11,14,[START_REF] Hulet | [END_REF].

In conclusion, we have presented a general treatment for the AL of BQPs in an interacting BEC subjected to a random potential with finite-range correlations. We have calculated the Lyapunov exponents for a 1D speckle potential and we have shown that the localization is strongest when σ R ∼ ξ. We have found that the localization length can be smaller than the size of the BEC for experimentally accessible parameters. We expect that the AL of BQPs could be observed directly, for instance as a broadening of the resonance lines in Bragg spectroscopy, a well mastered technique in gaseous BECs [START_REF] Stamper-Kurn | [END_REF].

Figure 1 :

 1 Figure 1: (color online) Density plot of the Lyapunov exponent of the BQPs for a 1D speckle potential.

Figure 2 :

 2 Figure 2: (color online) Lyapunov exponent of the BQPs for a 1D speckle potential. The lines correspond to Eq. (13) and the points to numerical results for ξ = 1.25 × 10 -5 L, VR = 0.075µ and σR = 3/2ξ (solid red line), σR = 3.7ξ (dashed blue line), σR = 0.4ξ (dash-dotted green line). Inset: localization maximum versus the ratio of the correlation length of the disorder to the healing length of the BEC (the dash-dotted and dashed lines correspond to the limits σR ≪ ξ and σR ≫ ξ respectively; see text).

We thank M. Lewenstein, G. Shlyapnikov, S. Stringari, and W. Zwerger for stimulating discussions during the Workshop on Quantum Gases at the Institut Henri Poincaré -Centre Emile Borel. This work was supported by the French DGA, MENRT and ANR, and the ESF QUDEDIS program. The Atom Optics group at LCFIO is a member of the Institut Francilien de Recherche sur les Atomes Froids (IFRAF).


Philip Feinsilver

René Schott

Finite-Dimensional Calculus

Keywords: inversion formulas, orthofermions, Heisenberg-Weyl algebra, polynomial systems

We develop finite-dimensional calculus using matrices.

The truncated Heisenberg-Weyl algebra is called a TAA algebra after Tekin, Aydin, and Arik who formulated it in terms of orthofermions. The matrix approach is used to implement our method of polynomial inversion in one-variable and multivariable settings. Here we establish notations, present some algebraic developments, and discuss the univariate case in detail.

Introduction

The TAA algebra is generated by an operator a and its adjoint a * . Let N = a * a. Then the defining commutation rule is [a, N] = a. Taking adjoints gives the complementary rule [N, a * ] = a * . In other words, we assume that N and a generate the two-dimensional Lie algebra of the affine group. Writing it out aa * a -a * aa = a [START_REF] Feinsilver | The affine group and generalized Gegenbauer polynomials[END_REF] we see that the commutation rule aa * -a * a = 1 of the Heisenberg-Weyl algebra has been multiplied by a on the right. This modification is enough to yield finite-dimensional representations, including the truncated HW-algebra given by the operators X =multiplication by x and D = d/dx acting on polynomials of a given bounded degree. Writing matrices for these operators we will see that they obey equation [START_REF] Feinsilver | The affine group and generalized Gegenbauer polynomials[END_REF].

This article may be thought of as realizing Rota's idea of "Finite Operator Calculus". There is work of P. R. Vein along similar lines [START_REF] Vein | Identities among certain triangular matrices[END_REF][START_REF] Vein | Matrices which generate families of polynomials and associated infinite series[END_REF]. The main feature here is that, in fact, the operator calculus is done on finite-dimensional spaces and can be carried out explicitly using matrices. The approach in this paper is based on algebraic properties of the operators and includes indications for the multivariable case. The one-variable case is dual to that presented in [START_REF] Feinsilver | Inversion of analytic functions via canonical polynomials: a matrix approach[END_REF].

Orthofermion formulation

Orthofermions with regard to connections with the HW algebra have been studied in [START_REF] Tekin | A new interpretation of the orthofermion algebra[END_REF]. Start with a set of operators {c 1 , . . . , c p }, with p a positive integer and form the star-algebra generated by the {c i } modulo the following relations

c i c j = 0 c i c * j + δ ij p k=1 c * k c k = δ ij 1 (2)
where 1 is the identity operator. Setting Π = 1 -p k=1 c * k c k (as in [START_REF] Mostafazadeh | On the Representation Theory of Orthofermions and Orthosupersymmetric Realization of Parasupersymmetry and Fractional Supersymmetry[END_REF]) we can write this last relation as c i c * j = δ ij Π where one readily shows from the defining relations (2) that Π 2 = Π, i.e., Π is a projection as suggested by the notation.

From the defining relations, we see that Πc k = c k and from the second relation of eq. ( 2) follows the useful relation

c i c * j c k = δ ij c k (3)
Within the orthofermion algebra, following [START_REF] Tekin | A new interpretation of the orthofermion algebra[END_REF], modifying slightly their formulation, we set

a = c 1 + p k=2 k c * k-1 c k a † = c * 1 + p k=2 c * k c k-1
Using equation (3), we get

aa † -a † a = 1 -(p + 1) c * p c p
which then yields the relation corresponding to equation (1) of the TAA algebra.

One-variable calculus with matrices

Restricting the differentiation operator to the finite-dimensional space of polynomials of degree less than or equal to p is no problem. Use the standard basis {1, x, x 2 , . . . , x p }. For p = 4, we have

D =       
0 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4 0 0 0 0 0

      
with the extension to general p following the same pattern. However, multiplication by x must be cut off. If we define X x i = x i+1 for i < p and X x p = 0, we no longer have the relation DX -XD = 1. But we still have DXD -XDD = D as for the TAA algebra. The matrix of X has the form, for p = 4, e.g.,

X =        0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0       
Note that Xp+1 = 0. To keep in line with the powers of x, we label the basis elements starting from 0. So let e k denote the column vector with the only nonzero entry equal to 1 in the (k + 1) st position. The vacuum state is Ω = e 0 , satisfying DΩ = 0. And Xk Ω = e k , for 1 ≤ k ≤ p. As expected, these are raising and lowering operators satisfying

Xe k = e k+1 θ kp , De k = k e k-1
where θ ij = 1 if i < j, zero otherwise.

With the inner product e n , e m = δ nm n! , we indeed have D * = X.

Let E ij denote the standard unit matrices with all but one entry equal to zero, (E ij ) kl = δ ik δ jl , 1 ≤ i, j, k, l ≤ p+1. The connection with orthofermions is given by the (p + 1) × (p + 1) matrix realization ĉi = E 1 i+1

for 1 ≤ i ≤ p. We have the orthofermion relations and particularly for this realization ĉ * i ĉj = E i+1 j+1 Note that Π = E 11 and that the star-algebra generated by the ĉi is the full matrix algebra.

As long as X never multiplies the power x p , the matrix implementation agrees with usual calculus. The TAA relation formulates this algebraically.

The following theorem shows that D and X not only do not generate a Heisenberg algebra, but, in fact, are as far as possible from doing so. Theorem 3.1 For p > 0, let D and X be (p+1)×(p+1) matrices defined by

D = p k=1 k E k k+1 , X = p k=1 E k+1 k . Then the Lie algebra generated by { X, D} is sl(p + 1) . Proof: For convenience set n = p + 1. First we have H = [ D, X] = -p E nn + p k=1 E kk Set ξ 1 = X, η 1 = D, and H 1 = H. For 2 ≤ k ≤ n, let ξ k = H( ←-- ad X) k , and η k = (ad D) k H, where (ad A)B = [A, B] and A( ←-- ad B) = [A, B].
Then it is easily checked by induction that

ξ k = -n E n n-k+1 and η k = a k ξ † k
for nonzero constants a k , the † denoting matrix transpose. Thus, we have

E in and E ni for 1 ≤ i ≤ p. Noting that [E in , E nj ] = E ij if i = j,
we have all of the off-diagonal E ′ s. And

H k = [η k , ξ k ] = -na k (E n-k+1 n-k+1 -E nn )
fill out the Cartan elements of sl(n).

Examples

Here we look at some important operators for the case p = 4.

Example 3.2

The number operator is XD. We have

X D =             0 0 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4            
And, in general, we calculate

X D = p k=1 k E k+1 k+1
as we want.

Example 3.3 The Ornstein-Uhlenbeck operator , XD -D 2 , has the form

            0 0 0 0 0 0 1 -2 0 0 0 0 2 -6 0 0 0 0 3 -12 0 0 0 0 4             Example 3.4
The translation operator, e D looks like

            1 1 1 1 1 0 1 2 3 4 0 0 1 3 6 0 0 0 1 4 0 0 0 0 1            
with columns given by binomial coefficients, corresponding to the action x → x + 1 on the basis polynomials x i .

Example 3.5 The Gegenbauer operator , see, e.g. [START_REF] Feinsilver | The affine group and generalized Gegenbauer polynomials[END_REF], (XD + α) 2 -D 2 , has the matrix form

            α 2 0 -2 0 0 0 (1 + α) 2 0 -6 0 0 0 (2 + α) 2 0 - 12
0 0 0 (3 + α) 2 0 0 0 0 0 (4 + α) 2            
with corresponding Gegenbauer polynomials as eigenfunctions.

Canonical variables and polynomials

Basic to our approach is the use of canonical variables which are functions of X and D obeying the HW relations on an infinite-dimensional space, which restricts to the TAA relation on spaces of polynomials in x of a given bounded degree.

Let us review the basic construction and notations. We will discuss the general, multivariable, case.

Note. We use our usual convention of summing over repeated Greek indices, irrespective of position.

Given V : C N → C N , V (z) = (V 1 (z 1 , . . . , z N ), . . . , V N (z 1 , . . . , z N ))
holomorphic in a neighborhood of the origin, satisfying V (0) = 0, we construct a corresponding abelian family of dual vector fields. Corresponding to the operators X i of multiplication by x i , we have the partial differentiation operators, D i . In this context, a function of x = (x 1 , . . . , x N ), f (x), is identified with f (X)1, the operator of multiplication by f (X) acting on the vacuum state 1, with D i 1 = 0, for all 1 ≤ i ≤ N. We define operators V (D) = (V 1 (D 1 , . . . , D N ), . . . , V N (D 1 , . . . , D N )). These are our canonical lowering operators, corresponding to differentiation.

Denoting the Jacobian ∂V i ∂z j by V ′ (z), let W (z) = (V ′ (z)) -1 , be the inverse (matrix inverse) Jacobian. Then the boson commutation relations extend to

[V (D), X i ] = ∂V ∂D i . Now define the operators Y i = X µ W µi (D)
These are our canonical raising operators, corresponding to multiplication by X i . The canonical system of raising and lowering operators {Y j }, {V i },

1 ≤ i, j ≤ N indeed satisfy [V i , Y j ] = δ ij 1.
The essential feature, which has to be checked, is that, [Y i , Y j ] = [V i , V j ] = 0. Notice that exchanging D with X is a formal Fourier transformation and turns the variables Y i into the vector fields Ỹi = W (x) µi ∂/∂x µ . Thus, the Y i are dual vector fields.

Notation. We complement the standard notations used along with V and W , letting U denote the inverse function to V . I.e., U

• V = V • U = id. Explicitly: U(V (z)) = z.
Observe that since W = V ′-1 , we have W (z) = U ′ (V (z)). In other words, converting from z to V acting on functions of the canonical variables Y i , we have

X = Y U ′ (V ) -1 .
The main formula (cf. [3, p. 185, eq. ( 1)]) is

exp(α µ Y µ ) 1 = exp(x µ U µ (α)) = n≥0 α n n! y n (x)
the n denoting a multi-index (n 1 , . . . , n N ). This expansion defines the canonical polynomials:

y n (x) = Y n 1.

Multivariable calculus with matrices

Here we extend Section 3 to N variables. For matrices, A, B, the tensor product A ⊗ B denotes the Kronecker product of the two matrices. That is, if A is n × n, and B is m × m, then A ⊗ B is nm × nm with entries formed by replacing each entry a ij in A with the block matrix a ij B. For products of more than two matrices, we conventionally associate to the left.

For a fixed p, we have (p + 1) × (p + 1) matrices D and X. Let I denote the (p + 1) × (p + 1) identity matrix. Then we set

Dj = I ⊗ I ⊗ • • • ⊗ D ⊗ I • • • ⊗ I ( D in the j th spot) Xj = I ⊗ I ⊗ • • • ⊗ X ⊗ I • • • ⊗ I ( X in the j th spot)
Then Dj and Xj will satisfy the TAA relations while [ Dj , Xi ] = [ Xj , Xi ] = [ Dj , Di ] = 0 for i = j.

Canonical calculus with matrices

First consider the case N = 1. We have a function V (z) analytic in a neighborhood of the origin in C, normalized to V (0) = 0, V ′ (0) = 0. Let W (z) = 1/V ′ (z) have the Taylor expansion

W (z) = w 0 + w 1 z + • • • + w k z k + • • • The corresponding canonical variable is Y = XW (D), satisfying [V (D), Y ] = 1.
The canonical basis polynomials are y n (x) = Y n 1, n ≥ 0. Fix the order p. Let Ŵ = W ( D). Then we employ the algebra generated by the operators V = V ( D) and Ŷ = X Ŵ . Note, e.g., that since Dp+1 = 0, the operators V and Ŵ are polynomials in D. Similarly, since Xp+1 = 0, the polynomials y n ( X) are truncated if n > p. However, for n ≤ p, the correspondence between the polynomials y n (x) and vectors ŷn = y n ( X)e 0 is exact. Namely, the vector ŷn gives the coefficients of the polynomial y n (x). The reason this works is that up to order p, the operator X never acts on a power of x greater than p.

Examples

Example 6.1 A basic example is given by

V (z) = e z -1 , U(v) = log(1 + v) so W (z) = e -z , Y = Xe -D .
We can easily calculate

y n (x) = x(x -1) • • • (x -n + 1) . For p = 4, we get Ŷ =             0 0 0 0 0 1 -1 1 -1 1 0 1 -2 3 -4 0 0 1 -3 6 0 0 0 1 -4             and Ŷ 2 =             0 0 0 0 0 -1 2 -4 8 -15 1 -3 8 -20 43 0 1 -5 18 -46      0 0 0 0 0 24 -
           
with the first column giving the coefficients of the corresponding polynomial y n , where we can see the truncation beginning in this last.

With

U(v) = log(1 + v), the relation X = Y U ′ (V ) -1 reads X = Y + Y V or xy n = y n+1 + ny n yielding the recurrence y n+1 = (x -n)y n .
Example 6.2 Another interesting example is the Gaussian with drift,

V (z) = αz -z 2 /2 , U(v) = α - √ α 2 -2v
the minus sign taken in U(v) to have U(0) = 0. Then W (z) = 1 α -z , and

Ŷ =             0 0 0 0 0 α -1 α -2 2 α -3 6 α -4 24 α -5 0 α -1 2 α -2 6 α -3 24 α -4 0 0 α -1 3 α -2 12 α -3 0 0 0 α -1 4 α -2            
Powers of Ŷ yield the canonical polynomials, the first few of which are

y 1 = x α y 2 = x α 3 + x 2 α 2 y 3 = 3 x α 5 + 3 x 2 α 4 + x 3 α 3 y 4 = 15 x α 7 + 15
x 2 α 6 + 6

x 3 α 5 +

x 4 α 4 y 5 = 105 x α 9 + 105

x 2 α 8 + 45

x 3 α 7 + 10

x 4 α 6 + x 5 α 5
These are a scaled variation of Bessel polynomials.

In this case

U ′ (V ) -1 = α 1 - 2V α 2 1/2
. Thus, expanding and rearranging the relation

X = Y U ′ (V ) -1 , αY = X+αY V α 2 + 1 2 V 2 α 4 + 1 2 V 3 α 6 + 5 8 V 4 α 8 + 7 8
V 5 α 10 + 21 16

V 6 α 12 + 33 16

V 7 α 14 + . . . which translates to α y n+1 = xy n + n α y n + n(n -1) 2α 3 y n-1 + n(n -1)(n -2) 2α 5 y n-2 + . . . = xy n + n α y n + n k=2 n k (2k -3)!! α 2k-1 y n-k+1
Example 6.3 For our final example in this section, we mention the reflected LambertW function, which we denote W to avoid confusion with our W . Take V (z) = ze -z , see [4, p. 110]. Then U(v) = -W(-v). We find Y = Xe D (I -D) -1 and with p = 7 the corresponding matrix

Ŷ =                       0
                     
One can show that y n = x(x + n) n-1 and that the relation X = Y U ′ (V ) -1 leads to the recurrence

y n+1 = (x + 2n)y n + n k=2 n k (k -1) k-1 y n-k+1 .

Matrix expansions

Using the unit matrices E ij of size (p + 1) × (p + 1), we can write formulas for the basic operators. We have used in Theorem 3. For N > 1, the matrices for D j and X j provide the operators Ŷj = Xµ W µi ( D) as matrices. Repeated multiplication on the vacuum vector e 0 yields exactly the coefficients of the polynomials y n up to order p. Here, order p means that no variable x i appears to a power higher than p.

Summary and prospects

In this article, the one-variable case of the matrix approach has been presented in some detail, along with the basic theory for the multivariate case.

The TAA algebra conveniently replaces the HW algebra in the finite-dimensional setting. The connection with orthofermions is interesting and clarifies the underlying structure.

In Part II of this series of articles, we will look at multivariate systems from the point of view of dual vector fields. We will present and use corresponding matrix formulations. In Part III, we expect to have some C and/or Java code implementing the techniques of Parts I and II. For Part IV, the plan is to develop useful rigorous estimates for establishing the order of truncation required for applications.

E

  j E k k+j with (k) j = k(k + 1) • • • (k + j -1)denoting the rising factorial. Multiplying by X gives Ŷ = 1≤k≤p 1≤k+j≤p+1 k+1 k+j (k) j w j


Clemence Petit

Sylvain Meille
email: sylvain.meille@insa-lyon.fr

Eric Maire

Laurent Gremillard

Jérôme Adrien

Grace Y Lau

Antoni P Tomsia

Clémence Petit

Fracture behavior of robocast HA/beta-TCP scaffolds studied by X-ray tomography and finite element modeling

Keywords: Calcium phosphate, porous ceramics, X-ray tomography, mechanical behavior, finite element modeling

come   L'archive ouverte pluridisciplinaire

Introduction

Calcium phosphate ceramics, especially hydroxyapatite (HA) and β-tricalcium phosphate (β-TCP) have been widely used for bone substitution [START_REF] Bohner | Calcium orthophosphates in medicine: from ceramics to calcium phosphate cements[END_REF][START_REF] Legeros | Biphasic calcium phosphate bioceramics: preparation, properties and applications[END_REF][START_REF] Bohner | Technological issues for the development of more efficient calcium phosphate bone cements: A critical assessment[END_REF][START_REF] Stevens | Biomaterials for bone tissue engineering[END_REF] because they exhibit a unique combination of properties: biocompatibility, osteoconductivity, osteointegration, bioresorption, and bioactivity [START_REF] Wagoner Johnson | A review of the mechanical behavior of CaP and CaP/polymer composites for applications in bone replacement and repair[END_REF]. HA/β-TCP composites (commonly named biphasic calcium phosphates [BCP]) have been developed [START_REF] Daculsi | Biphasic calcium phosphate concept applied to artificial bone, implant coating and injectable bone substitute[END_REF] with the aim of optimizing material in terms of mechanical properties and resorption in the human body. An interconnected porosity is necessary to promote bone regeneration and vascularization around and inside the scaffolds, as well as to allow circulation of fluids in the scaffold. Several studies have shown that scaffolds need interconnected macropores with a size higher than 100 µm and micropores with a size between 10 and 100 µm to improve bone regeneration [START_REF] Bignon | Effect of micro-and macroporosity of bone substitutes on their mechanical properties and cellular response[END_REF][START_REF] Hing | Microporosity enhances bioactivity of synthetic bone graft substitutes[END_REF][START_REF] Hulbert | Tissue reaction to three ceramics of porous and non-porous structures[END_REF][START_REF] Liu | Influence of porosity and pore size on the compressive strength of porous hydroxyapatite ceramic[END_REF]. Various fabrication techniques have been used to fabricate porous calcium phosphate samples: the replica technique [START_REF] Jo | Highly porous hydroxyapatite scaffolds with elongated pores using stretched polymeric sponges as novel template[END_REF], sacrificial templates [START_REF] Chevalier | Fabrication of porous substrates: A review of processes using pore forming agents in the biomaterial fields[END_REF], and direct foaming [START_REF] Ahn | Production of highly porous triphasic calcium phosphate scaffolds with excellent in vitro bioactivity using vacuum-assisted foaming of ceramic suspension (VFC) technique[END_REF]. Additive manufacturing (AM) has gained interest because it can provide samples with complex shapes and an intricate internal architecture, allowing control of mechanical properties [START_REF] Chartier | Stereolithography of structural complex ceramic parts[END_REF][START_REF] Farzadi | Effect of layer printing delay on mechanical properties and dimensional accuracy of 3D printed porous prototypes in bone tissue engineering[END_REF]. Among the AM techniques available to process ceramics, robocasting is a robotic deposition of a ceramic ink (viscous slurry) through a nozzle following a computer-aided model. Robocasting can be distinguished from the other AM processes by the use of aqueous ceramic inks with a low amount of organic additives (typically less than 1wt%) [START_REF] Saiz | Preparation of porous hydroxyapatite scaffolds[END_REF][START_REF] Dellinger | Robotic deposition of model hydroxyapatite scaffolds with multiple architectures and multi-scale porosity for bone tissue engineering[END_REF]. The inks are also highly viscous thus with a high solid loading. These features facilitate the densification during the heat treatment. The ceramic ink has to flow through the nozzle during the extrusion step and set after deposition [START_REF] Dellinger | Robotic deposition of model hydroxyapatite scaffolds with multiple architectures and multi-scale porosity for bone tissue engineering[END_REF]. The ink's characteristics (e.g. high solid loading, concentration of organic additives such as dispersant) must be carefully controlled to avoid some specific defects. For example, if the solid loading is too low, the shrinkage will be too important during drying and cracks will appear [START_REF] Dellinger | Robotic deposition of model hydroxyapatite scaffolds with multiple architectures and multi-scale porosity for bone tissue engineering[END_REF][START_REF] Guo | Aggregation effects on the compressive flow properties and drying behavior of colloidal silica suspensions[END_REF]. After deposition, there is also a risk of bending of the rods [START_REF] Franco | Direct write assembly of calcium phosphate scaffolds using a water-based hydrogel[END_REF]. A too high solid loading will impose a high pressure for extrusion. This can cause an inhomogeneous re-distribution of liquid and ceramic particles into the ink [START_REF] Guo | Aggregation effects on the compressive flow properties and drying behavior of colloidal silica suspensions[END_REF][START_REF] Smay | Colloidal inks for directed assembly of 3-D periodic structures[END_REF]. Thus the solid loading will be different throughout the structure. Thus, several studies have been dedicated to the optimization of this fabrication process for ceramics and the formulation of ceramic inks [START_REF] Saiz | Preparation of porous hydroxyapatite scaffolds[END_REF][START_REF] Dellinger | Robotic deposition of model hydroxyapatite scaffolds with multiple architectures and multi-scale porosity for bone tissue engineering[END_REF][START_REF] Duoss | Sol-gel inks for direct-write assembly of functional oxides[END_REF]. For example, Houmard et al. [START_REF] Houmard | On structural, mechanical and biodegradation properties of HA-TCP robocast scaffolds[END_REF] fabricated BCP samples with a wide range of volume fraction and pore size. Other authors focused more on the mechanical properties of scaffolds processed by robocasting, such as Miranda et al., who reported the propagation of long cracks along the compression direction for HA and β-TCP samples [START_REF] Miranda | Mechanical properties of calcium phosphate fabricated by robocasting[END_REF][START_REF] Miranda | Fracture modes under uniaxial compression in hydroxyapatite scaffolds fabricated by robocasting[END_REF]. The aim of this research is to perform a fundamental and qualitative study of the fracture process in uniaxial compression of BCP scaffolds by means of X-ray tomography. This nondestructive technique provides three-dimensional (3D) images of the samples. It has been extensively used to characterize highly porous materials [START_REF] Maire | X-Ray tomography applied to the characterization of highly porous materials[END_REF]. In or ex situ mechanical tests can also be performed by adding special rigs inside the tomograph. We have recently used this approach to characterize the uniaxial compression behavior of a porous β-TCP sample with a stochastic porosity generated by a sacrificial template method [START_REF] Petit | Mechanical behavior of a β-TCP ceramic with a random porosity: study of the fracture path with X-ray tomography[END_REF]. Additionally, the tomographic images can be used as input files for finite element (FE) modeling [START_REF] Petit | Cellular solids studied by X-ray tomography and finite element modeling -a review[END_REF]. Meshing the 3D images allows performing calculations based on the real architecture of the specimens. In the present study, porous BCP samples were processed via robocasting and scanned by Xray tomography. Model materials were designed by adding artificial defects in the struts thanks to pore formers. Our aim was to generate defects having a controlled size in the structure. The location of these defects in the initial microstructure was determined by highresolution scans obtained with the so-called "local tomography" mode. Then an ex situ compression test was performed to follow the crack propagation by standard low resolution X-ray tomography. Finally, high-resolution images of the initial samples were processed to create an FE model and to check the influence of artificial defects in the struts on the mechanical behavior of the scaffolds.

Materials and methods

Sample fabrication

Cellular samples with a solid phase containing a volume fraction of 80% HA and 20% β-TCP were fabricated following a protocol thoroughly described by Miranda et al. [START_REF] Miranda | Sintering and robocasting of β-tricalcium phosphate scaffolds for orthopaedic applications[END_REF]. A 20 vol.% Pluronic F127 solution in water was prepared and cooled to approximately 4°C. Adequate quantities of HA (Science Applications Industries, France) and of β-TCP (Fluka, Switzerland) were mixed with this solution to form an aqueous paste with a solid loading of 30 vol.%, fluid at low temperature and solid at room temperature thanks to the addition of Pluronic. Two vol.% of polymethacrylate (PMMA) spheres (Diakon CLG356, Lucite International, France) were added as porogens to create artificial defects inside the solid struts. The PMMA beads were previously sieved and the spheres with a diameter from 45 to 100 µm were selected. This particular size range was chosen because it allows having controlled defects in the rods which can be detected by X-ray tomography with resolution of few microns, as used in this work. This size range can also be distinguished from the residual micropores with a size of few microns [START_REF] Franco | Direct write assembly of calcium phosphate scaffolds using a water-based hydrogel[END_REF][START_REF] Miranda | Mechanical properties of calcium phosphate fabricated by robocasting[END_REF]. The ceramic ink and the PMMA beads were mixed in a jar. The mixture of ceramic ink and PMMA beads was introduced in the conical nozzle of the robocasting device (3-D Inks, Stillwater, USA) controlled by a software program (Robocad 3.0, 3-D Inks, Stillwater, USA). The nozzle had a diameter of 150 µm and was displaced to deposit the ceramic ink on a plane with a constant deposition speed of 20 mm.s -1 . The spacing between two rods (from rod center to rod center) was fixed to 300 µm. After finishing the deposition of several parallel rods, the nozzle was displaced in the vertical direction and another layer of rods was deposited in a direction perpendicular to the previously deposited rods. The process was continued until all the layers were deposited. The deposition was made in an oil bath to avoid collapse of the layers during the process. After deposition, the samples were air-dried for 3 days. Finally, a two-step thermal treatment took place. The first step consisted of a debinding treatment at 550°C during 2 hours to remove the PMMA spheres (artificial defects) and organic binders. The sintering occurred during the second step at 1100°C during 4 hours. During sintering, it was expected that the large artificial defects would remain unsintered while the small pores between solid particles would close. The heating rate was 1°C/min until 550°C and 2°C/min from 550 to 1100°C. The cooling rate was 5°C/min. The overall approach described in this paper (i.e., ex situ compression test, image processing steps, and FE modeling) was applied to one sample. This is justified by the difficulty to reproduce this timeconsuming approach to a large number of samples. For this reason, our goal was to develop a qualitative approach about the influence of the defects on the macroscopic behavior. The results presented in this paper concern a sample with dimensions of 3.58 mm (in direction 3, which was the layer-by-layer direction during fabrication and the loading direction during compression), and 4.62 mm and 4.58 mm (length and width, directions 1 and 2, respectively).

Structural characterization

Different characterization methods were used. (1) Geometric relative density was calculated from the weight and the dimensions of the sample, considering a density of 3.16 g.cm -3 and 3.14 g.cm -3 for the fully dense HA and β-TCP, respectively. (2) One of the faces another sample was carefully polished to 1 µm with diamond paste to observe the solid phase. After polishing, the sample was cleaned by ultrasonic bath to remove the polishing debris in the pores. After polishing, the sample was etched by a thermal treatment at 1000°C over 1 hour. Then the polished and etched face was observed with scanning electron microscope (SEM), (ZEISS, SUPRA55 VP, Germany) at a voltage of 1 kV without metallization. (3) The phase composition was identified by X-ray diffraction (XRD, Siemens D-500, Munich, Germany) using Cu Kα radiation (λ=0.15406 nm). The acquisitions were carried out at a scanning rate of 1.8°/min in the 2θ range of 20-80°. The HAp, β-TCP and α-TCP were searched by comparison with the JCPDS cards 09-0432, 09-0169 and 09-348 respectively. The weight ratio of the phases was also determined. (4) Mercury intrusion porosimetry (Auto Pore IV, Micromeritics, USA) was used to determine the interconnection pore size distribution by intrusion of liquid mercury into the pores with a maximal pressure of 400 MPa.

X-ray tomography

A laboratory X-ray tomograph (v/Tome/X, Phoenix, Germany) described in detail by Buffière et al. [START_REF] Buffière | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF] was used to scan the whole sample in its initial state at a "low" resolution of 15 µm. The tomograph was operated at a voltage of 80 kV and an intensity of 280 µA with no filtering on incident X-rays. The 3D final images were then obtained by a standard filtered back projection algorithm implemented in the software coupled to the tomography [START_REF] Kak | Principles of computerized tomographic imaging[END_REF]. To observe the finer details of the microstructure in this initial state, local tomography was used [START_REF] Zhang | Local tomography study of the fracture of an ERG metal foam[END_REF]. For this, the sample was placed near the X-ray source; the detector's field of view was then smaller than the size of the sample. This resulted in only a part of the sample (the cylinder of reconstruction) being irradiated by X-rays under all the viewing angles. Some authors have proved that reconstruction can be performed in this local mode [START_REF] Youssef | Finite element modelling of the actual structure of cellular materials determined by X-ray tomography[END_REF][START_REF] Toda | Quantitative assessment of microstructure and its effects on compression behavior of aluminum foams via high-resolution synchrotron X-ray tomography[END_REF]. It is particularly well adapted to porous samples due to their low density [START_REF] Zhang | Local tomography study of the fracture of an ERG metal foam[END_REF]. The reconstructed tomogram has a lower signal to noise ratio than in conventional tomography, but a high resolution (voxel size of 3 µm) reconstruction of the irradiated part of the sample can be obtained. To scan another part of the sample, it was necessary only to physically displace the specimen above the rotating stage. Finally, we gathered eight high-resolution images, each representing a part of the sample, to have a representation of the whole sample. To complete the characterization by X-ray tomography in the initial state, small pieces of rods were additionally scanned with a second higher resolution tomograph described by Etiemble et al. [START_REF] Etiemble | 3D morphological analysis of copper foams as current collectors for Li-ion batteries by means of X-ray tomography[END_REF] (EasyTom Nano, Rx Solutions, France) with a voxel size of 0.7 µm. The rods were scanned with a tension of 60 kV and an intensity of 150 µA without using local tomography. An ex situ uniaxial compression test was finally performed to follow the crack initiation and propagation in the sample at low resolution (voxel size of 15 µm). The sample was incrementally loaded outside of the tomograph thanks to a BOSE Electroforce 3200 compression machine (Prairie Valley, USA) equipped with a 200 N load cell and a capacitive displacement sensor, and scanned in the tomograph after the application of different levels of strain. To test the samples at higher loads (above 200 N), an Instron 8502 (High Wycombe, UK) universal hydraulic machine equipped with a 5000 N load cell and a linear variable differential transformer (LVDT) displacement sensor was used. The test was displacement controlled with a speed of 0.03 mm.min -1 and was interrupted several times to scan the sample. To avoid contact damage, the two faces of the sample in contact with the loading platens were embedded with paraffin wax for the two testing conditions. Additionally, during the first loading, the sample was unloaded and reloaded six times to calculate the Young's modulus from the unloading slopes.

Image processing procedure

Image processing of the 3D reconstructed tomograms was performed using ImageJ software [35]. The 8 high-resolution images were first processed to improve the image quality through adjustment of contrast and brightness and the application of a median filter over neighbors of one voxel (Fig. 1a). This figure shows that the artificial defects are indeed present in the sintered rods. Then the images were thresholded to separate the macropores, the ceramic phase, and the artificial defects. Each of these phases had a different grey level (Fig. 1b): 0 for the macro-pores, 128 for the artificial defects, and 255 for the ceramic phase. At this stage, it could have been possible to stitch the thresholded images together to obtain a high-resolution image of the whole sample. However, the result would have been a large image which could not be easily used for finite element modelling. It appears as more relevant to generate an image of the initial sample at a low resolution (i.e., the one used for the ex situ compression test: 15 µm) and containing the information brought by the high-resolution images (i.e., the small defects in the solid struts). For this purpose, these high-resolution thresholded images were scaled to change their voxel size from 3 to 15 µm (Fig. 1c). The scaling operation preserved the segmentation made at high resolution but blurred the image (as shown in Fig. 1c). To avoid this blurring effect, a second thresholding was performed on the scaled images, using a conservative value of the grey level to make sure that no artificial defect was lost during this step (Fig. 1d). The smallest artificial defects have a typical diameter of 3 voxels at the low-resolution scale. Therefore, they could possibly be lost without this precaution. It must be stated, however, that this procedure clearly results in a non-physical slight increase in the size and number of the artificial defects. Finally, the 8 images were stitched two by two until the image of the whole sample was created (not shown in Fig. 1). The final image of the whole sample contained the information of the presence of the artificial defects in different locations inside the sample, rather than actually describing their morphology in detail. This image was used as an input for the FE model (see section 2.7).

Quantitative analysis of the image

The high-resolution thresholded images (i.e., those with a voxel size of 3 µm) were used to perform different morphological characterizations, whose measurement principles are explained in detail by Maire et al. [START_REF] Maire | Characterization of the morphology of cellular ceramics by 3D image processing of X-ray tomography[END_REF]. It was performed by applying plugins (i.e., Java programs implemented in ImageJ) to the images. A first plugin aimed at obtaining the relative density of the sample by counting the number of voxels belonging to the solid phase (the white voxels) and dividing it by the total number of voxels. Plots of the relative density in each slice as a function of the location of these slices were also created. A second plugin measured the characteristic size of each phase (rod size and pore size). This was made by applying a series of 3D erosion and dilation operation with an increasing size to the phase of interest. This measurement can be compared to sieving of a powder [START_REF] Maire | Characterization of the morphology of cellular ceramics by 3D image processing of X-ray tomography[END_REF]. The results were presented as the pores' or rods' sizes distributions in the sample.

Determination of the Young's modulus by nanoindentation

Nanoindentation was performed by an instrumented nanoindentation system (Agilent G200, Agilent Technologies, USA). A Berkovich tip indented polished rods with a speed of 10 nm.s -1 , a dwell time of 10 s, and a maximum load of 100 mN. An indentation depth of 2 µm was chosen to have an indent size sufficiently large to be representative of the entire rod's elastic properties (i.e., the ceramic phase and the potential residual micropores). A series of 10 indents was made with one indent per rod. The Young's modulus was calculated from the force/displacement curve using the Oliver-Pharr method [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF]. A value of 28 ± 2 GPa was found and used in the FE model (see section 2.7). This value is in agreement with values of Young's modulus obtained by nanoindentation for porous HA and β-TCP samples with the same volume fraction of micropores considering a rule of mixtures for BCP components [START_REF] Wang | Influence of sintering temperatures on hardness and Young's modulus of tricalcium phosphate bioceramic by nanoindentation technique[END_REF][START_REF] He | Mechanical behavior of porous hydroxyapatite[END_REF].

FE model

The image of the initial sample previously fabricated (see section 2.4) was used as an input to build the FE model. The rods (i.e., the ceramic and the artificial defects) were meshed with 1124000 quadratic tetrahedra having an average edge size of 69 µm using Avizo 8 software [START_REF] Caty | Modelling the properties of closed-cell cellular materials from tomography images using finite shell elements[END_REF]. To differentiate the mesh elements corresponding to the ceramic phase from those corresponding to the artificial defects, a special Java program was developed to create an Abaqus input file where the elements were classified according to their corresponding phase. This program associated the grey levels of the initial image (i.e., 128 for the artificial defects and 255 for the ceramic phase) to each element of the mesh. Two element sets were considered in the simulation: a first one for the elements corresponding to the ceramic and a second one for the elements corresponding to the artificial defects. The macropores between rods were not meshed. A relatively coarse mesh was chosen to ensure a reasonable calculation time, with an average size of 69 µm, larger than the mean defect size comprised between 35 and 84 µm after sintering. An elastic behavior was assumed for the two phases, as this calculation aimed at evaluating the local stresses in the sample. Calculations including a failure criterion were beyond the scope of this paper. The average Young's modulus of the rods was determined by nanoindentation, as previously explained. For the elements corresponding to artificial defects, an empirical Young's modulus of 1 MPa was considered.

The mesh was coarse because the elements' average size of 69 µm was higher than the defects' size. Thus the elements associated to artificial defects correspond to zones containing mainly artificial defects but also some ceramic phase. A Poisson's ratio of 0.28 was considered for the two set of elements [START_REF] Grenoble | Elastic properties of hard tissues and apatites[END_REF]. Different negative displacements equivalent to the strains reached during the different loadings in the ex situ test were applied to the nodes of one side of the sample whereas the nodes on the other side were displacement constrained. A reference calculation where all the elements were assigned the behavior of the ceramic phase (i.e., as if the rods did not contain any artificial defect) was also performed. Comparison between the results given by the two calculations could bring information about the influence of the artificial defects on the macroscopic behavior. Fig. 3 gives the relative density distribution measured with the high-resolution images along the sections perpendicular to the three main directions of the sample shown in Fig. 2. Due to the spatial resolution of the technique, micropores resulting from incomplete sintering are not detected and this measure assumes fully dense rods. The relative density calculated from ImageJ is 0.60 ± 0.01, meaning a volume fraction of macropores in the whole sample of 0.40 ± 0.01. The evolution of relative densities in the cross sections in directions 1 and 2 between a highest value of 0.85 and a lowest of 0.30 is explained by the alternating orthogonal rods. The cross sections with the maximum relative density along direction 1 correspond to sections crossing the rods orientated in direction 2 and the less dense ones to sections crossing the sample between these rods (Fig. 3). Along directions 1 and 2, a higher density is calculated in the first and last cross sections corresponding to the sides of the samples. The higher density is due to the deposition done at the end of the rods. To go from one rod to another, the nozzle does not stop the deposition as shown in Fig. 2b in the plane perpendicular to the deposition. Along direction 3, the variations of the relative density are less logically important. Fig. 4a shows the rod size distribution measured by ImageJ with the 3D granulometry plugin. The rod diameter lies between 119 and 161 µm with a peak for a diameter of 147 µm and a smaller peak for 126 µm. The pore size distribution measured by 3D image analysis (Fig. 4b) is characterized by a main peak between 84 and 133 µm with a maximum frequency at 120 µm. This peak is related to the macroporosity between the rods. The pore size distribution also

Results

Architecture of the sample

Characterization of the microstructure in the solid rods

Local tomography (Fig. 5a) gives details on the microstructure of the rods: intentionally added artificial defects (some are highlighted in a red frame) and imperfect boundaries (indicated by an arrow) between two orthogonal layers. As seen in Fig. 5b, the small artificial defects are randomly located in the rods. Due to this random location, they are connected or not to the macropores located between the rods. Some of them have a non-spherical shape (Fig. 5a): these pores with a more elongated shape are probably formed by the burning out of two adjacent PMMA beads.

5: (a) Slice of a high-resolution image (voxel size of 3 µm) from local tomography showing artificial defects and incomplete bonding between the different layers (arrows), (b) 3D visualization of the artificial defects from local tomography images

To investigate in more detail the microstructure of the rods, a polished cross section was prepared and observed with SEM. Residual micropores due to incomplete sintering are revealed in the rod (Fig. 6a). This residual porosity is also noticed with specific very highresolution tomography by a difference of grey level values in the rods (Fig. 6b). This tomographic image has a voxel size of 0.7 µm and was acquired using a small piece of sample (with dimensions of approximately 1 x 1 x 0.7 mm 3 ). Moreover, the incomplete sintering at the junction between the two orthogonal rods (already seen in local tomography images in Fig. 5a) is visible in this image.

6: (a) SEM image of a cross section of a rod (magnification x 1000) and (b) an X-ray tomography high resolution image (voxel size of 0.7 µm) of several rods

The XRD pattern of the scaffold is shown in Fig. 7. Only HA and -TCP phases are detected and no trace of -TCP is found. HA/β-TCP weight ratio of 19/81 is obtained. This value is close to the targeted one (20/80). The size of the residual pores left by incomplete sintering measured through SEM images is 2.6 ± 1.1 µm. In addition, the interconnection pore size distribution of the robocast HA/β-TCP sample was measured by mercury intrusion porosimetry (Fig. 8). The peak between 47 and 114 µm with a maximum value of 91 µm corresponds to the interconnections through the macropores and the open artificial defects. The second peak between 0.3 and 2.4 µm with a maximum value of 1.3 µm is linked to the intrusion of mercury into the residual micropores left by incomplete sintering.

Fig. 8: Pore size distribution measured by mercury porosimetry

Mercury porosimetry reveals a total volume fraction in the sample of micropores with a size lower than 15 µm of 0.26 ± 0.02. This value is coherent with the difference between the volume fraction of pores obtained by geometric measurement (total porous volume, 0.64 ± 0.02) and by image analysis of the tomographic volume (macroporous volume, 0.4 ± 0.01). Consequently, the respective volume fractions of each type of pores in the whole sample are 0.4 for macropores, 0.26 for residual micropores and 0.012 for artificial defects. These values give a total volume fraction of pores in the whole sample of 0.67 which is coherent with the geometric value of 0.64. The difference can be explained by some errors in the measurement (resolution of the tomographic volume used, presence of some micropores which were not impregnated by mercury during porosimetry measurement). In the rods, the volume fraction of pores is 0.02 for the artificial defects and 0.43 for the micropores. on the beginning of the stress/strain curve is also shown in Fig. 9b. The shape of the curve is discussed in the next section. The mechanical behavior is characterized by a linear elastic part followed by a plateau before brittle failure. The compression strength (i.e., the maximum stress reached during the test) is 14.2 MPa. Table 1 presents the values of compression strength for the sample and values from literature which will be commented in the Discussion section. Average on the four unloading curves 2.31 ± 0.33

Ex situ compression tests

First unloading curve (Fig. 8) 1.84

Fourth unloading curve (Fig. 8)

2.65

Average on the four loading curves 1.47 ± 0.19 Fig. 9 presents a loading and two unloading curves for the sample at low strains used to calculate the Young's modulus. The stiffness is higher during the unloading cycle compared with the loading cycle (Fig. 9), and reaches the highest value at the onset of unloading. The stiffness during unloading increases with the cycle number. Loading/unloading hysteresis is also visible in the curve as well as an irreversible residual strain at the end of the unloading steps. The effective Young's modulus was calculated from the loading and unloading parts of the loading/unloading cycles (Fig. 10). The different values obtained are summarized in Table 1. A higher value of Young's modulus (2.3 GPa) is obtained compared with the value obtained with the loading curve (1.5 GPa). The Young's modulus increases as the loading/unloading cycle number increases. The reason for these differences is discussed in the next section. In the case of the Young's modulus, a comparison with the literature is very difficult because this parameter is highly dependent on the relative density, the fabrication process, and the volume fraction of the two phases (HA and β-TCP). In a previous study, we obtained a Young's modulus of 1.38 GPa (calculated with unloading curves) for a randomly porous β-TCP with a relative density of 0.29 [START_REF] Petit | Mechanical behavior of a β-TCP ceramic with a random porosity: study of the fracture path with X-ray tomography[END_REF]. The second slice (Fig. 12) located inside the sample shows the crack propagation in direction 2. The crack reaches this section of the sample after the third loading (see the red arrows in Fig. 12d). At the end of the test, the crack separates the sample into two parts. In addition, the insight in Fig. 11c shows the crack initiates in a rod where an artificial defect is present. A second defect is also present close to the first one in the left. These two close defects can act as stress concentrators and initiate the crack. However, this image does not allow for concluding exactly about the precise location of the crack initiation because the crack could also be initiated from the macropores located below. Section 3.4 details the results of the FE calculations. on the crack initiation are presented in Fig. 13c (tomographic slice), 13d (reference calculation without the artificial defects) and 12e (calculation with the artificial defects). Without the artificial defects (Fig. 13d), high stresses are predicted in the horizontal rods and especially near the joints with the orthogonal rods. This is consistent with the experimental observations, which show the propagation of the crack in the horizontal rods. The stress distribution in the section is logical because it is controlled by the periodic structure. When including the artificial defects, the stress distribution appears to be different (Fig. 13e). Different stress concentrations are observed in the vicinity of the artificial defects. In the zone where the crack initiates, the calculation without the artificial defect does not predict higher stress between the two macropores (see the yellow frame in Fig. 13d). When the artificial defects are taken into account, stress concentrations appear around the defect (see the yellow frame in Fig. 13e). A comparison between Fig. 13d and 13e shows that the stress concentrations are present along the crack path. The second artificial defect present in the left also creates stress concentration in its vicinity (see the purple frames in Fig. 13d, and13e).

Finally, with the artificial defects in the model, an increase of stress is noted in a horizontal rod where no artificial defect is present (compare the bottom part of the images in Fig. 13d and13e).

Discussion

Characterization of the robocast samples

The cellular ceramic sample considered here shows a square prism periodic structure. It has a high volume fraction of pores (0.64). However, the macropores, which characterize the macroscopic structure, only represent 0.4 of the whole sample. The difference between these two values represents the fraction of micropores and artificial defects in the solid rods of the sample. The sintered rods are not fully densified, as testified by an internal volume fraction of micropores of 0.43. This can be related to the relatively low sintering temperature of 1100°C leading to a limited densification of the ceramic, as already noted by several authors [START_REF] Chu | Research of phase transformation induced biodegradable properties on hydroxyapatite and tricalcium phosphate based bioceramic[END_REF][START_REF] Prokopiev | Dependence of the mechanical properties of sintered hydroxyapatite on the sintering temperature[END_REF]. However, this temperature was chosen to avoid the phase transformation from β to α-TCP which generally occurs around 1150°C [START_REF] Ryu | An improvement in sintering property of β-tricalcium phosphate by addition of calcium pyrophosphate[END_REF][START_REF] Kivrak | Synthesis of calcium hydroxyapatite-tricalcium phosphate (HA-TCP) composite bioceramic powders and their sintering behavior[END_REF]. A small amount of α-TCP phase has already been detected in similar HA/β-TCP scaffolds sintered at 1140°C [START_REF] Houmard | On structural, mechanical and biodegradation properties of HA-TCP robocast scaffolds[END_REF]. It is known that α-TCP can have a detrimental effect on the mechanical properties [START_REF] Chu | Research of phase transformation induced biodegradable properties on hydroxyapatite and tricalcium phosphate based bioceramic[END_REF][START_REF] Ryu | An improvement in sintering property of β-tricalcium phosphate by addition of calcium pyrophosphate[END_REF]. HA/β-TCP weight ratio of 84/16 was found in the sintered sample whereas a value of 80/20 was targeted. These close values show that the robocasting and sintering did not affect the sample's chemistry. During the robocasting process, the syringe diameter was 150 μm. Considering the solid loading of the ceramic ink of 30 vol% and the relative density of the sintered rods of 57%, the volume shrinkage is 46% and the linear shrinkage is 13%. Thus, the expected rod diameter after shrinkage is 130 μm. This is consistent with the measured rod diameter distribution between 119 and 161 µm. The highest values probably correspond to rod overlapping between two adjacent rods. It is important to note that we neglect a possible elastic movement of the ink when it is deposited. A rod spacing (from center to center of the rods) of 300 µm was chosen for the ink's deposition. Considering the nozzle diameter of 150 µm, it means that during the process, the space between two rods was 150 µm. Since the rods shrink in every dimension during sintering, the macropore size is also expected to decrease in the same proportion: thus the expected macropore size is also 130 µm (for 150 µm during deposition). The experimental macropore size is between 84 and 133 µm. This may be explained by an elastic "displacement" of the ink when it gets out of the nozzle and by the rod collapse during drying. This phenomenon is shown in Fig. 6b, with an interpenetration of the rods between adjacent layers and of non-perfect bonding between rods. The use of local tomography allows the observation of these types of defects due to fabrication processes (Fig. 5a and6b) even if the overall structure seems to be regular (Fig. 2a). Local tomography also shows the artificial defects having a size between 35 and 84 µm. After shrinkage, the expected artificial defects' size range is 39-87 µm due to the PMMA beads' size range (45-100 µm). The robocast samples are often observed by microscopy, which does not provide a complete characterization [START_REF] Dellinger | Robotic deposition of model hydroxyapatite scaffolds with multiple architectures and multi-scale porosity for bone tissue engineering[END_REF][START_REF] Franco | Direct write assembly of calcium phosphate scaffolds using a water-based hydrogel[END_REF]. This study shows the importance of a 3D characterization of the scaffolds, especially with local tomography.

Mechanical behavior in compression

The stress/strain curve of the sample in uniaxial compression is characterized by a first elastic domain followed by a stress plateau. Non-linearity is observed in the elastic part of the loading curves with an increase of the axial stiffness (Fig. 9b) as the applied stress increases. This has already been observed during compression tests of porous ceramics [START_REF] Bruno | Porous microcracked ceramics under compression: Micromechanical model of non-linear behavior[END_REF] or rocks [START_REF] David | Sliding crack model for nonlinearity and hysteresis in the uniaxial stress-strain curve of rock[END_REF]. The increase of stiffness with increasing load has been attributed to a progressive closure of microcracks favorably oriented relative to the loading direction. High-resolution images show the presence of such defects, such as imperfect bonding between the layers (Fig. 5a and6b). The non-linearity can also be explained by the removal of the samples after each loading cycle for scanning. Crack re-opening possibly occurs during sample handling between scans.

Additionally, the two faces in contact with the compression platens were embedded with paraffin wax to avoid local brittle fracture near the faces. Therefore, the increase of stiffness can also be due to a first compression of the wax before the "real" compression of the specimen. Then, at the onset of unloading, the stress/strain curve is linear: no mechanism associated with non-linearity is active. The axial stiffness slowly decreases during unloading, possibly due to the opening of the closed cracks and/or sliding of the closed cracks' faces [START_REF] Bruno | Porous microcracked ceramics under compression: Micromechanical model of non-linear behavior[END_REF][START_REF] David | Sliding crack model for nonlinearity and hysteresis in the uniaxial stress-strain curve of rock[END_REF]. This axial stiffening of the unloading curve can explain the higher value obtained for the effective Young's modulus calculated during unloading as compared to loading (Table 1). Finally, the increase of the stiffness at the onset of the unloading curves with the loading/unloading cycle number (Fig. 10) can be linked to an increased crack closure at higher applied strains. At high applied strains, the elastic domain is followed by a stress plateau: it usually corresponds to the damageable behavior commonly observed for highly porous ceramics with a stochastic microstructure. In this latter case, it is characterized by the propagation of multiple small cracks in the solid phase of stochastic porous ceramics [START_REF] Pecqueux | Influence of microporosity and macroporosity on the mechanical properties of biphasic calcium phosphate bioceramics: modelling and experiment[END_REF]. In this study, although the robocast sample has a total volume fraction of pores of 0.64, a behavior characteristic of dense brittle ceramic is observed with the propagation of a main crack parallel to the loading direction. The rupture is not perfectly brittle because the ex situ compression test shows that the crack propagates step by step along the sample. This specific behavior has already been observed by Miranda et al. [START_REF] Miranda | Fracture modes under uniaxial compression in hydroxyapatite scaffolds fabricated by robocasting[END_REF] with robocast HA scaffolds. Because of the high volume fraction of pores, we could have expected a damageable behavior with the propagation of multiple small cracks. It seems that the observed brittle-like behavior is rather explained by the volume fraction of macropores (0.4) significantly lower than the total porous volume. A transition from a brittle-like to a cellular-like behavior with the increase of volume fraction of pores has been observed by Meille et al. [START_REF] Meille | Mechanical properties of porous ceramics in compression: on the transition between elastic, brittle, and cellular behavior[END_REF] for a macroporous alumina with a stochastic microstructure. In this latter study, the transition between the two behaviors occurred at a volume fraction of pores of around 0.5. It seems therefore than the brittle or damageable character in uniaxial compression is mainly controlled by the macroporous volume fraction. It confirms observations made in our previous study about a macroporous β-TCP sample [START_REF] Petit | Mechanical behavior of a β-TCP ceramic with a random porosity: study of the fracture path with X-ray tomography[END_REF]. This sample having a volume fraction of macropores of 0.5 exhibited a damageable behavior. Specific use of two-scale X-ray tomography in the study of cellular ceramics Local tomography allows to have a new insight on the rods' structure and the sample's mechanical behavior. They allowed observing the artificial defects in the rods which have an influence on the crack initiation. In addition to the artificial defects, other defects (e.g., imperfect bonding between layers) are detected. The presence of defects in structures made by AM and their influence on the mechanical behavior have already been noticed in metallic samples [START_REF] Petit | CoCrMo cellular structures made by Electron Beam Melting studied by local tomography and finite element modelling[END_REF]. This means that crack initiation and propagation could be delayed if the defects due to the fabrication process could be removed. Therefore, further work is needed to improve the robocasting process and obtain defect-free specimens. The high-resolution images were also processed to be used as an input for FE modeling. In this paper, an original two-scale method was developed to take into account the details of the microstructure revealed by local tomography in a model of the material at a lower resolution. This avoided the use of very large images and long calculation times. Our approach allowed having a lightweight mesh of the macrostructure including information about a finer scale. The calculated local stress distribution is modified by the presence of the artificial defects (see Fig. 13d ande). A high stress concentration is present around the two main defects where the crack initiates during the experimental test (Fig. 13e). However, another stress concentration is visible below a macropore where no artificial defect is present (see the bottom part in Fig. 13e) and no crack initiation occurs at this location. This stress concentration may be related to the imperfect bonding between adjacent layers which can be distinguished in the tomographic slice (see the bottom part of Fig. 13c). This means that the link between the stress distributions predicted by FE models and the crack path is not obvious. To perform calculations within a reasonable time, we meshed the 3D images with tetrahedra with an average edge size of 69 µm. This resulted in a relatively coarse mesh with some approximations in the geometrical description of the sample, especially around the artificial defects. A refinement of the mesh or the use of higher resolution images could bring more accurate results. The simple elastic simulation presented here could also be completed by calculations with a more realistic behavior (for instance, damageable behavior) for the ceramic material.

Conclusions

This paper focuses on the analysis of the fracture behavior in compression of a robocast biphasic calcium phosphate scaffold. This sample consist of different layers of parallel rods, the rods of one layer being perpendicular to the rods of the adjacent layer. They contain three types of pores: a macroporosity of hundreds of μm between the layers defining the architecture of the specimens, artificial defects with a size 30-80 μm intentionally added in the rods, and residual pores due to incomplete sintering with a size of a few microns. The ex situ compression test reveals a brittle-like behavior. A crack propagates along the loading direction (i.e., perpendicular to the tensile stress created by Poisson's effect inside the sample) and cuts the sample through its height. Over the same time, it propagates along the specimen width and finally opens the sample in two parts. The spherical artificial defects intentionally added are shown to have an influence on the crack initiation, despite the presence of stress concentration zone at the interconnection between rods. The residual microporosity in the rods could not be observed with local tomography. Therefore, their role in the mechanical behavior is not known. X-ray tomography was used to characterize the sample at two scales. The macrostructure was imaged with a low-resolution volume whereas local tomography provided high-resolution images of the rods and especially the artificial defects. These latter images were processed to build a FE modeling having a moderate number of elements and including the presence of the artificial defects. They were included as a "second phase" different from the ceramic with a specific behavior. The maximal principal stresses obtained by this FE modeling confirm the influence of the artificial defects with an associated important stress concentration. However, other calculations implementing a more complex material behavior or with a finer mesh size would be necessary to understand better the respective role of the microstructure and the macrostructure on the scaffold's properties.

Fig. 1 :

 1 (a) Local tomography image slice of part of the sample (voxel size: 3 µm), (b) the same slice after the image processing steps separating the three phases (voxel size: 3 µm), (c) the slice after scaling (voxel size: 15 µm), and (d) the scaled slice after the second thresholding (voxel size: 15 µm)

Fig. 2

 2 Fig.2shows a 3D reconstruction and a SEM image of a part of the sample showing its architecture. The alternation of orthogonal rods from one layer to an adjacent one is clearly visible (see Fig.2a).

  35 and 84 µm, corresponding to the artificial defects added by PMMA beads in the rods.

Fig. 3 :

 3 Fig. 3: Relative density distributions along the three main directions calculated by 3D image analysis

Fig. 7 :

 7 Fig. 7: XRD pattern of the sintered scaffold. HA (*) and β-TCP (○) are identified with comparison with JCPDS cards 09-0432 and 09-0169 respectively.

Fig. 9a presents

  Fig.9apresents the stress/strain curves obtained for the sample during the ex situ test. The black curve represents the entire stress/strain curve, until the rupture of the sample. The other curves correspond to four loadings interrupted to scan the sample at different strains. An inset

Fig. 9 :

 9 Fig. 9: (a) Compression stress/strain curves corresponding to the successive loadings and to the test until rupture (black curve) and (b) a magnified view of the final stress/strain curve's beginning showing the increase of stiffness in the linear part

Fig. 10 :

 10 Fig. 10: Loading and unloading cycles undergone by the sample. The loading curve corresponds to the first loading curve of the Fig. 8a.

Fig. 11 :Fig. 12 :

 1112 Fig. 11: (a) 3D reconstruction of the sample showing the location of a vertical slice (red frame) which is presented after at different steps of the ex situ compression test: (b) at the initial state, (c) after the second loading (with a magnified view on the crack initiation), (d) after the third loading, and (e) after the fourth loading (red arrows show the crack initiation and propagation)

Fig. 13 :

 13 Contours of the calculated maximal principal stress for the calculation (a) without and (b) with the artificial defects, (c) magnified view of the slice from tomographic image on the crack initiation, and (d) and (e) magnified views of the contours where the crack initiates. In (d) and (e), a yellow arrow points out the stress field around the artificial defect crossed bythe crack, and a yellow frame shows the stress field around the second artificial defect. The yellow arrow and frame are reported in (c) to illustrate the exact location of the defects.

Table 1 :

 1

Values of effective Young's modulus calculated from the loading/unloading curves

Calculation

Effective Young's modulus (GPa)

Acknowledgements

Research reported in this publication was supported by the National Institute of Dental and Craniofacial Research of the National Institutes of Health under award number 1R01DE015633 and by the French Ministry of Education and Research (grant for C. Petit PhD, 2012-2015).


A logic identifying isomorphic propositions

Alejandro Díaz-Caro, Gilles Dowek

To cite this version:

Alejandro Díaz-Caro, Gilles Dowek. A logic identifying isomorphic propositions. 2018. hal-01109104v5

The paper has been withdrawn since we have found an error in the proof of strong normalisation and a counter-example for such a property: Let

y : τ z : τ ⇒ τ ⇒ τ δ = λx τ .π τ ⇒τ ((zy) × x)x : τ ⇒ τ Ω = δ(δy)
Then Ω → π τ ⇒τ ((zy) × (δy))(δy)

π τ ⇒τ ((z × δ)y)(δy)

π

  τ ⇒τ ⇒τ (z × δ)y(δy) π τ ⇒τ ⇒τ (z × δ)(δy)y π τ ⇒τ ((z × δ)(δy))y π τ ⇒τ ((z(δy) × Ω)y 1