EP3262466A1

EP3262466A1 - Method for numerical calculation of the diffraction of a structure

Info

Publication number: EP3262466A1
Application number: EP16714972.3A
Authority: EP
Inventors: Wolfgang IFF; Alexandre Tishchenko
Original assignee: Centre National de la Recherche Scientifique CNRS; Universite Jean Monnet Saint Etienne
Current assignee: Centre National de la Recherche Scientifique CNRS; Universite Jean Monnet Saint Etienne
Priority date: 2015-02-24
Filing date: 2016-02-22
Publication date: 2018-01-03
Also published as: WO2016135406A1; FR3033063A1; US20180046087A1; FR3033063B1

Abstract

The invention relates to a method for numerical calculation of the diffraction of a structure (1), wherein a numerical model defines the dielectric permittivity at each point, comprising the steps of: - defining a numerical modelling of a diffraction model of a structure, including: - splitting the numerical model of dielectric permittivity into N numerical models of dielectric permittivity; - determining the diffraction model of said layers starting from the numerical model of said layer; - numerically calculating the diffraction of the structure, including: - an initial iteration during which: - a wave propagation is performed in the direction of layers 1 to N, - and a wave propagation is performed in the direction of layers N to 1, - subsequent iterations are performed, each subsequent iteration of index k including: - a wave propagation in the direction of layers 2 to N: - a wave propagation in the direction of layers N-1 to 1.

Description

METHOD FOR DIGITIALLY CALCULATING THE DIFFRACTION OF A

STRUCTURE

The invention relates to the diffusion / diffraction of an electromagnetic wave by complex material structures, and in particular the modeling of diffusion / diffraction properties and the numerical calculation of diffusion / diffraction for heterogeneous structures and of considerable thickness compared to at the wavelength.

In many technical fields, such as photolithography or diffractive optical elements, it is essential to be able to accurately model and calculate the response of a component to an incident light beam. The modeling aims to calculate the diffraction of a structure characterized by the spatial distribution of the dielectric permittivity defined along different axes at each point of the structure.

The electromagnetic modeling of a structure is carried out in the state of the art to allow either an approximate resolution method, a statistical resolution method, or an exact resolution method.

The approximate resolution methods make it possible to obtain a relatively fast diffraction calculation result. However, such methods have insufficient accuracy to obtain a satisfactory calculation result.

The exact resolution methods of the state of the art are very slow and require a very large capacity of digital memory, which limits their application to optical elements of low complexity. The exact resolution methods are therefore most often used only as a means of spot verification of modeling carried out by an approximate resolution method.

Among the exact resolution methods known from the state of the art, mention may be made in particular of the method of calculation by means of the S diffusion matrices. For this, the structure to be modeled is decomposed into a set of N contiguous layers and M orders. of diffraction. The diffusion matrix Si of each index layer i is calculated and makes it possible to express the outgoing field amplitudes as a function of the amplitudes incident on the two faces of this layer. It can thus be noted that an outgoing amplitude of a layer constitutes an incidental amplitude for an adjacent layer.

A matrix S for the component is formed by combining all the matrices Si of the different layers.

The modeling then makes it possible to determine the outgoing field magnitudes of the structure formed from the set of layers as a function of the incident field amplitudes on this structure:

However, this combination of the matrices Si is not a simple product but a complex calculation. The calculation time is proportional to M ³ and generally linear in N.

An alternative method to the S-matrix method is known from the original Bremmer publication and the following Sluijter development for calculating the transmission and reflection of a system of superimposed uniform layers by calculating the transmission and the layer-to-layer reflection of a plane wave from one edge to the other of the multilayer followed by the same calculation in the opposite direction with memorization of the intermediate amplitudes, this round trip being repeated iteratively and the results of the successive iterations being summed until the sum converges. This method of propagation back and forth through the structure has been extended to diffracting structures composed of microstructured layers laterally, each step in a round trip calculating here the amplitude of the modes of each layer obtained by the method described in the document. DM Pai and KA Awada, "Analysis of dielectric gratings of arbitrary profiles and thicknesses", J. Opt. Soc. Am. A 8, 755-762. The validity of this method has been evaluated in the publication M. Nevière and F. Montiel, "Deep gratings: a combination of the differential theory and the multiple reflection series", Opt. Common. 108, 1-7 which concludes that this explicit summation of the iterative results only converges for structures where the modulation of refractive index or interfaces is weak and represents a perturbation of a basic structure.

EP2302360 discloses a method of modeling the diffraction properties of periodic microscopic structures and a method of calculating associated diffraction. The method expresses the problem by an integral volume form of a vector field replacing the electric field. The vector field is obtained from the electric field by a base change, so as to present a continuity to the limits of the material. Convolutions are performed on the vector field using convolution operators, according to Laurent's finite series. It is thus possible to produce matrix products by means of fast Fourier transformations. An operator of convolution and Base change is configured to transform the vector field to the desired electric field, via a base change depending on the material and geometry properties of the periodic structure.

