Lionel Gélébart
email: lionel.gelebart@cea.fr

FFT-based simulations of heterogeneous conducting materials with combined non-uniform Neumann, Periodic and Dirichlet boundary conditions

Keywords: heterogeneous material, FFT, Periodic, Neumann, Dirichlet, thermal conductivity

Because of their simplicity, efficiency and ability for parallelism, FFT-based methods are very attractive in the context of numerical periodic homogenization, especially when compared to standard FE codes used in the same context. The purpose of the paper is to go beyond the use of periodic Boundary Conditions (BC), but keeping the advantages of FFT-based implementations. The present paper focuses on conductivity problems, considered as a first step towards mechanical problems. The proposed implementation is highly flexible, allowing to apply non-uniform loadings, to choose between periodic, Neumann and Dirichlet BC for each face of the unit-cell (of couple of faces for periodic BC), and to choose between different types of Finite Differences schemes (two types are considered here). The implementation relies on the use of Discrete Trigonometric Transforms (i.e. sine and cosine transforms) and their relation with the Discrete Fourier Transform on 4 times extended signals. The use of DTTs implicitly considers the signal to be anti-symmetric or symmetric with respect to each domain boundary. A direct relationship exists between the choice of a DTT (among 16 available) used in a given direction and the type of symmetry assumed on each boundary. Symmetry and antisymmetry assumptions are respectively related to Dirichlet BC and anti-symmetry. The implementation has been precisely validated for various types of loading (uniform or non-uniform, full or combined, Dirichlet, Neumann and periodic BC) from a direct comparison with the same FE simulations.

-Introduction

Solving mechanical problems defined on heterogeneous unit-cells, submitted to periodic BC, with FFTbased algorithm was proposed initially by Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. Its application to conducting materials is straightforward, and even simpler as the size of the tensors is decreased compared to mechanics (i.e vectors replaced by scalar, second order tensors by vectors etc…). Because of their simplicity, efficiency and ability for parallelism, they have become very popular and various improvements have been proposed since the original proposition. The reader is invited to read the review papers of Schneider [START_REF] Schneider | A review of nonlinear fft-based computational homogenization methods[END_REF], Lucarini [START_REF] Lucarini | Fft based approaches in micromechanics: fundamentals, methods and applications[END_REF] or Gierden [START_REF] Gierden | A review of fe-fft-based two-scale methods for computational modeling of microstructure evolution and macroscopic material behavior[END_REF], for an overview of the different advances.

An important limitation of current FFT-based implementations is their restriction to periodic BC. To circumvent that problem, the unit-cell can be embedded in a buffer whose properties can be adjusted to simulate Neumann BC (see for example [START_REF] Chen | Analysis of the damage initiation in a sic/sic composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF] for a simulation on a tube), or combined to a modified algorithm, to simulate Dirichlet BC [START_REF] Gélébart | A modified fft-based solver for the mechanical simulation of heterogeneous materials with dirichlet boundary conditions[END_REF]. Using a buffer zone is also proposed by To [START_REF] To | Fourier transform approach to nonperiodic boundary value problems in porous conductive media[END_REF] [START_REF] Quy | Fourier transform approach to homogenization of frequencydependent heat transfer in porous media[END_REF] with a rather different method, using form factors. Very recently, the question was addressed differently for conducting materials by Monchiet [START_REF] Monchiet | Fft based iterative schemes for composite conductors with uniform boundary conditions[END_REF]: instead of using a buffer, the unit-cell and the loading is duplicated and symmetrized consistently with the type of BC (Neumann or Dirichlet). The advantage is that a classical FFT-based implementation can be used 'as it is', the drawback is that the unit-cell is multiplied by 8 in 3D. In addition, their proposition is limited to uniform BC. Another approach was proposed recently by Grimm-Strele [START_REF] Grimm-Strele | Fft-based homogenization with mixed uniform boundary conditions[END_REF], for mechanical problems, in order to apply orthogonal mixed uniform BC proposed by Pahr [START_REF] Pahr | Influence of boundary conditions on computed apparent elastic properties of cancellous bone[END_REF]. This approach extends the previous works of Wiegmann [START_REF] Wiegmann | Fast poisson, fast helmholtz and fast linear elastostatic solvers on rectangular parallelepipeds[END_REF] to heterogeneous materials: in that case, fields are implicitly symmetrized with the use of discrete Sine and Cosine transforms, similarly to the use of the discrete Fourier Transform. However, the implementation was limited to uniform loadings and to a limited number of loading cases. Note that a similar approach was also recently followed by Wise [START_REF] Wise | Pseudospectral time-domain (pstd) methods for the wave equation: Realizing boundary conditions with discrete sine and cosine transforms[END_REF] for solving the wave propagation equation in a homogeneous medium with various types of BC and a Pseudo Spectral Time Derivative solver. Following the same approach for heterogeneous conducting materials, the aim of the paper is to propose a highly flexible implementation, allowing to apply non-uniform loadings, to choose between periodic, Neumann and Dirichlet BC for each face of the unit-cell (of couple of faces for periodic BC), and to choose between different types of Finite Differences schemes (two types are considered here).

The position of any point 𝑀 is given by a vector 𝒙, whose coordinates (𝑥 , 𝑥 , 𝑥 ) are related to an orthonormal coordinate system (𝑂, 𝒆 𝟏 , 𝒆 𝟐 , 𝒆 𝟑 ). The unit cell consists of a parallelepiped of dimensions 𝐿 , 𝐿 and 𝐿 . The unit-cell is bounded by 6 faces denoted 𝑆 and 𝑆 (𝑑 = 1: 3): 𝑆 with an outer normal vector 𝑵 𝒅 𝟎 = -𝒆 𝒅 , positioned at (𝑥 = 0), and 𝑆 positioned at 𝑥 = 𝐿 with a normal 𝑵 𝒅 𝟏 = -𝒆 𝒅 . The field equations governing the conduction problem are reported below: div(𝒒) = 0 𝒒 = -𝐾. 𝛁𝑻 𝑇 = 𝑇 * + 𝑇 [START_REF]Cast3m[END_REF] with 𝒒 the heat flux, 𝑇 the temperature and 𝐾 the heterogeneous thermal conductivity. Note that an arbitrary given field 𝑇 * is introduced, so that the unknown field is the temperature fluctuation 𝑇 .

Different boundary conditions (BC) can be associated to each face (or couple of faces). For periodic BC on faces 𝑆 and 𝑆 , 𝑇 and 𝒒. 𝒆 𝒅 have to be 𝐿 -periodic (same value on opposite points). For Dirichlet BC on a given surface 𝑆 , the fluctuation 𝑇 is null on 𝑆 (so that the value of 𝑇 corresponds to 𝑇 * ). For the implementation has to be flexible enough to account for 125 possible combinations (250 in fact as two discretization schemes will be considered). This flexibility comes from the fact that all these boundary conditions can be regarded as periodic boundary conditions on 2 or 4 times larger unit-cells using appropriate symmetries. At the centre of the proposed implementation, the FFT-based algorithm used with periodic BC is briefly summarized below.

Solving problem (1) with periodic boundary conditions, using a regular spatial discretization of the fields over the unit-cell, can be done with Discrete Fourier Transform (with the notation 𝐷𝐹𝑇(𝑓) = 𝑓 or simply 𝑓 ) as follows:

⎩ ⎪ ⎨ ⎪ ⎧ 𝒒 = -𝐾. (𝛁 𝐃 𝑻 * + 𝛁 𝐃 𝑻 ) 𝑝 = div 𝒒 -𝑘 𝛁 𝐃 𝑻 div 𝛁 𝐃 𝑻 = - 𝑝 𝑘 ⟺ ⎩ ⎪ ⎨ ⎪ ⎧ 𝒒 = -𝐾. 𝛁 𝐃 𝑻 * + 𝛁 𝐃 𝑻 𝑝 = div 𝒒 -𝑘 𝛁 𝐃 𝑻 T = 𝑝k ‖𝝃 𝑫 ‖ (≐ 𝐺 (𝑝)) (2)
In the present formulation (2), all the variables are discrete fields, div and 𝛁 𝐃 stand for discrete divergence and gradient operators associated to the choice of a finite discretization D. With this choice comes the definition of modified wave-vectors 𝝃 𝐃 , so that dıv (𝒒) = i𝒒 . 𝝃 𝐃 and 𝛁 𝐃 𝑻 = i𝑇 𝝃 𝐃 . Their expressions can be found in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] or [START_REF] Schneider | A review of nonlinear fft-based computational homogenization methods[END_REF] for different finite discretization schemes. Formulation (2) allows for the definition of a fix-point algorithm: first step, given 𝑇 evaluate 𝑝 in real space, second step, apply the discrete Green operator (i.e. evaluate 𝑇 from a back and forth in Fourier space), back to first step. The choice of the numerical parameter 𝑘 affects the convergence of the fix-point algorithm [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF], however, using an acceleration convergence procedure (the Anderson's acceleration as proposed in [START_REF] Chen | Analysis of the damage initiation in a sic/sic composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF] and implemented in [8]) drastically accelerates the convergence and reduces its sensitivity to 𝑘 .

After introducing the geometric discretization and the two Finite Differences (FD) schemes used in the present paper, the second step of the algorithm (apply Green operator T = 𝐺 (𝑝)) is first detailed, using Discrete Trigonometric Transforms to account for different types of symmetries in the evaluation of the discrete Green operator, 𝐺 . Then, the first step (evaluate 𝑝 in equation ( 2)) is described, using symmetries to extend the size of the unit-cell and evaluate discrete derivatives in real space. Finally, after a section devoted to additional remarks, comes the final validation with various different loadings demonstrating the flexibility of the implementation. Two Finite Differences (FD) schemes are used: the "Octahedral" scheme corresponds to the classical first order centered derivation scheme (points involved in the derivation are placed on an octahedron), with the notation 𝑂 replacing D in (2), and the "Hexahedral" scheme, strictly equivalent to the first order derivation obtained using hexahedral finite elements with reduced integration [START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF], with notation 𝐻 (points involved in the derivation are placed on an hexahedron). These different schemes are represented in Figure 1, on which important remarks are made:

-Unit-cell discretization and Finite Differences schemes

-For the 𝑂-scheme, all the quantities (𝑇 , 𝛁 𝐃 𝑻, 𝒒, 𝑝) are defined on the grid of nodes, whereas for the 𝐻-scheme, the quantities (𝑇 , 𝑝) and (𝛁 𝐃 𝑻, 𝒒) are defined on the grid of nodes and the grid of centres, respectively. -The evaluation of the heat flux 𝒒 = -𝐾. 𝛁 𝐃 𝑻 supposes 𝐾 and 𝛁 𝐃 𝑻 defined on the same grid.

So, as a consequence of the previous point: for the H-scheme, 𝐾 is evaluated on the grid of centres and for the 𝑂-scheme, 𝐾 is evaluated on the grid of nodes. Hence, the discretization of the material properties is slightly different depending on the FD scheme. -For the 𝑂-scheme, evaluation of derivatives at the boundary requires the knowledge of quantities outside of the domain. For the 𝐻-scheme, only the derivation of 𝛁 𝐃 𝑻 or 𝒒, defined at centers, requires quantities defined outside of the domain. In both cases, these quantities are defined according to the symmetry conditions associated to the choice of BC (see next sections). Using periodic boundary conditions in the 3 directions, with 𝑇 defined on the range [0: 𝑛 -1] : , the modified wave-vectors for the two schemes, 𝝃 𝑶 and 𝝃 𝑯 , can be found in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] or [START_REF] Schneider | A review of nonlinear fft-based computational homogenization methods[END_REF].

(0,0) (1,0) (0,1) (0,0) (1,0) (0,1)
𝛁 𝐇 𝑻 div H (𝒒)

