Simulation of mechanical properties of epoxy-based chemically amplified resist by coarse-grained molecular dynamics

Polymer 53 (2012) 4834e4842

Contents lists available at SciVerse ScienceDirect

Polymer journal homepage: www.elsevier.com/locate/polymer

Simulation of mechanical properties of epoxy-based chemically amplified resist by coarse-grained molecular dynamics Hiromasa Yagyu a, *, Yoshikazu Hirai b, Akio Uesugi b, Yoshihide Makino b, Koji Sugano b, Toshiyuki Tsuchiya b, Osamu Tabata b a b

R&D Division, Mitsuboshi Belting Ltd, 4-1-21, Hamazoe-dori, Nagata-ku, Kobe 653-0024, Japan Department of Micro Engineering, Kyoto University, Japan

a r t i c l e i n f o

a b s t r a c t

Article history: Received 30 March 2012 Received in revised form 18 August 2012 Accepted 21 August 2012 Available online 28 August 2012

The mechanical properties of an epoxy-based chemically amplified resists with various cross-linking ratios were simulated using a newly developed coarse-grained molecular dynamics simulation that employs a beadespring model. Models with the different cross-linking ratios were created in the molecular dynamics calculation step and uniaxial elongation simulations were performed. The results reveal that the simulated elastic modulus of the resist modeled by the KremereGrest model with an extended angle bending potential depends on the cross-linking ratio, its dependency exhibits good agreement with that determined by nanoindentation tests. Ó 2012 Elsevier Ltd. All rights reserved.

Keywords: Beadespring model Elongation Elastic modulus

1. Introduction Epoxy-based chemically amplified resists (hereafter, resists) were originally developed as masks for substrates such as silicon and glass during etching, but they have been attracting attention as three-dimensional microstructural polymers for microstructuring microelectromechanical systems (MEMS). Fabricated permanent polymer structures have been utilized in various MEMS, including microfluidic devices, because of their low cost and superior mechanical, optical, and chemical properties. An interesting approach along this direction is the utilization of a semi-crosslinked resist proposed by Hirai et al. [1,2] and its application to a microfluidic device. They investigated the relationship between the developer permeability and the cross-linking ratio by varying photolithography parameters such as the UV dose and the postexposure-bake (PEB) temperature and time, and found that a cross-linking ratio between 0.5 and 0.75 is suitable for fabricating centimeter-long embedded microchannels in a negative resist. Based on this approach, a technique for fabricating embedded microchannels in a resist was proposed. They also demonstrated that a semi-cross-linked resist functions as a nanofiltration

* Corresponding author. Tel.: þ81 78 682 3381; fax: þ81 78 687 7516. E-mail addresses: [email protected], [email protected] (H. Yagyu). 0032-3861/$ e see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.polymer.2012.08.050

membrane with nanopores [3]. Thus, utilization of a semi-crosslinked resist promises a rapid and effective technique for fabricating complex microchannel networks in microfluidic devices. This technique does not suffer from limited developer diffusion from the open end [2]. The functionalization of microfluidic devices with a nanofiltration membrane can be realized without using an additional materials and complex microfabrication processes, making it highly attractive to biophysicists and bioengineers because filtration is an important step in applications such as tissue engineering and regenerative medicine [4]. To design microfluidic devices that have filtration membranes with the semi-cross-linked resist, it is essential to determine the dependences of its mechanical properties and the diffusion properties of various solvents through the membrane on the crosslinking ratio of the resist. As a first step toward this goal, we performed nanoindentation tests to measure the elastic moduli of resists with various cross-linking ratios, the results revealed that the elastic modulus of the resist depends on the cross-linking ratio [5]. As the next step, the relationship between the pore size, which dominates diffusion, and the cross-linking ratio needs to be clarified. However, it is difficult to experimentally determine the pore size to an accuracy of a few nanometers to a few 10s of nanometers. We found that the molecular dynamics method is useful for clarifying the cross-linked structures of polymers [6,7]. The

H. Yagyu et al. / Polymer 53 (2012) 4834e4842

molecular dynamics method is also expected to be able to calculate the elastic modulus and pore size of a semi-cross-linked resist by modeling cross-linked polymer structures in the glassy region. However, a fully atomistic molecular dynamics simulation cannot be performed on a general purpose computer since the model with large number of cross-links was required for analyzing the dependency of the elastic modulus and the pore size on the number of cross-links in the complex threedimensional molecular structure of a resist. The difficulty in fully atomistic model is related to procedure of cross-linking process (e.g. avoiding generated excessive force at bond by compulsory reaction in the modeling). Varhney et al. [8] presented an atomistic model based on molecular modeling for cross-linked epoxy polymer. They proposed a multistep relaxation procedure for relaxing the molecular topology during crosslinking, and calculated several material properties of the model.


 U LJ rij ¼

8 > > <

"( 43 ij

> > : 0;

sij rij

!12 

sij rij

!6 )

 

12  6 # sij ;  c c

sij

2. Simulation method 2.1. Beadespring model The coarse-grained molecular dynamics simulation employs the beadespring polymer-chain model proposed by Kremer and Grest [12]. This model utilizes a polymer chain that consists of i beads connected by the following potential energy U



 U ¼ U FENE rij þ U LJ rij

(1)

where rij is the distance between beads i and j. The bond potential energy UFENE(rij) and non-bond potential energy ULJ(rij) are given by

  2 
 rij 1 U FENE rij ¼ kR20 ln 1  2 R0


 rij  c 