The document 'New fast and memory-sparing method for rigorous electromagnetic analysis of 2D periodic dielectric structures', by Ms Shcherbakov and Tishchenko, published in Journal of Quantitative Spectroscopy & Radiative Transfer on pages 158-171, describes another method of exact resolution of a modeling of electromagnetic properties of a structure. This method is particularly suitable for dielectric structures having a periodicity in the plane, such as diffraction gratings.

The modeling method splits the structure into different flat layers parallel to an XY plane in a Cartesian XYZ coordinate system. Each layer has a respective thickness in the Z direction. A two-dimensional diffraction grating is modeled by a periodic variation of the dielectric permittivity of a layer according to two different directions defined by vectors of the network in the XY plane.

The resolution process transforms the wave equation into an implicit integral equation. The resolution process then converts the wave equation into reciprocal space along the X and Y axes. Because of the division of the structure into layers of the same thicknesses, the integral formulation can be expressed as equations. matrices, corresponding to sums of harmonics and expressing the stacking of layers in the form of sums. The components of the normal and tangent electric field are then separated. A block-Toplitz form of the matrices can be established without requiring matrix inversions. The resolution then consists in performing a matrix inversion of a matrix that can be expressed as products of block-diagonal matrices and Toplitz-block matrices. Matrix inversion is in practice computed by matrix multiplications, with a method of solving iterative linear equations of GMRES type, rather than by a direct calculation of matrix inversion. The block-Toplitz form makes it possible to carry out calculations by fast Fourier transformation, with a calculation time then substantially proportional to M, and a digital memory utilization substantially proportional to M also.

The calculation time obtained with this resolution method is proportional to N for simple and not very thick structures. However, the calculation time and the memory occupancy increase rapidly with the thickness of the layers of the structure to be modeled. A large number of iterations is then necessary to obtain a convergence of the iterative method, which results in a significant increase in computing time and digital memory. necessary. In addition, the resolution method is poorly suited for structures with heterogeneities of very different layers in the structure, which calls into question certain assumptions of the resolution method. Furthermore, the calculation remains demanding in the amount of digital memory used, the resolution data of the entire structure to be stored throughout the calculation. This amount of computing memory required greatly limits the thickness of the structures that can be modeled.

The document US6898537 describes a method for the numerical calculation of the diffraction of a structure. A numerical model of the structure defines the dielectric permittivity at each point of it. The method defines a numerical modeling of a diffraction pattern of the structure including splitting the model into a number N of digital models of respective superimposed planar layers. The method determines the diffraction pattern of each of the layers from the digital dielectric permittivity model of this layer. The method describes the numerical computation of the diffraction of the structure, based on the propagation of a wave across the layers in a given direction.

The document published by David Windt entitled "IMD-Software for modeling the optical properties of multilayer film" ds Omputers in Physics., Volume 12, No. 4, ^{January 1,} 1998, page 360, describes a digital calculation process features naked optics multilayer structure, of which a numerical model defines the dielectric permittivity at each point

The document EP1804126 describes a method for the numerical calculation of the diffraction of a structure, a numerical model of which defines the dielectric permittivity at each point. This method comprises the definition of a numerical modeling of a diffraction pattern of the structure. This definition of numerical modeling includes the splitting of the digital dielectric permittivity model of the structure into a number N of digital dielectric permittivity models of superimposed flat layers.

The invention aims to solve one or more of these disadvantages. The invention thus relates to a method of numerical calculation of the diffraction of a structure of which a numerical model defines the dielectric permittivity at each point, as defined in the appended claims.

Other characteristics and advantages of the invention will emerge clearly from the description which is given hereinafter, by way of indication and in no way limitative, with reference to the appended drawings, in which: FIG 1 is a sectional view of an example of structure to be modeled split into different layers;

FIG. 2 is a diagrammatic sectional view of an example of a simplified structure of the OLED type, to be modeled, divided into different layers;

FIG. 3 diagrammatically represents a system configured to model and numerically calculate the diffraction of a structure;

FIGS. 4 and 6 schematically represent various parameters calculated during a first phase of a calculation method;

FIGS. 5 and 7 schematically represent various parameters calculated during a subsequent phase of the calculation method;

FIGS. 8 and 9 schematically represent various parameters calculated during a first iteration of a calculation method according to another variant;

FIGS. 10 and 11 represent schematically different parameters calculated during subsequent iterations of a calculation method according to this other variant;

FIGS. 12 and 13 schematically represent various parameters calculated during a first iteration of a calculation method according to yet another variant;

FIG. 14 illustrates a logic diagram of an example of a method for the numerical calculation of the diffraction of a structure.

Figure 1 is a sectional view of an example of a structure 1 to be modeled. The structure 1 here has a parallelepiped structure with a flat upper face and a lower face on which electromagnetic waves can be incident. The structure 1 here comprises different zones, for example made of different materials having distinct dielectric permittivity distributions.

A numerical model of the dielectric permittivity of structure 1 at each point (for a given resolution of the numerical model) is known. The dielectric permittivity of structure 1 is therefore determined at each of its points, or determinable by a law (distribution, for example) in each of its points.

