Combined Modeling of the Optical Anisotropy of Porous Thin Films

: In this article, a combined approach for studying the optical anisotropy of porous thin films obtained by the glancing angle deposition is presented. This approach combines modeling on the atomistic and continuum levels. First, thin films clusters are obtained using the full ‐ atomistic molecular dynamics simulation of the deposition process. Then, these clusters are represented as a medium with anisotropic pores, the shapes parameters of which are determined using the Monte Carlo based method. The difference in the main components of the refractive index is calculated in the framework of the anisotropic Bruggeman effective medium theory. The presented approach is tested and validated by comparing the analytical and simulation results for the model problems, and then is applied to silicon dioxide thin films. It is found that the maximum difference between the main components of the refractive index is 0.035 in a film deposited at an angle of 80°. The simulation results agree with the experimental data reported in the literature.


Introduction
The deposition in a vacuum is the most widely used technique for thin film production [1]. The structural and optical properties of the deposited films depend on the deposition condition, in particular, on the angle between the direction of atoms flow and normal to the substrate surface [2]. Deposition at large angles, when atoms move almost parallel to the substrate, leads to the formation of glancing angle deposited (GLAD) films consisting of the separate nanostructures with different shapes and dimensions [3]. Due to the low reflectance and anisotropy of the refractive index, GLAD films are widely used in optical instruments [4][5][6][7][8].
To describe the anisotropic properties of the deposited films, the anisotropic Bruggeman effective medium approach is usually used [6,[9][10][11][12][13][14]. The Bruggeman formalism was initially developed for the calculation of the effective dielectric constant of a medium with randomly distributed spherical inclusions [15]. The values of the dielectric constant of medium and inclusions are different. In Reference [16], the approach was generalized for the case of the ellipsoidal inclusions. The anisotropy of a medium with inclusions was taken into account through depolarization factors depending on the ellipsoids shape parameters.
The voids between separated nanostructures in GLAD films can be considered as inclusions of different shapes and dimensions in the films' matter. Within the framework of the Bruggeman's formalism, these inclusions are represented by a set of identical ellipsoids oriented in the same direction. Since the experimental determination of the shape parameters of ellipsoids is still a challenge, these parameters are used to fit experimental results within effective medium models [6,[10][11][12][13][14]. On the other hand, shape parameters can be calculated through the Cartezian coordinates of atoms in the clusters obtained by the simulation of the thin film deposition process. Currently, the different methods of the atomistic level, including the ballistic approach [17,18], molecular dynamics (MD) [19,20], and kinetic Monte Carlo (kMC) [14] method were applied to study the deposition process, structural and mechanical properties of GLAD films. In Reference [14], shape parameters were defined for the titanium dioxide GLAD films using the structure tensor, calculated by the integral of density gradient [21]. The main components of this tensor correspond to the averaged inverse shape parameters. The planar birefringence, important in GLAD films applications for polarization optic [22,23], was calculated using the shape parameters and was compared with the experimental data from Reference [14].
In the present work, a combined approach for modeling of the optical anisotropy of porous thin films is presented. First, the MD-based procedure developed earlier in References [24][25][26] is used to simulate the high-energy deposition of silicon dioxide thin films. Further, the averaged pore shape parameters determining the anisotropy properties of porous thin films are calculated using the Monte Carlo based method. The method is quite general and is characterized by high computational efficiency. It can be applied to any porous atomistic clusters, obtained using MD, kMC or ballistic approach and consisting of up to 10 7 -10 8 atoms. Finally, the values of the difference between the main components of the refractive index ∆n are calculated in the framework of the anisotropic Bruggeman effective medium approach. The approach is validated by comparing the analytical and simulation results for model problems. The obtained results agree with experimental data for silicon dioxide thin films.

Simulation Method
This section presents a combined approach to calculating ∆n. We start from the description of the MD procedure that is used in the simulation of thin film growth. Then a brief overview of the anisotropic Bruggeman effective medium theory as applied to the porous atomistic clusters is given. In the final subsection, an original Monte Carlo based approach to the calculation of the averaged shape parameters is described.

Molecular Dynamics Simulation of Thin Film Deposition Process
The simulation of thin film growth is the most time-consuming step in the ∆n calculation. Nevertheless, the use of parallel computations allows one to increase the dimensions of clusters in the full-atomistic MD simulation up to several tens of nanometers. An increase in cluster dimensions is especially significant in the case of GLAD films consisting of separate nanoscale structures.
The geometry of the deposition process is shown in Figure 1. The simulation area includes parts of the substrate, film and gas phase above growing film. An increase in the deposition angle α leads to an increase in the porosity of the film and anisotropy of the structure. In particular, deposition at α = 80° leads to the formation of an anisotropic structure with separated slanted columns ( Figure 1). For this reason, the dependence of ∆n on α is investigated in a wide range of α from 0° to 80°. In this work, we use clusters that were simulated in our previous works [25,26], except for the cluster deposited at an angle α = 50°.
All silicon dioxide film clusters were obtained using the MD-based step by step method, as described in Reference [27]. The energy of interatomic interactions was calculated in the frame of the DESIL force field [27]: where qi(j) charge of i(j)-th atom, qO = −0.5qSi = −0.65e, Aij and Bij, parameters of Lennard-Jones potential for the van der Waals interaction, rij − interatomic distance, ASiO = 4.6 × 10 −8 kJ•nm 12 /mol, ASiSi = AOO = 1.5×10 −6 kJ•nm 12 /mol, BSiO = 4.2 × 10 −3 kJ•nm 6 /mol, BSiSi = BOO = 5 × 10 −5 kJ•nm 6 /mol. MD simulation is performed using the GROMACS package [28]. The leap-frog algorithm is used for the integration of particles motion equations [29]. Horizontal dimensions of the substrate are 20 nm and 60 nm, the vertical thickness of the substrate is taken equal to 6 nm (Figure 1). At each injection step, silicon and oxygen atoms with the stoichiometric proportion of 1:2 are inserted randomly at the top of the simulation box. Then injected atoms move toward the substrate and form new sub-layers on the surface of the growing film. The maximum value of the deposited film thicknesses is about 45 nm, which corresponds to about 4000 injection steps. The initial values of the silicon and oxygen atoms kinetic energies are E(Si) = 10 eV and E(O) = 0.1 eV, which corresponds to the high-energy deposition processes like ion beam sputtering (IBS) [3]. At each injection step, the NVT (constant number of particles, volume and temperature) ensemble is used. The vertical dimension of the simulation box is increased by 0.01 nm after each injection step in order to compensate for the growth of film thickness. The horizontal dimensions of the box remain unchanged. The Berendsen thermostat [30] is applied to keep the simulation cluster temperature, T = 300 K, constant. The duration of one injection step is 6 ps, and the time step of MD modeling is 0.5 fs. The simulation was carried out using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University [31].

Bruggeman Effective Medium Theory for the Porous Structures
In the framework of the anisotropic Bruggeman effective medium approach, the pores in the structure are approximated by a set of ellipsoids [9]. These ellipsoids are randomly distributed over the film volume and have the same shape parameters ax, ay, az and the same orientations ( Figure 2). The original method for the calculation of shape parameters is described in Section 2.3. After determining the values of ax, ay and az, the depolarization factors Lx(y,z) are calculated as follows [9,16]: where x ' is the integration parameter. These factors depend only on the ratios of ax, ay and az. Indeed, using the substitution x ' = ax 2 t, we obtain for Lx: where t is the dimensionless variable. Similar equations can be obtained for others depolarization factors.
The main components of the dielectric constant tensor εx(y,z) of the medium consisting of k materials are calculated as follows [9,16]: , where fi = Vi/V and εi are the volume fraction and dielectric constant of the i-th material. The components of the refractive index nx(y,z) are calculated using the Maxwell relation between refractive index and dielectric constant n 2 = ε.
For the case of porous films Equation (4) is simplified as follows: where εdf = 2.22 [32] is the dielectric constant of dense silicon dioxide film, f1 is the free volume fraction of film. Value of f1 is calculated using the film density dependence on the deposition angle (α): In Equation (6), it is assumed that a normally deposited film (α = 0) has no voids and the density of the dense fraction in films, deposited under different α, is equal to (0).

Calculation of the Averaged Shape Parameters of Pores in Deposited Film
A sample of the growing film structure is shown in the right part of Figure 3. The elongated pores are formed between columns and are highlighted in dark blue. It is seen that the pores are slanted approximately in the same direction.
In the framework of the anisotropic Bruggeman effective medium approach, it is necessary to calculate the averaged pores dimensions ax, ay and az along the axes of the coordinate system. Furthermore, these parameters are considered as ellipsoids shape parameters in calculating the depolarization factors, see Equations (2) and (3). The coordinate system can be defined by the plane of incidence (x,z) and the plane (x,y) slanted to the substrate plane at an angle β [9]. The Y axis is directed parallel to the substrate plane and perpendicular to the incidence plane, due to the symmetry of the deposition process geometry (Figure 3). The Z axis should be oriented along the direction in which the pores are elongated. This direction is given by the angle βm at which the ratio az/ax reaches its maximum.  The procedure for calculating the ax(y,z) and nx(y,z) values consists of the following steps:  The initial value of β angle is chosen.  Set of N points with random coordinates (xi, yi, zi) inside the cluster is generated. For each of the Nin points located inside the pores, the pore dimensions axi, ayi, and azi are calculated, as shown in The averaged pores dimensions ax(y,z), are calculated as follows: When calculating Equation (7), there is no need to determine the volume of every individual pore, which ensures high computational efficiency of the algorithm.  Step 2 and 3 are repeated with growing Nin values. The procedure is finished when ax(y,z)(Nin) are varying insignificantly with a further increase in Nin. The ratio az/ax is calculated.  The values of Lx(y,z) and nx(y,z) are calculated using Equations (3) and (4). The difference of the main components of the refractive index ∆n is calculated.  The value of β is changed, and steps 2-5 are repeated.  The ax(y,z)(Nin) values calculated using β, corresponding to the maximum of az/ax ratio, are the shape parameters. Let us prove that Equation (7) gives the ax(y,z) values averaged over the pores volume. Indeed, the volume element dV(az), including all of the points giving the same value of az, is defined in Figure 3: In the frame of the Monte Carlo method [33] and in the Nin →∞ limit, the next condition is satisfied: (9) where N(az) is the number of points, giving az value. Thus, the summation in (7) can be replaced by integration: The values of ay(z) are calculated in a similar way. For the case of ellipsoidal pores of different dimensions, oriented in the same direction, Equation (10) gives values of ax(y,z) which are one and a half times higher than the averaged values of the ellipsoids semi-axes (Appendix A). However, this difference does not matter, since only ratios of ax, ay and az are required in Equation (3). The pore boundary search algorithm takes into account that every atom in the clusters is considered as a hard sphere with a certain radius R. The values of the atomic radii are discussed in the text above Table 1, in Section 3. Thus, the pore surface consists of fragments of spherical surfaces centered on the atoms. To calculate the ax(y,z)I values that are required in Equation (7), the next algorithm is used:  The coordinate system is centered in the probe point, which is randomly chosen.  The distance from each atom to the X axis is equal to .  If rA is less than atomic radius R, the distance from the probe point to the point of intersection of the X axis with the sphere of radius R, centered on the atom, is calculated as follows:  The minimum value of x2A is equal to x2i.  The x1i value is calculated similarly. Then the axi value is calculated as axi = x2i -x1i.  The y1(2)i and z1 (2)

Results and Discussion
To validate the method for calculating shape parameters, two model geometries of the pores are considered: A spherical pore and a set of ellipsoidal pores of various dimensions (right side of Figure  5). Pores are created in the atomistic cluster of silicon dioxide film obtained by high-energy deposition. Horizontal dimensions of the cluster are 20 nm and 60 nm, the vertical dimension is 30 nm, and the directions of the coordinate axes can be seen in Figure 3. Atoms whose centers are inside the pores are removed.
The diameter of the spherical pore D is taken equal to 6 nm. For a continuous medium, the integral on the right-hand side of Equation (10) can be calculated analytically. This gives ax = ay = az = 3D/4 = 4.5 nm. The results of numerical simulation are shown in the plots in the upper part of Figure  5. Values of ax(y,z) converge to approximately 4 nm with an increase in the number of probe points. The small difference between the results of analytical and numerical solutions is due to the difference between the pore surface and the spherical surface. Indeed, the pore surface is formed by the fragments of spheres centered on the atoms closest to the center of the pore ( Figure 5, right part). Since these fragments are partially located inside a sphere of diameter D, the calculated shape parameters are smaller than the analytical parameters for a continuous medium.
The shape parameters are also calculated for the three ellipsoidal pores. The principal axes of all ellipsoids are directed along the coordinate system axes (Figure 3), the value of β is equal to 30°.
Ellipsoids are placed inside a cluster so as to avoid their overlapping. To diverse the shape of the ellipsoids, values of their semi-axes are taken equal to ax1;2;3 = 8; 7; 4 nm, ay1;2;3 = 4; 5; 5 nm, az1;2;3 = 2; 3; 6 nm. The dependences of shape parameters ratios on the number of probe points are shown in the lower part of Figure 5. These ratios converge fast to the theoretical values with an increase in the number of probe points. Some differences between the numerical and analytical results are explained in the same way as for the spherical pore. To apply the developed approach to the silicon dioxide films, the values of the atomic radii of silicon and oxygen atoms should be determined. As follows from Equation (5), the components of the dielectric constant and refractive index tensors depend on the free volume fraction of the film. The free volume is defined as volume outside spheres centered on the atoms of the cluster. Therefore, the free volume is fully determined by the Cartesian coordinates of the atoms and values of the atomic radii. On the other hand, the dependence of free volume fraction on the deposition angle is determined by Equation (6). Thus, the value of the atomic radii should be chosen so as to reproduce this dependence.
Earlier in the porosity calculation [34] the van der Waals radii of oxygen and silicon atoms were taken equal to RvdW(O) = 0.152 nm [35] and RvdW(Si) = 0.21 nm [36]. However, these radii leads to f1(α) values significantly higher than those obtained using Equation (6), see the bottom row in Table 1. This is due to the large number of small pores appearing between adjacent atoms if atomic spheres do not cover the entire volume of the film. Thus, the values of the atomic radii should be increased in comparison with the values of van der Waals radii. It was found that radii equal to 0.21 nm for The convergence of the pore shapes parameters in the thin film clusters with an increase in the number of probe points N is shown in Figure 6. A value of N is approximately ten thousand, which is enough to calculate the values of ax(y,z) . This calculation takes about ten minutes for clusters consisting of 10 6 atoms. Since the algorithm has linear scaling with the number of atoms, it takes about one hour for clusters with 10 7 atoms. As can be seen in Figure 6, an approximate estimation of ratios of ax(y,z) values can be done much faster. In any case the calculation of the of ax(y,z) values requires much less time than the atomistic simulation of the deposition process using MD or kMC methods. Besides, the algorithm can be easily parallelized since the probe points in different parts of the clusters can be chosen independently. Dependencies of the shape parameters and depolarization factors on the orientation of the coordinate system axes are shown in Figure 7. Due to the symmetry of the deposited structure, the differences in these values achieve maximum when the Z axis is directed almost parallel to the slanted columns (see the right part of Figure 7). This orientation of axes ensures the maximum value of az close to 4 nm (see the left part of Figure 7). On the contrary, with this axes' orientation, the value of ax reaches its minimum, since the X axis is directed almost perpendicular to the elongated pores. Value of ay changes insignificantly with an increase of rotation angle β.
The dependencies of the depolarization factors on the angle β are explained in the same way. In accordance with the theoretical results [9,16], the sum of all factors is constant and equal to 1. At a rotation angle βm, the factor Lx is maximum, since the ratio of ax/az is minimal, see Equation (3). The column tilt angle  can be approximately calculated as 90° − βm (see the right part of Figure 7), which gives about 50° for . The tangential rule [17] predicts  = 53° for the deposition angle  = 70°, which is close to our value of 50°. The results of the refractive index calculation are shown in Figure 8 and in Table 2 and are compared with experimental data known from the literature. For α = 50° the ∆n value is less than 0.003, i.e., is at the level of anisotropy of the normally deposited SiO2 films [32]. The dependences of the refractive index on the rotation angle β for the deposition angles  = 60°, 70°, 80° are similar. The difference ∆n = nz -nx reaches its maximum at a certain value of the rotation angle β = βm. The value of βm increases with increasing of . This is explained by an increase of the tilt angle of the separated columns in GLAD films with the increase in the deposition angle [3]. Since the pores are elongated along these columns, the angle between the substrate plane and Z semi-axis of the ellipsoids also increases.
As expected, the ∆n value monotonically increases with increasing the deposition angle, due to the increase in the film porosity. The difference between the nz and nx values reaches a maximum value of 0.035 when  is equal to 80°, and βm is approximately 55°. The calculated values of ∆n are less than experimental ones for all of the deposition angles ( Table  2). This can be explained as follows. As it is mentioned in Section 2.2, dimensions of the open pores in the vertical direction are limited by the upper boundary of the film. This leads to the underestimation of the az/ax ratio calculated for the coordinate system orientation corresponding to the rotation angle βm (Figure 7). To take into account, the contribution of the large elongated pores to the az/ax ratio and ∆n, the simulation with clusters of larger dimensions should be performed.
The calculated values of the refractive index n = (nx + ny + nz)/3 are slightly higher than the experimental ones (Table 2). A possible reason for this is an underestimation of the porosity of the clusters deposited in our MD simulation runs. This also leads to the underestimation of the calculated ∆n values compared to the experimental data, see the paragraph above.
The error in the refractive index in the experiments cited in Table 2 is about 0.005 [6]. When modeling, various sets of random numbers were used in the Monte Carlo method. This led to a difference in the refractive index values less than 0.002. The convergence of the calculated values with an increase in the number of probe points was ensured in our simulation experiments, see plots in Figure 6.
To summarize, the reached level of the simulation accuracy is close to that reported in Reference [14]. To increase the accuracy of the developed approach, it is necessary to increase the dimensions of simulation clusters. This requires the use of more high-performance computational recourses. It should be emphasized that the model does not have any adjustable parameters used to reproduce the experimental results.

Conclusion
In this study, the combined modeling of the optical anisotropy of porous silicon dioxide thin films was carried out. Atomistic clusters representing the structure of the films are obtained using the full-atomistic MD simulation of the deposition process. The pores in films structure are approximated by a set of ellipsoids having the same shape parameters and the same orientation. These parameters are calculated in the framework of the Monte Carlo based method, and are further used in the anisotropic Bruggeman effective medium theory to determine the difference of the main components of the refractive index ∆n.
The ∆n values are calculated for the high-energy deposited clusters. Deposition angles vary from 0 to 80, which correspond to the normal and glancing angle deposition, respectively. Normally deposited film is considered as a homogeneous structure without void volume. The increase in the deposition angle leads to the decrease of the film density, due to pore formation. The maximum value of the difference of the refractive index main components is equal to 0.035. While the reached level of the simulation accuracy is close to that reported earlier [14], the calculated values of ∆n are still somewhat lower than the experimental ones. Possible reasons for this are discussed, and it is pointed out that even larger atomistic clusters should be considered to increase the accuracy of simulation experiments.
The developed method of combined modeling of the anisotropic properties of atomistic clusters can be applied to other oxide films. For the calculation of ∆n values, information is required about the deposition conditions-including the energy of the incoming atoms, the value of the deposition angle, and the substrate temperature.

Conflicts of Interest:
The authors declare no conflict of interest.