(2,0)

Grid of nodes (for 𝑇 and 𝑝, using 𝑂and 𝐻-schemes, and 𝛁 𝐎 𝑻 and 𝒒, using 𝑂-scheme) Grid of centres (for 𝛁 𝑯 𝑻 and 𝒒, using 𝐻-scheme) Defining problem (1) on a 𝑀 times larger unit-cell and keeping the same grid discretization (i.e. with the range [0: 𝑀𝑛 -1] : ), the modified wave-vectors will be noted as 𝝃 𝑫𝑴 (in practice 𝑀 will be 2 or 4).

-Application of the Green Operator with symmetry conditions

Solving div 𝛁 𝐃 𝑻 =is equivalent to write 𝑇 = 𝐺 (𝑝), with 𝐺 a discrete Green operator (nonlocal). Evaluating the Green operator can be done straightforwardly with DFT if periodic BC are prescribed (see equation ( 2)). The aim of the present section is to build the discrete Green operator with non-periodic BC. It makes use of Discrete Trigonometric Transform (DTT), which can be related to standard DFT on larger unit-cells.

Note that below, DTTs gathers sine and cosine discrete transforms, and DTs, for discrete transforms, gathers the DTTs and the DFT.

Notations for discrete symmetry types and their associated DTTs

For a 1D discrete signal, the definitions of symmetry extensions outside of the domain depends: on the type of symmetry, which can be A (for Anti-symmetry) or S (for Symmetry), but also on the position of the (anti-)symmetry axis, whether located on an ending point, notation W (for Whole), or a halfelement after, notation H (for Half). Assuming that the location of the symmetry axis has the same location type T (whether W or H), the notation TXTY specifies the symmetry type with X and Y (whether A or S) respectively for the left (first) and right (or last) point of the discrete signal. In that case, the number of different combinations is of 8. As an example, the notation WAWS assumes an extension of the discrete signal with anti-symmetry (A) on the left side and symmetry (S) on the right side, the symmetry axis being located on the ending points (W). In addition, it must be noted that each type of symmetry extension corresponds to the definition of a given DTT (see [START_REF] Frigo | The design and implementation of fftw3[END_REF] and [START_REF] Wise | Pseudospectral time-domain (pstd) methods for the wave equation: Realizing boundary conditions with discrete sine and cosine transforms[END_REF]) and can also be associated to a type of boundary condition (Dirichlet or Neumann) as seen in equation [START_REF] Lucarini | Fft based approaches in micromechanics: fundamentals, methods and applications[END_REF]. Table 1 gathers all these informations. Finally, in case of 2D (or 3D) discrete fields, the symmetry is related to a given direction 𝑑 (𝑑 = 1: 3, in 3D) with the notation TX TY .

In addition, various graphical representations are given in appendix to visualize these symmetries and accommodate with the notation.

Discrete symmetry type (TXTY) Discrete Trigo. Transf. [ 1 : The table gathers the different discrete symmetry types (T=W or H for the location of symmetry axis, X=S or A for the type of symmetry on the left side, and Y=S or A for the right side), see explanations in the next for a more detailed description of the notation), the corresponding type of DT according to the names defined in the fftw library [START_REF] Frigo | The design and implementation of fftw3[END_REF], and the associated types of boundary conditions, see equation [START_REF] Lucarini | Fft based approaches in micromechanics: fundamentals, methods and applications[END_REF].

-DFT

The definition of the 1D DFT 𝐹 , of a sequence of 𝑛 points 𝑓, associated to the definition of the wavenumber 𝜉 , is:

𝐹 = 𝑓 exp (-i𝑘𝜉 ) , 𝜉 = 2𝜋𝑗 𝑛 , 𝑘 ∈ [0: 𝑛 -1] (3)
Extending this definition to 3D consists of applying 1D FFTs, separately in each direction, so that the definition, together with the definition of wave-vectors 𝝃 𝒋 , becomes:

𝐹 𝒌 = 𝑓 𝒋 exp (-i𝒌. 𝝃 ) 𝒋∈[ : ] : , 𝝃 𝒋 = 2𝜋𝑗 𝑛 , 2𝜋𝑗 𝑛 , 2𝜋𝑗 𝑛 , 𝒌 ∈ [0: 𝑛 -1] : (4)
As proposed in the previous section, we use the notation 𝐹 for the DFT of a signal with range [0: 𝑀𝑛 -1]. Note that the definition of the wave number 𝜉 , is deduced from the definition of 𝜉 as follows:

𝐹 , = 𝑓 exp (-i𝑘𝜉 , ) , 𝜉 , = 2𝜋𝑗 𝑀𝑛 , 𝑘 ∈ [0: 𝑀𝑛 -1] (5)

-DCT1 (WSWS)

The definition of the Discrete Cosine Transform of type 1, or DCT1, is associated to the extension of the signal defined on the range [0: 𝑛] to the range [0: 2𝑛 -1] that is assumed periodic (see appendix). Filling the extension is done according to symmetry assumptions on both sides of the initial domain. For DCT1, the left point (𝑖 = 0), belongs to the symmetry axis (W) and Symmetry (S) is assumed, this explain the first part of WSWS. The right point (𝑖 = 𝑛) has the same properties which explain the second part of WSWS. Note that W stands for "whole" when the symmetry axis crosses the "whole" point. H stands for "half" when the axis crosses between two points. With the symmetry WSWS, the extension is filled according to the relation:

𝑓 = 𝑓 for 𝑗 ∈ [0: 𝑛].
Applying the DFT definition (5) on the twice (𝑀 = 2) extended signal (with range [0: 2𝑛 -1] for 𝑗 and 𝑘) and using the relation 𝑓 = 𝑓 leads (after simple manipulations) to the definition of the DCT1 (as used in the library fftw [START_REF] Frigo | The design and implementation of fftw3[END_REF]):

𝐹 , = 𝑓 + 𝑓 (-1) + 2 𝑓 cos ( 𝑘𝜋𝑗 𝑁 -1 ) = 𝐹 , 𝑘 ∈ [0: 𝑛], (6)
with 𝑁 the number of points in the initial sequence (i.e. 𝑁 = 𝑛 + 1). DCT1 is performed on 𝑓 , for the complete range 𝑗 ∈ [0: 𝑛].

-DST1 (WAWA)

On the difference to the previous case, the twice extended signal is filled with anti-symmetry assumptions on left and right points (see appendix), leading to the relation: 𝑓 = -𝑓 for 𝑗 ∈ [0: 𝑛]. This condition necessitates that the signal vanishes on extreme points: 𝑓 = 𝑓 = 0. Following the same procedure as for DCT1 leads to the definition of DST1 (as used in fftw [START_REF] Frigo | The design and implementation of fftw3[END_REF]):

𝐹 , = -i. 2 𝑓 * sin 𝑘𝜋(𝑗 + 1) 𝑁 + 1 = -i𝐹 , 𝑘 ∈ [1: 𝑛 -1], (7)
with 𝑁 the number of non-vanishing points in the initial sequence (i.e. 𝑁 = 𝑛 -1). Actually, in that case, DST1 is performed on 𝑓 , for a reduced range 𝑗 ∈ [1: 𝑛 -1], excluding 𝑓 and 𝑓 necessarily set to 0. Equation ( 7) uses a renumbering with 𝑓 * = 𝑓 for 𝑗 ∈ [0: 𝑁 -1](= [0: 𝑛 -2]).

-DCT3 (WSWA)

In the previous cases, the symmetry conditions were the same on both sides. Now, the conditions are symmetry on left point and anti-symmetry on right point. Hence, the signal must be extended 4 times, from the range [0: 𝑛] to the range [0: 4𝑛 -1], to obtain a periodic signal consistent with the symmetry conditions (see appendix). For each quarter of the extended signal we use the relations deduced from the symmetry conditions: 𝑓 = -𝑓 , 𝑓 = -𝑓 and 𝑓 = 𝑓 , for 𝑗 ∈ [0: 𝑛]. Note that the first condition necessitates 𝑓 = 0. Following the same procedure as before, but with a 4 times extension (i.e. 𝑀 = 4), leads to the definition of DCT3 (as used in fftw [START_REF] Frigo | The design and implementation of fftw3[END_REF]):

𝐹 , = 2 𝑓 + 𝑓 cos 𝜋(𝑘 + 1/2)𝑗 𝑛 = 2𝐹 , 𝑘 ∈ [0: 𝑛 -1]. (8)
Here, DCT3 is performed on 𝑓 for a reduced range 𝑗 ∈ [0: 𝑛 -1], excluding 𝑓 necessarily set to 0.

-DST3 (WAWS)

DST3 is defined using the same procedure as DCT3, with the symmetry relations, 𝑓 = 𝑓 , 𝑓 = -𝑓 and 𝑓 = -𝑓 , for 𝑗 ∈ [0: 𝑛]. The last relation, together with the periodicity condition over the extended domain (𝑓 = 𝑓 ), necessitates 𝑓 = 0. Similar to DST3, the definition of DCT3 reads:

𝐹 , = -2i 𝑓 * (-1) + 2 𝑓 * sin 𝜋(𝑗 + 1) 𝑘 + 1 2 𝑛 = -2i𝐹 , 𝑘 ∈ [0: 𝑛 -1]. (9)
DST3 is performed on 𝑓 for a reduced range 𝑗 ∈ [1: 𝑛], excluding 𝑓 necessarily set to 0. Equation ( 9) uses a renumbering with 𝑓 * = 𝑓 for 𝑗 ∈ [0: 𝑛 -1].

-Relations between the symmetries of 𝑓 and its derivatives

For the sake of clarity, 2D fields are considered (adaptation to 3D or 1D fields being straightforward). Notation TX TY used below is detailed at the beginning of section 3. If a field 𝑓 has the symmetries TX TY /TX TY , for direction 1/direction 2, with X and Y the type of symmetry (A or S) and T the location of the symmetry (W or H, for field defined at nodes or at centers), then, the symmetry is inverted in the direction of derivation and kept in the other direction. In other words, 𝑑𝑓/𝑑𝑥 has the symmetry TX TY /TX TY and 𝑑𝑓/𝑑𝑥 the symmetry TX TY /TX TY , with the notation A = S and S = A.

For example, assuming a field 𝑓 with the symmetry WSWA/WAWA (i.e. a field defined at nodes (W), Symmetric (S) and Anti-symmetric (A) on the left and right sides in direction 1, and Anti-symmetric on both sides in direction 2), the derivative 𝑑𝑓/𝑑𝑥 has the symmetry WAWS/WAWA and 𝑑𝑓/𝑑𝑥 has the symmetry WSWA/WSWS.

Note that above, the location T is not modified by the derivative, which corresponds to a Finite Difference scheme for which the quantity and its derivative is evaluated at the same physical point (the nodes here), such as the 𝑂-scheme used in the present paper. Alternatively, the derivative can be evaluated at centers, between the nodes where the quantities are given (see Figure 1), or conversely at nodes if quantities are given at centers. This is the case of the 𝐻-scheme, for which the derivation 'inverts' location T in T (with the notation W = H and H = W).

A general rule, useful for a flexible implementation, is proposed below. If a 3D field, 𝑓, has the symmetries TX TY /TX TY /TX TY , its derivative in direction d , 𝑑𝑓/𝑑𝑥 , exhibits :

-inverse symmetry types in direction 𝑑 ( X and Y ) and no inversion in the other directions, -for 𝐻-scheme, an inverse symmetry location 𝑇 (no inversion for the 𝑂-scheme).

-Relations between symmetries and non-periodic BC

