Angel Leon

Anaïs Barasinski

Emmanuelle Abisset-Chavanne
email: emmanuelle.abisset-chavanne@ec-nantes.fr

Elias Cueto
email: ecueto@unizar.es

Anais Barasinski
email: anais.barasinski@ec-nantes.fr

Francisco Chinesta
email: francisco.chinesta@ensam.edu

Wavelet-based multiscale proper generalized decomposition

Keywords:

come

Introduction

Model Order Reduction -MOR -techniques allow nowadays solving, under real-time constraints, complex models. Intense research activities allowed reaching at present a certain maturity in the domain of model order reduction. Among the numerous references, the interested reader can refer to some review papers and books [START_REF] Chinesta | A short review in model order reduction based on Proper Generalized Decomposition[END_REF][START_REF] Chinesta | Parametric PGD based computational vademecum for efficient design, optimization and control[END_REF][START_REF]Separated Representations and PGD Based Model Reduction: Fundamentals and Applications[END_REF][START_REF] Chinesta | Model order reduction[END_REF], covering three major MOR technologies: POD (Proper Orthogonal Decomposition), RB (Reduced Basis), and PGD (Proper Generalized Decomposition).

Proper Orthogonal Decomposition (POD) is a general technique for extracting the most significant characteristics of a system's behavior and representing them in a set of "POD basis vectors". These basis vectors then provide an efficient (typically low-dimensional) representation of the key system behavior, which proves useful in a variety of ways. The most common use is to project the system-governing equations onto the reduced-order subspace defined by the POD basis vectors. This yields an explicit POD reduced model that can be solved in place of the original system. The POD basis can also provide a low-dimensional description in which to perform parametric interpolation, infill missing or "gappy" data, perform model adaptation, or define hyper-reduction procedures [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF]. There is an extensive literature and POD has seen broad application across fields. Some review of POD and its applications can be found in [START_REF] Volkwein | Model Reduction Using Proper Orthogonal Decomposition[END_REF][START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF].

Another family of model reduction techniques lies in the use of Reduced Basis constructed by combining a greedy algorithm and "a posteriori" error indicators. As for the POD, the Reduced Basis method requires some amount of offline work, but then the reduced basis model can be used online for solving different models with control of the solution accuracy, because of the availability of error bounds. When the error is unacceptably high, the reduced basis can be enriched by invoking a greedy adaption strategy [START_REF] Patera | Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: application to transport and continuum mechanics[END_REF].

