Luv Sharma

Ron H J Peerlings
email: r.h.j.peerlings@tue.nl

Franz Roters

Marc G D Geers

FFT-based interface decohesion modelling by a nonlocal interphase

Keywords: FFT-based spectral methods, Interfaces, Damage, Nonlocality

In this paper, two nonlocal approaches to incorporate interface damage in fast Fourier transform (FFT) based spectral methods are analysed. In FFT based methods, the discretisation is generally non-conforming to the interfaces and hence interface elements cannot be used. This limitation is remedied using the interfacial band concept, i.e., an interphase region of a finite thickness is used to capture the response of a physical sharp interface. Mesh dependency due to localisation in the softening interphase is avoided by applying established regularisation strategies, integral based nonlocal averaging or gradient based nonlocal damage, which render the interphase nonlocal. Application of these regularisation techniques within the interphase sub-domain in a one dimensional FFT framework is explored. The effectiveness of both approaches in terms of capturing the physical fracture energy, computational aspects and ease of implementation is evaluated. The integral model is found to give more regularised solutions and thus a better approximation of the fracture energy.

Introduction

Polycrystalline materials at a microscopic level show clear heterogeneous deformation patterns. This heterogeneity arises from the locally fluctuating mechanical properties of different phases and differences in lattice orientations between different grains.

