Grain Boundary Resistance in Copper Interconnects from an Atomistic Model to a Neural Network

TL;DR

This paper investigates copper grain boundary resistivity using atomistic models and develops a neural network-based compact model, providing insights into orientation effects and resistivity distribution in large structures.

Contribution

It introduces a methodology combining atomistic tight binding, EAM, and NEGF for resistivity modeling, and constructs a neural network-based predictive model.

Findings

Resistivity distribution in three-grain structures is approximately normal.

Methodology validated with <5nm grain boundaries showing 6.4% deviation.

Neural network model effectively predicts grain boundary resistivity.

Abstract

Orientation effects on the resistivity of copper grain boundaries are studied systematically with two different atomistic tight binding methods. A methodology is developed to model the resistivity of grain boundaries using the Embedded Atom Model, tight binding methods and non-equilibrum Green's functions (NEGF). The methodology is validated against first principles calculations for small, ultra-thin body grain boundaries (<5nm) with 6.4% deviation in the resistivity. A statistical ensemble of 600 large, random structures with grains is studied. For structures with three grains, it is found that the distribution of resistivities is close to normal. Finally, a compact model for grain boundary resistivity is constructed based on a neural network.

Tables

Table 1. Table 1: Specific resistance for different CSL ( Σ ​ N Σ 𝑁 \Sigma N ) calculated by TB, EH and DFT.

	Specific resistance CSL ${\gamma }^R\left({10}^{-12}\Omega c{m}^2\right)$
CSL GB	${\gamma }_{DFT}$	${\gamma }_{EH}$	${\gamma }_{TB}$	Experimental	Other References
$\Sigma 3$	0.156	0.173	0.158	0.170  Lu et al. (2004)	0.202  Kim et al. (2010)
					0.155  Zhou et al. (2010)
					0.158  César et al. (2014)
					0.148  Zhang et al. (2007)
$\Sigma 5$	1.759	1.934	2.240		1.885  Kim et al. (2010)
					1.49    César et al. (2014)
$\Sigma 9$	1.82	1.72	2.14		1.75    César et al. (2014)
$\Sigma 11$	0.64	0.57	0.71		0.75    César et al. (2014)
$\Sigma 13$	2.01	1.72	2.09		2.41    César et al. (2014)

Table 2. Table 2: TB parameters for Cu following the notation on ref  Hegde  et al. ( 2014 )

Parameter Name	Value	Parameter Name	Value
$VBO$	3.6540	$I\_D\_D\_\Delta$	-0.08
$E\_S$	-4.5236	$q\_D\_D\_\sigma$	4.8355
$E\_Px$	-0.1458	$q\_D\_D\_\Pi$	4.7528
$E\_Py$	-0.1458	$q\_D\_D\_\Delta$	4.2950
$E\_Pz$	-0.1458	$I\_S\_S\_\sigma$	0.4
$E\_Dxy$	-4.3034	$I\_S\_P\_\sigma$	0.4457
$E\_Dyz$	-4.3034	$I\_S\_D\_\sigma$	-0.36819
$E\_Dxz$	-4.3034	$I\_P\_P\_\sigma$	1.5605
$E\_D{z}^2$	-4.3034	$I\_P\_D\_\sigma$	-0.2532
$E\_D{x}^2\_y^2$	-4.3034	$I\_P\_P\_\Pi$	-0.1348
$V\_S\_S\_\sigma$	-0.9588	$I\_P\_D\_\Pi$	0.0135
$V\_S\_P\_\sigma$	1.4063	$q\_S\_S\_\sigma$	2.20333
$V\_S\_D\_\sigma$	-0.1841	$q\_S\_P\_\sigma$	2.6554
$V\_P\_P\_\sigma$	1.4025	$q\_S\_D\_\sigma$	0.2495
$V\_P\_P\_\Pi$	-0.5730	$q\_P\_P\_\sigma$	1.5905
$V\_P\_D\_\sigma$	-0.4607	$q\_P\_P\_\Pi$	2.9059
$V\_P\_D\_\Pi$	0.3373	$q\_P\_D\_\sigma$	3.8124
$V\_D\_D\_\sigma$	-0.3709	$q\_P\_D\_\Pi$	3.9330
$V\_D\_D\_\Pi$	0.2760	$p\_S\_S\_\sigma$	1.3692
$V\_D\_D\_\Delta$	-0.0735	$p\_S\_P\_\sigma$	2.8794
$I\_D\_D\_\sigma$	-0.15	$p\_S\_D\_\sigma$	3.94296
$I\_D\_D\_\Pi$	-0.2498	$p\_P\_P\_\sigma$	5.5023
$p\_P\_P\_\Pi$	0.536231	$p\_P\_D\_\sigma$	-1
$p\_P\_D\_\Pi$	-1	$p\_D\_D\_\sigma$	-0.83723
$p\_D\_D\_\Pi$	0.66507	$p\_D\_D\_\Delta$	4.8475
$R_0\_inter$	0.25526	$R_0\_intra$	0.25526

