Numerical study of pattern formation in compliant surfaces scraped by a rigid tip

The emergence of surface patterns on the surfaces of compliant materials subject to plowing wear is a complex problem which can be quantitatively characterized, e.g., on polymer surfaces scraped by an atomic force microscope (AFM) tip. Here we explore the applicability of a phenomenological model recently introduced to describe this phenomenon. Based on the competition between the viscoplastic indentation and the elastic shear stress caused by the tip, the model is able to reproduce the wavy features ( ripples ) observed when the tip is scanned along a series of parallel lines. For low values of the driving velocity v and the spacing b between scan lines, the existence of dotted areas formed by variously oriented pit alignments is observed. Moreover, coexistence of rippled with dotted domains is also observed at suitable parameter values. The formation process of the ripples is also described in detail. The amplitude, period, and orientation of these features are estimated numerically for different values of v and b parameters. We have also revisited the formation of the wavy patterns formed when a single line is scanned, and derived an equation which correctly describes their period and depth, and the static friction as well. This equation is not applicable when several lines are scanned one after the other and the ripples emerge as result of a cooperative process which involves the scanning of several lines. DOI:


I. INTRODUCTION
In 1992, Leung and Goh [1] published a pioneering work on the "orientational ordering of polymers by atomic force microscope tip-surface interaction." As stated in the abstract, they observed that, by action of an atomic force microscope (AFM) tip "nanometer-size structures are induced, resulting in a pattern that is periodic and is oriented perpendicular to the scan direction." Such wavy patterns, later called ripples, appear as the outcome of the tip motion, which reshapes the polymer surface locally. They have been observed and experimentally studied in depth by several authors afterwards [2][3][4][5][6][7][8][9][10][11]. However, these structures have been elusive to a theoretical understanding for years. This relates to the fact that the complex plowing process investigated in the aforementioned references involves concepts of physics, chemistry, and material science in a rather intermixed way which makes traditional approaches to the problem quite challenging.
In a recent investigation [12] and in a preliminary work limited to the scanning of single lines [9] we introduced a mesoscopic dynamic model which allowed us to reproduce the main features of the observed patterns. The model parameters were fitted to experiments performed over polystyrene films using a commercial AFM operated in contact mode and standard silicon probes. There it was showed that the surface patterns originated from the plowing process are ultimately due to the competition between two basic mechanisms: (i) viscoplastic indentation of the tip in the polymer surface caused by a constant normal force and (ii) lateral shear of the same surface caused by the tip elastically driven (at constant velocity) along a raster scan pattern.
As a result, the compliant surface evolves in time governed by two different rates: The indentation rate, which determines the change in the surface profile and, correspondingly, in the tip-surface interaction, and the driving velocity, which controls the time increase of the elastic energy stored in the cantilever supporting the tip. The onset of ripple patterns suggests that during indentation the tip tends to remain at a given position on the surface while the spring force increases till an elastic instability is reached and the tip slips into a new equilibrium position on the surface. This assumption is confirmed by the stick-slip behavior (sawtooth time variation of the spring force) observed experimentally [2,3,5]. Guided by a formal analogy with the Prandtl model for atomic-scale friction with variable potentials [13], the tip can be modeled as a point particle on a surface potential resembling the evolution of the polymer surface morphology. This is the basic idea of our approach.
After reviewing the used model, introduced in [9] and [12], we will rediscuss the simplest case of single scan lines, where we are able to offer an analytical interpretation beyond the numerical results presented so far. Then we will focus on the description and understanding of the ripple emergence process. At this point, it will become clear how the ripple emergence from parallel scan lines is a genuine three-dimensional (3D) process resulting from the correlation between consecutive scans. We will see that in case of strong correlation (short distance between scan lines and/or slow sliding velocity) the ripples are replaced by various arrangements of single pits. To end, we will characterize the ripples at different parameters in an experimentally accessible space.

