Tridimensional multiphysics model for the study of photo-induced thermal effects in arbitrary nanostructures

In the present paper, we detail the implementation of a numerical scheme based on the Finite Element Method (FEM) dedicated to a tridimensional investigation of photo-induced thermal effects in arbitrary nano-structures. The distribution of Joule losses resulting from the scattering of an incident wave by an arbitrary object embedded in a multilayered media is used as source of a conductive thermal transient problem. It is shown that an appropriate and rigorous formulation of the FEM consists in reducing the electromagnetic scattering problem to a radiative one whose sources are localized inside the scatterer. This approach makes the calculation very tractable. Its advantage compared to other existing methods lies in its complete independence towards the geometric, optical and thermal properties of both the scatterer and the medium in which it lies. Among the wide range of domain of application of this numerical scheme, we illustrate its relevance when applied to two typical cases of laser damage of optical components in high power applications. [DOI: 10.2971/jeos.2011.11037]


INTRODUCTION
On the one hand, a precise control and enhancement of the absorption of metallic nano-structures can lead to useful applications like hyperthermia of tumors [1,2], targeted destruction of microorganisms [3], triggered drugs or genes release [4], chemical vapor deposition [5], guiding [6] and mixing [7] of fluidic flow, bio-imaging [8]. . .On the other hand, this very same effect has proven to be harmful and costly in high-power applications where local scatterers and absorbers can enhance the local electric field and induce strong local temperature rise leading [9,10] to the damage of expensive large optical components.But in any case, there has recently been a growing interest in photo-induced thermal effects in various domains where metallic nano-structures are at stake.
The rigorous study of this multi-physics effect can be decomposed in two phases: (i) solving of the vector Maxwell equations in order to determine the distribution of the electromagnetic field scattered by the lossy defect(s), then (ii) solving the heat equation whose sources depend on the former Joule loss distribution.The electromagnetic problem is very well known and has been extensively treated for simple geometries.One can refer to [11] for a detailed review of the most widely used methods.However, in order to take into account realistically and rigorously the whole physical complexity of the problem (i.e. an arbitrarily shaped, possibly graded index and/or anisotropic scatterer as well as the structure in which it is embedded), a very general method is required.Both the Finite-Difference Time-Domain (FDTD) [12,13] and the Finite Element Method (FEM) [14,15,16] allow the computation of electromagnetic vector fields and are adapted to such arbitrary geometries.The FDTD is a powerful method which operates in the time domain and thus appears to be less adapted to the monochromatic effects investigated in this paper.The FEM is a widely spread method adapted to solve partial differential equations.Its use is particularly adapted to multiphysics models, since different physical phenomena represented by distinct partial differential equations can be solved on the same tetrahedral meshing of the tri-dimensional (3D) computational domain.In this paper, we address the following multi-physics problem: (i) solving the vector Helmholtz propagation equation in frequency domain and (ii) solving a transient heat equation whose 3D distribution of sources is deduced from the former electromagnetic problem.
In a first part, we show some theoretical and practical details about the implementation of the FEM.First of all, an original formulation already introduced in the case of periodic [17,18] problems is adapted to the case of scattering cases.This problem amounts to looking for the vector field scattered by one or several defects of arbitrary opto-geometric properties embedded in an arbitrary multilayered stack enlighten by a plane wave of arbitrary incident angle and polarization in the harmonic domain.One of the main difficulties in this physically unbounded problem is the rigorous truncation of the computation domain whereas the sources of the incoming wave are at infinity.This is overcome through a rigorous treatment of the plane wave sources thanks to an equivalence of the scatter-ing problem with a radiation one whose sources are localized inside the scatterer itself.
We deduce from the complex electric field map the distribution of Joule losses in the lossy regions of the problem, which are used as a source for a problem of heat transfer by conduction assuming thermal temperature-dependant parameters.This later problem is projected on the same mesh than the former problem but one relies on the use of nodal elements.The validity of the electromagnetic treatment is limited by the length of the temporal laser pulse, i.e. as long as the spectral dispersion of the pulse can be neglected.The thermal part is adapted to pulses in the nanosecond range and longer where thermal effects overcome avalanche and photoionization [19].
In a second part, a convergence study of the electromagnetic part is presented by comparing the FEM solution to the closed form solution obtained from a Mie code [20].Eventually, two examples of application of the method are provided, in the field of high power laser interaction with optical materials.Though the issue of laser damage has been known for a long time [21], its intrinsic multifactorial nature makes it extremely difficult to understand.The presented examples illustrate the crucial relevance of taking into account as precisely as possible the full optical, thermal and geometric characteristics of the possible scattering defect(s) in an optical component.In the first example, two absorbing defects embedded in a homogeneous medium are considered and the influence of their closeness on the laser induced damage threshold is investigated.In the second example, an ellipsoidal defect in a KDP lattice is considered.This two cases exhibit drastically different temperature responses depending on the geometric configuration.