rij >c

Recently, Nouri and Ziaei-Rad [9] reported cross-linked epoxy polymer model using fully atomistic molecular dynamics simulation with Graphical Processing Units (GPUs). Since the model has sufficient scale in order to evaluate the effect of cross-linking ratio on the property of the model, the relation of mechanical property and cross-linking ratio was investigated. However, although these large scale models can also be calculated using the general parallel molecular dynamics code and a multi core processor as well as GPUs, there are some difficulty in use of the system easily, and total calculation including become expensive. To overcome this problem, coarse-grained molecular dynamics simulations using the KremereGrest model (beadespring model) have been performed [10,11]. In the molecular models of coarsegrained molecular dynamics simulations, a single bead represents several monomer units. Although simulations cannot account for the chemical effects of chain linkages in the model, the chain dynamics during deformation on a microscale and over a long time scale can be simulated by employing an appropriate potential energy between the beads, such as bond and non-bond potentials [12]. The beadespring model has been used to model non-equilibrium situations such as shear [13] and elongation flow [14] by employing the above viewpoint [15]. Although some studies [10e15] have used the beadespring model to investigate polymer melts, the beadespring model has not been used to evaluate a cross-linked resist. Furthermore, the dependencies of the elastic modulus and pore size of the polymer on the crosslinking ratio have not been researched. We propose the quantitative simulation focusing on mechanical property of the glassy cross-linked polymer model using simple beads-spring models based on these prior studies. In this study, we modeled an epoxy-based chemically amplified resist using a simple beadespring model that includes the angle bending potential. To investigate the mechanical properties of the resist, uniaxial elongation was studied by performing a coarsegrained molecular dynamics simulation. The elastic modulus of the resist models was then calculated from the elongation simulation results and the calculated dependence of the elastic modulus on the cross-linking ratio was compared with that measured by nanoindentation tests.

4835

(2)

(3)

where k is the finitely extensible nonlinear elastic (FENE) spring constant, R0 is the maximum extension of the spring, sij and 3 ij are the LennardeJones parameters between beads i and j, and c is the cutoff distance. The bond potential UFENE(rij) is the FENE potential, it is harmonic at its minimum but bonds cannot be stretched beyond a maximum length determined by R0. The non-bond potential ULJ(rij) describes the non-bond interaction between beads separated by a distance rij and it is given by the LenardeJones potential. In this study, we extend the basic KremereGrest model described above by including the angle bending potential UANGLE to simulate the deformation dynamics of a cross-linked glassy polymer. The following angle bending potential is a function of the angle q and it is defined for three connected beads (i, j, and k).




 U ¼ U FENE rij þ U ANGLE rijk þ U LJ rij

(4)

2
 1  U ANGLE rijk ¼ kq qijk  q0 2

(5)

where

cos qijk ¼



 rj  ri $ rj  rk rij rjk

(6)

kq is the angle bending spring constant, q0 is the equilibrium angle, and ri is the position of bead i. The bending angle qijk is defined as the inner angle and is 180 for a linear molecular chain. The effect of the angle bending potential is discussed in Section 4.5. The time evolution of the bead at position ri is calculated using the following Langevin equation:

m

d2 ri vU dr ¼   G i þ Wi ðtÞ vri dt dt 2

(7)

where m is the mass of the bead, U is the total potential energy of the system, G is the friction constant, and Wi(t) is the Gaussian

4836

H. Yagyu et al. / Polymer 53 (2012) 4834e4842

white noise, which is generated according to the following equation:

Wi ðtÞWj ðt 0 Þ ¼ 2kB TmGdij I dðt  t 0 Þ

(8)

where d is the delta function, T is the temperature of the system, kB is the Boltzmann constant, and I is the identity matrix. The calculation is performed using reduced units for which s, 3 , and m are used as units of length, energy, and mass, respectively. The unit of time is s (¼s(m/3 )1/2) and the unit of temperature is 3 /kB. In this study, the friction constant G was set to 0.5s1 [12]. The equation of motion is solved using the velocity Verlet algorithm with a time step Δt ¼ 0.01s. The simulation was performed using periodic boundary conditions for all walls of the cell. The density of the system was set to 0.85ms3 under the condition of constant temperature and density (NVT ensemble). When evaluating the glass transition temperature of the model (see Section 2.4), the simulation was performed using the constant-pressure condition of Andersen (NPT ensemble). All the calculations were performed using the coarse-grained molecular dynamics program OCTA/ COGNAC [16]. 2.2. Resist model Fig. 1(a) shows the chemical structure of the monomer before exposure to a conventional epoxy-based chemically amplified resist such as SU-8Ô (MicroChem Co.) [17,18] or TMMRÔ (Tokyo Ohka Kogyo Co. Ltd.) [19]. The resist is based on bisphenol A diglycidyl ether (BADCE), a common industrial reagent. Four BADCE were connected by methylene group. Each molecule has eight epoxy groups on average. This is averaged structure of the resist used for experiment in Section 3. A single molecule consists of epoxy, phenyl, and isopropylidene groups and the molecules are connected between the phenyl groups with methylene group. The resist consists mainly of an epoxy resin and a photoacid generator. Photoacid (Lewis acid, designatad HþX) is photochemically produced in a solid photoresist by UV absorption of the photoacid generator. The chemical structure of catalyst is not shown in the manuscript because the structure of catalyst does not influenced the calculation and modeling in this simulation. The photoacid acts as a catalyst and it reacts with the epoxy group to open the epoxy ring in the subsequent cross-linking reaction that occurs during the PEB. Exposed resist contains acid catalyst, whereas unexposed resist does not. The epoxy groups with opening rings have two reactive bonds and react with another epoxy group, ether, in the same or different molecule. It is necessary to perform a PEB because