Separated representations, at the heart of the so-called Proper Generalized Decomposition methods, are considered when solving at-hand partial differential equations by employing procedures based on the separation of variables. Then they were considered in quantum chemistry for approximating multidimensional quantum wave-function. In the 1980s, Pierre Ladevèze proposed the use of space-time separated representations of transient solutions involved in strongly nonlinear models, defining a non-incremental integration procedure [START_REF] Ladevèze | The large time increment method for the analyze of structures with nonlinear constitutive relation described by internal variables[END_REF][START_REF] Ladevèze | Nonlinear Computational Structural Mechanics. New Approaches and Non-incremental Methods of Calculation[END_REF]. Later, separated representations were employed for solving multidimensional models suffering the so-called curse of dimensionality [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representation[END_REF], as well as in the context of stochastic modeling [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF]. Then, they were extended to the separation of space coordinates, making possible the solution to models defined in degenerated domains, e.g., plate and shells [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF], as well as for addressing parametric models where model parameters were considered as model extra-coordinates, making possible the offline calculation of the parametric solution, which can be viewed as a metamodel or a computational vademecum to be used online for real-time simulation, optimization, inverse analysis, and simulation-based control [START_REF] Chinesta | Parametric PGD based computational vademecum for efficient design, optimization and control[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF].

Separated representations

Within the PGD framework, four kinds of separated representations have been widely considered.

(i) Space-time separated representations that allowed the construction of efficient incremental and non-incremental integrators.

Within the standard finite element method, a space-time solution u(x, t), x ∈ ⊂ R 3 and

t ∈ I ⊂ R, of a transient problem is approximated from u(x, t) ≈ M i=1 u(x i , t)N i (x) ( 1
)
where M is the number of nodes employed for interpolating the unknown field, located at positions x i , and N i (x) the so-called shape functions. Because of the interpolative property of the shape functions, the approximation coefficients correspond to the nodal value of the approximated field, u(x i , t). Thus, in general, when solving a nonlinear problem, at least a linear system of size M must be solved at each time step. When considering P time steps (P can reach several millions), the complexity grows very fast. When considering POD-based model order reduction, the solution is projected into the reduced basis composed of functions {φ 1 (x), • • • , φ R (x)} extracted from some collected snapshots of the problem solution, with in general R M, and consequently the solution approximation reads

u(x, t) ≈ R i=1 ξ i (t)φ i (x) (2)
which requires the resolution of linear systems of size R instead of the ones of size M characteristic of finite element solutions. The use of a reduced basis implies in many cases impressive computing-time savings. Approximations (1) or (2) imply a finite sum of time-dependent coefficients and space functions. The last are assumed known; they consist of the usual finite element shape functions or the modes extracted by applying, for example, Proper Orthogonal Decomposition -POD. A step forward could consist in assuming space functions to be also unknown and in computing both time and space functions on the fly. In this case, the approximation reads

u(x, t) ≈ N i=1 T i (t)X i (x) (3)
Because both functions involved in approximation (3) are unknown, it defines a nonlinear problem whose solution requires an appropriate linearization strategy. The interested reader can refer to [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF] and the references therein for additional details on the separated representation constructor.

Expression (3) evidences that the solution procedure requires the resolution of about N problems, with N M and N ∼ R (in fact a bit more because of the nonlinearity induced by separated representations) involving the space coordinates (in general three-dimensional -3D -and whose associated discrete systems are of size M) for computing the space functions X i (x) and about N one-dimensional -1D -problems for calculating the time functions T i (t). Due to the fact that the computing cost related to the solution of 1D problems is negligible with respect to the solution of 3D problems, the resulting computational complexity reduces drastically, scaling with N instead of P.

(ii) Space separation allowed addressing multi-physics problems defined in degenerated geometries in which at least one of its dimensions remains much smaller that the other ones (e.g., beams, plates, shells, laminates...) or processes involving additive layers (e.g., automated tape placement, 3D printing, or additive manufacturing).

If domain can be decomposed as = x × y × z , the solution u(x, y, z) could be approximated by using the separated representation

u(x, y, z) ≈ N i=1 X i (x)Y i (y)Z i (z) (4)
which allows calculating the 3D solution from a sequence of 1D problems. For some geometries, as the ones associated with plates or shells, in-plane-out-of-plane separated representation becomes specially suitable,

u(x, y, z) ≈ N i=1 X i (x, y)Z i (z) ( 5
)
where the 3D complexity is reduced to the one characteristic of 2D problems, the ones related to the calculation of in-plane functions X i (x, y).

(iii) Space-time-parameter separated representations allowed constructing the so-called computational vademecums (also known as abacus, virtual charts, nomograms...) efficiently considered for multiple purposes: simulation, optimization, inverse analysis, uncertainty propagation and simulation-based control, all of them under the real-time constraint.

When the unknown field involves space, time, and a series of parameters μ 1 , . . . , μ Q , its associated separated repre-

sentation reads u(x, t, μ 1 , . . . , μ Q ) ≈ N i=1 X i (x)T i (t) Q j=1 j i (μ j ) (6)
(iv) Separated representations of intrinsically multidimensional models involving differential operators applying on time, space, and a series of conformation coordinates c 1 , . . . , c C . In this case, the solution is approximated according to

u(x, t, c 1 , . . . , c C ) ≈ N i=1 X i (x)T i (t) C j=1 C j i (c j ) (7)

Paper motivation

This work focuses on the separated representations previously discussed and contributes in two different ways. First, in general the calculation of functions involved in the finite sum representations is performed on a given approximation basis. Thus, the possible multi-scale character of the individual modes (functions involved in the finite sum) and/or the reconstructed solution (the finite sum itself) is not addressed.

Even if today it is accepted within the PGD community that both the number of terms in the finite sum, and the number of nodes used for discretizing each function involved in the former are of major relevance for ensuring accuracy, and even if different error estimators have been proposed in this context [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Nadal | A separated representation of an error indicator for the mesh refinement process under the Proper Generalized Decomposition framework[END_REF], the multi-scale nature of both, the terms in the sum and the sum itself, remain almost unexplored. An adaptive progressive PGD strategy was proposed in [START_REF] Falcó | Proper generalized decomposition for nonlinear convex problems in tensor Banach spaces[END_REF].

This work considers the use of a naturally multi-scale approach, a wavelet-based approximation technique, for approximating the different functions involved in the separated representation. Then, recombining the different scales involved in this approximation allows one to capture the multi-scale character of the solution. Even if, in our knowledge, this kind of approximations was employed in a Galerkin setting, as referred to later, its consideration within the PGD framework remains unexplored.

The next section revisits the PGD constructor and the wavelet-based multi-resolution analysis that will be applied in Section 3 within the PGD constructor described above. Then Section 4 will illustrate some potential applications of the proposed numerical technique.

Methods

PGD constructor at a glance

Consider the solution to the Poisson equation u(x, y) = f (x, y) [START_REF] Patera | Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations[END_REF] in a two-dimensional rectangular domain = x × y = (0, L) × (0, H). We specify homogeneous Dirichlet boundary conditions for the unknown field u(x, y), i.e. u(x, y) vanishes at the domain boundary . Furthermore, we assume that the source term f is constant over the domain .

For all suitable test functions u * , the weighted residual form of (8) reads

x × y u * ∂ 2 u ∂x 2 + ∂ 2 u ∂ y 2 -f dx d y = 0 (9)
Our goal is to obtain a PGD approximate solution to [START_REF] Patera | Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations[END_REF] in the separated form

u(x, y) = N i=1 X i (x) • Y i (y) ( 10
)
We shall do so by computing each term of the expansion one at a time, thus enriching the PGD approximation until a suitable convergence criterion is satisfied. Thus, at each enrichment step n (n ≥ 1), we have already computed the n -1 first terms of the PGD approximation

u n-1 (x, y) = n-1 i=1 X i (x) • Y i (y) (11)
We now wish to compute the next term

X n (x) • Y n (y) to obtain the enriched PGD solution u n (x, y) = u n-1 (x, y) + X n (x) • Y n (y) = n-1 i=1 X i (x) • Y i (y) + X n (x) • Y n (y) (12)
Both functions X n (x) and Y n (y) are unknown at the current enrichment step n, and they appear in the form of a product.

The resulting problem is thus non-linear, and a suitable iterative scheme is required. We shall use the index p to denote a particular iteration. At enrichment step n, the PGD approximation u n,p obtained at iteration p thus reads

u n,p (x, y) = u n-1 (x, y) + X p n (x) • Y p n (y) (13)
and the algorithm proceeds by (i) calculating

X p n (x) from Y p-1 n (y), then (ii) updating Y p n (y) from the just-computed X p n (x);
and finally (iii) checking the convergence by evaluating

X p n (x) • Y p n (y) -X p-1 n (x) • Y p-1 n (y) .
In the first step, by introducing the approximate

u n,p (x, y) = u n-1 (x, y) + X p n (x) • Y p-1 n
(y) [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF] as well as the weight function

u * (x, y) = X * n (x) • Y p-1 n (y) ( 15
)
into the problem integral form (9), after integrating over y , one finds

x X * n • α x d 2 X p n dx 2 + β x X p n dx = - x X * n • n-1 i=1 γ x i d 2 X i dx 2 + δ x i X i dx + x X * n • ξ x dx (16)
where the coefficients result from the integration over y of functions depending on the y-coordinate. For more details, the interested reader can refer to [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF]. The above one-dimensional problem is then discretized by using a standard mesh-based discretization technique, as for example the finite element method.

Having thus computed X

u n,p (x, y) = n-1 i=1 X i (x) • Y i (y) + X p n (x) • Y p n (y) ( 17
)
which finally leads to the one-dimensional problem involving the unknown field Y p n (y):

y Y * n • α y d 2 Y p n d y 2 + β y Y p n d y = - y Y * n • n-1 i=1 γ y i d 2 Y i d y 2 + δ y i Y i d y + y Y * n • ξ y d y ( 18
)
It is important to realize that the original two-dimensional Poisson equation defined over = x × y has been transformed within the PGD framework into a series of decoupled one-dimensional problems formulated in x and y . The procedure can be generalized when considering D-dimensional problems. As previously indicated, the functions X n (x) and Y n (y) are approximated using standard finite element approximations or orthogonal polynomials (spectral approximations). Thus, in first approximation one could expect that the accuracy of the computed solution should depend on the number of terms involved in the finite sum N and the number of nodes M (or the polynomial degree) considered for approximating the different functions involved in the finite sum, X n (x) and Y n (y) in the previous case study.

In [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF], the authors considered the Poisson equation ( 8) on a two-dimensional rectangular domain

= x × y = (0, 2) × (0, 1) with f = 1, which has as exact solution u ex (x, y) = m,n odd 64 π 4 4n 2 + m 2 sin mπx 2 sin (nπy) (19)
Functions X i (x) and Y i (y) were sought on a uniform one-dimensional grid with M points. The solution convergence (with respect to the exact one) was evaluated by using the error

E M (u N M ) = 1 0 2 0 u ex (x, y) -u N M (x, y) 2 dx d y ( 20
)
where integrals were performed numerically and u N M refers to the solution obtained using a finite sum composed of N terms and the involved one-dimensional functions were approximated by using M nodes. Fig. 1 depicts the PGD error E M (u N M ) as a function of the number of enrichment steps N for different numbers M of discretization points. This figure proves that both M and N are relevant to ensure convergence, and that both must be considered together, because, for a given mesh, and independently of the number of terms considered in the finite sum, the error reaches a plateau. In the other sense, even a very fine mesh is unable to capture the solution if the number of terms in the finite sum is insufficient.

Wavelet-based multi-resolution analysis and adaptive approximations

Multi-resolution analysis is based on the construction of a series of embedded subspaces

V j ⊂ V j+1 , ∀ j ∈ Z, i.e. {∅} • • • ⊂ V j-1 ⊂ V j ⊂ V j+1 • • • ⊂ L 2 (R) (21)
where each subspace is spanned by the integer translation of a single function, the scaling function φ(x) [START_REF] Gopalakrishnan | Wavelet Methods for Dynamical Problems with Application to Metallic, Composite, and Nano-Composite Structures[END_REF].

If φ(x) ∈ V 0 , the functions φ 0,k = φ(xk) constitute a basis of V 0 , and the functions φ j,k (x) = 2 j/2 φ(2 j xk) define a basis of subspace V j . Thus, the projection of the function f (x) in V j reads

P j f (x) = k=+∞ k=-∞ c j,k φ j,k (x) (22) The fact that V 0 ⊂ V 1 allows writing φ(x) = k=+∞ k=-∞ a k φ(2x -k) (23)
and because V j ⊂ V j+1 we can define the orthogonal complement W j , such that [START_REF] Kumar | A three-step wavelet Galerkin method for parabolic and hyperbolic partial differential equations[END_REF] Spaces W j are also spanned by the integer translation of a single function, the wavelet function ψ(x), defining the bases at each scale:

V j+1 = V j ⊕ W j (24) with V j ⊥ W j . An important consequence is that j∈Z W j = L 2 (R) (25) or V 0 ⊕ j∈N W j = L 2 (R)
ψ j,k = 2 j/2 ψ(2 j x -k). Because W 0 ⊂ V 1 , we can write ψ(x) = k=+∞ k=-∞ b k φ(2x -k) (27)
In what follows we are considering Daubechies wavelets [START_REF] Daubechies | Orthonormal basis of compactly supported wavelets[END_REF] for approximating the different functions involved in the separated representation of the problem solution. This choice is motivated by their compact support and regularity (approximation consistency). In opposition to usual approximation techniques, its use allows space and frequency localization, the latter because of its multiscale character and the former due to its compact support (in opposition to approximations based on Fourier polynomials). The approximation is fully defined by the coefficients a k in Eq. ( 23). The coefficients b k depend on

a k from b k = (-1) k a 1-k .
The coefficients a k result from the following conditions:

(i) the normalization condition, ∞ -∞ φ(x) dx = 1 (28)
which, taking into account Eq. ( 23), leads to (iii) if we consider a choice with N coefficients, the previous conditions only provide N /2 + 1 equations. Thus, other N /2 -1 conditions should be added in order to compute the N coefficients. One common route consists in enforcing that the scaling function exactly represents polynomials of order M = N /2, i.e.

f (x) = a 0 + a 1 x + • • • + a M-1 x M-1 (32)
from which, by using

f (x) = k=∞ k=-∞ c k φ(x -k) ( 33
)
and the fact that the wavelet function ψ(x) is orthogonal to the scaling translates φ(xk), it results

f (x), ψ(x) = k=∞ k=-∞ c k φ(x -k), ψ(x) = 0 (34) with g(x), h(x) = ∞ -∞ g(x) h(x) dx (35)
By replacing f (x) by its expression (32) into Eq. (34) leads to

a 0 ψ(x), 1 + a 1 ψ(x), x + a M-1 ψ(x), x M-1 = 0 (36)
which applies for any value of coefficients a 0 , • • • , a M-1 . Thus, by choosing a l = 1 and a m = 0, ∀m = l, it results ψ(x),

x l = 0, l = 0, 1, • • • , M -1 (37)
implying that the M first moments of the wavelet function vanish. This condition can be proved to be equivalent to

k=∞ k=-∞ (-1) k a k k l = 0, l = 0, 1, • • • , M -1 (38)
Now, all the conditions above allow determining the N coefficients that determine completely the approximation.

In order to discretize the weak form of a given problem, one must have access to the scaling and wavelet function and their derivatives at any point (in particular at the integration points considered for evaluating the integrals involved in the problem weak form). For that purpose, we start by considering

φ(x) = a 0 φ(2x) + a 1 φ(2x -1) + • • • + a N -1 φ(2x -N + 1) (39)
or, taking into account that the support of the scaling function is

[0, N], it results ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ φ(0) = a 0 φ(0) φ(1) = a 0 φ(2) + a 1 φ(1) + a 2 φ(0) φ(2) = a 0 φ(4) + a 1 φ(3) + a 2 φ(2) + a 3 φ(1) + a 4 φ(0) . . . ( 40
)
which defines an eigenvalue problem giving access to the value of the scaling function at the integers. Now, in order to obtain its value at the so-called dyadic points, we make use of

φ(x/2) = k=∞ k=-∞ a K φ(x -k) (41)
and from those we can compute φ(x/4) and so on. Now, for the derivatives, we proceed in a similar manner, by taking the derivative of Eq. (39), which allows us, by following the same rationale and solving the associated eigenproblem, to calculate the derivative at the integer points, and from this the one at the dyadic points, and so on.

Approximations based on wavelets have been successfully considered for discretizing partial differential equations, most of the time by using Galerkin formulations [START_REF] Gopalakrishnan | Wavelet Methods for Dynamical Problems with Application to Metallic, Composite, and Nano-Composite Structures[END_REF][START_REF] Avudainayagam | Wavelet-Galerkin method for integro-differential equations[END_REF][START_REF] Jones | The wavelet-Galerkin method for solving PDE's with spatially dependent variables[END_REF][START_REF] Kumar | A three-step wavelet Galerkin method for parabolic and hyperbolic partial differential equations[END_REF][START_REF] Mahmoudi | Wavelet Galerkin method for numerical solution of nonlinear integral equations[END_REF].

If, for a while, we consider the one-dimensional problem involving the unknown field u(x) and similarly for x = 0.85, with a small enough constant. When using wavelet-based approximations, the unknown field can be expressed at the lowest level by using the scaling functions

d dx K (x) du dx = 0 (42) with x ∈ [0, 1], u(x = 0) = 0, u(x = 1) = 1 and the model parameter K (x) defined from K (x) = K max if x ∈ [0.
u 0 (x) = k c 0,k φ(x -k) (45)
The problem solution (described later) consists of the set of coefficients c 0,k . The solution at the next level is written by combining both the scaling and wavelet function at the lowest level, i.e.

u 1 (x) = k c 0,k φ(x -k) + k d 0,k ψ(x -k) (46)
whose solution consists now of the coefficients c 0,k and d 0,k .

Because of the multi-resolution property of the employed approximation, one expects that the coefficients d 0,k become higher at locations k where the lowest-level approximation u 0 (x) based on the use of the lowest level scaling functions is not able to represent the solution with sufficient accuracy. Thus, we could define a threshold value κ, allowing us to identify two sets S - 1 and S + 1 , composed respectively by integers k such that d 0,k ≤ κ and d 0,k > κ.

Now, the next approximation level u 2 (x) is defined by adding to the previous one u 1 (x) only the wavelet functions related to points in S + 1 and the ones having a natural neighbor in S + 1 . The resulting extended set is denoted by S+ 1 , and now the approximation u 2 (x) reads

u 2 (x) = u 1 (x) + k∈ S+ 1 d 1,k ψ 1,k (x) (47)
After calculating coefficients d 1,k , they are classified in sets S - 2 and S + 2 , and the adaptive approximation continues. Fig. 2 depicts the evolution of the approximation basis on the left (with the red line being the computed solution), where the points at the lowest level (blue ones) are related to the uniform lowest-level approximation (based on the use of the scaling functions). Then, the second level, red points, are the ones in S+ 1 , and so on. On the right (Fig. 2), the different contributions to the solution are depicted: u 1 (x) (blue curve) as well as the adaptive contributions (for example, the red curve represents u 2 (x)u 1 (x), whose support is defined by S+ 1 ). In the next section, the wavelet approximations will be considered for discretizing PDEs within the separated representation format characteristic of PGD methods.

Multi-scale Proper Generalized Decomposition

For the sake of simplicity, we consider the heat equation

∇ (K (x, y)∇u(x, y)) = f (x, y) (48)
where K (x, y) represents the material conductivity, assumed isotropic.

Problem (48) is defined in = x × y = (0, L) × (0, H) and where without loss of generality the homogeneous Dirichlet boundary condition is enforced on the domain boundary ≡ ∂ .

The weak form related to Eq. (48) reads

x × y ∇u * • (K ∇u(x, y) + f (x, y)) dx d y = 0 (49) whose solution, as described in Section 1, is searched under the separated form

u n (x, y) = n i=1 X i (x) • Y i (y) (50)
which implies the alternative resolution of two one-dimensional problems, the first one for calculating X i (x), and the second one for calculating Y i (y), similarly to Eqs. ( 16) and [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF].

When considering wavelet-based approximations, the functions X i (x) and Y i (y) at the enrichment iteration n and the nonlinear iteration p, involved in the one-dimensional problems of the PGD separated representation constructor, are approximated according to:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ X p n (x) = 2 j 0 /2 k X p,n j 0 ,k φ 2 j 0 x -k + J j= j 0 2 j/2 k X p,n j,k ψ 2 j x -k Y p n (y) = 2 j 0 /2 k Y p,n j 0 ,k φ 2 j 0 y -k + J j= j 0 2 j/2 k Y p,n j,k ψ 2 j y -k X * n (x) = 2 j 0 /2 k X * ,n j 0 ,k φ 2 j 0 x -k + J j= j 0 2 j/2 k X * ,n j,k ψ 2 j x -k Y * n (y) = 2 j 0 /2 k Y * ,n j 0 ,k φ 2 j 0 y -k + J j= j 0 2 j/2 k Y * ,n j,k ψ 2 j y -k X i (x) = 2 j 0 /2 k X i j 0 ,k φ 2 j 0 x -k + J j= j 0 2 j/2 k X i j,k ψ 2 j x -k Y i (y) = 2 j 0 /2 k Y i j 0 ,k φ 2 j 0 y -k + J j= j 0 2 j/2 k Y i j,k ψ 2 j y -k ( 51
)
where j 0 refers to the lowest level and where the coefficients affecting the different scaling and wavelet functions are the unknowns.

It is important to note that, even if the approximates (51), as written, involve all the translations k at each level j, the multi-resolution analysis allows considering, as previously discussed, only the wavelet functions located in the regions in which they contribute to the solution improvement. Moreover, even if the largest scale is assumed to be the same for all the functions involved in the separated representation, J for all them, in general each function at each level could involve a different number of scales.

As previously described, all these functions, as well as their derivatives, can be computed at dyadic points, until approaching sufficiently the considered point. Thus, by injecting the previous approximates into the weak form of the one-dimensional problems resulting from the PGD discretization of Eq. (49) and using an integration quadrature, the approximation coefficients can be obtained.

It is important to mention that Eq. (49) also involves the conductivity K (x, y), which should be separated. For this purpose, we proceed as described in our former works (e.g., [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF]) by constructing its separated representation

K (x, y) ≈ K i=1 K x i (x)K y i (y) (52)
The separation could be performed by using a standard SVD -Singular Value Decomposition. Another possibility consists in using the PGD (see Chapter 3 in [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations[END_REF]), which consists in solving the problem

K (x, y) -K (x, y) = 0 (53) with K (x, y) = K i=1 K x i (x)K y i (y) (54)
whose integral form reads

x × y K * (x, y) K (x, t) - K i=1 K x i (x)K y i (y) dx d y = 0 (55)
Again, this equation is solved by using the PGD constructor and now, because no regularity is needed in the representation of the conductivity, the Haar's wavelet is chosen, i.e. used for approximating functions involved in the conductivity separated representation K x i (x) and K y i (y).

The main issue when considering the Daubechies and Haar representations of the temperature field u(x, y) and the conductivity K (x, y) is that the integrals in Eq. ( 49) involve the product of two or three scaling and wavelet functions, of different nature (Daubechies and Haar) at different scales. This difficulty has been circumvented by considering the technique proposed in [START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF] and summarized in Appendix A.

After solving the problem, that is, after computing the coefficients affecting the different scaling and wavelet functions appearing in the approximation of functions involved in the solution separated representation, one realizes that each mode n ≤ N, X n (x) • Y n (y), involves different levels. The Galerkin PGD constructor does not allow associating modes with levels; each mode contains many levels of description, and each level of description is present in most of the PGD modes.

However, in order to extract the multi-scale features of the solution, after having solved the problem, one could proceed with re-ordering the modes by enforcing the first mode to contain the lowest-level contribution (zero level), the second mode the first level contribution, and so on.

To better describe the re-ordering procedure, we consider the solution approximation

u N (x, y) = N i=1 X i (x)Y i (y) (56) with ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ X i (x) = 2 j 0 /2 k X i j 0 ,k φ 2 j 0 x -k + J j= j 0 2 j/2 k X i j,k ψ 2 j x -k Y i (y) = 2 j 0 /2 k Y i j 0 ,k φ 2 j 0 y -k + J j= j 0 2 j/2 k Y i j,k ψ 2 j y -k (57)
Solution ( 56) can be rewritten in the form

u N (x, y) = J i=0 Z i (x, y) (58)
with

Z 0 (x, y) = N i=1 2 j 0 /2 k X i j 0 ,k φ 2 j 0 x -k • 2 j 0 /2 k Y i j 0 ,k φ 2 j 0 x -k (59) Z 1 (x, y) = N i=1 2 j 0 /2 k X i j 0 ,k φ 2 j 0 x -k • 2 j 0 /2 k Y i j 0 ,k ψ 2 j 0 y -k + 2 j 0 /2 k X i j 0 ,k ψ 2 j 0 x -k • 2 j 0 /2 k Y i j 0 ,k φ 2 j 0 y -k (60)
and so on. Here, each term Z i (x, y) contains a level of description, from the coarsest one to the finest one.

It is important to note that each term Z i (x, y) has a separated representation, a finite sum of functional products, with one of the involved functions depending on the x-coordinate, and the other one on the y-coordinate.

Numerical examples

The Poisson problem: convergence analysis

In order to check the proposed strategy, we consider the problem discussed in Section 1, Eq. ( 8), with f = 1 and homogeneous boundary conditions whose exact solution, as previously discussed, is available. Fig. 3 depicts the reconstructed solution and the associated modes in each direction. It was noticed that the approximation of functions involved in higher modes (e.g., X N (x) and Y N (y)) required higher levels. The difference between the reconstructed and the exact solution is depicted in Fig. 4. The convergence analysis considered increasing the number of models involved in the separated representation for different levels in the mode approximation and the inverse, increasing the number of levels for a fixed number of modes. The error is calculated using the L 2 -norm. As expected, a strategy increasing either the number of modes or the number of levels rapidly reaches a plateau, as Fig. 5 attests to. Only by increasing simultaneously both of them, the convergence rate can be maintained. Fig. 6 compares the first three PGD models X i (x) • Y i (y) and the re-ordered multi-scale modes, referred to as WPGD modes, Z i (x, y), according to the strategy described in the previous section.

Finally, Fig. 7 illustrates the contribution of different levels to each PGD mode. It can be noticed that the first PGD modes involve mainly lower levels, whereas higher levels contribute mainly to the last PGD modes.

Multi-scale analysis

In this case study, we consider a thermal problem (heat transfer equation) in the domain depicted in Fig. 8, consisting of two phases with different thermal conductivities (both of them assumed isotropic). As discussed before, Daubechies wavelets with N = 6 are considered for approximating the functions composing the different solution modes, i.e. the functions X i (x) and Y i (y), while the separated representation of the conductivity uses a Haar wavelet representation.

First, the conductivity is approximated at a fine enough level by using a Haar wavelet representation. Then, the solution is approximated at different levels using Daubechies wavelets. When considering fine scales, as fine as the microstructure representation, the solution becomes locally representative of the exact solution, whereas when considering coarser representations with respect to the microstructure length scale, the latter is considered in an averaged sense in the low-scale solution approximations. Fig. 9 represents the solution when considering different representation levels of the solution.

Finally, Figs. 10 and 11 depict respectively the PGD modes X i (x) • Y i (y) and their associated multi-scale re-ordered modes Z i (x, y) described in the previous section, where the multi-scale character is pointed out, with higher frequencies appreciated at higher modes.

Parametric solutions

Until now, separated representations were associated with the space coordinates x and y. However, as discussed in Section 1, one appealing feature of separated representations concerns the solution to parametric models. In [START_REF] Chinesta | Parametric PGD based computational vademecum for efficient design, optimization and control[END_REF], the solution to parametric models within the PGD framework was widely described. Thus, if one is interested in calculating the temperature in the domain for any value of the material conductivity K ∈ I ⊂ R + , the solution separated representation reads

X i (x) • Y i (y). u(x, K ) ≈ N i=1 X i (x) C i (K ) (61)
which, injected in the problem integral form, allows calculating the functions X i (x) and C i (K ). However, such a procedure is too intrusive to be considered for calculating a parametric solution when using an external simulation code able to provide the temperature field for a given value of the conductivity, that is, u(x; K ). The last approach is based on an appropriate sampling of the parametric domain by choosing particular values of the conductivity grouped in the set K = {K 1 , K 2 , • • • , K K } and solving the problem for all these values in order to obtain u(x; K i ), ∀i ∈ K. Then, all these solutions can be adequately interpolated for defining the solution for any conductivity value K ∈ [K 1 , K K ]. This is the key idea behind surrogate models, response surfaces, etc.

If we imagine for a while that we solved the problem at hand for the minimum and maximum values of the conductivities, leading to u(x, K min ) and u(x, K max ), the conductivity for any other value of the conductivity K ∈ [K min , K max ] can be linearly approximated according to

u(x, K ) = u(x, K min )N 1 (K ) + u(x, K max )N 2 (K ) (62)
where the approximation functions in the parametric space N 1 (K ) and N 2 (K ) read

N 1 (K ) = K -K min K max -K min (63) and N 2 (K ) = K max -K K max -K min (64)
Finer approximations require finer samplings and the subsequent interpolation. The main difficulty is how choosing these points to keep the accuracy under control. In [START_REF] Borzacchiello | Non-intrusive sparse subspace learning for parametrized problems[END_REF], the authors proposed a non-intrusive sparse subspace learning approach using hierarchical approximation bases.

Another possibility consists in taking advantage of the multi-resolution property of wavelet-based approximations for controlling the sampling richness, that is, the number of elements in the set K. Thus, one can control the convergence by analyzing the value of the associated coefficients related to the wavelet functions. This appealing property is absent in usual polynomial approximations. Using a wavelet representation of the parametric evolution employing spline-wavelets, φ (x) = has the advantage of recovering coefficients that correspond to the value of the approximated function at its location. Then, as soon as the solution is computed at two consecutive levels by approximating the parametric functions by using the scaling functions φ j,k and φ j+1,k , respectively, their difference represents the wavelet contribution that can be expressed in the wavelet basis ψ j,k , and whose higher coefficients indicate the locations at which sampling must be refined.

This procedure allowed us to reduce significantly the number of sampling points, of some orders of magnitude in some of our numerical experiments.

Conclusions

In this work, we proved that separated representations involved in PGD solution procedures can be combined with wavelet-based functional approximations for extracting the multi-scale behavior present in the solutions.

Moreover, we proved that the multi-resolution property inherent to wavelet representations allows defining simple adaptive procedures with multiple applications: model refinement or efficient sampling for calculating parametric solutions.

The connection between multi-scale solutions in complex microstructures and homogenization constitutes a work in progress.

Appendix A. Calculating integrals involving three functions at the same and at different levels

We consider the integral where d1, d2 and d3 refers to the derivative order. Integral (66) can be expressed as

d 1 ,d 2 ,d 3 l,m = 2 d 1 2 d 2 2 d 3 ∞ -∞ k a k φ d 1 (2x -k)
= 2 d 1 2 d 2 k a k r a r ∞ -∞ φ d 1 (z) φ d 2 (z + k -2m -r) φ d 3 (z + k -l) = 2 d 1 2 d 2 k a k r a r 2m+r-k,l-k ( 72
)
which results again in an eigenproblem.

pn

  (x), we are now ready to proceed with the second step of iteration p, which updates Y p n (y) from the just-computed X p n (x). The procedure exactly mirrors what we have done above. Indeed, we simply exchange the roles played by all relevant functions of x and y, starting from the solution approximate

Fig. 1 .

 1 Fig. 1. PGD error E M (u N M ) as a function of the number of enrichment steps N for different numbers M of discretization points.

  the orthogonality condition between φ(x) and its integer translates into∞ -∞ φ(x)φ(xl) dx = δ 0l (30)with δ the Kroenecker delta function. Taking into account Eq. (23), the previous equation results in k=∞ k=-∞ a k a k+2l = 2δ 0l , ∀l ∈ Z (31)

Fig. 2 .

 2 Fig. 2. Adaptive approximation strategy.

Fig. 3 .

 3 Fig. 3. Reconstructed solution and involved modes in each direction.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Difference between the reconstructed solution and the exact one for different number of models considered in the solution separated representation.

Fig. 6 .

 6 Fig. 6. PGD modes X i (x) • Y i (y) versus WPGD modes Z i (x, y).

Fig. 7 .

 7 Fig. 7. Contribution of the different levels to each PGD mode.

Fig. 8 .

 8 Fig. 8. Domain exhibiting two phases with different thermal conductivities.

Fig. 9 .

 9 Fig. 9. Reconstructed solution for different approximation levels.

Fig. 10 .

 10 Fig. 10. PGD modes: X i (x) • Y i (y).

Fig. 11 .

 11 Fig. 11. Multi-scale modes Z i (x, y) constructed from the wavelet-based PGD separated representation.

d 1 ,d 2 φ

 12 d1 (x) φ d2 (xl) φ d3 (xm) dx (66)

r a r φ d 2 φ

 2 (2x -2lr) p a p φ d 3 (2x -2mp) dx (67)By considering the change of variable y = 2xk, the preceding expression readsd 1 ,d 2 ,d 3 d 1 (y) φ d 2 (y + k -2lr) φ d 3 (y + k -2mp) d y (68)which can be rewritten asd 1 ,d 2 ,d 3 rank-deficient eigenproblem.The rank deficiency is circumvented by adding the required number of extra equations[START_REF] Latto | The evaluation of connection coefficients of compactly supported wavelets[END_REF]. Now we consider that the integral involves functions at different scales,d 1 ,d 2 ,d 3 d1 (y) φ d2 (ym) φ d3 (2yl) d y (70)which can be rewritten asd 1 ,d 2 ,d 3 d 1 φ d 1 (2yk) r a r 2 d 2 φ d 2 (2y -2mr) φ d 3 (2yl)


Cédric Ray
email: cedric.ray@univ-lyon1.fr

Mathieu Caillau

Christian Jonin

Emmanuel Benichou

Christophe Moulin

Estelle Salmon

Melissa E Maldonado

Anderson S L Gomes

Virginie Monnier

Emmanuelle Laurenceau

Jean-Louis Leclercq

Yann Chevolot

Thierry Delair

Pierre- François Brevet

Quadratic Nonlinear Optics to Assess the Morphology of Riboflavin doped Chitosan for Eco-friendly Lithography

Keywords: Nonlinear optics, second harmonic generation, green lithography, riboflavin, chitosan

We report the use of the Second Harmonic Generation response from a riboflavin doped chitosan film as a characterization method of the film morphology. This film is of particular interest in the development of new and bio-sourced material for eco-friendly UV lithography.

The method allows us to determine how riboflavin is distributed as a function of film depth in the sample. This possibility is of importance in order to have a better understanding of the riboflavin influence in chitosan films during the lithography process. On the contrary, linear optical techniques provide no information beyond the mere confirmation of the riboflavin presence.

Introduction

Photolithography can nowadays achieve designs with high 2D spatial resolution for a wealth of industrial applications, notably in the field of nanoelectronics. This possibility results from the technical developments realized over the years in nanofabrication in the decrease in wavelength or the development and optimization of materials (masks, resists) and devices (lens, wavelength sources, exposure tools …) for each technological node [START_REF] Ronse | Optical lithography-a historical perspective[END_REF].

Photolithography is based on the use of photoresists generally composed of a polymer, a solvent, a photoactive component as a photo-acid generator, a quencher and different additives. This polymer solution is coated on the substrate leading to a thin film then patterned through photolithography. After development, it forms a temporary mask that protects selected areas of the underlying substrate allowing for the micro/nanopatterning (by etching or depositing matter locally) of the substrates for the fabrication of the devices.

However, with a growing concern on environmental and health issues, new green photoresists are required. In this context, chitosan as a natural polymer, sub-product of the chitin shells of shrimps and other crustaceans, offers interesting possibilities. Its doping with riboflavin, a natural dye, seems promising.

Chitosan is a polysaccharide constituted of D-glucosamine and N-acetyl-D-glucosamine units linked by a β-(1→4) glycosidic bond (Figure 1). This compound is naturally present in insects and microorganisms but can be also industrially prepared from chitin by treatment with an alkaline solution [START_REF] Rinaudo | Chitin and chitosan: Properties and applications[END_REF]. Chitin is the second most abundant biopolymer on Earth. It is extracted from wastes of seafood industry like crustacean shells from crabs and shrimps in particular.

Chitosan is nontoxic, biocompatible and biodegradable At acidic pH, chitosan is soluble in aqueous solutions. Riboflavin or vitamin B2 (Figure 1) is a natural and hydrosoluble dye present in fruits, vegetables and animals [START_REF] Sheraz | Photo, thermal and chemical degradation of riboflavin[END_REF]. It is known for its photosensitivity and used as a photodegradation agent of organic pollutants in water treatment [START_REF] Escalada | Photosensitized degradation in water of the phenolic pesticides bromoxynil and dichlorophen in the presence of riboflavin, as a model of their natural photodecomposition in the environment[END_REF]. Riboflavin-doped photosensitive chitosan films can be used for an eco-friendly and water-developable photolithography.

In photolithography, the polymer solution is spin-coated on a substrate prior to illumination. Therefore, novel characterization methods need to be proposed to determine the physical properties of the polymer thin film in view of the targeted eco-friendly lithography. In particular, the morphology of the film must be determined since the riboflavin distribution 3/20 within the film is of utmost importance. Indeed, as the solvent is evaporated from the film, the solubility limit of riboflavin (around 2×10 -4 mol.L -1 ) may be reached leading to phase segregation whereas a homogeneous distribution is desired. Linear optical methods like UVvisible absorption and fluorescence spectroscopies are among the simplest methods and indeed provide valuable information. Chitosan, although presenting a weak fluorescence, has no prominent features in the visible part of the spectrum as opposed to riboflavin which possesses both strong absorption and fluorescence bands. Yet, linear optical methods are known to lack longitudinal resolution along the direction of propagation and are therefore of limited use when dopant depth profiles are sought at micrometer resolution or better.

In order to improve the longitudinal resolution of the optical characterization, we propose in this work to use different complementary nonlinear optical methods. First, we have applied the two photon excited fluorescence (TPEF) method which permits to obtain a threedimensional resolution owing to the multiphoton excitation process [START_REF] So | Two-Photon Excitation Fluorescence Microscopy[END_REF][START_REF] Denk | Two-photon laser scanning fluorescence microscopy[END_REF][START_REF] Xu | Measurement of two-photon excitation cross sections of molecular fluorophores with data from 690 to 1050 nm[END_REF]. Indeed, because two photons are simultaneously required to achieve material excitation, the nonlinear optical process can only occur at the focus of the beam [START_REF] Kaiser | Two-Photon Excitation in CaF2:Eu2+[END_REF]. Thus, with a typical focal volume of 10 µm 3 , a depth resolution of few micrometers in the plane perpendicular to the propagation direction can be achieved. However, the linear and nonlinear methods have a similar resolution limited by diffraction and are not so well suited to investigate the film interfaces.

Therefore, we propose to also perform Second Harmonic Generation (SHG) measurements at the air-film interface [START_REF] Shen | Surface properties probed by second-harmonic and sum-frequency generation[END_REF]. The SHG process whereby two photons at the fundamental frequency are converted into one photon at the harmonic frequency, is highly sensitive to the symmetry of the material and in particular, is forbidden in media with inversion symmetry within the dipole approximation, that is within the volume of the film. Hence, the method is inherently surface sensitive and should allow us to go further in assessing any presence or accumulation of riboflavin at the air-film interface. The surface sensitivity is here determined by the extent of the centrosymmetry breaking into the film. The interface is a region where the centrosymmetry is broken whereas the volume is centrosymmetric by nature. Similarly to liquids, this interface region where the centrosymmetry is broken should not exceed a few molecular layers, i.e. should be of the order of a nanometer. With the combination of TPEF and SHS, to determine the riboflavin distribution in the film volume, and SHG, to determine its presence or absence at the film interfaces, one can thus retrieve a full picture of the 4/20 riboflavin spatial distribution. In fact, after the first few atomic or molecular layers, any centrosymmetric media will exhibit no coherent SHG signal. Only the symmetry breaking due a media change at the surface, can give any signal for a centrosymmetric media.

From the molecular structure of riboflavin which contains delocalized electrons in an asymmetric environment (Figure 1), we can expect that the molecule exhibits a significant first hyperpolarizability, namely its cross-section for SHG, and thus a significant second harmonic intensity. The first hyperpolarizability of a molecular compound can be determined by the Second Harmonic Scattering (SHS) method [START_REF] Clays | Hyper-Rayleigh scattering in solution[END_REF][START_REF] Dadap | Second-Harmonic Rayleigh Scattering from a Sphere of Centrosymmetric Material[END_REF][START_REF] Revillod | Hyper Rayleigh Scattering of Molecular Aggregates[END_REF], also referred as incoherent SHG. In this case, the harmonic light results from the scattering by a molecular suspension in a liquid due to orientation fluctuations. In this paper, we therefore report the SHG response from a riboflavin doped chitosan film deposited on a silicon substrate supporting a 200 nm silica layer. We first determined the film linear optical properties and the riboflavin first hyperpolarizability when dispersed in bulk solution using the SHS technique. Then, the TPEF and SHG intensities of riboflavin in thin chitosan films were recorded. The 2D spatial distribution of the molecules in the film was characterized by nonlinear imaging. Finally, surface SHG measurements were performed in order to confirm the absence of riboflavin at the air-film interface.

Experimental Section Chemistry

Riboflavin was purchased from Sigma-Aldrich (> 98 %). Chitosan was purchased from Mahtani Chitosan PVT, Ltd (India) and was characterized by its molecular weight and its degree of acetylation (DA) corresponding to the number of N-acetyl-D-glucosamine units along a chain of chitosan. In our case, chitosan was extracted from squid pens with a DA of 2 % and a molar mass of 569.9 kg.mol -1 .

A riboflavin solution concentrated at 2.2×10 -4 mol.L -1 in deionized water (ultra-pure, 18 MΩ.cm) was prepared. Chitosan was dissolved in this solution at a concentration of 0.9 % (w/v). To enable its total dissolution, acetic acid (Fluka Analytical, glacial acetic acid) was added at a final concentration of 0.05 mol.L -1 . Riboflavin doped chitosan films were realized by spin-coating (5000 rd.min -1 , 3000 rd.s -2 during 30 s) on n-doped (100) oriented Si substrates with a 200 nm thick SiO2 layer deposited by plasma enhanced chemical vapor 5/20 deposition (PECVD). After spin-coating, the resulting film was heated at 100 °C for 1 min. The riboflavin-doped chitosan film thickness was 180 nm ± 5 nm. Film thickness was measured by ellipsometry with a spectroscopic ellipsometer UVISEL ™ from Horiba using the classical formula dispersion as modelled in the Horiba software. The model was adapted to chitosan by adjusting the refractive index between 1.5-1.6 and the extinction coefficient to 0 at 630 nm following Nosal et al. [START_REF] Nosal | UV-vis-infrared optical and AFM study of spin-cast chitosan films[END_REF].

Exposure

Films were exposed to UV radiations through a photomask constituted of open zones (lines of 60 µm) in quartz and dark zones (squares of 440 µm x 440 µm) in chromium-iron oxide.

Exposure was realized during 1 minute on a MBJB4 UV SUSS Micro Tec UV lamp at 365 nm, 405 nm and 435 nm with a power density of 20 mW.cm -2 at 405 nm and a lamp power of 200 W. Patterns were revealed in the riboflavin-doped chitosan films after development in deionized water for 60 s. Non exposed films were stored in clean room in the dark after preparation [START_REF] Caillau | Sub-micron lines patterning into silica using water developable chitosan bioresist films for eco-friendly positive tone e-beam and UV lithography[END_REF].

Linear Optics

Ultraviolet-Visible absorption spectra were collected in a quartz cell with 1 cm path length using a double beam UV mc2 SAFAS spectrometer. One-photon excited fluorescence spectra were performed with a FLS920 fluorimeter from Edinburgh Photonics using quartz cells with 1 cm path length. A 450 W continuous xenon arc lamp was used as an excitation source.

Nonlinear optics

The nonlinear imaging experimental setup was based on a frequency doubled femtosecond Er-doped fiber laser source (Menlo Lasers, C-Fiber 780). The laser provided pulses with a duration of about 100 fs at a repetition rate of 100 MHz. The fundamental wavelength was set at 780 nm (spectral width of 13 nm) and an average power of about 65 mW was measured at the laser exit. The light beam was focused by a standard microscope objective (magnification X20, NA 0.75) mounted onto a microscope stand (Nikon Eclipse TE2000-U with 3D positioning in inverted geometry) holding the sample. Galvanometric mirrors allowed for the beam to be swept to acquire non-linear images of various samples. The nonlinear optical signals, SHG and TPEF, were collected into an optical fiber feeding a spectrometer after going through various color filters. These filters were placed on the emission light beam path to 6/20 avoid any unwanted fundamental or SHG light and also to separate TPEF and SHG signals. A cooled photomultiplier tube was placed after the spectrometer.

Hyperpolarisability measurements were performed in a SHS setup as described in Duboisset et al [START_REF] Duboisset | A bottom-up approach to build the hyperpolarizability of peptides and proteins from their amino acids[END_REF]. Briefly, incoherent second harmonic emission is generated in a solution illuminated by a femtosecond laser focused by a microscope objective and measured, after collection by a microscope objective, on a photomultiplier installed on a spectrometer.

The detailed description of the SHG setup can be found elsewhere [START_REF] Martin-Gassin | Compression Induced Chirality in Dense Molecular Films at the Air-Water Interface Probed by Second Harmonic Generation[END_REF]. Briefly, in our configuration, a fundamental 780 nm incoming wave was impinging onto the film surface with a 70° angle with controlled linear polarization. The harmonic outgoing wave at 390 nm was then collected in specular reflection through a spectrometer and a CCD camera [START_REF] Martin-Gassin | Compression Induced Chirality in Dense Molecular Films at the Air-Water Interface Probed by Second Harmonic Generation[END_REF].

Results and Discussion

Riboflavin properties in solution

In order to investigate the properties of Riboflavin-doped Chitosan films we first discuss the optical properties of the Chitosan/Riboflavin solutions used to produce the films. Within the UV-visible spectral range, the optical features in absorption and fluorescence are essentially provided by riboflavin with a secondary contribution from chitosan. The UV-visible absorption spectrum of the initial riboflavin solution used in the film fabrication is shown in Figure 2. The absorbance of the 0.9% chitosan solution is negligible in the visible spectrum and is only present in the UV range. On the opposite, the 2.2×10 -4 mol.L -1 riboflavin solution exhibits four peaks, two in the visible and two in the UV range albeit at higher wavelengths as compared to chitosan.

Because the lithography process entails a UV light exposure of the film, a comparison is made with an exposed riboflavin + chitosan solution. Exposure was performed at 365 nm, 405 nm and 435 nm. The absorption spectrum of the exposed solution differs quite dramatically from the non-exposed riboflavin. In fact, its absorption is very similar to that of lumichrome, the reaction product of the riboflavin photodegradation [START_REF] Koziol | Absorption spectra of riboflavin, lumiflavin, and lumichrome in organic solvents[END_REF].

As shown in Figure 3, fluorescence spectra were also recorded under single photon excitation at 377 nm. Chitosan is known to present fluorescence [START_REF] Revillod | Multipolar Contributions to the Second Harmonic Response from Mixed DiA-SDS Molecular Aggregates[END_REF]. However, in our experiments, its fluorescence was very weak. The observed fluorescence spectrum therefore closely matches that of the Riboflavin and the mixture solutions. Since our chitosan solution yields a much 7/20 lower fluorescence signal [START_REF] Huang | Enhanced fluorescence of chitosan based on size change of micelles and application to directly selective detecting Fe3+ in humanserum[END_REF] compared to riboflavin solution, the riboflavin and chitosan mixture solution closely matches the riboflavin spectrum.

We estimated also the first hyperpolarizability of riboflavin. The SHG scattered intensity was collected in a right angle geometry for different dilutions of a riboflavin aqueous solution.

When normalized to the intensity of the neat solvent these samples exhibit a linear behavior and the slope p is directly given by:

𝑝 = 〈𝛽 2 〉 𝑁 𝑤 〈𝛽 𝑤 2 〉
where Nw=55.56 M is the number of water molecules per unit volume and 〈𝛽 𝑤 2 〉 and 〈𝛽 2 〉 are the square of the first hyperpolarizability, for water and riboflavin respectively, in the laboratory frame averaged over all space orientation in accordance with the liquid phase isotropy. The measurement was conducted at room temperature with a fundamental wavelength of 780 nm leading to the estimation of a value of √〈𝛽 2 〉 = (2.5±0.5)×10 -28 esu using the 〈𝛽 𝑤 2 〉 water reference given by Duboisset et al. [START_REF] Duboisset | A bottom-up approach to build the hyperpolarizability of peptides and proteins from their amino acids[END_REF]. This value, as expected, is significant if compared to 4-(4-dihexadecylaminostyryl)-N-methylpyridinium iodide, a compound well known for its large hyperpolarizability of (17.6±0.5)×10 -28 esu [START_REF] Revillod | Multipolar Contributions to the Second Harmonic Response from Mixed DiA-SDS Molecular Aggregates[END_REF]. Thus, the detection of riboflavin with SHG was readily accessible although the experimental conditions were not optimal.

Riboflavin-chitosan film properties

First, the reflectance properties of the film were determined. Figure 4 shows that the film linear optical properties derive from chitosan only since the riboflavin content is too low.

Note that in Figure 4, the reflectance spectra clearly exhibit the interference fringes arising from the thin film.

Then, the nonlinear optical properties of the film were studied with TPEF. Figure 5 shows the TPEF spectra recorded for a pure chitosan film as well as a mixed riboflavin-chitosan film. Two photons, with a wavelength of 780 nm, were used and absorbed. This excitation can be compared to an irradiation at 390 nm (half the initial wavelength) which is close to the excitation used to acquire the one photon fluorescence spectra (377 nm). The fluorescence maximum was measured at a wavelength of about 500 nm with a shoulder at shorter wavelengths, around 450 nm.

8/20

On the opposite, the pure chitosan film showed almost no two-photon excited fluorescence signal. Interestingly, both films presented a SHG response located at half the 780 nm excitation wavelength. However, the SHG intensity seemed unaffected by the riboflavin presence and seems to therefore arise mainly from the pure chitosan film on the SiO2 substrate.

In order to characterize the distribution of riboflavin within the chitosan film, depth profiles

were recorded for different wavelengths on the 2D pattern imprinted on the riboflavin doped chitosan film. Two distinct areas of the sample were selected to analyze the exposed and non-exposed regions, respectively identified with a red or a blue square on Figure 6 (right panel).

It is clear from Figure 6 that the non-exposed region provided less fluorescence than the exposed one.

In Figure 6, the neat SiO2 substrate, designated as 'reference' and shown in black, is obtained in absence of riboflavin. The shifted position is due to the uncertainty when positioning different samples. This depth profile measurement resulted from the convolution of the beam focus and the film depth. Hence a TPEF maximum indicated that at the resolution imposed by the beam waist (few micrometers) riboflavin is distributed rather homogeneously within the film.

Figure 7 provides the depth profile recorded for SHS at the same areas as in Figure 6. The results are very similar to those obtained with TPEF and confirm that riboflavin was rather homogeneously distributed in the film. However, the larger contrast between the exposed and non-exposed areas obtained with SHG suggests that this method provides better contrasted images of the sample after UV exposure. The presence of riboflavin clusters with a size of tens of microns is also visible on the right panel of figure 7 (bright white signature).

This feature probably arises for undissolved riboflavin crystals although riboflavin reaches its maximum of solubility under the conditions used. The solution was however not filtered.

Despite the rather homogeneous distribution of riboflavin molecules in the bulk of the film, it was impossible to ensure the absence of riboflavin accumulation at the surface of the film (first few molecular layers) with the profiles obtained. This feature can be accessed with surface SHG known to provide a molecular information on the first few molecular layers of 9/20 centrosymmetric films like chitosan [START_REF] Shen | Optical Second Harmonic Generation at Interfaces[END_REF][START_REF]Liquid Interfaces Probed by Second-Harmonic and Sum-Frequency Spectroscopy[END_REF][START_REF] Volkov | Liquid-Liquid Interfaces Theory and Methods[END_REF]. We emphasize here that for a centrosymmetric film, similarly to liquids, it is assumed that this breaking of centrosymmetry pertains over a few molecular diameters only, i.e. over distances of the order of 1 nm. Any breaking over much larger distances, like micrometers, would be observed with the volume profile. Indeed, the SHG signal is only coming from the few molecular layers where the surface is not centrosymmetric. Moreover, the laser interaction with the bulk of the centrosymmetric film does not give any SHG signal.

We performed a polarization analysis of the surface SHG intensity collected from the air-film interface in presence and absence of riboflavin where the input fundamental beam angle of polarization was rotated for two perpendicular angles of the analyzer, respectively S and P.

Hence, two plots were recorded, the SHG intensities so-called IP and IS, for chitosan films and chitosan/riboflavin films under non-exposed conditions (Figure 8). Data obtained for other film compositions show very similar trends and are not shown. The plots all exhibit a two-lobe pattern for the P-polarized graph (squares) and a four-lobe pattern for the S-polarized graph (circles) of the SH intensity as a function of the input polarization angle, as expected from theory. These plots were analyzed using the standard simplified form of the SHG intensity as a function of the input polarization angle 𝛾 [START_REF] Bruyère | Reversibility of the supramolecular chirality of bridged binaphtol derivatives at the airwater interface[END_REF]:

𝐼 𝑃 ∝ |𝑎 ⋅ cos 2 (𝛾) + 𝑏 ⋅ sin 2 (𝛾)| 2 (2a) 𝐼 𝑆 ∝ |𝑐 ⋅ sin(2𝛾)| 2 (2b)
where a, b and c are three parameters depending on the geometrical configuration, the optical indices at 800 and 400 nm and the three non-vanishing independent tensor elements of the isotropic air-film interface, namely 𝜒 𝑆,𝑋𝑋𝑍 (2) , 𝜒 𝑆,𝑍𝑋𝑋

and 𝜒 𝑆,𝑍𝑍𝑍 (2) . The adjustment procedure for all different films composition led to the numerical values reported in Table 1. The recorded SHG intensity was dominated by the 𝐼 𝑃 component whereas the 𝐼 𝑆 (𝛾)-component was much weaker as shown on Figure 8. These results were identical for samples after and before exposition as well as for with pure chitosan films (left of Figure 8). Once the parameters are renormalized with respect to 𝑐, it is obvious that the parameters 𝑎 and 𝑏 are rather constant indicating that riboflavin does not accumulate at the surface. These results clearly demonstrate the absence of riboflavin in the first few molecular layers of the sample.

It is also clear that one observes a large increase of the 𝑎 coefficient upon the chitosan deposition on the surface of the neat substrate, as expected. Indeed the 𝑎 coefficient is proportional to 𝐼 𝑃 , that is the P polarized intensity at both fundamental and harmonic frequency. This large increase is the signature of drastic changes at the surface when chitosan is deposited on the Si/SiO2 bare substrate. This polarization configuration is the most sensitive one to these changes. In order to control further that the presence of riboflavin at the surface is independent of the substrate, we also analyzed the same films deposited on a glass substrate.

Conclusion and perspectives

We have reported the use of the nonlinear optical techniques of TPEF and SHG to investigate the distribution of riboflavin in a chitosan thin film deposited on a substrate. The results show that riboflavin does not accumulate at the air-film interface. Moreover, the SHG technique demonstrate a better characterization of the chitosan/riboflavin film interface than the TPEF technique. The nonlinear optical techniques used in this study are therefore highly complementary and efficient in determining the riboflavin depth distribution of the doped chitosan films. They should thus help improving eco-friendly lithography processes (especially UV exposure photo-mechanism in presence of dyes) by a better polymer film characterization.

As we have mentioned, some samples contained riboflavin clusters. Thus it would be interesting to decrease the riboflavin concentration in the chitosan solutions to avoid precipitation. The first lead would be to precisely determine the riboflavin solubility in chitosan in order to decrease its concentration while preserving the lithographic properties.

Another important parameter to investigate is the degree of acetylation (DA) of the chitosan.

Indeed, we used low DA (2%) chitosan but higher DA ones (up to 50%) may affect the riboflavin film solubilization through a synergic interaction between the acetamide group and riboflavin. 16/20 17/20 18/20 19/20

Glass

20/20

Graphical Abstract

Figure 1 :

 1 Figure 1 : Chemical structures of chitosan and riboflavin.

Figure 2 :

 2 Figure2: Absorbance spectra of the UV-exposed and non-exposed mixture of the 2.2×10 -4 mol.L -1 riboflavin and 0.9 % chitosan solutions compared to non-exposed solutions of 2.2×10 -4 mol.L -1 riboflavin and 0.9 % chitosan solution.

Figure 3 :

 3 Figure 3: One Photon Excitation Fluorescence (OPEF) spectra of pure riboflavin, pure chitosan and riboflavinchitosan mixture solutions excited at 377 nm.

Figure 4 :

 4 Figure 4: Reflectance of the silica substrate and the silica substrate coated with chitosan and a mixed riboflavinchitosan film.

Figure 5 :

 5 Figure 5 : TPEF spectra for a pure chitosan film and a mixed riboflavin-chitosan film with a 780 nm wavelength excitation. An SHG peak is observed in both samples at 390 nm.

Figure 6 :

 6 Figure 6 : (a:left) TPEF signal integrated over the 410 -650 nm spectral region as a function of the film depth for two different areas of the sample and the reference (SiO2 substrate). (b:right) Image of the exposed sample as obtained at the depth giving the maximum TPEF intensity (z=0 µm).

Figure 7 :

 7 Figure 7: (a:left) SHS (second harmonic scattering) recorded at 390 nm as a function of the film depth for the same areas of the previous figure. (b:right) Image of the exposed sample obtained at the depth giving the maximum intensity of SHS (z=0 µm).

Figure 8 :

 8 Figure 8: Surface SHG intensity as a function of the fundamental polarization angle for a chitosan film on Si/SiO2 substrate (a:left) and a non-exposed riboflavin-doped chitosan film on Si/SiO2 substrate (b:right). The black squares and red circles curves correspond respectively to the P-and S-harmonic intensities (IP and IS) with their corresponding fits.

Table 1 :

 1 Adjustment parameters and ratios normalized to the c coefficient, for the polarization resolved SHG intensity plots obtained using Eqs. (2a), (2b) for films on Si/SiO2 substrate.

	Si/SiO2 substrate	a	a/c	b/c
	Substrate	183.6±0.6	2.05±0.01	-2.03 ± 0.01
	Substrate + Chitosan	224.6±0.3	5.81±0.05	-1.14 ± 0.03
	Substrate + Chitosan + Riboflavin	212.9±0.2	5.90±0.05	-1.18 ± 0.02
	Substrate + Chitosan + Riboflavin (insolated) 221.1±0.4	5.62±0.04	-1.09 ± 0.04

Table 2 :

 2 Adjustment parameters and ratios normalized to the c coefficient, for the polarization resolved SHG intensity plots obtained using Eqs. (2a), (2b) for films on glass substrate.

Table 2 .

 2 As for the Si/SiO2 substrate, the 𝑎 coefficient exhibits an increase upon the deposition of a chitosan film, though the relative increase is much smaller. After normalization by the 𝑐 coefficient, one can note that the results are similar to the Si/SiO2 substrate indicating no clear sign of riboflavin accumulation within the first few layers of the film interface. These results also illustrate the great sensitivity of surface SHG[START_REF] Bruyère | Reversibility of the supramolecular chirality of bridged binaphtol derivatives at the airwater interface[END_REF] to probe these film interfaces.

	The fit parameters

Acknowledgements

This work is supported by the LABEX iMUST (ANR-10-LABX-0064) of Université de Lyon within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). This work is performed thanks to the Centre of Nano-Optics Nanoptec of the University Claude Bernard Lyon 1 and the France-Brazil CAPES-COFECUB (Sv852-15) program is also acknowledged.

NanoLyon platform is acknowledged for equipment and support. Green lithography project is supported by Pulsalys. MC also acknowledges the EEA doctoral school for a doctoral scholarship.