FFT based spectral solver was originally introduced to model the mechanical behaviour of composite microstructures [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. Since then it has emerged as a promising tool for modelling the micromechanical response of polycrystalline materials [START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF][START_REF] Eisenlohr | A spectral method solution to crystal elasto-viscoplasticity at finite strains[END_REF]. A comparison of different FFT formulations and solution approaches in a crystal plasticity constitutive framework [START_REF] Roters | Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: theory, experiments, applications[END_REF] was presented in [START_REF] Shanthraj | Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials[END_REF]. Recently, an FEM perspective on an FFT based spectral formulation for small strain non-linear material behaviour was given in [START_REF] Zeman | A finite element perspective on nonlinear FFT-based micromechanical simulations[END_REF] and its extension to a finite strain setting was presented in [START_REF] De Geus | Finite strain FFT-based non-linear solvers made simple[END_REF]. Alongside such improvements, much effort has gone into making the method suitable for various applications. The computational efficiency of FFT methods makes them attractive to solve multi-field problems, for e.g. a nonlocal crystal plasticity formulation [START_REF] Lebensohn | Numerical implementation of non-local polycrystal plasticity using fast fourier transforms[END_REF], ferroelectric switching [START_REF] Vidyasagar | Predicting the effective response of bulk polycrystalline ferroelectric ceramics via improved spectral phase field methods[END_REF], etc.

Polycrystalline microstructures show a combination of intergranular and intragranular damage. Gradient based nonlocal damage simulations in a polycrystalline material were performed using an FFT approach in [START_REF] Boeff | Formulation of nonlocal damage models based on spectral methods for application to complex microstructures[END_REF] without differentiating between bulk damage and interface decohesion. Interfaces like phase boundaries and grain boundaries are particularly susceptible to damage initiation and propagation. For example, damage nucleation in Dual Phase steel has been attributed to interface decohesion at ferrite-martensite phase boundaries and martensite cracking along prior austenite grain boundaries [START_REF] Archie | Micro-damage initiation in ferrite-martensite DP microstructures: a statistical characterization of crystallographic and chemical parameters[END_REF]. In our contribution, we focus on the computational treatment of interface decohesion and its interaction with the (undamaged) bulk behaviour in FFT solvers.

In finite element (FE) based solution schemes, such problems are generally tackled using cohesive zones [START_REF] Xu | Void nucleation by inclusion debonding in a crystal matrix[END_REF][START_REF] Remmers | The cohesive band model: a cohesive surface formulation with stress triaxiality[END_REF], see Fig. 1 (a 1 , a 2 , a 3 ). Interfaces are identified a priori and interface elements are introduced at these physical locations. This is possible since the FE discretisation is made to conform to the sharp interfaces. The opening behaviour of these interface elements is governed by a traction-separation law (TSL) that is assigned to them. FFT based spectral methods, however, rely on a regular grid and a method to include sub-dimensional (e.g. planar in 3D) elements cannot be used. Composite pixel based approaches, for example [START_REF] Merkert | Voxel-based fast solution of the Lippmann-Schwinger equation with smooth material interfaces[END_REF], which make use of the interface normal to assign homogenised mechanical properties to grid points near the interfaces, are a step in that direction. But at this point it is not clear how they will perform in situations of damaging interfaces that are inclined or curved. The high mechanical contrast combined with inherent pixelation effects of FFT may cause stresses relayed across damaging interfaces, which is undesirable for the modelling of cracks. This limitation is dealt with using the idea of an interphase (volumetric) band, as depicted in Fig. 1 (b 1 , b 2 , b 3 ). The material points in the vicinity of the sharp interface are identified and furnished with a damage constitutive behaviour (softening) and corresponding kinematics-in addition to the deformation mechanisms attributed to them as a part of the respective grains they belong to. We aim to have multiple points across the thickness of the interphase. We expect this to reduce the pixelation effect and ensure that the two crack faces are fully decoupled. This approach allows capturing the interfacial mechanics and still benefits from the computational (memory and speed) and implementation related advantages of FFT based methods. Such an approach to represent interfaces as interphase has been used previously in the literature. Hsueh-Hung et al. [START_REF] Fu | Computational description of nanocrystalline deformation based on crystal plasticity[END_REF] used it to assign different material behaviour to grain boundary regions than the bulk. Their approach was physically motivated and directed towards understanding metal plasticity in the nanocrystalline regime. Clayton et al. [START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF] used this approach within a finite element implementation of phase field damage to understand the competition between intergranular and intragranular damage.

In volumetric dissipation based models, one specifies a dissipation density and an internal length scale to the model-which for the present case is the band thickness, l. The total dissipation resulting from the volumetric damage process must equal the physical dissipation of the real sharp interface (or cohesive zone model). Accordingly, the dissipation density has to scale inversely with the band thickness. A straightforward inversely proportional relationship results if the entire interphase band damages uniformly. This is not guaranteed in a conventional local softening model due to localisation and lacking objectivity with respect to the grid spacing. Various regularisation strategies have been used in the literature to remedy this. They can broadly be classified into two categories-integral based averaging [START_REF] Bažant | Nonlocal continuum damage, localization instability and convergence[END_REF] and gradient damage based regularisation [START_REF] De Borst | Some future directions in computational failure mechanics[END_REF][START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF][START_REF] Geers | Finite strain logarithmic hyperelasto-plasticity with softening: a strongly non-local implicit gradient framework[END_REF][START_REF] Shanthraj | Elasto-viscoplastic phase field modelling of anisotropic cleavage fracture[END_REF][START_REF] Moës | Coupling local and non-local damage evolutions with the thick level set model[END_REF].

The objective of this paper is to analyse the applicability of these two regularisation approaches in an FFT method to model the softening interphases. We compare the effectiveness of the two approaches for the idealized case of one dimensional (1D) bar with a single interphase band. The constitutive response of the underlying (bulk) material is assumed to be linear elastic. This simple case allows us to systematically study the effect that the nonlocality has on the damage evolution-with emphasis on the desired objectivity of the predicted fracture energy with respect to the interphase band thickness. The organisation of this paper is as follows. In "Problem statement" section the 1D problem considered in our study is stated and the constitutive model and kinematics are described. Next, the gradient damage and integral averaging approaches are discussed in "Gradient based nonlocal damage" and "Integral based nonlocal damage" sections, respectively.

0 x 1 = -l 2 x 2 = l 2 -L 2 L 2 x Sharp Interface

Interphase band

Fig. 2 1D bar with an interphase band (bottom) representation of a central sharp interface (top)

In "Results and discussion" section, results of both the approaches are discussed. A conclusion and outlook is presented in "Conclusion and outlook" section.

Problem statement

Consider a 1D bar of length L as shown in Fig. 2. Let us assume a sharp interface at the center (x = 0) of the bar. To capture the response of this interface using an FFT based spectral scheme, the sharp interface is substituted by a volumetric interphase band of thickness l spanning the region from

x 1 = -l 2 to x 2 = l 2 .
The elastic properties are taken to be homogeneous, with Young's modulus defined as E. To model the kinematics of decohesion, a damage strain field, ε d (x), is introduced for x ∈ Ω i where Ω i = (x 1 , x 2 ) is the sub-domain representing the interphase region in the total domain Ω = (-L/2, L/2). Since only the interphase region accommodates the damage strain, the damage constitutive model is only used in Ω i . In Ω i it is superimposed on the elastic strain ε e (which exists on Ω) as,

ε = ε e + ε d .
(

A local damage field φ l (x), for x ∈ Ω i , is introduced. This field represents the material integrity. It is initialised to unity, representing a perfectly bonded state. When the damage evolves it drops towards zero, thereby reducing the capacity of a material point to transmit stresses. A damage constitutive model exploiting only the local field φ l is ill-posed and provides non-physical results due to a pathological tendency to localisation [START_REF] Bažant | Nonlocal continuum damage, localization instability and convergence[END_REF]. Therefore, a nonlocal counterpart field, φ nl , is defined also for x ∈ Ω i . The bar is subjected to an overall monotonic uniaxial tensile loading at a constant rate ε. This loading rate is related to the local strain rate field as

ε = 1 L + L 2 -L 2 ε dx, (2)
The stress at equilibrium has to satisfy,

∂ σ ∂ x = 0. ( 3
)

Constitutive model

The constitutive model relating the stress and the elastic strain for all material points at

x ∈ Ω reads σ = Eε e . ( 4
)
The evolution law for ε d is taken as

εd = εd0 |σ | σ c φ nl 2 -1 2 sgn (σ ) , (5)
where σ c φ nl 2 is the effective opening resistance, with σ c the critical stress of the model; εd0 is the reference opening rate. The Macaulay brackets, defined as

• = (• + |•|) /2
, make sure that the damage strain does not evolve until σ reaches the initial critical stress σ c . The nonlocal damage φ nl , in Eq. ( 5) is obtained by appropriately regularizing the local damage field φ l defined as

φ l = 1 - ε d ε f . ( 6
)
The critical damage strain ε f governs the softening slope. It is noted that the rate dependent form depicted in Eq. ( 5) is motivated from the sub-scale mechanics (micro-void formation, plasticity etc.).

Analytical solution for the rate-independent case

Since the applied loading is monotonic tension, one can rewrite Eq. ( 5) beyond the damage threshold as

σ = σ c φ nl 2 1 + εd εd0 1 2 . ( 7
)
In order to retrieve an analytical solution, a special case is considered. It is assumed that the nonlocal field coincides with the local field, which is furthermore assumed to be uniform within the interphase band. In this case we have, also using Eq. ( 6),

σ = σ c 1 - ε d ε f 2 1 + εd εd0 1 2 . ( 8
)
In the rate-independent limit, implying εd << εd0 , (

Equation ( 8) yields

σ = σ c 1 - ε d ε f 2 . ( 10
)
Due to the equilibrium condition in this 1D setting, a spatially constant stress field results. The average strain ε for any instant after the damage has initiated can be obtained by inverting Eq. ( 10) to give ε d as a function of σ and integrating the elastic and damage contributions to the strain field:

ε = 1 L L 0 ε e dx + x 2 x 1 ε d dx = σ E + 1 - σ σ c 1 2 ε f l L . ( 11
)
Thus, the overall stress-strain response for this rate-independent uniform damage case can be written in piecewise form as

ε = ⎧ ⎪ ⎨ ⎪ ⎩ σ E , if ε ≤ σ c E σ E + 1 -σ σ c 1 2 ε f l L , otherwise. (12)
The dissipation per unit volume due to the damaging process can be computed analytically as

g c = ε f 0 σ dε d = σ c ε f 0 1 - ε d ε f 2 dε d = 1 3 σ c ε f , ( 13
)
where use has been made of Eq. ( 10).

Relating the dissipation throughout the interphase volume to the expected interfacial dissipation G c , we obtain

g c l = G c . ( 14
)
Accordingly, the product of the damage constitutive parameters σ c and ε f should depend inversely on the interphase band thickness. This equality, together with Eq. ( 13), indicates how the strain to failure needs to be scaled with respect to the interphase band thickness l, in order to recover the strength and dissipation of a cohesive zone, for arbitrary l:

ε f = 3G c σ c l . ( 15
)

Gradient based nonlocal damage

In the differential form of nonlocality, the nonlocal damage field is determined as the solution of a differential equation. Here, the classical gradient based damage equation is considered to retrieve the nonlocal counterpart of φ l , which reads in one dimension:

∂ ∂ x λ 2 (x) ∂ ∂ x φ nl + φ l -φ nl = 0, (16)
where λ(x) is the nonlocal length scale. In principle, the damage variables φ nl and φ l are meaningful only for x ∈ Ω i , and therefore Eq. ( 16) for the nonlocal field φ nl should be solved on the sub-domain Ω i only, with flux free boundary conditions (also referred as interface condition) on ∂Ω i , i.e.

λ 2 ∂φ nl ∂ x = 0 on x ∈ ∂Ω i , ( 17
)
where ∂Ω i = {x 1 , x 2 }. Since FFT solvers are based on global and periodic shape functions with a regular discretisation grid, the problem has to be solved on the full domain instead, still respecting Eq. ( 17) at least approximately. This may be achieved by extending Eq. ( 16) to the entire domain Ω, but with a piecewise constant λ, equal to the desired value λ in inside the band and a much smaller value λ out outside the band. In that case, the interface condition at the edges of the interphase reads:

λ 2 ∂φ nl ∂ x x - = λ 2 ∂φ nl ∂ x x + for x ∈ ∂Ω i , ( 18
)
which corresponds with Eq. ( 16). In the limit of λ in /λ out → ∞ this condition enforces Eq. [START_REF] Bažant | Nonlocal continuum damage, localization instability and convergence[END_REF]. In practice, a finite, but large, contrast λ in /λ out is used and Eq. ( 17), the flux free boundary condition, is satisfied only approximately. Phase field approaches also provide a differential equation based continuum solution for evolution of damage. In principle, a phase field model based on a Ginzburg-Landau approach can be used as well, for example [START_REF] Shanthraj | Elasto-viscoplastic phase field modelling of anisotropic cleavage fracture[END_REF][START_REF] Shanthraj | A phase field model for damage in elasto-viscoplastic materials[END_REF]. The resulting governing equations are similar to models based on the micromorphic approach, with an evolving length scale [START_REF] Poh | Localizing gradient damage model with decreasing interactions[END_REF]. In fact, these approaches were introduced to avoid spurious spreading of the damage field (away from the nucleation point), which was observed when using a linear phase field potential in the former and a decreasing length scale in the latter.

In the present approach, based on the classical gradient based damage, the issue of compact support of damage field is avoided by using the length scale contrast to restrict the damage field to the interphase band. In order to recover physically meaningful results, which are objective with respect to the thickness of the interphase band, the interphase band needs to damage uniformly, even for bands that are not extremely thin. From this perspective, the classical gradient damage model provides better regularisation characteristics than these phase field approaches. If a phase field approach would be considered instead, models based on quadratic or double-well potential would suit best for the good regularisation characteristics needed in our application.

The numerical solution of the boundary value problem requires solving Eqs. ( 3) and ( 16) simultaneously. Note that these equations are coupled through the constitutive Eq. ( 5). Both differential equations are discretised using the FFT scheme. The coupled system of equations is solved using a staggered iterative scheme as detailed in an FEM context in [START_REF] Shanthraj | Elasto-viscoplastic phase field modelling of anisotropic cleavage fracture[END_REF]. The discretisation and residual evaluation of the mechanical equilibrium differential equation is well documented in [START_REF] Shanthraj | Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials[END_REF]. Here, we only outline the residual formulation of the damage Eq. [START_REF] Clayton | Phase field modeling of directional fracture in anisotropic polycrystals[END_REF].

In order to avoid taking the Fourier transform of the length scale field, the non-linear residual is constructed by utilizing the fixed-point concept. We start by splitting λ 2 (x) into a homogeneous and a fluctuation field as

λ 2 (x) = λ2 + λ2 (x). ( 19
)
Substituting Eq. ( 19) in Eq. ( 16) and taking the Fourier transform (F ) we obtain

Fφ nl -λ2 F ∂ 2 φ nl ∂ x 2 = Fφ l + F ∂ ∂ x λ2 ∂ ∂ x φ nl , ( 20
)
which can also be written as,

Fφ nl -λ2 F ∂ 2 φ nl ∂ x 2 = Fφ l + F ∂ ∂ x λ2 F -1 F ∂ ∂ x φ nl . ( 21
)
Now using the property of the Fourier transform applied to differential operators, i.e, F ∂ n f ∂ x n = (ι k) n f , the above equation can be rewritten as

φnl -λ2 (ι k) 2 φnl = φl + (ι k)F λ2 F -1 (ι k) φnl . ( 22
)
After some rearrangements and taking the inverse Fourier transform (F -1 ), Eq. ( 22) is rewritten in residual form as

R d = φ nl -F -1 φl + (ι k)F λ2 F -1 (ι k) φnl 1 -λ2 (ι k) 2 = 0. ( 23
)
Initialising φ nl ≡ 1, Eq. ( 23) is solved using the matrix free Newton solver available in the PETSc library [START_REF] Balay | Petsc users manual revision 3.7[END_REF]. The frequencies used in the discrete Fourier transforms in this work correspond to the modifications presented in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF].

Integral based nonlocal damage

For the integral based nonlocal damage case, a weighted moving average is used to calculate the nonlocal field φ nl from the local damage φ l , i.e. [START_REF] Bažant | Nonlocal continuum damage, localization instability and convergence[END_REF][START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF].

φ nl (x) = Ω i Ψ (y; x) φ l (y) dΩ Ω i Ψ (y; x) dΩ , ( 24
)
where y is the position of the infinitesimal volume dΩ ∈ Ω i . The weighting function Ψ (y; x) that is commonly used is the Gaussian distribution, which reads

Ψ (y; x) = exp -(y -x) 2 2λ 2 , ( 25
)
where x is the position vector of the point at which the distribution is centred and λ is the characteristic length, which determines the distance along which Ψ decays to zero. The denominator in Eq. ( 24) normalizes the weighting function, which ensures that for a homogeneous φ l (x), the nonlocal field φ nl calculated using Eq. ( 24) coincides with its local counterpart.

In our numerical setting, Eq. ( 24) is approximated by a cell averaging method, which can be written as

φ nl i = j φ l j Ψ ((Δ x(j -i)) 2 ) j Ψ ((Δ x(j -i)) 2 ) , ( 26
)
where i and j represent material points in the interphase sub-domain and Δ x is the grid spacing. Equation ( 26) can be implemented by looping over all the points in Ω i . A more efficient way is to store the normalised weights in Eq. ( 26) in a matrix and implement the convolution via a matrix-vector product.

Results and discussion

In this section, we present and discuss the results of two numerical studies. The first is on the effect of the length scale contrast method for the gradient damage model. Next, we compare this method with the integral averaging model. The comparison is made based on their performance in providing a delocalised damage field and thus giving a numerically calculated fracture energy G n c in close agreement with two ideal solutions: the rate-independent analytical solution derived above (denoted by 'Rate indep.') and a ratedependent uniform damage strain solution denoted by 'Rate dep.'. The latter is easy to obtain in our 1D numerical setting: it only requires running simulations of the underlying Fig. 3 Nonlocal damage profiles for various length scale (λ in /λ out ) 2 contrast at an instant when global stress σ has dropped to 0.007σ c local damage model (φ nl = φ l ) without imperfection. Given the piecewise uniform nature of this problem, the numerical solution for its calculation does not show localisation and thus does not need any regularisation.

Model parameters used

The numerical studies are performed on the 1D bar of Fig. 2, discretised by 1000 uniformly spaced Fourier grid points. Young's modulus of the material is E = 100 GPa and a critical stress σ c = 1 MPa is used in all the simulations performed. The elastic properties throughout the bar are assumed to be same. The bar is loaded at an average strain rate ε = 10 -5 s -1 for 5 s in time steps of Δt = 0.01 s. Since the considered properties are uniform, we introduce a small imperfection by means of a 1% reduction of the critical stress σ c for the grid points in the interval x ∈ [-l/20, l/20] within the interphase band, in order to trigger localisation.

Effect of length scale contrast in the gradient based model

In this part, we assess the performance of enforcing the interface condition (Eq. ( 17)) by the length scale contrast method. For this study, apart from the parameters described above, an interphase band thickness of l = 0.05L, reference strain rate εd0 = 0.1 s -1 and the final damage strain ε f = 0.001 were used. The nonlocal length scale used inside the band was λ in = l while outside it, the value λ out were varied to investigate the influence of the contrast (λ in /λ out ).

Figure 3 depicts the distribution of nonlocal damage φ nl inside the interphase band and the region outside, to which it has spread significantly. Due to symmetry, only half of the band (and its surroundings) is shown. The profiles are plotted for the cases (λ in /λ out ) 2 = 10 1 , 10 2 , 10 3 , 10 4 . The nonlocal damage variable field extends more than 1.5 times the interphase band thickness into the bulk region on either side of the interphase for (λ in /λ out ) 2 = 10 1 . This is due to the poor approximation of Eq. ( 17) by Eq. ( 18) with this value of length scale contrast. However, this extension of the interphase damage into the Fig. 4 Effect of the interface condition at the edges of the interphase on the global stress-strain response bulk drops systematically for higher contrasts. The maximum contrast that the numerical solver for the damage equation could handle was (λ in /λ out ) 2 = 10 4 . For this contrast, we observe that the zone of spread of interphase damage into the bulk has reduced to less than 0.2 l. Figure 3 also shows that for high contrast ratios, the non-uniformity of φ nl caused by the imperfection decreases and the damage profile becomes increasingly uniform inside the band. Consistent with Eq. ( 18), the slope of φ nl at the two edges of the interphase tends to zero as the length scale contrast is increased. For the lower ratios, the interaction with the undamaged region next to the band results in a slower evolution of φ nl near the edge of the band compared to its centre.

The effect of this slowing down is observable in the averaged stress-strain response, see Fig. 4. As a result of it, the initial softening slope is less steep. This in turn results in an over-prediction of the amount of dissipation during the damage process as compared with that for the cases with a sufficient contrast, i.e. (λ in /λ out ) 2 = 10 3 -10 4 .

Figure 5 (left) shows the calculated fracture energy G n c relative to the one predicted by Eq. ( 14). Remember that G c was derived for the rate-independent limit. Hence, we should expect the rate-dependent, computed G n c to approach G c , but not necessarily to asymptote to it. A much closer approximation of the rate dependent ideal solution (dashed curve) can be expected, as it captures the additional dissipation due to rate effects. Indeed, the ratio G n c /G c drops to less than 1.1 for the highest length scale contrast applied. This higher accuracy however entails significantly more iterations (Fig. 5 (right)) required to solve Eq. ( 23): this number goes up by a factor of 10 as the contrast is increased from (λ in /λ out ) 2 = 10 1 to 10 4 . This is essentially due to the d(λ 2 d dx ) dx structure of Eq. ( 16), which brings first order gradients of φ nl into the damage Eq. ( 16).

( λ in λ out ) 2 G n c G c 10 1 10 2 10 3 10 4 0 1 2 •10 4 ( λ in λ out ) 2

Iteration count

Fig. 5 The reduction of fracture energy error (left) and increase in linear iterations (right) as a function of length scale contrast

Comparison of gradient and integral nonlocal damage models

The analysis in "Analytical solution for the rate-independent case" section showed that for a uniform damage strain distribution in the rate-independent limit, choosing ε f according to Eq. ( 15), a constant fracture energy can be obtained, irrespective of the interphase band thickness l. In this section, we want to test the same hypothesis but in a numerical setting, where localisation is triggered by the presence of an imperfection (as discussed above). This is a test of the effectiveness of the two types of nonlocality in regularizing the damage (and damage strain rates) and thus facilitating the use of simple scaling Eq. ( 15) for the interphase band. Different interphase band thicknesses (l/L = 0.010, 0.015, 0.020, 0.025, 0.030) are studied. εd0 = 1.0 s -1 and ε f = 5 × 10 -5 L/l in combination with the parameters in "Model parameters used" section are used. For both models, three different values of the nonlocal length scale are used: λ/l = 0.5, 1.0, and 1.5, where for the gradient damage based approach λ in = λ and λ out is varied such that a constant length scale contrast (λ in /λ out ) 2 = 10 4 is obtained.

Figures 6 and7 show the damage profiles obtained towards the end of the damage process for the gradient and integral model, respectively. The damage distribution obtained for the gradient damage model, shown in Fig. 6, is increasingly uniform inside the interphase band as the length scale λ is increased, for all the band thicknesses l. However, a slight gradient remains even for λ/l = 1.5. This is due to the imperfect insulation of the interphase damage from the surrounding (undamaged) bulk, as discussed in "Effect of length scale contrast in the gradient based model" section. The integral nonlocal model, in Fig. 7, shows a nearly constant damage inside the band for λ/l ≥ 1.0-again for all l/L. For the smallest length scale considered, λ/l = 0.5, some of the non-uniformity introduced by the imperfection remains.

A more discriminating result is the damage strain rate εd , which is shown in Figs. 8 and9-again towards the end of the damage process and for the gradient and integral approach respectively. Here we do observe a significant amount of non-uniformity in all cases, demonstrating that the damage strain does localise late in the damage evolution process, triggered by the imperfection. The effect is however significantly milder for the integral model (Fig. 9) compared with the gradient model (Fig. 8). Particularly for the integral

ε E ε E ε
Fig. 10 The overall stress-strain response of the gradient damage model for different interphase band thicknesses and nonlocal length scales the applied average strain ε for all cases considered. Also shown (in black) is the underlying rate-independent response as given (analytically) by Eq. [START_REF] Xu | Void nucleation by inclusion debonding in a crystal matrix[END_REF]. For both types of nonlocality, and for all l/L shown, the computed stress-average strain responses for l/L = 1.0 and 1.5 practically coincide. They furthermore deviate only slightly from the analytical solution. This shows that for these values of the nonlocal length scale both regularisation approaches are effective in rendering the global response objective with respect to the interphase band thickness, despite the non-uniformity observed in φ nl (Figs. 6,7) and, particularly, the damage strain rate (Figs. 8,9). The smallest value of λ, λ = 0.5 l, clearly is insufficient, as the softening responses computed for it are systematically steeper than for larger values (and than the analytical solution).

For a more quantitative assessment of the objectivity of the computed response with respect to the interphase band thickness adopted, the variation in the calculated fracture energy G n c , normalised by the ideal input fracture energy G c , is shown for the gradient damage method and the integral averaging method in Fig. 12. In the rate-independent limit, one would like to observe this ratio to be constant at unity. Given the fact that the numerical model is rate-dependent, however, a slight deviation should be anticipatedwhich ideally would be constant. To illustrate the effect of the rate dependence, numerical solutions for the rate dependent case with piecewise constant damage are included in both diagrams as black dashed curves. Note that these are slightly above G n c /G c = 1 and show a downward trend, which is due to the fact that post-peak the strain rate in the interphase band varies with the interphase band thickness. The fracture energy computed by the gradient damage approach show a discrepancy with both reference values (i.e. rate-independent and rate-dependent) on the order of 10%. For λ/l = 0.5, the fracture energy is under-estimated, whereas the larger values of λ result in an over-estimation. In each of these cases, the deviation furthermore varies with l/L. The integral nonlocal model performs much better. For λ/l = 0.5, which we already observed to be too small to obtain a uniform damage in Fig. 7, it under-estimates G c by up to 10% (Fig. 12). But for λ/l ≥ 1.0, it over-estimates G c by ≤ 5%, in a way which is consistent with the rate-dependent reference solution: the difference between the regularised solutions for λ/l = 1.0 and 1.5 and the piecewise constant reference solution is less than 1%. For practical purposes, this is more than satisfactory and the slight trend with l/L, due to the rate-sensitivity, should also pose no problem.

Conclusion and outlook

In this paper, the issue of nonlocality associated with a method to incorporate interface decohesion at polycrystalline interfaces approximated by an FFT based spectral method was discussed. Interfaces were approximated as interphase bands. The softening nature of decohesion required the use of nonlocality within the interphase domain. The applicability and performance of integral and gradient damage based nonlocal averaging methods have been discussed.

For the gradient damage approach, it was found that the mechanical response and dissipation depends on the accuracy with which the flux boundary condition at the edges of the interphase domain can be enforced. In order to restrict the damage to the interphase only, a flux free condition at the interphase boundaries was used. Since FFT solvers do not allow solving an equation on an irregular domain, in the gradient approach, a contrast in the nonlocal length scale was used to approximate the flux free boundary condition at the interphase edges. The error in the fracture energy reduces upon increasing this length scale contrast. However, a high contrast renders the gradient based approach highly heterogeneous and entails numerical problems.

The implementation of the integral approach on the other hand was straightforward and effective. It only requires storage of the nonlocal weights in a matrix that can be optimized using sparse storage. Regularisation length scale values equal to or greater than the interphase band thickness were found to give accurate predictions for the fracture energy, largely independent of the (arbitrary) interphase band thickness l. The slight (≤ 1%) variation which remains is due to the fact that the scaling of ε f according to Eq. ( 15) does not take into account the strain rate sensitivity of the damage model. From the current study it is very clear that integral approach offers more advantage from the computational efficiency point of view. We expect the same advantages to carry over in multidimensional cases on periodic microstructures. Nevertheless, this still remains to be tested.

We wish to emphasize that, although it enables a rigorous and transparent comparison, the 1D problem considered here may not reveal all complexities that could be encountered in two or three dimensions. For instance, care should be taken that the nonlocality introduced to homogenise the damage across the interphase band thickness does not affect the propagation of damage along the band in an unrealistic manner. One possible way to avoid this is by introducing anisotropy in the nonlocal averaging. In the gradient damage approach this can be achieved by having a tensorial form for the nonlocal length scale, while for the integral approach an orientation dependent averaging kernel can be used. Modelling decohesion of polycrystalline interfaces will require a proper treatment of triple junctions to avoid the issue of non-unique interface normals. Furthermore, a method to couple interface damage with the bulk damage in voxel based models [START_REF] Nguyen | A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF] can also be explored to model more complex crack patterns-kinking and branching of cracks into the bulk. These issues are currently under investigation and will comprise future works.

Fig. 1

 1 Fig. 1 Methodology for modelling a polycrystalline interface in finite element method (a 1 -a 3 ) and FFT based solvers (b 1 -b 3 ). a 1 shows a polycrystal with interface elements embedded at sharp boundaries. a 2 depicts an interface conformal finite element mesh with cohesive zones that open during loading according to the traction (σ )-separation (δ) law shown in a 3 with characteristic fracture energy G c . b 1 shows the same polycrystal with interfacial bands representing the sharp interfaces. An example constitutive behaviour for the material points (within the volumetric interphase) assigned dissipation density g c is shown in b 3

Fig. 6 5 Fig. 7 Fig. 8 Fig. 9

 65789 Fig. 6 Nonlocal damage φ nl distribution at ε E σc = 4 of the gradient damage model for different interphase thicknesses and nonlocal length scales

E

Fig. 11 Fig. 12

 1112 Fig.11The overall stress-strain response of the integral averaging model for different interphase band thicknesses and nonlocal length scales

Acknowledgements

The authors would like to acknowledge the fruitful discussions with Prof. Bob Svendsen of RWTH Aachen University and Max Planck Institut für Eisenforschung, Düsseldorf, Germany. The work presented is part of the research done under the project S22.2.1349b in the framework of the research program Materials innovation institute M2i.

Funding

This research is supported by Tata Steel Europe through the Materials innovation institute (M2i) and Netherlands Organisation for Scientific Research (NWO), under the grant number STW 13358.

Availability of data and materials

Not applicable.

Abbreviations

FFT: fast Fourier transform; FE: finite element; FEM: finite element method; TSL: traction separation law.

Authors' contributions

LS carried out most of the study and drafted the manuscript. LS and PS developed the methodology and implemented it. RHJP, FR and MGDG conceived of the study and also participated in its design and coordination in addition to critically reviewing the manuscript. All authors read and approved the final manuscript.

Author details 1 Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands, 2 Microstructure Physics and Alloy Design, Max-Planck Institut für Eisenforschung, Max-Planck-Str. 1, 40237 Düsseldorf, Germany.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Mathilde Perray

Djenebou Traore

Lucas Riegel

Daniela Rojas Castro

Bruno Spire

Marion Mora

Adam Yattassaye

Marie Préau

Pathways of women peer researchers living with HIV in the evaluation of the Gundo-So programme

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

Gundo-So

DISCUSSION

• These results highlighted an activist commitment that was built up through several similar stages among the different facilitators. In particular, through their involvement in the community-based clinics, the increase in skills and knowledge through this commitment and the possibility of professionalization that follows were noted.

RESULTS

METHODS

•

  A community-based program • Supporting women living with HIV (WLHIV) in Mali regarding serostatus disclosure management • ANRS-12373/Gundo-So was a randomized trial implemented in 6 community-based clinics in Mali in 2019/2021 Community programs allow for the development of both psychological empowerment of community workers and organisational empowerment Pathways of women peer researchers living with HIV in the evaluation of the Gundo-So programme Perray, M., Traore, D., Riegel, L., Rojas Castro,D., Spire, B., Mora, M., Yattassaye, A., Préau, M


C Bourges

Stéphane Chevalier
email: stephane.chevalier@u-bordeaux.fr

Jeremie Maire
email: jeremie.maire@u-bordeaux.fr

Alain Sommier

Christophe Pradere

S Dilhaire

Mid-infrared spectroscopic

explanation and prediction of the spectroscopic behavior of the thermotransmittance coefficient in the mid-IR range

Main text

Although microsystems are increasingly subject to thermal challenges, there is no reliable contactless way to measure the temperature within a device for all types of materials.

In most cases, surface measurements performed by detecting the reflection of a light beam [START_REF] Norris | Femtosecond pump-probe nondestructive examination of materials (invited)[END_REF][START_REF] Anufriev | Quasi-Ballistic Heat Conduction due to Lévy Phonon Flights in Silicon Nanowires[END_REF][START_REF] Yazawa | Fast Transient and Steady State Thermal Imaging of CMOS Integrated Circuit Chips Considering Package Thermal Boundaries[END_REF][START_REF] Cahill | Thermometry and Thermal Transport in Micro/Nanoscale Solid-State Devices and Structures[END_REF] or the infrared (IR) [START_REF] Cernuschi | In-plane thermal diffusivity evaluation by infrared thermography[END_REF][START_REF] Groz | Thermal resistance field estimations from IR thermography using multiscale Bayesian inference[END_REF][START_REF] Breitenstein | Microscopic lock-in thermography investigation of leakage sites in integrated circuits[END_REF] emission from the surface are coupled to thermal models and inverse methods to estimate the temperature in buried layers and thermal properties such as thermal diffusivity [START_REF] Bedoya | Measurement of in-plane thermal diffusivity of solids moving at constant velocity using laser spot infrared thermography[END_REF] or thermal contact resistance between layers [START_REF] Chassain | Bayesian Estimation of Thermal Properties Using Periodically Pulsed Photothermal Radiometry: A Focus on Interfacial Thermal Resistances between Layers[END_REF]. However, these approaches are often developed for opaque surfaces and quantitative measurements in semitransparent media proves more challenging. In these materials, the volume of the sample contributes to the proper emission measured in IR thermography and an apparent emissivity [START_REF] Gieseler | Apparent emissivity measurement of semi-transparent materials part 2: Theoretical concept[END_REF] or emittance [START_REF] Jeon | Simultaneous measurement of emittance, transmittance, and reflectance of semitransparent materials at elevated temperature[END_REF] needs to be defined. Furthermore, the camera needs to be calibrated and background noise from the environment need to be filtered out [START_REF] Gieseler | Apparent emissivity measurement of semi-transparent materials part 1: Experimental realization[END_REF]. Although IR thermography is therefore not yet widely used, several works have either focused on measuring the emittance [START_REF] Gieseler | Apparent emissivity measurement of semi-transparent materials part 1: Experimental realization[END_REF][START_REF] Rozenbaum | A spectroscopic method to measure the spectral emissivity of semi-transparent materials up to high temperature[END_REF][START_REF] Adibekyan | Review of PTB Measurements on Emissivity, Reflectivity and Transmissivity of Semitransparent Fiber-Reinforced Plastic Composites[END_REF] or on the thermal diffusivity of the material. Despite these studies, developments of contactless measurement of the 3D temperature field in semitransparent media remain necessary. This is especially true in the IR spectral range, particularly from 2 to 5 µm. With the strong development of microsystems semi-transparent in that spectral range, such as microfluidic chips, electronic devices, energy generation and storage devices like micro-batteries and supercapacitors, there is an increasing need for operando temperature measurement within these systems. In a previous work, we have demonstrated the possibility to overcome the challenges inherent to semi-transparent media without the need for This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689

The origin of the thermotransmittance coefficient, 𝜅(𝜆), lies in the variation of the optical properties with temperature, and in particular the reflections at interfaces and the absorption in the material. The Lorentz oscillators model describes the response of the bound electrons to an electric field and allows the calculation of the relative permittivity and therefore the refractive index, based on electrons collectively acting as damped oscillators. The resulting temperature dependence of the refractive index then translates to a temperature dependence of the reflection and transmission. It can therefore appropriately model dielectric materials'

properties in the mid-IR range to obtain the wavelength dependent thermotransmittance coefficient. In our previous work, the wavelength was selected without consideration for the wavelength dependence of the thermotransmittance coefficient [START_REF] Bourgès | Infrared thermotransmittance-based temperature field measurements in semitransparent media[END_REF], which strongly limits the performance of this technique. However, another study focusing on a microfluidic chip highlighted the importance of the choice of wavelength for the simultaneous imaging of temperature and concentration [START_REF] Kakuta | Simultaneous imaging of temperature and concentration of ethanol-water mixtures in microchannel using near-infrared dual-wavelength absorption technique[END_REF]. Thus, it clearly appears that a phenomenological model to predict the thermotransmittance coefficient, and its extensive measurement in the mid-IR range would be of fundamental importance to strengthen the temperature measurements in semitransparent dielectric materials.

Therefore, this work aims at answering these needs by reporting the measurements of 𝜅 as a function of wavelength in the mid IR range, and by proposing a Lorentz oscillators model with thermodependent parameters to predict the value of 𝜅 as a function of wavelength. The experimental setup relies on an IR camera to detect the variations of transmission from an IR source through a Borofloat glass wafer heated at different temperatures [START_REF] Bourgès | Infrared thermotransmittance-based temperature field measurements in semitransparent media[END_REF]. We then introduce a semi-empirical expression of the transmittance based on Lorentz oscillators and compare the results to the transmission spectrum at ambient temperature. Finally, we introduce the temperature dependence in the model and fit the free parameters to the thermotransmittance coefficients measured at different wavelengths to obtain the thermotransmittance coefficient This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689 spectrum between 3100 nm and 4600 nm and provide measurements to justify the choice of measurement wavelength.

The measurement process consists in using IR camera to detect the intensity of an IR beam that is transmitted through the sample of interest. As the temperature of the sample changes, so do the optical properties of the sample and the corresponding variation in the intensity of the transmitted beam is recorded; it can be expressed as follows:

𝛥𝛤(𝑡, ) 𝛤 0 () = 𝜅()𝛥𝑇(𝑡),
Eq. 1

where 𝛤 0 () is the transmittance at ambient temperature 𝑇 0 , 𝛥𝛤(𝑡, ) = 𝛤(𝑇, 𝑡, 𝜆) -𝛤 0 (𝜆) is the variation in the transmitted signal as a function of temperature at time t, 𝛥𝑇(𝑡) = 𝑇(𝑡) -𝑇 0 is the variation of temperature as a function of the same time 𝑡. Note that since the optical properties depend on the wavelength , so does the thermotransmittance coefficient and therefore the measured relative transmittance variation.

To measure the thermotransmittance signal, we use the setup shown in Figure 1. The illumination wavelength is selected using a diffraction grating monochromator (Bentham Instruments, TMc300) from a stabilized broadband infrared lamp (IR-Si217), before passing through a beam expander to adjust the diameter of the IR beam to match the sample dimensions.

The transmitted IR beam is then detected using an infrared camera (FLIR SC7000), equipped with an indium-antimonide (InSb) detector of 512×640 pixels, a pixel pitch of 15 µm, and a 50 mm focal length lens. The spectral range of the camera is λ ∈ [2.5 -5.5] µm, and the optical resolution is 107 µm/pixel in our experimental configuration.

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689 m 2 /s [START_REF] Williamson | A Compilation of thermal property data for computer heat-conduction calculations[END_REF][START_REF] Churchill | Correlating equations for laminar and turbulent free convection from a vertical plate[END_REF]. The procedure to obtain the thermal properties of the sample is detailed in our previous work [START_REF] Bourgès | Infrared thermotransmittance-based temperature field measurements in semitransparent media[END_REF]. The measured temperature variation field Δ𝑇(𝑟) after demodulating the thermotransmittance signal at frequency 𝑓 𝑇 is shown in Fig. 1c.

The first step in comprehending the temperature dependency of transmission though a semitransparent sample consists in modeling the optical properties and calibrating the setup.

The optical behavior of the materials can then be predicted and the most thermally sensitive wavelengths can be identified to optimize experimental studies.

To determine the total transmittance of a slide of thickness 𝐿, we need to consider (i) the reflections at the sample interfaces, (ii) the absorption in the material and (iii) the contribution This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689

of multiple reflections that occur within the slide [START_REF] Hecht | Optics: International Edition[END_REF]. The total transmission coefficient, Γtot, is obtained by summing the contributions from multiple reflections as given in Equation 2:

Γ 𝑡𝑜𝑡 (𝜆) = (1 -𝑅(𝜆)) 2 × 𝑒 -α(𝜆)𝐿 × ∑ 𝑅(𝜆) 2(𝑞-1) ∞ 𝑞=1 𝑒 -2(𝑞-1)α(𝜆)𝐿 , Eq. 2
where the reflection coefficient R and the absorption  are given by the following equations:

𝑅(𝜔) = (𝑛 ′ (𝜔) -1) 2 + 𝑛 ′′2 (𝜔) (𝑛 ′ (𝜔) + 1) 2 + 𝑛 ′′2 (𝜔) , Eq. 3 α(𝜔) = 2ω𝑛 ′′ (𝜔) 𝑐 = 4π𝑛 ′′ (𝜔) λ ,
Eq. 4

where n' and n'' stand for the real and imaginary part of the refractive index 𝑛 ̃(ω) = 𝑛 ′ (𝜔) + 𝑖𝑛 ′′ (𝜔) and where the angular frequency is given by 𝜔 = 2𝜋𝑐/𝜆 , and c is the speed of light in vacuum (m/s).

Both parts can be calculated from the relative permittivity ε r ̃(ω) = ε ′ (𝜔) + iε ′′ (𝜔) with the following set of equations:

𝑛 ′ (𝜔) = √ √ε ′2 (𝜔) + ε ′′2 (𝜔) + ε ′ (𝜔) 2 ,
Eq. 5

𝑛 ′′ (𝜔) = √ √ε ′2 (𝜔) + ε ′′2 (𝜔) -ε ′ (𝜔) 2 .
Eq. 6

In our modeling approach, we assume that the dielectric properties of the transparent medium are well-described by a classical Lorentz oscillators' behavior that is characterized by the total permittivity. This complex value is the sum the dielectric constant 𝜀 ∞ and the contributions of This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689

several resonances processes, which can be modeled as a sum of oscillators. Each oscillator is then characterized by its resonance and plasma frequencies 𝜔 𝑂 and 𝜔 𝑝 , as well as a damping factor 𝛾. The resulting expression for the permittivity is given in Equation 7:

ε 𝑟 ̃(ω) = ε ∞ + ∑ ω 𝑝,𝑗 2 ω 0,𝑗 2 -ω 2 + 𝑖ωγ 𝑗 𝑗 Eq. 7
These parameters are not known a priori, and transmission measurements are required to set them. We use the setup shown in Fig. 1 at ambient temperature and sweep the wavelength of the monochromator over the detection range of our camera. The resulting transmission and reflection spectra are shown in Fig. 2. Several features are visible, such as three dips of varying depth in the transmission and a progressive decrease of the transmission at large wavelengths.

Therefore, we choose the resonant frequencies of the oscillators so that they correspond to the three dips as well as one chosen arbitrarily at a wavelength beyond 5000 nm to model the progressive decrease in transmission. The chosen values are reported in Table 1 and are in good agreement with the literature [START_REF] Davis | Quantitative infrared spectroscopic measurement of hydroxyl concentrations in silica glass[END_REF][START_REF] Tan | Determination of refractive index of silica glass for infrared wavelengths by IR spectroscopy[END_REF]. Although our values agree, note that the bound electrons resonances do not necessarily explain the physical or chemical origin of these bands, which have been explained at length in previous works, such as OH groups in oxide glasses [START_REF] Navarra | OH-related infrared absorption bands in oxide glasses[END_REF] or B2O3 and Na2O/K2O dopants [START_REF] Izawa | Optical attenuation in pure and doped fused silica in the ir wavelength region[END_REF].

Wavelength (nm) 2750 3620 3960 5500

Resonant frequency 𝜔 0 (rad/s) 6.85 × 10 14 5.20 × 10 14 4.76 × 10 14 3.43 × 10 14

Table 1. Identified resonant frequencies

Next, we perform an optimization to estimate all the remaining unknowns, namely the plasma frequencies, the damping rates and the dielectric constant. The latter was found to be 𝜀 ∞ = This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The result of the fit is shown in Fig. 2 alongside the measured spectrum. The approach of assimilating spring-bound electrons in the Lorentz oscillators model with several resonance frequencies gives results in good agreement with the experimental values, as the position, width and depth of the dips in transmission are well reproduced. However, the discrepancy observed between the measurement and the model for the transmission peak around 3300 nm suggests that other types of contributions are missing from the model such as that of free electrons [START_REF] Tompkins | Handbook of Ellipsometry[END_REF],

PLEASE CITE THIS ARTICLE AS

or contributions that could be accounted for by a semi-quantic dielectric function model [START_REF] Sousa Meneses | Contribution of Semi-Quantum Dielectric Function Models to the Analysis of Infrared Spectra[END_REF].

Another possibility is that the visible and near infrared absorption peaks are not considered in the model, which could have a slight impact on the peak at 3300 nm. Several more complex models have been proposed [START_REF] Kitamura | Optical constants of silica glass from extreme This is the author's peer reviewed[END_REF][START_REF] Efimov | IR absorption spectra of vitreous silica and silicate glasses: The nature of bands in the 1300 to 5000 cm-1 region[END_REF][START_REF] Sousa Meneses | Causal Voigt profile for modeling reflectivity spectra of glasses[END_REF][START_REF] Sousa Meneses | Structure and lattice dynamics of binary lead silicate glasses investigated by infrared spectroscopy[END_REF] that account for various phenomena, and some of the spectral features were attributed to the fabrication process. More complex models have been developed in the literature to account for all the phenomena Nonetheless, despite the highly non-monotonous transmission spectrum of our sample, the absolute error, calculated as the difference between the measured and modeled value, remains below 10%. The best fit therefore This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Once we have identified all the parameters at ambient temperature, we investigate the thermal dependency of the refractive index and its impact on the transmittance of the material. For simplicity, we assume that the resonant frequencies, plasma frequencies and damping rates are first-order functions of temperature given by Eq. 8-10.

PLEASE CITE THIS ARTICLE AS

ω 0 (𝑇) = ω 0 (𝑇 0 ) × (1 + α ω 0 Δ𝑇)

Eq. 8

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689

ω 𝑝 (𝑇) = ω 𝑝 (𝑇 0 ) × (1 + α ω 𝑝 Δ𝑇)
Eq. 9

γ(𝑇) = γ(𝑇 0 ) × (1 + α γ Δ𝑇)

Eq. 10

where ∆𝑇 = 𝑇 -𝑇 0 is the temperature variation (K) and 𝛼 𝜔 0 , 𝛼 𝜔 𝑝 and 𝛼 𝛾 the coefficients of thermal variation of the resonant frequencies, plasma frequencies and damping rates, respectively (K -1 ). The values at 𝑇 = 𝑇 0 correspond to the room-temperature values obtained in the previous section. For the plasma frequencies, the temperature dependence is primarily caused by the thermal expansion of the material that affects the electron density ne(T) [START_REF] Block | Tracking ultrafast hot-electron diffusion in space and time by ultrafast thermomodulation microscopy[END_REF].

To estimate the temperature coefficients for each resonance identified in the previous section, we need the experimental values of the thermotransmittance coefficient. We extract the thermotransmittance coefficient at each wavelength of interest using the method described in our previous work from the thermal fields shown in Fig. 1b,c. The variation of transmittance as a function of temperature is presented in Fig. 3a for the wavelength 𝜆 = 3300 nm and clearly shows a linear dependence. A linear regression then provides the thermotransmittance coefficient. Note that the higher dispersion in the thermotransmittance signal for large temperature variations correspond to pixels close to the contact area with the Peltier module used to modulate the temperature, where there is uncertainty about the exact position of the edge. We repeat the same procedure at selected wavelengths from 3100 nm to 4600 nm and plot the values of the thermotransmittance coefficient as a function of the wavelength.

The temperature-dependent model is then fitted to the thermotransmittance signal as a function of wavelength with a simplex minimization algorithm [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF]. The resulting temperature coefficients are shown in Table 3. The modeled thermotransmittance as a function of wavelength is plotted in the range from 3000 nm to 5000 nm in Fig. 3b alongside the experimental values. The good agreement between the experimental values and the best fit from the model attests that our model is a good approximation to calculate the thermo-optical This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689

PLEASE CITE THIS ARTICLE AS

Since the intensity of the transmittance variation with temperature is represented by the absolute value of the thermotransmittance coefficient, we see from Fig. 3 that the relative variation of transmission as a function of temperature will increase for longer wavelengths in the range investigated in this study. More specifically, the thermotransmittance scales with the absorbance and an illumination wavelength matching an absorbance peakequivalent to a transmittance dipgives the highest value of the thermotransmittance coefficient. To elucidate the reason behind this, we first introduce a first order temperature dependence in the reflectance and in the absorption coefficient, similarly to what has been done for the optical parameters.

After injecting these new expressions in Eq. 2, we calculate the temperature dependent relative variation of transmittance, that is our thermotransmittance signal. The resulting expression can then be linearized to obtain

𝛥𝛤 𝛤 0 ≈ ( 2𝑅 0 𝜅 𝑟 1 -𝑅 0 -𝛼𝜅 𝛼 𝐿) 𝛥𝑇, Eq. 11
where 𝛤 0 is the initial transmittance at ambient temperature, 𝛥𝛤 is the variation of transmittance for a variation of temperature 𝛥𝑇, 𝜅 𝑟 and 𝜅 𝛼 are phenomenon-specific thermotransmittance coefficients from the reflectance and absorbance, with 𝜅 𝑟 expressed in K -1 and 𝜅 𝛼 expressed in K -1 . We can see from this expression that by measuring samples with identical properties but different thicknesses, it will be possible to accurately measure the relative contributions of reflectance and absorbance to the measured signal, which will be investigated in a future work.

Therefore, measurements will strongly vary based on the thickness of the sample. All the values of the thermotransmittance coefficient are given here in K -1 , and are thus valid for a given thickness of the materialthe thickness dependence of this coefficient will be the subject of a further study.

These observations do not necessarily mean that longer wavelengths are more appropriate to perform thermotransmittance measurements, as the response of the whole acquisition chain has to be considered (see its spectral transfer function in Fig. S5 of the This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689 supplementary material). The measured signal variation indeed combines the relative variation due to temperature changes, but also the absolute value of the transmitted signal. The product of these two values, shown in Fig. 3c, shows a large dispersion in our particular sample but, although several wavelengths, both short and long, display high values, no single one stands out as the best choice for measurements. The measured intensity variation also depends on the transfer function of the whole acquisition chain, among which the sensitivity of the detector.

The total signal output by the camera as a function of wavelength, given in Fig. 3c, clearly highlights the best wavelengths to perform thermotransmittance measurements with our experimental setup, here around 4000 nm and 4400 nm for this particular sample. It is interesting to note the decrease in signal between the two curves of Fig. 3c around 4200 nm, which corresponds to the absorption peak of the CO2 present in the air. One of the interests of these spectroscopic measurements therefore lies in finding the wavelength at which the signal to noise ratio will be optimal.

In conclusion, we developed a semi-empirical model to estimate the temperature-dependent optical properties of dielectric materials, such as glass, in the mid-IR spectral range. The model is combined with thermotransmittance measurements that consist in measuring the temperaturedependent transmission through the material, to obtain the transmission and reflection spectra, as well as the thermotransmittance coefficient. This approach completes our previous work in which we have demonstrated the ability of thermotransmittance measurements to extract the temperature fields and thermal diffusivity of semi-transparent materials in the mid-IR spectral range. This work sheds light on how the thermotransmittance signal varies with the source wavelength and acquisition chain sensitivity.

More specifically, we have introduced here a formalism based on Lorentz oscillators and shown that a fit to the experimental data could help us choose an ideal wavelength for This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689 thermotransmittance measurements. Even in the restricted spectral range of study in this work, the thermotransmittance vary by as much as a factor three, showing that there are potentially wavelengths at which some materials are sensitive and others not. This discrepancy can allow to distinguish the origin of the thermotransmittance signal at specific wavelengths in heterogeneous structures. Derivating our model further can also help predict the dependence over other parameters, such as the sample thickness. With this thickness dependence, it becomes possible to quantitatively assess the proportion of the thermotransmittance signal stemming from reflectance and that coming from the absorbance of the material, which will be detailed in a future work. Overall, this work expands our previous demonstration of thermotransmittance to spectroscopic measurements that combine with a simple semi-empirical model to optimize temperature imaging in semi-transparent media.

Figure 1 .

 1 Figure 1. a) Thermotransmittance measurement system. The light source, comprising an IR lamp and a monochromator, shines through a beam expander and through the sample, and is detected by an IR camera. The camera acquisition is synchronized to an optical chopper to alternatively acquire images with and without the IR beam. insert: The glass sample mounted on the annular Peltier module. b, c) Amplitude field of the modeled and measured temperatures.

  DOI: 10.1063/5.0176689 provides a basis for introducing the temperature dependence in the coefficients of the semiempirical model.

