JP2001156140A - Wire em simulation method - Google Patents

Abstract

PROBLEM TO BE SOLVED: To accurately predict a wire lifetime by taking in advance a facet face generated in EM of an actual wire into a surface element of a void surface. SOLUTION: A wire EM simulation method comprises the steps of a face direction determining means 257 for providing a face direction in each surface element of a void surface of a wires, a means 256 for introducing surface energy in each face direction thereof, and a means 255 for calculating the chemical potential in each face direction from this surface energy. A void shape change is analyzed in figures by an atomic flux component in a void region generated by a current, and an atom flux component by a slope of the chemical potential.

Description

DETAILED DESCRIPTION OF THE INVENTION

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a method for simulating a failure caused by electromigration that occurs when a current flows through a wiring, and particularly to a method for predicting a failure caused by the occurrence of voids.

[0002]

2. Description of the Related Art With the miniaturization of semiconductor elements forming a semiconductor integrated circuit (LSI), the width of wiring in a chip has been reduced, and the structure of crystal grains forming wiring has changed from polycrystalline to polycrystalline bamboo. Was. In the wiring having such a bamboo-like grain boundary structure, since the number of grain boundaries as a fast diffusion path is reduced, the electromigration of the wiring (hereinafter, referred to as “electromigration”)
EM). However, this wiring has slit-like voids and has a small number of wirings that are cut quickly, as described in, for example, the article described in Applied Physics Letter, Vol. 61, pp. 3121, 1992. It is known that the variation in the service life is large. That is, a phenomenon in which the shape of the void correlates with the EM lifetime has been reported.

In order to obtain a guideline for improving the EM life, it is necessary to clarify the factors governing the void shape, and a simulation is performed for that purpose.
The simulation considers the influence of the current density distribution near the void and the temperature field. At the same time, the void surface tension is taken into account.

After inputting the initial void shape, the current density and temperature distribution are obtained by the finite element method, and the movement of the void is calculated by the difference method using the surface diffusion and surface tension of the wiring metal. The above-mentioned current heat transfer and surface diffusion calculations are repeated to simulate the growth of voids. Hereinafter, the technical contents of the above documents will be briefly described.

The shape deformation uses surface diffusion in consideration of surface tension and atomic flow (EM term) by electric current. The flux Js on the surface is expressed by equation (1). Then, the position change speed vn of the void surface in the normal direction is expressed by equation (2).

[0006]

Where δs is the thickness of the surface diffusion layer, Ds is the surface diffusion coefficient, Ω is the atomic volume, k is the Boltzmann constant, T is the temperature of the system, Z ^* is the effective charge, ρ is the resistivity, and js is the current density. , Γ is surface tension, κ is curvature, and s is coordinates along the surface.

FIG. 7 shows a void and a mesh.
Explain how the void surface moves. In FIG. 7, flux I of adjacent elements a and b by EM
se (a) and Ise (b) are each expressed by equation (3).

Here, js (a) and js (b) are calculated by the finite element method using the current densities at the elements a and b. From the mass conservation, the moving speed vne in the normal direction of the void surface of the node 2 is derived as in equation (4).

Here, sa and sb are nodes 1 and 2 and node 2
And a distance of 3. The direction of vne is perpendicular to the line connecting nodes 1 and 3. The distance une at which the node 2 travels at the speed vne during the time step Δt is represented by Expression (5).

Next, the flux due to the surface tension is calculated. First, the curvature of each node is calculated using the positions of adjacent nodes. In this case, a circle passing through three points is obtained, and is obtained as the reciprocal of the radius. The curvature is positive when the void surface is concave. The flux per unit film thickness using the curvature as the driving force can be calculated for element a and element b as in equation (6). Using these fluxes, the movement distance of the surface at the node 2 (during the time step Δt) can be calculated as in equation (7).

[0012]

Here, it is necessary to make certain assumptions at the nodes at both ends of the void surface. That is, since these nodes are located on the side wall of the wiring, they move along the side wall. In the equilibrium state, the angle between the void surface and the wiring side wall is assumed to be 90 degrees.

In the above, the processing has been described for the growth of voids without considering facets. However, by continuously modifying the processing, the conventional technique reproduces the growth of voids in a facet shape (a shape having a contour of a facet surface). are doing.

In order to reproduce the shape evolution of facet-like voids observed in the experiment, the angle dependence was introduced into the surface diffusion coefficient Ds as shown in equation (8).