a

O

O

O

O

Table 1 Parameters of resist model. LennardeJones length sij LennardeJones minimum energy 3 ij LennardeJones cutoff distance c FENE spring constant k Maximum extension distance R0 Angle bending spring constant kq Angle bending equilibrium angle q0

1.0 1.0 2.0s 303 /s2 1.5s 0, 10, 25, 503 /rad2 109.5

the reaction kinetics of the cross-linking mechanism is very slow at ambient temperature. Heating the resist above its glass transition temperature increases its molecular motion, thereby assisting the cross-linking reaction. Consequently, most cross-linking reactions occur during the PEB [20]. In this study, the resist was modeled using a simple beade spring model, as shown in Fig. 1(b). The resist monomer was defined as branched linear polymer, and the model used same parameters (e.g., weight, bond length, and interaction) for all beads in order to realize a fast calculation because this study focused on simulating mechanical property of the resists using the simple model. However, to realize cross-linking during the molecular dynamics simulation, beads at the end of the molecule have a reactive function. The resist monomer model consists of 20 beads. The equilibrium bond length r0 between the beads was set to s. The spring constant k and the equilibrium bond length R0 in the FENE bond potential were defined as 303 /s2 and 1.5s, respectively. There is a tradeoff between the maximum spring constant of bond potential and time step Δt. These parameters are the values at which no crossing of the beads occurred, and the dynamics of polymer can be simulated [10e12]. The equilibrium angle q0 in the angle bending potential was set to 109.5 , which is the equilibrium angle of a carbon network consisting of linear polymer chains [21]. To determine a suitable parameter for the epoxy resist, four angle spring constants kqq ¼ 0, 10, 25, 503 /rad2 were evaluated (kq ¼ 0 is the case for when there is no angle bending potential) The parameters in the LeonardeJones potential were sij ¼ 1 and 3 ij ¼ 1 for all pairs of beads. The cutoff distance parameter c for the non-bond interaction was set to 2.0s because a glassy polymer chain model requires attractive components in the LenardeJones potential [12,22,23]. Table 1 lists the resist model parameters. 2.3. Modeling of cross-linking process Fig. 2 shows the modeling flow of the epoxy-based chemically amplified resist. The resist monomer was divided into three groups:

b BS1

O

O

O

O

BS2 HC

HC

HC

HC CH

CH

CH

CH

BS3 BS3

O

O

O

O

BS1 O

O

O

O

Fig. 1. (a) Chemical structure of monomer molecule before exposure and (b) beadespring model of a conventional epoxy-based chemically amplified resist.

H. Yagyu et al. / Polymer 53 (2012) 4834e4842

Fig. 2. Modeling flow of epoxy-based chemically amplified resist.

(1) an epoxy group, (2) a phenyl group connected to methylene, (3) isopropylidene groups and phenyl groups, which were denoted by (1) BS1, (2) BS2, and (3) BS3 beads, respectively. The methylene group was neglected in this model because the beads in the beade spring model have many degrees of freedom and methylene enables the connection between phenyl groups to bend, which is not possible if no intermonomer is present. The catalyst was defined as an S1 bead and it was finally removed in the model. S1 beads are reactive beads that can react with BS1 beads. Table 2 shows the beads employed in the model. In this study, the crosslinked structure of the resist was produced in the molecular dynamics simulation by the following steps: Step 0: Hundred monomer models without the angle bending potential were positioned at random positions in the cell to realize a density of 0.85ms3 and molecular dynamics calculations were performed for an NVT ensemble at a temperature of T ¼ 1.03 /kB for 105 time steps (103s) as the relaxation step. After the relaxation step, the mean square radius of gyration hR2g i of the monomer was 2.5s2. To obtain the equilibrium state, the simulation was performed until the chains moved by at least 2hR2g i. This was realized by monitoring the mean square displacement (MSD) of the molecule during the calculations, which is the general procedure for a polymer model [12]. Step 1: To create a BS1eS1 bond, 50e900 S1 beads of the catalyst were placed at random positions in the cell and molecular dynamics calculations were performed for an NVT ensemble at a temperature of T ¼ 1.03 /kB for 105 time steps (103s). In this calculation, a reaction occurred when two beads approached a within 1.1s in the molecular dynamics simulation using the reaction function of the simulation. This step modeled the UV irradiation step in the experiment. The maximum number of bonds between BS1 and S1 beads was set to 1. BS1

Table 2 Type of beads used in the model. Types

Colors

Functions in the model

BS1 BS2 BS3 S1 BS11

Red Blue Gray Yellow Green

Epoxy group, to be reacted with catalyst Phenyl group connected to methylene Isopropylidene groups, phenyl groups Catalyst with one reactive bond, finally removed Activated epoxy group with two reactive bond

4837