Set up of the problem and notations
We denote by x, y and z the unit vectors of the axes of an orthogonal coordinate system Oxyz.In this paper, for the sake of simplicity, the materials are assumed to be isotropic and therefore are optically characterized by their relative permittivity ε and relative permeability µ.Note that the relative reluctivity, i.e the inverse of permeability, is denoted ν.The relative permittivity and relative permeability are represented by complex valued functions which allows the study of lossy materials.The optical scattering system that we are addressing in this paper can be split into the following regions as suggested in Fig. 1: • The superstrate (z > z 0 ) is supposed to be homogeneous, isotropic and lossless, and therefore characterized by its relative permittivity ε + and its relative permeability µ + (= 1/ν + ) and we denote k + := k 0 ε + µ + , where k 0 := ω/c = 2 π/λ.Its thermal properties are denoted ρ + , C + p and κ + (respectively density, specific heat capacity at constant pressure and thermal conductivity) and are functions of the temperature T.
• The multilayered stack (z N < z < z 0 ) is made of N layers which are supposed to be homogeneous and isotropic, and therefore characterized by their relative permittivity ε n , their relative permeability µ n (= 1/ν n ) and their thickness e n .We denote k n := k 0 √ ε n µ n for n integer between 1 and N. Their thermal properties are denoted ρ n , C n p and κ n and are functions of the absolute temperature T.
• The scattering region (z s < z < z s−1 ) is an heterogeneous region constituted of the layer indexed s (ε s , µ s ) and the scatterer itself.The permittivity, permeability and thermal properties of the scatterer are denoted ε d (x, y, z), µ d (x, y, z), ρ d , C d p and κ d , where the superscript d stands for defect.The method described in this paper does work irrespective of whether the scatterers are homogeneous: The permittivity ε d and permeability µ d can vary continuously (gradient index scatterer) or discontinuously (step index scatterer).This scattering region is thus characterized by the scalar fields defined by part ε s (x, y, z) and µ s (x, y, z) (= 1/ν s (x, y, z)), where ε s = ε d and µ s = µ d (resp.ε s = ε s and µ s = µ s ) inside (resp.outside) the defect.The thermal properties are denoted accordingly ρ s , C s p and κ s and are functions of the temperature T.
• The substrate (z < z N ) is supposed to be homogeneous and isotropic and therefore characterized by its relative permittivity ε − and its relative permeability µ − (= 1/ν − ) and we denote k − := k 0 ε − µ − .Its thermal properties are denoted ρ − , C − p and κ − and are functions of the temperature T.
Let us emphasize the fact that the method principles remain unchanged in the case of several scatterers made of distinct geometry and/or material.It is of importance to note that this approach (see Ref. [18]) allows to study diffractive effects as well and what follows can be applied to the study of thermal effect in diffractive optics [22,23].

Electromagnetic problem
We only deal with time-harmonic fields, consequently electric and magnetic fields are represented by the complex vector fields E and H, with a time dependance in exp(−i ω t).The incident field on this structure is denoted: with and where The scattering problem that we address in this paper is therefore to find the solution of Maxwell equations in harmonic regime i.e. the unique solution (E, H) of: such that the scattered field satisfies a so-called Outgoing Waves Condition (OWC [24]).
One can choose to calculate arbitrarily E, since H can be deduced from Eq. (4a).The scattering problem amounts to looking for the unique solution E of the so-called vectorial Helmholtz propagation equation, deduced from Eqs. (4a,4b): such that the scattered field satisfies an OWC.