Here, Ds and min are the diffusion coefficient on the surface with the slowest surface diffusion, and the dimensionless numbers A and z are
The degree of anisotropy and the symmetry (number) of the crystal. Angle φ
Is the angle between the tangent to the void surface and the longitudinal direction of the wiring, and the angle φ0 is the angle representing the relative arrangement of the wiring and the crystal.

1980 Surface Science, 97
Vol. 73 (Surface Science 97
According to the paper of (1980) 73-87), it is reported that the surface diffusion coefficient differs by two or more digits depending on the surface. Based on this report, the relationship of the magnitude of the diffusion coefficient of the plane orientation was calculated as Ds (110)> Ds (111)> Ds (100).
And From equation (8), the azimuth at which the diffusion coefficient becomes maximum is z
And it is considered from the above assumption that there is a (110) plane along the direction of the maximum diffusion coefficient.

FIG. 8 shows the result of a calculation incorporating the growth of facet-like voids in this manner. FIG. 8 is a simulation result showing a temporal change of the void shape. Here, in addition to the above assumptions, a constant growth rate (B) is assumed on the void surface. The solid arrow in the figure indicates the direction of the surface having the fastest diffusion coefficient. Each condition in FIG. 8 is as follows.

FIG. 8A: A = 10, z = 3, φ0 =
0, B = 0.36 μm ² / h FIG. 8 (b); A = 10, z = 1, φ0 = 0, B = 0.3
6 μm ² / h FIG. 8C; A = 10, z = 1, φ0 = 0, B = 0.0
72 μm ² / h FIG. 8D; A = 10, z = 1, φ0 = −π / 4, B =
0.072 μm ² / h FIG. 8A shows that at z = 3, a crystal grain having three {110} plane groups matches that when the paper surface matches (111). The void grows from the first semicircle to a shape close to a hexagon (half). However, each side is assumed (110)
It is slightly off the plane.

FIGS. 8B, 8C and 8D show the case where z = 1.
This is an arrangement of crystal grains in which only one (011 bar) plane is perpendicular to the paper surface and the paper surface is the (110) plane. In these cases, the simulation results appear to be along the (11 bar 1) plane.

[0021]

However, in the prior art as described above, a plane deviated from the assumed plane orientation appears, and the growth of voids cannot be accurately reproduced. As a result, the life prediction becomes inaccurate. Further, the calculation in the case where there is a grain boundary on the void surface cannot reflect the properties of the grain boundary. Therefore, the prediction of the void shape due to the properties of the grain boundaries becomes inaccurate, and the prediction of the life becomes inaccurate.

This is semi-quantitative since the orientation of the crystal in the shape change of the void surface is considered only in the order of the size, depending on the plane orientation dependence of the diffusion coefficient Ds. This is because the dependence cannot be reproduced quantitatively.

The angle dependence of the diffusion coefficient Ds is given by
It is taken as the square of the cos of z (φ + φ0) as in equation (8). Therefore, the angle at which Ds becomes the maximum is (φ = nπ / z-φ0, n is an arbitrary integer) around D
Since s takes a value close to the maximum value, the angle of the growing facet surface tends to shift. Furthermore, since only the direction of the largest diffusion coefficient Ds is meaningful, other plane orientations are not physically defined. Therefore, even if the shape is reproduced, the other plane orientations are only estimates and ambiguities arise in the interpretation of the plane orientations. In addition, the surface energy is not taken into account, and the mechanism does not include a mechanism in which an original low index surface is preferentially generated.

An object of the present invention is to solve the above problems,
An object of the present invention is to provide a simulation method capable of easily performing a numerical analysis of a shape change of a void formed by a facet surface, and to improve the wiring reliability of an LSI.

[0025]

For this purpose, according to the wiring EM simulation method of the present invention, a plane orientation is provided for each surface element on the void surface of the wiring, and a chemical potential is calculated from the surface energy of the plane direction. The flux of the atoms constituting the wiring is obtained from the gradient of the chemical potential, and the time change of the void surface is numerically analyzed. Here, among the void surfaces, the chemical potential is determined such that elements on a surface having a small surface energy are large and elements on a surface having a large surface energy are small.

According to the wiring EM simulation method of the present invention, in analyzing the time change of the void surface, an electric field analyzing means in the wiring, a temperature field analyzing means, and a diffusion and a drift of atoms constituting the wiring are analyzed. Means.

According to the wiring EM simulation method of the present invention, a step of determining a plane direction to have a plane direction for each surface element on the void surface of the wiring, a step of deriving a surface energy for each plane direction, Calculating a chemical potential for each plane orientation from energy.

