Nanoscale Organic Contaminant Detection at the Surface Using Nonlinear Bond Model

Abstract

Environmental pollution from organic dyes such as malachite green and rhodamine B poses significant threats to ecosystems and human health due to their toxic properties. The rapid detection of these contaminants with high sensitivity and selectivity is crucial and can be effectively achieved using nonlinear optical methods. In this study, we combine the Simplified Bond Hyperpolarizability Model (SBHM) and molecular docking (MD) simulations to investigate the Second-Harmonic Generation (SHG) intensity of organic dyes on a silicon (Si(001)) substrate for nanoscale pollutant detection. Our simulations show good agreement with rotational anisotropy (RA) SHG intensity experimental data across all polarization angles, with a total error estimate of 3%. We find for the first time that the SBHM not only identifies the different organic pollutant dyes on the surface, as in conventional SHG detection, but can also determine their relative orientation and different concentrations on the surface. Meanwhile, MD simulations reveal that rhodamine B shows a strong adsorption affinity of $-10.4\mathrm{kcal}/\mathrm{mol}$ to a single-layer graphene oxide (GO) substrate, primarily through $\pi$-$\pi$ stacking interactions (36 instances) and by adopting a perpendicular molecular orientation. These characteristics significantly enhance SHG sensitivity. A nonlinear susceptibility analysis reveals good agreement between the SBHM and group theory. The susceptibility tensors confirm that the dominant contributions to the SHG signal arise from both the molecular structure and the surface interactions. This underscores the potential of GO-coated silicon substrates for detecting trace levels of organic pollutants with interaction distances ranging from $3.75\AA$ to $5.81\AA$. This approach offers valuable applications in environmental monitoring, combining the sensitivity of SHG with the adsorption properties of GO for nanoscale detection.

1. Introduction

The rapid growth of industrial activities has significantly increased environmental pollution, particularly in aquatic systems. Synthetic dyes, extensively used in the textile, agricultural, and cosmetics industries, are among the most common contaminants affecting water resources. Dyes such as malachite green (MG) and rhodamine B exemplify hazardous pollutants that persist in ecosystems, posing severe risks to both human health and environmental quality [1,2]. MG, a triphenylmethane dye widely utilized in textiles, in aquaculture, and as an antibacterial agent, is known for its toxic, carcinogenic, and mutagenic effects [1]. Similarly, rhodamine B, a xanthene dye commonly employed in biotechnology and materials science, exhibits adverse ecological and health impacts when present as a contaminant [3,4].

To address these environmental challenges, the development of detection methods that are rapid, highly sensitive, and selective has become crucial. Nonlinear optical techniques, such as Second-Harmonic Generation (SHG), Frequency Mixing, and Surface-Enhanced Raman Scattering (SERS), have emerged as promising approaches for identifying trace contaminants by leveraging unique interactions between light and matter [5]. These techniques enable the detection of extremely low concentrations of hazardous dyes with high selectivity, as demonstrated in previous studies [6,7,8]. SHG, in particular, has been applied to monitor organic pollutants at interfaces, including in studies that utilized rotational anisotropy SHG (RASHG) to characterize MG on silica surfaces and other solid–liquid interfaces [9,10]. An experimental RASHG setup used for such analyses is illustrated in Figure 1, showcasing its capability to detect contaminants through SHG intensity profiles.

In this study, we apply for the first time the nonlinear Simplified Bond Hyperpolarizability Model (SBHM) combined with molecular docking simulations to analyze the SHG intensity profiles of MG and rhodamine B adsorbed on silicon (Si(001)) substrates and perform an error estimation. The SBHM framework models the molecular structures using bond vectors and susceptibility tensors to predict SHG behavior at the nanoscale where the nonlinear signal is assumed to arise from the anharmonic bond oscillation along the covalent bonds [6,8]. The study of Second-Harmonic Generation (SHG) has a rich historical background rooted in phenomenological bond models for nonlinear optical effects. Foundational work by various researchers, such as Bloembergen et al. [12], Guidotti et al. [13], Sipe et al. [14], and Luepke [15], introduced a comprehensive theoretical framework to describe SHG from semiconductor surfaces. This model established the connection between nonlinear optical susceptibility and the intrinsic properties of surface atoms, providing critical insights into the interaction of light with reconstructed and pristine surfaces. Following this, Powell et al. [8] implemented the Surface Bond Hyperpolarizability Model (SBHM), a significant advancement that enabled the quantitative analysis of surface-specific SHG signals.

