FDTD Simulations of Shell Scattering in Au@SiO2 Core–Shell Nanorods with SERS Activity for Sensory Purposes

The article describes the results of Finite-Difference Time-Domain (FDTD) mathematical modeling of electromagnetic field parameters near the surfaces of core–shell gold-based nanorods in the Au@SiO2 system. Three excitation linewidths (λ = 532, 632.8, and 785 nm) were used for theoretical experiments. Electric field parameters for Au nanorods, Au@SiO2 nanorods, and hollow SiO2 shells have been calculated and evaluated. The correlations between electric field calculated parameters with nanorod morphology and shell size parameters have been clarified. The optical properties of nanoobjects have been simulated and discussed. The highest maximum calculated value of the electric field tension was E = 7.34 V/m. The enhancement coefficient was |E/E0|4 = 3.15 × 107 and was obtained on a rod with a SiO2 shell with dimensional parameters of height 70 nm, rod width 20 nm, and shell thickness 20 nm. As a result, a flexible simulation algorithm has been developed for the simulation of electric field parameters in each component of the Au@SiO2 system. The developed simulation algorithm will be applicable in the future for any other calculations of optical parameters in any similar component of the core–shell system.


Introduction
Gold nanoparticles (NPs) are currently widely used for modern sensory purposes [1][2][3]. Such NPs can be modified with various compounds and shells to obtain additional modalities. Gold NPs can be modified with antibodies [4] for targeted delivery [5], shells for encapsulating drugs [6], fluorescent and Raman labels to perform labeled fluorescence, and Raman spectroscopy [7]. Gold NPs are also used in the latest applications in physics and medicine, which represent the studies on the delivery of multimodal NPs to a target and simultaneously perform therapeutic, sensory, and targeted delivery functions [8]. In this case, the NP often has a shell that performs various functions of the structure. A coating shell can enhance colloidal stability, decrease toxicity, and allow for further functionalization of the NP to form a theranostic complex. Inorganic silica (SiO 2 ) is widely used as a capsuling material due to its chemical inertness, optical transparency, and excellent shell thickness control [9]. Since Au NPs and Au@SiO 2 NPs can act as energy coupling between free electrons and fluorophores, such structures can be applied for sensory purposes. The problem is to find the optimal shape of NPs and shells, in which the values of surface plasmon resonance will have the highest values. The former can be interpreted in the sphere of sensory functions of these NPs. Most of the papers focus mainly on forms of NPs [10] or represent the results of core-shell NPs morphological and optical properties investigation [11][12][13][14]. However, such works, especially those related to purely theoretical modeling, are not sufficient, especially those with a large comparative sample within geometrically identical NPs and their shells. There are a huge number of studies on surfacefunctionalized NPs dealing with their application to models of delivery systems, in which various parameters of multilayer NPs were indicated or slightly varied [15,16]. Optical

