A Quantum Chemistry Approach Based on the Analogy with π-System in Polymers for a Rapid Estimation of the Resonance Wavelength of Nanoparticle Systems

In this paper, the Variational Method based on the Hückel Theory is applied to NPs chain and aggregate systems in order to estimate the energy of the plasmon and, in turn, the resonance wavelength shift, which is caused by the interaction of adjacent NPs. This method is based on the analogies of NPs dipole interactions and the π-system in molecules. Differently from the Hartree-Fock method that is a self-consistent model, in this approach, the input data that this method requires is the dimer energy shift with respect to single NPs. This enables us to acquire a simultaneous estimation of the wavefunctions of the NPs system as well as the expectation energy value of every kind of NPs system. The main advantage of this approach is the rapid response and ease of application to every kind of geometries and spacing from the linear chain to clusters, without the necessity of a time-consuming calculation. The results obtained with this model are closely aligned to related literature and open the way to further development of this methodology for investigating other properties of NPs systems.


Introduction
The understanding of light-matter interaction at a nanoscale has been developed rapidly in recent years. Various models [1][2][3] have been proposed to understand and optimize several applications such as Surface Enhanced Raman Scattering (SERS) [4], photo-emission [5], fluorescence [6], Nanoparticle Enhanced Laser Induced Breakdown Spectroscopy (NELIBS) [7], and more. It is well known that, when an incident radiation illuminates a metal nanoparticle (NP), it induces polarization of the NPs conduction electrons and, if the radiation frequency is in resonance with the plasmonic frequency of the particle itself, a strong electromagnetic field is induced on the border of the NPs [8]. If the radiation with adequate frequency interacts with a set of NPs placed at a distance, which is smaller than the diameter of the NP, a strong field enhancement is generated in the gap between the NPs, i.e., the so-called "hot spot" [9,10]. In this case, an effective variation of the Coulomb restoring force that exists between the positive metallic lattice of the NPs and the oscillating electrons occurs because of the interaction of adjacent NPs dipoles. This results in a notable change of the system's resonance energy. Several experiments and models have been performed in order to quantify this phenomenon and, thus, determine the resonance wavelength of the NPs system [11,12]. Although various compositions of multi-NPs have been investigated [13,14], due to the mathematical approach complexity and the Let us consider a general wavefunction of a linear system of n equally sized and spaced NPs given by the following linear combination. ψ k = c 1 ϕ 1 + c 2 ϕ 2 + c 3 ϕ 3 + . . . + c n ϕ n (1) where ϕ is the wavefunction associated with the dipole oscillation of an excited individual NP. If we use n different basis for the linear combination, n different wavefunctions will be obtained, with each one describing a different mode. In this case, we will consider only longitudinal modes, restricting our attention to those modes which are susceptible when the incident radiation is perpendicular to the NPs cluster.
(Ω 11 − ω) (Ω 12 where (a) the Coulomb's integral Ω ii = τ ϕ iω ϕ i dτ = ω 0 that is the energy of the single particle, (b) the interaction integrals Ω ij = Ω ji = τ ϕ iω ϕ j dτ = τ ϕ jω ϕ i dτ = β ≈ −∆ω dim is equal to the shift of energy of the dimer compared to the single NP when i and j represent adjacent NPs and is null for non-adjacent particles because we assume that non-adjacent dipoles do not interact between them. In the case of a linear chain of NPs, Ω ij = β when i − j = 1 and is equal to 0 when i − j > 1; (c) the overlapping integral S ij = S ji = τ ϕ i ϕ j dτ = τ ϕ j ϕ i dτ ≈ 0 because the NPs are far enough to neglect the contribution of the efficient wavefunctions overlapping.
Under these conditions, for the system of equation to be null, the secular determinant must be numerically equal to 0. As an example, in the case of a linear chain of n NPs, the following condition must be assumed.
Some other examples are reported in the supplementary materials Figure S1. The solution of this determinant for n particles, i.e., n wavefunctions, gives n energy modes, including some that are bright modes (ω < ω 0 ) while the rest are dark modes (ω > ω 0 ). Some other examples are reported in the supplementary materials Figure S1. The solution of this determinant for n particles, i.e., n wavefunctions, gives n energy modes, including some that are bright modes ( < 0) while the rest are dark modes ( > 0).
A representative sketch in the case of a chain of four equal NPs is shown in Figure 1. As a result of this calculation, we obtain several wavefunctions describing several dipole configurations, as shown in the example reported in Figure 1.
It is worth underlining that, with the present approach, whenever the plasmon energy of the single NP and the shift of energy of the dimer with respect to the single NP is known, the calculation of the secular determinant of Equation (5) is a minor and immediate operation. Thus, this method shows a general trend of the NPs system where the effect of NPs distances is included in the interaction integral whose value can be obtained by using an independent theoretical or experimental approach. The exact wavefunction of the single NPs is not required for the computation of the energy in agreement with the Hückel theory and so the method can be applied virtually to all the geometries of the single element as the effect of the properties of the single NPs in the complex system, which is contained in the input data 0 and dim. As is well known, dielectric polarization plays an important role in the plasmon resonant energy. For this reason, the input data should be selected properly, taking into account the dielectric medium where the NPs system is embedded.
Dimeric systems have been intensively studied in the last decade because the optical properties of this kind of system can be determined accurately using different experimental and theoretical approaches [21][22][23][24][25][26][27][28][29][30]. For instance, we reported in Table 1 some data retrieved from published works on Au and Ag NPs dimers. As a result of this calculation, we obtain several wavefunctions describing several dipole configurations, as shown in the example reported in Figure 1.
It is worth underlining that, with the present approach, whenever the plasmon energy of the single NP and the shift of energy of the dimer with respect to the single NP is known, the calculation of the secular determinant of Equation (5) is a minor and immediate operation. Thus, this method shows a general trend of the NPs system where the effect of NPs distances is included in the interaction integral whose value can be obtained by using an independent theoretical or experimental approach. The exact wavefunction of the single NPs is not required for the computation of the energy in agreement with the Hückel theory and so the method can be applied virtually to all the geometries of the single element as the effect of the properties of the single NPs in the complex system, which is contained in the input data ω 0 and ∆ω dim . As is well known, dielectric polarization plays an important role in the plasmon resonant energy. For this reason, the input data should be selected properly, taking into account the dielectric medium where the NPs system is embedded.
Dimeric systems have been intensively studied in the last decade because the optical properties of this kind of system can be determined accurately using different experimental and theoretical approaches [21][22][23][24][25][26][27][28][29][30]. For instance, we reported in Table 1 some data retrieved from published works on Au and Ag NPs dimers.

NPs Chain Systems
Let us consider the case of a chain of AuNPs and one of AgNPs, where the plasmon energy ω 0 and the dimer shift ∆ω dim are ω 0 = 2.375 eV, β = −0.28 eV and ω 0 = 3.10 eV and β = −0.52 eV, respectively, for AuNPs of 40 nm and AgNPs of 50 nm [22]. Results of the VM are reported in Figure 2, which shows an increase of the wavelength as a function of the number of NPs in the chain, which, in turn, results in the red-shift of the resonance wavelength of the lowest energy mode, obtained with the solution of the secular determinant of Equation (5).  Considering the broadening of each state due to the quantum limit of the investigated system [34], it should be expected that all these states merge together in a wide band. Consequently, the system of NPs can be suitably polarized along a broad range of incident electromagnetic field frequencies generating a band of resonance including all the bright modes within the limits of possible configurations, which the dipoles in the chain can assume. Since the VM does not allow taking into account the broadening of the energy levels resulting from the quantum limits and a detailed investigation on the broadening would be required [35,36], with the only aim to show how the broadening of the states may form a band of resonance, in Figure 3, we also report the convolution of the modes, assuming a Gaussian shape with a broadening (FWHM, Full Width at Half Maximum) of 80 nm in agreement with the broadening of the SPR (Surface Plasmon Resonance) of the single AuNPs of 40 nm [34]. Even though, this is an arbitrary condition for a chain system and As expected, the wavelength of the longitudinal modes shifts to higher values as the number of NPs in the system increases. This moves toward a plateau for n → ∞. Although, with the present model, we are not investigating the behavior of the Hot Spot [31,32] generated in the NPs system. However, it is possible to estimate just the number of Hot Spots assuming they are equal to the number of coherent interactions between adjacent NPs.
For each mode, the number of available coherent dipole interactions, i.e., potential hot spots, is m = n-k, where n is the number of NPs constituting the system and k is the mode number, assuming k = 1 for the lowest energy mode. This observation makes it clear that the maximum number of available coherent dipole interactions is obtained for the mode with the lowest energy. This mode represents the dipole configuration, which allows the maximum electromagnetic field enhancement. On the contrary, the number of sites where the dipole assumes the opposite orientation corresponds to the number N of nodes of the wavefunction representing the specific mode, so that N = n − 1 − m = k − 1. Consequently, in order to polarize a given mode, it is necessary that the number of constructive interactions between the dipoles is greater than the number of dipoles with a destructive interaction. This means that the number of coherent dipole interactions has to be greater than the number of nodes, m > N (see Figure 1). This observation suggests that, in a bright mode m > N, consequently, ω < ω 0, and the corresponding wavefunction represents a system with a net total dipole. On the contrary, when m < N, then ω > ω 0 , and the wavefunction describes a system where no effective dipole exists, which is a dark mode [33]. In the case of the dark mode, the total wavefunction presents more nodes than coherent interactions and, therefore, in agreement with the general knowledge, they cannot be excited.
Considering the broadening of each state due to the quantum limit of the investigated system [34], it should be expected that all these states merge together in a wide band. Consequently, the system of NPs can be suitably polarized along a broad range of incident electromagnetic field frequencies generating a band of resonance including all the bright modes within the limits of possible configurations, which the dipoles in the chain can assume. Since the VM does not allow taking into account the broadening of the energy levels resulting from the quantum limits and a detailed investigation on the broadening would be required [35,36], with the only aim to show how the broadening of the states may form a band of resonance, in Figure 3, we also report the convolution of the modes, assuming a Gaussian shape with a broadening (FWHM, Full Width at Half Maximum) of 80 nm in agreement with the broadening of the SPR (Surface Plasmon Resonance) of the single AuNPs of 40 nm [34]. Even though, this is an arbitrary condition for a chain system and it underestimates the real broadening, it gives an idea of the correlation between the states of the system and the peculiarities of the band of resonance generally observed in the experimental applications. Considering the broadening of each state due to the quantum limit of the investigated system [34], it should be expected that all these states merge together in a wide band. Consequently, the system of NPs can be suitably polarized along a broad range of incident electromagnetic field frequencies generating a band of resonance including all the bright modes within the limits of possible configurations, which the dipoles in the chain can assume. Since the VM does not allow taking into account the broadening of the energy levels resulting from the quantum limits and a detailed investigation on the broadening would be required [35,36], with the only aim to show how the broadening of the states may form a band of resonance, in Figure 3, we also report the convolution of the modes, assuming a Gaussian shape with a broadening (FWHM, Full Width at Half Maximum) of 80 nm in agreement with the broadening of the SPR (Surface Plasmon Resonance) of the single AuNPs of 40 nm [34]. Even though, this is an arbitrary condition for a chain system and it underestimates the real broadening, it gives an idea of the correlation between the states of the system and the peculiarities of the band of resonance generally observed in the experimental applications.  This phenomenon underlines the most important feature of the NPs systems: the possibility to adapt the configuration of their dipoles to the incident electromagnetic radiation, which allows a sort of cooperation in reaching resonance with the incoming perturbation. Clearly, this effect is not achievable with more rigid systems, e.g., the atomic ones, where resonance occurs only at a specific frequency. This flexibility of the NPs systems in reaching resonance with the perturbation occurs as an effect of a system that borders strongly quantized systems such as atoms and macroscopic (i.e., not quantized) systems such as particles.
It is interesting to compare the results of the VM with other results based on different approaches reported in literature, as shown in Figure 4. In this case, we compare the lowest energy mode obtained with VM with different properties reported in literature and related to the resonance energy of the plasmonic system. Observation of Figure 4 indicates an impressively high level of agreement of the results of VM with those obtained with Electrodynamic Simulation, S, [22] and Finit Integration Technique, FIT, [26], as well as with experimental results, E, [24]. The values of ω 0 and β used for comparisons of Figure 4 are ω 0 = 2.34 eV, β = −0.57 eV for the comparison with experimental results (EM), ω 0 = 2.36 eV, β = −0.33 eV for the comparison with Electrodynamics simulation (MS) and ω 0 = 2.37 eV, β = −0.27 eV for the comparison with the Finit Integration Technique (MFIT). The agreement of the proposed method with literature data shows the similarity of the plasmonic characteristics of NPs system with a conjugated dienic system, where this theoretical approach is considered one of the most accurate ways to determine the energy of π-bonding. This phenomenon underlines the most important feature of the NPs systems: the possibility to adapt the configuration of their dipoles to the incident electromagnetic radiation, which allows a sort of cooperation in reaching resonance with the incoming perturbation. Clearly, this effect is not achievable with more rigid systems, e.g., the atomic ones, where resonance occurs only at a specific frequency. This flexibility of the NPs systems in reaching resonance with the perturbation occurs as an effect of a system that borders strongly quantized systems such as atoms and macroscopic (i.e., not quantized) systems such as particles.
It is interesting to compare the results of the VM with other results based on different approaches reported in literature, as shown in Figure 4. In this case, we compare the lowest energy mode obtained with VM with different properties reported in literature and related to the resonance energy of the plasmonic system. Observation of Figure 4 indicates an impressively high level of agreement of the results of VM with those obtained with Electrodynamic Simulation, S, [22] and Finit Integration Technique, FIT, [26], as well as with experimental results, E, [24]. The values of 0 and  used for comparisons of Figure

NPs Cluster Systems
One of the most attracting results that can be obtained with the VM is the determination of resonance energy of a more complex system than a dimmer and, as mentioned in the introduction, the VM allows us to study the interaction between NPs in a different geometrical configuration. As

NPs Cluster Systems
One of the most attracting results that can be obtained with the VM is the determination of resonance energy of a more complex system than a dimmer and, as mentioned in the introduction, the VM allows us to study the interaction between NPs in a different geometrical configuration. As an example, let us consider the heptamer in three geometries: linear chain, ring chain, and cluster. The energy diagram obtained for these systems is reported in Figure 5. Same as in the previous section, we considered a system of Au-NPs with ω 0 = 2.375 eV and β = −0.28 eV. We can note that the linear Nanomaterials 2019, 9, 929 8 of 16 chain and ring chain had a similar energy shift, 1.86 eV and 1.88 eV, respectively. However, while in the linear chain, there was no degeneracy of the modes, the ring chain presented a degeneracy equal to 2 of all the modes except the one with the higher energy with k = 7. This effect was due to the geometrical symmetry of the ring system. When comparing the chain systems with the heptamer cluster, a further redshift was observed as a result of the strong coupling with the NP placed in the center in the heptamer cluster. It interacts with all the surrounding NPs, which increases the number of coherent interactions. The latter observation introduces a general feature of 2D systems, which is an increase of the number of bright modes with respect to the dark modes. In this frame, the heptamer shows a general condition where the NPs packing (i.e., the ratio between NP diameter and interparticle mean distances) is optimized, which allows the central NP to have a full interaction with several other units, since it will be discussed below in the section about the 2D lattice of NPs. an example, let us consider the heptamer in three geometries: linear chain, ring chain, and cluster. The energy diagram obtained for these systems is reported in Figure 5. Same as in the previous section, we considered a system of Au-NPs with 0 = 2.375 eV and  = −0.28 eV. We can note that the linear chain and ring chain had a similar energy shift, 1.86 eV and 1.88 eV, respectively. However, while in the linear chain, there was no degeneracy of the modes, the ring chain presented a degeneracy equal to 2 of all the modes except the one with the higher energy with k = 7. This effect was due to the geometrical symmetry of the ring system. When comparing the chain systems with the heptamer cluster, a further redshift was observed as a result of the strong coupling with the NP placed in the center in the heptamer cluster. It interacts with all the surrounding NPs, which increases the number of coherent interactions. The latter observation introduces a general feature of 2D systems, which is an increase of the number of bright modes with respect to the dark modes. In this frame, the heptamer shows a general condition where the NPs packing (i.e., the ratio between NP diameter and interparticle mean distances) is optimized, which allows the central NP to have a full interaction with several other units, since it will be discussed below in the section about the 2D lattice of NPs. As a second case, we investigated the effect of ordered 2D NPs clusters. It was necessary to take into account that there were different types of interaction between adjacent NPs and particles placed on the diagonal direction. Therefore, for the latter ones, we used a different value of the interaction integral, which is ' =  /2, according to geometric considerations (see supporting information). As a result of the dipole's interaction, a new energy diagram was obtained and is reported in Figure 6A. The figure shows that the energy of bright modes decreases with the increasing of NPs in the reticular cluster and that the number of bright modes is always greater than the one of dark modes. As mentioned above, this is due to the increasing number of coherent dipole interactions. For this reason, and because the total energy of the levels must be conserved, the dark modes notably increase their energy. Once again, for  < 0, the energy stabilization and, in turn, the wavelength increasing is dependent on the number of interactions that every NP is able to perform. The ordered reticular system of NPs does not represent the best situation since the packing factor can be improved by decreasing the average distances between the NPs. This effect is demonstrated in Figure 6B where the energy diagram of the 2D lattice is reported. The figure clearly shows a further decrease of the energy of about 0.1 eV corresponding in a wavelength red-shift of about 40 nm when the chain composing the 2D-lattice is shifted with respect to the adjacent chains in order to optimize the packing of the NPs. As a second case, we investigated the effect of ordered 2D NPs clusters. It was necessary to take into account that there were different types of interaction between adjacent NPs and particles placed on the diagonal direction. Therefore, for the latter ones, we used a different value of the interaction integral, which is β = β/ √ 2, according to geometric considerations (see supporting information). As a result of the dipole's interaction, a new energy diagram was obtained and is reported in Figure 6A. The figure shows that the energy of bright modes decreases with the increasing of NPs in the reticular cluster and that the number of bright modes is always greater than the one of dark modes. As mentioned above, this is due to the increasing number of coherent dipole interactions. For this reason, and because the total energy of the levels must be conserved, the dark modes notably increase their energy. Once again, for ω < ω 0 , the energy stabilization and, in turn, the wavelength increasing is dependent on the number of interactions that every NP is able to perform. The ordered reticular system of NPs does not represent the best situation since the packing factor can be improved by decreasing the average distances between the NPs. This effect is demonstrated in Figure 6B where the energy diagram of the 2D lattice is reported. The figure clearly shows a further decrease of the energy of about 0.1 eV corresponding in a wavelength red-shift of about 40 nm when the chain composing the 2D-lattice is shifted with respect to the adjacent chains in order to optimize the packing of the NPs.
In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at ω 0 = 2.375 eV, β = −0.28 eV. In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV.  In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV.  In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV.  In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV.  In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV.  In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV.  In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV.  In Table 2, for a fast comparison of resonance energy of different cluster geometries, the energy and the corresponding wavelength of the lowest energy mode (k = 1), that, as mentioned above, corresponds to the main peak of the resonance band, is reported at 0 = 2.375 eV,  = −0.28 eV. Lastly, in Figure 7, we compared the case of triangular, square, and heptamer cluster with some results published in Reference [37]. For all the examined structures, there is good agreement of the Lastly, in Figure 7, we compared the case of triangular, square, and heptamer cluster with some results published in Reference [37]. For all the examined structures, there is good agreement of the resonance wavelength obtained with Generalized Mie Theory (GMT) and the wavelength Lastly, in Figure 7, we compared the case of triangular, square, and heptamer cluster with some results published in Reference [37]. For all the examined structures, there is good agreement of the resonance wavelength obtained with Generalized Mie Theory (GMT) and the wavelength corresponding to the lowest energy mode, as calculated applying the VM.
Lastly, in Figure 7, we compared the case of triangular, square, and heptamer cluster with some results published in Reference [37]. For all the examined structures, there is good agreement of the resonance wavelength obtained with Generalized Mie Theory (GMT) and the wavelength corresponding to the lowest energy mode, as calculated applying the VM. Lastly, in Figure 7, we compared the case of triangular, square, and heptamer cluster with some results published in Reference [37]. For all the examined structures, there is good agreement of the resonance wavelength obtained with Generalized Mie Theory (GMT) and the wavelength corresponding to the lowest energy mode, as calculated applying the VM.

Figure 7.
Comparison of results of the present work and Reference [37] in the case of triangle, square, and heptamer clusters. Input data have been estimated from Reference [37] and are 0 = 2.15 eV and  = −0.28 eV.

Effect of the Surface
The effect of the metallic surface on the chain plasmon resonance is of great interest and has been investigated by several authors (see [20,[38][39][40][41]). Two approaches may be feasible in the context of the VM. Applying the first approach, we can directly couple the surface plasmon energy independently with all the single NPs constituting the chain. From the mathematical point of view, it means adding a further equation in the system as well as adding a second term, , that takes into account the energy shift due to the interaction of an individual NP with the metallic surface. In any case, this approach would take into account only the main energy shift due to a strong coupling of the metal surface and the NPs. On the contrary, the surface can interact with the NPs chain, which induces two main groups of modes, with one representing the most effective coupling and the other one's surface interacting destructively with the modes of the chain [20]. In order to better represent this phenomenon, we propose coupling directly each chain wavefunction with the metallic surface. Consequently, each single mode determined, as discussed in the previous paragraph, is coupled with the metal plasma energy of the surface where the chain or cluster is placed. This obtains a 2 × 2 matrix for each mode. Thus, the determinant becomes: where k is the resonance energy of the mode k of the chain, S is the plasma energy of the metal, and the interaction integral  is the shift in energy due to the interaction of the chain with the surface. A sketch of the level formation is shown in Figure 8 for the chain system and for the chain on the metallic surface. For large numbers of NPs (N > 10) in the chain, the interaction with the surface generates three groups of levels. The first one concerns the influence of NPs system on a surface where the energy is beyond the plasma energy of the metal. It is clearly a set of dark modes. The other two groups represent two sets of modes of the chain on the metallic surface, where the interaction of the surface with the chain decreases the energy of both low and high energy modes of the insulated chain.  [37] in the case of triangle, square, and heptamer clusters. Input data have been estimated from Reference [37] and are ω 0 = 2.15 eV and β = −0.28 eV.

Effect of the Surface
The effect of the metallic surface on the chain plasmon resonance is of great interest and has been investigated by several authors (see [20,[38][39][40][41]). Two approaches may be feasible in the context of the VM. Applying the first approach, we can directly couple the surface plasmon energy independently with all the single NPs constituting the chain. From the mathematical point of view, it means adding a further equation in the system as well as adding a second term, α, that takes into account the energy shift due to the interaction of an individual NP with the metallic surface. In any case, this approach would take into account only the main energy shift due to a strong coupling of the metal surface and the NPs. On the contrary, the surface can interact with the NPs chain, which induces two main groups of modes, with one representing the most effective coupling and the other one's surface interacting destructively with the modes of the chain [20]. In order to better represent this phenomenon, we propose coupling directly each chain wavefunction with the metallic surface. Consequently, each single mode determined, as discussed in the previous paragraph, is coupled with the metal plasma energy of the surface where the chain or cluster is placed. This obtains a 2 × 2 matrix for each mode. Thus, the determinant becomes: ( where ω k is the resonance energy of the mode k of the chain, ω S is the plasma energy of the metal, and the interaction integral α is the shift in energy due to the interaction of the chain with the surface. A sketch of the level formation is shown in Figure 8 for the chain system and for the chain on the metallic surface. 2.38 and  = −0.28 eV. By a formal point of view, it would be more correct to use the energy of the metal surface plasmon (6.4 eV) and the corresponding interaction integral, but, as from the practical point of view, the energy of the metal is only the baseline for the calculation. We decided to use the plasma frequency (8.9 eV) and the corresponding interaction integral, which is more easy to estimate from literature data.  For large numbers of NPs (N > 10) in the chain, the interaction with the surface generates three groups of levels. The first one concerns the influence of NPs system on a surface where the energy is beyond the plasma energy of the metal. It is clearly a set of dark modes. The other two groups represent two sets of modes of the chain on the metallic surface, where the interaction of the surface with the chain decreases the energy of both low and high energy modes of the insulated chain. Moreover, it is worth underlining that the dark modes of the insulated chain become bright under the effect of interaction with the surface. It clearly explains the nature of the second peak when a chain of NPs is placed on the metallic surface [39,40]. This issue is very interesting as the introduction of the surface in the investigated system allows n new coherent dipole interactions due to the interaction of every single particle constituting the chain and the surface plasmon. In this view, the number of coherent dipole interactions where hot spots can be formed becomes m s = 2n − k for the positive interaction with the surface while the number of nodes remains the same as the insulated chain. On the contrary, for the destructive interaction of the chain with the surface, the number of nodes becomes N s = n + k − 1, while the number of coherent dipole interactions remains the same as the insulated chain. This effect is shown in Figure 9 in the case of a four NPs chain.  Results are shown in Figure 10, which depicts the trend of Peak I and II, interpreted as the lowest and highest energy mode of the bright band, as a function of the chain length. It is observable that, while the group of modes of Peak II is slightly blue-shifted with respect to the free chain, Peak I is strongly red-shifted, which allows resonance until the near IR. The two peaks represent two different ways in which the metallic surface interacts with the chain: Peak I is obtained by a further stabilization of the system due to the localization of the metal plasmon with opposite charges with respect to the particle charge in contact with the surface. This stabilization of the charges of the entire system decreases the plasmonic energy. On the contrary, Peak II is the result of the stabilization of high energy modes of the chain due to the effect of the coupling of these modes with the surface plasmon of the metallic substrate. As the interaction of the dipole in the chain for these high energy modes increases, the energy of the system leads to an increasing number of elements in the chain. The resonance wavelength is slightly blue-shifted.
The same approach can be applied to other geometries, as shown in Figure 11, for the case of the heptamer cluster. As a general result, the coupling of the NPs cluster with the surface induces an In order to explore the potentiality of this method, we calculated the resonance energy of the modes related to a chain of gold NPs on a gold surface, assuming the input data as follows: ω s = 8.9 eV, α = −2.6 eV, with ω k calculated using the VM, as described in the previous paragraph with ω 0 = 2.38 and β = −0.28 eV. By a formal point of view, it would be more correct to use the energy of the metal surface plasmon (6.4 eV) and the corresponding interaction integral, but, as from the practical point of view, the energy of the metal is only the baseline for the calculation. We decided to use the plasma frequency (8.9 eV) and the corresponding interaction integral, which is more easy to estimate from literature data.
Results are shown in Figure 10, which depicts the trend of Peak I and II, interpreted as the lowest and highest energy mode of the bright band, as a function of the chain length. It is observable that, while the group of modes of Peak II is slightly blue-shifted with respect to the free chain, Peak I is strongly red-shifted, which allows resonance until the near IR. The two peaks represent two different ways in which the metallic surface interacts with the chain: Peak I is obtained by a further stabilization of the system due to the localization of the metal plasmon with opposite charges with respect to the particle charge in contact with the surface. This stabilization of the charges of the entire system decreases the plasmonic energy. On the contrary, Peak II is the result of the stabilization of high energy modes of the chain due to the effect of the coupling of these modes with the surface plasmon of the metallic substrate. As the interaction of the dipole in the chain for these high energy modes increases, the energy of the system leads to an increasing number of elements in the chain. The resonance wavelength is slightly blue-shifted.  Results are shown in Figure 10, which depicts the trend of Peak I and II, interpreted as the lowest and highest energy mode of the bright band, as a function of the chain length. It is observable that, while the group of modes of Peak II is slightly blue-shifted with respect to the free chain, Peak I is strongly red-shifted, which allows resonance until the near IR. The two peaks represent two different ways in which the metallic surface interacts with the chain: Peak I is obtained by a further stabilization of the system due to the localization of the metal plasmon with opposite charges with respect to the particle charge in contact with the surface. This stabilization of the charges of the entire system decreases the plasmonic energy. On the contrary, Peak II is the result of the stabilization of high energy modes of the chain due to the effect of the coupling of these modes with the surface plasmon of the metallic substrate. As the interaction of the dipole in the chain for these high energy modes increases, the energy of the system leads to an increasing number of elements in the chain. The resonance wavelength is slightly blue-shifted.
The same approach can be applied to other geometries, as shown in Figure 11, for the case of the heptamer cluster. As a general result, the coupling of the NPs cluster with the surface induces an The same approach can be applied to other geometries, as shown in Figure 11, for the case of the heptamer cluster. As a general result, the coupling of the NPs cluster with the surface induces an evident red-shift of the modes and, consequently, some of those levels that were dark modes become bright modes. This is a further enlargement of the resonance band and a new peak in the left side of the band. It is analogical to what has been discussed in case of the chain coupled to the metallic surface. If the heptamer is tested with the VM at the experimental conditions used for the chain system discussed previously, the maximum of the resonance peak is expected at 1380 nm. Moreover, a resonance mode at the same energy of the single NPs is formed as well as a narrow band of levels at a high energy, as already observed in the case of the chain. evident red-shift of the modes and, consequently, some of those levels that were dark modes become bright modes. This is a further enlargement of the resonance band and a new peak in the left side of the band. It is analogical to what has been discussed in case of the chain coupled to the metallic surface. If the heptamer is tested with the VM at the experimental conditions used for the chain system discussed previously, the maximum of the resonance peak is expected at 1380 nm. Moreover, a resonance mode at the same energy of the single NPs is formed as well as a narrow band of levels at a high energy, as already observed in the case of the chain.

Conclusions
In conclusion, the most important advantage of the proposed method is the rapid response and ease of application to all kinds of chain and cluster NPs systems, if the basic interactions of the elements of the system are known. This means that this approach is unable to determine the

Conclusions
In conclusion, the most important advantage of the proposed method is the rapid response and ease of application to all kinds of chain and cluster NPs systems, if the basic interactions of the elements of the system are known. This means that this approach is unable to determine the dependence of the energy shift due to the interparticle distance and particle size. However, it enables us to estimate how a given interaction changes as a result of an increasing number of conjugated NPs. This is possible without having to resort to a time-consuming calculation routine. As a consequence of a higher number of dipoles involved constructively in the state with lower energy, the corresponding wavelength represents the wavelength of the peak of the system's resonance band. On the contrary, the whole band, with this quantum chemistry approach, is the result of the convolution of all the set of modes that are obtained by solving the secular determinant. The latter observation suggests that the incident electromagnetic radiation can polarize the NPs dipoles in a wide range of wavelengths.
Lastly, the aim of the application of this methodology is to support, by means of a rapid overview of the resonance wavelength characteristics of the system, the realization of more complex models as well as of experiments at nanoscale. Moreover, the present approach allows us to investigate the effect of the geometrical distribution of the NPs system, as well as the effect of substrates in an immediate way and this can be useful for selecting the laser source in laser-based analytical techniques as well as for selecting optimal substrates for the analysis. Although, in the present study, we focused our efforts on introducing a new model for interpreting the plasmon resonance in complex NPs system. In a future work, the use of the obtained wavefunctions will be investigated for estimating other properties of the plasmonic resonance from a more applicative point of view.
Author Contributions: All the authors contributed in the same way to the present work.
Funding: This research did no receive external funding.