A first aspect of the invention aims to define a numerical model of the diffraction of structure 1 for the application of incident electromagnetic waves on the upper face and / or the lower face. The invention aims in particular to define such a numerical model to allow a diffraction calculation whose amount of digital memory used is proportional to a number M of diffraction orders considered. The invention may be implemented by a digital processing system 2 illustrated in FIG. 3. Such a digital processing system 2 may for example include a computing device 21 (for example a server equipped with an operating system and appropriate computing applications), a storage device 22 of the digital dielectric permittivity model of the structure 1, and a storage device 23 for calculation results. The storage device 23 may, for example, store digital models (detailed below) of diffraction of different layers of the structure 1 or the digital diffraction model of the whole of the structure 1.

An example of a numerical modeling method according to the first aspect of the invention may be the following. The numerical model of dielectric permittivity of structure 1 is first split into a number N of digital dielectric permittivity models. Each of these digital dielectric permittivity models corresponds to a respective flat layer of the structure 1, the different layers being superimposed in a direction perpendicular to the upper and lower faces of the structure 1. Each of these layers is subsequently identified by its index i, the index i increasing between the lower face and the upper face of the structure 1 with values between 1 and N, as illustrated in FIG.

For the sake of simplification, the different layers here have the same thickness. The invention could, however, also be applied to a splitting of the structure 1 into layers of different thicknesses. For example, Figure 2 is a sectional view of another example of a structure 1 to be modeled. Figure 2 corresponds to a simplified diagram of an OLED type structure. The different layers of the digital model are here chosen according to the functions of the different layers of the structure 1, these different layers typically having different thicknesses. The numerical models of the different layers of the structure 1 may therefore for example include a digital model of an air layer 11, a digital model of a glass layer 12, a digital model of a diffusion layer 13, a digital model of a transparent electrode 14, a digital model of a polymer layer 15, a digital model of an emitting layer 16, a digital model of a polymer layer 17 and a numerical model of a metal electrode 18.

The number of layers N is chosen so that each layer is sufficiently fine so that a diffraction pattern of each layer can be determined from its numerical model of dielectric permittivity, and so that the diffraction pattern of each layer can be calculated with a numeric memory occupancy less than or equal to K * M. advantageously, this number of layers N is chosen so that each layer is sufficiently fine so that a diffraction pattern of each layer can be calculated in a time less than or equal to K * M * log (M). K is a factor independent of M, typically a constant.

The determination of the diffraction model of each layer can for example be implemented by the method called GSM (for Generalized Source Method in English) thereafter, and described in the document 'New fast and memory-sparing method for rigorous electromagnetic analysis of 2D periodic dielectric structures, Shcherbakov and Tishchenko, published in Journal of Quantitative Spectroscopy & Radiative Transfer on pages 158-171.

The diffraction model of each layer i can for example be noted as a linear operator ϋ such that:

with f and b _M the amplitudes of the M orders of diffraction of the respective incident waves on the faces of the layer i, and f _t and b _t the amplitudes of the M orders diffracted by the layer i. The operator U can be obtained by the generalized source method (GSM) described in AA Shcherbakov, AV Tishchenko, "New fast and memory-sparing method for rigorous electromagnetic analysis of 2D periodic dielectric structures", J. Quant. Spectrosc. Rad. Transfer 1 13, 158-171 (2012) or by an analytical formulation for very thin layers with respect to the wavelength described in AV Tishchenko, "Analytical solutions of 2D grating diffraction: GSM versus Rayleigh hypothesis," Proc. SPIE 5249 p. 683-694 (2004) or may even be the S matrix in the case of single layers where S is of the Töplitz or diagonal block type. According to a second aspect of the invention, a diffraction calculation of at least one incident wave is made from the diffraction patterns of the N layers. Thus, after the determination of the diffraction patterns of the N layers, diffraction calculations can be made for different incident waves from these diffraction patterns.

According to a first diffraction calculation variant, a propagation calculation is made in parallel by applying an incident wave on the layer 1 and a propagation calculation by applying an incident wave on the N layer. digital memory can be further optimized that with this first variant (here a memory occupancy proportional to 2M ^* (N + 1)), it nevertheless proves to be particularly suitable for being implemented by parallel computing means, for example processor systems or multiple graphics cards.

For the propagation calculation based on the application of an incident wave on layer 1, an initial iteration (or iteration of index 0) of the process is illustrated with reference to FIG. 4. M incident diffraction orders whose amplitudes are contained in the vector / ₀ ° are applied to the diffraction pattern of the layer 1, on its outer face. The M transmitted diffraction orders whose amplitudes are contained in the vector are calculated with this diffraction model, and the M reflected orders b ° are calculated and stored.

Then, for each index layer i between 2 and N, the M transmitted orders calculated in the diffraction model of the index layer i are applied. We calculate the M orders transmitted f °, and we calculate and memorize the M orders reflected b °. For the layer of index N, the M transmitted orders f ° are also memorized. The memorized elements are illustrated inside dashed circles in Figure 4.

During this initial iteration, the operator U is applied in practice in the following manner, in the absence of incidence b _N ° _{+ 1} :

For the propagation computation based on the application of the M incident orders on the N layer, an initial iteration of the process is illustrated with reference to FIG.