Figure 2 .

 2 Figure 2. Transmittance and reflectance spectra. a) Spectral transmittance and reflectance of the Borofloat wafer measured experimentally (solid lines) and modeled (dashed lines). The three resonances corresponding to the dips in transmittance are indicated with the dashed lines. b) Absolute error for both the transmittance and reflectance spectra.

4 Table 3 .Figure 3 .

 433 Figure 3. a) Relative variation of the transmittance as a function of the calculated temperature. b) Extracted thermotransmittance coefficient for different wavelengths and comparison to the temperature-dependent optical model output. c) Total transmitted signal arriving on the detector for a unitary temperature increase (circles, left) and corresponding signal weighed by the transfer function of the whole acquisition chain (triangles, right).

DOI: 10.1063/5.0176689

	2.6 ± 0.3, again in agreement with the literature [31], and the remaining parameter values are
	given in Table 2. The details of the optimization and sensitivities to the different parameters are
	given in Supplementary Materials.
	Wavelength 𝜆 (nm)	Plasma pulsation 𝜔 𝑝 (rad/s)	Damping rate 𝛾 (s -1 )
	2750	3.1 ± 0.3 × 10 12	3.3 ± 0.8 × 10 13
	3620	4.8 ± 0.2 × 10 12	2.5 ± 0.3 × 10 13
	3960	4.2 ± 0.5 × 10 12	3.8 ± 0.8 × 10 13
	5500	1.4 ± 0.4 × 10 13	1.2 ± 0.6 × 10 12

Table 2 .

 2 Best fit values of the plasma frequencies and the damping rates.

Acknowledgments and statements

This work was supported by the French National Research Agency (ANR) (Grant ANR-22-CE50-0015). This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0176689

The authors declare no conflict of interest.

Supplementary Material

The Supplementary Material presents a sensitivity analysis of the model and the fit parameters.