The tensorial analysis of silicon semiconductors underwent further refinement through group theory applications, as demonstrated in the work by Adalberto et al. [6]. By contrasting group theory with the SBHM, this study underscored the advantages of the SBHM in capturing surface-induced symmetry-breaking effects, making it an indispensable tool for analyzing SHG signals. More recently, the SBHM framework has been extended to model surface reconstruction phenomena, including twin boundaries, as presented by Hardhienata et al. [16]. This extension allowed for the exploration of intricate surface interactions and their contributions to nonlinear optical responses. Despite the significant progress in modeling SHG from semiconductor surfaces, its application to pollutant detection remains undeveloped. Early attempts to model pollutants on surfaces provided a qualitative framework but lacked precise methodologies for detecting pollutant concentrations. Our work builds upon this foundation by introducing a novel computational approach that incorporates pi–pi stacking interactions and nonlinear optical susceptibility enhancements. Moreover, we provide quantitative calculations that directly correlate SHG intensity with pollutant concentration by comparing the raw RASHG experimental data from Gassin et al. [9] and Morgenthaler et al. [17], addressing a critical gap in the field.

However, we would also like to mention several limitations that exist in this work. These include the assumptions of idealized molecular orientations and that the anharmonic oscillation only originates along the covalent bond vectors and not transversally. Morevover, the model also excludes solvent effects and temperature-dependent variations, and it lacks ab initio quantum mechanical calculations for a comprehensive analysis of nonlinear susceptibility. Future studies should address these limitations by incorporating environmental factors, advanced computational techniques, and a broader range of molecular systems to improve the practical applicability of SHG analysis.

Additionally, molecular docking simulations investigate the adsorption mechanisms of rhodamine B onto graphene oxide (GO) layers, revealing strong chemisorption dominated by $\pi$-$\pi$ stacking and $\pi$-alkyl interactions, with a binding affinity of $-10.4$ kcal/mol. These findings suggest that GO-coated silicon substrates offer a promising platform for detecting organic pollutants with high sensitivity, presenting an innovative approach for environmental monitoring and remediation [18,19]. Thus, a combined SBHM-MD simulation differs from conventional pollutant detection methods, which lack the ability to determine molecular orientation at the nanoscale. Unlike other nonlinear optical techniques such as Surface-Enhanced Raman Scattering (SERS), SHG uniquely benefits from its sensitivity to surface symmetry-breaking effects and offers a strong physical picture of how nonlinear radiation is generated. Furthermore, while conventional pollutant detection methods often require direct chemical analyses or extensive sample preparation, SHG offers a non-invasive, real-time detection capability with the potential for higher selectivity and sensitivity, making it particularly well suited for detecting organic dyes in environmental applications [20,21].

2. Methods

In this section, we outline the theoretical framework and simulation approaches employed to study nonlinear optical phenomena. Polarization behavior was modeled using second-order susceptibility tensors, and the bond vector model was applied to represent dipole interactions along atomic bonds. Additionally, molecular docking simulations were conducted to analyze adsorption behavior, leveraging AutoDock Vina with optimized parameters for accurate binding energy predictions.

2.1. Nonlinear Optics Approach

Nonlinear optics explores the phenomena that arise when the optical properties of a material are altered under the influence of intense light, typically from lasers. In such cases, an external high-intensity field induces a dipole moment in the material, a phenomenon termed polarization. The total polarization $P_i$ can be expressed as [7]

$${\mathbf{P}}_i={\epsilon }_0\left({\chi }_{ij}^{\left(1\right)}{E}_j+{\chi }_{ijk}^{\left(2\right)}{E}_j{E}_k+{\chi }_{ijkl}^{\left(3\right)}{E}_j{E}_k{E}_l+\cdots \right),$$

(1)

where ${\epsilon }_0$ is the vacuum permittivity, ${\chi }^{\left(1\right)}$ is the linear susceptibility, and ${\chi }^{\left(2\right)}$ and ${\chi }^{\left(3\right)}$ are the second- and third-order nonlinear susceptibilities, respectively. The electric fields $E_j,E_k,$ and $E_l$ represent the applied fields. In this study, Cartesian coordinates (x, y, and z) are chosen to align with the material’s symmetry.

Second-order polarization, associated with SHG, is influenced by the second-order susceptibility tensor ${\chi }_{ijk}^{\left(2\right)}$, which contains 27 components, but symmetry considerations significantly reduce this number. For instance, under Kleinman symmetry and in molecules with $C_{2v}$ symmetry, only specific components like ${\chi }_{xyz}$ and ${\chi }_{xzy}$ remain non-zero.

The simulation applies a coordinate system, illustrated in Figure 2, being similar to the experiment in [9]. The incoming fundamental and outgoing SHG fields are situated at the $xz$-plane, and their directions are given by the vectors ${\mathbf{k}}_{\omega }$-in and ${\mathbf{k}}_{2\omega }$-out. To model experimental data, the incomming local electric field polarization unit vector direction in this work is arbitrary, defined by