From a scattering problem to a radiative one with localized sources
According to Fig. 1, the scalar relative permittivity ε and reluctivity ν fields associated to the studied scattering structure can be written using complex-valued functions defined by part and taking into account the notations adopted in Sec.2.1: with ζ = {ε, ν}, z 0 = 0 and z n = − ∑ n l=1 e l for 1 ≤ n ≤ N.
It is now convenient to introduce two functions defined by part ε 1 and ν 1 corresponding to the associated multilayered case (i.e. the same stack without any scatterer) constant over Ox and Oy : with ζ = {ε, ν}.
We denote by E 0 the restriction of E inc to the superstrate region: We are now in a position to define more explicitly the vectorial scattering problem that we are dealing with in this article.It amounts to looking for the unique vector field E solution of: In order to reduce this scattering problem to a radiative one, an intermediary vector field denoted E 1 is necessary and is defined as the unique solution of: The vector field E 1 corresponds to an ancillary problem associated to the general vectorial case of a multilayered stack which can be calculated independently.A general analytical derivation of such a problem can be found for instance in [16].Thus E 1 is from now on considered as a known vector field.It is now apropos to introduce the unknown vector field E s 2 , simply defined as the difference between E and E 1 , which can finally be calculated thanks to the FEM: It is of importance to note that the presence of the superscript s is not fortuitous: As a difference between two scattered fields (Eq.( 11)), E s 2 also satisfies an OWC which is of prime importance in our formulation.By taking into account these new definitions, Eq. ( 9) can be written: where the right-hand member is a vector field which can be interpreted as a known vectorial source term −S 1 (x, y, z) whose support is localized inside the scatterer itself.To prove it, let us introduce the null term defined in Eq. ( 10) and make the use of the linearity of M , which leads to:

Weak form and truncation of physically unbounded regions
The weak form is obtained by multiplying scalarly Eq. ( 9) by weighted vectors E chosen among the ensemble of fields of Integrating by part Eq. ( 14) and making the use of the Green-Ostrogradsky theorem lead to: where n refers to the exterior unit vector normal to the surface The first term of this sum concerns the volume behavior of the unknown vector field whereas the right-hand term can be classically used to set boundary conditions (Dirichlet, Neumann or so-called quasi-periodic Bloch-Floquet conditions).
The solution E d 2 of the weak form associated to the diffraction problem, expressed in its previously defined equivalent radiative form at Eq. ( 12), is the element of L 2 (curl) such that: A set of Perfectly Matched Layers are used in order to truncate the substrate, the superstrate and each layer of the multilayered medium (see [25] for practical implementation of PML adapted to the FEM).In the examples of Sec.3, rectangular PML have been chosen since a cartesian orthogonal coordinate system was chosen at the first place, but spherical PML could be implemented regardless of the principles of the method.For the sake of clarity, the rectangular PMLs set adapted to a bulk are represented in Fig. 2. Since the proposed unknown E s 2 satisfies an OWC, this set of boundary conditions is perfectly sane: E s 2 is radiated from the diffractive element towards the infinite regions of the problem and decays exponentially inside the PMLs.The total field associated to the scattering problem E is deduced at once from Eq. (11).To solve this part of the problem, tetrahedral quadratic edge elements [26,27] (or Whitney 1-form) were used together with the direct solver PARDISO.The size of each element size in each domain is set to λ/(N m e{ε}).The mesh parameter N m represents the approximate number of elements by which a period of the electric field propagating in a media of permittivity ε is sampled.Compared to other FEM formulations for electromagnetic scattering or diffraction problems, the present rigorous approach exhibits several advantages.For instance, locating the sources of the electromagnetic problem inside the PML [28] or on the exterior bounds of the PML leads to the necessity to mesh the PML finely for accuracy.It should also be pointed out that the PMLs can be set as close to the scatterer as desired.Having the sources inside the domain of interest allows to gradually loosen the mesh in the PML without significant loss of accuracy.These last two practical aspects concerning the computation cell make the problem very tractable and represent two considerable assets when tackling large 3D problems.

