In this section, we extend ωFQFμ (Giovannini et al., 2022) to describe the optical properties of bimetallic systems. To this end, we follow the derivation reported in (Giovannini et al., 2019b) and (Giovannini et al., 2022). In ωFQFμ, the fluctuations of the electronic density of a plasmonic substrate composed of N atoms emerging as a response to the external electric field are represented by a set of N discrete fluctuating complex-valued charges q _i and N fluctuating complex-valued dipoles μ _i , which are located at atomic positions. The equation of motion of charges q originates from the continuity equation (Giovannini et al., 2019b) (atomic units are used throughout the paper): d q i d t = ∑ i j A i j n j 〈 p j 〉 i ⋅ l ̂ j i − n i 〈 p i 〉 j ⋅ l ̂ i j = ∑ i j A j n j 〈 p j 〉 i ⋅ l ̂ j i − A i n i 〈 p i 〉 j ⋅ l ̂ i j (1) where n _i is the electron density on atom i, ⟨ p i ⟩ j is the momentum of the electron associated with atom i averaged over all the trajectories towards atom j. l ̂ i j = − l ̂ j i is the unit vector connecting i and j. A _ij is the effective area dividing ith and jth atoms, and its value modulates the charge exchange between each atom pair. By following (Giovannini et al., 2019b; Giovannini et al., 2020; Giovannini et al., 2022), A _ij is approximated as an atomic parameter (A _i ).The electron momentum p _i can be estimated by means of the Drude model of conductance (Bade, 1957) as follows: d p i d t = E i t − p i τ i (2) where E _i (t) is the electric field acting on the ith atom and τ _i is the relaxation time associated with scattering events. Assuming the external uniform monochromatic electric field E ^ext(t) = E ^ext(ω) exp (−iωt) oscillates at frequency ω, we can reformulate Eqs 1, 2 in the frequency domain, i.e.: − i ω q i = ∑ j A j n j 〈 E j ω 〉 i ⋅ l ̂ j i 1 / τ j − i ω − A i n i 〈 E i ω 〉 j ⋅ l ̂ i j 1 / τ i − i ω (3) Following the original derivation reported in (Giovannini et al., 2019b), we can write: ⟨ E i ω ⟩ j ⋅ l ̂ i j ≈ ϕ j − ϕ i l i j (4) where ϕ _i is the chemical potential of atom i and l _ij = l _ji is the distance between ith and jth atoms. By this, Eq. 3 becomes: − i ω q i = ∑ j A j n j 1 − f j i l i j 1 / τ j − i ω + A i n i 1 − f i j l i j 1 / τ i − i ω ϕ i − ϕ j l i j (5) where quantum tunneling effects are expressed in terms of a Fermi-like function f _ij (l _ij ), which reads: f i j l i j = 1 1 + exp − d i j l i j s i j ⋅ l i j 0 − 1 (6) where l i j 0 is the equilibrium distance between atoms i and j, and the parameters d _ij and s _ij determine the shape of the damping function. Notice that the inclusion of a phenomenological description of quantum tunneling is needed to correctly describe the optical response of plasmonic subnanometer junctions and hot-spots (Esteban et al., 2012; Esteban et al., 2015; Scholl et al., 2013; Giovannini et al., 2019b; Bonatti et al., 2022; Giovannini et al., 2022). Indeed, in the case of plasmonic nanoaggregates, charge transfer plasmons (CT) may occur. In ωFQFμ, the charge exchange among different NPs is governed by Drude and tunneling mechanisms, and the total system charge is conserved (see Eq. 5), allowing for a physically consistent description of such plasmonic modes. To make the notation more compact, the following terms are introduced (Lafiosca et al., 2021): w i ω = 2 n i 1 / τ i − i ω (7) K i j = A i 1 − f i j l i j l i j (8) z i ω = − i ω w i ω (9) Dividing Eq. 5 by w _i (ω), we obtain: z i ω q i = 1 2 ∑ j K i j + w j ω w i ω K j i ϕ i − ϕ j (10) where the chemical potential ϕ _i can be written as: ϕ i = ∑ k T i k q q q k + ∑ k T i k q μ μ k + V i ext (11) The first and the second terms in Eq. 11 are the electric potential generated by the charges and the dipoles of the system, mediated by the interaction tensors T ^qq and T ^qμ (see (Giovannini et al., 2019a) for their definitions). V i ext is the electric potential associated with the external electric field. By plugging Eq. 11 into Eq. 10, the linear equations ruling charge evolution are obtained, i.e.: Z ω q = 1 2 K ̄ + H ̄ ω T q q q + 1 2 K ̄ + H ̄ ω T q μ μ + 1 2 K ̄ + H ̄ ω V (12) where the following matrices are introduced: H i j ω = w j ω w i ω K j i (13) K ̄ i j = K i j − ∑ k K i k δ i j (14) Z i j = z i δ i j (15) H ̄ i j = H i j − ∑ k H i k δ i j (16) The electric dipoles μ _i , which are introduced to properly model interband contributions, are obtained by solving the following set of linear equations: μ i = α i ω E i ext ω + E i q ω + E i μ ω (17) where α _i (ω) is the frequency-dependent polarizability of the ith atom, related to interband terms, which multiplies the total electric field acting on the ith dipole, composed by the external field E i ext , the field generated by the fluctuating charges, E i q , and the field generated by the fluctuating dipoles, E i μ . If the system is homogeneous, α _i is calculated as reported in (Giovannini et al., 2022): α i ω ≡ α A IB ω = ε A IB ω 4 π n A (18) where ε A IB is the interband contribution to the frequency-dependent permittivity function and n _A is the number density of atoms of the chemical element A. ɛ ^IB is calculated by subtracting the Drude term from the frequency-dependent permittivity.