Incident orders ° ₊₁ are applied to the diffraction model of the N layer on its outer face. The transmitted commands c _N ° are calculated with this diffraction model, and the reflex orders g _N ° are calculated and stored.

Then, for each layer of index i between N-1 and 1, the transmitted orders calculated c ₊₁ are applied to the diffraction model of the index layer i. The orders transmitted c ° are calculated, and the reflected orders g ° are calculated and stored. For the index layer 1, the transmitted commands c ° are also memorized. The memorized elements are illustrated inside dashed circles in Figure 6. During this initial iteration, the operator U is practically applied in the following way, in the absence of incidence:

The output amplitudes of the initial iteration are calculated by adding the amplitudes of the incident orders on an external face and the amplitudes of the orders reflected by this same face:

A vector of basic amplitudes is shown, forming a basic solution containing the readings inside the structure:

with all i between 1 and N-1.

This initial iteration requires a computation time proportional to the number N of layers and requires a memory occupation proportional to 2M ^* (N + 1).

This basic solution v ° does not correspond to the solution retained later. The solution chosen is determined following further iterations of the diffraction calculation detailed below.

The vector v of the amplitudes of the final solution containing the amplitudes of all the diffraction orders of all the layers is determined from the basic solution v ° by successive iterations. Each iteration is here followed by a convergence test, to determine if the solution vector of the iteration is a final solution.

Each iteration can be expressed by means of an operator P in the form ν * = Ρ.ν ", with k the index of the iteration.In principle, the solution sought adds the set of vectors of the iterations. v searched vector solution is then the solution of the implicit equation v = v ° + P.

This implicit equation is true for any linear electromagnetic system and results from the application of the Ewald-Oseen theorem described in the document Born, Max; Wolf, Emil (1999), Principles of Optics (7th ed.), Cambridge: Cambridge University Press, § 2.4. During each iteration of index k> 0, upward propagation is carried out between the layers 2 and N of each vector g * ^"1 added to each vector, and a downward propagation of each vector is carried out. added to each vector c ₊₁ , between the layers N-1 and 1.

Thus, for the upward propagation, the operator U is practically applied in the following manner, as illustrated in FIG. 5:

with f = 0.

We then memorize and.

For downward propagation, the operator U is practically applied in the following manner, as illustrated in FIG.

with c _N ^k = 0.

We then memorize and . Each iteration requires a computation time proportional to the number N of layers and requires a memory occupation proportional to 2M * (N + 1).

The final vector solution v is the solution of the implicit equation v = v ° + P.Y. A convergence test is carried out at the end of each iteration. The convergence test is part of a method of solving a system of linear equations defined by the implicit equation above.

The resolution of the implicit equation is for example based on the GMRES method, usually used to iteratively obtain a numerical solution of a system of linear equations without matrix inversion. The GMRES method is notably described in the document 'Iterative Methods for Sparse Linear Systems' by M. Saad, second edition of 'Society for Industrial and Applied Mathematics', published in 2003 (ISBN 978-0-89871-534-7). The GMRES method usually seeks to solve the system of linear equations of the type:

A.x = b

Matrix A is assumed to be invertible and of size (m x m). Moreover, we assume that b is normed, i.e., || b || = 1, || - || representing here the Euclidean norm.

The n-th space of Krylov for this problem is defined as:

K _n = {Vect} {b, Ab, A ² .b, A < ⁿ -> .b}

Where Vect corresponds to the generated vector subspace.

The GMRES method gives an approximation of the exact solution of Ax = b by the vector x _n GK _n which minimizes the norm of the residue: || Ax _n - b ||.

To guarantee the linearly independent character of the vectors b, Ab, A ^n_1 .b, we use the Arnoldi method to find orthonormal vectors qi, q2,. . . , q _n which constitute a base of K _n . Thus, the vector x _n GK _n can be written Xn = Q _n y _n with y _n GR _n , and Q _n a size matrix (mxn) formed of

The Arnoldi method also produces an upper Hessenberg matrix H _n of size (n + 1) .x _n with

A.Qn = Qn + 1. H _B

Since Q _n is orthogonal, we have

|| A.Xn - b || = || H _B .y _n - β.βι ||

Where ei = (1, 0,0, ..., 0) is the first vector of the canonical base of R _n + i, and β = || bA.xo ||, with xo an initialization vector (by zero example). Thus, x _n can be found by minimizing the norm of the residue r _n = H _n . y _n - β.βι

Each iteration of the GMRES algorithm usually includes:

perform a step of the Arnoldi algorithm;

find y _n that minimizes || r _n || ;

calculate x _n = Q _n ..y _n ;

repeat as long as the residue is larger than a quantity (called tolerance) chosen arbitrarily at the beginning of the algorithm. For the application of the GMRES method to the convergence test of the invention, the operation (x-P.x) is substituted for the multiplication operation A.x mentioned above.