Computation of losses and absorption cross-section
The classical dimensionless absorptivity σ can be obtained through the computation of the following ratio: where The numerator in Eq. ( 17) clarifies losses in watts inside the considered lossy scatterer and are computed by integrating the Joule effect losses density Q h (x, y, z) (the superscript h stands for harmonic) over the volume V of the lossy scatterer.The denominator normalizes these losses to the incident power, i.e. the time-averaged incident Poynting vector flux across the surface S resulting from the projection of the scatterer edges on a plane parallel to Oxy in the superstrate, whose normal oriented along decreasing values of z is denoted n.Since E 0 is nothing but the plane wave defined at Eqs. (8,2,3), this last term is equal to (A 2 e ε 0 /µ 0 S)/(2 cos(θ 0 )), where s denotes the surface of S in m 2 .Volumes and normal to surfaces being explicitly defined, the absorptivity σ is instantly available once E computed and interpolated between mesh nodes.Numerical examples and comparisons with an independent code are given in Sec. 3.

Temporal domain and thermal coupling
On the one hand, the thermal coupling is possible provided that the incident pulse has a negligible spectral dispersion, which is generally the case for nanosecond pulse lasers addressed in this paper.On the other hand, the plane wave approximation is fully justified since the size of the scatterers is insignificant compared to the waist the incident laser beam.Indeed the temporal pulse shape must be considered when tackling laser-damage effects [29].Let us consider for instance a laser pulse with a gaussian temporal irradiance profile I(t) given in W.m −2 .The pulse is classically defined by its fluence F in J.m −2 , its temporal width at half-height t 1 − t 0 in seconds and we are numerically spanning the temporal domain [t min , t max ].Note that the presented method remains identical for an arbitrary temporal pulse shape.The fluence F of the pulse is defined as the temporal integral of the gaussian irradiance I(t) profile as follows: where τ = (t 1 + t 0 )/2 and δ = (t 1 − t 0 )/(2 2 ln(2)).Let us assume that the electromagnetic problem was solved with an amplitude of electric field corresponding to the maximum irradiance incident on the structure I max , which is worth according to Eq. ( 19): Finally, since losses Q h in harmonic domain are directly proportional to the incident irradiance I, the spatio-temporal distribution of Joule losses are expressed as follows:

Weak form of the heat equation
The complete derivation of the thermal transient problem using finite elements is outside of the scope of this paper and one can find useful details and references about its practical implementation in [30], as well as comparisons between different iterative solvers.The main goal of this part is to strengthens the rigorousness of the coupling between electromagnetic and thermal calculations by taking into account the distribution of the electromagnetic losses Q t , as well as the temperature dependance of the constitutive materials of the system.The classical conductive heat transfer equation in solids that governs the spatio-temporal evolution of the absolute temperature T(x, y, z, t) in K, or Fourier's law: where the density ρ in kg.m −3 , the specific heat capacity at constant pressure C p in J • (kg • K) −1 and the thermal conductivity κ in W • (m • K) −1 are piecewise temperature dependent functions.Thus, these functions can be defined by part using the notations of Sec.2.1 as follows: with ζ = {ρ, C p , κ}.Following the Galerkin method of weighted residual, i.e. the weighted functions are chosen among the shape functions, and integrating by part leads to the following weak form of Eq. ( 22): Note that convective and radiative thermal contributions are note taken into account in this formulation.Considering the scalar nature of T, the temperature T(x, y, z, t) is approximated using nodal elements of the form: Finally, inserting the decomposition of Eq. ( 25) in Eq. ( 24) leads to final matrix system.The stat-of-the-art methods used to solve the resulting differential algebraic equations are fully described in [31,32].