Equations

$$G=\frac{1}{R}=\frac{2e^2}{h}\int T(E_f,\mathbf{k})d^{2}k,$$

$${\gamma }^E=\frac{E_{slab}-N{E}_0}{A},$$

$$\rho (\alpha ,\beta ,\gamma )=$$

$$-6.15{\beta }^2+13.41\gamma -3.91\beta \gamma -5.18{\gamma }^2$$

Full text

I INTRODUCTION

Due to the aggressive downscaling of logic devices, interconnects have reached the nanoscale, making quantum effects important. According to the roadmap provided by ITRS, interconnects are expected to reach sizes of 10 to 30 nm in the next decade ITR (2014). Previous work by Graham et al. Graham et al. (2010) demonstrates that surface scattering and grain boundary (GB) scattering play major roles in the resistance of structures smaller than 50 nm. Earlier works based on semi-empirical parameters have described polycrystalline films and surface scattering Fuchs and Mott (1938); Mayadas (1969) for macroscopic systems, but the fact that those models require fitting parameters for each experimental setup limits the scope of their applications. The ultra-scaled interconnects suggested by the roadmap require better descriptions of orientation and confinement effects to correctly model scattering in wires. Recently, first-principles calculations have been used to describe the resistance of a single grain boundary by making use of non-equilibrium Green’s function with Density Functional Theory (DFT-NEGF) formalism César et al. (2014). The results demonstrate a strong correlation between resistance and the geometry of the grain boundary , and show agreement with both experimental Kim et al. (2010) and other theoretical work Feldman et al. (2010); Zhou et al. (2010); Zhang et al. (2007). However, the studied structures are limited to relatively small sizes containing single grain boundaries and less than a few hundred atoms because of the computational burden required to perform DFT-NEGF calculations.

The purpose of this manuscript is to introduce an atomistic model that describes the resistivity due to grain boundary effects for realistic copper interconnects as projected by the ITRS roadmap ITR (2014) without depending on any phenomenological parameter. Even though the atomistic model is much faster than an ab initio method, parametric models have the advantage of easily providing a quantitative value of resistivity. Therefore, a compact model which reduces the computation time is generated by making use of a neural network that is based on large statistical sample. The rest of the manuscript has been organized as follows. Section  II presents the main characteristics of the atomistic models and benchmarks tight binding parameters against first principles calculations for a copper FCC structure. Section  III constructs single grain boundaries based on coincident site lattice (CSL) and validates their electronic properties against an ab initio method. Section IV describes grain boundary effects on copper interconnects using a system of three grains of 10 nm length simulated with an atomistic method which is benchmarked in the previous sections and quantifies the effect of misorientation. Section V proposes a compact model based on three different algorithms and finds that a neural network approach best matches the results obtained from the atomistic methods, allowing the results to be generalized to any grain boundary system configuration with a total length of 30 nm. Section VI presents a summary of this work.

II Description of Tight Binding Models