As mentioned in the previous sub-sections, the WA symmetry condition is associated to a null value at the boundary: 𝑇 = 𝑇 = 0 for WAWA, 𝑇 = 0 for WAWS, and 𝑇 = 0 for WSWA. Hence, using the DTT associated to the anti-symmetry conditions, supposes implicitly a null Dirichlet BC that is 𝑇 = 0.

According to the previous subsection, if a 1D field 𝑇 has a WS symmetry condition then 𝑔 = 𝑑𝑇 /𝑑𝑥 has a WA or HA symmetry (whether the 𝑂or 𝐻-scheme is used). As noticed just before, the WA symmetry is associated to a null value at the boundary. For the HA symmetry, the relation is (for the left point) 𝑔 = -𝑔 . However, located at centers 𝑔 and 𝑔 are defined on both side of the boundary so that the value at the middle point, that falls on the boundary, is evaluated by 𝑔 / = (𝑔 + 𝑔 )/2 = 0. For 3D fields, applying a WS symmetry in direction 𝑑 is equivalent to a null Neumann condition that is 𝛁 𝑫 𝑻 . 𝑵 𝒅 = 0 (with 𝑵 𝒅 the outer normal to the boundary in direction 𝑑).

In brief, applying null Neumann or Dirichlet BC on 𝑇 can be done assuming appropriate symmetries for 𝑇 as follows: The connections established in the previous sub-sections between various DTs and DFT can be summarized as follows:

S
𝑆 (𝑝̂ ) = 𝛼𝑆 (𝑝) (11)
𝑆 is a selection operator for arrays in real space. The range of 𝑝 is [0: 𝑛] and the range of 𝑆 (𝑝) is [0: 𝑁 -1], with 𝑁 the number of selected points. The points selected by 𝑆 for each DT (see previous sub-sections) are gathered in Table 2.

𝑆 is a selection operator for arrays in Fourier space and 4 times extended unit-cells. The range of 𝑝̂ is [0: 4𝑛 -1] and 𝑆 (𝑝̂ ) is an array of range [0: 𝑁 -1], with 𝑁 the number of selected points. The points selected by 𝑆 for each DT are gathered in Table 2. Note that in the previous subsections, selection is made for different extension factor : 1 (trivially, for DFT), 2 (for DCT1 and DST1) or 4 (for DCT3 and DST3). To simplify, selection is performed for the same extension factor 𝑀 = 4 (i.e. in the same range [0: 4𝑛 -1]), assuming 4 periodic repetitions for DFT, and 2 for DCT1 and DST1). Then the selection is done by multiplying by 2 (for DCT1 and DST1) or 4 (for DFT) the ranges selected for 𝑀 = 2 or 𝑀 = 1, as done in Table 2.

The factor 𝛼 is not reported in Table 2 as its effect will be annihilated by a factor 1/𝛼 in a further step (equation ( 13)). It can be easily deduced from previous sub-sections.

Once a selected range of Fourier coefficients 𝑆 (𝑝̂ ) has been evaluated from DT applications (equation ( 11)), the Green operator 𝐺 is applied for the same selection range:

𝑆 𝑇 = 𝑆 𝐺 𝑝̂ = 𝑆 𝐺 𝑆 (𝑝̂ ) = 𝑆 (𝑝̂ ) 𝑘 ‖𝑆 (𝜉 )‖ (12)
where 𝜉 are the modified wave-numbers evaluated for a 4 times extended domain.

Equation ( 11) is now applied to deduce 𝑆 𝑇 from 𝑆 𝑇 according to:

𝑆 𝑇 = 𝑆 𝑇 𝛼 (13)
Then, inverse DT can be applied to 𝑆 𝑇 to obtain 𝑆 𝑇 .

Finally, the points which have not been selected by 𝑆 in [0: 𝑛], that corresponds to points located at the boundaries, must be filled: in case of Dirichlet BC on 0 and/or 𝑛, then 𝑇 = 0 and/or 𝑇 = 0, in case of periodic BC, then 𝑇 = 𝑇 .

For the sake of simplicity, the presentation was done in 1D. However, the extension to 3D is rather straightforward and relies on the fact that the 3D transform corresponds to a succession of 1D transforms in each direction (with appropriate DT corresponding to the BC in each direction). The selection operators 𝑆 and 𝑆 now operates on 3D arrays with the selection rules proposed in Table 2 for each direction. In equation ( 12), the wave-number 𝜉 in 1D is replaced by the wave-vector 𝝃 𝑫𝟒 , whose definition is given in [19] [15].

In terms of implementation, it must be emphasized that if 4 times extended fields are considered in Fourier space, only the values selected by 𝑆 appear in equations (11)(12)(13), and each quantity of type 𝑆 (𝑋 ) corresponds to a 3D array of range [0: 𝑁 -1] : , with 𝑁 the number of points used in DT (see Table 2). Note that the array 𝑆 (𝝃 𝑫𝟒 ) in equation ( 12) corresponds to a selection of wavevectors that can be evaluated once, at the beginning of the iterative fix-point algorithm.

To summarize, assuming that 𝑝 and 𝑇 have the same symmetries, the application of the (non local) Green operator 𝑇 = 𝐺 (𝑝) is rather simple and flexible as soon as the selection operators 𝑆 and 𝑆 are correctly initialized in agreement with Table 2. The four steps are:

-DT on a selected range (𝑆 ) of 𝑝 (→ equation ( 11)), -Application of the Green Operator, for selected (𝑆 ) wave-vectors (→ equation ( 12)), -Inverse DT (of equation ( 13)), to obtain 𝑇 on the selected range (𝑆 ), -Adjust non-selected points (for Dirichlet and Periodic BC). 2 : Range selections in real (𝑆 ) and Fourier space (𝑆 ) for various Discrete Transforms

-Evaluation of 𝑝 (see equation (2)) with non-uniform BC

The fix-point algorithm described in equation ( 2) (left), summarizes in two steps: evaluate 𝑝(𝑇 ), apply the Green operator (𝑇 = 𝐺 (𝑝)). The second step was the purpose of the previous section. The evaluation of 𝑝, described below, is done in real space using finite differences with the 𝑂or 𝐻-scheme. Recall that 𝑝 is defined on the grid of nodes, as well as 𝑇 , with the range [0: 𝑛].

-Evaluation with non-uniform Dirichlet or periodic BC and null Neumann BC

The definition of 𝑝 is given in equation ( 2) and it is divided in two parts 𝑝 = 𝑝 -𝑝 with 𝑝 = div (𝒒) and 𝑝 = 𝑘 div 𝛁 𝑫 𝑻 where the norm of 𝑝 is used as convergence criterion. The important steps of the evaluation are reported below.

Step 1: Evaluating 𝛁 𝑫 𝑻 The derivation of 𝑇 using 𝐻-scheme, naturally gives the complete field of 𝛁 𝑫 𝑻 defined at centres (range [0: 𝑛 -1] : ). On the contrary, using the 𝑂-scheme, 𝛁 𝑫 𝑻 being defined at nodes (range [0: 𝑛 ] : ), it requires the knowledge of 𝑇 in the 6 planes (indices 𝑖 = -1 and 𝑖 = 𝑛 + 1, for each direction 𝑑) surrounding the unit-cell. The value of 𝑇 for these additional points are evaluated according to the symmetry conditions, chosen in agreement with the choice of BC (see section 3.7).

Step 2: Evaluating div 𝛁 𝑫 𝑻 div 𝛁 𝑫 𝑻 is defined at nodes (range [0: 𝑛 ] : ). Using the 𝑂-scheme, 𝛁 𝑫 𝑻 being defined at nodes (with the same range [0: 𝑛 ] : ), it requires the knowledge of 𝛁 𝑫 𝑻 in the 6 surrounding planes (of nodes). Similarly, using the 𝐻-scheme requires the knowledge of 𝛁 𝑫 𝑻 in the 6 surrounding planes (of centres). In both cases, the symmetry conditions of 𝛁 𝑫 𝑻 , deduced from the symmetry of 𝑇 in section 3.6, are used.

Step 3: Evaluating 𝒒 = -𝐾. (𝛁 𝐃 𝑻 * + 𝛁 𝐃 𝑻 ).

The evaluation of 𝒒 is very simple, and emphasis is put on the term 𝛁 𝐃 𝑻 * . Defined during the initialization, it allows for non-uniform Dirichlet or periodic BC.

For Dirichlet BC on both sides (i.e. WAWA symmetry), in a given direction 𝑑, 2D heterogeneous temperatures fields are imposed on the two opposite faces of the unit-cell, for example 𝑇 * [0,0: 𝑛 , 0: 𝑛 ] and 𝑇 * [𝑛 , 0: 𝑛 , 0: 𝑛 ] if 𝑑 = 1. Any field compatible with these boundary conditions can be used for 𝑇 * , a simple linear interpolation is chosen. The same linear interpolation is used if Dirichlet BC are applied on 1, 2 or 3 couples of opposite faces.

For Dirichlet BC on one side and Neumann BC on the other side (i.e. WAWS or WSWA symmetry): the same linear interpolation is used, with the value of 𝑇 * on the Neumann side arbitrarily chosen equal to 𝑇 * on the Dirichlet side.

Finally, note that for the 𝐻-scheme, the definition of 𝑇 * on the range [0: 𝑛 , 0: 𝑛 , 0: 𝑛 ] is sufficient to evaluate 𝛁 𝐃 𝑻 * on the grid of centres. For the 𝑂-scheme, 𝑇 * must be defined on the range [-1: 𝑛 + 1, -1: 𝑛 + 1, -1: 𝑛 + 1] to evaluate 𝛁 𝐃 𝑻 * on the grid of nodes. The linear function used for the interpolation inside the unit-cell is kept for extrapolation at points located outside of the unit-cell.

For periodic BC, gradient fields 𝛁 𝐃 𝑻 * as proposed for 2-sides Dirichlet BC, can be prescribed. As the temperature fluctuation 𝑇 is periodic, its 1D-average gradient in a given direction 𝑑 is a null 2D field, so that the 1D-average gradient of the temperature 𝑇, in direction 𝑑, is equal to the 1D-average gradient of 𝑻 * ,in direction 𝑑.

Note that the specific case of uniform 𝛁 𝐃 𝑻 * , is extensively used in the context of first order homogenization, with both Dirichlet and periodic BC.

Step 4: Evaluation of div (𝒒)

The same considerations as those described for step 2 (evaluation of div 𝛁 𝑫 𝑻 ) apply to step 4. However, it is worth noting that the symmetry conditions used for 𝛁 𝑫 𝑻 are also used for 𝒒. Hence, a WS symmetry condition on 𝑇 , corresponding to a null-Neumann BC on 𝛁 𝑫 𝑻 (i.e. 𝛁 𝑫 𝑻 . 𝒏 𝒅 = 0), also corresponds to a null-Neumann BC on 𝒒 (i.e. 𝒒. 𝒏 𝒅 = 0).

-Adding non-zero Neumann BC

As demonstrated above, using the DT framework and the set of equation ( 2) only allows to apply null-Neumann BC. The only way to introduce non-zero Neumann BC is to define an appropriate volume source field 𝑟 * null everywhere except at, or in the inner neighbourhood of, the boundary. In practice, 𝑟 * is defined by 𝑟 * = div (𝒒 * ) with a flux 𝒒 * , null everywhere in the unit-cell except at, or in the outer neighbourhood of, the boundary (precisions will be given below for the definition of 𝒒 * , consistently with the applied fluxes). The definition of 𝑝 is then modified according to:

𝒒 = -𝐾. 𝛁 𝐃 𝑻 * + 𝛁 𝐃 𝑻 + 𝒒 * 𝑝 = div 𝒒 -𝑘 𝛁 𝐃 𝑻 (14)
As a consequence, as soon as 𝒒 * is correctly initialized, the minor modification consists in adding 𝒒 * to -𝐾. 𝛁 𝐃 𝑻 * + 𝛁 𝐃 𝑻 step 3. Details about the initialization of 𝒒 * are given below, for a Neumann BC applied in direction 1, on the left side, so that the outer normal vector 𝑵 is -𝒆 𝟏 and (𝒒. 𝒆 𝟏 ) = (𝑞 ) = -𝑄 * , with 𝑄 * an applied outcoming heat flux (a 2D field in the most general case).

For the 𝑂-scheme, 𝒒 being defined at nodes, the condition (𝑞 ) = -𝑄 * is prescribed at nodes and the range of the applied flux 𝑄 * , is [0: 𝑛 , 0: 𝑛 ]. As div is applied to 𝒒 * , it must also be defined on planes of nodes surrounding the unit-cell, so the range of 𝒒 * is [-1: 𝑛 ] : . 𝒒 * is null everywhere except, for the BC considered here, on the two first planes of nodes for which :

𝑞 * [-1,0: 𝑛 , 0: 𝑛 ] = -2𝑄 * [0: 𝑛 , 0: 𝑛 ] (first plane out of the unit -cell) 𝑞 * [0,0: 𝑛 , 0: 𝑛 ] = -𝑄 * [0: 𝑛 , 0: 𝑛 ] (boundary plane)

Assuming, as in 4.1 -step 4, that 𝒒 = -𝐾. 𝛁 𝐃 𝑻 * + 𝛁 𝐃 𝑻 has the same symmetries as 𝛁 𝐃 𝑻 , then a null-Neumann BC applies to 𝒒 in direction 1 so that (𝑞 ) = 0. Hence, the choice (𝑞 * ) = -𝑄 * made for the boundary plane (equation ( 15 As div is applied to 𝒒 * , it must also be defined on planes of centres surrounding the unit-cell, so the range of 𝒒 * is [-1: 𝑛 ] : . 𝒒 * is null everywhere except, for the BC considered here, on the first plane of centres (outside the unit-cell) for which :

𝑞 * [-1,0: 𝑛 -1,0: 𝑛 -1] = -2𝑄 * [0: 𝑛 -1,0: 𝑛 -1] (16)
With this choice, the flux at the boundary, located in the middle of the two centres consists of (𝑞 ) = Here, attention must be paid to the points located at the edges of the domain [-1, -1: 𝑛 , -1: 𝑛 ], which are not considered in equation [START_REF] Quy | Fourier transform approach to homogenization of frequencydependent heat transfer in porous media[END_REF]. They must be filled with the symmetry conditions associated to 𝛁 𝐃 𝑻 in direction 2 and 3. This remark holds for the 𝐻-scheme. For the 𝑂-scheme, the points at the edges are not used to evaluate div (𝒒 * ) so that applying these symmetry conditions is not mandatory. The proposed implementation is then able to simulate any of the configurations given by these inputs. Considering the choice of BC, it must be emphasized that 5 × 5 × 5 = 125 types of BC can be prescribed. This proof of concept relies on the matlab language, and the toolbox matlab-dtts proposed by Treeby [START_REF] Treeby | Matlab discrete trigonometric transform library[END_REF], who implemented an interface with the 8 Discrete Sine and Cosine Transforms available in the library fftw [START_REF] Frigo | The design and implementation of fftw3[END_REF].

The next step will consist of implementing this method in the massively parallel AMITEX code [8].

-Convergence criterion

The convergence criterion, associated to the equilibrium equation div (𝒒) = 0, reads:

𝜀 = ‖div (𝒒)‖ 𝒒 𝑑𝑥 < 𝜀 (17)
with 𝑑𝑥 the voxel size. The iterative algorithm is interrupted as soon as this criterion is satisfied.

-Convergence Acceleration

The convergence acceleration technique proposed by Anderson [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF], implemented in the code AMITEX [8] to accelerate the FFT-based fix-point algorithm [START_REF] Chen | Analysis of the damage initiation in a sic/sic composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF], is also implemented here with a depth of 3 and an acceleration every 3 iterations. Briefly, it consists of saving in a buffer the last 4 couples of solution and residual fields, and to propose every 3 iterations, a new solution field deduced from the saved fields (see [START_REF] Chen | Analysis of the damage initiation in a sic/sic composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography[END_REF], for a more detailed description).

-Discussion

In the present paper, choice is made to evaluate 𝑝 = div 𝒒 -𝑘 𝛁 𝐃 𝑻 in real space and then perform a DT, apply the Green operator in Fourier space and inverse DT to obtain 𝑇 . Another, and more usual way, is to evaluate 𝒒, perform the DTs, apply the Green operator (different from the one presented here) and inverse DTs to obtain 𝛁 𝐃 𝑻 .

In the first implementation, DT/inverse DT are performed on a scalar field, for the second, they are performed on vector fields and the number of DT/inverse DT is multiplied by 3.

The arguments for reducing the number of DT/iDT are the following:

-The 𝐻and 𝑂-schemes proposed here are very 'local' (8 or 6 points are involved in the derivation operators) and it can be more efficient to apply derivatives in real space, instead of applying additional transforms and inverse transforms. -This point is reinforced when using Discrete Sine and Cosine Transforms, whose implementation in the library fftw is, according to the authors, less efficient than DFT. -It is also reinforced when dealing with a massively parallel implementation (the next target after this proof of concept). Actually, with data distributed over a 2D pencil decomposition as in [8], the 3D transform requires a transposition of data to perform the 1D transforms in each direction. This transposition corresponds to a "ALLTOALL" communication between all the processes that comes with an additional communication cost.

-Validation

As simulations performed with the 𝐻-scheme are equivalent to Finite Element simulations with linear hexaedral finite elements (Hex8R) [START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF], a quantitative cross validation can be made by comparisons of strictly equivalent simulations performed with the FE and FFT-based codes. For that purpose a cubic unit-cell (dimension 1 × 1 × 1) is used, with a non-centred spherical inclusion (radius 0.3, centre at (0.15,0.4,0.6), see Figure 2). The unit-cell is simple but non-symmetric, with the inclusion crossing the boundary. This choice allows to demonstrate that the method is not limited to unit-cells with inner symmetries and that it is capable to account for material discontinuities located at the boundary. To enforce strong heterogeneities a high contrast is chosen for the thermal conductivity (matrix 1, inclusion 10 ). The unit-cell is discretized with 32x32x32 voxels. A high accuracy is prescribed for this validation (convergence criterion 𝜀 = 10 ).

Note that all the simulations below have been performed with a single flexible implementation that allows to choose between 125 types of symmetries and 2 finite difference schemes. A selection of loadings of increasing complexity is proposed below. Table 3 gathers the different loadings used in the following sub-sections. Table 3 : Description of the loading cases used for the validation of the method. The unit-cell and the faces 𝑆 are defined in Figure 2. For the "periodic" case, 𝑇 periodic (equal on opposite faces). For the "non-uniform" case, 𝑇 and 𝑄 * are defined in section 6.3. For the other cases, applied 𝑮 * or 𝑸 * (see definition in section 6.1) are set to one, leading to the BC reported in the present table.

Neumann

  BC on a given surface 𝑆 , the normal heat flux 𝒒. 𝑵 𝒅 𝒋 must be equal to an arbitrary applied flux 𝑄 * . The standard implementation of FFT-based methods for the simulation of heterogeneous materials only considers periodic BC on the 3 couples of faces. The proposed implementation allows for every combinations of BC on the 6 different faces. With 5 combinations per couple of opposite faces (Periodic, Dirichlet/Dirichlet, Dirichlet/Neumann, Neumann/Dirichlet and Neumann/Neumann)

Figure 1 ,

 1 Figure1, illustrates the definitions used for the discretization of the unit-cell. The grid of nodes touches the boundaries of the unit-cell, whereas the grid of centres, translated by half a voxel do not touch the boundary. The size of the grid of centres is 𝑛 𝑛 𝑛 with the range [0: 𝑛 -1] : , and the size of the

Figure 1 :

 1 Figure 1 : Illustration of the unit-cell discretization: definition of the ranges of the grids of corners and centres, locations of the different quantities depending on the Finite Difference scheme (𝑂 or 𝐻), and illustration of the 𝑂 and 𝐻 schemes.

4 •

 4 𝛁 𝐎 𝑻 , div (𝒒) Definitions for this 2D grid: • Dimensions are 𝑛 1 = 𝑛 2 = Range of the grid of nodes : [0: 𝑛 1 , 0: 𝑛 2 ] • Range of the grid of centres : [0: 𝑛 1 -1,0: 𝑛 2 -1]

10 ) 3 . 8 -

 1038 symmetry for 𝑇 ⟹ null Neumann BC (𝛁 𝑫 𝑻 . 𝑵 𝒅 = 𝟎) A symmetry for 𝑇 ⟹ null Dirichlet BC (𝑻 = 𝟎) (Green operator with null-Dirichlet, null-Neumann and Periodic BC This section describes step by step the application of the Green operator 𝑇 = 𝐺 (𝑝), assuming that 𝑝 and 𝑇 are 1D fields with the same range [0: 𝑛], and the same symmetries (i.e. BC).

  )) allows satisfying the expected Neumann BC (𝑞 ) = (𝑞 ) + (𝑞 * ) = -𝑄 * . The choice 𝑞 * = -2𝑄 * made for the first plane out of the unit-cell, is done to respect the condition that the derivative evaluated at the boundary with the 𝑂-scheme, * = * [ , : , : ] * [ , : , : ] = * [ , : , : ] , should be equal to the derivative evaluated from nodes located inside the cell, * = * [ , : , : ] * [ , : , : ] = * [ , : , : ] = * . For the 𝐻-scheme, 𝒒 being defined at centres, the condition (𝑞 ) = -𝑄 * is prescribed at the middle points of two centres separated by the unit-cell boundary. Hence, the range of the applied flux 𝑄 * is [0: 𝑛 -1,0: 𝑛 -1].