If the system is heterogeneous, the link between the macroscopic, bulk dynamical permittivity ɛ ^IB(ω) and the microscopic, atomistic frequency-dependent interband polarizability α _i (ω) is not straightforward. Let us focus on bimetallic systems, composed of two chemical elements A and B. We assume the interband polarizability of the ith atom belonging to the A moiety to be a function of the local composition of the system, i.e., α A , i IB,alloy ( N A , i , N B , i ) , where N _A,i and N _B,i are the number of nearest neighbors of A and B type, respectively. The simplest approach is to model α A , i IB, alloy as the weighted average of α A IB and α B IB . Therefore, the interband frequency-dependent polarizability of the ith atom of species A in the alloy composed by materials A and B is defined as follows: α i ω ≡ α A , i IB,alloy ω = N A , i + 1 α A IB ω + N B , i α B IB ω N A , i + N B , i + 1 − 1 (19) Notice that, possible alternative approaches, e.g., to resort to arithmetic averages have been investigated. However, the weighted harmonic mean indeed gives the best numerical results (see Supplementary Figure S1 and Sec. S1 in the Supplementary Materials – SM).

Starting from Eq. 17, the final expression for the electric dipoles becomes (Giovannini et al., 2022): μ i = α i ω E i ext + ∑ k T i k μ q q k + ∑ k T i k μ μ μ k (20) where the field generated by the multipoles is mediated by the interaction tensors T ^μq and T ^μμ, defined according to (Giovannini et al., 2019a).

