3D Simulation of Nano-Imprint Lithography

A proof of concept study of the feasibility of fully three-dimensional (3D) time-dependent simulation of nano-imprint lithography of polymer melt, where the polymer is treated as a structured liquid, has been presented. Considering the flow physics of the polymer as a structured liquid, we have followed the line initiated by de Gennes, using a Molecular Stress Function model of the Doi and Edwards type. We have used a 3D Lagrangian Galerkin finite element methods implemented on a parallel computer architecture. In a Lagrangian techniques, the node point follows the particle movement, allowing for the movement of free surfaces or interfaces. We have extended the method to handle the dynamic movement of the contact line between the polymer melt and stamp during mold filling.


Introduction
In the recent years, a considerable effort has been made in the development of polymer-based micro-and nano-fabrication techniques for applications in micro-electromechanical systems (MEMS). Techniques such as electron beam lithography (EBL) and atomic force microscopy lithography (AFML) have proven successful in transferring nano-patterns with line width of the order of few tenths of nanometers. However, the cost and the deficient processability for mass production have caused the emergence of promising techniques that can retain or improve the resolution of EBL and AFML and circumvent their inherent shortcomings. One of these techniques is thermal nanoimprint lithography (henceforth simply referred to as nanoimprint lithography or NIL) [1] offering resolutions below 10 nm.
Nano-imprint lithography (NIL) is essentially a two-step process. In the first step, a resist film (a thermoplastic polymer) is cast onto a hard substrate and heated above its glass transition temperature. A rigid stamp (or mold) containing some well-defined structures patterned on its surface is brought into contact with the film and pressed in so that the negative replica of the stamp will be transferred. The pressure load is held for a period of time and then the stamp is cooled down until its temperature is below the glass transition temperature of the polymer. This ensures the mechanical stability of the film and then demolding can be carried out.
The NIL operational window [2] covers pressures up to 200 bars, temperatures of about 200°C (it has been experimentally found that operation temperatures of 70-90°C above the glass transition temperature of the polymer are suitable for the imprint), and stamp sizes of 2 cm 9 2 cm. The fact that the patterning can be done over large areas is particularly attractive in connection with the manufacturing of high-volume data storage and highspeed data-processing components.
Nano-imprint lithography (NIL) can be regarded as a particular case of hot embossing lithography (HEL), but with the specific characteristic that the magnitude of the thickness of the samples is comparable to the height of the indents in the mould. Hence, the patterning is not a surface modulation process, as the HEL, but demands the transport of a large amount of the polymer from areas in contact with the stamp into the cavities. In most cases, NIL is an isothermal process where both stamp and sample have identical temperatures.
Many related NIL publications have appeared in the last decade. Only few of them have consistently addressed the investigation of polymer flow during imprint. Different experimental approaches have shown that the flow mechanism of polymers into nano-structures [3,4] is not trivial, and this has been attributed to a variety of causes such as the type of flow, the stamp geometry, the viscoelastic properties of the polymer, as well as the role played by surface tension effects. The latter includes attraction forces of different nature between the stamp and the polymer [5,6].
The experimental work has often been followed by attempts to simulate the flow during imprint mainly for the cases of periodic squared arrays or squared cavities. Simulations have been carried out using finite element methods, and all the flow problems were solved in a two-dimensional framework. Many of the numerical investigations adopted a general scope focusing in the understanding of the flow into micro/submicrometer structures.
Most of the early approaches treated the polymer either as a purely viscous fluid [7,8] (sometimes using generalized Newtonian models) or as purely elastic materials [9]. In a few cases, the surface tension of the polymer was taken into consideration but interfacial tension between polymer and stamp was not accounted for. To our knowledge, only two approaches considered viscoelasticity in the material. Hirai et al. [10] modeled the polymer flow using a linear viscoelastic constitutive equation and established a qualitative comparison with their experiments. The only study with relevance to the NIL or HEL process, applying a nonlinear constitutive equation is that of Eriksson et al. [11]. Eriksson et al. [11] addressed both experimentally and numerically a compression molding on the millimeter scale combined with the filling of periodic squared arrays at micrometer length scales. Probably due to the applied nonlinear constitutive equation, they obtained a quantitative agreement between the filling profiles measured in the experiments and those provided by numerical simulations. The latter approach used an integral constitutive equation of the K-BKZ type [12,13], within a Molecular Stress Function (MSF) approach [14], to describe the dynamics of the polymer.