(

  𝑞 ) + (𝑞 * ) , with (𝑞 * ) = * [ * and (𝑞 ) = 0 (as a null-Neumann BC is imposed to 𝒒 , see section 4.1 -step 4). Hence, the expected Neumann BC, (𝑞 ) = -𝑄 * , is satisfied.

  The definition of 𝒒 * proposed here for a non-zero Neumann BC on face 𝑆 can be easily extended to any other face. Submitted to a non-zero Neumann BC5 -Additional remarks 5.1 -A proof of conceptThe implementation made of the FFT-based method used for the validation tests, is a proof of concept demonstrating its flexibility. The user defined inputs are:-the discrete thermal conductivity field defined on the whole unit-cell, -for each couple of opposite faces, the choice of BC (Dirichlet/Dirichlet, Dirichlet/Neumann, Neumann/Dirichlet, Neumann/Neumann or Periodic BC corresponding to the 5 possible symmetry conditions WAWA, WAWS, WSWA, WSWS and P), -the fields 𝑇 * and 𝒒 * associated respectively to the Dirichlet and Neumann BC, -the FD scheme chosen between the 𝑂and the 𝐻-scheme.

Figure 2 : 1 𝑞

 21 Figure 2 : Thermal conductivity field corresponding to a resistive spherical inclusion of radius 0.3, centred at (0.15, 0.4, 0.6) in a 1x1x1 unit-cell. The six different faces 𝑆 , used to apply various loadings with either periodic, Dirichlet or Neumann BC are reported.

6. 1 -

 1 Uniform loadings with full Periodic, Dirichlet and Neumann BC. The loadings considered here are given by, for all 𝑥 at the boundary of the unit-cell, 𝑇(𝒙) = 𝑮 * . 𝒙 + 𝑇 (𝒙), 𝑇 submitted to null-Dirichlet or periodic BC, or 𝒒(𝒙). 𝑵(𝒙) = 𝑸 * . 𝑵(𝒙) for Neumann BC. They are called 'uniform' in the sense that 𝑮 * and 𝑸 * are uniforms over the unit-cell. Such loadings are commonly used in the homogenization context as they provide respectively, an upper bound, an of discrete signals and corresponding Discrete Trigonometric Transforms : graphical representations. The figures below give different examples of symmetries for discrete signals, symmetries that are implicitly assumed when using the corresponding Discrete Trigonometric Transforms.

Figure 5 :

 5 Figure 5 : Discrete signal with symmetry HAHA (upper) or HSHS (lower). The symmetrized signal has a 2n period. HAHA corresponds to DST2 and HSHS to DCT2 (according to the fftw definition used in the present paper).

Figure 6 :

 6 Figure 6 : Discrete signal with symmetry WAWA (upper) or WSWS (lower). The symmetrized signal has a 2n period. WAWA corresponds to DST1 and WSWS to DCT1 (according to the fftw definition used in the present paper).

Figure 7 :

 7 Figure 7 : Discrete signal with symmetry WAWS. The symmetrized signal has a 4n period. WAWS corresponds to DCT3 (according to the fftw definition used in the present paper).

Acknowledgments

The author would like to thank Marc Josien from CEA-Cadarache for fruitful discussions and Olivier Fandeur from CEA-Saclay who introduced reduced integration elements for thermal problems in the FE code CAST3M, required for the cross validation between FE and FFT codes.

estimate and a lower bound, for the homogenized conductivity. It is easily demonstrated that 𝛁𝑻 = 𝑮 * for periodic and Dirichlet BC and 𝒒 = 𝑸 * for Neumann BC. In the simulations below the loadings are applied with 𝐺 * = 𝐺 * = 𝐺 * = 𝑄 * = 𝑄 * = 𝑄 * = 1.

Finite element simulations are performed with the FE code CAST3M [START_REF]Cast3m[END_REF], using Hex8R elements. Table 4 gathers the average heat flux and temperature gradient for the three FE simulations and the relative error (definition in the table legend) between FE and FFT-based simulation. As expected, the relative error is very low with a maximum value of less than 10 for Neumann simulation. This is consistent with the low value of the criterion of 10 used in the FFT-based algorithm. The precision on the average temperature gradient is null for Dirichlet and periodic BC: actually this quantity is strictly imposed by the algorithm. On the contrary, for Neumann BC, the applied flux is not strictly applied by the algorithm and the precision of the average flux depends on the convergence criterion. At last, the maximum relative error is by far the error on ∇𝑇 for the Neumann BC. It is related to the fact that, on face 𝑆 , a constant heat flux is applied directly to the inclusion whose conductivity is very low. That induces very large (and heterogeneous) temperature gradient as observed in the table (when comparing ∇𝑇 to the other components). The higher number of iterations also reveals that the Neumann BC simulation (with the low conductive inclusion crossing the boundary) is more delicate.

It.

Factor ∇𝑇 ∇𝑇 ∇𝑇 q q q Periodic 47 1.000 1.000 . On the opposite face 𝑆 , the flux is null except on a square of size 0.5 × 0.5, centred on the face. The comparison with FE simulations, on the average gradient and heat flux in Table 4, reveals errors lower than 10 . As expected from equation ( 16) and the chosen BC, the average flux in direction 1, 𝑞 is equal to 𝑄 * /4 (𝑄 * being applied on a surface of size 0.5 × 0.5). The non-vanishing values for 𝑞 and 𝑞 is explained once again by the coupling arising when using Dirichlet BC (see equation( 16) ).

For this non-uniform case, a relative error on the complete flux fields is also evaluated with 𝐸 = ‖𝑞 -𝑞 ‖ /‖𝑞 ‖ . The values, gathered in Table 5 for different convergence criteria, confirms, now on full fields, that the error is very low and reduces with the convergence criteria. For the 𝑂-scheme, without any appropriate FD code that could provide a rigorous comparison with FFT-based results, validation is then reduced to a comparison between FFT-based simulations performed with the 𝑂and 𝐻-schemes. Average heat flux and temperature gradient are compared quantitatively as a function of the spatial resolution in Table 6. As expected, the difference decreases when increasing the spatial resolution, both simulations converge towards the solution. The largest difference is observed with Neumann BC, especially when applied to the face crossing the low conducting inclusion. For all the simulations, it has been observed that the 𝐻-scheme converge faster than the 𝑂-scheme, and the evolution of the differences in Table 6 is mostly due to the 𝑂-scheme. As an example, for the full-Neumann BC simulation, ∇𝑇 evolves between -20.72 and -20.61 when increasing the resolution from 32 to 128, for the 𝐻-scheme, and between -16.30 and -19.56 for the 𝑂scheme.

To complete qualitatively this validation, the heat flux field 𝑞 evaluated for the non-uniform loading with the 𝑂-scheme, is displayed in Figure 4. It compares very well with Figure 3 The first point concerns the use of discrete sine and cosine transforms, whose definition can be linked to the definition of discrete Fourier transform applied to a 4 times extended signal, using appropriate symmetries. Hence, the periodic Green operator can be applied to this extended signal for an appropriate selection of wave vectors. Besides, the discrete sine or cosine transforms have also to be performed on an appropriate selection of grid points. With these two selection operators, the discrete Green operator, associated to any set of symmetry conditions, can be easily applied.

The second point is the term 𝑝 (see equation ( 2)), evaluated in real space using discrete derivation schemes and taking into account the appropriate symmetries, it also introduces the non-uniform loading through the quantities 𝑇 * and 𝒒 * whose definition is associated to Dirichlet (and periodic) and Neumann BC, respectively.

This flexible implementation has been validated from a cross-comparison between FFT-based simulations, using the 𝐻-scheme for discrete derivation, and a FE simulations, using linear hexahedral elements with reduced integration. As expected from these strictly equivalent simulations [START_REF] Schneider | Fft-based homogenization for microstructures discretized by linear hexahedral elements[END_REF], the simulations provide the same results, up to the precision of the FFT-based iterative algorithm.

𝑞 𝑞

The implementation also proposes two different types of Finite Difference schemes: the 𝐻-scheme, for which the derivation operator changes the support from a grid of nodes to a grid of centres, or viceversa, and the 𝑂-scheme, for which all quantities are defined on the same grid (of nodes). It should pave the way to implement other FD schemes.

From the simulations performed with a low conductivity inclusion (contrast 1000 with the matrix) crossed by the boundary, and a very high precision (𝜀 = 10 ), various remarks can be made: the number of iterations is very reasonable (between 46 and 152) and it is higher when Neumann BC are applied to the face crossed by the inclusion (in agreement with intuition). The present implementation appears much more efficient than using a buffer zone to apply Dirichlet BC [START_REF] Gélébart | A modified fft-based solver for the mechanical simulation of heterogeneous materials with dirichlet boundary conditions[END_REF]. Limited to unit-cells with parallelepiped shapes, the present method is not able to deal with volumes with arbitrary shapes and in that case using a buffer zone [START_REF] Gélébart | A modified fft-based solver for the mechanical simulation of heterogeneous materials with dirichlet boundary conditions[END_REF] [17] [START_REF] Quy | Fourier transform approach to homogenization of frequencydependent heat transfer in porous media[END_REF] remains an interesting solution. It also extends the capabilities of the method proposed recently [START_REF] Monchiet | Fft based iterative schemes for composite conductors with uniform boundary conditions[END_REF] to non-uniform loadings as well as Dirichlet/Neumann BC on opposite faces.

As future prospects, the present implementation proposed in the context of conducting materials is a first step towards mechanical problems. Besides, the present proof of concept implemented through the matlab language, will have to be introduced in the FFT-based simulation code AMITEX [8] to take benefit of its massively parallel implementation.


Antonio Imbruglia
email: antonio.imbruglia@st.com

Daniela Cancila
email: daniela.cancila@cea.fr

Marina Settembre
email: msettembre@fub.it

A da Us er J o ur na l 5G Communication and Security in Connected Vehicles

Keywords: Antonio Imbruglia, Daniela Cancila, ECS-SRIA, 5G, Cybersecurity, Safety, V2X

5G communication and security in connected vehicles

Introduction

The exponential development of IoT systems, electric embedded devices and cyber-physical systems (CPS), with stronger intelligence and interconnection requirements, takes advantage of the development of 5G networks. In spite of new emerging application scenarios, new thorny issues in cyber-security should be considered due to an increased attack surface and evolving threat landscape. These systems are going to operate in a more dynamic and open environment respect to traditional systems, with increasing needs to exchange information with other systems. Consider electric and autonomous vehicles, where the connection with both other vehicles and the environment (e.g. infrastructure, roads, traffic lights pedestrians) can improve road safety, reduce environmental impact, and provide a new driving experience. However, if susceptible to cyber-attacks, these systems could have catastrophic consequence. Moreover, in traditional systems, specialized operators manage the system. For example, in the railway application domain, a team of experts manages the system (including train maintenance). In this new scenario, the actors involved in the management can be more heterogenous, belonging to different organizations and not always so specialized to recognize the nominal or degraded behavior of the system, following for example a cyberattack.

The paper, without claiming to be exhaustive, offers insights and reflections on the very broad topic that integrates innovative devices, 5G and cybersecurity by illustrating the main European directions indicated by both the Strategic Research and Innovation Agenda for Electronic Component and Systems (ECS-SRIA), [START_REF]Electronic Components and Systems -Strategic Research and Innovation Agenda[END_REF] and some insights into the evolution of 5G network standards. For the sake of concreteness, some considerations will be reanalyzed in Vehicles to X (V2X) case studies.

ECS-SRIA

The Strategic Research and Innovation Agenda for Electronic Component and Systems (ECS-SRIA) aims at describing the "Major Challenges, and the necessary Research & Development & Innovation efforts to tackle them, in micro-and nanoelectronics for smart systems integration all the way up to embedded and cyber physical systems, and System of Systems" [START_REF]Electronic Components and Systems -Strategic Research and Innovation Agenda[END_REF].

ECS-SRIA is developed by three industrial associations: AENEAS [2], Inside Industry Association [3] and EPOSS [START_REF]EPoSS (smart-systems-integration[END_REF], with the support of a European team of experts in the discipline coming from European industries, Research and Technology Organizations and Academics in Europe. ECS-SRIA is based on three main areas (see Figure 1):  Key application areas, such as Mobility, (Green colour)  Fundamental Technologies, such as Embedded software and Beyond, (Blue colour)  Cross-Sectional technologies, such as cybersecurity, safety, connectivity, (purple colour).

Figure 1: ECS-SRIA 2022 structure [1]

V ol um e XX, N um ber XX, Mo n th 2022 A da Us er J o ur na l Some of the main ECS-SRIA objectives address:

 To boost industrial competitiveness through interdisciplinary technology innovations [START_REF]Electronic Components and Systems -Strategic Research and Innovation Agenda[END_REF]  To ensure EU digital autonomy through secure, safe and reliable ECS supporting key European application domains [1]  To establish and strengthen sustainable and resilient ECS value chains supporting the Green Deal [START_REF]Electronic Components and Systems -Strategic Research and Innovation Agenda[END_REF]. The following sections discuss some insights related to the chapters of the ECS SRIA in which the authors are more active.

The Cyber-Security and Privacy Challenge

The new generation of CPS are closed to the non-expert human via societal applications. For example, CPS may be customized to individual, including an embedded combination of heterogeneous subsystems of different quality, artificial Intelligence (AI) and connectivity functionality. In this context, ECS-SRIA identifies the cyber-attacks risks and the potential leakage of sensitive data, highlighting a research direction towards a "robust root of trust system, with unique identification enabling security without interruption from the hardware level right up to the applications, including AI involved in the accomplishment of the system's mission in dynamic unknown environments" [START_REF]Electronic Components and Systems -Strategic Research and Innovation Agenda[END_REF]. A particular emphasis is placed on trustworthiness and the hardware security of such systems. Many organizations in Europe, including the IC and Digital System Division, at CEA LIST, and STMicroelectronics, are working in this direction.

The challenge also highlights the need to achieve a common recognized certification scheme under the composition of modular trusted hardware and trusted software. Electric and autonomous vehicles, for example, include several heterogeneous sensors, actuators and embedded devices, that can change or updated over time. Evolution of certification dealing with system changing in some parts over time is a challenge both in terms of cost and security, and, hence, competitiveness.

Ensuring privacy is another relevant theme in the challenge. This focus includes not only privacy-by-design approach, but also of quantum-safe cryptography modules everywhere in the system. Finally, the challenge addresses "ensuring both security and safety" properties. In the case of autonomous vehicle, enhanced connectivity can be useful for driving safety, but, if susceptible to cyber security attacks, may have catastrophic consequences. Ensuring both safety and security requirements has been also capital during the Covid-19 pandemic. Then, in order to achieve a greater level of trustworthiness, it is necessary to properly manage safety and cyber-security issues in the same system.

The Connectivity Challenge

The main ECS SRIA challenges to ensure European leadership in terms of connectivity technologies as well as associated hardware technology supporting the development of connectivity solutions can be summarized as follows:



Innovation Aspects in the 3GPP 5G Standard Evolution

The innovative aspects of 5G networks are not only evolutionary in nature, basically characterized by performance improvements (e.g. capacity, mobility management, connection density, spectral efficiency, latency, energy efficiency), but also present many revolutionary elements (e.g., network softwerization and programmability, network slicing, cloud native approach, more flexible and expanded use of frequencies, artificial intelligence, new approach to cybersecurity). 5G can be considered as a multiplicity of dedicated, flexible, and intelligent networks open to new actors to connect anything and to dynamically serve different verticals, exploiting the concept of slicing. In such a complex scenario the standardization, regulatory and institutional bodies, vendors, and operators should face many challenges to ensure security, safety, privacy, net neutrality and transparency for all users [6]. The third-generation partnership project, 3GPP [6], which is the reference standard for mobile communications, identifies three phases in 5G evolution [7]. During phase 1, the first set of 5G specifications on the new radio have been completed both for non-standalone solutions, when the new radio interacts with an existing 4G/LTE core network, and for the standalone solution, when the new radio interconnects within a 5G core network. The new radio truly represents a step forward from previous generations of radios both in terms of flexibility, throughput, latency, and reliability to meet even mission critical service requirements, exploiting a frequency spectrum that expands from low frequencies below 1 GHz to mmWave, flexible subcarrier spacing up-to 400MHz Channel bandwidths for a single-component (CC) carrier, 3D beamforming for improvement of spectral efficiency, dynamic time-division duplexing (TDD). New radio plays an important role in the autonomous driving application scenario that we will consider in the Section 4. In Phase 2, 3GPP introduced the new concept of Service Base Architecture (SBA) for the 5G core network (Release 16) [8], aiming at providing unprecedented flexibility agilityrespect to the traditional architecture. In Phase 3, Release 17 and beyond, new 5G enhancements are considered both on the radio part (e.g., MIMO, positioning, side-link) and on the core network and slicing, extending it to the access part [START_REF]System architecture for the 5G System (5GS)[END_REF]. New security features have been introduced compared to previous generations, such as a stronger cryptographic algorithm for 256-bit encryption, better air interface security, user privacy protection and enhanced roaming security, but SBA is a completely new concept and leads to completely new security challenges. The SBA is composed of services. Each Network Function (NF) can be regarded as a service [6]. In the upper green part of the Fig. 2 are represented the NFs of the control plane, while the lower part represents the user plane and access networks. It is out of the scope of the present paper entering into details of each network function, (a detailed description can be found in [START_REF]System architecture for the 5G System (5GS)[END_REF]), but for improving the understanding the NFs are, schematically and not rigorously, grouped by functionality.

Fig.2 SBA architecture view, [6]

NFs are self-contained, independent and can be connected to a service bus, using common Internet protocol, [START_REF]Service Based Architecture in 5G[END_REF][START_REF] Borsatti | A Network Slicing Architecture for Mission Critical Communications[END_REF][START_REF]Threat Landscape for 5G Networks[END_REF]. Each NF is a software running on virtual machines or containers efficiently deployed as a Virtual Network Function (VNF) by exploiting cloud computing. NFV and SDN allow the implementation and automation of customized services on a fully programmable platform. Several NFs can be combined to form a logical block, named slice, addressing a specific purpose, with specific QoS and latency features. 3GPP has already standardized some specific slices as Machine Type (mMTC), enhanced Mobile Broadband (eMBB), Ultra Reliable Low Latency (URLLC) and Vehicle-to-X (V2X) communications, but many others can be defined as well. It is undoubtable that the integration of network slicing, NFs, NFV, SDN allows a powerful, flexible, fast, and dynamic service deployment, but its management is complex, and the security should be properly addressed,[6]. 5G threat landscape is complex and continuously evolving, [START_REF]Threat Landscape for 5G Networks[END_REF]. Potential source of attacks (threat actors) could come from end devices, untrusted networks, roaming networks, internet application service providers or insider attack. These attacks constitute different threats to network assets, that can be schematically categorized as loss of availability, confidentiality, integrity, and control [START_REF]Security in 5G Specifications[END_REF].

V2X scenario: communication, safety, and security issues

The automotive industry is at the center of a real revolution. 5G with its low latency, higher bandwidth and great flexibility can provide many previously unachievable features to automotive sector. In recent years vehicles have been equipped with an increasing number of electronic Advanced Driver Assistance Systems (ADAS) developed to increase the level of safety and driving comfort, exploiting four types of sensors: radar, lidar, camera and ultrasonic shortrange sensors, [START_REF] Duncan | The way forward and opportunities towards autonomous driving[END_REF]. However, these sensors have limitations. They are unable to see a pedestrian around the corner or to alert drivers to hazards or slowdowns before they meet on the road: they do not communicate. To complement ADAS systems there are systems called Vehicle to X (V2X) overcoming intrinsic limitation of ADAS systems and, hence, providing features of 360° vision, Non-Light of Sight view, extended range, and the ability to communicate between vehicles, [START_REF]V2X: white paper[END_REF]. There are basically two types of communication for connected vehicles. Dedicated short-range communication (DRSC) is the first standard of V2X technology, based on the IEEE 802.11p standard and operating at 5.9 GHz. It provides the ability for vehicles to communicate with other vehicles and infrastructure around them, exchanging basic safety messages to prevent collisions, but it is not a system for long range communications (about 300 m). Subsequently, 3GPP has introduced the cellular V2X (C-V2X) starting with Release 14, which in turn was based on Release 12 related to Device to Device (D2D) communications. 3GPP specifies four types of C-V2X communication: Vehicle to Vehicle (V2V), Vehicle to Infrastructure (V2I), Vehicle to Network (V2N), and Vehicle to Pedestrians (V2P). The 3GPP's C-V2X standard supports two complementary communication methods: a cellular network-based communication, that uses the air LTE Uu interface and uses the cellular communication bands for long-range transmissions to connect to the network and thus to services; and a direct or side-link that uses the PC5 interface for short-range (less than 1 km) V2V and V2I communications. 3GPP Release 14 defines the foundations of C-V2X communication for basic security message exchange. New features are added in Release 15, but a real step forward is possible with Release16 with the use of 5G New Radio. 5G NR V2X provides lower latency, ultra-reliable communication, and high data rate useful for addressing challenging autonomous driving requirements (e.g. high throughput sensor sharing, intent/trajectory sharing, real time local updates to build detailed maps and share them, coordinated driving etc.), [START_REF]V2X: white paper[END_REF]. In the connected vehicles scenario some critical challenges and possible C-V2X solutions are:  high relative speed of two vehicles driving in opposite directions causes a Doppler shift and frequency offset. C-V2X provides an improved signal design and additional reference signal symbols for better channel estimation.  high vehicle density can lead to radio resource congestion. C-V2X can adopt algorithms that detect available resources, selecting the least congested ones.

 Loss of out-of-coverage synchronization. C-V2X is inherently a synchronous system using GNSS. An application scenario for connected vehicles with challenging requirements on the radio component is platooning, where the lead vehicle of a convoy relays information to the vehicles behind it. Platooning offers many advantages for long-distance travel since vehicles can drive at a constant speed with a short distance in between, reduced air resistance and, hence, a reduction in fuel consumption and CO2 emissions. Moreover, platooning improves safety by, reducing response times and minimizing the risk or, at least, the impact of accidents. ETSI (European TLC Standards Institute) define, the requirements in terms of end-to-end latency, reliability, and data rate, for enhanced V2X scenarios, (e.g. platooning, as driving with high or full automation, remote driving, sharing large amounts of sensor data). It results that 5G NR V2X represents the most suitable solution for Advanced Safety Automated Driving scenarios, as reported in [START_REF]V2X: white paper[END_REF], but the increased connectivity and automation expose them to several cyber threats. [START_REF]ENISA good practices for security of smart cars[END_REF]. ITU-T Recommendation X.1372 provides security guidelines for different V2X scenarios, some depicted in Fig. 3,[18]. Specifically, ITU-T Recommendation X.1372 identifies cybersecurity threats, security requirements to mitigate these threats and describes possible implementations of secure V2X communications. Security issues deal with: i) identification, authenticity to authorize access to services and information, ii) message integrity to ensure that information is accurate and reliable, iii) availability of services and information, iv) confidentiality & privacy of users, their data and their actions from eavesdropping and exploitation and v) non-repudiation and accountability of the source of information. In the case of identification, authenticity, and integrity the main threats to be considered are message manipulation by an attacker providing false information, credential manipulation allowing access to information without authorization, manipulation of sensor data and information, causing traffic congestion, accidents etc., and interception of message and malicious reply in place of the authorized user. In the case of availability, the main threats consist of i) jamming and distributed denial of service (DDoS) attacks caused malicious service requests congesting the channel capacity and so impacting the reliability or availability of C-V2X services, ii) timing attack aiming at delaying the delivery of a safety message to other vehicles and iii) hacking of sensors providing incorrect values determining transient or permanent faults. For confidentiality and privacy issues, the main threats are due to eavesdropping of V2V, V2I, V2P, and V2D messages or leaking of personally identifiable information as identity, position, actions, trajectories of a user of the V2X service. Finally, for non-repudiation case main threats may occur if an attacker manipulates the certification database or accesses the private key to a certificate without authorization. The analysis of possible mitigations to the above-mentioned threats for C-V2X is out of the scope of the presen paper.