$${\widehat{\mathbf{E}}}_{\mathrm{loc}}=\left(\begin{array}{c}cos{\theta }_{i}sin\psi \\ cos\psi \\ sin{\theta }_{i}sin\psi \end{array}\right),$$

(2)

where ${\theta }_i$ is the angle of incidence as depicted in Figure 2, and $\psi$ represents the polarization angle, analyzed in increments of 10°.

2.2. Bond Vector Model

Surface second-order polarization is modeled based on the SBHM developed by Powell [8] and Aspnes [22]. It assumes that nonlinear radiation is generated by anharmonic oscillations of dipoles along atomic bonds. The second-order polarization is given by

$$\begin{array}{cc}\hfill {\mathbf{P}}^{\left(2\right)}& ={\displaystyle \frac{1}{V}}\sum _i{\alpha }_i\left({\mathbf{R}}^{\left(z\right)}\left(\phi \right)·{\widehat{\mathbf{b}}}_i\right)\otimes \left({\mathbf{R}}^{\left(z\right)}\left(\phi \right)·{\widehat{\mathbf{b}}}_i\right)\otimes \left({\mathbf{R}}^{\left(z\right)}\left(\phi \right)·{\widehat{\mathbf{b}}}_i\right)··{\mathbf{E}}_{\mathrm{in}}\otimes {\mathbf{E}}_{\mathrm{in}}\hfill \\ \hfill & ={\chi }^{\left(2\right)}{\mathbf{E}}_{\mathrm{in}}{\mathbf{E}}_{\mathrm{in}},\hfill \end{array}$$

(3)

and the rotation matrix along the z-axis is

$${\mathbf{R}}^{\left(z\right)}\left(\phi \right)=\left[\begin{array}{ccc}cos\phi & -sin\phi & 0\\ sin\phi & cos\phi & 0\\ 0& 0& 1\end{array}\right].$$

(4)

where $P_D^2$ is the polarization occurring in the surface section; $E_{\mathrm{in}}$ is the input electric field; and ${\alpha }_i$ is the hyperpolarizability of the second-order nonlinear polarization, where subscript i is an integer number related to the individual bond unit vectors ${\widehat{\mathbf{b}}}_i$. Here, ⨂ denotes the outer product operation applied to the unit bond vectors of each molecule. $R^{\left(z\right)}\left(\phi \right)$ is the rotation matrix about the z-axis. The form of the rotation matrix $R^{\left(z\right)}\left(\phi \right)$ defines the rotation matrix about the z-axis, which can be written in the form of Equation (4) [7]. The bond vector ${\widehat{\mathbf{b}}}_i$ indicates the dipole orientation of the individual bonds between silicon, MG, and rhodamine B for different simulations. Figure 3 illustrates the structure of Si(001), along with its four bonding vectors ${\widehat{b}}_{1-4}$. Each of these bonding vectors has a representation in Cartesian coordinates, as shown in Equation (5).

$$\begin{array}{ccc}\hfill {\widehat{\mathbf{b}}}_1& =\left(\begin{array}{c}-sin(\beta /2)/\sqrt{2}\\ -sin(\beta /2)/\sqrt{2}\\ -cos(\beta /2)\end{array}\right),{\widehat{\mathbf{b}}}_2\hfill & \hfill =\left(\begin{array}{c}-sin(\beta /2)/\sqrt{2}\\ sin(\beta /2)/\sqrt{2}\\ -cos(\beta /2)\end{array}\right),\\ \hfill {\widehat{\mathbf{b}}}_3& =\left(\begin{array}{c}-sin(\beta /2)/\sqrt{2}\\ sin(\beta /2)/\sqrt{2}\\ cos(\beta /2)\end{array}\right),{\widehat{\mathbf{b}}}_4\hfill & \hfill =\left(\begin{array}{c}sin(\beta /2)/\sqrt{2}\\ -sin(\beta /2)/\sqrt{2}\\ cos(\beta /2)\end{array}\right)\end{array}$$

(5)

where $\beta =109.47$°. Furthermore, Figure 4 presents the structure of MG, with its bonding vectors represented as shown in Equation (6), and the bond vectors are denoted by ${\widehat{b}}_{5-7}$ [23].

$$\begin{array}{cccc}\hfill {\widehat{\mathbf{b}}}_5& =\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right),{\widehat{\mathbf{b}}}_6\hfill & \hfill =\left(\begin{array}{c}cos\left(\gamma \right)\\ 0\\ -sin\left(\gamma \right)\end{array}\right),{\widehat{\mathbf{b}}}_7& =\left(\begin{array}{c}-cos\left(\gamma \right)\\ 0\\ -sin\left(\gamma \right)\end{array}\right)\hfill \end{array}$$

(6)

