Numerical modeling and analysis of plasmonic flying head for rotary near-field lithography technology

Rotary near-field lithography (RNFL) technology provides a route to overcome the diffraction limit with a high throughput and low cost for nanomanufacturing. Utilizing the advantage of the passive flying of a plasmonic head, RNFL can achieve a 10 m/s processing speed with a perfect near-field condition at dozens of nanometers. The flying performance of the plasmonic flying head (PFH) is the pivotal issue in the system. The linewidth has a strong correlation with the near-field gap, and the manufacturing uniformity is directly influenced by the dynamic performance. A more serious issue is that the unexpected contact between the PFH and substrate will result in system failure. Therefore, it is important to model and analyze the flying process of the PFH at the system level. In this study, a novel full-coupled suspension-PFH-air-substrate (SPAS) model that integrates a six-degree of freedom suspension-PFH dynamics, PFH-air-substrate air bearing lubrication, and substrate vibration, is established. The pressure distribution of the air bearing is governed by the molecular gas lubrication equation that is solved by the finite element method (FEM) with a local pressure gradient based adaptive mesh refinement algorithm using the COMSOL Multiphysics software. Based on this model, three designs of the air bearing surface are chosen to study the static, dynamic, and load/unload performance to verify whether it satisfies the design requirements of RNFL. Finally, a PFH analysis solver SKLY.app is developed based on the proposed model.


Introduction
Optical lithography has been extensively performed in the semiconductor industry for decades. Currently, to sustain Moore's law [1], the critical dimension is a 10-nm node or even 7 nm node according to the International Technology Roadmap for Semiconductors (ITRS); this leads to a substantial increase in the cost [2]. In particular, the fabrication of the high-quality lithography masks is expensive and difficult with multiple-lithography. Focused ion beam (FIB) and electron-beam lithography (EBL) offer the ability to attain the resolution of less than 10 nm; however, the low throughput of these methods owing to their slow scanning ability, restrains their applications for nanofabrication on a large scale [3,4]. Near-field lithography (NFL) is a newly emerging technique that has been explored for realizing ultrahigh resolution beyond the Rayleigh diffraction limit with a specially designed plasmonic focusing lens and without image aberration under an extremely short work distance [5−8]. Lee et al. [9] used an active nanogap control implementation in which a metallic sharp ridged nanometer-scale aperture was fabricated on a solid immersion lens (SIL), and the nanometer-scale aperture excited the surface plasmon polaritons (SPPs) and focused the light on the center of the SIL. In 2011, the Zhang's group [10,11] achieved maskless lithography with a resolution of 22 nm by rotary near-field lithography (RNFL) on an inorganic thermal photoresist TeO x film. Recently, we performed RNFL on a new molecular glass chemically amplified resist (CAR) that has a higher sensitivity and resolution than a thermal photoresist and can be easily applied in the industry without additional cost [12,13]. The system scheme is illustrated in Fig. 1. The plasmonic lens is fabricated on a transparent slider and flies above a substrate coated with the photoresist layer. An array of plasmonic lens, forming the plasmonic flying head (PFH), generates a self-adaptive near-field gap nf ( ( , )) d x y of tens of nanometers between the PFH and high speed revolving substrate that is insensitive to the environmental vibration. The plasmonic lens can focus the incident laser beam to a spot size of sub-20 nm with enhanced field intensity by exciting the surface plasmon polaritons.
The flying performance of the PFH is key to RFNL [12]. Owing to the exponential decay of the surface plasmon in the near-field range, the near-field gap of the PFH will directly determine the laser light absorption of the photoresist that further influences the lithography linewidth. In addition, the dynamic performance of the PFH will affect the manufacturing uniformity. Furthermore, it is difficult to fly on the organic CAR because of its poor mechanical properties; specifically, the unexpected contact between the PFH and substrate will result in system failure. Therefore, modeling and analysis of the flying of the PFH is significant for a reliable RNFL design. The principle of PFH is similar to that of the magnetic head in hard disk drives (HDDs). In the last few decades, the head-disk interface (HDI) has been extensively studied numerically and experimentally along with the development of the HDD industry [14−17]. A small gap of the order of 5-20 nm, which is much smaller than the mean free path of the air molecule, is maintained between the disk and slider in an HDD; therefore, the continuum model is no longer valid in this region [18]. The generalized Reynolds equation is introduced by modifying the traditional Reynolds equation with the slip theory [19,20] or linearized Boltzmann equation [21] to solve the problem. Among these, Fukui and Kaneko's molecular gas lubrication (MGL) equation is widely used in the HDD industry because of its good convergence to the Boltzmann equation [22,23]. The finite difference [24,25], finite volume [26,27], and finite element (FE) methods [28,29] are implemented to solve the generalized Reynolds equation. For the FE method, use of a structured mesh is not required, so that it is highly suitable for a complex geometry and particularly for local adaptive refinement in regions where the solution displays large gradients of the solution parameter or the slider geometry changes. The slider equilibrium attitude is determined by the air pressure and applied suspension forces. The coupling between the two sets of equations can be solved by two approaches. The solution can be found either by the transient solution approach until the steady-state conditions are found, or by directly solving the steady coupled equations. Several models have been established to analyze the dynamic performance of the magnetic head flying [30,31]. However, the PFH in an RNFL system has a different system setup and design goals, so a new model is required.
In this paper, we proposed a novel full-coupled suspension-PFH-air-substrate (SPAS) model that integrates a six-degrees of freedom (6-DOF) suspensionhead dynamics, head-air-substrate air bearing lubrication, and substrate vibration for RFNL. The pressure distribution of the air bearing is governed by the MGL equation that is solved by the FE method with a local pressure gradient-based adaptive mesh refinement algorithm using the COMSOL Multiphysics software. The Newton-Raphson method is utilized to solve the slider equilibrium attitude by neglecting the inertia term and damping term of the proposed model. Next, three designs of the air bearing morphology are studied to examine the static, dynamic, and load/ unload performance to verify whether it satisfies the requirements of the PFH in RNFL. Finally, a PFH analysis solver, SKLY.app is developed based on the proposed model.