Conclusions

In this paper, without claiming to be exhaustive, some insights from the ECS-SRIA and the evolution on 5G networks have been discussed focusing mainly on connectivity, cybersecurity, and safety issues. The cellular-V2X scenario has been analysed as a relevant use case where connectivity, safety and security have a key interwork.

  Strengthening the EU connectivity technology portfolio to maintain leadership, secure sovereignty and offer an independent supply chain.

	 Investigate innovative connectivity technology (new
	spectrum or medium) and new approaches to improving
	existing connectivity technology to maintain the EU's
	long-term leadership
	 Autonomous	interoperability	translation	for
	communication protocol, data encoding, compression,
	security and information semantics.
	 Architectures and reference implementations of
	interoperable, secure, scalable, smart and evolvable IoT
	and SoS connectivity
	 Network virtualisation enabling run-time engineering,
	deployment and management of edge and cloud network
	architectures
	Finally, secure communication and control by powerful
	computation system applying AI are fundamental.

Acknowledgment

This paper has been originated from a lecture the authors gave in the framework of the course 0921SIC10, organized by AMES Society -AEIT (Italian Association of Electrical, Electronics, Automation, Information and Communication Technology). The authors wish to thank AMES -AEIT for supporting this activity.


Claire Guyot

Jean-Paul Barnes