The two tight binding methods used in this study are an environmental orthogonal tight binding model (TB) Hegde et al. (2014) and a non-orthogonal tight binding method based on the Extended Hückel (EH) model Hoffmann (1988). The TB model has an orthogonal basis with an interaction radius up to the second nearest neighbor (2NN). However, it requires a large number of parameters to include strain effects (48 parameters for copper). In comparison, the EH model has a non-orthogonal basis with a larger interaction radius up the third nearest neighbor (3NN). It requires a smaller number of parameters than the TB method (11 parameters for copper).

Existing parameters for the TB model Hegde et al. (2014) fail when used in highly distorted atomic systems such as GB. Due to the exponential dependence of the inter-atomic coupling on the bond length, the inter-atomic matrix elements corresponding to bond lengths with a 5% or greater distortion generate unphysical results. The problem is solved by obtaining a new parametrization with additional constraints on the inter-atomic coupling. This new parameter set is summarized in TABLE 2 in Appendix A. The parameters for the EH model are taken from literature Cerdá and Soria (2000). Both EH parameters and the new TB parameters show a good match for the Cu unit cell when compared against an ab initio method as shown in Fig. 1. The ab initio result, used as a reference, is obtained by density functional method with a Perdew-Burke-Ernzerhof version of the generalized gradient approximation (GGA PBE) exchange- correlation functional Brandbyge et al. (2002). An energy cutoff of 150 Ry is used and the Brillouin zone is sampled with a 10$\times$10$\times$10 mesh. An FCC copper lattice with a lattice constant of 0.361 nm, as reported experimentally Straumanis et al. (1969), is considered.

III Coincident Site Lattice Grain Boundaries

To validate the tight binding models, the effects of GB scattering were studied for a single coincident site lattice grain boundary. The simplest GB configurations are obtained by a rotation of one of the grains until its lattice vector becomes coincident with the vector of the unrotated lattice  Bristowe and Crocker (1978) as shown in Fig. 2. A fairly small number of atoms ($<400$) is required to construct these systems, which allows the tight binding models to be benchmarked against a first principles calculation as implemented in the ATK package Brandbyge et al. (2002).

CSLs are labeled by $\Sigma N$, where $N$ corresponds to the ratio of the CSL unit cell size to the standard unit cell size. In this work, the CSL GB are generated with GBSTUDIO GBStudio (2014) and relaxed using an ab initio method. The relaxation is carried out with GGA PBE exchange-correlation functional. A Double Zeta polarized basis set is used for copper atoms with an energy cutoff is 150 Ry and the Brillouin zone sampled with a 4$\times$4$\times$1 mesh, until all atomic forces on each ion are less than $10^{-5}$ eV/Å.

Once the ionic relaxation is completed, the transmission spectrum for CSL structures is calculated by the recursive Green’s function method Lake et al. (1997) implemented in NEMO5 Steiger et al. (2011) in an energy range between -2 and 2 eV around the Fermi level with a Brillouin zone sampled with a 30$\times$30$\times$1 mesh.

The integrated transmission spectra in the k space obtained by the tight binding methods are compared against the spectrum obtained by the ab initio method with a similar basis set, energy cutoff and Brillouin mesh as is used in the ionic relaxation. The results in Fig. 3 show that the EH method captures the main features of DFT not only at the Fermi energy ($E_{f}$), but also over a large energy window. On other hand, while the transmission spectrum calculated by TB also shows reasonable agreement with DFT around the Fermi window, it fails to describe the ab initio transmission spectrum for energies away from the $E_{f}$.

Subsequently, the resistance for the CSL GB in the ballistic limit is obtained based on the Landauer formalism assuming a low bias condition Supriyo Datta (1997) as:

[TABLE]