II. MODEL
Within the conceptual frame of the Prandtl model, the tip dynamics of the AFM tip apex is described by the equation where is a damping coefficient, r 0 is the two-dimensional vector giving the tip position on the XY plane defined by the initially flat polymer surface, and U tot is the total potential energy, which includes two terms: The elastic potential of the cantilever spring, U elas , and the tip-surface interaction energy U sur . In Eq. (1) we assume that the dynamics of the tip is fully dissipative, thus the inertial term is completely neglected. The tip is thus pinned at the location corresponding to the minimum of the variable total potential during the stick phase. When the minimum turns into an unstable equilibrium point the tip slips towards a new minimum energy position guided by the gradient of U tot .
As in the Prandtl model [14][15][16][17] the cantilever spring energy is given by where R c is the support position andk is the (anisotropic) effective spring constant or stiffness of the system. For consistency with the experiments described in [12] it is indeed necessary to distinguish between two distinct values k x and k y (corresponding to longitudinal and transverse deformations with respect to the scan direction). The surface interaction energy landscape explored by the tip corresponds to a continuous field which evolves in time due to the penetrating action of the tip itself. This assumption is parametrized by an interaction potential which resembles the morphology of the evolving surface, given by Here z( r, t ) is the surface height at a given point, (ρ), with ρ = | r − r 0 |, is a geometric profile corresponding to the indentation print [ Fig. 1(a)], N (z 0 ) is the indentation rate, and is the conversion factor between z and the potential energy. We assume that it depends on the surface height at the tip position, z( r 0 ) = z 0 , as In this way, the indentation process is assumed to saturate as the penetration depth increases, Figure 1(b), with a crossover depth λ that can be estimated from the experiments, as shown in [9]. To keep the discussion as general as possible we assume that the indentation print profile is described by a combination of Gaussian functions: where σ determines the width of the indentation print, as seen in Fig. 1(a), and the parameter α ≈ 0.2 is chosen in a way such that ∞ 0 2πρ (ρ)dρ = 0, i.e., to guarantee that the volume below the surface is conserved, as expected in the case of plastic deformation.
Under these hypotheses, the shape and time evolution of the typical patterns observed on polymer surfaces indented by AFM are well reproduced [12]. In line with [9,12] the numeric simulations have been run with the values σ = 50 nm for the parameter defining the print width, N = 100 nm/s for the indentation rate, and λ = 1.8 nm for the crossover depth. The energy conversion factor involved in Eq. (3) is estimated as = 6 × 10 3 nN so that the energy indentation rate N is 6 × 10 −4 nJ/s and the crossover energy λ is about 10 −5 nJ. For the anisotropic spring constants we set k x = 2 nN/nm and k y = 10 nN/nm. The support moves at a velocity v of a few μm/s, with v = 10 μm/s the value taken in most of the simulations. Scanning occurs from left to right and from bottom to top in a series of parallel horizontal lines. The distance between consecutive scan lines, b, is varied between 5 and 100 nm (b = 20 nm is our reference value). Finally, the damping coefficient = 5 × 10 −4 kg/s is chosen to guarantee an overdamped dynamics and set a time constant suitable for numerical simulations reproducing the experimental observations.