Olivier Renault

Tony Maindron

Surface analysis insight

Keywords: Organic electronic, correlative analysis, ToF-SIMS, XPS, degradation

With the development in the early 2000s of new cluster ion sputtering sources, a reliable analysis of organic multilayers and interfaces has become possible. Nowadays, ToF-SIMS and XPS depth profiling are routinely used to investigate organic stacks. However, analysis beams can cause bond scission or beam-induced degradations that accumulate in buried layers. We developed a correlative protocol that minimises damages related to analysis beams. It uses a shallow angle bevel crater fabricated inside the ToF-SIMS analysis chamber. With this preparation method, the in-depth information is displayed over the surface of the bevel crater. XPS profiles and high-resolution spectra paired with ToF-SIMS images enable an easy identification of the organic layers and complete understanding of their chemistry. The reduction of beam-induced degradation is achieved by minimising the acquisition times, therefore beam exposure on materials. Finally, an important advantage of this preparation method is that the analysis can be performed on exactly the same spot by multiple techniques. Several ToF-SIMS and XPS acquisitions can be carried out with various parameters (investigation of backscattered Arn + cluster ion fragments, tandem MS imaging, …) as well as analysis with other techniques that possess limitations in spatial resolution and/or inaptitude to probe buried layers such as Raman or AFM.

Introduction

When compared to inorganic semiconductors, organic semiconductors offer various specific advantages such as flexibility and transparency of the devices, low material and production costs [START_REF] Sekitani | Stretchable, Large-area Organic Electronics[END_REF][START_REF] Yin | Highly Flexible Fabric-Based Organic Light-Emitting Devices for Conformal Wearable Displays[END_REF][START_REF] Abbel | Roll-to-Roll Fabrication of Solution Processed Electronics[END_REF] . For the past three decades, these materials have been introduced in different applications: from organic light-emitting diodes (OLEDs) to organic photovoltaics (OPVs) and organic field-effect transistors (OFET). Huge advances in the performances and availability of lab-based characterisation techniques enabled the development and commercialisation of these devices over the last decades. In particular, surface characterisation techniques, such as X-rayand ultraviolet photoelectron spectroscopy (XPS, UPS) and time-of-flight secondary ion mass spectrometry (ToF-SIMS), provide chemical information that is essential to understand degradation mechanisms and improve the device performances [START_REF] De Moraes | Influence of the applied charge on the electrochemical degradation in green phosphorescent organic light emitting diodes[END_REF][START_REF] Lin | Migration of small molecules during the degradation of organic light-emitting diodes[END_REF][START_REF] Wang | Photochemical deterioration of the organic/metal contacts in organic optoelectronic devices[END_REF][START_REF] Nakayama | Complete Demonstration of the Valence Electronic Structure Inside a Practical Organic Solar Cell Probed by Low Energy Photoemission[END_REF][START_REF] Jeon | Improved lifetime in organic solar cells using a bilayer cathode of organic interlayer/Al[END_REF] . However, because of their surface sensitivity, important information concerning the bulk, buried layers and interfaces are unreachable. To analyse thin film multilayers, a sputter ion gun is therefore necessary. Monoatomic ion sputtering is inadequate for the analysis of organic layers as the high energy per atom tends to induce a large amount of fragmentation in the bombarded molecules [START_REF] Postawa | Sputtering simulations of organic overlayers on metal substrates by monoatomic and clusters projectiles[END_REF] . With the development of new cluster ion sources in the early 2000s [START_REF] Yamada | Materials processing by gas cluster ion beams[END_REF][START_REF] Weibel | A C60 Primary Ion Beam System for Time of Flight Secondary Ion Mass Spectrometry: Its Development and Secondary Ion Yield Characteristics[END_REF] , a reliable analysis of these organic buried layers and interfaces became possible.

Nowadays, it is possible to carry out combined XPS and ToF-SIMS analysis to obtain a full understanding of the chemistry of organic multilayers by depth profiling. However, there are some limitations to conventional depth profiling when analysing organic semiconductors. Indeed, the analysis beams used in XPS and ToF-SIMS may damage organic molecules. Highenergy ion bombardment or X-ray exposure can cause bond scission or beam-induced degradation [START_REF] Postawa | Sputtering simulations of organic overlayers on metal substrates by monoatomic and clusters projectiles[END_REF][START_REF] Mclaren | Analysis induced reduction of a polyelectrolyte[END_REF][START_REF] Morgan | Revisiting degradation in the XPS analysis of polymers[END_REF] that accumulate in buried layers when depth profiling. Here, we introduce a correlative protocol to analyse organic semiconductors multilayers in a reproducible and reliable way. This protocol serves to minimise the degradation induced by both analysis beams when investigating buried layers. Furthermore, this approach enables an easy correlation of ToF-SIMS and XPS information as the analysis is performed at the same location of the sample, achieved using a bevel crater where the depth information is displayed over the bevel surface.

Bevel crater

In the last decades, different methods to prepare a bevel crater have been investigated. From ultra-low angle microtomy to FIB cross-section or in-situ bevel crater, these methods have demonstrated their ability to provide chemical information when combined with XPS and/or ToF-SIMS [START_REF] Barnes | Focused ion beam time-of-flight secondary ion mass spectroscopy tomography of through-silicon vias for 3D integration[END_REF][START_REF] Fiducia | Understanding the role of selenium in defect passivation for highly efficient selenium-alloyed cadmium telluride solar cells[END_REF][START_REF] Greunz | Cryo ultra-low-angle microtomy for XPS-depth profiling of organic coatings[END_REF][START_REF] Otto | In Situ Investigation of Lithium Metal-Solid Electrolyte Anode Interfaces with ToF-SIMS[END_REF][START_REF] Iida | A new approach for determining accurate chemical distributions using in-situ GCIB cross-section imaging[END_REF] . However, mechanical methods, ion milling and FIB cross-section damage organic semiconductors and destroy valuable molecular information. The development of argon gas cluster ion beam (Ar-GCIB) in the early 2000s enabled the sputtering of organic semiconductors while preserving molecular information. With an energy shared by every atom forming the cluster, the actual energy per atom is very low (few eV/atom) in comparison with conventional monoatomic sputtering (few keV/atom) [START_REF] Matsuo | Size effect in cluster collision on solid surfaces[END_REF][START_REF] Ninomiya | Molecular depth profiling of multilayer structures of organic semiconductor materials by secondary ion mass spectrometry with large argon cluster ion beams[END_REF] . Moreover, due to this low energy per atom, the sputtering direction will be mostly oriented horizontally instead of vertically in the case of monoatomic sputtering [START_REF] Postawa | Sputtering simulations of organic overlayers on metal substrates by monoatomic and clusters projectiles[END_REF] , inducing a lower implantation of the argon atoms inside the organic layers.

Different papers reported the use of an Ar-GCIB to fabricate bevel structures. Mihara et al. demonstrated the combination of monoatomic gun with cluster sputtering to form a bevel crater [START_REF] Mihara | Embedding-Free Method for Preparation of Cross-Sections of Organic Materials for Micro Chemical Analysis Using Gas Cluster Ion Beam Sputtering[END_REF] . First, a mask of damaged organic semiconductors is created with a high-energy monoatomic beam, creating a region where the sputter yield has two orders of magnitude less than pristine organic molecules. The surface is then sputtered with the argon gas cluster ion beam enabling the fabrication of a well-controlled bevel crater as the mask protect underlying material from sputtering. A fine control of the dimension of the crater is achieved but the angle of the bevel crater cannot be modified as the angle of the sputter gun in the analysis chamber is fixed (45°). Another method relies on Zalar rotation while sputtering to form a circularly symmetrical crater, exposing the different layers from the center to the edges of the crater [START_REF] Yun | Fabrication of Bevel Structures by Means of the Spatial Gradient of Ion Dose in Ar Gas Cluster Ion Beam and Its Unique Characteristics[END_REF] . This approach is easy to implement and available on many instruments. However, it does not provide fine control over the dimensions and angle of the bevel crater, which needs to be extremely small for amplifying ultra-thin buried layers.

A bevel crater with total control over its dimensions and angle can be fabricated varying the fluence of an ion beam over the surface of a sample as demonstrated by Mao et al. with a sputtering of C60 + 23 . To fabricate this crater, a FIB-like cross section is operated inside the ToF-SIMS analysis chamber by varying the fluence of the sputtering beam over the surface. Therefore, the in-depth information is converted into a lateral information over the bevel surface (figure 1); hence converting a hundred nanometres thick organic stack into a distance of a few hundred micrometres over its surface. Depending on the dimension of the bevel crater, the grazing angle of the bevel crater (up to 0.005°) can lead to an amplification factor of the organic layers up to 10 000. For instance, a thin layer of 5 nm depth would appear with a dimension of about 50 µm on the bevel surface for a 1000 µm-wide crater. This amplification of the apparent thickness of organic layers enables an easy correlative analysis between ToF-SIMS and XPS, enabling analysis over the same area for both techniques. Figure 1. 3D scheme of the bevel crater fabrication process. By applying an increasing ion fluence along the y-direction of the raster area, a bevel-shape crater is formed. Here, an argon-cluster ion beam is used as the sputter gun with an energy of 5 keV and a cluster size of 7500 atoms of argon. In this case, the dimension of the crater is 800 x 800 µm² and the maximum depth is about 100 nm, leading to an angle of 0.007° and an amplification coefficient of the organic layers of about 8700.

ToF-SIMS and XPS measurements on organic light-emitting diode multilayer stacks

An example of correlative ToF-SIMS and XPS measurements on a bevel crater is presented on organic light-emitting diode (OLED) multilayers.

The studied green top-emission OLED consists of a five-layer electrically-doped structure sandwiched between an anode of Al : Cu / TiN and a 15 nm-thick silver cathode (figure 2). 1,4,5,8,9,11-Hexaazatriphenylenehexacarbonitrile (HATCN) was used as the 1 nm-thick hole injection layer (HIL). The hole transport layer (HTL) of 30 nm consists of a matrix of 2,2',7,7'-Tetra(N,N-di-p-tolyl)amino-9,9-spirobifluorene (STTB) doped with 4% of 2,3,5,6-Tetrafluoro-7,7,8,8-tetracyanoquinodimethane (F4TCNQ). The electron blocking layer (EBL) consists of 6 nm of N,N′-Bis(naphthalen-1-yl)-N,N′-bis(phenyl)benzidine (NPB). In the emissive layer (EML) of the OLED device, we used 9,9′-(4,4′-(Phenylphosphoryl)bis-(4,1 phenylene))bis(9Hcarbazole) (BCPO) as host and tris[2-(p-tolyl)pyridine]iridium(III) (Ir(mppy)3) as triplet dopant, with a doping of 15%. 2,2′,2"-(1,3,5-Benzinetriyl)-tris(1-phenyl-1-H-benzimidazole) (TPBi) is used as the 5 nm-thick hole blocking layer (HBL). Finally, 4,7-Diphenyl-1,10phenanthroline (BPhen) and Ca were used as host and dopant of the 20 nm-thick electron transport layer (ETL) respectively, with a doping of 4% [START_REF] Shelhammer | Doping effects and stability of calcium in organic electron-transport materials[END_REF] . The devices are then encapsulated with 80 nm-thick SiO and 25-nm thick Al2O3. All organic semiconductors, cathode and SiO layers were deposited by physical vapor deposition while Al2O3 was deposited by low temperature atomic layer deposition [START_REF] Maindron | Defect analysis in low temperature atomic layer deposited Al2O3 and physical vapor deposited SiO barrier films and combination of both to achieve high quality moisture barriers[END_REF] . Before introducing the sample mounted on the metallic sample holder with double-sided conductive tape in the ToF-SIMS analysis chamber (TOF-V from IONTOF GmbH), the encapsulation layers and the cathode are peeled off using an adhesive tape under protective atmosphere (N2) [START_REF] Langer | Advanced chemical characterization of organic electronic materials and devices[END_REF] . The bevel crater is then formed in the organic layers using the method described previously. The sputtering is realized with Ar7500 + ions at an energy of 5 keV with a current around 60 pA. The size of the studied bevel crater is 800 x 800 µm² for a maximum depth of about 100 nm (thickness of the organic stack). As mentioned previously, the layers are extended by the lateral spreading enabled by the bevel crater. Here, the amplification factor is about 8700: for instance, the 5 nm-thick TPBi layer will be converted into a distance of about 43.5 nm on the bevel surface. Finally, the size of the crater has been optimised in order to obtain high intensity with minimum analysis beam exposure (in particular for XPS measurements) as well as fabricating a bevel crater with a controlled raster.