Here, in the wiring EM simulation method of the present invention, a step of determining an initial void surface shape, a step of generating a mesh in the void, a step of calculating a current density distribution in the void region, Calculating the temperature distribution in the void region.

The wiring EM simulation method according to the present invention includes a step of obtaining an atomic flux from a current in the wiring, and a step of determining whether or not the void is exposed on the wiring surface. When the voids are exposed on the wiring surface by flowing the gas, the atomic flux in the void region generated by the current is determined, and the change in the void shape is numerically analyzed with the flux component due to the chemical potential gradient.

As described above, according to the present invention, the EM
Since the facet surface generated in step (1) is previously incorporated into the surface element of the void surface, the wiring life can be accurately predicted.

[0031]

DESCRIPTION OF THE PREFERRED EMBODIMENTS The present invention predicts a failure caused by voids in a method for simulating a failure caused by electromigration that occurs when a current flows through a wiring. In particular, in a wiring made of a polycrystalline Al alloy, voids often consist of facet planes having a low index plane orientation, and the present invention provides a method for simulating a shape change of a void holding a facet plane. .

A plane orientation is provided for each surface element of the void surface, and the surface energy and the surface diffusion coefficient obtained for each plane orientation are used. A void surface deformation simulation method characterized by incorporating a flux of a chemical potential gradient due to surface energy so as to increase a surface element having a small surface energy and reduce a surface element having a large surface energy.

An embodiment of the present invention will be described below with reference to the drawings. FIG. 1 shows an outline of an apparatus used in the simulation of the present invention and its means.

As shown in FIG. 1, an input device 1 such as a keyboard and a data processing device 2 operating under program control
And a storage device 3 for storing information, and an output device 4 such as a display device or a printing device.

The storage device 3 comprises a shape storage unit 31, a mesh information storage unit 32, a temperature distribution storage unit 33, a potential distribution storage unit 34, and an atomic flow rate storage unit 35. The data processing device 2 includes an initial shape input unit 21, a mesh generation unit 22, an electric field analysis unit 23, a temperature field analysis unit 24, a diffusion drift analysis unit 25, and a shape development unit 26.

Here, the initial shape input means 21 inputs, from the input device 1, the shape of the wiring and the surrounding insulating film, physical properties corresponding to the material, boundary conditions for performing the wiring failure simulation, and the like. Then, it is stored in the shape / boundary condition storage unit 31. The mesh generation unit 22 extracts shape data for each of the wiring and insulating film in the shape / boundary condition storage unit 31, generates a mesh, and stores the mesh information in the mesh information storage unit 32. The electric field analysis unit 23 refers to the mesh information of the wiring from the mesh information storage unit 32 and the boundary condition from the shape / boundary condition storage unit 31 and uses the analysis means such as the finite element method or the finite volume method to calculate the Poisson equation. The potential at each node is stored in the potential distribution storage unit 34. The temperature field analysis unit 24 refers to the mesh information of the wiring and the insulating film from the mesh information storage unit 32, and
Referring to the boundary conditions related to the thermal analysis from the boundary condition storage unit 31, solving the heat conduction equation using analysis means such as the finite element method or the finite volume method, calculating the temperature at each node,
It is stored in the temperature distribution storage unit 33.

Then, the diffusion drift analysis means 25 refers to the mesh information of the grain boundary, the interface, and the void surface from the mesh information storage unit 32, and calculates the diffusion coefficient and the diffusion coefficient at each node on the grain boundary, the interface, and the void surface. The driving force based on the EM, the concentration, and the surface energy is calculated, the flux based on each driving force is obtained, and the flux is stored in the atomic velocity storage unit 35.

The shape developing means 26 includes an atomic velocity storage unit 35
, The balance of the flux at each surface element is calculated, the amount of movement of the surface element is determined, the position of the node forming the surface is changed, and the shape / boundary condition storage unit 31 is stored. Further, the shape is output to the output device 4.

Next, the structures of the diffusion drift calculation means 25 and the shape development means 26 will be described in detail with reference to FIG.

As shown in FIG. 2, the atomic velocity storage unit 35
In more detail, the surface diffusion coefficient storage unit 331, the surface energy storage unit 332, the crystal grain data storage unit 333,
It comprises a grain boundary diffusion flux storage unit 334, an EM flux storage unit 335, and a surface flux storage unit 336.