The coupled charge-dipole equations can finally be written in a compact notation as follows: 1 2 K ̄ + H ̄ ω T q q − Z ω 1 2 K ̄ + H ̄ ω T q μ T μ q T μ μ − Z IB ω q μ = 1 2 K ̄ + H ̄ ω V ext − E ext (21) where Z ^IB is a diagonal matrix, of which the elements Z i IB ( ω ) = − 1 / α i ( ω ) are defined according to Eq. 19. Notice that in the case of homogenous systems, i.e., A = B, the standard ωFQFμ equations presented by Giovannini et al. (2022) are recovered. From the values of charges and dipoles, the complex polarizability α ̄ can be computed, from which the absorption cross section σ ^abs can be calculated: α ̄ ω k l = ∂ d ̄ k ω ∂ E l ω = ∑ i q i , k ω ⋅ k i E l ω − ∑ i μ i , k ω (22) σ abs = 4 π 3 c ω t r I m α ̄ ω (23) In the previous equations d is the total complex dipole moment, i runs over NP atoms, k represents x, y, z positions of the ith atom, and l runs over x, y, z directions of the external electric field. c is the speed of light and Im ( α ̄ ) is the imaginary part of the complex polarizability α ̄ .

The developed ωFQFμ approach is challenged to reproduce the optical response of several Ag-Au nanostructures. In particular, we first consider Ag-Au spherical nanoalloys with diameter D = 5.2 nm (4,347 atoms) and Ag-Au nanorods with D = 2.5 nm and length L = 10 nm (2,560 atoms), which are generated by using Atomic Simulation Environment (ASE) Python module v. 3.17 (Larsen et al., 2017). A lattice constant of 4.08 Å (Haynes, 2014) and a Face-Centered Cubic (FCC) packing are exploited. For both spheres and nanorods, ten alloy compositions are considered by starting from pure Ag structures and increasing the percentage of Au atoms, which randomly replace Ag atoms (from 0% to 100%, with a constant step of 10%). To gain statistical significance, for each Au percentage, ten nanoalloy structures are generated by randomly replacing the proper fraction of Ag with Au, similarly to the strategy followed in previous studies (Sørensen et al., 2021). As a further example, core-shell spherical systems, which are constituted by an inner Au sphere (D = 5.0 nm, 3,851 atoms) surrounded by an outer Ag shell (D = 6.25 nm, 3,698 atoms), are investigated.

ωFQFμ equations are solved by using a stand-alone Fortran 95 package. ωFQFμ parameters (see Eqs 5, 6, 18) for Ag and Au atoms are taken from (Giovannini et al., 2022). In particular, Drude and interband parameters are recovered from Ag and Au experimental permittivity functions (Etchegoin et al., 2006). The parameters d _ij = 12.0 and s _ij = 0.95 entering Eq. 6 are obtained by fitting ωFQFμ results for i = Ag and j = Au (and vice versa) on reference time-dependent density functional theory + tight binding (TD-DFT + TB) (Asadi-Aghbolaghi et al., 2020) absorption spectra of Ag-Au clusters (see Sec. S2 in the SM).

The modeling of the dependence of the absorption cross-section σ ^abs of a spherical NP on the Au concentration has received much attention in the literature (Papavassiliou, 1976; Link et al., 1999; Rioux et al., 2014; Ristig et al., 2015). In particular, a simple linear combination of Ag and Au permittivity functions combined with the Mie theory cannot reproduce the linear behavior of the plasmon red-shift upon increasing Au concentration which is experimentally observed (Link et al., 1999). This issue can be solved by following the strategy proposed by Rioux et al. (2014), who have developed an analytical and rigorous model to predict a composition-dependent complex dielectric function of Ag-Au alloys, based on critical point analysis of the band-structure of Ag and Au. By coupling such modeling to the Mie theory, the correct trend is followed for spherical nanoalloys. It is worth noting that such a strategy might be exploited within the more refined, yet continuum, BEM method, thus overcoming the limitations of this approach in describing the plasmonic properties of nanoalloys. However, while the methodology proposed by Rioux et al. (2014) is rigorous from the theoretical point of view, its extension to different alloy compositions, and to alloys of different chemical nature rather than Au and Ag is far from trivial, due to the high (about 30) number of parameters that would need to be fitted to obtain the composition dependent dielectric function.