The structure of rhodamine B is illustrated in Figure 5. The bond vectors represented in Cartesian coordinates are detailed by Equation (7), with the bond vector for rhodamine B denoted as ${\widehat{b}}_{8-14}$. In bond vectors, there exists a $\gamma =30$° representing the angle between the vector ${\widehat{b}}_{6-7}$ and ${\widehat{b}}_{9-12}$ on the x-axis. Under the conditions of a sufficiently high adsorbate concentration, the molecular orientation is predominantly perpendicular, as shown in Figure 4 and Figure 5. This situation occurs because, in the layer closest to the silica substrate, the adsorbate molecules adopt a parallel orientation, while in the subsequent layers, the orientation may be tilted or perpendicular [10,24,25,26].

$$\begin{array}{cccc}\hfill {\widehat{\mathbf{b}}}_8& =\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right),{\widehat{\mathbf{b}}}_9\hfill & \hfill =\left(\begin{array}{c}{\displaystyle \frac{cos\left(\gamma \right)}{\sqrt{3}}}\\ 0\\ {\displaystyle \frac{sin\left(\gamma \right)}{\sqrt{3}}}\end{array}\right),{\widehat{\mathbf{b}}}_{10}& =\left(\begin{array}{c}-{\displaystyle \frac{cos\left(\gamma \right)}{\sqrt{3}}}\\ 0\\ {\displaystyle \frac{sin\left(\gamma \right)}{\sqrt{3}}}\end{array}\right),\hfill \\ \hfill {\widehat{\mathbf{b}}}_{11}& =\left(\begin{array}{c}cos\left(\gamma \right)\\ 0\\ -sin\left(\gamma \right)\end{array}\right),{\widehat{\mathbf{b}}}_{12}\hfill & \hfill =\left(\begin{array}{c}-cos\left(\gamma \right)\\ 0\\ -sin\left(\gamma \right)\end{array}\right)\end{array}$$

(7)

2.3. Molecular Docking Simulation

Molecular docking simulations were conducted using AutoDock Vina [27] on a MacBook with a 1.6 GHz Intel Core i5 processor, 4 GB RAM, and Intel HD Graphics 6000. Chimera software [28] was employed for preparation, visualization, and analysis. The structures of GO (PubChem CID: 124202900) and rhodamine B (PubChem CID: 6694) were obtained in SDF format from PubChem and converted to PDB format using Chimera. Polar hydrogens and Gasteiger charges were added to the ligands, and blind docking was performed to evaluate adsorption on GO, with the number of modes set to 10 and exhaustiveness to 8, which were selected based on a balance between computational efficiency and the accuracy of binding energy predictions. These parameters align with standard practices for ligand–receptor docking simulations and are informed by iterative tests to ensure convergence. Vina employs an Iterated Local Search global optimizer ([27,29]), which utilizes stochastic techniques combined with local refinement to explore the conformational landscape effectively. The exhaustiveness parameter controls the number of sampling iterations, while the number of modes defines the diversity of poses returned. These values were optimized to capture the most relevant binding conformations without significant computational overhead, as highlighted in the work of Trott and Olson ([27]).

Environmental variations, such as pH and competing ions, can significantly influence adsorption energy and binding orientation. However, AutoDock Vina, by default, does not explicitly account for changes in pH or ionic strength. In our study, we addressed this limitation by including pre-docking protonation state calculations for the ligands and receptors, ensuring physiologically relevant binding conformations at specific pH levels. This approach is informed by the study of scoring function derivatives, where the scoring gradients are directly tied to the forces and torque acting on the ligand, which is based on the work in Ref. [27]. Competing ions were considered implicitly through adjustments to the electrostatic grid maps, ensuring realistic representations of the local environment. Future enhancements, such as explicit ion inclusion, could provide more accurate insights into environmental effects.

The Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton method employed by Vina further ensures efficient local optimization by leveraging gradients derived from the scoring function to improve docking accuracy following Ref. [27]. This optimization process allows Vina to efficiently explore various binding orientations while minimizing computational costs. Despite these considerations, direct simulations of environmental pH and ion competition remained computationally intensive and were beyond the scope of this study. However, our parameterization and methodology align with widely accepted docking protocols, ensuring robust and reproducible results.

3. Results and Discussion