The diffusion drift calculator 25 includes a grain boundary diffusion flux calculator 251 and an EM flux calculator 25.
2. It comprises a surface flux calculating means 253, a surface diffusion coefficient searching means 254, a chemical potential calculating means 255, a surface energy searching means 256, and a plane orientation determining means 257.

The grain boundary diffusion flux calculation means 251 takes out the temperature at each node on the grain boundary with reference to the temperature distribution storage unit 33 and calculates the grain boundary diffusion coefficient. Further, a concentration gradient is calculated from a grain boundary diffusion coefficient and an atomic concentration at an adjacent node on the grain boundary, a flux of the grain boundary diffusion is calculated, and stored in the grain boundary diffusion flux storage unit 334.

The EM flux calculation means 252 takes out the temperature at each node on the grain boundary with reference to the temperature distribution storage section 33 and calculates the grain boundary diffusion coefficient. Further, the current density is calculated from the grain boundary diffusion coefficient and the temperature, and the potential at an adjacent node on the grain boundary, and the EM flux is calculated and stored in the EM flux storage unit.

The surface flux calculation means 253 extracts the gradient of the chemical potential along the surface from the chemical potential calculation means 255 and the surface diffusion coefficient from the surface diffusion coefficient search means 254, calculates the surface flux, and calculates the surface flux. To be stored.

Then, the chemical potential calculation means 255
Determines the orientation of each surface element along the void surface using the plane orientation determining means 257, extracts the surface energy according to the plane orientation information from the surface energy searching means 256, obtains the chemical potential, and determines the chemical potential along the surface. The calculated chemical potential gradient is calculated and passed to the surface flux calculating means 253.

Next, the procedure of the wiring fault simulation of the present invention will be described with reference to FIG. FIG. 3 is a flowchart for explaining an embodiment of the wiring failure simulation method according to the present invention. That is, after the initial shape is input, the current density and temperature distribution are obtained by the finite element method, and the wiring based on electromigration is used to calculate the grain boundary and the surface diffusion amount of the atoms on the crystal grain boundaries and the void surface based on the current density and temperature distribution It is a figure showing a flow chart of failure simulation. Hereinafter, each step will be described in order.

Step A1: First, the initial structure according to the wiring shape is input from the input device 1 by the initial shape input means 21, and physical values and boundary conditions corresponding to the wiring shape and material are stored in the shape / boundary condition storage unit 31. To be stored.

Step A2: Next, a mesh is generated for the initial structure by the mesh generation means 22 and stored in the mesh information storage unit 32.

Step A3: Next, the electric field analyzing means 23
Thus, the Poisson equation is solved using analysis means such as the finite element method or the finite volume method, and the potential at each node is stored in the potential distribution storage unit 34.

Step A4: Next, the temperature analysis means 24
By using the same mesh as that for calculating the current density distribution, the temperature distribution is obtained using the Joule heat generated by the obtained current as a heat source, and stored in the temperature distribution storage unit 33.

Step A5: Next, the atomic flux at each node is obtained by the diffusion drift analysis means 25.
The flux in each element is stored in the atomic velocity storage unit 3
5 is stored. Further, the shape evolving means 26 changes the atomic density in the grain boundary at the grain boundary, moves the surface element in consideration of the entrance and exit of the atomic flux on the surface, and advances the shape deformation of the void.

Step A6: Next, the number of repetitions is updated, and it is determined in step A7 whether the required number of repetitions has been reached. If not, step A2
The process returns to the generation of the mesh, and the process is ended if the mesh has been reached.

Next, the deformation of the void shape will be described with reference to the detailed flowchart of the deformation of the void shape shown in FIG.

Step B1: First, potential gradients are obtained from potentials at both ends of all surface elements or grain boundary elements by the EM flux calculating means 252 shown in FIG. 2, and the effective charge Z ^* , the charge e of the electron, the electric field The EM flux J _E on the surface or on the grain boundary is determined from (= potential gradient) E and the diffusion coefficient corresponding to the temperature of the node. Thus, J _E = D (T) / kT · Z ^* · eE.

Step B2: Next, since the subsequent processing differs depending on the surface or the grain boundary, a judgment is made. The processing in the subsequent steps is performed on the entire surface element or the entire grain boundary element.

Step B3: In the case of the surface, the surface diffusion coefficient is determined in advance from the interatomic potential using the molecular dynamics method (MD) based on the EM flux and taking into account the plane orientation, and is tabulated. Then, the surface diffusion coefficient is obtained by the surface diffusion coefficient search means 254, and the surface elements are moved according to the increase and decrease of the atoms by the EM flux.

