Roughness of surface decorated with randomly distributed pillars

We have presented a quantitative analysis of roughness of planar surfaces decorated with randomly distributed, cylindrical pillars, disks, or cavities. We have described the roughness in terms of the surface power spectral density (PSD). First, we have derived a general equation for the PSD of such surfaces. Then, we have found the PSD for the special case of statistically isotropic, circular areas. We have demonstrated that the PSD provides quantitative information on the radius of the circular area, dimensions of the pillar, and surface coverage. We have also discussed the numerical method of extracting the parameters from experimental PSD data obtained from discrete Fourier transform of surface scanning measurements.

Finally, we have discussed our results in the context of potential application of the PSD for the quantitative determination of the pillar layer parameters.

Analytical Considerations
Fourier transform of pillar decorated surface. Let us consider a planar area A in the Cartesian coordinate system shown in Fig. 1, enclosed by the curve r s (α), where α is the polar angle of the vector r s . The area is decorated with N identical, circular, randomly distributed pillars of the height h and is otherwise smooth. We can uniquely describe the system geometry by specifying the curve r s (α) and set of position vectors r j = (x j , y j ), j = 1, …, N, where x j and y j are the coordinates of the j-th pillar's base center. Unless stated otherwise, for the sake of notation simplicity, we have normalized all dimensional parameters and variables by the pillar radius a to a power matching their physical dimensions, to make them dimensionless. In what follows, to neglect boundary effects, we have assumed that linear dimensions of the area A are many times larger than the pillar radius, i.e., that |r s (α)| ≫ 1 for any value of α.
Once we know the position vectors of all pillars, we can quantitatively describe the height of the pillar decorated surface. For that, we have used the function of height: Equation (2) describes the function of height of the surface decorated with a single pillar located at r j . Mathematically, the single pillar system can be generated by translation of the pillar, originally located at the origin, by the vector r j . In Fig. 2a, we have presented the pillar at its original position. The function of height for this system has the canonical form Please note that this function, because of axial symmetry of the system, depends on the modulus r of the vector r only.
The Fourier transform of z(r, r 1 where Z j (q, r j ) is the Fourier transform of z j (r, r j ) over A and q is the wave vector of the modulus q = 2π/λ. Here, λ denotes the dimensionless wavelength. We can exploit the Fourier transform translation theorem to rewrite the single pillar transform as where Z 0 (q) is the Fourier transform of z 0 (r) over A and i is the imaginary unit. Please note that the function Z 0 (q) depends on the wavenumber q only. It is easy to show that the transform equals A 0 2 0 1 Figure 1. Scheme of planar area A in Cartesian coordinates, decorated with randomly distributed pillars. The area is enclosed by the curve r s (α) 18 . where J 1 (q) is the first order Bessel function of the first kind. In Eq. (6), as well as in the whole paper, we have assumed that q > 0. Thus, in the very special case of pillar at the origin, its Fourier transform is a purely real function. In Fig. 2b, we have presented the function Z 0 (q) in the reciprocal space. In general, however, the transform Z j is a complex, highly variable function of the vector q. In Fig. 2c and d, we have presented the real and imaginary parts of the function, respectively, for the vector r j = (1, 1). Substituting Eqs (5) and (6) into Eq. (4), we obtain the Fourier transform of pillar decorated surface in the form j N 1 N j 1 1 PSD of pillar decorated surface. The PSD of a finite surface area equals 1 where the angle brackets denote ensemble averaging. Substituting Eq. (7) into Eq. (8) and using the trigonometric form of complex numbers, we get Please note that the available pillar positions are determined by the shape of the area A. Therefore, C(q, N) is dependent both on the area's size and, implicitly, on its shape.
A simple transformation of Eq. (9) yields   1 2 which is independent of the pillar height h and, as we will see later on, is relatively weakly dependent on the pillar number N.
As we can see, the function T(q, N) represents an ensemble averaged sum over all ordered pillar pairs. We can distinguish two subsets of the pairs. The first contains all the uncorrelated pillar pairs at a distance r lm larger than some characteristic distance r c . This distance depends on the distribution of pillar positions on the surface. In systems of the structure consistent with the model of random sequential adsorption (RSA) 7 , e.g., the distance changes with surface coverage and, in practice, is always less than five 10  Let us first consider correlated pillar pairs. Again, we can distinguish two subsets of the pairs. The first contains all the ordered pairs with the first pillar located in the inner part A i of the area A, separated from the boundary r s (α) by the distance r c (see Fig. 3a). The contribution from these pillar pairs to the first term on the RHS of Eq. (13) equals where r and r′ are the position vectors of the first and second pillars located, respectively, in the area A i and circular area A c = π r c 2 centered at the point r, n(r) is the ensemble averaged pillar areal number density at the point r, and n(r′, r) is the ensemble averaged pillar areal number density at the point r′ given a pillar at r.
Because of statistical homogeneity and local, at the scale of r c , statistical isotropy of the system, the pillar number densities are equal n(r) = θ/π and n(r′, r) = g(ρ, θ) θ/π, where θ = π N/A is the pillar surface coverage, ρ is the length of the vector ρ = r′ -r, and g(ρ, θ) is the radial distribution function. Please note that this function equals zero for ρ ≤ 2, therefore we can integrate over the entire area A c without violating the condition of pillar non-overlapping. Thus, for any r, the pillar number density n(r′, r) is described by the same function of ρ and θ only. In this sense, the function is independent of r. The same is true for the argument of the cosine function appearing in Eq. (14). Therefore, we can separate the integration over A i and A c . Considering this and using polar coordinates, we can rewrite Eq. (14) as where β and φ are the phase angle of the vector q and polar angle of the vector ρ, respectively. We can express the inner integral in a closed-form 11 , Eq. 9.1.18; 12 , Eq. 10.9.1 to get is the zeroth order Bessel function of the first kind.
Next, let us consider the second subset of correlated, ordered pillar pairs with the first pillar located in the outer part A o of the area A, close to the curve r s (α) (see Fig. 3b). Their contribution to the first term on the RHS of Eq. (13) equals where A e (r) denotes the external part of the disk A c , extending outside the area A. Using the same reasoning as for the area A i , we can express Eq. (17) as Here, ρ 0 (r) denotes the minimum value of |ρ|, while φ 0 (ρ, r) and φ 1 (ρ, r) are the minimum and maximum values of the polar angle φ (see Fig. 3b): where γ is the polar angle of the vector r and, from the law of cosines, The angles α i , appearing in Eq. (21), denote the minimum and maximum polar angles of the vector r s for the given area A e and value of ρ. Summing up the contributions from both subsets of correlated pillar pairs, we get the first term on the RHS of Eq. (13) to be The second term on the RHS of Eq. (13) represents the contribution from uncorrelated pillar pairs. We can calculate it as a difference between the integral over all possible positions of the ordered pillar pairs and integral over the positions of correlated pairs at close distances: n n q r r q r r r r r r q r r r r r r q r r r r or, in terms of surface coverage, Please note that Eq. (23) is correct whether or not we neglect the correlations between close positions of pillars, as long as we do so in all integrals. Therefore, for the sake of simplicity, we have neglected the correlations and assumed that everywhere in Eq. (23) g(ρ, θ) = 1, which is the correct value for uncorrelated pillar pairs. That yields  Summing up the contributions from correlated and uncorrelated pillar pairs, we get is the total correlation function 13 , p.72 and We can split the middle integral on the RHS of Eq. (25) into two pieces, the part over 0 ≤ ρ < 2 and the part over 2 ≤ ρ ≤ r c . In the former, h(ρ, θ) = −1 and therefore 11 , Eq. 9.1.30; 12 , Eq. 10.6.6 Let us comment on the physical meaning of the five terms of the sum appearing on the RHS of Eq. (31). The first, constant term represents the contribution from single pillars. The second term represents the contribution from all ordered pillar pairs with neglected overlapping and correlations between pillar positions. Thus, in this integral we have taken into account also the pairs of overlapping pillars and we have neglected deviations of the pillar number density at the distance ρ > 2. The third component is a correction term for pillar overlapping. The fourth component is a correction term for correlations between pillar positions. As we can see in Eq. (28), we have calculated the two correction terms by the integration over the radius r c of the disk A c . As a consequence, we have integrated over the second pillar position extending outside the total integration area A, which leads to an overestimation of the two terms. To compensate the overestimation, we have subtracted the fifth term, which is the integral over the outer area A o and external part of the disk, A e . Thus, the fifth component is a correction term for the third and fourth components.
If, for any α, r s (α) ≫ r c , then A o ≪ A and the last term of the sum on the RHS of Eq. (31) becomes negligibly small. Then, in the case of large surface area, Eq. (31) simplifies to If, on the other hand, θ ≪1, the total correlation function at ρ > 2 tends to zero. Then, in the case of low surface coverage, the integral I c (q, θ) → 0 and we get Please note that the equations derived in this section are valid also for randomly distributed cylindrical disks and cavities.

PSD of pillar decorated circular area. Equations (31)-(35) are valid for surfaces of various shape. In
what follows, we have restricted ourselves to a disk of a radius r s , centered at the origin. Then, we can derive a closed-form expression for the second term of the sum on the RHS of Eq. (31), representing the contribution from all pillar pairs as an explicit function of r s :  where the integrals I c (q, θ) and I e (q, θ, r s ) are given by Eqs (29) and (44), respectively. Please note that this function is independent of wave-vector direction, as we have expected for the system of circular symmetry. Thus, in general, to get the reduced PSD, we need to calculate the integrals I c (q, θ) and I e (q, θ, r s ). As we can see from Eqs (29)    respectively.

Numerical Verification
In principle, we can use Eq. (45) to determine important parameters of pillar layers. Specifically, we can find the pillar radius and height, disk radius, and surface coverage. For that, we have first to determine the function of height of the pillar layer of interest by, e.g., scanning its surface. Then, we need to numerically calculate the discrete Fourier transform of the function of height. There are a few open-source packages available for this purpose. Perhaps the most common is the package FFTW 14,15 . In our earlier paper, we described the application of this package in calculating the PSD of spherical particle monolayer 16 . Input data to the programs should usually be in the form of a rectangular array of z-coordinates of the scanned rectangular area A s . Therefore, to ensure a circular shape of the pillar decorated surface, we need to assign the height z = 0 to all points located outside the analyzed disk area A. Finally, we have to fit the dimensional version of Eq. (45), i.e.,  to the numerical results of discrete Fourier transform. In Eq. (50), Q = q/a, H = ha, R s = r s a, and S = Ca 4 are the dimensional wave number, pillar height, disk radius, and PSD, respectively. Please note that, to take into account the non-circular shape of the scanned area, we have to replace the factor θ/π, appearing in the dimensional equation for the PSD, with the factor θR s 2 /A s . In general, the application of Eq. (50) requires explicit analytical expressions for the integrals I c (q, θ) and I e (q, θ, r s ) that are 2D and 3D, highly variable functions of wavenumber. To find the expressions for a specific radial distribution function, we have to numerically calculate the integrals for a number of θ and r s values, over a range of wavenumber. Then, we can use the least-square method to fit analytical expressions to the numerical results. This procedure is out of the scope of this paper. We are going to present numerical results of such calculations for pillar layers consistent with the RSA distribution in a future report. Our preliminary results suggest, however, that the simplified form of Eq. (50) with the integrals I c (q, θ) and I e (q, θ, r s ) neglected, is sufficient to determine the parameters a and R s with a relative standard error less than 1%, at any surface coverage and for any size of analyzed surface as long as r s ≥ r c . To achieve this accuracy, the pillar number N must be large enough and therefore, at low surface coverage, ensemble averaging can be necessary. Figure 5 presents a comparison of PSD calculated numerically with the library FFTW (points) and fits of Eq. (51) to the numerical data (lines). We have derived the numerical results from a virtual system of 160 pillars of the radius a = 0.5 μm and height H = 1 μm, randomly distributed over a square area A s = 2500 μm 2 , created according to the RSA model. Thus, the pillar surface coverage in this system is low and equals θ = 0.05. We have determined the height of the system at 512 2 grid points equidistantly distributed over the area A s . To calculate the PSD for a circular area, we have set z = 0 at all grid points outside the disk of the radius R s = 25 μm, inscribed in the area A s . At this value of R s , corresponding to r s = 50, neglecting the integral I e (q, θ, r s ) results in the PSD relative error less than 1%, as we have found in our preliminary calculations. Then, using FFTW library, we have calculated discrete Fourier transform of the data. Finally, we have calculated powers of all wave-vectors and summed them over narrow wavenumber intervals to calculate the discrete values of PSD. To back determine the system parameters, we have fitted Eq. (51) to the numerical data over low-and high-wavenumber ranges where the PSD shows numerous, sharp minima determined by the Bessel functions appearing in Eqs. (50) and (51). From the general properties of the function J 1 (x) we can deduce that in the limit of small wavenumber, where Qa ≪ 1, the PSD is determined by the term (J 1 (QR s )/Qa) 2 with positions of minima dependent on R s . Thus, fitting the data in the range of low-wavenumber, we can find the disk radius. In the limit of large wavenumer, on the other hand, the PSD is determined by the factor (J 1 (Qa)/Q) 2 with zeros dependent on a. Therefore, fitting the data in the Figure 5. PSD of a circular area decorated with randomly distributed pillars. Points denote numerical results derived, through discrete Fourier transform of 512 2 height measurements, from a virtual pillar system with the parameters a = 0.5 μm, H = 1 μm, A s = 2500 μm 2 , R s = 25 μm, and θ = 0.05. Lines represent leastsquare fits to the numerical data in the low-, mid-, and high-wavenumber ranges 20 . In the lower and upper ranges, for Q < 7000 cm −1 and Q > 10 5 cm −1 , respectively, we have fitted Eq. (51). In the mid-range, for 3.6 × 10 4 cm −1 < Q < 6.6 × 10 4 cm −1 , we have used Eq. (52). high-wavenumber range, we can find the pillar radius. In our system, we have obtained the following radii and standard errors: a = 0.5007 ± 0.0007 μm and R s = 25.93 ± 0.21 μm. The fitting results are quite consistent with the original values used in our RSA simulation. The PSD simplicity in the limiting wavenumber ranges has also a negative consequence: the parameters H and θ are strongly correlated and their standard errors are usually very large. Moreover, sampling of the PSD in the low-wavenumber range is poor. In the high-wavenumber range, on the other hand, the calculated PSD of the rapidly changing function of height can be aliased because the discrete Fourier transform of such function is not bandwidth limited 17 , p. 494. Therefore, to determine the pillar height and surface coverage, we have to carry out the fitting procedure over a mid-wavenumber range. In our system, we have chosen the wavenumber interval 3.6 × 10 4 cm −1 < Q < 6.6 × 10 4 cm −1 . We have fixed the radii at a = 0.5007 μm and R s = 25.93 μm, according to the results of fitting in the lower and upper ranges. Also, to eliminate the PSD suppression with the increase in Q, we have fitted the reduced form of Eq. (51), i.e.,

Conclusion
Our results suggest that the PSD can be a valuable source of quantitative information on parameters of planar surfaces decorated with randomly distributed, cylindrical pillars, disks, or cavities. In the case of statistically isotropic, circular surface area, the PSD provides us with the radius of the area, dimensions of the pillar, and surface coverage. The analysis presented here provides means to design and produce surfaces of controlled roughness. To better exploit the results, we need to numerically calculate the integrals appearing in the equation for the PSD and to approximate them with analytical, 2D and 3D functions of the wavenumber, pillar surface coverage, and disk radius. That will be a subject of our future investigation.

Data Availability
All data generated and analysed during the current study are available from the corresponding author on reasonable request.