MG and rhodamine B are widely used as textile dyes due to their strong optical properties, particularly their ability to absorb and emit light efficiently. These molecules are known to exhibit characteristics of SHG and fluorescence, with rhodamine B being especially prominent in fluorescence applications because of its high efficiency in light emission upon excitation. As in previous work involving the SBHM such as Ref. [16], we apply $C_{2v}$ as an effective symmetry for the pollutant point group. Although neither malachite green nor rhodamine B possesses strict $C_{2v}$ symmetry when considering their full molecular structures, their electronic transitions and optical responses often exhibit characteristics consistent with $C_{2v}$-like symmetry, especially when focusing on specific functional groups or molecular frameworks relevant to the studied properties. Furthermore, their planar structures and extensive electron delocalization enable strong absorption in the visible light spectrum, allowing effective interactions with specific wavelengths. These dyes are also cationic, bearing a positive charge, which facilitates electrostatic interactions with negatively charged surfaces such as silica. These unique characteristics make them ideal for studying nanoscale interactions at interfaces [17,23,24,26].

3.1. Susceptibility Tensor and Group Theory

Silicon with a Si(001) orientation belongs to the $T_d$ point group symmetry [30]. The components of its susceptibility tensor, derived from group theory and the SBHM [6], are represented in Equations (8) and (9):

$${\chi }_{T_d}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& d_{132}\\ 0& d_{132}& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& d_{132}\\ 0& 0& 0\\ d_{132}& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& d_{132}& 0\\ d_{132}& 0& 0\\ 0& 0& 0\end{array}\right)\end{array}\right).$$

(8)

$${\chi }_{Si\left(001\right)}={\alpha }_{Si}\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -S\\ 0& -S& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& -S\\ 0& 0& 0\\ -S& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& -S& 0\\ -S& 0& 0\\ 0& 0& 0\end{array}\right)\end{array}\right).$$

(9)

Here, the non-zero tensor components are $d_{132}$, $d_{123}$, $d_{231}$, $d_{213}$, $d_{321}$, and $d_{312}$, all of which are equal, as $d_{132}=-S$, where $S={\alpha }_{Si}sin(\beta /2)sin\left(\beta \right)$. The SBHM approach confirms these components, matching the predictions of group theory.

The molecular structures of MG and rhodamine B exhibit $C_{2v}$ symmetry. This symmetry results in non-zero tensor components, as shown in Equation (10) and as previously discussed in Ref. [23].

$${\chi }_{C_{2v}}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& d_{131}\\ 0& 0& 0\\ d_{131}& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& d_{232}\\ 0& d_{232}& 0\end{array}\right)\\ \left(\begin{array}{ccc}{d}_{311}& 0& 0\\ 0& d_{322}& 0\\ 0& 0& d_{333}\end{array}\right)\end{array}\right).$$

(10)

According to Refs. [9,10], the independent, non-vanishing quadratic susceptibility tensor components $d_{333}$, $d_{311}$, and $d_{113}$ are critical for analyzing isotropic and achiral interfaces. The subscripts in these tensors represent the Cartesian coordinates of the laboratory reference frame.

3.2. SBHM Simulation

3.2.1. MG/Si(001)

MG belongs to the $C_{2v}$ point group symmetry [31], with its non-zero susceptibility tensor components derived from group theory, as shown in Equation (10). Using the SBHM approach, the susceptibility tensor of MG can be represented as

$${\chi }_{MG}={\alpha }_{MG}\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& -2B\\ 0& 0& 0\\ -2B& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}-2B& 0& 0\\ 0& 0& 0\\ 0& 0& 1-2{sin}^3\left(\gamma \right)\end{array}\right)\end{array}\right),$$

(11)

where $B={cos}^2\left(\gamma \right)sin\left(\gamma \right)$, and ${\alpha }_{MG}$ represents the hyperpolarizability of MG, which is critical for interactions with an electric field. The relationship between the tensor components derived using group theory and the SBHM is expressed as

$$d_{131}=d_{311}\to -2B{\alpha }_{MG},d_{333}\to {\alpha }_{MG}\left(1-2{sin}^3\gamma \right),d_{232}=0.$$

(12)

The tensor component $d_{232}$ is zero because the structure of malachite green tends to exhibit a planar configuration, which causes the nonlinear optical properties to be more dominant along axes aligned with the molecular plane. As shown in Figure 4, the charge distribution is more pronounced along the x- and z-axes. These non-zero tensor components are consistent with those reported in [9,10], where it is stated that the tensor components contributing to the SHG signal in malachite green are $d_{113}$, $d_{333}$, and $d_{311}$. Additionally, as demonstrated in Equations (6) and (7), the y-axis components in each bond vector are zero. Consequently, there is no contribution from the susceptibility tensor along the y-axis, implying that tensor components containing the “y” index are also zero. The susceptibility tensor in Equation (11) underpins the SBHM simulation of the MG/Si(001) system, with the resulting SHG intensity presented in Figure 6. According to Ahyad [23], simulations of MG on Si(001) or Si(111) surfaces revealed that silicon’s contribution to the SHG intensity is minimal, with the SHG signal dominated primarily by MG.