Polymer Melt as a Structured Liquid
The physics in the NIL is the isothermal flow of polymer melts. Polymers are structured systems where order phenomena develops the stresses in the material during flow (Pierre-Gilles de Gennes [15]). The reptation theory by Doi and Edwards [16] was the first contribution toward an exact formulation for the flow of polymers. Since then progress has been scarce. The 'interchain pressure' concept by Marrucci and Ianniruberto [17] is currently the only theoretical approach capable of accurately predicting the published homogeneous flow data for molar mass distributed as well as structurally well-defined polymer melts [18,19,20,21,22,23,24]. Just recently, Wagner et al. [22,23] and Rasmussen et al. [24] have suggested models capable of modeling the flow of polymer melts in general. These models are based on the 'interchain pressure' concept by Marrucci and Ianniruberto [17], and a Molecular Stress Function constitutive model [14]. The latter is a generalization of the ideas and model of Doi and Edwards [16].
The general form of MSF constitutive model is written as a memory-weighted time integral over a strain tensor and the square of the molecular stress function, f; Here, M(tt 0 ) is the linear viscoelastic memory function and the strain, . . . h i, is the independent alignment tensor from the Doi-Edwards reptation theory. The theoretical basis is the idea of a tube segment of unit length and orientation given by the unit vector u ¼ ðu 1 ; u 2 ; u 3 Þ. In the stress-free state, u is deformed into E Á u in the current state. The components of the macroscopic displacement gradient tensor, E, are given by E ij ðx; t; t 0 Þ ¼ ox i =ox j 0 ; i = 1, 2, 3 and j = 1, 2, 3. x 0 ¼ ðx 0 1 ; x 0 2 ; x 0 3 Þ are the coordinates of a given particle in the stress-free reference state (time t 0 ), displaced to coordinates x ¼ ðx 1 ; x 2 ; x 3 Þ in the current state (time t). The angular brackets, . . . h i, denote an average over a unit sphere . . . h i ¼ 1=ð4pÞ R juj¼1 du. The variables in the displacement gradient tensor, x and x 0 , are indexed by a Cartesian coordinate system attached to the particle. This corresponds to the use of Lagrangian variables.
Please notice, the flow physics depends on the choice of the molecular stress function, f [22,24].