Once the bevel crater is fabricated, ToF-SIMS measurements are first performed as the sample is already present inside the analysis chamber. An image of the surface is acquired with Bi3 + primary ions with an energy of 15 keV over an area of 1.1 x 1.1 mm². A wide area is analysed to identify the beginning of the bevel surface. The images were acquired in both polarities with a dose limit of 4.10 11 ion/cm², ensuring static ToF-SIMS conditions. To each pixel corresponds a mass spectrum, therefore, the identification of the different layers is based on the spatial position of the characteristic molecular ions over the bevel surface as presented in figure 3. a). It is also possible to access to an intensity profile of each characteristic peak (figure 3. b)). To differentiate the bevel surface from the sample surface, common contamination ions of the surface such as Na + , K + or Cs + as well as ions from the diffusion of cathode or encapsulation layer into the organic layers can be used. Here, the surface is identified by the contamination ion Cs + remaining in the analysis chamber from previous depth profiling. Its profile confirms that the crater is about 800 µm-long. The exposed distance of the BPhen : Ca layer on the bevel crater is shorter than expected as a non-negligible thickness is removed when the encapsulation layers and the cathode are peeled off [START_REF] Shin | Degradation mechanism of organic light-emitting device investigated by scanning photoelectron microscopy coupled with peel-off technique[END_REF] . The sample is then transferred to the XPS analysis chamber (PHI5000 VersaProbe II from Physical Electronics Inc.) under an inert atmosphere (N2, P>Patm) using a transfer vessel. The sample was mounted with double-sided conductive tape onto the metallic sample holder. First, the bevel crater is identified by the difference of image contrast between the organic layers and the substrate on the secondary electron image (figure 4. a)). No chemical contrast amongst the different organic layers is observed as each molecule is mainly composed of carbon and nitrogen. To identify the positions of the layers over the 800 µm-long bevel crater, a profile along a line (linescan) is acquired with a pass energy of 93.3 eV. All XPS measurements are performed using a monochromatized Al Kα source (1486.6 eV) with a beam size of 50 µm and no charge compensation. As the atomic concentration of carbon and nitrogen is very similar in each layer, the organic layers are identified thanks to the profiles of both the characteristic inorganic dopants and substrate. Calcium, iridium and titanium are present in the electron transport layer, emissive layer and the substrate respectively. The positions of the other layers are extrapolated from the dopants and substrate profiles (figure 4. b)). For instance, the TPBi layer, sandwiched between the electron transport layer and the emissive layer, is located at the crossing of the calcium and iridium profiles in the linescan. Afterwards, survey spectra are acquired in each layer. It confirms the correct location of the layers and the elements presents at each position. Finally, high-resolution core-level spectra of and band region spectra are acquired to access the chemical environment of each element. High-resolution spectra and valence band region spectra are acquired with a pass energy of 23.5 eV and 46.95 eV resulting in an overall energy resolution (X-ray bandwidth included) of 0.67 and 0.83 eV respectively.

Valence band region spectra are smoothened using linear Savitzky-Golay smoothing with 5 datapoints. The high-resolution spectra and valence band region spectrum of the emissive layer are presented in figure 4. c). A comparison with reference spectra has been performed on each spectra. However, only the comparison of the valence band region is presented. Indeed, the valence band region of organic compounds acts as a fingerprint [START_REF] Drummond | Practical polymer film characterisation using high performance XPS methods[END_REF][START_REF] Pireaux | X-ray photoemission study of core-electron relaxation energies and valence-band formation of the linear alkanes. II. Solid-phase measurements[END_REF] . Slight chemical modifications may not be observed on core-level spectra but they may be visible in the valence band region. Here, except a shift to lower binding energies (≈ 0.7 eV), the valence band regions are similar regarding the energy separation and relative intensities of its characteristic spectral components. This shift is also observed for all core-level and can be explained by band bending in multilayer structures. Therefore, it is concluded that the fabrication of the bevel crater and the subsequent ToF-SIMS measurements have not induced any detectable damage to the organic molecules.

What are the benefits for the analysis of organic multilayers with this method?

This protocol serves multiple goals. First, the use of a bevel crater provides an easy correlation of ToF-SIMS and XPS results on the same area by locating and individually analysing each layer over the bevel surface, as demonstrated previously.

The second objective concerns the degradations induced by analysis beams. Indeed, modifications of the molecule chemistry and bond scission can occur under illumination and charge neutralization. Here, the described method is not employing any neutralizer whatsoever. In ToF-SIMS, once the Bin q+ primary ions reach the surface, the cluster will dissociate and penetrate the layers over a distance of about 10-20 nm [START_REF] Muramoto | Exploring the Surface Sensitivity of TOF-Secondary Ion Mass Spectrometry by Measuring the Implantation and Sampling Depths of Bin and C60 Ions in Organic Films[END_REF] . The resulting collision cascade into the organic semiconductors can cause bond scission and fragment organic molecules [START_REF] Postawa | Sputtering simulations of organic overlayers on metal substrates by monoatomic and clusters projectiles[END_REF] . In XPS, X-ray-induced damages on organic semiconductors have been well documented. Changes in the chemical environment of elements and the ratio of atomic concentration can be observed after an hour of X-ray exposure [START_REF] Mclaren | Analysis induced reduction of a polyelectrolyte[END_REF][START_REF] Morgan | Revisiting degradation in the XPS analysis of polymers[END_REF] . Moreover, X-rays probe about 9 nm depth into organic semiconductors. As these characterisation techniques are surface sensitive, depth profiling is mandatory to access in-depth information. However, analysis beams induce damages over a depth of about 10 nm, thus buried layers may already be modified before their analysis. With the described preparation method, the in-depth information is displayed over the surface of the bevel without previous analysis beam exposure, preserving the chemical information of buried layers. Furthermore, acquisition times, and therefore analysis beam flux, are reduced in comparison with depth profiling. In ToF-SIMS, images in both polarities are acquired with a primary ion dose of 4.10 11 ions/cm², ensuring static conditions. In XPS, even though several measurements are performed, the total X-ray exposure time varies from 60 min to 120 min. First, a profile along an 850 µm-length line to identify each layer lasts about 70 min for 85 pixels, resulting in an effective time per pixel of 50 s. Then, survey spectra are acquired in each identified layers for about 3 min. Finally, high-resolution core-level spectra and valence band region are performed for about 50 min to 110 min depending on the number of elements to acquire.

Finally, several measurements can be performed on a bevel crater with different ToF-SIMS and XPS parameters or with other characterisation techniques. Firstly, changing incident energy or primary ion type induces major changes in the ToF-SIMS mass spectrum. With a higher primary ion incident energy, more fragmentation peaks will be present in the mass spectrum. This phenomenon is also observed for smaller cluster or monoatomic primary ions such as Ga + and Bi1 + . Argon cluster Arn + can also be used as primary ions to probe physical information by studying the backscattering ratio of small argon cluster [START_REF] Mochiji | Mass spectrometric analysis of the dissociation of argon cluster ions in collision with several kinds of metal[END_REF][START_REF] Chundak | Probing the Surface Curie Temperature of Ferroelectric P(VDF-ran-TrFE) Copolymers by Argon Gas Cluster Ion Scattering[END_REF][START_REF] Poleunis | Temperature Dependence of Arn+ Cluster Backscattering from Polymer Surfaces: a New Method to Determine the Surface Glass Transition Temperature[END_REF] . Tandem MS imaging of the bevel crater can be used to resolve ambiguities [START_REF] Fisher | A New Method and Mass Spectrometer Design for TOF-SIMS Parallel Imaging MS/MS[END_REF] or help identify products of degradation [START_REF] Iida | Time-of-flight secondary ion tandem mass spectrometry depth profiling of organic light-emitting diode devices for elucidating the degradation process[END_REF] . Secondly, the large dimensions and the grazing angle of the crater enables performing other measurements as well, for instance using characterisation techniques with lower spatial resolution such as Raman [START_REF] Giorgetti | A study of the degradation of poly(3-octylthiophene)-based light emitting diodes by Surface Enhanced Raman Scattering[END_REF] or techniques with limitations to probe buried layers such as AFM [START_REF] Kolosov | Direct observation of structural changes in organic light emitting devices during degradation[END_REF] .

Conclusions

In summary, this paper presents a correlative approach to analyse organic multilayers with minimal beam damage using a bevel crater. Fabricated by varying the fluence of the Ar cluster ions along the y-direction of the raster area over the surface of the sample, the bevel structure exposes the organic layers over its surface. The grazing angle (~0.007°) and the dimensions (800 µm long for a depth of 100 nm) of the bevel crater enable an amplification coefficient of the organic stack of about 8700. An easy correlation is achieved by coupling ToF-SIMS images with XPS profiles and high-resolution spectra. Profiles along the crater provide depth information by locating the organic layers over the bevel surface. The correlation of individual analysis of organic layers by ToF-SIMS and XPS offers a more complete understanding of the chemistry of organic molecules. A reduction of beam-induced degradation is achieved by minimising the acquisition times, therefore analysis flux. Moreover, compared to depth profiling, the buried layers displayed over the bevel surface have not been previously exposed to analysis beams, ensuring the preservation of their chemical information. Several measurements can be performed on such structure: variation of ToF-SIMS and XPS parameters as well as analysis with other characterisation techniques with limitation in spatial resolution and /or inaptitude to probe buried layers are possible. Finally, this method can be applied to various organic multilayers samples and may be of particular interest for non-homogenous samples such as dye-sensitized solar cells or polymer light-emitting electrochemical cells, where it is critical to analyse the same area of the sample to avoid variations due to sample heterogeneity.

Figure 2 .

 2 Figure 2. a) Schematic diagram of the device structure and b) utilized organic molecules in this study.

Figure 3 .

 3 Figure 3. a) Secondary ion image of the surface of the bevel crater and identified layers based on the position of characteristic molecular ions, Ti + and Ca + ions. The surface is identified by the surface contamination ion Cs + ; b) mean intensity profile of the characteristic molecular ions over the selected area (yellow rectangle).

Figure 4 .

 4 Figure 4. a) secondary electron image of the bevel crater superposed with its complementary image in ToF-SIMS, the XPS linescan (white) and the analysis spot of each layer (coloured dots), b) Ca 2p, Ir 4f and Ti 2p profiles over the surface of the bevel crater along a line and identification of the position of each layer and c) high resolution core-level spectra and valence band region of BCPO : Ir(mppy)3 layer.

Acknowledgments

This work, carried out on the Platform for NanoCharacterisation (PFNC), was supported by the "Recherches Technologiques de Base" program of the French National Research Agency (ANR).