Suspension-PFH-air-substrate system modeling
The PFH in RNFL utilizes the concept of the air bearing head in the HDD. An air bearing surface (ABS) with a suitably designed morphology exerts an aerodynamic force to support the head to fly on a nanoscale distance with a high-speed rotary substrate that can provide a near-field condition for the lithography. An actual PFH system in RNFL has several components such as the main suspension, gimbal, limiter, head, substrate, as well as ambient air. The different parameters of these components will result in different static head postures and system dynamic performances that are the key lithography quality evaluation indices of RNFL. Therefore, performing a numerical simulation of the PFH system is usually necessary because it can be used to study the different factors influencing the static and dynamic performances, and furthermore to assist in the design of the head and suspension parameters to satisfy the demand of NFL.
We established a fully coupled suspension-PFH-airsubstrate (SPAS) model. Figure 2 shows the schematic of the different states of the suspension-head-substrate system and its corresponding dynamic model. The 6-DOF model divides the suspension into three parts, namely, the tab, local suspension, and main suspension, and the acting force between these parts is simplified by the spring and damping. The interaction between the suspension and head consists of three parts: the dimple/head contact force, gimbal force/moments, and limiter contact force, when the dimple separation is larger than the limiter gap. The head has three DOF with corresponding posture parameters of near-field gap, pitch angle, and roll angle. The aerodynamics force, head/substrate contact force, and intermolecular force act on the head/substrate interface (HSI).
To model the PFH system, we establish the following 6-DOF dynamic motion equations that can describe the complete process of the PFH system: are the generalized mass and displacement vectors, respectively, of the tab, local suspension, main suspension, and head, which are the corresponding 6-DOF parameters in the PFH system. Combining the different states, we can obtain the stiffness matrix, damping matrix, and external force vector, where i k and The Reynolds equation is usually used to calculate the air pressure distribution at the HSI. However, because of the significant rarefied effect when the spacing between the head and substrate is extremely small, such that it is much less than the mean free path of air molecule, the continuum model is not applicable. The Reynolds equation has been modified in the last few decades to calculate the pressure distribution in such severe conditions. The modified equation can be expressed in the following dimensionless generalized form: where H is the dimensionless air bearing thickness;  is the bearing number defined as where  is the ambient gas viscosity, U is the dimensionless gas velocity, and  is the squeeze number defined as  is the angular velocity; and Q is a flow rate coefficient to account for the gaseous rarefaction effect which has been provided by Fukui and Kaneko's model derived from the linearized Boltzmann equation [21]. After obtaining the pressure distribution, the aerodynamics force and moments can be calculated by where a p is the atmospheric pressure, air F x and air F y are the coordinates of the air bearing force center, and G x and G y are the coordinates of the gravity center.
For the head/substrate contact force, Chang et al. [32] proved that for extremely smooth surfaces the contact is mostly elastic throughout. Therefore, here we use the Greenwood and Williamson (GW) model [33], the basic elastic rough surface contact model, to calculate the head/substrate contact force.
where  csd E is the equivalent stiffness of the head and substrate;  , r, and  are the density of the asperities, asperity radius of the curvature, and standard deviation of the asperities heights, respectively; ( ) s  is the dimensionless asperity height distribution function; and csd F x , csd F y are the coordinates of the contact force center.
For the intermolecular force, with the expression derived from a plane-plane interaction, the total force and moments owing to the intermolecular interaction between the head and substrate are obtained by In the load/unload process, the contact states at the dimple/head, limiter/suspension, and tab/ramp interfaces will change, causing the PFH system to have several states. The limiter contact is defined as a rigid contact. The limiter force is defined as where li k is the limiter stiffness and gap z is the initial distance between the suspension and limiter. The Hertzian contact model is used to calculate the dimple/head contact force and ramp/tab contact force. The contact forces are defined as where  cls E is the equivalent elastic modulus of the dimple and head, and cls R and crt R are the equivalent radius of the curvatures of the dimple and tab, respectively.