where $G$ is the conductance, $R$ is the resistance, $e$ is the elementary charge, $h$ is Planck’s constant and $T(E_{f},{\bf k})$ is the transmission for a particular wave vector ${\bf k}$ at the Fermi energy. The Fermi levels in Figs. 1 and 3 are calculated at the leads of the device self consistently for DFT and non-self consistently for tight binding models. In this last case, the Fermi level is obtained by integrating over the DOS from $-\infty$ to $E_{f}$ until this value becomes equal to the total number of electrons at a zero temperature approximation Valencia et al. (2016). Following Ref. César et al. (2014) the specific resistances of the CSL grain boundaries are obtained by $\gamma^{R}=(R-R_{B})\,A$, where $R$ is the resistance of the configuration that contains the GB, $R_{B}$ is the resistance of the perfect bulk copper, and A is the grain cross section. The specific resistances for those CSL configurations are calculated by TB and EH and compared to DFT as shown in Table 1

The results in the Table 1 and Fig. 4 show less than 10.4 $\%$ difference in the specific resistance between EH and DFT, and less than $11.2\%$ between TB and DFT. Thus the atomistic methods (TB and EH) are able to describe copper interconnects with reasonable accuracy. These methods are chosen to study GB systems with $10^{3}$ to $10^{4}$ atoms because they require significantly fewer computer resources than the ab initio calculations Valencia et al. (2016).

Only non ab initio methods are capable of relaxing structures of this size ($\gg 10^{3}$ atoms), therefore a force field potential method based on an Embedded Atom Model (EAM) is used. The relaxation is performed using LAMMPS software package Plimpton (1995) with an EAM potential constructed by Y. Mishin et al. that is fitted to first principles calculations to correctly describe grain boundaries and point defects in copper  Mishin et al. (2001). .

 

The accuracy of this approach is determined by comparing the formation energy for CSL GBs obtained by ab initio and the EAM method. The formation energy $\gamma^{E}$ is defined as follows:

[TABLE]

where $E_{slab}$ is the total energy of a slab configuration that contain a CSL GB, $N$ is the number of atoms in the CSL GB, $E_{0}$ is the energy of a single atom of bulk copper and $A$ the cross sectional area. The ionic relaxation carried out by ab initio methods used the plane wave DFT package (VASP) Kresse and Furthmüller (1996) and a PBE GGA exchange-correlation functional. The plane wave energy cutoff is 500 eV and the Brillouin zone is sampled with a 4$\times$4$\times$1 mesh, until all atomic forces on each ion are less than $10^{-5}$ eV/Å.

Comparison of the relaxation energy, computed using the EAM potential, with the DFT result (see Fig. 5), shows that the difference is less than 7% with for all CSL orientations except the $\Sigma_{11}$, which shows a larger error of 20 %. These results indicate that the EAM potential calculation is an acceptable method to relax the grain boundary structures with the benefit of reduced computational burden, compared to DFT.

IV Specific Resistance for Grains of 10 nm length

Based on the prediction of the ITRS roadmap that interconnects will reach 10 to 30 nm length in the coming years ITR (2014), a set of copper thin films of 30 nm is constructed and modeled by tight binding methods as described in Section II. The copper interconnects are formed by three grains of 10 nm length. Each grain is constructed with a super cell growing in the ${[}110{]}$ orientation with a lattice constant of 0.361 nm which has the highest conductance Hegde et al. (2014), as reported experimentally Straumanis et al. (1969). In order to quantify the effect of GB orientation on the specific resistances for copper interconnects, two different types of GBs are generated by Voronoi diagrams Rycroft et al. (2006). These GB types are based on the rotation direction of the middle grain shown as “Tilt” and “Twist” GBs respectively, which generates two boundaries as shown in Fig. 6 a) and b). In order to have a lower impact on the specific resistivity due to the electrode setup, three grains are modeled in this work. In both configurations, only the middle GB is initially rotated then a periodic boundary condition is applied in the ${[}001{]}$ direction for the ionic relaxation and the electronic transport calculation. Therefore, atomic surface roughness is present in the structures as a result of the relaxation. Additionally it is assumed that each configuration shown in the Fig. 6 and 8 is connected to a pristine source and drain lead oriented in the [110] direction, whose atoms are fixed during the ionic relaxation.

 