The vector space W _k = (v ₀ v _l ... vj, || || = 1 of Krylov is formed using ΙΝΓ ^v ° ^{Vo = 1}

The asterisk ^* after a vector indicates its transposition with the conjugation.

The application of an iterative step v _k - Pv _k makes it possible to create a new vector v _{k + 1} and to enlarge the Krylov space up to W _{k + 1} .

The normalization operation of the vectors v _k defines the Hessenberg matrix H _ik :

such that v, - Py _k = v _k ^* _{+ l} .H _k v _k - Pv _k = v; ₊₁ H _k

The Hessenberg matrix H _k is represented in the form H _k = Q _{k + l} .R _k Where R _k is a right upper matrix, ie a matrix in which the extra-diagonal terms under the diagonal are zero, and where Q _{k + l} is a rotation / reflection matrix with \\ Q _{k + l} \\ = 1

The solution x _k is sought in the space W _k as a linear superposition of the vectors \ _k taken with the coefficients y _k :

k-l

7 = 0

We show that the best choice of y for the solution x _k is:

In this case the error standard is minimal:

The error standard is calculated at each iterative algorithm step of the GMRES method. The calculated error standard is compared to a predefined threshold. Yes the calculated error standard is greater than the threshold, a new propagation iteration is calculated. If the calculated error standard is below the threshold, it is determined that the calculated solution is sufficiently close to the final solution to interrupt the propagation iterations.

Due to the criterion of choice of the number of layers N (determined so that a diffraction model of each layer can be calculated from its numerical model of dielectric permittivity with a digital memory occupancy less than or equal to K ^* M), the GMRES method is applicable even with computing devices having limited memory resources, and the maximum thickness of the structure that can be modeled is greatly increased.

In a final step, after verifying that the convergence criterion was reached for the iteration of index k:

we find the vector solution between the layers to calculate the fields in the layers

for all i between 1 and N-1.

we thus find the vector solution at the exit of the structure A convergence test based on a GMRES method has been detailed here, but other methods for solving linear equation systems can be used for a convergence test applied to an implicit equation, such as the stabilized biconjugate gradient (BCGS) method. . During the final iteration, an upward propagation is carried out between the layers 2 and N of each vector _M added to each vector f _t _ _y , and a downward propagation of each vector b _{M is} added to each vector c _M , between layers N-1 and 1. The advantage is to find the vector solution between layers.

Thus, for the rising propagation, the operator U is practically applied in the manner illustrated in FIG. with / L = 0 and with i between 2 and N.

We then memorize and f _N.

For downward propagation, the operator U is practically applied as illustrated in FIG.

with c _N = 0 and for all i between N-1 and 1.

We then memorize g _N and.

This final iteration requires a computation time proportional to the number N of layers and requires a memory occupancy proportional to 4M. Then we find the vector solution at the exit of the structure According to a modification of the algorithm, the amplitudes reflected on the external faces of the initial iteration are calculated:

f ^~ SN

b ° = b?

The base vector v ° is enlarged by also including the amplitudes f ° and c °. At each iteration of index k> 0, we add the amplitudes f ^ and c with the vector v * to find an enlarged vector. The final solution found by such an iterative method includes the amplitudes f _N and c _x . The amplitudes / and b at the output of the structure are found by addition with the reflected amplitudes of the initial iteration: This modification results in a larger memory occupancy, proportional to (N + 1) for each iteration, hence we do not need the final iteration, the total number of iterations is therefore reduced. According to a second variant, the application of a propagation in one direction is sequentially carried out, then the orders reflected in the opposite direction are applied. This second variant makes it possible to optimize the use of the digital memory with a proportional occupation to M ^* (N + 1).

During an initial iteration (or iteration of index 0) of the process illustrated with reference to FIG. 8, for example propagation of the incident diffraction orders on layer 1 is carried out.

The vector of the incident orders / ₀ ° is applied to the diffraction model of the layer 1, on its external face. The transmitted commands are calculated with this diffraction model, and the reflected orders b ° are calculated and stored.

Then, for each layer of index i between 2 and N, we apply the transmitted orders calculated in the diffraction model of the index layer i. The orders transmitted f ° are calculated, and the amplitudes of the reflected orders b ° are calculated and stored. For the index layer N, the amplitudes of the transmitted commands f ° are also memorized. The memorized elements are illustrated inside dashed circles in Figure 8.

During this initial iteration, the operator U is applied in practice in the following manner, in the absence of incidence b ° _{+ 1} :

The initial iteration is continued by the application of the incident orders ° ₊₁ to the diffraction model of the layer N, on its external face. The transmitted commands c _N ° are calculated with this diffraction model, and the orders reflected g _N ° are calculated and memorized.

Then, for each layer of index i between N-1 and 1, we add the calculated transmitted orders c ₊₁ and the previously calculated reflex orders b ° _{+ 1} in the diffraction model of the index layer i.

The amplitude of the orders transmitted c ° is calculated, and the amplitudes of the reflected orders g ° are calculated and stored. For the index layer 1, the amplitudes of the transmitted commands c ° are also memorized. The memorized elements are illustrated inside dashed circles in FIG. 9. In order to optimize the quantity of digital memory used, after calculating a vector transmitted orders c ° and a vector g _t , it is possible to release the digital memory occupied by the amplitudes of orders b ° _{+ 1} calculated previously.

During this initial iteration, the operator U is practically applied in the following way, in the absence of incidence: The output amplitudes of the initial iteration are calculated by adding between the amplitudes of the incident orders on an external face and the amplitudes of the orders reflected by this face:

We denote by v ° a base vector, forming a basic solution with the amplitudes inside the structure: v ° = fc °} with all i between 1 and N-1.