Computed ωFQFμ values for σ ^abs of a spherical NP with different Au concentrations are reported in (see Figure 1). In all cases, spectra are characterized by a main peak (associated with the dipolar mode, see Supplementary Figure S2 in the SM), which redshifts and lowers in intensity as the Au concentration increases. By plotting the PRF (in nm) as a function of the percentage of Au (see Figure 1C), a linear trend is observed (R ² ∼ 1.00). This behavior perfectly follows Vegard’s law, which reads: λ V egard x = 1 − x λ A g + x λ A u (24)

where x is the percentage of Au atoms. Such a linear dependence of the absorption wavelength on Au percentage is reported by many experimental works (Papavassiliou, 1976; Link et al., 1999; Rioux et al., 2014; Ristig et al., 2015). The slope (2.06 nm/% Au) of the line fitted on ωFQFμ results is in good agreement with experimental data (Link et al., 1999) (∼ 1.35 nm/% Au). Furthermore, by moving from pure Ag to pure Au NPs, the intensity of the plasmon band exponentially decreases (R ² = 0.99, see Figure 1D), thus perfectly reproducing experimental findings (Link et al., 1999).

We note that, remarkably, other classical atomistic approaches may fail at reproducing the experimentally reported linear trend in wavelength and the exponentially decreasing intensity (Sørensen et al., 2021). Differently, ωFQFμ can correctly reproduce all the plasmonic properties of these structures, which is a key result of this paper.

We now move to discuss the plasmonic properties of Ag-Au alloy nanorods. Also in this case, we study the dependence of the absorption frequency on the Au concentration, which varies from 0% (pure Ag nanorod) to 100% (pure Au nanorod), with a constant step of 10%. ωFQFμ absorption cross sections, which are reported in Figure 2B, are clearly dominated by an intense peak which redshifts and decreases in intensity as the gold concentration increases. This peak is associated with the longitudinal dipolar mode (see also Supplementary Figure S3 in the SM). From a comparison with spherical NPs (see Figure 1B), we note that the dipolar peak moves to lower energies, independently of the chemical composition of the alloy nanorod. Also, in this case, additional bands can be appreciated in the region between 350 and 500 nm. There, octupolar plasmons arise (Bonatti et al., 2020), together with transversal dipolar modes. The associated bands associated are graphically depicted as an inset in Figure 2B, and redshift and lower in intensity as the Au concentration increases.

The PRFs (in nm) of both the longitudinal and transversal dipolar modes as a function of the Au percentage are graphically depicted in Figure 2C. The transversal mode follows the linear trend (circles, R ² = 0.99) predicted by Vegard’s law (see Eq. 24), while for the longitudinal plasmon, an evident deviation from the linear regime is observed, especially for Au concentrations larger than 40%. Our findings are in agreement with the experimental measurements reported by Bok et al. (2009), where for the transversal peak a linear trend is measured. Also, the experimental slope (1.47 nm/%Au) correlates significantly well with our calculations (1.92 nm/%Au). Differently, the longitudinal absorption wavelength is reported to be rather independent of the chemical composition for Au percentages > 40%, while a slight blueshift is expected by increasing the Ag concentration Bok et al. (2009). Note that a DDA description is able to reproduce the composition dependence of the longitudinal mode but fails at describing the linear trend of the transversal peak (Bok et al., 2009). This is probably related to the exploited linear combination of Ag and Au permittivity functions. Remarkably, ωFQFμ correctly models both behaviors.

These findings can be explained by considering that the longitudinal dipolar mode falls at energies that are reasonably far from interband transition regions of both Ag and Au. In this regime, Drude free electrons play a dominant role (Link et al., 1999; Bok et al., 2009) and the good agreement with experiments reported by both DDA (Bok et al., 2009) and ωFQFμ show that the permittivity function in this region can be reasonably approximated as a linear combination of Ag and Au Drude contributions. Differently, the transversal plasmon mode absorbs at much higher energies, where interband effects cannot be neglected (Kreibig et al., 1995). This proves the reliability of ωFQFμ, which, differently from other models (Draine and Flatau, 1994), describes the plasmonic resonance by means of both electric charges and dipoles, modeling the intraband and interband mechanisms, respectively.