The “Tilt” GBs are generated by a rotation of the middle grain with respect to the [001] direction by an angle $\theta$ in a range between 0 and $\pi/2$. Each grain is formed by a supercell of 10 nm length ($L$) in the transport direction ${[}110{]}$, 10 nm width ($W$) in the ${[}\bar{1}01{]}$ direction and 0.361 nm thickness ($T$) in the periodic direction ${[}001{]}$ as shown in Fig. 6 a) and c).

The “Twist” GBs are generated by a rotation of the middle grain with respect to the [$\bar{1}$11] direction by an angle $\theta$ in a range between 0 and $\pi/2$. The rotation is applied in the same direction as the periodicity, therefore thicker grains are constructed to ensure the grains overlap after rotation. In this configuration setup each grain is formed by a supercell of 10 nm length ($L$) in the transport direction ${[}110{]}$, 3 nm width ($W$) in the ${[}\bar{1}01{]}$ direction and 3 nm thickness ($T$) in the periodic direction [001] as shown in Fig. 6 b) and d).

It is important to clarify that after any rotation for “Tilt” or “Twist” GB the [110] direction is no longer the transport direction for that grain. Similarly, the rotation angle corresponds to the initial value, but this value will be slightly modified after relaxing the structure.

The specific resistance for “Tilt” and “Twist” GBs for different orientations are obtained by a procedure similar to that described in Section III as $\rho=R\times A$, where $R$ is obtained by Eq. (1) and each configuration is relaxed by an EAM potential. In order to compare the specific resistivity for “Tilt” and “Twist” GBs for different angles $\theta$, the “Tilt” GBs values are normalized such that “Tilt” and “Twist” GBs are calculated over the same cross sectional area. Those values are plotted in Fig. 7. In both systems, specific resistance increases with an increase in the angle, until the angle reaches $\pi/6$, and then becomes almost constant, although the “Tilt” GB shows a reduction after $\pi/3$. The specific resistance dependence for “Twist” GBs shows more noise than for “Tilt” GBs, because “Twist” structure has more points per unit area where the grain boundaries intersect (see Figs. 6c, d), which leads to a higher number of dislocations. Differences between TB and EH, especially pronounced for “Twist” GBs, are due to the fact that the TB model does not correctly describe strained systems, where the atoms are coupled by distances much smaller than the bulk bond length. In particular, for the “Twist” system with rotation angles such as 6, 8 and 70 degrees, there are many atoms with a small distance between nearest neighbors which results in unphysical peaks in the specific resistance dependence (see Fig. 7b) when it is calculated by the TB method.

In order to understand and create a compact model to predict how specific resistance changes as a result of GB orientation, a sample set of 600 configurations are generated. Each GB is constructed with three grains and each of them is rotated with an angle ($\alpha,\beta,\gamma$) in a range between 0 to 180 degrees parallel to the GB boundary. The dimensions of the GB are similar to those used for “Tilt” GB with thickness, width and length equal to 0.5 nm, 3 nm and 10 nm respectively as shown in Fig. 8. A periodic boundary condition in the [001] direction is imposed.

The specific resistance for these samples is calculated with the EH method because it is more reliable over angle rotations than the TB method. Making use of the results obtained from these samples, a boxplot for $\alpha$ and $\gamma$ in a range between [math] to $180$ degrees and a constant angle $\beta$ is plotted in Fig. 9 which shows a symmetry in the specific resistance in a range between 0 to 90 degrees and 90 to 180 degrees. This observation is confirmed by a statistic nonparametric Kolmogorov-Smirnov test Conover (1999) which compares the distribution function for the group of samples in a range between 0 to 90 degrees against those between 90 to 180 degrees and finds that both groups of samples are drawn from an equivalent, continuous distribution. A p-value of 0.16 is obtained for the Kolmogorov-Smirnov test, confirming that there is no difference between the specific resistance distributions for both cases with a confidence of 95$\%$. The symmetry in the specific resistance is due to the fact that the crystal symmetry of copper is not totally disrupted by the structural relaxation. The probability distribution for the three different angles ($\alpha$, $\beta$ and $\gamma$) in a range between 90 to 180 degrees is plotted in Fig. 10.