This initial iteration requires a computation time proportional to the number N of layers. Due to the release of the digital memory occupied by the computed orders b ° _{+ 1} , this initial iteration requires a memory occupation proportional to M * (N + 1).

This basic solution v ° does not correspond to the solution retained later. The chosen solution is determined following further detailed iterations thereafter.

As for the first variant, the final solution vector v is determined from the base solution v ° by successive iterations. Each iteration is followed by a convergence test as detailed above, to determine if the solution vector of the iteration is a final solution.

Each iteration can be expressed by means of an operator P in the form Y ^k = p. \ ^{K ~ l} , with k the index of the iteration. The vector solution v is the solution of the implicit equation v = v ° + P. During each iteration of index k> 0, an upward propagation is carried out between the layers 2 and N of each vector, by taking f = 0, added to each vector g * ^"1 and each vector b resulting is stored. Then, downward propagation of each vector c ₊₁ added to each vector b _{+ 1} between the N-1 and 1 layers is performed, taking c _N ^k = 0, and each resulting vector g is stored.

Thus, for upward propagation, the operator U is applied in practice as follows, as illustrated in FIG.

with / * = 0.

We then memorize .

After each calculation of b. , we release the digital memory occupied by

For the downward propagation, in practice the operator U is applied as follows, as illustrated in FIG. with ci = 0.

We then memorize g. .

After each calculation of g. , we release the digital memory occupied by

Each iteration requires a computation time proportional to the number N of layers and requires a memory occupation proportional to M ^* (N + 1).

As for the first variant, a convergence test is performed at the end of each iteration to solve the implicit equation v = v ° + P.Y. The convergence test also uses a method of solving a system of linear equations to search for the solution of this implicit equation without matrix inversion.

During the final iteration, upward propagation is carried out between the layers 2 and N of each vector g _M added to each vector f _t _ _y , and a downward propagation of each vector b _{i + 1 is} added to each vector c _{i + 1} , between the layers N-1 and 1. For upward propagation, the operator U is practically applied in the following manner, as illustrated in FIG.

with f _x - = 0.

We then memorize b _i .

After each calculation of b _i , the digital memory occupied by

For the downward propagation, in practice the operator U is applied as follows, as illustrated in FIG. with c _N = 0.

We then memorize g _l .

After each calculation of g _t, the digital memory is freed occupied by b _M ^■

This final iteration requires a computation time proportional to the number N of layers and requires a memory occupation proportional to M * (N + 1). Then we find the vector solution at the exit of the structure

According to a modification of the algorithm, the base vector v ° is enlarged by including the amplitudes f ° and c °. At each iteration of index k> 0, we add the amplitudes f ^ and c with the vector v * to find an enlarged vector. The amplitudes / and b at the output of the structure are found by addition with the reflected amplitudes of the initial iteration: This modification results in a larger memory occupancy, proportional to (N + 1) for each iteration, hence we do not need the final iteration, the total number of iterations is therefore reduced. According to a third variant, a propagation calculation is carried out in parallel by applying an incident wave on the layer 1 and a propagation calculation by applying an incident wave on the N layer. The use of the digital memory remains the same as in the first variant (a memory occupancy proportional to 2M ^* (N + 1)), the propagation computation in the layers is more parallelizable, which proves to be particularly suitable to be implemented for a large number of parallel computing means (for example at least 4, or with a number of parallel computing means at least equal to a quarter of the number of layers N), for example systems with multiple processors or graphics cards. For simplicity, we will now detail an example in which one carries out as many calculations in parallel as of layers N.

For the calculation of propagation based on the application of an incident wave on layer 1, an initial iteration (or iteration of index 0) of the process is illustrated with reference to FIG. 12.

M incident diffraction orders whose amplitudes are contained in the vector / ₀ ° are applied to the diffraction model of the layer 1, on its outer face. The M transmitted diffraction orders whose amplitudes are contained in the vector are computed with this diffraction model and stored, and the M orders b b ° are calculated and stored:

For this initial iteration, only the transmitted orders b ° and the reflex orders g _N ° are computed. For all other amplitudes, we assume orders transmitted:

b- =

for any i between 2 and N.

For the propagation calculation based on the application of the M incident orders on the N layer, an initial operation of the process is illustrated with reference to FIG. 13.

Incident orders ° ₊₁ are applied to the diffraction model of the N layer on its outer face. The orders transmitted here are calculated with this diffraction pattern and stored, and the reflected orders g _N are calculated and stored.

For all other amplitudes, we assume that orders are reflected:

0 A for all i between 1 and N-1

We denote by v ° a basic vector, forming a basic solution with the amplitudes ° and b °

for all i between 1 and N-1.

This initial iteration requires a calculation time of a layer since the calculation is parallelized and requires a memory occupancy proportional to 4M.