ωFQFμ is not limited to the description of nanoalloys, but can generally describe bimetallic nanostructures. To showcase its applicability, we consider pure core-shell (Au@Ag-0%) spherical NPs with a diameter D of 6.25 nm (3,698 atoms, see Figure 3A), constituted by a Au core (D _Au = 5 nm), and an external Ag shell (external D _Ag = 6.25 nm, D _Au /D _Ag = 0.8). In addition, to show the flexibility of ωFQFμ, the effect of alloying such structures, which has been experimentally realized (Ristig et al., 2015; Blommaerts et al., 2019), is also taken into account, by randomly replacing 18% of the core Au atoms with Ag, and viceversa for the Ag shell (Au@Ag-18%, see Figure 3A).

Computed ωFQFμ absorption cross sections in the visible range for both core-shell systems are reported in Figure 3B, together with reference spectra for Ag and Au spherical NPs with D = 6.25 nm. Similarly to Figure 1B, Ag and Au reference absorption spectra are characterized by a main peak, located at about 360 nm and 570 nm for Ag and Au, respectively. The spectrum of the pure Au@Ag-0% core-shell is instead dominated by two peaks A and B, in agreement with (Peña-Rodríguez and Pal, 2011a; Chen et al., 2012; Gao et al., 2014), which fall at about 356 nm (A) and 517 nm (B). To deeply investigate the nature of the plasmons associated with these two bands, we graphically report in Figure 3C, left panel, the electron densities at the excitation energies for pure Au@Ag-0% core-shell (on a slice in the xy plane). As it can be appreciated, both peaks are associated with a global dipolar plasmon. However, they differ for the electron distribution in both the Au core and the interface between the two regions (see yellow dashed line in Figure 3C). Indeed, for peak A a huge charge accumulation at the interface is reported, while for peak B, the boundary between the Ag shell and the Au core is only partially marked. This is due to the fact that our model allows for charge exchange between the two layers, which is for instance not accounted for by other classical descriptions (Szántó et al., 2021). However, such a charge transfer is limited by the different tunneling barriers between the two metals, as introduced in Eq. 6, and by the different chemical nature (chemical hardness, polarizability, Drude parameters, …). As a consequence, charge accumulation at the Ag-Au interface is expected.

Similarly, the Au@Ag-18% absorption cross section (see Figure 3B) is characterized by two main peaks A′ and B′, which are blue- and red-shifted as compared to their counterparts in Au@Ag-0% system. Such a behavior can be justified by considering that by increasing the alloying degree up to 50%, a fully alloyed Ag-Au spherical NP is obtained. Therefore, in this situation, a single band is expected at about 450 nm (see also Figure 1B). The plasmonic nature of the electron densities computed at the A′ and B′ frequencies is graphically reported in Figure 3C, right panel. A global dipolar plasmon is observed. However, by also comparing with Au@Ag-0%, charge accumulation at the core-shell boundary is still present, but less pronounced, in particular for B′. This confirms the above speculation: indeed, by increasing the alloying degree, a decrease in the potential barrier at the interface is obtained.

To conclude the discussion on Au@Ag core-shell systems, we compare our results with experimental absorption cross sections reported by Peña-Rodríguez and Pal (2011b) (see inset in Figure 3B), who studied the plasmonic response of Au@Ag core-shell spherical NPs with D = 44.2 nm and characterized by D _Au /D _Ag = 0.81 (comparable with our simulations). As it can be noticed, the experimental spectrum is characterized by two main peaks at about 400 nm and 497 nm. This is in very good agreement with our calculations, thus further demonstrating the reliability of our model. In fact, only a slight discrepancy in relative intensities, probably due to different experimental conditions (size of the NP, solvent, …), is observed.