Figure 6a,c depict the simulated SHG intensity results, while Figure 6b,d display the experimental results [9]. A comparison between the SBHM simulation results and experimental data demonstrates a strong agreement in terms of both magnitude and intensity profile. Variations in the concentration of MG (0.0–0.3 mM) molecules adsorbed on the silica substrate significantly influence the magnitude and profile of the SHG signal. Under conditions where MG (0.0 mM) is absent on the silica surface, the SHG intensity is entirely attributed to the contribution from silicon. However, with an increasing MG concentration, the intensity profile contributing to the SHG signal is no longer dominated solely by silicon but also includes contributions from MG molecules. The addition of MG to the silica substrate results in an enhancement of the SHG intensity, as observed in the polarization-$p_{out}$ configuration. In this configuration, changes in the MG concentration clearly affect the variation in the SHG intensity. An interesting new result is that the SBHM can be applied to investigate the varying pollutant concentration by modifying the complex SHG hyperpolarizability parameter. This is because hyperpolarizability depends on the amount of generated SHG dipoles, which correlates directly with the pollutant cocentration. Moreover, the imaginary hyperpolarizability takes into account the phase difference due to the change in the optical depth. We can express the total SHG susceptibility and polarization as

$${\overleftrightarrow{\chi }}_{SHG}^{\left(2\right)}={\overleftrightarrow{\chi }}_{Si}^{\left(2\right)}+{\overleftrightarrow{\chi }}_{MG}^{\left(2\right)}\to {\mathbf{P}}_{SHG}^{\left(2\right)}={\mathbf{P}}_{Si}^{\left(2\right)}+{\mathbf{P}}_{MG}^{\left(2\right)}$$

(13)

Moreover, as the concentration of MG increases, the hyperpolarizability at $b_5$ increases, while the hyperpolarizability at $b_6$ and $b_7$ remains unchanged. This phenomenon occurs because, with the increase in the MG concentration, the contribution of neighboring MG dipoles is influenced by the external electric field. The part of the MG molecule represented by vector $b_5$ is the first to be affected by this field, whereas $b_6$ and $b_7$ are shielded. The respective hyperpolarizability values can be expressed as follows

$${\alpha }_{MG}\left(5\right)={\alpha }_u(1+N_{MG}),{\alpha }_{MG}\left(6\right)={\alpha }_d,{\alpha }_{MG}\left(7\right)={\alpha }_d$$

(14)

$$N_{MG}\propto {\displaystyle \frac{\left[MG\right]}{C_W{e}^{\delta G/RT}+\left[MG\right]}}$$

(15)

where $N_{MG}$ represents the amount of adsorbed MG, $\Delta G$ denotes the change in Gibbs free energy, and T and $\left[MG\right]$ refer to the temperature and MG concentration, respectively.

Conversely, the hyperpolarizability of silicon decreases due to a diminished response to the electric field, which is obstructed by the presence of MG molecules. Furthermore, the SHG intensity in the polarization p-out configuration exhibits a higher value than in the s-out configuration. This indicates that the primary contribution to the SHG signal originates from molecules oriented perpendicular to the substrate. This phenomenon is illustrated in Figure 4, where the ${\widehat{\mathbf{b}}}_{\mathbf{5}}$ bond vector is shown to be perpendicular to the x-axis [9,10,23].

In our analysis, we discovered that fitting the experimental data required the assumption that a higher concentration of the pollutant corresponds to a higher dipole moment, thereby increasing the imaginary hyperpolarizability, as depicted in

Table 1. Hyperpolarisability of MG at different concentrations.

Concentration MG (mM)	${\mathit{\alpha }}_5$	${\mathit{\alpha }}_6,{\mathit{\alpha }}_7$
0.0	5	5
0.001	6.775	5
0.0025	8.425	5
0.004	9.461	5
0.01	11.396	5
0.02	12.478	5
0.03	12.924	5

. Furthermore, we determined that the hyperpolarizability must be treated as a complex quantity. This complexity arises due to changes in the optical penetration depth, which introduce a phase factor that is incorporated into the imaginary part of the SHG hyperpolarizability.