beads were defined as cross-link point in the next step. Theses cross-link points can be numerically determined without using S1 beads. However, using S1 beads randomly dispersed in the cell, the cross-linking points were obtained at positions selected in 3D space during the molecular dynamics calculation. Moreover, the effect of dispersion of catalyst in the resist material on the property can be considered using S1 beads. Step 2: Reacted BS1 beads were converted to BS11 beads. BS11 beads are defined as beads that can react with other BS11 beads to form a cross-link. Subsequently, scission of the BS1eS1 bond was performed at an NVT ensemble temperature of T ¼ 1.03 /kB for 105 time steps (103s) in the molecular dynamics simulation. After the calculation, S1 beads were removed from the cell. Step 3: The BS11eBS11 cross-link was formed at an NVT ensemble temperature of T ¼ 1.03 /kB for 105 time steps (103s) in the molecular dynamics simulation. The parameters about the reaction between BS11 bead to another BS11 are same as that of step 1. After creating the cross-linked structure, the conformation of the chain was evaluated and the angle bending potential was added to the three connected beads in the model using a program developed by the authors. The maximum number of bonds for a BS11 bead was set to 3. Since the BS11 bead in the initial structure has one connection with a BS2 or BS3 bead (see Fig. 1(b)), there are a maximum of two activators in this step. After the calculation, the creation of uniform dispersion of crosslinks was checked by the snapshot of the model. Step 4: Finally, the relaxation process was performed for 2  107 time steps (¼2  105s) using cell size that gave a density of 0.85ms3 and a temperature of T ¼ 0.253 /kB. This step took almost four days of CPU time on a computer with a 2.93-GHz Xeon processor. The temperature in the relaxation process was set to that in the glassy region defined by the evaluated glass transition temperature of the model. The details of glass transition temperature calculation are discussed in Section 2.4. It is assumed that the relaxation for 2  107 time steps is insufficient to relax completely the glassy cross-linked polymer chain if based on MSD of the molecule. But, we carried out relaxation process for 2  107 time steps, which can be obtained initial structure required for estimating elastic modulus in elongation simulation because the models are equilibrium state before the cross-linking process. In addition, the cross-linking points are set at a random position in step 1. The dynamics of the polymer model at the glass state created by the same relaxation process is reported [22]. The effect of relaxation time on the properties will be discussed in future work. Three-dimensional periodic boundary conditions were employed in all simulations. According to above modeling steps, models using 0e900 S1 beads of the catalyst were produced. This resulted in models with different cross-linking ratios. In addition, we produced a model with kq ¼ 0, 10, 25, 503 /rad2 in the angle bending potential and the relationship between the cross-linking ratio and the elastic modulus from the elongation simulation were investigated in this study. In order to evaluate the properties of the resist of one cross-linking ratio, one initial configuration was prepared. We created and estimated many models with different cross-linking ratio for improvement of the reliability of calculations.

2.4. Glass transition temperature calculation method The glass transition temperature Tg is an important parameter when determining the mechanical properties of polymer materials. Since non-cross-linked epoxy-based chemically amplified resists such as SU-8Ô (MicroChem Co.) and TMMRÔ (Tokyo Ohka Kogyo

4838

H. Yagyu et al. / Polymer 53 (2012) 4834e4842

Co. Ltd.) have glass transition temperatures that are about 50  C above room temperature [24], an elongation simulation should be performed below the glass transition temperature of the resist model to simulate elongation behavior in an experiment at room temperature. In this study, Tg was calculated by two methods. The first method is the conventional method used in molecular dynamics simulations of polymers. The glass transition temperatures of the resist model were calculated from the relationship between the volume and the temperature using periodic boundaries and a constant pressure in the Andersen method (NPT ensemble) [25], where the pressure is set to 1.03 /s3. The calculations were performed in the range of temperature 2.0e0.13 /kB at 0.13 /kB interval (20 temperatures for one model) for 105 time steps (103s). The mean volume was calculated from the data for last 3  104 time steps at each temperature. The glass transition temperatures Tg were estimated from the break points in plots of the volume against temperature ðTgNPT Þ. The second method determined the glass transition temperatures from the variation in the MSD for the calculation times for the NVT ensemble ðTgNVT Þ. The glass transition temperatures Tg were estimated from break points in plots of MSD against time [23]. 2.5. Elongation simulation method A uniaxial elongation simulation was performed using a Poisson ratio of 0.3, which is close to the experimental Poisson ratio for a cross-linked resist [26e28], and an elongation rate of 0.03s/s at a temperature below the glass transition temperature of the model, to clarify the relationship between the mechanical properties and the cross-linking ratio of the resist model. The normal stress S in the model was calculated from the microscopic virial tensor, P ¼ Pzz e (Pxx þ Pyy)/2, which x is the direction of elongation. The stresse strain curve of the model was obtained by the elongation simulation, and the elastic modulus was calculated from the gradient of the stressestrain curve below a strain of 0.05.

to obtain the required structure for tensile or bulge tests may seriously affect the measurement results since it may damage the samples or alter their cross-linking ratios. The nanoindentation tests were conducted using a hardness tester (DUH-W201, Shimadzu Corp.) with a Berkovich diamond indenter at room temperature, and the sample was indented with the maximum force of 15 mN with loading rate of 0.284 mN/s and holding time of 60 s [5]. In the nanoindentation test, the elastic modulus of the sample can be calculated by Eq. (9):

pffiffiffiffi Er ¼

p S

2b

pffiffiffiffiffi Ac

(9)

where S is the slope of the unloading curve and known as the contact stiffness of the material, b is a constant depending on the geometry of the indenter, Ac is the projected contact area, and Er is the reduced elastic modulus which is calculated from Eq. (10):