III. RESULTS
Although we are interested in scanning 2D regions, it is important to start showing the results corresponding to a single scan line. Figure 2 shows the patterns obtained for five velocities between 0.05 and 10 μm/s. The patterns result from the competition between the indentation rate and the scan velocity. At low velocities the tip indents the sample deeply, with the indentation depth limited by the saturation effect. In this case the elastic force required to overcome the surface-tip interaction force is large and the tip remains stuck inside the indented well for a long time before slipping into the new equilibrium position where the indentation process starts again. The result is a series of periodically spaced pits, Fig. 2(a). As the velocity increases, the separation and depth of the pits decrease, Fig. 2(b). Contiguous pits are independent from each other until the separation gets close to the pit width 5σ , Fig. 2(c). If the velocity increases further, the pits overlap, Fig. 2(d). The stick-slip process is suppressed at higher velocities, and instead of a sequence of pits a single small-amplitude groove is obtained, Fig. 2(e). In this case the tip does not have time enough to indent the surface and moves smoothly along the surface. This result is reminiscent of the transition from stick-slip to continuous motion observed in atomic-scale friction measurements [18].
These observations can be quantitatively substantiated as follows. Under pure indentation, v = 0, the surface deformation is given by Note that (0) −1, so the argument in the log function is always positive. The pinning force is given by | − ∂U sur /∂ρ|  and has a maximum at ρ σ . This force acts as a time increasing static friction force. When the cantilever moves from left to right, the tip barely follows it and rather indents the surface as long as the pulling elastic force is smaller than this force. When the elastic force overcomes the maximum pinning force the tip slips, the elastic force goes to zero and a new indentation process starts. The effective indentation time, t i , can be computed from the expression Then, the distance between contiguous pits is approximately given by vt i + σ , its depth by z(0, t i ), and the maximum lateral force (in the longitudinal direction) by k x (vt i − σ ). This is fully confirmed by our numerical simulations, as shown in Fig. 3.
Let us now move on to consider the scanning of 2D regions. Figure 4(a) shows the typical rippled pattern obtained when a square area is scratched with a separation distance b = 20 nm (0.4σ ) between contiguous scan lines. At first glance one may think that the observed ripples are the natural evolution of the pits observed in the single scan line, with each pit slightly extended to the right (in the x direction) and up (in the y direction) as the scan process is repeated along the series of parallel lines. However this is not the case. As shown by our simulations, ripple emergence is a complex collective phenomenon which builds up from several scan lines. Figure 4 We can compare the pit size and the spacing b. The width of the indentation print is about 5σ (Fig. 1), which accounts for 12.5b at b = 20nm. It shows that every point in the square area is visited, and thus indented, by the tip several times. Figures 4(b) to 4(e) show pictures of the evolution of the surface topography after 1, 5, 10, and 20 scan lines. It is clearly seen that ripples emerge in a cooperative process only after a few scan lines. The initial groove first deepens after each new scan and then breaks starting from the left into several drops which move towards the upper-right direction to create the finally observed rippled pattern. A movie of the process is available as Supplemental Material [19]. Figure 4(f) shows the topography of the surface resulting after an arbitrary scan line (line 77) and, in green dots, the trajectory followed by the tip in this line and the previous four ones. Each new trajectory is shifted up by b = 20 nm and by a length x to the right with respect to the previous one. As a result the ripple pattern propagates in the northeast direction, with an inclination angle given by tan φ = b/ x. Figure 5 shows the time evolution of the two components of the spring force acting on the tip: F x = k x (X c − x 0 ) and F y = k y (Y c − y 0 ). Both forces vary periodically with time and show the expected stick-slip dynamics that suggested the use   Fig. 4(f).
of a modified Prandtl model. The dynamics can be better understood combining this figure with Fig. 4(f). During the stick phase the tip moves slowly to the right (slightly reducing the growth of F x ) and a bit up (which, due to the high value of k y , significantly decreases F y ), while trapped by the ripple groove. When the spring force becomes high enough to unpin the tip, the latter slips into the next groove where it gets trapped again. In the chosen approximation the tip always follows the minimum of the time varying total potential energy. As a result, both components of the spring force oscillate in antiphase with a similar amplitude of about 250 nN.
We see that the ripple inclination is larger for larger b values. This was already observed in our previous simulations and experimental results in [12]. There it was reported the observation of ripples for b values between 0.2σ and 1.2σ and the increasing of the ripple tilt angle with b in full agreement with numerical results. We include here additional figures to show that arrays of pits rather than ripples are seen at small b values, Fig. 6(a). The pits are mainly aligned in the northwest direction, except in the lower part of the image, where they are northeast oriented. The extension of the lower part increases for higher values of b, till the pits merge into the ripples and the northeast alignments are suppressed for b 0.2σ , Fig. 6(b). If we increase b an uniform ripple pattern is observed; see Fig. 4(a) where b = 0.4σ . For larger b values as shown in Fig. 6(c), vertical ripples appear also as a border effect. At high enough b 2.5σ consecutive scan lines do not significantly interfere and we just get a series of horizontal grooves, Fig. 6(d).
In Figs. 6(a1), 6(a2), and 6(b1) we show enlarged areas of Figs. 6(a) and 6(b) to better appreciate the character of the boundaries between different phases. In Fig. 6(a1) two domain walls, soliton type [20], separating similar regions slightly displaced one from each other are observed. In Fig. 6(a2) are visible the transitions between three different patterns: Ripples in the lower part and two differently aligned sequences of pits above the ripples. In Fig. 6(b) ripples dominate the bottom part of the figure but chains of pits appear above. Figure 6(b1) shows the transition between the two configurations.
Next we present results obtained for different values of the scan velocity v, Fig. 7. An array of pits is observed at low v values, where indentation is deeper. As v increases, ripples If the ripple pattern is cut either along the x or y direction the surface profile is almost sinusoidal. Not only the ripple orientation is well defined, but also their amplitude and period, and it makes sense to study how these quantities change when some of the experimentally controlled parameters are varied. This is done in Fig. 8, where the values of these quantities are plotted as functions of v and b. Increasing either of these parameters results in a decrease of the ripple amplitude and an increase of their tilt angle. The ripple period along the x direction, λ x is much less influenced by v and b. Indeed λ x ≈ 200 nm, i.e., close to the value of 5σ = 250 nm for the amplitude of the tip print. A very different behavior is observed when single lines are scanned, where the distance between consecutive pits decreases with increasing v (Fig. 2). Thus, it is further confirmed that the ripple period results from the combined effect of tip geometry and the cooperative effects accompanying the tip dynamics. In Figs. 8(c) and 8(d) we also show the real value of the period of the rippled pattern, perpendicular to the ripple crests, λ ⊥ = λ x sin φ with φ the orientation angle of the ripple.