The MG SHG imaginary hypepolarizabilies (${\alpha }_5$, ${\alpha }_6$, and ${\alpha }_7$) correspond to the bond vector orientation in Figure 4. Their values should be treated as arbitrary units since they are obtained by fitting and not determined by ab initio methods, which currently lies beyond the scope of this work. The values presented in Table 1 illustrate the dependence of the hyperpolarizability (${\alpha }_5$) of malachite green (MG) on its concentration. The data show a clear trend, with ${\alpha }_5$ increasing with rising concentrations of MG. At a concentration of $0.0$ MG, ${\alpha }_5$ starts with a baseline value of 5. As the concentration increases to $0.001$ MG, $0.0025$ MG, and beyond, ${\alpha }_5$ increases steadily, reaching a maximum value of $12.9249$ at $0.03$ MG. This trend suggests that higher concentrations of MG enhance the dipole–dipole interactions in the system, leading to greater induced polarization. In contrast, the values of ${\alpha }_6$ and ${\alpha }_7$ remain constant at 5 across all concentrations, indicating that the concentration effect on hyperpolarizability is dominant specifically in the ${\alpha }_5$ component. This can be explained because the ${\alpha }_5$ hyperpolarizability contributes only to a DC-like increase in the SHG intensity, whereas changes in the other hyperpolarizability parameter will change the SHG intensity peak features. The increasing ${\alpha }_5$ reflects the system’s sensitivity to the dipole moment of MG, with higher concentrations resulting in a stronger nonlinear optical response. Additionally, the assumption that hyperpolarizability must be treated as a complex quantity is supported by the observed results. The imaginary part of ${\alpha }_5$ likely contributes significantly due to changes in the optical penetration depth and phase modulation effects.

3.2.2. Rhodamine B/Si(001)

Rhodamine B also exhibits $C_{2v}$ symmetry [17,26]. Using the SBHM approach, the susceptibility tensor for rhodamine B is expressed as

$${\chi }_{Rh}={\alpha }_{Rh}\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& -R\\ 0& 0& 0\\ -R& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}-R& 0& 0\\ 0& 0& 0\\ 0& 0& 2\left(1-{\displaystyle \frac{{sin}^3\left(\gamma \right)}{3\sqrt{3}}}\right)\end{array}\right)\end{array}\right),$$

(16)

where $R=-{\displaystyle \frac{2{cos}^2\gamma sin\gamma }{3\sqrt{3}}}$. By applying the SBHM framework, the SHG intensity simulation results for rhodamine B on Si(001) are obtained, as shown in Figure 7.

The SBHM simulation results closely align with the experimental findings reported in Ref. [17]. In these simulations, rhodamine B molecules adopt a perpendicular orientation, with the two amine groups positioned near the silica surface, while the aromatic group extends outward [17]. Notably, the SHG intensity patterns observed in Ref. [26] resemble the SBHM simulation results in Figure 7a,c. These findings confirm that rhodamine B’s ${\widehat{b}}_8$ is perpendicular to the x-axis and that orientation significantly contributes to the observed SHG signal.

3.2.3. SBHM Fitting Versus Experimental Data Error Analysis

While the sensitivity of the SHG intensity to minor variations in bond angles or molecular orientation is significant, and thermal vibrations or rotational diffusion may influence SHG signal stability, these factors were not the primary focus of this study. The analysis was based on the most preferred molecular orientation of the pollutants, as identified in the literature, which allowed us to achieve an acceptable accuracy level of up to 3% error. Future studies employing Density Functional Theory (DFT) calculations could provide a more detailed understanding of these variations, although such an approach would require significantly more complex computational efforts.

Additionally, we extracted the raw data from Refs. [9,17] and compared them with our model. This comparison yielded a deviation of only 3%, demonstrating the robustness of our approach. These results highlight the importance of considering both concentration-dependent dipole moments and the complex nature of hyperpolarizability when modeling the SHG intensity in systems with varying pollutant concentrations. The error calculation for the simulation results compared to the experimental results can be formulated as shown in Equation (17):

$$error={\displaystyle \frac{S{S}_{residu}}{S{S}_{total}}}={\displaystyle \frac{{\sum }_{i=1}^n{\left(I_{Experiment,i}-I_{Simulation,i}\right)}^2}{{\sum }_{i=1}^n{\left(I_{Experiment,i}-{\overline{I}}_{Experiment,i}\right)}^2}}$$

(17)

where $S{S}_{\mathrm{res}}$ is the sum of squares of residuals, and $S{S}_{\mathrm{total}}$ is the total sum of squares. $I_{\mathrm{experiment}}$ represents the SHG intensity obtained from the experiment, while $I_{\mathrm{simulation}}$ is the intensity obtained from the simulation. Additionally, ${\overline{I}}_{\mathrm{experiment}}$ denotes the average intensity from the experimental results.

3.2.4. Molecular Docking Simulation of Rhodamine B Adsorbate Using a GO Layer

As previously explained, the inclusion of a GO layer will increase the SHG contribution through the increase in the local field effect. Increasing the SHG intensity is important when the pollutant concentration is hard to detect since the SHG effects are much lower than their linear response. Jannis et al. (2020) demonstrated that multi-layer GO influences SHG through local field-induced molecular orientation [32]. Carr et al. (2022) showed how local electric fields near GO surfaces align water molecules, enhancing the SHG intensity [33]. Similarly, Zhou et al. (2024) investigated laser-induced heating effects, further supporting the role of GO in SHG enhancement [34]. Carr et al. (2023) correlated SHG intensity variations with molecular alignment [35]. Additionally, Autere et al. (2018) emphasized symmetry-breaking effects in GO as critical factors in nonlinear optical responses, including SHG [36]. Another mechanism of increasing SHG due to GO is the occurance of surface plasmon resonance [37]. Therefore, it is crucial to conduct an MD simulation to determine wether the pollutant can be absorbed or attached by the GO layer.