1 1  m2s 1  m2i ¼ þ Er Es Ei

(10)

where Es and ms are the elastic modulus and Poisson’s ratio for the sample. Ei and mi are the Poisson’s ratio and the elastic modulus of the indenter. For a Berkovich indenter, the material properties are often used as Ei ¼ 1140 GPa and mi ¼ 0.07. A commercial 30-mm-thick dry-film negative resist (TMMF S2030, Tokyo Ohka Kogyo Co., Ltd.) was used to fabricate the samples (TMMF is a film of TMMR). Samples with different crosslinking ratios were prepared by varying the UV dose and using PEB temperatures between 90 and 105  C [2,5]. In the experiment, the cross-linking ratio of the resists was determined by Fouriertransform infrared (FT-IR) spectroscopy measurements. FT-IR measurement is a common method for the characterization of organic materials. According to them, the measurement of the variation of the absorption by the epoxy group (910 cm1) during the UV exposure and the PEB provides a means for evaluating the relationship between the UV exposure, PEB conditions, and the cross-linking reaction [31].

3. Nanoindentation test A nanoindentation test was performed to experimentally evaluate the elastic modulus since it can better determine the mechanical properties of small volumes on a substrate than conventional mechanical tests such as tensile and bulge tests [29,30]. It has a further advantage in that it requires minimal preparation for samples with a semi-cross-linking ratio, the additional processing (e.g., etching) required for negative resist samples

4. Results and discussion 4.1. Cross-linking structure Fig. 3 shows snapshots of the non-cross-linked resist model (Nv ¼ 0) and the cross-linked resist model (Nv ¼ 900) after the relaxation step. The green (light gray in black and white) beads (in

Fig. 3. Snapshots of (a) non-cross-linked resist model and (b) cross-linked resist model with Nv ¼ 900.

H. Yagyu et al. / Polymer 53 (2012) 4834e4842 Table 3 Calculated cross-linking ratio of the model. M is the number of molecules in the system. Nv

M

Nb

Nc

n

0 50 100 200 300 400 450 500 550 600 650 700 750 800 900

100 74 52 10 7 3 1 1 1 1 1 1 1 1 1

0 50 100 200 300 400 450 500 550 600 650 700 750 789 800

0 44 93 193 295 395 442 491 539 591 638 693 742 780 788

0 0.055 0.116 0.241 0.369 0.494 0.553 0.614 0.674 0.739 0.798 0.867 0.928 0.975 0.985

the web vesion), which are the cross-linked beads, are uniformly dispersed in the cell. To investigate the cross-linked structure in detail, the cross-linking ratio of the resist model n, was calculated using the following equation:

v ¼

Nc  Nd ðNa þ Nb Þ=2

(11)

where Na is the number of end beads (BS1), Nb is the number of activated end beads (BS11), Nc is the number of cross-linked bonds (BS11eBS11), which was obtained from the simulation results, and Nd is the maximum number of bonds of the activated end beads. The calculated n represents the experimental cross-linking ratio, which is the ratio of epoxy groups with opening rings. Table 3 lists the calculated cross-linking ratios of the resist model with Nd ¼ 2 and Na þ Nb ¼ 800. The cross-linking ratio increased with increasing number of catalyst beads Nv, the molecule was a single molecule from n ¼ 0.494 (Nv ¼ 400). Fig. 4 shows the cross-linking ratio dependencies of the glass transition temperature Tg of the resist model with an angle spring constant kq ¼ 253 /rad2 in the angle bending potential. The glass transition temperature Tg of the non-cross-linked model (i.e., crosslinking ratio, n ¼ 0) was approximately 0.63 /kB. Tg increased with increasing cross-linking ratio and became constant value at a crosslinking ratio of about n ¼ 0.6. In the other model with spring

1.0

constants kq ¼ 0, 10, 503 /rad2 in the angle bending potential, the glass transition temperature Tg of the non-cross-linked model were determined to be 0.4 (kq ¼ 0), 0.5 (kq ¼ 103 /rad2), and 0.73 /kB (kq ¼ 503 /rad2). Therefore, elongation simulations were performed at temperatures of 0.253 /kB, which is sufficiently below the glass transition temperatures of all resist models to realize elongation behaviors of the glassy polymer. 4.2. Scaling of the model To compare the elastic modulus obtained by the simulation with experimental values, the length s and energy 3 were determined as follows. First, s was mapped using the mean square radius of gyration hR2g i of non-cross-linked molecules of the coarse-grained models and that of a fully atomistic molecular dynamics model with the GAFF (general AMBER force field) [32] for the bond potential and the LennardeJones potential for the non-bond potential. In this fully atomistic model, the charge of the atoms was calculated using density functional theory using the B3LYP/631Gd method with GAMESS [33] and Ewald summation was used to include long-range Columbic interaction between point charges [34]. An amorphous structure of fully atomistic resist model with 100 monomers was obtained using three relaxation steps. In the first step, the model was calculated using an NPT ensemble and a pressure of 100 MPa for 0.1 ns. The next step used an NPT ensemble and a pressure of 0.1 MPa for 0.1 ns. The final step was performed for an NVT ensemble and a pressure of 0.1 MPa. All the calculations were performed at a temperature of 300 K. From the final model, the mean square radius of gyration hR2giFL of the fully atomistic model was defined as 0.598 nm2. In addition, the mean square radius of gyration hR2giCG of the coarse-grained model after relaxation at the temperature of 0.253 /kB was calculated to be 2.55s2 for kq ¼ 0, 2.84s2 for kq ¼ 103 /rad2, 2.81s2 for kq ¼ 253 /rad2, and 2.88s2 for kq ¼ 503 /rad2. Comparison of these simulated values for the coarse-grained model with that of the fully atomistic model (hR2g iFL ¼ hR2g iCG) gives s ¼ 0.484 nm (kq ¼ 0), s ¼ 0.459 nm (kq ¼ 103 / rad2), s ¼ 0.461 nm (kq ¼ 253 /rad2), and s ¼ 0.456 nm (kq ¼ 503 / rad2). Next, in the elongation simulation, room temperature (300 K) was set to 0.253 /kB and 3 was calculated to be 3 ¼ 1.63  1020 J using kB ¼ 1.38  1023 J/K and the calculated temperatures of the models in the molecular dynamics simulation. Table 4 shows the calculated unit parameters of the coarse-grained resist models. 4.3. Elongation simulations Fig. 5(a) shows stressestrain curves obtained in the uniaxial elongation simulations of the resist model with an angle bending spring constant of kq ¼ 253 /rad2 in the angle bending potential at a temperature of 0.253 /kB. The stress increased with increasing cross-linking ratio. Fig. 5(b) shows stressestrain curves obtained in the uniaxial elongation simulation of the resist model with an angle bending spring constant of kq ¼ 0 in the angle bending potential at a temperature of 0.253 /kB. The stress was constant in all models in the low strain region below the experimental breaking strain of about 0.05 [24]. These results demonstrate that the angle bending potential strongly influences the stress response.