Numerical methodologies
The most challenging step in the numerical simulation process is to solve the generalized Reynolds equation. In static analysis, there are two problems: one is the forward problem to determine the pressure distribution for the given attitude parameters of the head, and the other is the inverse problem to determine the static attitude parameters for the given external forces and moments. In dynamic analysis, the problem is establish an approach to couple the transient generalized Reynolds equation and motion differential equations to simulate the system dynamic performance. With the Galerkin method, the partial differential equation (Eq. (5)) is multiplied by test function v and integrated over domain  to obtain Processing the equation with integration by parts, along with the boundary condition P = 1 over  , we can obtain the weak form, In this study, Eq. (13) was solved by using the commercial software COMSOL Multiphysics based on the FEM. In the simulation, the Newton-Raphson iterative method was chosen to solve the nonlinear problem, and the global system of equations obtained in each iteration step were solved by the Multifrontal Massively Parallel Sparse Direct Solver (MUMPS).
For the forward problem, the time-dependent terms in Eq. (5) were neglected to obtain the steady pressure distribution at the given head postures. For the inverse problem, we adopted the Newton-Raphson method to calculate the force equilibrium equation that is Eq. (1) but without the inertia terms and damping terms at a given pre-force and moment. We solved the forward problem in each iteration step to obtain the head 3-DOF attitudes. The contact force and intermolecular force were calculated at every element with the given near-field gap distribution. For the dynamic problem, we first calculated the substrate dynamic response, and then coupled the response with the 6-DOF suspension-substrate model to derive the head dynamic response. Although the load acting on the substrate is not axisymmetric, there is a minor nonaxisymmetrical effect on the substrate dynamic response. Equation (5) was calculated transiently with the backward difference method in the time direction, and the Runge-Kutta method was adopted to solve the Eq. (1). Simultaneously solving the equations, we can obtain the dynamic performance of the head. The relative discontinuity of the head geometry results in a large gradient of the pressure distribution; therefore, the meshes for the simulation were generated by the adaptive mesh refinement algorithm to achieve accurate results.