Per the Shapiro-Wilk test Conover (1999) with a p value of 0.15 and a 95% confidence, the specific resistance distribution follows a normal distribution with a mean and standard deviation equal to 31.7 $\times 10^{-12}\,\Omega\,cm^{2}$ and 2.8 $\times 10^{-12}\,\Omega\,cm^{2}$. The Q-Q plot in Fig. 10 b) shows that the specific resistance distribution is likely normal, although the left and right tails do not follow a normal distribution.

V Grain boundaries Modeled by a Neural Network

Atomistic models based on a tight binding approach can describe the effects of the GB orientation on specific resistance with the same accuracy as DFT methods, but with a much lower computational burden. However, the specific resistance calculated by atomistic models such as EH and TB for a combination of three grains of 10 nm length in the transport direction are still not as fast as conventional models such as the Fuchs-Sondheimer and Mayadas-Shatzkes models  Fuchs and Mott (1938); Mayadas (1969) which describe surface roughness and grain boundary effects respectively in copper interconnects. However these models require experimental input to fit some parameters which limits the transferability for different configurations. Therefore, compact models based on the statistical results obtained from an atomistic model described in Section IV are proposed to describe the scattering effects on grain boundaries for a system of 3 grains of 10 nm length. Three different algorithms are used to construct the compact models, including a polynomial fit, a nearest neighbor search model and a neural network as described in the following subsections. The inputs for the compact models are the orientation angles $\alpha$, $\beta$ and $\gamma$ and the output is the specific resistance of the GB $\rho(\alpha,\beta,\gamma)$. The compact models are trained with a random selection of 80% of the 600 samples plotted in the Fig. 10 and validated with the remaining 20% of the data.

Polynomial Fit

A polynomial fit of second order is carried out based on a least squares adjustment, obtaining the following parametric relationship between the misorientation angles $(\alpha,\beta,\gamma)$ and the specific resistivity:

[TABLE]

The expected values obtained from the model are compared against the remaining 20% of the atomistic data as show in the figure  11. The parametric fitting based on a polynomial approximation displays a poor match with the atomistic results with a 70% variability of the resistivity for the training dataset and a MSE equal to 13.94 $\times 10^{-12}\,\Omega\,cm^{2}$. This result shows that grain boundary effects cannot be modeled as a simple additive effect between each orientation. Therefore, a more complicated dependency exists between the resistivity and the orientation angles.

Nearest Neighbor Fitting

Since the polynomial fit provides a poor fitting for the specific resistance of a GB oriented by the angles $(\alpha,\beta,\gamma)$, a non-parametric model is explored based on a “Nearest-Neighbor” search which uses the “dsearchn” triangulation method implemented in Matlab’s optimization package MATLAB . The comparison between the expected specific resistance and the predicted specific resistance obtained from the testing data is plotted in Fig. 12. While this algorithm exhibits a mean square error for the specific resistance equal to 2.67 $\times 10^{-12}\,\Omega\,cm^{2}$ which is much lower than the error of the polynomial method, it does not support systems that have more than 3 degrees of freedom.

Neural Network Model

Finally, a compact model based on a Neural Network (NN) Hagan et al. (1996) algorithm is introduced. NN models have been widely used to model complex problems; in the TB approach, NN algorithms have been used to describe potential minimization Marim et al. (2003) and materials parametrization Bholoa et al. (2007). In this work, a multilayer neuronal network (MLN) is applied with a back-propagation algorithm Hagan et al. (1996) to quickly obtain the specific resistance of the GB. .

 

The neural network shown in Fig. 13 is achieved after testing different types of neural networks and varying the number of hidden layers. The final system is formed by an input layer, three hidden layers, and one output layer. The input layer p$=(\alpha,\beta,\gamma$) is represented by a row vector of dimension $3\times 1$. The hidden layer is composed of three inner layers $i$ with 10, 6, and 3 neurons, respectively; the weight Wi and bias bi vectors for a given layer $i$ are shown in Fig. 13.