To find the exact solution of the diffraction problem, iterations are necessary and are detailed later.

As for the first and the second variant, the vector v of the amplitudes of the final solution is determined from the basic solution v ° by successive iterations. Each iteration is followed by a convergence test as detailed above to determine whether the solution vector of the iteration can be considered as a final solution.

Each iteration can be expressed by means of an operator P in the form v * = pv ^{k ~} with k the index of the iteration. The vector solution v is the solution of the implicit equation v = v ° + .v

During each iteration of index k> 0, a propagation of each vector * ^"1 added to each vector b _{+ 1} to a layer i is carried out Such a calculation is carried out in parallel for each layer of index i between 1 and NOT : with = 0 and 6 ^ = 0. Each iteration requires a calculation time of a layer and requires a memory occupation proportional to 2M ^* (N-1).

As for the first two variants, a convergence test is performed at the end of each iteration to solve the implicit equation v = v ° + P. \. The convergence test also uses a method of solving a system of linear equations to search for the solution of this implicit equation without matrix inversion. During the final iteration, a vector propagation g _N _ _{x is carried out} towards the N layer: and a vector propagation c _x to layer 1

Such a calculation is carried out in parallel for the two layers 1 and N with g ^ = 0 and b ^ = 0.

This iteration requires a calculation time of a layer and requires a memory occupation proportional to 2M.

The vector solution at the output of the structure is found as the sum

According to a modification of the algorithm, the base vector v ° is enlarged by also including the amplitudes f ° and c °. At each iteration of index k> 0, we add the amplitudes f ^ and c with the vector v * to find an enlarged vector. The amplitudes / and b at the output of the structure are found by addition with the reflected amplitudes of the initial iteration: This modification results in a larger memory occupancy, proportional to (N + 1) for each iteration, hence we do not need the final iteration, the total number of iterations is therefore reduced. FIG. 14 schematically represents an example of a sequence of steps implemented in a method of numerical calculation of the diffraction of a structure.

In step 301, a numerical model of the dielectric permittivity of a structure is defined at each of its points.

In step 302, the digital dielectric permittivity model of the structure is divided into N digital dielectric permittivity models, each of these digital models each corresponding to a planar layer of the structure, these planar layers being superimposed in a perpendicular direction to the upper and lower faces of the structure. The layers, corresponding to the splitting of the dielectric permittivity model of the structure, are determined so that the diffraction pattern of each of these layers can be calculated from its numerical model of dielectric permittivity in a time less than or equal to K ^* M ^* log (M).

In step 303, a diffraction pattern is determined for each of the N layers from its digital dielectric permittivity model. The diffraction model of each layer i can for example be noted as an operator that

with f _t _ _x and b _{i + l} the incident amplitudes on the faces of the layer i, and f _i and b _; the amplitudes diffracted by the layer i.

In step 304, an initial propagation iteration of the incident orders is carried out through the diffraction patterns of the layers of the structure.

During step 305, a propagation iteration of the orders diffracted and calculated during the previous iteration is carried out through the diffraction models of the layers of the structure.

In step 306, a convergence test of the last diffracted order propagation iteration is carried out using a method of solving a system of linear equations without matrix inversion and applying it to the resolution of the implicit equation sought for convergence of iterations. If the condition of the convergence test is not fulfilled, step 305 is executed again for a new iteration of propagation of the diffracted orders. At the step 307, a final iteration of propagation of the orders diffracted and calculated during the previous iterations is carried out, through the diffraction models of the layers of the structure. We deduce the result of the diffraction of the structure.

The method of numerical computation of the diffraction of a structure detailed in the preceding examples can be applied in different technical fields, a non-exhaustive list of which is presented below by way of examples.

In diffracting optics, for the design of structures with micro- and nano-structures (micro and nano-structures whose characteristic dimension is smaller than 3 to 5 times the wavelength), the diffraction properties depend on the polarization and the usual scalar methods model them erroneously. The calculation methods according to the invention make it possible to simulate large sections of such periodic or non-periodic structures, or even the complete structure. Such a method of calculation can for example be applied to a refracto-diffracting micro-lens of minimum size and minimum weight used in a scanning system or optical scanning ultra-fast.

A calculation method according to the invention can also be applied in diffusing optics. Such a calculation method can in particular be applied for the optimization of the desired distribution of mono- or poly-chromatic light by a scattering layer. Such a diffusing layer typically comprises a host medium containing microspheres or micro-polyhedra of refractive index different from that of the host medium. The diffusing layer generally illuminates a screen uniformly, desirably illuminating a volume from a localized light source or part of a light wall. By making it possible to model exactly large sections of such scattering layers, the calculation method described takes into account the multiple diffraction effects in a layer and the coherence effects.

Such a calculation method can also be used in diffusive optics for optimizing a diffusing layer intended to extract efficiently the light produced, for example, by a light-emitting diode of the Electro-Luminescent Diode type (LED in English) or an organic electroluminescent diode.