Numerical illustrations
Among the various fields addressed in the introduction by the described method, we chose to illustrate the method with Laser-induced damage typical situations.Laser-induced damage in optical materials is a main issue for high-power density laser systems.The irreversible degradation and growth of damage of optical components under laser irradiation seriously limit the operational conditions of MegaJoule class laser for fusion ignition [33] as well as small and compact laser for airborne or space applications [34].In dielectric materials irradiated by nanosecond pulses, if the wavelength is far from the absorption bands, laser damage in optical materials is mainly initiated by nanometric defects that absorb the laser energy, inducing a fast temperature rise and a subsequent plasma that will lead to mechanical failure due to extreme temperatures and pressures [35].For surfaces, the initiating defects are often linked to polishing residues [36], fractures [37,38], or contamination [39].If the surface is coated, different additional limiting defects can be found in the multilayer system: Nodules [40], substrate contaminants [41], impurities coming from the coating chamber or from the source material [42,43], voids, grain boundaries. . .In the case of crystals, the origin of defect may be linked to the growth condition.Indeed since the growth requires certain additives, the precursor defects may be constituted of atomic impurities (such as Al, Cr, Fe, Si) [10] or dislocations [44].Structural defects that can form absorbing clusters can also be created during the crystal growth [45].The creation or photo-induced defects during laser operation has also been reported in different materials [46].
In order to refine the knowledge on the precursor defects and improve the understanding of the damage initiation process, models of laser-induced damage are needed [47].The development of such models has begun with the first studies on laser damage, for instance with the work of Hopper and Uhlmann [21] and has lead to more refined models in the recent years [48]- [51].Modeling laser damage in optical components requires the ability to deal with complex structures such as optical multilayer systems, which are often the limitation of an optical system, or diffractive optics which are more and more used in high powerful laser systems [52].Laser damage initiation by complex defects such as fractures [38] or nodules has been theoretically addressed with 2D or 3D FDTD electromagnetic calculations developed for a particular case [51].
The present numerical scheme relying on the use of unstructured mesh presents the advantage of being suitable to single or multiple arbitrary defects of any shape, in complex optical structures, by coupling electromagnetic to thermal calculations.

Convergence of the electromagnetic field as a function of the mesh refinement
In this section, values of absorptivity obtained with the FEM (Eq.( 17)) are compared to those obtained with a classical Mie scattering code.Let us consider a spherical defect made of hafnium (Hf of permittivity ε d = 7.04 + 25.30 i) merged in a hafnium dioxide (HfO 2 of permittivity ε s = ε − = ε + = 4) bulk and enlighten by a plane wave at 1064 nm.Fig. 3b represents the relative error between the absorption cross section obtained by FEM and using a Mie code [20].This case can be considered as a typical defect in a coating [53].
As depicted in Fig. 3a, the FEM calculation shows excellent numerical agreement with the Mie scattering code for various values of the radius r de f ect of the spherical defect.Furthermore, as shown on Fig. 3b, the convergence of the FEM scattering results is fast and very accurate values can be obtained as soon as N m ≥ 6.For instance, the relative error on the absorptivity σ with r de f ect = 150 nm is lower than 0.1% for N m = 8.Interestingly enough in terms of computational cost, note that this relative error is only 2.4% when N m is set to 1.The electromagnetic part of the presented method can produce very ac- curate results at low computational cost.Its interest towards a Mie code relies of course in its independence towards the geometry of the diffractive object and the structure in which it is embedded