Static analysis
As previous studies have shown, for an RNFL system (Fig. 1), the intensity profile of the modified Bull's eye plasmonic lens (PL), which is characterized by the full width half maximum (FWHM), emitted from the plasmonic lens on the photoresist expands with increase in the off-plane distance, but does not change much when the off-plane distance is in the range of 0-30 nm [34]. The exponential decay of the intensity with the increase in the off-plane distance also increases the demand for the sensitivity of the photoresist. However, the actual focus spot size is determined by not only the near-field gap but also the PL design. Ji et al. [35] shows a lower FWHM can be achieved by an optimal plasmonic lens design. It should be noted that the decrease in the near-field gap would generate a series of problems such as unexpected contact and wear, which would cause an early failure of the system. Therefore, the steady near-field gap at the plasmonic lens position PL ( ) d must be carefully designed. In this work, the design target of PL d is set as 0-30 nm. In addition, to ensure a manufacturing uniformity for the substrate, PL d is expected to be as robust as possible for different values of the substrate radius and preload forces. The ABS must have a relatively large area for the array of PLs to enable parallel writing and highthroughput.
Here we have used three ABS designs modified from HDD designs, as shown in Fig. 3(a). Head 1 is a simple design from the triple pad design. Head 2 is a five-pad head design in Ref. [36] designed in the computer mechanics laboratory (CML). Head 3 is a head of a commercial design used in the HDD industry. Different colors represent different recess depths. The three designs are chosen and evaluated to determine whether they satisfy the RNFL system requirement. A large area at the trailing edge of the array of the PLs is present in all the three designs. The position of the plasmonic lens is also shown in Fig. 3.
We first implemented the forward problem to test the simulation model. The operational parameters for the heads are fixed at a nominal nf d (defined as the clearance at the trailing edge center) of 10 nm, pitch of 150 rad, roll of 1 rad, and radius of 18 mm with a 0° skew. The rotation speed of the substrate disk is set as 5,400 rpm. Figures 3(b) and 3(c) show the adaptive meshes based on the local air pressure gradient and the pressure distributions of the three heads, respectively. The heads are mainly supported by the high-pressure peaks generated by the central trailing pads, maintaining the fly height of the air bearings. The sub-ambient pressure zone in the central regions helps to improve the entire air bearing stiffness and decrease the sensitivity of the substrate speed or other parameters. The side-pad pressure assists in maintaining the stiffness along the rolling direction.
The computation accuracy and efficiency of the adaptive mesh algorithm are also studied because they are important for the PFH analysis. Figure 4 shows the result, where L is the ABS length and l is the maximum size of the initial uniform triangle mesh. Mesh parameter L/l represents the initial mesh intensity. Different levels of adaptive mesh refinement are tested. With the increase in the initial mesh intensity, the accuracy will have a significant improvement accompanied with a slight increase of computing time. We find that the accuracy tends to be stable, whereas the computation time increases linearly with the increase in the initial mesh intensity. With the increase in the adaptive level, the computation accuracy will be improved evenly, though at the price of a longer computing time because the computing time increases much faster at larger mesh parameters. Hence, to ensure the accuracy of the model and make A preload is generated by the bending of the suspension, which mainly determines the static attitude. As shown in Fig. 5(a), heads 1, 2, and 3 are designed to provide a near-field condition for the plasmonic lens at ~30 nm, ~5 nm, and ~15 nm, respectively. The fluctuation in PL d is significant for head 1, which may induce a nonuniformity in the substrate at different radii. PL d of head 2 is proportional to the radius with a small slope. Head 3 has a minor fluctuation in the short radius range and then remains almost constant from 15-30 mm. The reason PL d of heads 1 and 3 does not monotonically increase with the radius is because of the strong negative pressure effect. For the two heads, when the linear speed is sufficiently high, the negative pressure is so strong that the head changes its flying height, pitch angle, and roll angle simultaneously to maintain its balance. In such cases, the flying height no longer exhibits a simple monotonic change with the linear speed and may decrease with the increase in the radius. Figure 5(b) illustrates that the preload has a significant influence on PL d ; with an increase in the preload, PL d decreases because the head must balance the load with the air bearing force by a smaller near-field gap. Head 2 shows a relatively low sensitivity to the preload. Figure 5(c) exhibits that PL d is proportional to the PSA, and head 1 is more sensitive to the PSA than the other two heads. The results of the PSA and preload analysis can help to control PL d according to different parameters of the suspension. Based on these calculation results, it can be concluded that head 2 is the best, whereas head 1 is the worst among the three designs in terms of the static performance.