Step B4: Next, the gradient of the chemical potential μ is obtained by differentiation along the void surface s by the chemical potential calculation means 255, and is stored in the surface flux storage unit 336 by the surface flux calculation means 253. Here, the surface flux due to the chemical potential gradient is represented by Js = −Ds / ΩkT · dμ / ds.

Step B5: Next, for each surface, the surface flux is calculated in proportion to the difference between the flux flowing into the surface element and the flux flowing out of the surface element according to the surface flux, and the amount of movement of the surface in the normal direction is calculated. Is determined, the surface is moved, and the processing on the surface is completed.

Step B6: In the case of a grain boundary, first, the atomic concentration in the grain boundary element is changed according to the EM flux calculated in step B1.

Step B7: Next, the flux due to the concentration gradient generated in step B6 is calculated by the grain boundary diffusion flux calculation means 251 and stored in the grain boundary diffusion flux storage section 334.

Step B8: Finally, the density at the grain boundaries is calculated based on the grain boundary diffusion flux at each grain boundary, and the process exits.

Next, a method of calculating a chemical potential for providing a driving force for surface diffusion will be described with reference to FIG. Chemical potential is the change in the amount of unit substance,
Since it is an energy change, when two facet surfaces are adjacent to each other, how the surface energy changes and the amount of the substance are calculated together with the surface shape deformation (movement of the facet surface).

In the surface elements a and b, the length of the surface element a is La, the length of the surface element b is Lb, and the angle formed on the material side is θ. Further, assuming that the surface energy per unit area is γa, γb, La × γa + L
b × γb. Now, suppose that the substance flows from another adjacent surface, and La moves by va and Lb moves by vb. La
The length (surface area) of La′Lb changes to Lb ′. From geometrical considerations, La ′ and Lb ′ are finally obtained as in equation (9). The area increase ΔS at that time is (1
It can be obtained from equation (0).

[0064]

As described above, the chemical potential accompanying this change can be calculated from the surface energy of the void and the increase in the area thereof. That is, in the case of the present invention, equation (9) and (1)
Equations (0) to (11) can be used.

The growth ratio between va and vb is represented by parameters α, β, and v, and va = αv and vb = βv.
Here, α ² + β ² = 1. These relationships are defined as (1
Substituting into equation (1) gives equation (12).

[0067]

Then, the chemical potential μ is obtained by limiting v to infinity in equation (12). In this way, the chemical potential μ becomes as shown in equation (13).
Note that α or β is determined so that the chemical potential μ is minimized. This is intended to locally minimize the facet surface energy. The gradient of the chemical potential μ along the surface becomes the driving force for the flow of matter along the surface.

Next, FIG. 6 shows the result of a simulation calculation of the shape deformation of the void. This is a view from the top of the wiring,
There are two crystal grains and one grain boundary between them. A void having a facet face in one crystal grain (left) was created in the initial state. The facet surface almost perpendicular to the wiring was (100), and the oblique facet surface was (111).
The current flows from left to right, the electrons flow from right to left, the current density is 2 × 10 ⁶ / cm ^{2 in} the void-free area, and the temperature is 250.
° C.

Each surface ((111), (110), (10
0)) surface energy and, of surface diffusion coefficient for each surface parameters (D _S and the activation energy Ea)
Is based on the EAM potential calculated using the parameter set described in the 1987 Journal of Materials Research, Vol. 5, page Voiles and Dow papers, 1991 Surface Science, 253
Vol. 334, published in Liu et al. Al
Table 1 shows the values for metals.

[0071]

[Table 1]

When the surface energy and the surface diffusion coefficient of another surface are required, the surface energy and the surface diffusion coefficient are determined by using the EAM potential and the molecular dynamics method, and the data is added to the table.

FIG. 6A shows the initial state, and FIG.
After h, FIG. 6C shows a simulation result of voids after 20 h. Generally, in the Al wiring system, the flow of atoms along the interface is small, and the body diffusion through the crystal grains is smaller than that and can be ignored. Therefore, atoms flow from right to left through the upper and lower interfaces. Further, the presence of the grain boundaries causes the atoms to flow to the upper interface, so that the atoms are depleted in the void portion. FIG. 6B shows the shape change after 10 hours. The atoms decrease from the entire void and the surface recedes, but the surface diffusion is fast and the (111) plane with low surface energy expands,
The (100) plane has a larger area, but hardly recedes. This is a result of the effect of the surface energy according to the present invention. FIG. 6 (c) shows the result after 20 hours. When the grain boundary and the (111) plane merge, the crystal grains are different. Surface appeared as a void surface.

[0074]

As described above, the present invention is a method for simulating a change in the shape of a void holding a facet surface, in which each surface element of the void surface has a plane orientation and is determined for each plane orientation. Use surface energy and surface diffusion coefficient. Also, the flux of the chemical potential gradient due to the surface energy is adopted so that the surface element having a small surface energy is increased and the surface element having a large surface energy is reduced.

According to the EM simulation method of the present invention, the facet growth of the void in the wiring can be accurately predicted. This is because, in the deformation of the void shape, 1) the surface diffusion coefficient and the surface energy are obtained in advance for each plane orientation, and 2) the surface energy and the surface diffusion coefficient depending on the plane orientation are incorporated in the flux.