FDTD Approach
We used an approach based on the FDTD method using the basic Yee algorithm to solve Maxwell's equations numerically. Yee proposed a spatial rectangular grid to discretize the selected computational domain [20]. The electric fields are located along the boundary of the Yee-described cube block [20], while the magnetic fields are located toward the center of the block. All grid components are spaced and independent of each other, thereby they will have different parameters at each point. The For non-magnetic materials, the first two curl-equations have the following form [21] [Equation (2)]: where → H, → E and → D denote magnetic, electric field, and offset field, respectively, and ε r (ω) = n 2 where n is refractive index. In three dimensions, the Maxwell equations have six components: E x , E y , E z and H x , H y , and H z . If we assume that the structure is infinite, for example, in the z plane and that the fields are independent of z, in particular, we obtain [Equation (3)]: ε r (ω, x, y, z) = ε r (ω, x, y,) The FDTD method uses finite differences as approximations to both spatial and temporal derivatives that appear in Maxwell's equations (in particular, Ampère's and Faraday's laws). To move from an analytical solution to a numerical one, we consider the Taylor series expansion of the function f(x) expanded with respect to the point x 0 with a shift of ±δ = 2. Next, we combine the expansion into a Taylor series with the value +δ and −δ, divide by the error δ, and create a finite difference scheme [Equation (4)]: If δ is small enough, a reasonable approximation to the derivative can be obtained by simply neglecting all subsequent terms in the series [Equation (5)]: The central difference provides an approximation of the derivative of the function at the point x 0 , to be precise, the function is chosen at the neighboring points x 0 + δ/2 and x 0 − δ/2 to approximate the vector function in Maxwell's equations. Since the smallest power of the ignored δ is of the second order, this means that the central difference is of the second-order precision. Let us move on to the main governing equations when constructing the FDTD algorithm for a three-dimensional grid. These equations [Equations (6) and (7)] are: where, ε-permittivity, µ-magnetic permeability, σ-electrical conductivity, and σ m -magnetic conductivity.
The three-dimensional mesh for modeling consists of six nodes [Equations (8)- (13)], which can be conventionally denoted by the indices m, n, and p.
Parts of a 3D cell are shown in Figure 1. This type of image is called a Yee cell. This cell consists of electric field nodes located along the edges of the cube and magnetic field nodes located on the faces. Parts of a 3D cell are shown in Figure 1. This type of image is called a Yee cell. Thi cell consists of electric field nodes located along the edges of the cube and magnetic field nodes located on the faces. Figure 1. Nodes in a 3D FDTD mesh as a Yee cell. In this image, not all nodes have the same indexes As shown here, the cube will consist of four Ex nodes, four Ey nodes, and four Ez nodes, i.e., th electric fields are located along the edges of the cube. The magnetic fields are on the faces of th cube and hence there will be two Hx nodes, two Hy nodes, and two Hz nodes.
With the arrangement of nodes shown in Figure 1, the components of Equations (6 and (7) can be expressed at the appropriate evaluation in Equations (14)- (19): These equations ignore instantaneous losses and currents. The time derivative o each electric field component is always determined by the spatial derivative of the tw magnetic field components and vice versa. In addition, the components of one field ar associated with two orthogonal distributions of the components of the other field. As ha been completed previously, the loss term can be approximated by the average of the field in two-times steps. Update equations for 3D meshes can be written by simply checkin the underlying equations in continuous representation. Update equations in the final form [Equations (20)- (25)] are as follows: As shown here, the cube will consist of four E x nodes, four E y nodes, and four E z nodes, i.e., the electric fields are located along the edges of the cube. The magnetic fields are on the faces of the cube and hence there will be two H x nodes, two H y nodes, and two H z nodes.
With the arrangement of nodes shown in Figure 1, the components of Equations (6) and (7) can be expressed at the appropriate evaluation in Equations (14)- (19): These equations ignore instantaneous losses and currents. The time derivative of each electric field component is always determined by the spatial derivative of the two magnetic field components and vice versa. In addition, the components of one field are associated with two orthogonal distributions of the components of the other field. As has been completed previously, the loss term can be approximated by the average of the field in two-times steps. Update equations for 3D meshes can be written by simply checking the underlying equations in continuous representation. Update equations in the final form [Equations (20)- (25)] are as follows: The coefficients in the update equations are assumed to be constant over time, but as in the case of our calculations, they can be functions of the position. In accordance with the accepted notation, and assuming a uniform grid (wherein, the grid in software simulation can be of different accuracy), in which ∆x = ∆y = ∆z = δ, the magnetic field update coefficients can be expressed as [Equations (26)-(31)]: For the electric-field update equations, the coefficients are in the following form [Equations (32)-(37)]: C eze m, n, p C ezh m, n, p + These coefficients can be solved with the Courant number, which is related to the convergence condition of this equation. For a uniform grid in three dimensions, the Courant limit must be equal to 1/ √ 3. We will not touch upon a rigorous derivation of this limit. Let us use a simple empirical explanation. To ensure stability, we need to make sure that the distance traveled in the real continuous world in three-time steps is less than the distance over which the grid can transmit information. The grid must have time to calculate the propagation of the same perturbation no slower than if this perturbation took place in real time. In our case, we use the most accurate standard grid (not counting the additional refinement grid).

Raman Scattering and Enhancement Factor
The Raman scattering effect is the result of inelastic scattering between a photon and the vibrational modes of a molecule. The power of the scattered Raman signal can be described by the expression [Equation (38)]: where N denotes the number of active scatterings within the propagated excitation, σ RS denotes the scattering cross section, I(ν L ) is the intensity of the incident beam at the frequency ν L , and v s is the frequency of the scattered Raman signal. This expression takes into account only the Stokes shift, in which case the Raman signal will be less than the frequency of the incident light [22]. However, in Raman spectroscopy, the scattering intensity can be put into a linear dependence on the intensity of the incident field, E 2 0 . The Stokes shift occurs when an incident photon interacts with a molecule in its ground vibrational state. Near the surface of scattering metal nanostructures, this process is accompanied by signal amplification and is known as SERS [23]. Since the magnitude of the field intensity on these surfaces is significantly increased, the intensity of Raman scattering can be related to the absolute value of the square E out on the surface of the NPs. Let us give, as an example, the equation for a metal sphere [Equation (39)] [24]: where g denotes a value that depends on the dielectric constant of the medium and the metal NP and θ is the angle between the incidence field vector and the vector directed to the position of the molecule on the surface. Peak amplification occurs when θ is 0 • or 180 • , which corresponds to the position on the axis of light propagation. In cases where g is large, the maximum gain approaches [Equation (40)] [24]: We take into account the contribution of the applied field to Raman scattering which, as noted above, induces a vibrational dipole in the molecule on the surface. Then, this dipole radiates, and in the approximation there is a probability that the emitted light is, as in theory, shifted by the vibrational frequency of the molecule. The first order approximation is to use an expression similar to Equations (14)- (19), except that it is evaluated at the Raman-Stokes bias frequency, so the following expression can be written in approximation [Equation (41)] [24] as follows: This expression is defined as the theoretical gain (EF) of the SERS. In the literature, this expression is called the fourth power of field enhancement on the surface of NPs |E| 4 . More details on the use of the coefficient |E| 4 can be found, for example, in these works [24,25].
This section shows the main essence of the use of EF SERS related directly to the interpretation of experimental data. In our work, however, a slightly different theoretical model of EF is used, as described in Section 2.3, not including the presence of the analyte and depending largely on the strength of the electric field and the incident wave.

Simulation Process
Modeling was performed using the Lumerical FDTD Solutions software package [v. 8.19.1584, Lumerical Inc. (Vancouver, BC, Canada)]. The simulation was carried out for rod-shaped Au NPs (NRs). The NRs were varied radii of 5, 10, 15, and 20 nm. NRs' lengths were taken as 70 and 80 nm. Moreover, NRs of the same dimensions were modeled using an encapsulating SiO 2 core-shell shell with shell thicknesses of 2, 3, 5, 10, 15, and 20 nm. Schematic illustrations of NRs FDTD simulations are illustrated in Figure 2.
We take into account the contribution of the applied field to Raman scattering which, as noted above, induces a vibrational dipole in the molecule on the surface. Then, this dipole radiates, and in the approximation there is a probability that the emitted light is, as in theory, shifted by the vibrational frequency of the molecule. The first order approximation is to use an expression similar to Equations (14)- (19), except that it is evaluated at the Raman-Stokes bias frequency, so the following expression can be written in approximation [Equation (41)] [24] as follows: This expression is defined as the theoretical gain (EF) of the SERS. In the literature, this expression is called the fourth power of field enhancement on the surface of NPs | | . More details on the use of the coefficient | | can be found, for example, in these works [24,25].
This section shows the main essence of the use of EF SERS related directly to the interpretation of experimental data. In our work, however, a slightly different theoretical model of EF is used, as described in Section 2.3, not including the presence of the analyte and depending largely on the strength of the electric field and the incident wave.

Simulation Process
Modeling was performed using the Lumerical FDTD Solutions software package  The stages of the simulation process were performed as follows: (1) The counting area, grid resolution, and boundary conditions were set. For the computational domain, a rectangular grid was used from the basic Yee algorithm in the Cartesian coordinate system. The main modeling quantities (material properties and object geometry, electric, and magnetic fields) were calculated separately at each grid point. The size of the computational area along the axes varied within different limits based on the size of the objects. To maintain accuracy, the meshing algorithm generated a smaller mesh with a high index (to maintain a constant number of mesh points per wavelength). The minimum grid spacing was set to 0.25 nm. Then, an additional refinement mesh was installed for the simulation. The size of the computational region of the additional grid was set by the grid step: dx, dy, and dz = 2.5 nm. We have chosen standard absorbing perfectly matched layer (PML) and boundary conditions designed to absorb incident light with minimal reflections. Their parameters were as follows: − Layers (for PML area discretization purposes) = 8. (2) A body with specified optical and geometrical parameters was placed inside the counting region. Next, the optical and geometric parameters of the samples were set. We used materials from a Lumerical database (Au, SiO 2 ) and changed their parameters (size, shape, and geometry) for modeling. The values of such a parameter as the real and imaginary parts of the permittivity, which depends on the radiation frequency, were taken into account. In our case, a plane wave equal to λ = 532, 632.8, and 785 nm. We also used theoretical ε (Au) values for modeling, taken from the Lumerical database (Table 1). (3) At the next step, the radiation source parameters were set for three wavelengths. In our study, a total scattered field source (TFSF) is used, which is often suitable for studying scattering by small particles illuminated by a plane wave. The TFSF source divides the computational domain into two separate domains: (a) the total field domain includes the sum of the incident field wave plus the scattered field, and (b) the scattered field domain includes only the scattered field. The TFSF source is an extended source. It is important to note that the physical field is a total field, and the division into the incident and scattered fields requires careful interpretation. For NPs in a homogeneous medium, the incident field is a p-polarized plane wave. We obtained the magnitude of the electric tension in the maximum values. We also calculated the array of electric field values in the region of plasmon generation defined by us as the integral sum ∑ k p E attenuation point of the propagated field. It should be noted that the position of the source of a plane ppolarized wave is important for calculating the maximum electric field value, especially under conditions of single measurements, but this is not decisive for the overall distribution of the dependence in a rather large sample. We have installed the radiation source close to the NR at 3/2 of its length.
(4) To provide final information, the monitor plane was set perpendicular to the xzplane, which gave us the final information about the value of the electric field as a function of position in space, in the form of a 2D slice. The use of field monitors in the frequency domain allowed us to collect a field profile in this region and provide simulation results in some spatial domains to the FDTD solver. The wave was polarized along the z axis, i.e., its direction corresponded to the normal vector to the monitor surface on the XZ axis.
(5) As the last step, the calculated values of the electric field were converted using the Lumerical program code (scripts) into the intensity values of the Raman and SERS according to the theoretical effective |E/E 0 | 4 or enhancement factor (EF) for SERS calculated in the xz-plane. In order to do this, we use TFSF sources and models in this area of the structure to find the maximum values of EF. Even though a SERS calculation is performed mostly for rough surfaces or dimers on which it is calculated at hot spots between gaps, we modeled it for sheathed NRs. Comparing such NRs with NRs without shells, we can check the LSPR presence in NPs by the presence of a dielectric shell. The theoretical gain has been calculated for NRs with/without shells. This parameter, in contrast to experimental measurements, can be calculated without conditions for the presence of the analyte.
Let us take into account the absence as such of a strict dependence between E as the magnitude of the electric field strength of localized surface plasmons, which for a TFSF source determines mainly the magnitude of the total plus scattered field, to the SERS scattering intensity and EF coefficient, which depend on the resonant absorption frequency of various gold NPs. It is possible that the integral sum of the electric fields distributed near the surface of the NRs will approximate this dependence. The integral sum of field ∑ k p E is defined in relative terms, as it depends on the NR radius. The SERS and EF values themselves are equivalent in our theoretical modeling, i.e., one value of SERS will correspond to only one value of EF. However, in experimental studies, there are inconsistencies between these values due to the lack of strict definitions [26].

NR without a SiO 2 Shell
As a result of the first part of the simulation, the values of the local maximum of the electric field were obtained as a function of SERS and EF for Au NRs with a thickness of 5, 10, 15, and 20 nm. Depending on the NR radius, three radiation wavelengths were used: 532, 632.8, and 785 nm. NRs of two lengths (70 and 80 nm) were applied in experiments. This is due to insignificant differences in the studied values when varying this parameter, since all parameters were calculated in the plane parallel to the NR length (xz). In this line of modeling, we did not calculate the integral sum of the fields, which is mostly comparative in nature with a large sample. The data are summarized in Table 2. It can be seen that the NRs give low values of electric field strength, mainly increasing with NR radius growth. A similar increase also occurs in the analysis of identical NRs with an increase in the radiation wavelength. The main results for NR with a 70 nm and 80 nm length for different excitation sources are shown in Figure 3.
The SERS intensity and EF appeared to increase with the NR radius, which was manifested significantly at a 785 excitation wavelength. That fact correlates with [27] data and depends on plasmon excitation maximum near 785 nm. 20 3.55 6.29 6.74 520 530 3600 26.6 27.5 1320 It can be seen that the NRs give low values of electric field strength, mainly increasing with NR radius growth. A similar increase also occurs in the analysis of identical NRs with an increase in the radiation wavelength. The main results for NR with a 70 nm and 80 nm length for different excitation sources are shown in Figure 3. The SERS intensity and EF appeared to increase with the NR radius, which was manifested significantly at a 785 excitation wavelength. That fact correlates with [27] data and depends on plasmon excitation maximum near 785 nm.

NR with a SiO2 Shell
The second series of modeling experiments were carried out with the Au@SiO2 complex using a xz monitor plane. Below are the data obtained for Au NRs with radii of 5, 10, 15, and 20 nm and a SiO2 shell of various thicknesses (2, 3, 5, 10, 15, and 20 nm). The same parameters were studied for NRs without a shell, to which the study of the integral sum of fields was added. The effect of different shell thicknesses, as expected, introduced a nonlinear character into the quantities under study (Table 3).

NR with a SiO 2 Shell
The second series of modeling experiments were carried out with the Au@SiO 2 complex using a xz monitor plane. Below are the data obtained for Au NRs with radii of 5, 10, 15, and 20 nm and a SiO 2 shell of various thicknesses (2, 3, 5, 10, 15, and 20 nm). The same parameters were studied for NRs without a shell, to which the study of the integral sum of fields was added. The effect of different shell thicknesses, as expected, introduced a nonlinear character into the quantities under study (Table 3).  When analyzing the results obtained from the electric field simulation for NRs with a SiO 2 shell, it can be seen that the electric field values for three different excitation wavelengths (532, 632.8, and 785 nm) showed a similar distribution as the "bare" NRs. An increase occurred for relatively thin NRs (5 and 10 nm radii), smoothing out the differences at 632.8 and 785 nm for thick NRs (15 and 20 nm radii). For SERS and EF, the values were approximately equal for the same radii at 532 and 632.8 nm for thin NRs (5 and 10 nm radii) increasing in difference in favor of 532 with radius growth. Moreover, SERS/EF values appeared to become higher at the 785 nm excitation wavelength. Analyzing the influence of the SiO 2 shell, we did not observe a clear monotonous character: electric field and SERS signal intensities fell at medium values of SiO 2 thickness. First, in comparison with bare NRs, shelled NRs excited at wavelengths of 532 and 632.8 nm had lower electric field strengths for thin-coated NRs and higher values for thick NRs coated with a thick SiO 2 shell. This can be explained by the fact that the thickness of the SiO 2 shell in the case of NRs excited by laser radiation of the aforementioned wavelengths affects the refractive index by changing the extinction coefficient. This can explain the relatively low values of maximum electric field strength at an average shell thickness. For the 532 nm excitation wavelength, the largest values of the electric field almost always corresponded to the relatively largest values of SERS. No such correspondence was observed for the 632.8 and 785 nm excitation wavelengths. This can be explained by the fact that such sizes of sheathed gold rods can fall into absorption peaks in the long wavelength region and contribute to SERS, but they contribute less to the scattered field. For the NRs irradiated by the 785 nm laser emission, an increase in the thickness of the shell leads to an increase in the difference between the maximum values of the electric field and the relatively maximum SERS. This is consistent with the fact that such sizes of sheathed gold NRs fall into absorption peaks in the long wavelength region and contribute to SERS, but do not contribute to scattering.

A Single (hollow) SiO 2 Shell
The most interesting case was the simulation of the silica shell scattering without the NR (single, or hollow NR). The quantitative contribution was determined for the scattering of the silica shell in the overall pattern of changes for the simulated parameters of the electric field. The obtained results show that silica shells contribute to the scattering of the overall system at all excitation wavelengths used in our simulation. These results (Table 4) correlate with [23,28].
Based on the results of the paper, it was found that the silicon shell, although it contributes to scattering, cannot be used as an independent medium for a significant increase or decrease in the scattering of structures adsorbed on it. In addition, as the shell thickness increases, the scattering decreases, from which we can conclude that the overall efficiency of Au@SiO 2 sensors decreases.

Calculation of the Integral Sum for a Single (Hollow) SiO 2 and NR Au@SiO 2
The maximum values of the electric field integral sums almost always correspond either to the maximum values of E, which is logically explained by the fact that one of the values taken will be included in the largest, or to the high values of SERS, which can be explained by the fact that SERS itself is nothing more than the sum of the peaks of the Raman signal shift.
The simulations of the electrical field integral sum distributed near the surface of the NRs and shell demonstrate strong E presence in the system. At 632.8 and 785 nm excitation wavelengths (Table 5), electrical field strength values are higher than at 532 nm (for both 70 and 80 nm NR lengths). For 5 and 10 nm inner shell radii, E decreases in Au@SiO 2, but in SiO 2 , E did not show any obvious correlations with NR or shell size parameters.

Conclusions
As a result, a flexible simulation algorithm based on the FDTD method has been developed for the simulation of electric field parameters in each component of the Au@SiO 2 system. It has been found that the SiO 2 shell on Au NR, although it contributes to the total system scattering, demonstrates low values. Some cases linearly increased from the thickness of the SiO 2 shell at the laser excitation wavelength of λ = 532 nm and λ = 632.6 nm, because for λ = 785 nm laser excitation there is an absorption peak. The highest value of the E = 6.29 V/m was revealed for nanorods without a shell, excited by a monochromatic wave λ = 785 nm using dimensional characteristics: length 80 nm and radius 20 nm. The parameter |E/E 0 | 4 = 1.41 × 10 7 was obtained for a rod with almost the same parameters, except for the length (l = 70 nm). On the hollow SiO 2 shells, the highest E value (E = 2.81 V/m) was obtained with the simulation parameters of length 70 nm, inner radius 20 nm, and shell thickness 20 nm, upon excitation with a wavelength of λ = 785 nm. The highest value for the |E/E 0 | 4 = 1.8 × 10 5 was obtained with the simulation parameters of length 70 nm, inner radius 5 nm, and shell thickness 15 nm, with excitation at a wavelength of λ = 785 nm. For the SiO 2 coated rod, the highest maximum values of the electric field strength (E = 7.34 V/m), as well as the |E/E 0 | 4 = 3.15 × 10 7 ), when excited by laser radiation λ = 785 nm, with dimensional parameters were length 70 nm, rod width 20 nm, and shell thickness 20 nm. As we can see, this simulation object had the highest values of all. As a result, the electric field parameters of Au NR, Au@SiO 2 , and a single (hollow) SiO 2 shell were calculated and estimated. The results of the work can be used in methods for the synthesis of core-shell NPs.