Dynamic analysis
In the operation of RNFL, the system will experience 450 Friction 6(4): 443-456 (2018) | https://mc03.manuscriptcentral.com/friction different shock events from the environment. Such shocks may directly cause damage to the PFH and photoresist owing to the severe contact between the head and substrate surface. In addition, the PL d fluctuations during the shock events result in an instability of the near-field condition that determines the width of the lithography. Therefore, it is required to rapidly and accurately simulate the operating shock response to obtain a better design of the PFH system. In our study, a shock is modeled as a half sine acceleration wave defined by its peak amplitude and pulse width.
In Fig. 6, we display the free substrate and 6-DOF suspension-PFH system shock response, without the air bearing effect, to a 100G shock of different pulse widths and power spectra. The substrate FE model mesh and response at the outer diameter (OD) and middle diameter (MD) are shown in Fig. 6(a). It has been shown in various studies that the shock response of a rotating substrate to an axisymmetric shock is primarily composed of the first axisymmetric (umbrella) mode, first radial mode, and first axisymmetricradial coupled mode. From the power spectra of the substrate shock in Fig. 6(c), we can observe that the corresponding mode frequencies are 12.4 kHz, 74.02 kHz, and 188.64 kHz that are all excited in the 0.2 ms pulse case. Figure 6(b) displays the displacements of different components of the PFH system and dimple-head contact force during the shock. We can see that the motion of the head from  the green lines that follow the suspension but with a local vibration during the shock, which is owing to the change in the dimple contact state. The primary mode of the suspension at all the pulse widths is the main suspension mode with a frequency of 170 Hz, as shown in power spectra in Fig. 6(d). With the increase in the shock pulse, the peak value of the displacement increases. Next, we consider the case where the acceleration shock is applied to the entire SPAS system uniformly.
Here, a positive shock is defined as one that causes the substrate to initially move toward the head, and both move in the same direction. The suspension is usually designed to be more flexible than the substrate; therefore, the head will move far away from the substrate and the clearance between the head and substrate will increase under the shock. However, in a negative shock, the substrate is initially followed by the slider, but because of the different stiffness, the head substrate clearance will first decrease. Both positive and negative shocks are studied in our following simulations. Figure 7 exhibits the response of the full-coupled SPAS model to a 100-G positive shock. Head 3 is chosen to study the shock pulse effect. The results are shown in Figs. 7(a)-7(c). First, comparing with the PFH response without air bearing in Fig. 6(b), the head displacement in the SPAS model is quite small because the air bearing assists the PFH to sustain a steady near field condition. The suspension being more flexible than the substrate, the PFH will move far away from the substrate under a positive shock, and thus, the clearance between the PFH and substrate will increase as shown in Fig. 7(b). With the increase in the pulse width, the maximum displacement of all the components during the shock increases, resulting in a relatively large fluctuation in the post-shock. The fluctuation of PL d shows the same result. The large head slap owing to the separation and snap back may cause failure of the PFH. There are three frequency components in all the pulse cases, namely, 3.25 kHz, 6.54 kHz, and 12.4 kHz, observed in the power spectra in Fig. 7(c). The former two frequencies are the constrained modes of the suspension system, and the latter is the first substrate umbrella induced mode. Figures 7(d)-7(f) show the shock responses of the three ABS designs at the same pulse width. We can notice that the suspension components of the three designs exhibit similar responses. From the PL d response in Fig. 7(e), we can observed that head 2 has the highest stiffness to sustain a relative steady near-field gap, whereas head 1 has the least stiffness. This is because head 2 has a relatively high negative pressure that tends to improve the dynamic performance. The power spectra of the three ABS designs also show the three frequency components. Figure 8 presents the response of the full-coupled SPAS model to a 100-G negative shock. Figures 8(a)-8(c) show the head 3 responses at different pulse widths. We can expect that the suspension will come close to the substrate under a negative shock and bounce back owing to the high air bearing force. However, the head will contact the substrate, resulting in the failure in the large pulse width case, as shown in Figs. 8(a) and 8(b), because the inertia load of the substrate overcomes the air bearing force. The power spectra exhibits a constrained mode and substrate mode. The high stiffness performance of head 2 is also indicated in Fig. 8(e).