A calculation method according to the invention can also be applied in microelectronics for the optimization of different processes of photolithography, in particular when the characteristic dimension of the patterns at the level of the reticle (the mask) and / or at the level of the silicon wafer is of the order of magnitude or smaller than the projection wavelength. This is particularly the case for the technological nodes of 45, 30 nm and below the wavelengths of excimer KrF (248 nm) and ArF (193 nm) lasers. One of the processes to be performed at the level of the reticle is generally optical proximity correction (OPC or Optical Proximity Correction in English). Such a correction OPC assumes the exact calculation of the transmission of the reticle at variable optical incidence, the reticle including patterns of different depths of materials such as SiO2, chromium, MoSiO, TaO.

The usual transmission calculation methods for such structures are scalar methods with Euristic patches based on accurate local calculations. The use of these patches makes it possible to preserve the speed of computation of the scalar methods but their field of validity is limited. A calculation method according to the invention makes it possible to push back very far the boundary between the domain of structures which can be modeled accurately and that of approximately calculable structures.

Another process in microelectronics photolithography is that of modeling the latent image in a layer of photoresist deposited on a substrate. The topography and surface composition result from a number of previous technological processes, which significantly affects the incident light power distribution by reflection diffraction. The scalar methods of the state of the art are faced with great difficulties because the structure below the photoresist layer where the projected latent image is formed can be very thick and very complex. A calculation method according to the invention makes it possible to take into account the actual structure implemented.

Claims

A method of numerical calculation of the diffraction of a structure (1) of which a numerical model defines the dielectric permittivity at each point, comprising the steps of:

-Definition of a numerical modeling of a diffraction model of a structure, including:

the division of the digital dielectric permittivity model of the structure into a number N of digital dielectric permittivity models of respective flat layers superimposed in one direction and ordered according to an index i between 1 and N in said direction;

the determination of the diffraction model of each of said layers from the digital dielectric permittivity model of this layer, the number N of layers being determined so that the digital memory occupancy for the diffraction calculation by the diffraction model of each layer is less than K ^* M, with M a number of diffraction orders of this layer and K a factor independent of M;

-digital calculation of the diffraction of the structure, including:

an initial iteration during which one realizes:

a wave propagation in the direction of layers 1 to N, including:

the calculation of a reflected wave and of a wave transmitted by application of an incident wave to the diffraction model of the index layer 1;

for a layer of index i comprised between 2 and N, calculating a reflected wave and a wave transmitted by applying an incident wave to the diffraction model of the layer of index i, this incident wave being transmitted wave calculated for the index layer i-1;

memorizing the reflected waves calculated for each of the layers of index i and the transmitted wave calculated for the layer of index N;

and a wave propagation in the direction of the N-to-1 layers, including:

a reflected wave and a transmitted wave are calculated by applying an incident wave to the diffraction model of the N index layer;

for a layer of index i between N-1 and 1, a reflected wave and a transmitted wave are calculated by applying an incident wave to the diffraction model of the index layer i, this incident wave including the transmitted wave calculated for the index layer i + 1; memorizing the reflected waves calculated for each of the layers of index i and the transmitted wave calculated for the layer of index N;

subsequent iterations, each subsequent iteration of index k including:

a wave propagation in the direction of layers 2 to N, including:

for a layer of index i lying between 2 and N, calculating a reflected wave and a wave transmitted by applying an incident wave to the diffraction model of the layer of index i, this incident wave including transmitted wave calculated for the index layer i-1 during this iteration and including the calculated reflected wave for the index layer i-1 during the last propagation in the opposite direction;

a wave propagation in the direction of the N-1 layers towards 1, including:

for a layer of index i lying between N-1 and 1, calculating a reflected wave and a wave transmitted by applying an incident wave to the diffraction model of the layer of index i, this incident wave including the transmitted wave calculated for the index layer i + 1 during this iteration and including the reflected wave calculated for the index layer i + 1 during the last propagation in the opposite direction;

memorizing the reflected waves calculated for each of the layers of index i and the transmitted wave calculated for the layer of index 1;

a convergence test of the solution Vk formed of the reflected and transmitted waves calculated during the iteration of index k, by applying an iterative method of solving systems of linear equations to the resolution of the implicit equation of convergence Vk = P.Vk, with P an operator corresponding to the transformation of the solution of any solution Vk-i into a solution Vk at the iteration k.

2. Numerical calculation method according to claim 1, wherein said iterative method of solving systems of linear equations is a GMRES type method.

The numerical calculation method according to claim 1, wherein said iterative method of solving systems of linear equations is a stabilized biconjugate gradient method.

4. A numerical calculation method according to any one of the preceding claims, wherein said determination of the diffraction pattern of each of said layer is implemented by a GSM method.

The numerical calculation method according to any one of the preceding claims, wherein N is at least 5.

6. Numerical calculation method according to any one of the preceding claims, wherein the calculation of propagation in the direction of the layers 2 to N and the propagation calculation in the direction of the layers N-1 to 1 are carried out in parallel.

7. A method of numerical calculation according to any one of claims 1 to 5, wherein the computation of propagation in the direction of the layers 2 to N and the propagation calculation in the direction of the layers N-1 to 1 are carried out sequentially.

The digital computation method according to any one of claims 1 to 5, wherein the computation of propagation comprises calculations in parallel for different layers.