Flow Modeling of Polymer Melts in NIL
The modeling of the NIL process requires time dependency in the applied method. A scant effort has been made in the development of numerical methods capable of handling time-dependent flow of integral constitutive equations [25,26,27]. To perform a realistic modeling, it will (at least) require a fully three dimensional as well as free surface viscoelastic flow. The method of Wapperom [26] and the one by Rasmussen [27] both apply Lagrangian particle variables. Therefore, both methods can handle (large) movement of moving surfaces or interfaces without extra effort. The movement can be free or specified. Only the approach by Rasmussen [27] has been numerically formulated in three-dimensional (3D), although the step from 2D to 3D in most cases is a minor problem.
The major concern in fully three-dimensional computations is the immense increase in the number of unknown variables. Three-dimensional time-dependent computations in most cases require an efficient numerical formulation as well as code parallelization. Therefore, we will use the recent method by Marín et al. [28]. This method is a Lagrangian finite element method, based on a proven convergent Galerkins principle [29]. In a Lagrangian method, the node point follows the particle movement. The particle variables and the pressure field are approximated by quadratic (tetrahedral) and linear interpolation functions, respectively. The method is third-order accurate in both space and time, and is numerically stable for all chosen time steps sizes. The equations are solved by a robust Newton-Raphson iterative scheme and the code is fully parallelized.
The use of (Lagrangian) variables indexed by a coordinate system attached to the particle enables the numerical method to handle arbitrary large movement of the material freely moving surfaces or interfaces. The modeling of a NIL, as well as any other mold filling problem, is complicated by the presence of dynamic movement of the contact line between the melt and stamp during filling. It is therefore essential to apply an easy handling of this. Here, we treat the contact of the particles as follows: the time and position of the contact to a solid surface of the particles in nodal points (in the finite element discretization) was calculated with an explicit second order prediction. Any particle (e.g., node point) is then attached (e.g., sticking) to the solid mold, starting from the time (and position) of contact.
The finite element flow solver used here has been adapted to simulate the filling of an arbitrary stamp morphology. We will consider the example depicted in Fig. 1 where a stamp containing a surface patterned with periodic squared sinus cavities aligned with each other and equally spaced is pressed against a polymeric film of thickness T. Exploiting the symmetry, the simulation domain can be reduced to the region delimited by the blue line.
For simplicity, we will only use a molecular stress function f = 1. The modeling is then based on a liquid described constitutively by the independent alignment strain tensor from the Doi-Edwards reptation theory only. This corresponds to the exact flow properties of polymers melts when the 'Interchain Pressure' is relaxed. Currently, a relation between the 'Interchain Pressure' relaxation and the linear viscoelastic relaxation has not been established, although the 'Interchain Pressure' seems to relax at time scale in the size of the Rouse time. Further, we have applied a memory function as a continuous BSW spectrum [30] given by Mðt À t 0 Þ ¼ ðg 0 ð1 À nÞ=k 2 Þððt À t 0 Þ=kÞ Àðnþ1Þ Cðn þ 1; ðt À t 0 Þ=kÞ: CðÁ; ÁÞ is the incomplete gamma function, g 0 the zero-shear viscosity and k the maximal relaxation time constant. The BSW represent the linear viscoelastic dynamic of a monodisperse melt and is easily extended to broadly distributed polymers. In all computations, we use a Currie [31] approximation of the independent alignment tensor from the Doi-Edwards reptation theory.
A non-dimensionalization of the the stress (e.g., equation of motion), using non-dimensional variables as x Ã i ¼ x i =T; t * = t/k and r Ã ij ¼ r ij =G 0 N , can be applied. G 0 N is the elastic plateau modulus for the polymer and T the initial thickness of the substrate. Introducing these non-dimensional variables, the mass conservation as well as the momentum equations will only contain one non-dimensional parameter, n. We have used a fixed value of n = 0.07 here.
The initial unstructured finite element mesh is shown in Fig. 2 (top and left) along with the shape of the cavity (denoted by the solid line). Non-slip boundary conditions are specified on the bottom plane of the domain. The three lateral planes are treated as symmetry planes. The top plane contains two differentiated regions; the semicircular one that is treated as a free surface and the secondary region formed by the nodes initially in contact with the mold. The flow can be induced either moving the initial contact area at a given velocity or applying a pressure load on it. Here, we have applied a constant velocity. Figure 2 shows the progressive flow of the polymer in the cavity. The stamp is pressed in the polymer at constant velocity. As the fluid flows into the cavity, new nodal points get in contact with the stamp surface and they then follow the motion of the stamp. The face coloring is a measure of internal sample pressure. Pressures are normalized with the maximum pressure recorded during the Fig. 1 Illustration of the NIL process simulated with the present numerical method. The green domain represents the polymer cast onto a hard substrate and the grey domain represents the cavities patterned in the stamp. The dimensionless aspect ratio, relative to the initial thickness of the sample T, is H/T = 4. The separation between cavities is L/T = 20. The blue line denote the simulated domain simulation. Red and blue colors indicate high-and lowpressure areas, respectively. One has to notice; throughout any Lagrangian simulation, the mesh undergoes distortion. Therefore, several new meshes need to be regenerated in order to continue and complete the computations. Despite the difficulties in mesh generation, it gives the possibility to relocate elements where they are needed.
The actual movement of the free surface is similar to what commonly is observed experimentally as discussed in detail in Rowland et al. [8]. The compression of the melt between the hard substrate and the stamp forces an outward flow which results in the creation of a central suppression in the surface. Only in the final stage of the filling, a single peek flow appears due to the narrow gap [8].
Please notice, if there are significant differences in temperature between the stamp and the hard substrate a non-isothermal approach should be applied. The approach to non-isothermal flow commonly follows the assumption of polymer melts as a thermorheological simple material (Morland and Lee [32]), applying a pseudo time (Crochet and Naghdi [33]). The currently only existing numerical method for non-isothermal time-dependent flow of polymers as a structured liquid based on the pseudo time approach is the 3D Lagrangian finite element method by Marín et al. [34]. This numerical method is a non-isothermal extension of the numerical method used here.

Conclusion
To summarize, the theoretical exploitation of the NIL process for technical purposed seems to be feasible. It has been demonstrated that the time-dependent modeling of a NIL process, similar to hot embossing, is possible in fully three dimensions. A Lagrangian or particle finite element approach has been applied in the numerical simulations. As the Lagrangian formulation is the natural basis for structural-based constitutive models (e.g., integral models), it is capable of using a correct physical basis for the flow physics of polymer melts as a structured liquid.