As described above, according to the present invention, the wiring life of an actual wiring by EM can be accurately predicted.

[Brief description of the drawings]

FIG. 1 is a block diagram of an apparatus and means used in a simulation of the present invention.

FIG. 2 is a block diagram for explaining the above-described device and means in detail.

FIG. 3 is a flowchart illustrating a simulation method according to the present invention.

FIG. 4 is a flowchart illustrating a simulation method according to the present invention.

FIG. 5 is a schematic diagram for explaining a chemical potential deriving method used in the present invention.

FIG. 6 is a diagram showing a simulation result of the present invention.

FIG. 7 is a mesh diagram in a finite element method for explaining a conventional technique.

FIG. 8 is a simulation result of a conventional technique.

[Explanation of symbols]

 Reference Signs List 1 input device 2 data processing device 3 storage device 4 output device 21 initial shape input means 22 mesh generation means 23 electric field analysis means 24 temperature field analysis means 25 diffusion / drift analysis means 26 shape development means 31 shape / boundary condition storage unit 32 mesh Information storage unit 33 Temperature distribution storage unit 34 Potential distribution storage unit 35 Atomic flow velocity storage unit 251 Grain boundary diffusion flux calculation unit 252 EM flux calculation unit 253 Surface flux calculation unit 254 Surface diffusion coefficient search unit 255 Chemical potential calculation unit 256 Surface energy search Means 257 Surface orientation determining means 331 Surface diffusion coefficient storage section 332 Surface energy storage section 333 Crystal grain data storage section 334 Grain boundary flux storage section 335 EM flux storage section 336 Surface flux storage section

Claims (8)

[Claims] 1. A plane orientation is given to each surface element of a void surface of a wiring, a chemical potential is calculated from a surface energy of the plane orientation, and a flux of atoms constituting the wiring is obtained from a gradient of the chemical potential. A time change of the void surface is numerically analyzed;
M simulation method. 2. The change in surface energy when adjacent surface elements are displaced in the normal direction is divided by the change in volume before and after the displacement, and the minimum value of the change in surface energy per unit volume is defined as the chemical potential. 2. The wiring EM simulation method according to claim 1, wherein: 3. The wiring EM simulation method according to claim 1, wherein an adjustment is performed so that a value obtained by adding all surface energies of the respective surface elements is minimized. 4. An analysis of a time change of the void surface includes an electric field analyzing means in a wiring, a temperature field analyzing means, and a means for analyzing diffusion and drift of atoms constituting the wiring. The wiring EM simulation method according to claim 1, 2, or 3. 5. A step of determining a plane orientation for giving a plane orientation to each surface element of a void surface of a wiring, a step of deriving a surface energy for each plane orientation, and a chemical potential for each plane orientation based on the surface energy. 5. The wiring EM simulation method according to claim 1, further comprising the step of: 6. A step of determining an initial void surface shape, a step of generating a mesh in the void, a step of calculating a current density distribution in the void region, and a step of calculating a temperature distribution of the void region. 6. The wiring EM simulation method according to claim 4, wherein: 7. A step of obtaining a flux of atoms from a current in a wiring, and determining whether a diffusion path is diffusion of a surface of the wiring or diffusion of a grain boundary of a crystal.
7. The method according to claim 1, further comprising a step of calculating a flux of atoms proportional to a chemical potential gradient due to surface energy in the case of surface diffusion and a concentration gradient at a grain boundary in the case of grain boundary diffusion. The wiring EM simulation method according to claim 1. 8. When the void is formed on the side wall of the wiring by applying a current to the wiring, an atomic flux on the void surface due to the current is obtained, and a change in the void shape is numerically calculated using a flux component due to the chemical potential gradient. The wiring EM simulation method according to claim 1, wherein the method is performed.