Multiscale Mathematics Research Project
Description
Objectives of Research Project:
The goal of this research project is to develop mathematical and computational algorithms for adaptive modeling of important classes of physical systems in which multiple spatial and temporal scales are prominent features of the response. The basic approach is built around general methods of a posteriori estimation of modeling (and approximation) error and the control of error through systematic adaptive processes called Goals algorithms. The approach, in theory, can make possible the dynamic generation of multi-scale models capable of delivering results at preset levels of accuracy. Important applications of the mathematical and computational procedures are considered, including the analysis of complex multi-scale behavior encountered in nano-manufacturing of semi-conductor devices.
Figure: Close-up view of chip (left) and modeling challenges (right).
Step and Flash Imprint Lithography Process:
The application of interest in this project is the simulation of the polymerized network that is used in the Step and Flash Imprint Lithography process (SFIL). SFIL is a low cost, high resolution patterning technique that shows great promise as a new technology for nanofabrication and computer chip fabrication.
All Imprint lithography techniques are based on the same principle: a template with a prefabricated topography (usually generated using the electron beam technology) is pressed into a displaceable material which takes the desired shape. The material is then hardened and after removing the template the desired pattern is obtained. The SFIL process, however, notably distinguishes itself, as it uses no projection optics (commonly used in lithography processes), operates at room temperature, and is able of generating relief patterns with line widths smaller than 100 nm. The details of the SFIL process are illustrated in the figure below. An organic polymer layer, the transfer layer, is coated on a substrate, typically silicon. A low viscosity, photopolymerizable, organosilicon solution, the etch barrier, is then dispensed in the area to be imprinted. A transparent template bearing patterned relief structures is then aligned over the coated silicon substrate. The template is lowered into contact with the substrate, displacing the etch barrier and trapping the photopolymerizable liquid into the template relief. Ultraviolet radiation through the backside of the transparent template subsequently cures the etch barrier into a cross-linked polymer film. A fluorocarbon release layer on the template allows separation from the substrate, leaving an organosilicon replica of the template pattern on the coated substrate. A halogen etch is then used to break through the undisplaced residual layer of the etch barrier to expose the transfer layer. Finally, a reactive ion etch (RIE) is used to transfer the image through the transfer layer, amplifying the aspect ratio of the imprinted image. For the desired geometric precision of the device, control of the profile of the etch barrier is essential, but presents many technological challenges.
Figure: Step and Flash Imprint Lithography process.
Polymerization step:
In the SFIL process, the target etch barrier material is produced in two steps: (1) a solution of chemical constituents flows into the template relief (the template contours) as the template is pressed toward the transfer layer; (2) ultra-violet (UV) light is passed through the template (which is translucent quartz) into the mixture. The application of ultraviolet light causes chemical reactions between the chemical constituents which leads to the formation of chain-like macromolecules consisting of various repeated molecular units. The units are monomers and the macromolecules are polymer chains. This chemical process is referred to as polymerization.
The equations of chemical kinetics only provide the global species concentrations as functions of time (and the initial species distributions). To determine the molecular structure at the termination of the polymerization process, one must follow the likely chemical reactions that can occur between the molecular components. The possible conformations resulting from reactions of a given initial distribution of constituents can be generated through a Monte-Carlo-like algorithm that shall be referred to here as the kinetic Monte-Carlo process. The key to simulating the polymerization process in this way is to observe that the rate coefficients from chemical kinetics for each reaction is related to the classical Arrhenius law. Within this context, the Arrhenius law provides the probability of reaction occurrence. We introduce three-dimensional regular lattice L where the number of lattice sites is set equal to the estimated number of constituent molecules in a cube of the initial etch barrier mixture. There are five constituents in the SFIL process under study: the monomer M1, the monomer M2, the cross-linker XL, the initiator I, and possible voids V. The voids are introduced to allow diffusion of the constituents during the model process.
The first step of the algorithm consists in the population of the lattice. A molecular constituent is assigned to each cell as follows: (1) boundary cells are assigned either a template molecule or a transfer layer molecule, depending on the location of the part of the boundary where the cell is located; (2) each lattice site is visited in order and a uniform random number r is selected between 0 and 1. If r falls into an interval representing the fraction of constituent j, then the cell is assigned this constituent; (3) a random swapping procedure of the cells is used to further ``randomize'' the lattice. The figure (left) below illustrates the lattice placement process schematically.
The second step consists in applying the kinetic Monte-Carlo type algorithm for the polymerization of the new populated lattice. A schematic of this algorithm is presented in the figure (right) below. Initiation: if an initiator is randomly selected that is not a free radical, then it is made a free radical if the reaction is determined to occur. This is depicted in (a) from left to right. Propagation: if a free radical is randomly selected, then a random neighbor is selected. If a bond has not been formed, t hen a bond is formed if the reaction is determined to occur. This is depicted in (b) from left to right. Void diffusion: if an unreacted particle is randomly selected, then a random neighbor is selected. If that neighbor is a void, then the cell location of the void and the neighbor is switched. This is depicted in (c) from left to right.
Figure: Population and propagation steps in polymerization process.
At the conclusion of the kinetic Monte-Carlo process, the location of the site of each constituent and the connectivity of bonds forming the polymer chains is known. It is observed experimentally that upon completion of the polymerization process, a volume shrinkage of approximately 20 percent occurs upon removal of the quartz template. To account for this densification effect, bond potentials must be assigned to the polymerized etch barrier and a mathematical model must be formulated to describe the motion of the molecules due to the formation of the bonds. To this end, it is assumed that the bulk deformation of the polymer is a quasi-static process that can effectively be modeled using molecular statics.
Mechanics of the Polymer - Molecular Statics:
The goal of molecular statics is to minimize the global energy function of of the lattice structure. It is assumed here that the potentials are pair potentials and that there is no external loading. The covalent bonds between molecules are modeled here with harmonic spring potentials with spring stiffness k and unstretched length l (see table below). Weaker harmonic springs have also been used to describe van der Waal's interactions. The van der Waal's potentials are only assigned to the nearest neighbors and only 18 of the 26 neighbors are included (see figure below). Solutions of the minimization problem are obtained using the software packages PETSc and TAO. Figure below show the configuration of the polymer lattice before and after densification.
Figure: The bonding configuration for each particle. The green particles represent neighboring particles which are allowed to bond to the red; they can be convalent or van der Waal's bonds. Blue particles are not allowed to bond to the red particle.
| Monomer 1 | Monomer 2 | Cross-linker | Initiator | Substrate | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| k | l | k | l | k | l | k | l | k | l | |
| Monomer 1 | 1.10 | 5.40 | 0.80 | 5.50 | 1.20 | 5.70 | 0.90 | 4.80 | 1.00 | 5.34 |
| Monomer 2 | 0.80 | 5.50 | 1.00 | 5.20 | 0.70 | 4.90 | 1.30 | 5.90 | 1.00 | 5.34 |
| Cross-linker | 1.20 | 5.70 | 1.30 | 5.90 | 1.50 | 5.30 | 0.60 | 5.80 | 1.00 | 5.34 |
| Initiator | 0.90 | 4.80 | 0.70 | 4.90 | 0.50 | 5.80 | 0.60 | 5.10 | 1.00 | 5.34 |
| Substrate | 1.00 | 5.34 | 1.00 | 5.34 | 1.00 | 5.34 | 1.00 | 5.34 | - | - |
Table: Spring stiffness and unstretched length used for bonds between polymer molecules.
Figure: Molecular modeling of polymer feature before and after densification.
Continuum model:
A critical step in the multiscale modeling procedure developed here is the construction of a continuum model that represent events at larger scales than the base model, but are compatible to the base model in some sense. Generally, coarser-scale models may result from averaging the features of fine-scale models through various homogenization methods or ensemble-averaging techniques. In this case, a continuum model is chosen, but one whose corresponding constitutive equation coefficients are unknown. A scheme is described below that is designed to determine these unknown coefficients using ``virtual'' (numerical) experiments that are performed on a representative volume element (RVE) of polymer material.
The polymerization process described earlier does not involve any features that would lead to macroscale inhomogeneities or anisotropies, and the monomer and crosslinker constituents were designed to avoid or minimize rate effects. Thus, the assumption that the process leads to a material behaving as an amorphous elastomer seems reasonable. The macroscale model of the polymerized etch barrier is then that of a homogeneous, isotropic, hyperelastic material with a stored energy per unit volume relative to a reference configuration. We consider here the right Cauchy-Green deformation tensor C=FFT where F is the deformation gradient tensor of the motion. A fundamental condition on the stored energy function W is that it be form-invariant under all changes in the observer frame of reference, which, generally, means that W must depend on invariants of the deformation. One way of guaranteeing this invariance is to write W as a function of the principal invariants of C:
- A cube D of polymerized material, initially of size L0 with N particles, is generated using the kinetic Monte-Carlo algorithm, and is allowed to assume a relaxed (densified) equilibrium configuration. The cube D is the initial RVE.
- Uniform tractions are applied over opposite faces of the cube by assigning values of the net force on molecules residing on the near-planar boundaries of D.
- The resulting principal stretches λ1, λ2, λ3, are calculated by computing the length of the stretched RVE and by taking the ratio with the initial length of the relaxed RVE: λ1 = Lx/L0, etc. These lengths are calculated by averaging the positions of the particles on each face and computing the difference of these averages. Furthermore, the Jacobian J=det F = (λ1,λ2, λ3)1/2 is calculated from the RVE.
- The total energy ED in D is computed in the relaxed and stretched configurations.
-
The total energy density in D is equated to the continuum
stored energy density in D:
W = ED / V0 where V0 is the initial volume of D.
Figure: RVE testing procedure. Step 1: Relaxation of the polymer lattice. Step 2: Deformation (uniaxial and biaxial stretches). Step 3: Measurements of volume, energy, and principal stretches of material lines. Step 4: Calibration of model parameters by curve fitting.
Important aspects of this process are noted:
- The constants determined by this process should be independent of the size of the RVE; thus the dimension of the domain D should be increased until the computed values of the material parameters do not change.
- A single molecular model of the RVE represents only one realization of the polymeric structure. Thus, a large number of such realizations should be generated to determine the statistical variations of the material parameters with sufficient accuracy.
- The form of W assumed in the process should represent a stable characterization of the material, relatively insensitive to small perturbations in the molecular model.
Coupling approach:
This section is concerned with the development of a general scheme to systematically couple particle and continuum models in order to produce surrogates of the molecular base model. The strategy involves enforcing the displacement and derivative constraints between the two models. The constraints are achieved using Lagrange multipliers on a region of overlap between the continuum and particle models. This scheme is an adaptation of the Arlequin method used for coupling two continuum models with differing scales of finite element discretization (a global coarse mesh and a local fine mesh). Well-posedness of the method for a one-dimensional problem coupling a harmonic spring model with a linearly elastic rod has been shown. Here, we extend the formulation to the coupling of the three-dimensional, nonlinear polymer model and the nonlinear, elastic continuum model discussed previously (see figures below).
Figure: (Left) Example of a geometry for the arlequin method. The solid green represents the continuum model while the red, blue, and yellow particles represent monomer, cross-linker, and initiator, respectively. Notice that the zoomed shows the particle region, while the particles contained in the green define the overlap region. (right) Equilibrium configuration of an Arlequin approximation to a 51x51x51 polymer lattice. The mesh is colored by elements of the zz-component of the Cauchy stress while the slice is the interpolated zz-Cauchy stress.
Adaptive modeling:
The adaptive algorithm is driven by error estimates in quantities of interest so as to control the position of the interface between the continuum and particle models. The error e is defined here are the difference between the solution u0 of the surrogate model provided by the coupling method and the solution u of the full particle model. The optimal position of the interface is the one that reduces the error E = Q(u) - Q(u0) in the quantity of interest within some given tolerance. The error E is estimated using the framework of so-called goal-oriented error estimation. It can be shown that this error can be represented in terms of the solution p of the dual problem associated with the particle model:
The quantity R(u0; p) or an approximation of this quantity can be decomposed into local contributions on predefined cells (elements) and the adaptive algorithm consists in switching continuum cells to particle cells for which the contributions are the largest. We show in the figures below the successive adaptive steps in the case of a cubic 21x21x21 polymer lattice for which the relative error in the displacement at the center of the top face has been reduced to a tolerance of 6%.
Figure: (Left) Configurations chosen by the adaptive algorithm for a 21x21x21 polymer lattice with an Arlequin surrogate model. (Right) Exploded view of the configurations chosen by the adaptive algorithm. Red cells denote particle model regions, green cells the continuum model, and yellow cells denote the interface regions.
This description of the project has been adapted from the manuscript: P.T. Bauman, J.T. Oden, and S. Prudhomme, Adaptive multiscale modeling of polymeric materials: Arlequin coupling and goals algorithms, Computer Methods in Applied Mechanics and Engineering, Submitted (2008).