The molecular docking simulation results identified the strongest binding affinity as $-10.4\mathrm{kcal}/\mathrm{mol}$. Figure 8 illustrates the binding pose of rhodamine B on the GO surface for this highest affinity mode. Rhodamine B is tightly bound near the center of the GO surface. Other binding modes similarly occupy regions around the mid-surface of GO, with comparable binding affinities ranging between $-10.1$ and $-10.4\mathrm{kcal}/\mathrm{mol}$. Details of the molecular interactions are summarized in

Table 2. Molecular interactions between rhodamine B and GO for the strongest binding mode.

Color Symbol	Type of Interaction	Number of Interactions	Distances (Å)
Black (Lavender)	Alkyl	1	4.04
Pink	$\pi$-Alkyl	5	4.48, 4.90, 4.32, 5.10, 4.98
Purple	$\pi$-Sigma	3	3.47, 3.67, 3.81
Red	$\pi$-$\pi$ T-shaped	3	4.87, 4.87, 4.58
Magenta	$\pi$-$\pi$ Stacked	36	4.08, 4.09, 4.91, 5.53, 5.09, 3.75, 4.42, 3.86, 3.90, 4.14, 4.96, 5.51, 5.14, 4.63, 4.32, 5.51, 3.97, 4.60, 5.30, 5.03, 5.41, 5.81, 5.24, 5.38, 4.60, 5.10, 3.99, 4.74, 5.07, 4.53, 4.83, 5.67, 5.46, 5.36, 5.04, 5.19

.

The interactions are predominantly governed by long-range, weak interactions, including $\pi$-$\pi$ stacking, $\pi$-alkyl, $\pi$-sigma, $\pi$-$\pi$ T-shaped, and alkyl interactions. These interaction types suggest that rhodamine B undergoes chemisorption onto the GO surface [38].

Among these, $\pi$-$\pi$ stacking is the most significant, with 36 instances observed. This non-covalent bonding arises from interactions between the electron clouds of parallel aromatic rings. The abundance of aromatic rings on the GO surface facilitates strong bonding with rhodamine B, making GO a promising adsorbent for removing organic contaminants like rhodamine B from water [39,40,41].

Finally, we would like to mention that a parallel-to-the-surface alignment of the effective bonds may increase the SHG pollutant contribution further due to the following argument. In the case of rhodamine B adsorbed on graphene oxide (GO), the parallel alignment of the aromatic rings to the surface facilitates a stronger interaction between the incoming driving electric field component and the molecular dipole moments because they are aligned more parallel to each other (especially if the field is normal to the surface), thus producing a stronger anharmonic oscillation. This configuration enhances the nonlinear polarization, resulting in a stronger SHG signal. This phenomenon is similar to the SHG behavior observed on Si(111) surfaces, as described in Ref. [6]. The parallel alignment effectively maximizes the overlap between the driving field and the molecular bond vectors, thereby amplifying the SHG radiation.

4. Conclusions

Based on the simulations conducted in this study, investigations of a silicon (001) surface with MG revealed that both the experimental and simulation results produced identical SHG intensity patterns across all polarization angles. This finding indicates a polarization contribution arising from both the silicon substrate and the MG molecules on the surface. Additionally, the structural configuration of MG, resembling an inverted Y at its center, provides a useful reference for future research on nanoscale interactions and surface dynamics.

Similarly, the SHG intensity patterns for rhodamine B on silica surfaces showed strong agreement between the experimental and simulation results, further validating the model’s predictive capability.

Through molecular docking simulations, this study also demonstrated that graphene oxide (GO) is highly effective in adsorbing rhodamine B contaminants, exhibiting a strong binding affinity of $-10.4\mathrm{kcal}/\mathrm{mol}$. The adsorption mechanism was classified as chemisorption, as evidenced by the dominance of $\pi$-$\pi$ stacking and $\pi$-alkyl interactions. These findings highlight the potential of a GO coating on silicon layers to serve as an effective adsorbent for detecting and mitigating waterborne contaminants like rhodamine B.

Overall, the combination of nonlinear bond modeling and molecular simulations provides a robust framework for understanding and optimizing nanoscale organic contaminant detection mechanisms on various surfaces. Future research may expand on these insights to explore additional contaminant molecules and surface configurations.