Xiaoyu Xu

Zhuoxiang Ren

Zhi Qin

Hakeim Talleb

Shuai Yan

Application of PGD on parametric modeling of a piezoelectric energy harvester

Keywords: Proper generalized decomposition, Model reduction, Finite element method, Piezoelectricity, Energy harvesting

come

I. INTRODUCTION

HEN designing an energy harvester, it is crucial to investigate its performance under various conditions. Take the cantilevered piezoelectric EH for example, one usually seeks to calculate its output power when the excitation frequency and electrical loads connected between its electrodes vary. The conventional approach is modeling the device with numerical methods such as the FE method while taking the varying variables as parameters. Such parametric modeling leads to high dimensional problems that require a vast amount of computational resource to solve. In practical applications, appropriate methods must be introduced to alleviate the computational complexity.

Recently, researchers have developed various model reduction techniques, in either a posteriori [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF] or a priori [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] fashion. The former consists in extracting from available solutions a reduced basis which alleviates the burden in subsequent calculations. Algorithms like the singular value decomposition (SVD), principle component analysis (PCA) and proper orthogonal decomposition (POD) are popular for the purpose. The second category, in contrast, is normally more challenging as it searches a reduced basis without knowing the solution a priori. The proper generalized decomposition is one of the dominant methods under this framework [START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations: a primer[END_REF]. This work introduces applying the PGD method on the parametric model of a piezoelectric EH.

The essential assumption of the PGD method is that the numerical solution can be expressed as low rank tensors, or a sum of rank-one tensors [START_REF] Nouy | A low-rank methods for high-dimensional approximation and model order reduction[END_REF]. Each term in the sum is called a mode while the number of terms is referred to as the tensor rank [START_REF] Grasedyck | A literature survey of lowrank tensor approximation techniques[END_REF]. From a numerical perspective, each mode is a multi-dimensional array in which every dimension is related to a physical quantity e.g. spatial dimensions and frequency. In practice, modes are calculated successively, which is termed as the enrichment procedure. During each enrichment, simultaneously determining all vectors in the array leads to a nonlinear problem that can be effectively computed using the alternating direction method (ADM) which is sometimes named the fixed point iteration procedure [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF]. Advantage of the PGD approach is twofold. Firstly, instead of solving the high-dimensional problem through looping over each sample of the parameter (the conventional approach), it solves the problem with all samples of parameters at the same time, thus avoiding a large number of loops and reducing the CPU time. Secondly, the number of unknowns to be determined is generally dramatically reduced since PGD solutions actually form reduced basis of the conventional solutions. Indeed, for a vast range of problems, only moderate number of enrichments and few iterations within each enrichment are required. Naturally, less solution data needs to be stored which is crucial in real-time simulations [START_REF] Nirromandi | Real-time simulation of biological soft tissues : a PGD approach[END_REF]. The PGD method also relies on the separability of the problem. Separability of the right hand side (r.h.s.) is usually feasible because excitations and boundary conditions based on which the r.h.s. is formed are normally known. For parameter dependent cases, a posteriori techniques e.g. SVD can be used to obtain a separated representation. On the contrary, separating the coefficient matrix in a general case is still an open problem, albeit some literatures are now available. For instance, in [START_REF] Aguado | Advanced strategies for the separated formulation of problems in the proper generalized decomposition framework[END_REF] operator separating in both linear and nonlinear cases are discussed. In [START_REF] Modesto | Proper generalized decomposition for parameterized Helmholtz problems in heterogeneous and unbounded domains: Application to harbor agitation[END_REF], the same issue in high-dimensional cases are investigated. After separating, PGD formulations can be constructed using the Galerkin orthogonality criteria or residual minimization, depending on whether the operator is self-adjoint [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. Convergence accelerating methods and optimality of the solution (i.e. attaining the same accuracy with as fewer modes as possible) are also discussed in the same reference. A recently developed approach named ideal minimal residual method is reported in [START_REF] Billaud-Friess | A tensor approximation method based on ideal minimal residual formulations for the solution of high-dimensional problems[END_REF] for improving optimality in highdimensional cases. It is based on measuring the residual and the target tensors with appropriately chosen norms. As theoretical framework of the PGD becoming more and more mature, its applications for specific problems are also increasingly visible, such as uncertainty quantification in [START_REF] Chevreuil | A model order reduction based on proper generalized decomposition for the propagation of uncertainty in structural dynamics[END_REF], mechanics in [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and non-incremental methods of calculation[END_REF], magneto quasi-static problems in [START_REF] Henneron | Proper generalized decomposition method applied to solve 3D magnetoquasi-static field problems coupling with external electric circuits[END_REF] and nonlinear magnetic dynamic problems in [START_REF] Qin | A proper generalized decomposition based solver for nonlinear magneto-thermal problems[END_REF].

W

Although literatures introducing PGD for mechanics or electromagnetics are now available, those involving both mechanics and electromagnetics remain to be exploited. Such problems are particularly important in modeling energy harvesting devices in which energy conversion between mechanical and electromagnetic forms is considered, along with various working conditions and energy storing issues [START_REF] Talleb | Finite element modeling of a magnetoelectric laminate composites in considering nonlinear and load effects for energy harvesting[END_REF][START_REF] Talleb | Finite element modeling of a magnetoelectric energy transducer including the load effect[END_REF]. It is generally much more computationally expensive to solve these multiphysics problems because the strong coupling effect between physics can significantly increase the problem size, especially in large scale parametric FE models. Here, we propose to apply the PGD method on the FE model of a piezoelectric EH. Although based on appropriate simplifications and developments, analytical models (see e.g. [START_REF] Zhao | Elecroelastic modeling and experimental validation of piezoelectric energy harvesting from broadband random vibrations of cantilevered bimorphs[END_REF]) are less computationally intensive than FE models when dealing with devices of simple geometries and boundary conditions, they become less applicable when the model becomes more sophisticated. FE models, on the other hand, are advantageous due to their versatility in modeling complicated geometry and its capability in obtaining full field solutions under any boundary conditions. In [START_REF] Elvin | A coupled finite element -circuit simulation model for analyzing piezoelectric energy generators[END_REF], FE methods are utilized for analyzing piezoelectric energy generators using different commercial software modeling multiphysics fields and the electric circuits, individually. Dedicate loops are designed in order to pass results between software. In this paper, PGD is applied so that variables are decomposed into space, frequency and electrical load related components, and the original large-size system is solved through solving a sequence of smaller-size subsystems, thus reducing complexity of the problem. We also introduce a method to incorporate electrodes and electrical loads in the piezoelectric model so that separability of the problem is straightforward which enables getting around separation of the coefficient matrix. In this manner, advantages of the FE model over its analytical counterpart are preserved while the CPU time is significantly decreased. The paper is divided as follows. In section II, formulations of the multiphysics model, with excitation frequency and electrical load as parameters, are derived. In section III, how solution of the model can be approximated with low rank tensors using PGD is introduced. Implementing details are then presented. In Section IV, properties of the proposed method are revealed through a set of numerical illustrations. It finishes with general summaries and conclusions in section V.

II. FORMULATIONS OF THE MULTIPHYSICS MODEL

In what follows, quantities of the same uses are expressed in uniform representations. Specifically, lower-case letters are used for scalars; plain capitals for vectors; bold capitals for matrices; and calligraphic letters for tensors of order higher than two.

A. The model problem

As in most piezoelectric FE models, only linear piezoelectricity [START_REF]IEEE Standard on Piezoelectricity[END_REF] is considered here. Assume that we have a threedimensional domain Ω φ ∈ ℝ 3 A common practice of modeling piezoelectric materials is to choose Ω φ coinciding with Ω u which, however, does not appropriately account for the fringing effect, from an electrical modeling point of view. Here, domains are chosen so that Ω φ ∕ Ω u ≠ ∅, as shown in Fig. 1. Indeed, Ω φ ∕ Ω u represents an air domain surrounding the piezoelectric part. The strong form of the coupled problem (for detailed derivation refer to [START_REF] Qin | Finite element modeling of thermos-elastic attenuation in piezoelectric surface acoustic wave devices[END_REF] and references therein) consists of finding fields of mechanical displacements and electrical potentials subjecting to certain boundary conditions and excitations, as shown in the coupled mechanical equilibrium equation

c ijkl E u k,li + e kij φ ,ki + f j = ρ v üj in Ω u
(1.a)

u i = u ̅ on Γ u D (1.b) n i σ ij = t̅ on Γ u N (1.c)
and the coupled electrostatic equation Electrodes are deposited on the piezoelectric part and impose equipotential electrical conditions. Influence of the electrodes are mechanically negligible since they are very thin. Therefore, electrodes can be considered as surfaces and involves merely imposing the equipotential condition. Electrodes and electrical loads are introduced after the fields are approximated using FEs.

e kij u i,jk -ε ij S φ ,ij = 0 in Ω u (2.a) -ε ij S φ ,ij = 0 in Ω φ ∕ Ω u (2.b) φ = φ ̅ on Γ φ D (2.

B. Finite element approximation

To approximate the coupled problem using FEs, consider the following weak form. Find u i ∈ 𝒰 and φ ∈ 𝒱 such that

B 𝑢𝑢 (u′ j , u j ) Ω u + B 𝑢𝜑 (u′ j , φ) Ω u = ∀u ′ j ∈ 𝒰 ̂ (3)
-(u ′ j , f j ) Ω u + (u ′ j , t̅ ) Γ u N B 𝜑𝜑 (φ′ , φ) Ω φ + B 𝜑𝑢 (φ′ , u i ) Ω u = 0 ∀φ ′ ∈ 𝒱 ̂ (4)
where the integrals are defined as

B 𝑢𝑢 (u′ j , u j ) Ω u : = ∫ u′ j (c ijkl E u k,li -ρ v üj) Ω u dΩ (5.a) B 𝑢𝜑 (u′ j , φ) Ω u : = ∫ u′ j e kij φ ,ki Ω u dΩ (5.b) (u′ j , f j ) Ω u ∶= ∫ u ′ j f j Ω u dΩ (5.c) (u′ j , t̅ ) Γ u N ≔ ∫ u ′ j t̅ Γ u N dΓ (5.d) B 𝜑𝜑 (φ′ , φ) Ω φ ≔ ∫ φ ′ ε ij S φ ,ij Ω φ dΩ (5.e) B 𝜑𝑢 (φ′ , u i ) Ω u ≔ ∫ φ ′ e kij u i,jk Ω u dΩ (5.f)
The admissible spaces are

𝒰 ∶= {u j ∈ H 1 (Ω u ) | u j = u ̅ on Γ u D , j = 1,2,3} (6.a)
𝒰 ̂∶= {u′ j ∈ H 1 (Ω u ) | u′ j = 0 on Γ u D , j = 1,2,3} (6.b
)

𝒱 ∶= {φ ∈ H 1 (Ω φ ) | φ = φ ̅ on Γ φ D } (6.c) 𝒱 ̂∶= {φ′ ∈ H 1 (Ω φ ) | φ′ = 0 on Γ φ D } (6.d)
where H 1 is the first order Sobolev space.

The discrete counterparts of Eq.(3) and Eq.( 4) read: find u jh ∈ 𝒰 h and φ h ∈ 𝒱 h such that

B 𝑢𝑢 (u jh ′ , u jh ) Ω u + B 𝑢𝜑 (u jh ′ , φ h ) Ω u ∀u jh ′ ∈ 𝒰 ̂h (7) = -(u jh ′ , f j ) Ω u + (u jh ′ , t̅ ) Γ u N B 𝜑𝜑 (φ h ′ , φ h ) Ω φ + B 𝜑𝑢 (φ h ′ , u ih ′) Ω u = 0 ∀φ h ′ ∈ 𝒱 ̂h (8)
where 𝒰 h ⊂ 𝒰 , 𝒰 ̂h ⊂ 𝒰 ̂, 𝒱 h ⊂ 𝒱 , and 𝒱 ̂h ⊂ 𝒱 ̂ are finite dimensional spaces, with h characterizing the triangulation resolution. Within a given finite element 𝑒 of domain Ω u e (resp. Ω φ e ), the discretized mechanical (resp. electrical) test and trial functions can be constructed using the Galerkin method, as in Eq.( 9) .

u jh ′ = ∑ N m (u) n u m=1 u′ j m e (9)
u jh = ∑ N m (u) n u m=1 u j m e
and Eq.( 10)

φ h ′ = ∑ N m (φ) n φ m=1 φ′ m e ( 10
)
φ h = ∑ N m (φ) n φ m=1 φ m e
where n u and n φ are total number of nodes in the mechanical and electrical domain respectively. N m (u) and N m (φ) are shape functions and the same shape function N m (u) is used for all components. u′ j m e , u j m e , φ′ m e and φ m e are nodal values.

In the implementation, it is more convenient to arrange u jh ′s and φ h ′s in vectors. For instance, u jh ′s can be written as

U h = [ u 1h u 2h u 3h] T .

Re-writing the above equations in the compact form yields

U h ′ = 𝐍 u e U ′ , U h = 𝐍 u e U (11.
𝐍 u e = [N 1 (u) 𝐈 3 N 2 (u) 𝐈 3 ⋯ N n u (u) 𝐈 3 ] (12.e) N φ e = [N 1 (φ) N 2 (φ) ⋯ N n φ (φ) ] (12.f)
and 𝐈 3 is a 3 × 3 identity matrix.

The FE system can be obtained through substituting Eq.( 11) into Eq.( 7) and Eq.( 8) and integrating. Matrices and vectors are then calculated element-wise: for the α th element

e α ∈ Ω u 𝐌 u (α) U ̈α + 𝐊 u (α) U α + 𝐊 uφ (α) Φ α = F X (α) (13.a) 𝐊 uφ (α),T U α -𝐊 φ (α) Φ α = 0 (
𝐍 u e α = [N α1 (u) 𝐈 3 N α2 (u) 𝐈 3 N α3 (u) 𝐈 3 N α4 (u) 𝐈 3 ]
(15.c)

𝐁 u e αi = [ ∂N αi (u) ∂x 1 0 0 0 ∂N αi (u) ∂x 2 0 0 0 ∂N αi (u) ∂x 3 ∂N αi (u) ∂x 2 0 ∂N αi (u) ∂x 3 ∂N αi (u) ∂x 1 ∂N αi (u) ∂x 3 0 0 ∂N αi (u) ∂x 2 ∂N αi (u) ∂x 1 ] 𝑇 (15.d)
To model the structural damping, Rayleigh damping can be added which introduces the damping matrix as a fraction of the mass and stiffness matrix using α dM and β dK :

𝐂 u (α) = α dM 𝐌 u (α) + β dK 𝐊 u (α) (16)
The parameters α dM and β dK can be obtained from Eq.( 17) where ξ denotes the damping factor at resonant frequency ω 0 . Damping factors can be measured at the first two resonant frequencies with which Rayleigh damping are calculated and used for the whole frequency range, approximately [START_REF] Zhao | Elecroelastic modeling and experimental validation of piezoelectric energy harvesting from broadband random vibrations of cantilevered bimorphs[END_REF].

ξ = 1 2 ( α dM ω 0 + β dK ω 0 ) (17)
Adding Eq.( 16) into Eq.(13.a) leads to

𝐌 u (α) U ̈α e + 𝐂 u (α) U ̇α e + 𝐊 u (α) U α e + 𝐊 uφ (α) Φ α e = F X (α) (18)
Numerically, damping also prevents singularity when Eq.( 18) is solved near frequencies of resonance in the frequency domain. Before assembling elementary matrices into a global system that applies to the whole domain, modeling of the electrodes is addressed. . These dofs can be assigned with number , n u ̂ being the number of nodes on Γ u D . In assembling, positions of each elementary matrix and vector in the global system are determined by the mapping results. The global system can be written as

C. Modeling the electrodes and electrical loads

𝐌 u U ̈+ 𝐂 u U ̇+ 𝐊 u U + 𝐊 uφ Φ = F X (19.a) 𝐊 uφ T U -𝐊 φ Φ = 0 (19.b)
where

𝐌 u ∈ ℝ N u ×N u , 𝐂 u ∈ ℝ N u ×N u , 𝐊 u ∈ ℝ N u ×N u , 𝐊 uφ ∈ ℝ N u ×N φ ,𝐊 φ ∈ ℝ N φ ×N φ , and F X ∈ ℝ N u ×1 .
Eq.( 19) can be transformed from the time domain to the frequency domain. For a time-harmonic excitation at angular frequency ω, the transforming yields

(-ω 2 𝐌 u + jω𝐂 u + 𝐊 u )U + 𝐊 uφ Φ = F X (20.a) 𝐊 uφ T U -𝐊 φ Φ = 0 (20.b)
where j is the imaginary unit.

Incorporating electric circuits into Eq.( 20) can be achieved via describing each circuit element with the so called equivalent capacitance matrix [START_REF] Wang | A finite element circuit coupled simulation method for piezoelectric transducer[END_REF] and adding the resultant matrices into the dielectric matrix of Eq.(20.b). To be more precise, 𝐊 φ is replaced with 𝐊 φ + 𝐊 ̂φ in which the effective capacitance matrix 𝐊 ̂φ is initially null. For each circuit element x ̌ connected between the p th and q th electrode, define a vector V x ̌ that is of size N φ × 1 and its p th and q th component are respectively 1 and -1 while others are zero. Z ̃x ̌ equals to -Z x ̌ when x ̌ is a capacitor; while Z ̃x ̌= Z x ̌ when x ̌ is a resistor or an inductor, with Z x ̌ denoting the impedance of x ̌. The new global system can be expressed as

(-ω 2 𝐌 u + jω𝐂 u + 𝐊 u )U + 𝐊 uφ Φ = F X (22.a) 𝐊 uφ T U -(𝐊 φ + 𝐊 ̂φ)Φ = 0 (22.b)
For EH design, one of the most frequently used configuration is that composed of a piezoelectric bimorph connected with a resistor in series, as represented in Fig. 2. In this case, the electrodes 2 and 3 are grounded while the resistor is connected between electrodes 1 and 4. Consequently, the effective capacitance matrix can be simplified as -j ωr ⁄ 𝐊 ̆φ. This configuration is condidered in the following.

The above equation is multi-parametric, depending on ω and r. Usually, results of the FE model for a wide range of these parameters are of interest. For instance, one may seek the EH output for all combinations of ω in the N ω × 1 vector V ω ≔ [ω 1 ω 2 ⋯ ω N ω ] and r in the N r × 1 vector V r ≔ [ r 1 r 2 ⋯ r N r ], so as to find the optimal working condition. A conventional approach to solve the problems is to solve Eq.( 22) with two loops for ω and r, one comprised in the other. This approach leads to the model complexity increasing exponentially with the number of parameters. This approach takes a huge amount of time, especially when the total number of parameters is large.

III. PGD BASED MODEL ORDER REDUCTION

In this section, a novel approach that enables the complexity increase linearly with the number of parameters is introduced. Note that when Eq.( 22) is solved for values of ω and r, solutions can be arranged as three-way arrays 𝔘 and 𝔙, respectively for the displacements and potentials. Each dimension in the array relates to either the spatial coordinates, frequencies or resistors. The goal with the novel ap-proach is to approximate 𝔘 and 𝔙 in terms of low rank tensors in the canonical form [START_REF] Kolda | Tensor decomposition and applications[END_REF], as depicted in Fig. 3 where 𝔜 stands for 𝔘 or 𝔙. where U Xi ∈ ℂ N u ×1 , Φ Xi ∈ ℂ N φ ×1 , U ωi , Φ ωi ∈ ℂ N ω ×1 , and U ri , Φ ri ∈ ℂ N r ×1 . '⊗' denotes the Kronecker product. R is the total number of modes.

Note that decompositions in Eq.( 23) is the solution of a system obtained by expanding Eq.( 22) for all samples of the parameters. To obtain this system, Eq.( 22) needs to be expanded in tensor form. For the r.h.s., w.l.g., assume that identical mechanical excitations are applied in all cases of ω and r. Thus, r.h.s. in Eq.( 22.a) is expanded as F X ⊗ F ω ⊗ F r where F ω and F r , whose components are all equal to 1, are vectors of size N ω × 1 and N r × 1, respectively. Expanding matrices in Eq.( 22) in tensor form is straightforward as variables related to X, ω and r are individual terms in products, thanks to the method modeling electrodes and electrical loads. The resulting tensor comprises four blocks

𝔄 = [ 𝔄 11 𝔄 12 𝔄 21 𝔄 22 ] (24)
with where 𝐃 ω (resp. 𝐃 r ) equals to 𝐈 N ω (resp. 𝐈 N r ), the identity N ω × N ω (resp. N r × N r ) matrix. The others are defined as

𝔄 11 = -𝐌 u ⊗ 𝐀 ω ⊗ 𝐃 r + j𝐂 u ⊗ 𝐁 ω ⊗ 𝐃 r +𝐊 u ⊗ 𝐃 ω ⊗ 𝐃 r
𝐀 ω ≔ diag{ω 1 2 ω 2 2 ⋯ ω N ω 2 } ∈ ℝ N ω ×N ω (26.a) 𝐁 ω ≔ diag{ ω 1 ω 2 ⋯ ω N ω } ∈ ℝ N ω ×N ω (26.b) 𝐆 ω ∶= diag{1 ω 1 ⁄ 1 ω 2 ⁄ ⋯ 1 ω N ω ⁄ } ∈ ℝ N ω ×N ω (26.c) 𝐆 r ∶= diag{1 r 1 ⁄ 1 r 2 ⁄ ⋯ 1 r N r ⁄ } ∈ ℝ N r ×N r (26.d)
To keep notations neat, modes in Eq.( 23) are designated as

𝔵 i ≔ [𝔵 1 i , 𝔵 2 i ] T = [U Xi ⊗ U ωi ⊗ U ri , Φ Xi ⊗ Φ ωi ⊗ Φ ri ] T (27)
and the expanded r.h.s.s as

𝔟 ≔ [𝔟 1 , 𝔟 2 ] T = [F X ⊗ F ω ⊗ F r , 0] T (28)
Subscript 1 and 2 indicate association with mechanical and electrical variables, respectively. Components associated with X, ω and r are accessed with an additional subscript κ with κ = 1, 2 and 3, respectively. For instance, 𝔵 11 i represents U Xi , 𝔵 23 i represents Φ ri and 𝔟 12 represents F ω , and so forth.

To accesses matrices in the blocks of Eq.( 25), a superscript κ is added to 𝔄 so that 𝔄 ij κ depicts matrices in the block 𝔄 ij relating to the space, frequency and resistor when κ equals to respectively, 1, 2 and 3. The expanded system of Eq.( 22) is

𝔄 ∑ 𝔵 i R i=1 ≈ 𝔟 (29)
The PGD method suggests calculating each 𝔵 i successively.

Namely, in the m th enrichment for m = 1, ⋯ , R, 𝔵 m is calculated on the basis of sum of all its m -1 preceding terms.

Components 𝔵 λκ m in 𝔵 m are calculated individually through e.g. fixed point iterations.

B. PGD interpreted as an optimization problem

Firstly, define suitable subspaces 𝒮 λκ so that 𝔵 λκ m ∈ 𝒮 λκ .

Examples of such subspaces are 𝒮 11 ⊂ ℂ N u ×1 and 𝒮 23 ⊂ ℂ N r ×1 . As a result, 𝔵 λ m belongs 𝒮 λ ≔ 𝒮 λ1 × 𝒮 λ2 × 𝒮 λ3 and 𝔵 m ∈ 𝒮 ≔ {{𝔷 1 𝔷 2 } T ; 𝔷 λ ∈ 𝒮 λ for λ = 1,2}. Denote the residual 𝔟 λ -𝔄(λ, : )(𝔵 {m-1} + 𝔵 m ) by 𝔯 λ m where the "colon notation" in MATLAB is used to designate the λth row of 𝔄 in 𝔄(λ, : ). The new mode 𝔵 m can be calculated from Eq.( 30) which is sometimes referred to as the nearest Kronecker product problem [START_REF] Golub | Matrix computations[END_REF].

𝔵 m = arg min 𝔵 m * ‖𝔯 λ m (𝔵 m * )‖ 2 for λ = 1,2 (30)
Or find 𝔵 m ∈ 𝒮 such that for λ = 1,2

(𝔄(λ, : )𝔵 m , 𝔵 λ m ′ ) 𝒮 λ = (𝔯 λ m-1 , 𝔵 λ m ′ ) 𝒮 λ ∀ 𝔵 λ m ′ ∈ 𝒮 λ (31)
where (⋆,⋆) 𝒮 λ is the inner product on 𝒮 λ and it is calculated as (⋆,⋆) 𝒮 λ : = ∏ (⋆ κ ,⋆ κ ) 𝒮 λκ 3 κ=1 [START_REF] Kolda | Tensor decomposition and applications[END_REF]. Since the vector ⋆ κ is complex valued, the inner product (⋆ κ ,⋆ κ ) 𝒮 λκ is ⋆ κ H • ⋆ κ in which H denotes the conjugate transpose. When 𝔄 is Hermitian, the above equation can be reformulated as

𝔵 λ m ′ H 𝔄(λ, : )𝔵 m = 𝔵 λ m ′ H 𝔯 λ m-1 (32)
Unfortunately, in our case the operator is symmetric but non-Hermitian because matrices in Eq.( 25) are symmetric and real-valued, terms involving multiplication of this matrix with the imaginary unit make 𝔄 non-Hermitian. Consequent-ly, Eq.(30) needs to be calculated using, e.g. the least square method. This is achieved by replacing 𝔄 in Eq.( 32) by 𝔄 ̃, being 𝔄 ̃: = 𝔄 T 𝔄. However, it is not recommended because first, 𝔄 ̃ becomes dense rather than sparse as 𝔄 is. Second, convergence rate of the resultant PGD system is significantly slowed down [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. Nevertheless, as reported in [START_REF] Chevreuil | A model order reduction based on proper generalized decomposition for the propagation of uncertainty in structural dynamics[END_REF], the approach based on Eq.( 32) proves efficient to capture good approximations in many cases even 𝔄 is non-Hermitian. We consider Eq.(32) in the following.

C. Implementing details

Solving Eq.( 32) is a nonlinear problem as products of unknowns are involved. With the ADM, it is assumed that the two components other than the μ th in the m th mode are known when calculating for 𝔵 λμ m , as in Eq.( 33). Assumed known variables need to be initialized where for the first step they can be arbitrary (or simply taken as ones) and for the subsequent, solutions of the former enrichment can be taken.

Recall that vectors associated with X, ω and r are calculated when μ takes value from 1 to 3 while mechanical displacements and electrical potentials are calculated when λ takes 1 and 2. This process enforces orthogonality between the residual 𝔯 λ m and the subspace 𝒮 λ . Iterations are involved using the ADM. Denote the iteration counter as p and the value of 𝔵 λμ m calculated at the p th iteration as 𝔵 λμ m,p . Iterations are terminated when either a maximum number n In for p is attained or the following convergence criterion is met

max λ ∏ ‖𝔵 λμ m,p -𝔵 λμ m,p-1 ‖ F ‖𝔵 λμ m,p-1 ‖ F 3 μ=1 < ε In ( 35
)
where ‖⋆‖ F stands for the Frobenius norm and ε In is a predefined criterion. It is noted that calculating 𝔵 λ m,p as 𝔵 λ1 m,p ⊗ 𝔵 λ2 m,p ⊗ 𝔵 λ3 m,p and then comparing the difference between 𝔵 λ m,p and 𝔵 λ m,p-1 is less practical as it is memory intensive and reduces the speed. Indeed, the relation

‖𝔵 λ m,p ‖ F = ∏ ‖𝔵 λμ m,p ‖ F 3 μ=1
indicates that it has the same effect as what described above.

The enriched vectors associating with the space and frequency are normalized before they are added into the new mode. Actually, one can normalize any two in the three vectors. The purpose is to prevent cases where one vector tends to infinity while another tends to zero, which causes numerical overflow. Normally, both ϵ In and n p are set to ensure that iterations do not take too long.

The enrichment proceeds till m > 𝑅 , the maximum enrichment number or the following criterion is met

max λ ‖𝔵 λ m ‖ F ‖𝔵 λ 1 ‖ F < ε out (36)
which means the last mode is small enough compared with the first mode. Eq.( 36) is rudimentary whereas it is efficient. Deciding the values for ε out is problem dependent.

Once the enrichment process is finished, an approximation 𝔵 pgd with the R rank-one tensors

{𝔵 λ 1 | ⋯ | 𝔵 λ R } (denot- ed as 𝔵 λ {R}
) to the real solution can be recovered using Eq.( 23). To evaluate the accuracy of 𝔵 pgd , its relative error ε rel w.r.t. the solution obtained using the conventional approach 𝔵 con ≔ {𝔵 1 con 𝔵 2 con } T can be calculated.

ε rel = max λ ‖𝔵 λ pgd -𝔵 λ con ‖ F ‖𝔵 λ con ‖ F (37)

D. Updating strategy

In [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF] the above process is referred to as the progress Galerkin PGD. Its convergence can be accelerated via updating the results after each enrichment. The idea is to impose orthogonality between 𝔯 λ m and the subspace (𝒮 λ ) m ≔ × μ=1,2,3 (𝒮 λμ ) m , as expressed in Eq.(38).

𝔵 λ l H (𝔟 λ -𝔄(λ, :

)𝔵 {m} ) = 0, l = 1,2, … , m (38)
The ADM can be employed which consists in solving all m components 𝔵 λμ l while fixing the other two set of m components 𝔵 λκ l for κ ≠ μ and take μ = 1,2,3 in turn. However, it is found in our numerical tests that the convergence rate can be extremely slow, especially for a large m. In other words, increasing iteration numbers has insignificant effects on the improvement of convergence. Hence, only one iteration in the updating is considered here. It also worthies noting that when fixing all frequency and resistor related vectors and updating space related vectors, the size of the system to be solved is the square of m × (N u + N φ ) that quickly becomes too large to be tractable as m increases. Therefore, vectors of a large number of dofs (spatial components) are not updated. The implementing is summarized in the algorithm chart.

IV. NUMERICAL EXAMPLES

A. Model description

Consider a cantilevered piezoelectric EH that is fixed onto a base. Acceleration are applied on the base which makes the EH deform, acting as mechanical excitation for the EH. To model the device, it is more convenient to consider displacements of the EH w.r.t. the base. Thus, we assume that the EH is fixed on the left while subject to an acceleration in the vertical direction (Fig. 4a). Dimensions of the bimorph are adopted from [START_REF] Zhao | Elecroelastic modeling and experimental validation of piezoelectric energy harvesting from broadband random vibrations of cantilevered bimorphs[END_REF] where two PZT-5H layers of 27 × 3.2 × 0.258 mm 3 and a brass substrate of size 27×3.2×0.115mm 3 are connected in series. Their material properties are summarized in Table I.

Six hundred excitation frequencies equally spaced between 0 and 3kHz are considered. The acceleration is taken constant (9.81 m/s 2 ) for all frequencies. To compare simulation results with measurements reported in [START_REF] Zhao | Elecroelastic modeling and experimental validation of piezoelectric energy harvesting from broadband random vibrations of cantilevered bimorphs[END_REF], a set of six resistors are chosen (being 100 Ω , 1 kΩ , 10 kΩ , 91 kΩ , 500 kΩ and 910 kΩ). An area representing the air is added surrounding the EH. In the mesh, the air area is set to be relatively larger than the EH domain. A total of 76824 tetrahedral elements with 14783 vertices are generated (Fig. 4b). Not that only half of the domain needs to be calculated due to the vertical symmetric surface in the middle.

Algorithm for parametric studies of the EH 1 :

Initialize 𝔵 λμ 0 with λ = 1,2; μ = 1,2,3

:

for i out = 1 to R do

B. Simulation results

In the first test, parameters for the PGD approach are set so that the fixed point iteration stopping criterion is 10 -3 while that of the enrichment process is 10 -5 . Maximum number of iterations and enrichments are set to 10 and 100, respectively (influence of different criterions are addressed in the next section). The same model is also solved using the conventional approach in which one loop for all samples of the resistor and the other for frequency are included. Both simulations are carried out in MATLAB on a workstation of 32G RAM and an Intel Xeon 3.5 GHz processor. The total consumed CPU time, memory required to store the solutions for both approaches and relative error between them are presented in Table II. The conventional approach is carried out both in serial and parallel (the latter is referred to as P-Conventional where the parfor function in MATLAB is used). For the PGD approach, both PGD with and without updating are tested. The former is referred to as U-PGD. PGD and U-PGD undergoes the same amount of enrichments and have solutions of the same size. Solutions of the conventional approaches are taken as the reference for evaluating accuracy of the PGD solutions, based on Eq.(37).

As can be observed, parallel implementation of the conventional approach shortens the calculating time to less than half of the serial counterpart, meanwhile CPUs are fully charged throughout, unlike in other cases where CPUs are only slightly charged. On the other hand, both PGD simulations take less time than P-Conventional in which U-PGD consumes less than half CPU time as the P-Conventional does (that is less than one quarter of the serial conventional version). Differences in CPU time between the PGD and conventional approaches can be explained from their algorithm complexities. Define 𝒪 S as the complexity to solve the FE system Eq.( 22) for one combination of ω and r. Thus, complexity of the conventional approach is N ω × N r × 𝒪 S . For the PGD approach, the most time-consuming part is calculating the space associated components since N ω and N r are much smaller than N u and N φ . Its complexity can be estimated as ∑ n ̃in i R i=1

× 𝒪 S where n ̃in i is the number of iterations in the ith enrichment. Therefore, the ratio of complexities of the two approaches is

N ω × N r ∑ n ̃in i R i=1 ⁄
. Larger this ratio is more time shall be saved with the proposed approach. Another advantage with the PGD approach is that only space of orders of magnitude less is demanded to store the simulation results (in our case the ratio of required space for the conventional and PGD approach is (N u + N φ ) × N ω × N r versus (N u + N φ + 2N ω + 2N r ) × R ). It is also observed that the updating plays an important role in accelerating convergence rate of the PGD (see Fig. 5) which leads to more accurate solutions. In Fig. 6 and Fig. 7, space and frequency components for the first eight modes are depicted, respectively. As can be seen, only bending modes (none torsion nor other high order modal shapes) appear in the space components which is due to the fact that frequencies chosen in our example only cover the first two eigen frequencies of the bimorph. For the frequency components, most peak values locate near the eigen frequencies of the system. In Fig. 8 are presented impedance components for 𝑅 = 91𝑘Ω whose counterparts for other values show a similar trend. For the validation of our results obtained with 100 enrichments, voltages between electrode 1 and 4 are extracted from the PGD results. This voltage under different frequencies and resistors are shown in Fig. 9a. It is found that our simulated results compare favorably with the experimental results in the literature (c.f. Fig. 4a in [START_REF] Zhao | Elecroelastic modeling and experimental validation of piezoelectric energy harvesting from broadband random vibrations of cantilevered bimorphs[END_REF]). In Fig. 9b, the harvested power, or power consumed on the resistors are depicted. As can be seen, the resistor of 91kΩ is optimal.

C. Influence of the algorithm parameters

Additional simulations (U-PGD version) are carried out in order to investigate the influence of different parameters in the novel approach. Results are grouped and presented in tables, one for the maximum loop number and stopping criterion of the fixed point iterations (referred to as inner loops) and the other for enrichments (referred to as outer loops). Results are compared in terms of three items: the total iteration number n ite , the convergence indicator in the last enrichment 𝜀 conv calculated with Eq.(36), and the relative error 𝜀 rel obtained with Eq.(37). It is to be noted that only the order of magnitudes are of concern, because the convergence is not monotone and the values may locally oscillate.

In Table III, the maximum enrichment number is fixed at 100 and the convergence criterion is set to 10 -5 . It is found that increasing the maximum iteration number n IN does not significantly improve the accuracy while the complexity increases dramatically. It is recommended to take n IN as small as possible e.g. 5 or 10 whenever possible. However, it is more tricky with ε IN . When it is too large, solutions in step 4 of the algorithm can be far from the exact ones for which more enrichments are required to keep the approximation precision. However, if it is too small the results can be deteriorated, because more iterations are needed for small ε IN 's and it can diverge after certain iterations. Consequently, the solved modes are not necessarily more accurate albeit the larger number of iterations.

In Table IV, the maximum iteration number is 10 and the iteration convergence criterion is 10 -3 for all cases. Influences of the outer loop parameters are more evident. As suggested in Eq.( 23), a larger 𝑅 normally leads to a better approximation to the exact solution. Therefore, both smaller ε OUT and larger n OUT can hopefully lead to better approximations. However, it is important to realize that the limit of the approximation precision is decided by various factors, such as parameters of both inner and outer loops, and the updating strategy. As a result, the results are no longer improved to some limit by increasing the number of modes. In our case, a maximum number of 50 at 10 -5 is found to be appropriate.

V. SUMMARY AND CONCLUSIONS

A PGD based approach for the parametric study of EH numerical models is proposed in this paper. Compared with conventional approaches, the novel approach is capable to obtain numerical results of almost the same quality while at a much reduced complexity. Using the proposed approach, CPU time for solving parametric problems is effectively reduced along with less required space to store the results. In addition, the proposed approach can be readily implemented. Nevertheless, there is still much space for improvements of the proposed method such as effective methods to solve high dimensional problems in updating and more rigorous error estimators for convergence check. Additionally, solving the three dimensional spatial components dominates the simulation process which can be reduced by PGD based decomposition of the full spatial fields into lower dimensional ones [START_REF] Vidal | Proper generalized decomposition and layer-wise approach for the modeling of composite plate structures[END_REF] though extra efforts than in a pure structural problem are needed due to the piezoelectric coupling. Such problems offer tantalizing prospects for future extensions of this work.

  with Lipschitz boundary Γ φ = Γ φ D ∪ Γ φ N where Γ φ D and Γ φ N are the Dirichlet and Neumann electrical boundary conditions. The mechanical domain Ω u , which is a subdomain of Ω φ (Ω u ⊂ Ω φ ), has its Lipschitz boundary Γ u = Γ u D ∪ Γ u N , Γ u D and Γ u N being, respectively, the Dirichlet and Neumann mechanical boundary condition.

Fig. 1 .

 1 Fig. 1. Electric and mechanical domains of the multiphysics problem

  13.b) and for the β th element e β ∈ Ω φ ∕ Ω u -𝐊 φ (β) Φ β = 0 (13.c) where U α and Φ α (resp. Φ β ) represent all nodal values within the element e α (resp. e β ). Without losing generality (w.l.g.), assume that linear tetrahedral elements are used. Then, the vector Φ α equal to[φ is of size 4 × 1 and U α is of size 12 × 1, α1 to α4 being node numerations in element α. The mass matrix 𝐌 u (α) , stiffness matrices 𝐊 u (α) and 𝐊 φ (α) , coupling matrix 𝐊 uφ (α) and the mechanical load vector F X (α) are calculated as follows. β) = ∫ 𝐁 φ e α ,T [ε S ]𝐁 φ e α Ω e α,β dΩ (14.e) F X (α) = ∫ 𝐍 u e α ,T ρ{f} Ω e α dΩ + ∫ 𝐍 u e α ,T ρt̅ Γ u N ∩∂Ω e α dΓ (14.e) in which the matrices are defined as

  from m e + 1 to m e + n φ -∑ n e k m e k=1 -n φ ̂, n φ ̂ denoting the number of nodes on Γ φ D . Meanwhile, the mechanical mapping ℳ dof u can be such defined that the j th component of the dofs at the l th mesh node is assigned 3 × l + j -1 with j = 1,2,3. Hence, the total number of mechanical and electrical dofs N u and N φ are respectively, 3n u -n u ̂ and m e + n φ -n φ ̂-∑ n e k m e k=1

Fig. 2 .

 2 Fig. 2. Piezoelectric EH connected to electrical loads in series Denote 𝐊 ̆φx ̌ as the product V x ̌Vx ̌T (the incident matrix) for conciseness. Then the matrix -j ωZ ̃x ⁄ 𝐊 ̆φx ̌ is added to 𝐊 ̂φ. 𝐊 ̂φ = ∑ j ωZ ̃x ̌𝐊 ̆φx x ̌

Fig. 3 .

 3 Fig. 3. Approximation of the solution via low rank tensors

𝔄 22 =

 22 𝐊 φ ⊗ 𝐃 ω ⊗ 𝐃 r -j𝐊 ̆φ ⊗ 𝐆 ω ⊗ 𝐆 r(25.d)

  (33) into Eq.(32) and cancelling 𝔵 λμ m ′ on each side yields a linear system for 𝔵 λμ m

Fig. 4 .

 4 Fig. 4. Model configuration (a) the cantilevered piezoelectric EH, (b) mesh

Fig. 5 .

 5 Fig. 5. Convergence of enrichment for with and without updating

Fig. 6 .Fig. 7 .

 67 Fig. 6. Space components of the 1st to the 8th enrichment (displacements are depicted by the deformation from the original positions while potentials are represented by the color; the area for air is not shown)

Fig. 8 .

 8 Fig. 8. Impedance components for R = 91kΩ Fig. 9. Output of the energy harvester (a) generated voltage, (b) power consumed by the electrical load

TABLE III STUDY

 III ON PARAMETERS OF THE INNER LOOPS (N OUT = 100, Ε OUT = 10 -5 )

			𝑛 IN = 5			𝑛 IN = 10			𝑛 IN = 20			𝑛 IN = 30
	𝜀 IN	𝑛 ite	𝜀 conv × 10 -5	𝜀 rel × 10 -5	𝑛 ite	𝜀 conv × 10 -5	𝜀 rel × 10 -5	𝑛 ite	𝜀 conv × 10 -5	𝜀 rel × 10 -5	𝑛 ite	𝜀 conv × 10 -5	𝜀 rel × 10 -5
	10 -2	480	2.2	7.2	624	2.0	6.8	608	0.82	5.6	609	0.97	5.6
	10 -3	495	3.7	7.8	839	3.9	8.0	798	0.88	5.8	999	0.98	5.8
	10 -4	498	2.0	6.7	937	2.2	6.5	1052	0.95	6.1	1575	0.99	6.0
	10 -5	500	52	44	863	0.96	7.7	1119	0.97	5.6	1706	0.88	5.6

TABLE IV STUDY

 IV ON PARAMETERS OF THE OUTER LOOPS (N IN = 10, Ε IN = 10 -3 )

			𝑛 OUT = 30		𝑛 OUT = 50		𝑛 OUT = 100		𝑛 OUT = 200
	𝜀 OUT	𝑛 ite	𝜀 conv × 10 -2	𝜀 rel × 10 -2	𝑛 ite	𝜀 conv × 10 -3	𝜀 rel × 10 -3	𝑛 ite	𝜀 conv × 10 -4	𝜀 rel × 10 -4	𝑛 ite	𝜀 conv × 10 -5	𝜀 rel × 10 -4
	10 -2	228	0.72	4.8	228	7.2	48	228	72	480	228	720	480
	10 -3	254	5.7	5.4	417	2.7	1.9	425	8.2	18	425	82	18
	10 -4	254	5.7	5.4	417	2.7	1.9	623	0.92	3.6	623	9.2	3.6
	10 -5	254	5.7	5.4	417	2.7	1.9	835	0.17	0.69	948	0.82	0.58


Nicolas Goubet
email: nicolas.goubet@upmc.fr

P A Albouy

A Thompson

M P Pileni

Polymorphism in nanoparticle-based crystals depending upon their single or polycrystalline character

published or not. The documents may come

Polymorphism in nanoparticle-based crystals depending upon their single or polycrystalline character

Introduction

There is no doubt that the research on nanoparticle (NPs) selfassembly is coming to a new stage with the challenge of the growth of single-crystals of NPs large enough to allow easy manipulation. This step is determinant as the precise knowledge of the structure and symmetries of these objects open the door to a deeper study of their physical properties. 1-7 Indeed the arrangement of NPs coated with organic ligands give rise to a new class of functional materials where specific properties, different from those of the NP or of the corresponding bulk phase, can be expected. [8][START_REF] Pileni | Nanocrystals Forming Mesoscopic Structures[END_REF][10][11] It must be insisted that these single-crystals are periodic crystals as defined by the IUCr 12 and need to be distinguished from the recently introduced notion of a mesocrystal, where orientational order, and not translational order, is essential.

13- [START_REF] Gur | [END_REF] We shall keep to the term "colloidal crystal" instead of "superlattice" to designate them, albeit this last is often found in the literature; however it has a specific and quite different meaning in crystallography. A challenge for the development of technologies based on these materials is to achieve a precise control of the self-assembly process. The presence of the translational order expands the variety and complexity of the mesoscopic structures made of one or multiple components. [16][17][18][19][20][21][22][23][24][25] In case of spherical NPs with low size dispersion, face-centered (fcc) or body-centered cubic (bcc) lattices are most often observed. Hexagonal close packing (hcp) is more seldom reported although the difference in free energy compared to fcc is not expected to be large. [26][27][28] On the other hand, the free energy balance between fcc and bcc structures is mainly controlled by the gold volume fraction Φ Au and the ratio χ between the ligand length and the radius of the inorganic core. 29,30 It is generally accepted that for a given Φ Au value, the bcc structure is favored above some threshold for χ. Below this threshold, the NPs are assimilated to hard spheres with short-range interactions that self-assemble into the denser fcc packing. However, several studies report the unexpected observation of the bcc structure. Choi et al. have reported that oxidized PbS NPs with anisotropic ligand coverage self-assemble into the bcc structure, whereas PbS NPs with an isotropic coating pack on a fcc lattice. 31 Solvation processing of the ligand shell have shown to be a good way for the control of the bcc/fcc polymorphism control in assemblies of PbSe NPs. 32 Inclusion of either extra organic molecules or ligand excess can also induce the fcc to bcc transition in PbSe or Au NPs colloidal crystals. 33,34 These emerging parameters have thus to be taken into account. Following the same rationale, the involvement of orientational ordering between NPs in colloidal crystals, as encountered in mesocrystals, is still unclear. Indeed, contrary to semiconductors, Au NPs can be grown as single or polycrystals due to the easy generation of planar defects; well-ordered colloidal crystals can be grown by self-assembly in both cases.

In this paper, we elaborate sub-millimeter colloidal singlecrystals by self-assembling either gold truncated octahedral single-crystal NPs (thereafter denominated nanocrystals, NCs)

or polycrystalline decahedral and icosahedral NPs also known as Multiply Twinned Particles (MTP). 35 Their average size is ca.

5nm in both cases. The high X-ray scattering power of gold atoms combined with the brilliance of synchrotron sources allowed the acquisition of high-quality diffraction patterns from individual colloidal single-crystals. To the best of our knowledge, only a few studies are reported based on such a technique, mainly due to difficulties in growing and handling large and isolated samples. [36][37][38] It is thus possible to reconstruct in great detail a large portion of the reciprocal space associated with these objects. This approach allows to access important structural details as, for instance, intrinsic disorder, that are hard to reveal when samples are under the form of a crystalline powder as usually reported. [18][19][20]32,39 This is presently exemplified by the observation of diffuse scattering outside the Bragg reflections.

Experimental

Materials chlorotriphenylphosphine gold (I) and tert-butylamine borane were obtained from STEM chemicals. Dodecanethiol, ethanol were obtained from Sigma-Aldrich Chemical Co. Toluene (98%) from Riedel de Haen. Ethanol (99.8%) from Prolabo. All reagents were used as received without further purification.

Synthesis of 5,2 nm gold nanocrystals

The gold nanocrystal were synthesized by revisiting a method from the Stucky's groups. 40 The solution of gold precursor is made with 124 mg of chlorotriphenylphosphine gold (I) dissolved in 25 mL of toluene with 500µL of dodecanethiol. The reducing solution is made of 330 mg of tert-butylamine borane complex dissolved in 5 mL of toluene. The two solutions are heated up to 100°C in a silicone oil bath. The tertbutylamine borane solution is then quickly added to the gold precursor solution. The solution changes from colorless to dark red. After 5 minutes of reaction, the solution is let to cool down to room temperature. After evaporation of the solution, the gold nanocrystals are then washed with ethanol and redispersed in toluene. The resulting colloidal solution consists in a mixture of single crystals and polycrystals as shown in Fig. S1.

Formation of colloidal crystals with gold single crystals

Ethanol (3% v/v) is added to the washed colloidal solution at [Au]=2. 5 10 -6 M and 500µL are poured into a NMR tube with a glass capillary inside. In order to prevent toluene evaporation, the NMR tube is confined inside a larger tube with toluene at the bottom. After one week, some shinny reflections can be seen on the glass wall of the NMR tube that reveals presence of micrometric crystals.

Formation of colloidal crystals with gold polycrystals

The remaining colloidal solution is poured into a new NMR tube and confined in a larger tube with 2 mL ethanol at the bottom. The solution is let to destabilize for one week. The solution becomes colorless and shinny reflections appear on the wall of the tube and capillary.

Results and discussion

The colloidal crystals are obtained by self-assembly of gold NPs coated with dodecanethiol; these NPs are synthesized by degradation of the organometallic precursor triphenylphosphine gold (I) in toluene. 40 The as-synthesized NPs reveal a size distribution below 5% that is mandatory for the observation of 3D long-range order. However, their single or polycrystalline character could not be controlled using synthesis parameters and dark field pictures reveal a mixture of both types (see ESI). They were separated using the crystallinity segregation already reported in previous works. 39,41 However, contrary to these previous work, the NPs are too small to induce a spontaneous self-assembly at room temperature. 42 Consequently, the NPs mixture is redispersed into toluene and a first crystallization process is induced by adding a small amount of bad solvent (ethanol) to the solution (see Experimental Section). Solution evaporation is prevented by placing the solution in an enclosure saturated with toluene vapor. Addition of a polar solvent is preferred to induce ordered self-assembly compare to lowering the temperature, which may prevent crystallization as reported elsewhere.

42

After one week, the solution remains red-colored but shiny reflections can be seen on the glass walls and at the bottom of the tube. This is indicative for the growth of large NPs-based single-crystals (hereafter denominated type I colloidal crystals) as confirmed by the optical micrograph displayed in ESI. The solution with remaining free NPs is poured into another tube and placed in an enclosure saturated with ethanol vapor. The slow (one week) diffusion of this bad solvent into the colloidal solution induces NPs precipitation. The process is made visible by the progressive discoloration of the solution that starts from the top and propagates to the bottom, while some shiny reflections become visible at the bottom and on the glass wall. These last colloidal crystals are denominated type II. SEM pictures of the as-grown colloidal crystals show that their shape is mostly distorted octahedron for both types (Fig. 1 ) for a colloidal crystal grown from polycrystalline gold nanoparticles (type II). The planes are reconstructed based on 720 diffractions patterns obtained in the rotating crystal geometry. levels could be observed using the protein crystallography beamline PROXIMA1 of the synchrotron facility SOLEIL. Samples were glued at the tip of a glass capillary by a small amount of silicon grease and placed on the φ-rotation stage of the goniometer. 720 diffraction patterns were taken every 0.5° with an exposure time of 1s per frame. In order to access to both levels of periodicities, the measure was repeated for two different sample-to-detector distances and wavelengths without changing the crystal orientation: 850mm and 0.1042nm for the detection of smaller scattering vectors, and 196mm and 0.0867nm for larger ones. The choice of wavelengths is dictated by the necessity of minimizing fluorescence and absorption effects. A common orientation matrix can thus be used for both sets of data. This last is determined by the software CrysAlisPro 43 based on the reflections associated to the colloidal crystal periodicities. The subsequent reconstruction of the reciprocal space and data reduction are performed with the same software. Data relevant to the NPs packing are discussed first. We show representative reconstructed reciprocal planes in Fig. 3 and4 for type I and II respectively: the large numbers of well-defined Bragg reflections (50 to 100) is a sign for the good quality of both types of colloidal crystals. The ellipses and straight lines originate from the inactive inter-chips zones of the Pilatus 6M detector and of the beam-stop holder. The indexation of the Bragg reflections indicates that both colloidal crystals belong to the cubic system. Condition for reflection is h+k+l = 2n for type I, that points to a bcc direct lattice (fcc reciprocal lattice).

Conversely condition for reflection in type II is similar parity for h, k, l indices: this means that the direct lattice is fcc (bcc reciprocal lattice). Apparently forbidden reflections visible in Fig. 4-c (white hexagons) are in fact diffuse scattering features that are discussed below. Three samples of type I with better quality were selected for data analysis; they have very similar cubic parameters with an average value a = 8.02 ± 0.05nm. The two samples of type II with better quality that have been selected display somehow different cubic parameters, 9.79nm and 10.45nm; it is shown below to correlate with slightly different NP size. Diffraction features associated to the NP inner structure are observed at larger scattering vectors and discussed below (see Fig. 5).

The case of type I samples is first discussed. Patterns in Fig. 5 Unexpectedly, the colloidal crystals type I are distorted octahedra ({111} CC facets), which is not the more stable crystal habit for a bcc lattice (rhombic dodecahedron with {110} facets). 44 However, the co-axial alignment between NCs and colloidal crystals makes favorable the stabilization of colloidal crystal facets made with the larger NC facets, i.e. {111} NC facets. A close examination of the {111} NC diffraction line in Fig. 6-a further reveals well-defined secondary maxima; their observation is made possible by the very low NC size dispersion within the colloidal crystal. They are clearly distinguished in the radial scan plotted in Fig. 6-b. Fitting the line profile gives valuable information on the shape and size of the NCs. Based on HRTEM images, a NC is first modelled by a cube originally containing N 3 fcc unit cells and subsequently truncated so as to grow symmetric {111} NC facets. For this purpose, nodes of coordinates (u,v,w) such that ±𝑢 ± 𝑣 ± 𝑤 > 𝑚 are removed; particular values are 𝑚 = 𝑁 for cuboctahedron, 𝑚 = 3𝑁/4 for truncated octahedron and 𝑚 = 𝑁/2 for an octahedron. As noticed above, the statistical misorientation of NCs crystallographic axis with respect to those of the host colloidal crystal may be represented by a 2D-Gaussian function with half-width at half-maximum 10°; the Box-Muller method is most convenient to create orientations that follow this type of statistics. 45 Both levels of periodicities are incommensurate and the amplitudes scattered by the NCs present no phase relation at large angles; in these conditions the intensities corresponding to the accessible NC orientations simply add. The atomic scattering factor of gold is taken from the International Tables for Crystallography. [START_REF]International Tables for X-Ray Crystallography[END_REF] radius of 2.59 nm containing 4196 gold atoms (see Fig. 6-b, blue curve). The number of gold atoms within the NC is actually quite similar for both shapes. We now turn to the case of type II colloidal crystals composed of multiply twinned particles (MTP). We see in Fig. 5-d that the diffraction lines are now homogeneous which indicates a perfectly random orientation of the coherent domains; the corresponding radial scan for the scattered intensity is shown in Fig. 6-c. In order to determine the location and breadth of the main peak and of the shoulder, the profile has been adjusted by two Gaussian functions. The location of the main peak at 26.93nm -1 and corresponds to a cubic parameter 0.404nm for a (111) NP assignation. The broadening and shift towards smaller angles for the (200) NP diffraction line is typical for decahedral or icosahedral NPs that can be viewed as twinned assemblies of defective fcc tetrahedral subunits. [START_REF] Hall | [END_REF] For the sake of simplicity, we simply estimate a NP size p defined as the cube root of the coherent volume using the Scherrer formula 𝑝 = 𝐾𝜋/Δ𝑞 1/2 . Here Δ𝑞 1/2 = 1.04 𝑛𝑚 -1 is the (111) NP peak half-width at half-maximum and K is taken equal to 0.8894, as found in the literature for the (111) NP diffraction line in tetrahedral crystallites. 48 The value p=2.62nm gives an estimate of 18nm 3 for the coherent volume. We show below that the NPs are close to spheres of radius 2.74nm for this particular sample. It corresponds to 4.8 coherent volumes per particle close to the value of 5 expected for decahedral structures. However, no definite conclusion can be drawn for the proportion of decahedral and icosahedral NPs, as tetrahedrons in icosahedron are expected to be much smaller, and thus to give a much more diffuse signal. The intensities of Bragg reflection associated to the NP lattice are now analyzed for the two types of colloidal crystals. As mentioned above, an automatic data reduction is performed by the software CrysAlisPro; however, the listed intensities need to be inspected to reject data clearly impacted by the proximity of the inactive zones of the detector during acquisition. A reliability factor is defined as 𝑅 𝑓 = ∑||𝐹 𝑜𝑏𝑠 | -|𝐹 𝑐𝑎𝑙𝑐 || / ∑|𝐹 𝑜𝑏𝑠 | where F obs and F calc respectively refer to the observed and calculated amplitudes. A Debye-Waller-like correction 𝑒𝑥𝑝(-𝑈 𝑖𝑠𝑜 𝑞 2 /2) is further introduced; indeed its main purpose is to account for possible static shifts of NPs out of their equilibrium position (by an average amplitude √𝑈 𝑖𝑠𝑜 ) as temperature is shown below to have no visible impact on the Bragg intensities. Due to the shape factor of NPs, the intensity of Bragg reflections falls rapidly with increasing scattering vector and signals typically disappear into the background for q  8nm -1 for type I and q  6nm -1 for type II.

For the sake of simplicity only the 26 first reflections are included in the fitting procedure (q  5nm -1 for both colloidal crystals). We first consider the data analysis of type I colloidal crystals. Contrarily to the previous case, NPs diffract coherently in this scattering range and averages on diffracted amplitudes and not intensities have to be considered. They are calculated considering the different possible NCs orientation as explained above. Following the same lines, {111} NC truncated cubes are first tested. A best agreement is found for N=13 and m=9.5 (4794 atoms; R f = 0.17 and U iso = 0.02nm 2 ) that correspond to a distance of 5.32nm between {100} NC opposite facets and 4.49nm between {111} NC opposite facets. It is to notice that the distance between {100} NC opposite facets is slightly bigger than the one deduced previously. However a significantly better reliability factor (Rf = 0.13) is obtained with a spherical shape, as already seen in the treatment of the larger angle data (see Table 1). This probably reflects some rounding up of the NCs edges. The as-determined radius is 2.63nm close to 1. The corresponding amplitudes are expected similar for a perfectly spherical shape, which is not completely verified: this probably reflects a still present faceting. The Debye-Waller factor value corresponds to a rather small out-of-equilibrium average displacement: √𝑈 𝑖𝑠𝑜 = 0.14 𝑛𝑚. The minimum inter-particle gap is along the longer diagonal and its determination is rather model-dependent: 2.45nm for the truncated cube and 1.69nm for the sphere. Data for type II colloidal crystals are now discussed. In accordance with HRTEM images, a spherical particle shape is directly chosen. The procedure has been applied to the two above-mentioned best samples with slightly different cubic parameters. As an example, data for sample with cubic parameter 10.45nm are given in Table 2. A good agreement between measured and calculated amplitudes is found for a particle radius 2.74nm (R F =0.04); it is significantly better than for type I. The amplitudes of reflections with similar d-spacing are further rather close so that a spherical model is really pertinent. As expected, the particle radius for the second sample with smaller cubic parameter 9.79nm is found smaller 2.59nm (U iso = 0.1nm 2 ; Rf = 0.08); indeed the minimum interparticle gap found along the shorter diagonal is rather similar for both samples: 1.91nm and 1.74nm. We now discuss the possible origin of the difference in packing mode observed for NCs and MTPs. It has been shown in the literature that the presence of ligand excess can induce polymorphism in colloidal crystals. 33,34 However in the present study, the gold colloidal crystals are made with the same colloidal solution, i.e. same ligand excess. We mentioned in introduction that a large corpus of experimental data and theoretical work enlightens the importance of the gold volume fraction Φ Au and the ratio χ between the ligand length and the radius of the inorganic core; in particular a (Φ Au , χ) curve is determined that separate fcc and bcc as the most stable packing modes. 29,30 Indeed, Φ Au and χ as evaluated with the assumption of a spherical particle shape for both types of colloidal crystals and using the X-ray data, and are found very similar: Φ Au =0.3 and χ =0.59 (with the hypothesis of an extended length for a ligand with an alkyl chain of n carbons 0.12(n+1)=1.56nm). 49 The corresponding point in the (Φ Au , χ)

diagram lies within the domain for expected fcc packing. Unexpectedly, type I colloidal crystals exhibit a bcc lattice instead of the expected fcc which shows that the softness of the ligand shell alone is not enough to predict the packing. Indeed, the nanocrystallinity influences the packing process via two major effects. The first is related to the nanocrystal shape. Indeed, the parameters for the cube truncation (N=13 and m=9.5 as found by the fitting procedure) are very close to those of a truncated octahedron (m=9.75) that constitutes the Wigner-Seitz cell of the bcc structure; bcc would be the most efficient packing for NPs with such a shape. Another parameter that can favor such a face-to-face packing is the spatial distribution of the ligands on the NCs. Indeed, we know that the dominant interaction between the nanoparticles comes from van der Walls forces between the ligand bundles fixed on their surface. 50 Contrary to PbS NCs, 31 the potential functions for the S-Au bond given shows that the {100} and {111} facets have similar S-Au bonding behavior, so that the coverage is quite similar for both surfaces. 51 However, this is not affecting the bundle formation 51 and one can expect that the organic shell replicates the nanocrystal shape. In the case of single crystals, the truncated octahedron is the Wigner-Seitz cell of the bcc lattice. However in the case of polycrystals, fivefold symmetry of icosahedron and decahedron is incommensurate to the translational order of the colloidal crystals, which leads to polycrystals behaving like spheres. One can notice that the bcc/fcc polymorphism was not observed in In their study, the NCs and polycrystals self-assemble into fcc and fcc/hcp lattice mixture, respectively. This only modest influence of nanocrystallinity can be due to the method for NPs deposition, as evaporation of the solvent is known to induce compression of the NPs lattice and favor a compact structure.

4

An analysis of the Bragg reflection intensities only yields an average structure for the colloidal crystal (lattice parameters, average NP shape). Diffuse scattering results from any difference between the real crystal structure and the average one; 52 the most simple example is given by thermal vibrations that displace atoms from their average position, as it was observed in materials such as silicon, ice or tin polymorphs. [START_REF] Holt | [END_REF][54][55] Moreover, the presence of acoustic phonons in such NC colloidal crystal has been revealed in recent work by pumpprobe spectroscopy. 56,57 In the present case, diffuse lines extending between Bragg reflections are clearly visible on diffraction patterns from both colloidal crystals (see Fig. 3, 5 and 7). We show in Fig. 7 successive cuts in the reciprocal space with increasing l-index that allow a spatial reconstruction of the diffuse scattering. Animations are provided in ESI that displays this exploration in a vivid and precise way. In case of type I colloidal crystals with bcc packing, diffuse scattering takes the form of broad diffuse rods passing through the Bragg reflections and pointing in the eight {110} CC directions: this is schematically represented in Fig. 8. The observation of rods is typically correlated to some planar disorder that should affect the (110) CC denser planes in real space. [START_REF] Guinier | X-ray Diffraction in Crystals, Imperfect Crystals, and Amorphous Bodies[END_REF] No diffuse scattering is furthermore detected close to the center of the reciprocal space, which is in favor of a displacement disorder. This disorder can be static or dynamic, which is easily probed comparing similar diffraction planes recorded at different temperatures. For this purpose, X-ray data were collected on a same sample at room temperature and at 110K; cooling was obtained by gently blowing cold nitrogen gas on the sample using an Oxford Cryostream 700 cooling device devoted to protein crystallography (see ESI). No significant differences were observed due to this change in temperature, which indicates that the disorder is static. Indeed collective motion of NCs requires at least some thermal activation of the organic capping agent and this is not expected to occur for temperatures below typically 350K. It can be assumed that the truncated octahedrons interact more strongly through their facet-to-facet contacts within the denser {110} CC planes. In these conditions we propose the existence short-range orientational correlations between neighboring NPs within a given {110} CC plane; this is exemplified in Fig. 8-c where a perfect orientation of the central NP with respect to the colloidal crystal axis is chosen for the sake of simplicity. However the intensity distribution along the diffuse lines appears rather homogeneous so that these orientational correlations should not extend to neighboring NPs located on adjacent {110} CC planes (no explanation is proposed for this last point).

We now turn to the case of type II colloidal crystals. Diffuse scattering consists in well-defined lines that extend from each Bragg reflection in the four [111] CC directions, except for those directions that point toward the reciprocal space center. Apparently forbidden well-defined reflections are clearly visible on the plane of equation ℎ + 𝑘 + 𝑙 = 0 (see hexagons in Fig. 4-c). In fact they correspond to intersections at the points (h+1/3,k+1/3,l+1/3) of those lines that run along the [111] CC direction and that join a h,k,l Bragg reflection of the plane ℎ + 𝑘 + 𝑙 = -1 to the next (h+1,k+1,l+1) reflection of the planeℎ + 𝑘 + 𝑙 = 2. The presence of these apparently forbidden reflections was also reported in observations by electron diffraction of Ag and Au plate like crystals and by Xray diffraction on opals. [58][59][60] The diffuse lines may be interpreted as being due to stacking faults that are known to be frequent in fcc NPs assemblies. [START_REF] Warren | X-Ray Diffraction[END_REF] Due to the isotropic growth of the colloidal crystal in solution, no specific [111] faulting direction is favored and all four directions are equally represented. Indeed, by mapping the reciprocal space perpendicularly to the four compact directions, each reconstructed slice shows forbidden reflections (see ESI).

Conclusions

In summary, large (up to 50µm) and high-quality well-faceted NPs colloidal crystals could be grown in solution from either single-or polycrystalline particles; the average particle diameter was quite similar for both cases (ca. 5.2nm in diameter). Their high crystallographic quality allowed us to collect diffraction data from single samples, opening the door to detailed structural investigations. In particular, weaker diffuse scattering features outside Bragg reflections have been evidenced in both mesostructures and have been related to orientational or packing disorder. With similar gold volume fraction and the ratio of coating length to nanoparticle radius, the packing mode for these colloidal crystals exhibits a drastic sensitivity to the single or polycrystalline character of NPs: selfassembly of polycrystalline particles leads to a compact fcc structure with stacking faults distributed along the four [111] CC stacking axes, whereas single-crystal particles self-assemble in a bcc structure. In this last case, a very pronounced preferential orientation of the NCs crystallographic axis along those of the colloidal crystal is observed. This observation is related to the faceting of the single-crystal particles clearly observed by HRTEM. In a more general perspective, this study is a further demonstration for the possibility to record high-quality X-ray diffraction data (including diffuse scattering) from isolated NPs colloidal crystals. Indeed, a detailed knowledge of the internal structure of such materials is a prerequisite for the comprehension of their physical properties.

from EPFL, Switzerland for very helpful discussions on the use of the CrysAlisPro software. This work was carried out at the PROXIMA1 beamline at SOLEIL under proposal n° 20130947. Financial support was obtained from the advanced grant of the European Research Council under n° 267129.

Fig. 1

 1 Fig. 1 (a) View of a colloidal single-crystal made of singlecrystal gold nanoparticles (typeI), (b) View of a colloidal singlecrystal made of polycrystalline gold nanoparticles (typeII).

Fig. 2 (

 2 Fig. 2 (a and d) Low magnification TEM images (scale are 20 nm), (b and e) diameter distribution and (c and f) HRTEM images for particles obtained after dissolution of a colloidal crystal into a solvent and drop deposition on a carbon grid (scale bar are 2 nm). Images a and c (resp. d and f) correspond to particles obtained after dissolution of a few colloidal crystals of the first (resp. second) type.

Fig. 3 .

 3 Fig. 3. Reconstruction of the reciprocal planes of equation l = 0 (a), h+k = 0 (b) and h+k+l = 0 (c) in the smaller q-value zone (q max = 7 nm-1) for a colloidal crystal grown from single-crystal gold nanoparticles (type I). The planes are reconstructed based on 720 diffractions patterns obtained in the rotating crystal geometry.

Fig. 4

 4 Fig. 4 Reconstruction of the reciprocal planes of equation l = 0 (a), h+k = 0 (b) and h+k+l = 0 (c) in the smaller q-value zone (q max = 7 nm -1

Fig. 5

 5 Fig. 5 Reconstruction of the reciprocal planes of equation l = 0 (a), h+k = 0 (b) and h+k+l = 0 (c) in the larger q-value zone (q max = 42 nm -1 ) for a type I (a, b and c) and type II (d) colloidal single-crystal. The planes are reconstructed based on 720 diffractions patterns obtained in the rotating crystal geometry.

Fig. 6 (

 6 Fig. 6 (a) Oscillations surrounding an enhancement in the 111 Debye-Scherrer ring of a type I colloidal crystal. (b) Adjustment for the radial intensity profile by different particle shapes. For the sake of clarity the square root of the intensity is plotted in ordinate.

Fig. 7

 7 Fig. 7 Cuts of the reciprocal space for different elevations of index l are used for a reconstruction of the diffuse scattering. These snapshots are taken from an animation provided in ESI. (a) type I colloidal crystal, (b) type II colloidal crystal.

Fig. 8

 8 Fig. 8 Representation of the various features observed in the reciprocal space for the type I (a) and II (b) of colloidal crystals. The blue lines represent the diffuse scattering located between the Bragg nodes (red dots). Packing model of the single crystal (c) and polycrystals (d).

Table 1

 1 Comparison of measured and calculated (spherical particle shape) amplitudes for type I colloidal crystals. Indices marked with the same color are associated to reflections of similar d-spacing.

	h k l	│F obs │ │F calc │ h k l	│F obs │ │F calc │
	1 1 0 819	853	4 3 1 24	22
	2 0 0 186	145	5 1 0 12	22
	2 1 1 195	135	5 2 1 18	16
	2 2 0 262	190	4 4 0 24	28
	3 1 0 148	145	4 3 3 38	35
	2 2 2 77	71	5 3 0 30	35
	3 2 1 14	4	4 4 2 36	37
	4 0 0 60	45	6 0 0 29	37
	4 1 1 74	68	5 3 2 32	34
	3 3 0 68	68	6 1 1 25	34
	4 2 0 70	73	6 2 0 20	28
	3 3 2 74	62	5 4 1 18	21
	4 2 2 48	43	6 2 2 10	12

the previous determination. Reflections associated to similar d-spacing are pointed to by similar colors in Table

Table 2

 2 Comparison of measured and calculated (spherical particle shape) amplitudes for type II colloidal crystals. Indices marked with the same color are associated to reflections of similar d-spacing.

	h k l	│F obs │ │F calc │ h k l	│F obs │ │F calc │
	1 1 1 744	758	6 0 0 24	24
	2 0 0 509	492	4 4 2 23	24
	2 0 2 46	44	6 2 0 9	10
	1 1 3 147	136	5 3 3 4	0
	2 2 2 145	139	6 2 2 4	3
	4 0 0 86	90	4 4 4 7	10
	3 1 3 40	39	7 1 1 11	13
	4 2 0 25	23	5 5 1 11	13
	4 2 2 21	21	6 4 0 12	13
	5 1 1 32	36	6 4 2 12	13
	3 3 3 29	36	5 5 3 11	10
	4 4 0 38	36	7 3 1 10	10
	5 3 1 27	27	8 0 0 6	16

J. Name., 2013, 00, 1-3 | 3

J. Name., 2013, 00, 1-3 | 5

J. Name., 2013, 00, 1-3 | 7

J. Name., 2013, 00, 1-3 | 9

Acknowledgements

Nicolas Goubet and Pierre-Antoine Albouy thank Pierre Fertey from the CRYSTAL beamline at SOLEIL and Alla Arakcheeva