Tg NPT

0.9

4839

Tg NVT

0.8 0.7 0.6

Table 4 Calculated mean square radius of gyration hR2giCG, the unit of length s, the unit of energy 3 , and the calculated temperature T of the coarse-grained resist models.

0.5 0.0

0.2

0.4

0.6

0.8

1.0

Cross-linking ratio Fig. 4. Dependency of the glass transition temperature of the model (kq ¼ 253 /rad2) on the cross-linking ratio. TgNPT was calculated with NPT ensemble (C) and TgNVT was calculated with NVT ensemble (B).

kq [3 /rad2]

hR2g iCG [s2]

s [nm]

T [3 /kB]

3

0 10 25 50

2.55 2.84 2.81 2.88

0.484 0.459 0.461 0.456

0.247 0.245 0.241 0.247

1.68 1.72 1.68 1.66

[J]    

1020 1020 1020 1020

4840

H. Yagyu et al. / Polymer 53 (2012) 4834e4842

a

at the slower strain rate, rubber-like stress response is observed. The average strain rate l from the initial state to 20% strain state is 2  103s1 for elongation rate of 0.03s/s, while the relaxation time se for the non-cross-linked model at the temperature of T ¼ 1.03 /kB (>Tg) is estimated to be about 2  102s (1/se ¼ 5  103s1) from Rouse relaxation time [12]. At the temperature of T ¼ 0.253 /kB (
Fig. 5. Stressestrain curves for resist models with (a) kq ¼ 253 /rad2 and (b) kq ¼ 0 in the angle bending potential.

The elongation rate was estimated using the relation of

s ¼ s(m/3 )1/2. Where, the bead is approximated with phenyl group, and the mass m was defined as 78.11 amu. From this definition, the elongation rate of 0.03s/s (strain rate of 2  103s1) was calculated as 5.8  105 mm/min. This value is 105 times as fast as that of 5 mm/ min used in the conventional tensile test [17]. The elastic modulus of the model with the cross-linking ratio of 0.985 at the elongation rate of 3  105s/s (strain rate of 2  106s1), which is the same order rate as that of the experiment, was calculated to 4.3 GPa. This value was 12% smaller than that in the experimental tensile test, which was reported as about 4.9 GPa [17]. But, in this paper, the elastic modulus measured with nanoindentation in the experiment. In general, the elastic modulus measured with nanoindentation resulted in slightly large in comparison with that measured with tensile test, and the reason of this discrepancy was not clear. It is assumed that the reason is that the nanoindentation test was carried out under the condition of stress control and the calculation method of the elastic modulus differs from that with tensile test. In future work, we must calculate the elastic modulus using indentation model. In this paper, we used the elongation rate of 0.03s/s, which is the condition that the elastic modulus of the model with cross-linking ratio of 0.985 corresponds with the value of the sample of cross-linking ratio of 1 measured with nanoindentation. In the elongation simulation, the stress depends on the strain rate. If the strain rate is fast in comparison with the relaxation rate of the model, glass-like stress response is observed. In contrast,

4.5. Effect of angle potential In Section 4.3, the angle bending potential was shown to strongly influence the stress response. In this section, the role of the angle bending potential in the elongation simulation of the resist model is discussed by analysis of varying the bond and angle potentials of the model because the stress in this molecular dynamics simulation was calculated using virial stress, which is

8

6

E [GPa]

b

In the experiment, the elastic modulus in low cross-linking ratio region cannot be measured because the samples are too soft for measurement. Therefore, for analysis of variation of dynamics, it is important to simulate the properties in the region which cannot be evaluated with experiment, and evaluate the cross-linking ratio dependence. Fig. 6 shows the dependencies of the elastic modulus on the cross-linking ratio of the resist model with kq ¼ 0, 10, 25, and 503 /rad2 of the angle bending spring constant in the angle bending potential. The elastic modulus E of the model with kq ¼ 25 and 503 / rad2 increased when the cross-linking ratio n increased from about 0.3, the model with the larger angle bending spring constant induced a larger elastic modulus at a large cross-linking ratio (n ¼ 1). These results demonstrate that the elastic modulus of the model with a large angle bending spring constant in the angle bending potential is correlated with the cross-linking ratio. Fig. 7 compares the experimentally measured elastic modulus in the cross-linking ratio range 0.4e1.0 obtained from nanoindentation tests [5] with that obtained by the simulation. The tendency of the elastic modulus obtained by the simulation with kq ¼ 253 /rad2 was in good agreement with that observed in the experiment.

4

2

0 0

0.2

0.4

0.6

0.8

1

Cross-linking ratio Fig. 6. Dependencies of the elastic modulus on the cross-linking ratio of the resist model.

H. Yagyu et al. / Polymer 53 (2012) 4834e4842

Fig. 7. Dependencies of the elastic modulus on the cross-linking ratio in the resist model (kq ¼ 0, 10, 25, and 503 /rad2) and experiment (PEB 90, 95, 100, and 105  C).

defined as the sum of the kinetic and potential energies. In this section, the ratio of the potential energy to initial was plotted to evaluate the motion of beads during elongation. Fig. 8(a) shows the of the rate of increase of the angle bending potential UANGLE/UANGLE 0 resist model with n ¼ 0.985 (Nv ¼ 900) for kq ¼ 0, 10, 25, and 503 / is the angle bending rad2 during elongation simulation. UANGLE 0 potential energy at a strain of zero and UANGLE is the angle bending potential energy at a non-zero strain. The angle bending potential

a

4841

were the sum of the energy calculated using Eq. (5) using kq ¼ 253 / rad2 and the angles of all three connected beads of the model. Although the model with kq ¼ 0 shown in Fig. 8(a) does not have the angle bending potential, the angle bending potential of the model with kq ¼ 0 was calculated as a pseudopotential using Eq. (5) with kq ¼ 253 /rad2 to clarify the effect of the angle potential on the bead parameters such as the bond length during elongation. UANGLE/ of the model without angle bending potential (kq ¼ 0) UANGLE 0 increased linearly with increasing the strain. In contrast, the angle bending potential of the model with kq ¼ 10, 25, 503 /rad2 increased of the model with from the strain of about 0.05, but UANGLE/UANGLE 0 kq ¼ 25 and 503 /rad2 decreased up to the strain of about 0.05. These reductions are considered to be due to the average angle between the beads at a strain of zero is 107.4 for kq ¼ 253 /rad2 and 107.1 for kq ¼ 503 /rad2, which are small compared with the equilibrium angle of 109.5 . From Fig. 8(a), it was confirmed that the effect of the angle bending potential on the elongation becomes strong in the large strain region. In the elongation simulation of polymer, Yashiro et al. pointed out that the angle bending potential influenced the chain conformation and the stress of model [35]. The result in this paper was in agreement with their opinions. Fig. 8(b) of the shows the rate of increase of the bond potential UFENE/UFENE 0 resist model of n ¼ 0.985 (Nv ¼ 900) with kq ¼ 0, 10, 25, and 503 /rad2 in the angle bending potential during elongation simulation. UFENE 0 is the bond potential energy at a strain of zero and UFENE is the bond potential energy at a non-zero strain. The bond potentials were the sum of the energy calculated using Eq. (2) using the bond length of of the models increased with all connected beads. UFENE/UFENE 0 increasing the angle bending spring constant, and the model with larger angle bending spring constant showed larger bond potential energy. The stress in Fig. 5 is related to potential energy, but slightly differs from Fig. 8 because of including of kinetic energy and non-bond potential. When the angle bending potential energy was considered, the stress of the model increased in comparison with that of model with kq ¼ 0 because the initial total potential energy of the model increased. Based on these results, for the model with kq ¼ 0 in the angle bending potential, the bond length was assumed to be constant and the size of the cross-linked mesh surrounding by the cross-links was assumed not to change during elongation. On the other hand, for the model with a relatively large angle bending potential, the addition of angle increases the bond length because the angle of the beads was controlled to the equilibrium angle of 109.5 , during the elongation simulation. It was assumed that the relation of the bond potential and the angle bending potential is a tradeoff, and the relation of the variation of the bond length and the degree of constraint of angle depends on the elastic modulus of the resist. 5. Conclusions

Fig. 8. Dependencies of the rate of increase of (a) angle bending potential energy and (b) bond potential energy on the strain. The angle bending potential of the model with kq ¼ 0 and n ¼ 0.985 is the pseudopotential calculated from the angles between beads and kq ¼ 253 /rad2.

A coarse-grained molecular dynamics model for an epoxy-based chemically amplified resist was produced using a simple beade spring model. Models with different numbers of cross-linking ratios were used in an uniaxial elongation simulation. The uniaxial elongation simulations results reveal that at low strains below the experimental breaking strain of the resist, the elastic modulus obtained using the KremereGrest model extended angle bending potential depends on the cross-linking ratio. The dependency of the model with kq ¼ 253 /rad2 in the angle bending potential agreed well with that obtained by experiments. The fact indicated that the averaged structure of the resist molecules during deformation can be showed by the model with the angle bending spring constant of kq ¼ 253 /rad2. This confirms the validity of the developed simulation model and its applicability for epoxy-based chemically amplified resists.

4842

H. Yagyu et al. / Polymer 53 (2012) 4834e4842

The mechanical properties of the model depend on the potential energy (bond, angle, non-bond), temperature, and elongation rate. To calculate the elastic modulus of the resist, the simulation was carried out at enough low temperature in comparison with glass transition temperature of the model, and the elongation rate at which the stress in glassy region is realized. With respect to the potential energy, the angle bending potential was optimized since the bond and non-bond potentials used in this model are same as that of KremereGrest model except for cutoff distance of Lenarde Jones potential, which can be realized polymer dynamics. This optimization enabled the prediction of the elastic modulus measured with nanoindentation. In future research, we intend to use the developed beadespring model to predict the filtration ability of the resist (such as the effect of the molecule length of the solvent and the interaction between the resist and solvent) and the diffusion rate of short length molecules in the resist model such as water and organic solvents in the polymer. Acknowledgement This work was supported by JSPS Grant-inAid for Young Scientist (B) (23710148). The authors would like to thank N. Taneichi, M. Hagihara, and H. Shinbori of Tokyo Ohka Kogyo Company Ltd. for helpful discussions about resist materials and Prof. M. Hojo of Kyoto University for the many valuable discussions regarding the simulation results.

References [1] Hirai Y, Sugano K, Tsuchiya T, Tabata O. J Micromech Microeng 2010;20: 065005. [2] Hirai Y, Sugano K, Tsuchiya T, Tabata O. IEEE/ASME J Microelectromech Syst 2010;19:1058e69. [3] Hirai Y, Nakai Y, Makino Y, Sugano K, Tsuchiya T, Tabata O. Epoxy-based permeable membrane fabrication for 3D microfluidic device. In: Proceeding of (NEMS2011): 6th IEEE international conference on nano/micro engineered and molecular systems, Kaohsiung; 2011. p. 188e91.

[4] Napoli M, Eijkel JCT, Pennathur S. Lab Chip 2010;10:957e85. [5] Hirai Y, Uesugi A, Makino Y, Yagyu H, Sugano K, Tsuchiya T, et al. Negativephotoresist mechanical property for nano-filtration membrane embedded in microfluidics. In: Proceeding of transducers2011: 16th international conference on solid-state sensors, actuators and microsystems, Beijing; 2011. p. 2706e09. [6] Wu C, Xu W. Polymer 2006;47:6004e9. [7] Clancy TC, Frankland SJV, Hinkley JA, Gate TS. Polymer 2009;50:2736e42. [8] Varshney V, Patnaik SS, Roy AK, Farmer BL. Macromolecules 2008;41:6837e42. [9] Nouri N, Ziaei-Rad S. Macromolecules 2011;44:5481e9. [10] Yagyu H, Utsumi T. Comput Mater Sci 2009;46:286e92. [11] Grest GS, Kremer K. Macromolecules 1990;23:4994e5000. [12] Kremer K, Grest GS. J Chem Phys 1990;92:5057e86. [13] Kroger M, Makhloufi R. Phys Rev E 1996;53:2531e6. [14] Kroger M, Luap C, Muller R. Macromolecules 1996;30:526e39. [15] Kroger M. Phys Rep 2004;390:453e551. [16] Aoyagi T, Sawa F, Shoji T, Fukunaga H, Takimoto J, Doi M. Comput Phys Commun 2002;145:2676e9. [17] Lorenz H, Despont M, Fahmi N, LaBianca N, Renaud P, Vettiger P. J Micromech Microeng 1997;7:121e4. [18] Shaw JM, Gelorme JD, LaBianca NC, Conley WE, Holmes S. IBM J Res Dev 1997; 41:81e94. [19] Misumi K, Saito K, Yamanouchi A, Senzaki T, Honma H. J Photopolym Sci Technol 2006;19:57e62. [20] Zhang J, Chan-Park MB, Li CM. Sens Actuators B 2008;131:609e12. [21] Carraher Jr CE. Giant molecules. 2nd ed. New Jersey: John Wily & Sons; 2003. p. 33. [22] Bulacu M, van der Giessen E. Phys Rev E 2007;76:011807. [23] Morita H, Tanaka K, Kajiyama T, Nishi T, Doi M. Macromolecules 2006;39: 6233e7. [24] Feng R, Farris RJ. J Micromech Microeng 2003;13:80e8. [25] Andersen HC. J Chem Phys 1980;72:2384e93. [26] Luo C, Schneider TW, White RC, Currie J, Paranjape M. J Micromech Microeng 2003;13:129e33. [27] Feng R, Farris RJ. J Mater Sci 2002;37:4793e9. [28] Dellmann L, Roth S, Beuret C, Racine GA, Lorenz H, Despont M, et al. Sens Actuators A 1998;70:42e7. [29] Oliver WC, Pharr GM. J Mater Res 1992;7:1564e83. [30] Al-Halhouli AT, Kampen I, Krah T, Büttgenbach S. Microelectronic Eng 2008; 85:942e4. [31] Yoshino H, Matsumoto H. Jpn J Appl Phys 1992;31:4283e7. [32] Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. J ComputChem 2004;25: 1157e74. [33] Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, et al. J Comput Chem 1993;14:1347e63. [34] Allen MP, Tildesley DJ. Computer simulation of liquids. Oxford: Clarendon Press; 1987. [35] Yashiro K, Naito M, Ueno S, Feng J. Int J Mech Sci 2010;52:136e45.