IV. DISCUSSION AND CONCLUSIONS
We have shown the results of a rather simple model aimed to reproduce a complex surface structuring problem. The model allows for understanding many of the experimental observations made to date and can be used as a guide for predicting the appearance of different patterns under suitable working conditions. It is important to mention that the observed patterns were numerically found to be very robust against thermal noise and surface disorder, as was already reported in [12]. Here, we have numerically explored the existence of ripples when playing with two experimentally controllable model parameters: b and v. Dynamics and patterning result from the competition between two rates: The scan velocity v (which favors the slip) and indentation rate N (which favors the stick). In our case v/N ∼ 100. These parameters essentially define the two energy (force) relevant rates. On the one hand the driving velocity causes an increase of the elastic energy, and thus elastic force, in the system. On the other hand, the indentation rate gives an increasing of the pinning energy and thus of the pinning force. With our parameters the increase rates for the elastic energy (∼k x σ v) when compared to the energy indentation rate ( N = 6 × 10 −4 nJ/s) gives a ratio of about 5/3. This ratio can be used as a guide for seeking ripples at arbitrary parameter values.
We also assumed dissipative dynamics, which is guaranteed if (ω/γ ) 2 1, where ω = √ k x /m, γ = /m, and m corresponds to the effective mass of the tip relevant in the problem. For instance, in our case (ω/γ ) 2 < 0.01 for m < 10 −9 kg, which definitely holds for standard AFM probes.
Our general model could be easily refined to take into account specific experimental details. For instance, we have used a simple Gaussian profile for the indentation print, but this profile could be modified by introducing, e.g., pyramidal (Berkovich type) tips. In this context, the importance of tip face orientation has been suggested by recent experiments by Yan et al [11,21].
In Eq. (5) the indentation rate N depends on z 0 in a simple empirical way in order to get a logarithmic time increase of the indentation depth, similar to that observed in previous experiments (see Fig. 7 in [9]). For simplicity σ was chosen as constant. Thus the indentation print has a (constant) width of about 5σ in the whole indentation process. Different σ (z) and N (z) dependencies are required to describe the tip geometry more realistically, but the correct simulation of a σ (z) dependence is a great numerical challenge since it would require a very fine spatial discretization for U (x, y). In addition, when the contact is established σ → 0 and the pressure exerted by the tip diverges. Some cutoff should be introduced to numerically overcome both problems.
Another important parameter is the value of the driving spring constant. We have already shown in [12] that, due the shape and orientation of the cantilever, the elastic response of the system is anisotropic: Ripples do not appear if the sample is scanned in the direction of the cantilever axis. The lateral stiffness of the cantilever associated with its buckling (in the y direction) is much larger than the lateral stiffness associated with its torsion (along x) [22]. Our numerical simulations show that ripples appear at different k x values only if k y 3k x (in the chosen parameter space). Further increasing the value of k y does not change the results significantly, since the effect of this parameter is only limiting the mobility of the tip in the y direction in this case. On the other hand, the lower the value of k x is, the deeper is the tip penetration. The elastic force roughly behaves as k x vt and thus it grows more slowly at small k x values. It favous the production of patterns with pits. In the opposite case, high k x values weaken the stick-slip. In this case the tip follows closely the cantilever and ripples are not created.
To sum up, we have studied a quite simple and general model that provides all necessary ingredients to simulate the observed formation of surface structures on polystyrene surfaces scraped by a sharp tip. Polymers are a class of materials easily accessible to experimental validations on the nanoscale, but the application of the model to other compliant materials would be also interesting to test. Our work can be also considered as a possible method to design nanostructured materials using scanning probes, a topic of indubitable current interest [23]. In this regard, the surprising qualitative similarities observed with the problem of nanoscale pattern formation at surfaces by ion-beam irradiation is also worth mentioning (for a recent review on the topic, see [24]). 16 × 10 6 points). The tip position, r 0 (t ), can take any value in the XY plane.
Given the tip position and surface profile at time t the algorithm computes the new tip position at time t + dt, with dt the RK time step (fixed typically to 50 μs in our case, which at v = 10 μm/s gives an advance of 0.5 nm of the cantilever each time step). The spatial gradient involved in the equations is evaluated in the tip position from the discrete surface field z g with a centered next nearest neighbor finite difference approximation based on [25]. We observed that the use of a simpler linear nearest neighbor interpolation produces numerical instabilities. Once the new tip position is obtained, this result is used to update the surface z g by solving Eq. (4) with a forward Euler method. The spatial range of this update is limited to the nodes of the grid at distance less than 5σ from the tip position (about 10 6 nodes). At this distance, | (5σ )/ (0)| ≈ 0.002. We have checked that our numerical results do not change if we use a finer grid with 0.5 nm spaced nodes, smaller time steps, or a longer cutoff for updating the surface profile. In both cases the obtained surface field z( r, t ) and tip dynamics r 0 (t ) are basically the same. We have also numerically checked the robustness of the solutions against spatial disorder (starting from a reasonable rough surface [26]) and thermal fluctuations (in this case we solved the equations using a stochastic RK algorithm).