The MLN is implemented in the statistical software R making use of the package Neuralnet Fritsch et al. (2008). The value of the parameters Wi and bias bi are obtained by the gradient descent method Haykin (1999) which minimizes the mean square error of the output layer. In the NN, the functions fi represent logistic functions employed at each layer, except for the last layer f4 to which is applied a linear function.

The mean square error (MSE) prediction for specific resistance for this NN is 1.44$\times 10^{-12}\,\Omega\,cm^{2}$ . The results obtained for the testing data of the MLN are plotted in Fig. 14; the model shows good agreement for low values of specific resistance and larger variability for GB with a specific resistance over the range 29.0 - 39.0$\times 10^{-12}\,\Omega\,cm^{2}$.

Comparing the three methods described above, it is observed that the neural network method has a much smaller MSE than the other methods. It can also be generalized to describe more complicated configurations with different geometries and a number of grains not possible with non-parametric methods such as “Nearest Neighbor” or linear fitting.

VI SUMMARY

In summary, the effect of orientation on grain boundary resistance for copper interconnects is studied using two different atomistic tight binding methods (EH and TB). The transmission spectrum and specific resistance calculated by these methods are benchmarked for coincident site lattice single GB ($\Sigma N$) against first principles calculations. These results show that the EH method captures the main features of DFT in the Fermi window between -2 to 2 eV. On other hand, the transmission spectrum calculated by TB also shows reasonable agreement with DFT around the Fermi window, but fails to describe the ab initio transmission spectrum for energies away from the Fermi energy. Since the computational requirements for tight binding methods are also much smaller than for first principle calculations, the EH method is an effective way to describe the specific resistance of interconnects with lengths greater than 30 nm.

Orientation effects for “Tilt” and “Twist” GBs for copper interconnects of 30 nm length relaxed by a semi-classical EAM potential are also benchmarked against first principles. Rotations perpendicular to the transport direction have a larger effect on the specific resistance of the GB than rotations parallel to the transport direction. Statistical analysis of GB specific resistance shows that the inversion symmetry of copper is still manifested for the considered grain geometry.

Finally, statistical models based on three different algorithms are studied. The parametric model based on a polynomial fit of the misorientation angles $(\alpha,\beta,\gamma)$ shows a poor match with the test results from the atomistic model, confirming that a complex relationship exists between the specific resistance and the orientation angles. While the nearest neighbor model displays a better fit with an error of 2.67 $\times 10^{-12}\,\Omega\,cm^{2}$ , it can only support three degrees of freedom. Among the studied models, the compact model based on neural network is the best algorithm to describe the specific resistance with a MSE lower than 1.44 $\times 10^{-12}\,\Omega\,cm^{2}$. Additionally, the neural network can be used for systems with more than three degrees of freedom.

In this manuscript, the ballistic resistivity due to the grain boundary effect has been studied. While electron phonon scattering are reported to play an important role in copper resistivity at room temperature and when the grains are larger  Graham et al. (2010); Plombon et al. (2006), these effects have not been included in this work. Future work will use the neural network to generate a compact model that includes electron-phonon scattering in addition to grain boundary effects to describe the resistivity for copper interconnects.

ACKNOWLEDGMENTS

This work was supported by the FAME Center, one of six centres of STARnet, a Semiconductor Research Corporation program sponsored by MARCO and DARPA. Support by the US Department of Energy National Nuclear Security Administration under Grant No. DE-FC52-08NA28617 is acknowledged. The authors also acknowledge the staff and computing resources of both the Rosen Center for Advanced Computing (RCAC) at Purdue University and the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (award number ACI 1238993). Finally, the authors would like to thank Dr. Bozidar Novakovic and David Guzman for stimulating discussions about the topic.

VI.1 Appendix

Parameters for bulk copper with the environmental tight binding method (TB) are obtained by direct fitting bulk band structure Hegde et al. (2014), but additional constraints on the inter-atomic coupling are included during the parametrization process.