Load/unload analysis
The taking-off and landing of the PFH in RNFL is the key problem because of the peeling off the soft viscoelastic photoresist caused by the collision and friction. In a previous study [12], we have processed a transition zone for the taking-off to avoid the catastrophic collision between the plasmonic head and photoresist in the taking-off stage. The load/unload technology relies on the heads being lifted off the substrate onto a ramp. In comparison with the transition zone method, the load/unload method has a higher repeatability, a much higher shock resistance in the non-operational state, and the advantage of avoiding the problem of stiction that generates wear debris from the ABS and photoresist. However, the load/unload method must be studied to minimize the risk of the contacts in the operation process. In the simulation, the load/unload velocity is set as 50 mm/s and the substrate velocity is 5,400 rpm. Figure 9 shows the typical stages of the successful load/unload process of head 3. The force and displacements histories during the load process are displayed in Figs. 9 (a) and 9(b), respectively. They can be obviously divided into three stages. First, the dimple contact with the head causes a vibration of the head and tab, which can be illustrated by the ramp contact force and dimple contact force fluctuation. The displacement of each component also shows a fluctuation in the first stage. Second, the dimple contacts the head stably and the air bearing effect is gradually generated by the positive and negative pressures. Owing to the inertia of the main and local suspension, the flying height (FH) and air bearing force of the head initially experience a fluctuation and then tend to stabilize with the damping effect of the air bearing. In the final stage, the ramp contact force becomes nearly zero, indicating that the tab is separated from the ramp and head load successfully with a minimum near-field gap of 12 nm on the substrate. The results of the unload process are depicted in Figs. 9(c) and 9(d). The variations can be divided into four stages. In the first stage, the air bearing force decreases to zero and the FH steadily increases. The positive pressure resultant force of the bearing decreases and the suction force is almost a constant.
In the second stage, the dimple separates, reaching the maximum separation of the limiter distance. The bearing force changes from 0 to a negative value. In the third stage, the limiter is contacted and the contact force increases gradually. The air bearing force continues to decrease to a maximum negative value of −14 mN. Then the air bearing effect rapidly disappears, whereas the FH sharply increases. This occurs because the head is pulled by the suspension and limiter. In the last stage, the head strongly vibrates owing to the combined effects of the head inertia and dimple contact force. In addition, the PFH system experiences a vibration because of the suspension inertia. Figure 10 displays the load/unload performance comparison of the three different head designs. The force and minimum FH histories during the load process are shown in Figs. 10(a) and 10(b), respectively. There are two moments at which the head may contact the substrate. The first risk point is at the start of the second stage, where all the heads experience a vibration owing to the inertia of the suspension. Head 3 exhibits a better performance than the other two heads because of the rapid formation of the air bearing and its strong damping effect. The second point is in the process of the formation of the air bearing, during which the torque of the positive and negative pressures may cause the overturn of the heads. Therefore, a positive pitch torque is expected at the center location of the positive pressure, close to the air upstream, and the opposite is expected for the negative pressure. Head 2 undergoes a severe force  | https://mc03.manuscriptcentral.com/friction vibration during the end stage of the load process because of the high asperity force caused by the low FH. However, the FH has no significant fluctuation owing to the high stiffness. Figures 10(c) and 10(d) show the force and minimum FH histories during the unload process, respectively. Head 2 has a relatively high negative pressure that is advantageous for the stiffness, but makes the peeling off from the substrate difficult. It shows that head 3 has the lowest lift-off force of 5.8 mN compared with the values of 11.1 mN and 14.1 mN for the other two heads, implying that the unload performance of head 3 is the best among the three head designs.
Finally, based on the proposed the full-coupled SPAS model, a PFH analysis solver, SKLY.app (State Key Laboratory of Tribology + Fly) has been developed on the COMSOL Multiphysics Application Builder platform that can perform the analyses of the static, dynamic, and load/unload performance of the PFH or magnetic head in HDDs.

Conclusions
A novel full-coupled suspension-PFH-air-substrate (SPAS) model including the 6-DOF suspension-head dynamics, air bearing aerodynamics, and rotating substrate vibration for the flying performance analysis in RFNL, is proposed. The model exhibits a high accuracy and computing efficiency with a local pressure gradient-based adaptive mesh refinement algorithm. Three designs, heads 1, 2, and 3 of the ABS were investigated for the static, dynamic, and load/unload performance. Numerical simulations show that head 3 exhibits relatively good static and dynamic performances to achieve a ~15-nm near-field gap as well as has the best load/unload performance among the three ABS designs. Finally, a PFH analysis solver SKLY.app has been developed based on the proposed model. To achieve a finer lithography pattern, the design of ABS morphology with a near-field gap smaller than 10 nm and better load/unload performance must be further studied through the careful exploration of the dominant parameters.