Two Hf spheres in a HfO 2 bulk
In this section, we are dealing with multiple scatterers.With regard to the single scatterer case, the proposed method does not require any further adjustment: both the piecewise electromagnetic and thermal sources are now non null in the several scatterers.As in Sec.3.1, the scatterers are made of hafnium (Hf of permittivity ε d = 7.0364 + 25.2960 i), merged in a hafnium dioxide (HfO 2 of permittivity ε s = ε − = ε + = 4) bulk and enlighten by a plane wave of wavelength 1064 nm and whose electric field is polarized along the direction of alignment of the spheres.Fig. 4a represents the absorptivity σ 2 of a defect of radius r de f ect calculated in the two scatterers case (normalized by its single scatterer counterparts σ 1 calculated in Sec.3.1) as a function of the radius r de f ect and distance d between the two closest points the spheres.The ratio σ 2 /σ 1 exhibits a flat portion corresponding to the distances d large enough where the defects do not interact between each other.
A closer look at the point (r de f ect , d) = (120, 300) nm (right hand side dashed circle on Fig. 4a) confirms the very weak perturbation of the symmetry of the field interference pattern, as shown on the field maps |E y | in Fig. 4c and |H x | in Fig. 4e.As a consequence, the distribution of losses shown at Fig. 4g is very similar to the one scatterer case, which results in a ratio σ 2 /σ 1 close to 1.However, the mesh plot shown in Fig. 4a exhibits a sharp peak centered at (r de f ect , d) = (40, 2) nm.In this case, because of the transverse polarization of the electric field, the two close defects strongly interact, which results in a strong confining of the electric field shown in Fig. 4b: The maximum value of |E y | is 8 times stronger than the one presented on Fig. 4c.A similar asymmetry can be observed on the magnetic field distribution in Fig. 4d, which results in large off-centered losses plotted on Fig. 4f, up to 9 times the maximum losses reached by the case depicted in Fig. 4g.
We are now comparing the temperature distribution result- ing from the two-scatterers case to its single scatterer counterpart by making the use of the thermal coupling introduced in Sec.2.3.The thermal parameters are the following [54]: and the thermal conductivity, considered linearly dependant on the temperature [54] κ(T) = −0.016* T + 26.596 W • (m • K) −1 .Let us consider an incident gaussian temporal pulse of width at half-height t 1 − t 0 = 10 ns at λ = 1064 nm.The considered temporal range is [t min , t max ] = [0, 30] ns sampled with 1 ns steps.As shown in Fig. 5 the maximum of temperature reached in both 1-and 2-scatterers cases occurs 2 ns after the maximum pulse irradiance I max ≈ 1.16 × 10 11 W.m −2 .However, the temperature maximum reached in the 2-scatterers case is more than twice as high as in the 1-scatterer one (2240 K instead of 1150 K).
The hafnium fusion point is located at 2233 K [55].The model can then be used to study cooperation between defects in laser damage, as they may electromagnetically and thermally cooperate to induce damage [56].

Ellipsoidal defect in a KDP lattice
The second example we have chosen concerns KDP (Potassium Dihydrogen Phosphate) which is an uniaxial anisotropic crystal, used in laser applications for its nonlinear and electrooptical properties.In the case of nanosecond laser irradiation, laser damage initiation is induced by precursor defects that efficiently absorb the laser energy, inducing a fast temperature rise and a subsequent shock wave in this material [35,50].The nature and origin of these defects is however unclear and still under investigation [57].Recent experiments [58] and theoretical works [50] suggest that non-spherical absorbing defects could be involved: planar or ellipsoid-shaped.We report then in the section below a study based on the previously described model on the influence of the polarization and shape of the defects on the absorption efficiency and temperature rise inside the KDP lattice.
The configuration of this study is the following.An ellipsoidal defect geometrically characterized by its semi-axis lengths denotes a d , b d and c d respectively along Ox, Oy and Oz lies in an infinite bulk lattice of KDP.The dielectric permittivity of this defect is chosen according to [58]: ε d = 0.08 + 0.06 i.The permittivity of the KDP lattice is set to 35.We are considering in what follows KDP as an optically isotropic material, which is not a big assumption since generally accepted ordinary/extraordinary optical index values are n o = 1.53 and ne = 1.49 at 1064 nm.Moreover, the simulations have been done using an anisotropic material with these values.The results showed negligible discrepancies in terms of absorption and temperature reached compared those obtained with an isotropic KDP of optical index 1.5318.
Still according to the experimental set up described in [58], let us consider an incident gaussian temporal pulse of width at half-height t 1 − t 0 = 6.5 ns and fluence 11 J.cm −2 at λ = 1064 nm polarized along Ox.The considered temporal range is [t min , t max ] = [0, 25] ns sampled with 1 ns steps.
As for thermal parameters of the KDP lattice, the conductivity is set to κ = 1.9 W • (m • K) −1 and its thermal diffusivity to κ/(ρ.C p ) = 6.5 10 −7 m 2 .s−1 .Thermal parameters for the typical defect lying in a KDP lattice being unknown so far, let us consider a typical metallic case of thermal conductivity 50 W • (m • K) −1 and diffusivity 10 −4 m 2 .s−1 .Fig. 6 sums up the absorptivity and maximum temperature reached in the system for the three different types of elongated ellipsoids.For each case, one of the three semi-axis is spanned from 5 to 150 nm and the two other semi-axis are considered of equal lengths such that the volume of the resulting ellipsoid V e = 4/3 π a b c remains constant and equal to the volume of a sphere of 25 nm radius.As depicted in Fig. 6a, the system exhibits a strong resonance when the field is polarized in the direction of elongation of the ellipsoid, for small values of the semi-axis a.The ellipsoidal defect is then very flat and one can note that this strong resonance occurs when the incident electric field is orthogonal to very sharp edges: this resonance leads to an absorptivity at least 4 times stronger than for the other polarization case (see Fig. 6b) and to a maximum temperature increase 8.5 times more important (11400 K in the case of Fig. 6a instead of 1300 K in the case of Fig. 6b).Note that the shape of the absorptivities presented in Fig. 6b and Fig. 6c are different only because the apparent surface differs from case a) to case b).Otherwise, the amount of absorbed incident power in W (the numerator of Eq. ( 17)) is the same.Finally, the maximum temperature increase inside all these systems occurs between 1 and 2 ns after the maximum of the pulse.This short timescale can be explain by both the metallic order of magnitude chosen for the diffusivity and the size of the defect.Finally, it is of importance to emphasize that these numerical values are of the same order of magnitude than in the experimental results described in [58,59].

Concluding remarks
The main advantage of the proposed finite element based multiphysics model is its complete arbitrariness towards the number, geometric shape, dielectric and magnetic properties, embedding environment of the scatterers.Moreover, the fact that the model takes into account both the spatial distribution of the thermal sources inside the scatterer and the thermal dependency.This is due to a rigorous treatment of the electromagnetic sources which leads to a precise determination of the spatial distribution of the source losses even inside a very lossy material.PMLs allow to truncate the computational domain.Using second order edge elements provides enough accuracy to match the value obtained using a Mie code, even using a very coarse mesh inside the PML.Finally, the distribution of Joule losses calculated using this rigorous and precise electromagnetic approach are used as a source of a classical temporal variational formulation of the heat equation.
We finally applied the proposed method to the study of laser damage of optical components used in high power applications in order to stress the importance of taking into account all the details of a defect scattering in such structures.For instance, considering a cluster of defects rather than a single one brings out new critical defect parameters compared to previous study where only the diameter of the defect was considered.Thus, two non critical defects if considered isolated could actually lead to a damage when considered together and close.Considering the typical shape of a defect is also of prime importance, since an elongated defect can lead to an temperature rise almost 10 times superior to a spherical one with identical volume.
Test structures have been conceived at Institut Fresnel to test the ability of this model to match experimental cases and are under realization.Taking into account more closely the damage initiators with this rigorous multi-physics model can allow a better understanding of the inner mechanisms of laser damage.As detailed in the introduction, this method can be applied to many other nanophotonic applications involving photo-induced thermal effects.

G
FIG. 1 Scheme and notations of the studied multi-physics problem.

alongFIG. 2
FIG.2Sketch of the various PML used to rigorously truncate the unbounded domain.

FIG. 3
FIG. 3 (a)-Absorptivity as a function of the scatterer radius, calculated with a Mie code (circle) and the presented FEM (asterisk).(b)-Convergence of the FEM calculation as a function of the mesh refinement.

FIG. 5
FIG. 5 Temporal evolution of the maximal temperature reached in the 1-(dash-dotted blue line) and 2-scatterer(s) (plain red line) cases.The right axe and the dotted black line represent the irradiance I(t) for a pulse fluence of 0.124 J.cm −2 .

FIG. 6
FIG. 6 Absorptivity and maximum temperature reached as a function of the elongation along Ox (a), Oy (b) and Oz (c).