Parametric Characterization of Nonlinear Optical Susceptibilities in Four-Wave Mixing: Solvent and Molecular Structure Effects

1. Introduction

Semiclassical models are used to study the nonlinear optical properties in complex system [1,2,3,4]. It is very common to model the active system subject to interaction with electric fields as a set of two levels that can represent the optical transition and its transitions are evaluated through a Gaussian or Lorentzian distribution of inhomogeneities. This type of modeling has the advantage of simplicity since it does not require an internal structure of the states that imply different transitions. Within this framework, a theoretical study carried out by Hesabi et al. [5], focused on the study of the information flow for a field with lorentzian spectrum through the dynamics of a two-level system immersed in a fluctuating classical field interacting with its environment through dipolar interaction, revealed the existence of two work regimes corresponding to the broad and sharp spectrum of the field, in which the memory effects are governed by the energy difference between two levels or by the interaction energy. Sakiroglu [6] studied the third-order linear and nonlinear optical absorption coefficients induced by an electric field and the changes in the refractive index in a quantum well, represented by the Morse potential, fusing a two-level system. He found that the total optical absorption increases and resonant peaks shift to higher energy values with an increase in electric field and decrease in the well parameter that describes the confining potential range and controls its asymmetry, also a decrease in the change in total refractive index with the well structural parameter was observed. Mukamel and Abramavicius [7] performed theoretical calculations of polarization and nonlinear optical susceptibilities to simulate nonlinear chromophore spectroscopies using the many-body approach. Pishchalnikov [8] used a brownian oscillator model to include stochastic effects in his study for the simulation of linear optical response in photosynthetic pigments such as chlorophyll and bacteriochlorophyll, employing 49 vibronic modes for chlorophyll and 60 modes for bacteriochlorophyll; he calculated parameters of electron-phonon coupling and lower frequency modes, parameters that can be used to model the nonlinear optical response under light presence.

Paz et al. [9] have already carried out studies on the nonlinear optical responses in a saturation regime, incorporating a molecular structure in the two-level model through the insertion of spin-orbit residual terms of the total Hamiltonian of the system and including stochastic effects, observing modifications in the wave functions that describe the system, and allowing the presence of permanent dipole moments in molecules that present parity symmetry. Following this last line, other theoretical studies have now been considered to account molecular structure through a possible vibronic coupling, using both cumulant expansions in Gaussian distributions and without them [10,11,12,13]. Paz et al. [14] have reported, by constructing a kinetic model for the description of the spatial propagation of electromagnetic fields in four-wave mixing spectroscopy [15,16,17,18,19,20] in a two-level molecular model, that the nonlinear optical processes of absorption and scattering can be described by an analogy with the kinetic processes of a chemical reaction. They found that an apparent equilibrium constant can be defined that regulates the competition of the photonic processes involved, similar to Einstein's model of spontaneous emission, but extended to nonlinear optical processes and relaxation. Four-wave mixing has also been employed in frequency measurement and conversion, spatial signal processing, etc. [21,22,23,24,25,26]. The characterization of molecular systems of optoelectronic interest and quantum control with various technological applications have been the subject of study and analysis by different authors. Al Saidi et al. [27], have reported modifications in the sign and magnitude of the nonlinear refractive index, nonlinear absorption coefficient, and third-order nonlinear optical susceptibility with the chemical concentration of Giemsa dye, because of its high saturable absorption, which is useful for developing optical limiters. Cheng et al. [28], analyzed the optical properties of graphene to account for multiphotonic processes occurring in its band structure, reflected in the changes of electrical conductivity allowing to exploitation of the potential of this organic molecule in the development of photonic and optoelectronic devices such as: optical modulators, optical limiters, photodetectors, and saturable absorbers. Yao et al. [29], studied the optical properties of 2,9,16,23-phenoxyphthalocyanine and 2,9,16,23-phenoxyphthalocyanine-zinc in solution by nanosecond laser pulse excitation, finding that the optical limiting efficiency depends on the concentration, exhibiting good all-optical switching properties for potential application in signal processing. Siva et al. [30] have reported the synthesis and evaluation of nonlinear optical properties of a new metal-organic single crystal of bis[hexaacuacobalt(II) sulfate] 2-aminopyridinium monohydrate, finding that it possesses high nonlinearity, attributed to a saturable absorption process, an essential property for its use in devices of opto-electronic nature. Naik et al. [31] synthesized and crystallized two chalcones to investigate their structural, thermal, linear, and nonlinear optical properties. Among the third-order nonlinear optical parameters they determined were the absorption coefficient, nonlinear refractive index, susceptibility, and molecular hyperpolarizability. The result they obtained suggest that thiophene-chalcones could be used in frequency generators, switches, and optical limiters.

From the theoretical point of view, to reach all these above developments, we require more detailed modeling of the molecular structure with less pretensions in the formality of the electromagnetic field [32,33]. The inclusion of the electron-nuclear correlation in our study systems allows us to determine the sensitivity of the study parameters on the coherent and incoherent components of the optical responses associated with the dispersive and absorptive properties. Finally, the nonlinear optical responses show little sensitivity to the saturation order of the weak beam and to influences from both electron-nucleus correlation and solvent. More significant in the behavior of the optical response are the relaxation processes strictly controlled by the ${\mathrm{T}}_1\mathrm{/}{\mathrm{T}}_2$ratio. The importance of vibronic coupling in these studies is determined by the extension that can be made to molecular systems with permanent dipole moments without the need to resort explicitly to the rotating wave approximation.

2. Theory

A theoretical model involving the derivation of the optical Bloch equations through Liouville's formulation was developed. For this purpose, the Hamiltonians of the system, reservoir, and the two interaction potentials were considered. In the inclusion of the internal structure, adiabatic states of the chemical solute were constructed, considering the wave function as a linear combination of the electronic and vibrational states, and calculating on a new basis the adiabatic states responsible for the Bohr frequency and the permanent and transition dipole moments. Finally, it was also considered that the solvent induces shifts in the Bohr frequency of the molecular system, converting the conventional optical Bloch equations into stochastic ones. In this sense, our proposal incorporates from its base, two considerations: 1) the Bohr frequency shift to a time-dependent function, as a solvent effect in random terms, 2) the crossing of curves by residual Hamiltonian effects, which incorporates the electron-nuclear correlation.

2.1. Optical Bloch Equations in the Density Matrix Formalism

To describe the temporal dynamics of the systems, we start with the Liouville-von Newmann equation $\mathrm{i}\mathrm{\hslash }{\mathrm{\partial }}_{\mathrm{t}}\mathrm{\rho }(\mathrm{t})\mathrm{-}\mathrm{\left[}\mathrm{H}\mathrm{,}\mathrm{\rho }(\mathrm{t})\mathrm{\right]}\mathrm{=}0$, where $\mathrm{\rho }(\mathrm{t})$is the global density matrix, and where the Hamiltonian is defined as $\mathrm{H}\mathrm{=}{\displaystyle \mathrm{\sum }_{\mathrm{n}}{\mathrm{H}}_{\mathrm{n}}\mathrm{+}{\displaystyle \mathrm{\sum }_{\mathrm{m}\mathrm{\ne }\mathrm{n}}{\mathrm{H}}_{\mathrm{m}\mathrm{n}}}}$for $\mathrm{n}\mathrm{=}\mathrm{S}\mathrm{,}\mathrm{F}\mathrm{,}\mathrm{R}$(S: System, F:Field, R: Reservoir), with ${\mathrm{H}}_{\mathrm{m}\mathrm{n}}$represents the Hamiltonian of the interactions. In our model we do not incorporate the F-R interaction because it is out of resonance in the frequency absorption and we do not take the quantization of F. Considering the quantum theory of relaxation and applying the model to a two-level system, the equations of the dynamics of the density matrix components are finally obtained:

$${\mathrm{\partial }}_{\mathrm{t}}\mathrm{\rho }(\mathrm{t})\mathrm{=}\mathrm{A}(\mathrm{t})\mathrm{\rho }(\mathrm{t})\mathrm{+}{\mathrm{B}}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{x}}$$

(1)

$$\mathrm{where}\mathrm{\rho }(\mathrm{t})\mathrm{=}\mathrm{\left(}\begin{array}{c}{\mathrm{\rho }}_{\mathrm{D}}\\ {\mathrm{\rho }}_{\mathrm{b}\mathrm{a}}\\ {\mathrm{\rho }}_{\mathrm{a}\mathrm{b}}\end{array}\mathrm{\right)}\mathrm{;}{\mathrm{B}}_{\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{x}}\mathrm{=}\mathrm{\left(}\begin{array}{c}{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}\mathrm{.}}\mathrm{/}{\mathrm{T}}_1\\ 0\\ 0\end{array}\mathrm{\right)}\mathrm{,}$$

(2.a)

$$\mathrm{and}\mathrm{A}(\mathrm{t})\mathrm{=}\mathrm{\left(}\begin{array}{ccc}\mathrm{-}{\displaystyle \frac{1}{{\mathrm{T}}_1}}& \mathrm{-}{\displaystyle \frac{2\mathrm{i}}{\mathrm{\hslash }}}{\mathrm{V}}_{\mathrm{a}\mathrm{b}}& \mathrm{-}{\displaystyle \frac{2\mathrm{i}}{\mathrm{\hslash }}}{\mathrm{V}}_{\mathrm{a}\mathrm{b}}\\ \mathrm{-}{\displaystyle \frac{\mathrm{i}}{\mathrm{\hslash }}}{\mathrm{V}}_{\mathrm{b}\mathrm{a}}& \mathrm{-}\mathrm{\left(}\mathrm{i}{\mathrm{\omega }}_0\mathrm{+}{\displaystyle \frac{1}{{\mathrm{T}}_2}}\mathrm{\right)}\mathrm{+}{\displaystyle \frac{\mathrm{i}}{\mathrm{\hslash }}}\mathrm{\left(}{\mathrm{V}}_{\mathrm{a}\mathrm{a}}\mathrm{-}{\mathrm{V}}_{\mathrm{b}\mathrm{b}}\mathrm{\right)}& 0\\ {\displaystyle \frac{\mathrm{i}}{\mathrm{\hslash }}}{\mathrm{V}}_{\mathrm{a}\mathrm{b}}& 0& \mathrm{\left(}\mathrm{i}{\mathrm{\omega }}_0\mathrm{-}{\displaystyle \frac{1}{{\mathrm{T}}_2}}\mathrm{\right)}\mathrm{-}{\displaystyle \frac{\mathrm{i}}{\mathrm{\hslash }}}\mathrm{\left(}{\mathrm{V}}_{\mathrm{a}\mathrm{a}}\mathrm{-}{\mathrm{V}}_{\mathrm{b}\mathrm{b}}\mathrm{\right)}\end{array}\mathrm{\right)}$$

(2.b)

In eq. (1) we use the frequency domain to avoid possible future complications related to memory and coherence issues. For this purpose, we performed Fourier series expansions of the components of the density matrix. From these expressions, it can be seen how the components of the density matrix are related to each other in the frequency space. The coherent superposition of the eigenstates $\mathrm{|}\mathrm{b}\mathrm{\rangle }$ and $\mathrm{|}\mathrm{a}\mathrm{\rangle }$ induced by the radiative interaction oscillates at the frequencies of the incident beams and of the beam generated in the medium in which the molecular system is located; the component ${\mathrm{\rho }}_{\mathrm{a}\mathrm{b}}$ oscillates at the negative of the oscillation frequencies of ${\mathrm{\rho }}_{\mathrm{b}\mathrm{a}}$ and ${\mathrm{\rho }}_{\mathrm{D}}$ oscillates at a combination of both. Here, ${\mathrm{\omega }}_0$is the Bohr frequency, ${\mathrm{T}}_1\mathrm{,}{\mathrm{T}}_2$are the longitudinal and transversal relaxation times, respectively, and ${\mathrm{V}}_{\mathrm{a}\mathrm{b}}\mathrm{,}{\mathrm{V}}_{\mathrm{a}\mathrm{a}}\mathrm{,}{\mathrm{V}}_{\mathrm{b}\mathrm{b}}$are the nondiagonal, and diagonal potentials, respectively. The superscripts (eq) denote the equilibrium value of the population difference.

2.2. Vibronic Coupling and Stochastic Nature due to Solvent Presence

The interaction between nuclear and electronic motions in molecules gives rise to the crossing of two or more potential energy curves, a process associated with intramolecular coupling. The effect of coupling on the response of the molecular system of interest can be treated by considering wave functions on a new basis, composed of a linear combination of the product of electronic and vibrational wave functions [34,35,36,37]. For the consideration of vibronic coupling in the model, we use the previous results as described in ref. [9] which take into account for the energies of the adiabatic states, the reconfiguration of the Bohr frequency (${\tilde{\mathrm{\omega }}}_0$), and the new permanent and transition dipole moments (eqs.10-13) [13]. On the other hand, it was considered that the solvent, besides acting as a thermal reservoir, induces random shifts in the Bohr frequency of the molecular system under study. To address this situation, the natural frequency of the system is now considered to have a stochastic behavior of the form $\mathrm{\xi }(\mathrm{t})\mathrm{=}{\tilde{\mathrm{\omega }}}_0\mathrm{+}\mathrm{\sigma }(\mathrm{t})$. In this expression, all the stochasticity of ξ(t) has been introduced in the function σ(t). The latter function generates an inhomogeneous broadening of the upper energy level. Frequency fluctuations occur around the transition frequency ${\tilde{\mathrm{\omega }}}_0$, so one has $\mathrm{\langle }\mathrm{\sigma }(\mathrm{t})\mathrm{\rangle }\mathrm{=}0$. In solving this problem, an Ornstein-Uhlenbeck process (OUP) is considered, taking the reservoir as a colored noise source [38]. The noise considered has an intensity γ and an exponential correlation function with a decay rate τ, i.e., $\mathrm{\langle }\mathrm{\xi }(\mathrm{t})\mathrm{\xi }(\mathrm{t}\mathrm{\text{'}})\mathrm{\rangle }\mathrm{=}\mathrm{\gamma }\mathrm{\tau }\exp (\mathrm{-}\mathrm{\tau }\mathrm{\left|}\mathrm{t}\mathrm{-}\mathrm{t}\mathrm{\text{'}}\mathrm{\right|})$. In this distribution, the variance is given by the product $\mathrm{\gamma }\mathrm{\tau }$. Taking these considerations into account, by replacing the Bohr frequency of the system with the new time-dependent frequency ξ(t), the conventional optical Bloch equations are transformed into stochastic optical Bloch equations. It is important to note that ${\tilde{\mathrm{\omega }}}_0$ refers to the adiabatic states (A, B) by the insertion of the vibronic coupling.

2.3. Nonlinear Optical Susceptibilities: Stochastic and Intramolecular Considerations

Using eq.(2) for the diagonal and extra-diagonal frequency components of the reduced density matrix in the presence of the intramolecular coupling effects due to the internal molecular nature, and of the stochasticity due to the presence of the solvent and its strong collisions with the solute and the consequent shift in the deterministic frequency ${\tilde{\mathrm{\omega }}}_0$, we obtain the corresponding Fourier components at frequencies ${\mathrm{\omega }}_1\mathrm{,}{\mathrm{\omega }}_2\mathrm{,}{\mathrm{\omega }}_3$, necessary for the evaluation of the nonlinear polarizations induced in the medium. For this, we do: $\overrightarrow{\mathrm{E}}(\mathrm{t})\mathrm{=}{\displaystyle {\mathrm{\sum }}_{\mathrm{n}\mathrm{=}1}^3{\overrightarrow{\mathrm{E}}}_{\mathrm{n}}}({\mathrm{\omega }}_{\mathrm{n}})\exp (\mathrm{-}\mathrm{i}{\mathrm{\omega }}_{\mathrm{n}}\mathrm{t})$ where ${\overrightarrow{\mathrm{E}}}_{\mathrm{n}}({\mathrm{\omega }}_{\mathrm{n}})\mathrm{=}{\overrightarrow{\mathrm{E}}}_{\mathrm{n}0}\exp (\mathrm{i}{\overrightarrow{\mathrm{k}}}_{\mathrm{n}}\mathrm{.}\overrightarrow{\mathrm{r}}\mathrm{+}\mathrm{\varphi })$ and taking as interaction potential the one given by the electric dipole approximation ${\mathrm{V}}_{\mathrm{n}\mathrm{m}}\mathrm{=}\mathrm{-}{\overrightarrow{\mathrm{\mu }}}_{\mathrm{n}\mathrm{m}}\mathrm{.}\overrightarrow{\mathrm{E}}(\mathrm{t})$[39]. Considering the strong pump beam (${\overrightarrow{\mathrm{E}}}_1$) in saturation regimes, the probe beam (${\overrightarrow{\mathrm{E}}}_2$) in the second order and the emerging signal beam (${\overrightarrow{\mathrm{E}}}_3$) in the first order, under perturbative schemes, we have for the FWM signal:

$${\mathrm{\rho }}_{\mathrm{B}\mathrm{A}}({\mathrm{\omega }}_1)\mathrm{=}\mathrm{\left[}\mathrm{i}{\tilde{\mathrm{\Omega }}}_1\mathrm{-}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\displaystyle \frac{2\mathrm{i}{\mathrm{J}}^{\mathrm{*}}{\tilde{\mathrm{\Omega }}}_1{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}\mathrm{-}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\displaystyle \frac{2\mathrm{i}{\tilde{\mathrm{\Omega }}}_1^{\mathrm{*}}{\tilde{\mathrm{\Omega }}}_2{\tilde{\mathrm{\Omega }}}_3}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}\mathrm{\left(}{\mathrm{J}}^{\mathrm{*}}{\mathrm{u}}_{3\mathrm{,}\mathrm{-}1}\mathrm{+}\mathrm{J}{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}\mathrm{\right)}\mathrm{\right]}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_1}}$$

(3.a)

$${\mathrm{\rho }}_{\mathrm{B}\mathrm{A}}({\mathrm{\omega }}_2)\mathrm{=}\mathrm{\left[}\mathrm{i}{\tilde{\mathrm{\Omega }}}_2\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\displaystyle \frac{4\mathrm{i}{\tilde{\mathrm{\Omega }}}_1^2{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2{\tilde{\mathrm{\Omega }}}_3^{\mathrm{*}}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}1}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}\mathrm{-}{\displaystyle \frac{2\mathrm{i}\mathrm{J}{\tilde{\mathrm{\Omega }}}_2{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}\mathrm{-}{\displaystyle \frac{2\mathrm{i}\mathrm{J}{\tilde{\mathrm{\Omega }}}_1^2{\tilde{\mathrm{\Omega }}}_3^{\mathrm{*}}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}3}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}\mathrm{\right]}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_2}}$$

(3.b)

$${\mathrm{\rho }}_{\mathrm{B}\mathrm{A}}({\mathrm{\omega }}_3)\mathrm{=}\mathrm{\left[}\mathrm{i}{\tilde{\mathrm{\Omega }}}_3\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\displaystyle \frac{4\mathrm{i}{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2{\tilde{\mathrm{\Omega }}}_3}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}1}{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}\mathrm{-}{\displaystyle \frac{2\mathrm{i}{\mathrm{J}}^{\mathrm{*}}{\tilde{\mathrm{\Omega }}}_1^2{\tilde{\mathrm{\Omega }}}_2^{\mathrm{*}}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}\mathrm{-}{\displaystyle \frac{2\mathrm{i}{\mathrm{J}}^{\mathrm{*}}{\tilde{\mathrm{\Omega }}}_3{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2{\mathrm{u}}_{3\mathrm{,}\mathrm{-}1}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}\mathrm{\right]}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_3}}\mathrm{.}$$

(3.c)

We select the value of $\mathrm{m}\mathrm{=}0$ for the perturbative treatment at the first order for the probe beam and the value $\mathrm{m}\mathrm{=}1$, when considering the second order in the probe beam. In these expressions we have considered the validity of the rotating wave approximation RWA [40]; ${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}$ which represents the canonical population difference in the adiabatic basis in the presence of the saturation effects of the pump and probe beams. We can observe its tendency to the Boltzmann-type equilibrium value in the cases of the absence of radiation. We can notice that with the incorporation of the second order of probe ($\mathrm{m}\mathrm{=}1$), the expressions still have the same symmetry in their terms for the probe and signal beams; for $\mathrm{m}\mathrm{=}0$, the elements of the density matrix ${\mathrm{\rho }}_{\mathrm{B}\mathrm{A}}({\mathrm{\omega }}_2)$ and ${\mathrm{\rho }}_{\mathrm{B}\mathrm{A}}({\mathrm{\omega }}_3)$ maintain a symmetry of terms different from the previous one, and in that case, the pump-beam maintains a mathematical structure that only involves its beam in the photonic processes. We have defined: $\mathrm{J}\mathrm{=}{\mathrm{\Gamma }}_1\mathrm{+}2{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2{\mathrm{u}}_{3\mathrm{,}\mathrm{-}2}\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}2{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2\mathrm{/}{\mathrm{D}}_1$; ${\mathrm{u}}_{\mathrm{i}\mathrm{,}\mathrm{-}\mathrm{k}}\mathrm{=}1\mathrm{/}{\mathrm{D}}_{\mathrm{i}}\mathrm{+}1\mathrm{/}{\mathrm{D}}_{\mathrm{k}}^{\mathrm{*}}$; ${\mathrm{D}}_{\mathrm{n}}\mathrm{=}{\mathrm{T}}_2^{\mathrm{-}1}\mathrm{+}\mathrm{i}(\mathrm{\xi }(\mathrm{t})\mathrm{-}{\mathrm{\omega }}_{\mathrm{n}})$; ${\mathrm{\Gamma }}_{\mathrm{n}}\mathrm{=}{\mathrm{T}}_1^{\mathrm{-}1}\mathrm{-}\mathrm{i}\mathrm{n}\mathrm{\Delta }$; ${\tilde{\mathrm{\Omega }}}_{\mathrm{j}}$ is the Rabi frequency of the beam j in presence of the intramolecular coupling. Using eqs. (3), we evaluate the Fourier components of the complex nonlinear polarizations at different frequencies, according to $\overrightarrow{\mathrm{P}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}\mathrm{N}\mathrm{\left(}{\mathrm{\langle }{\mathrm{\langle }{\mathrm{\rho }}_{\mathrm{B}\mathrm{A}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{\rangle }}_{\mathrm{\xi }}{\overrightarrow{\tilde{\mathrm{\mu }}}}_{\mathrm{A}\mathrm{B}}\mathrm{\rangle }}_{\mathrm{\theta }}\mathrm{\right)}$. Here, the external bracket denotes average over the probability distributions of molecular orientations o the molecules, while the internal bracket implies the average over the set of realizations of the random variable in the range between $\mathrm{\xi }$ and $\mathrm{\xi }\mathrm{+}\mathrm{d}\mathrm{\xi }$, and N the concentration of active molecules interacting with the field. In the tensorial approximation, we express the polarization components of the form:

$$\overrightarrow{\mathrm{P}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\overrightarrow{\mathrm{P}}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{s}}^{\mathrm{C}\mathrm{,}\mathrm{I}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{+}{\overrightarrow{\mathrm{P}}}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{\mathrm{F}\mathrm{W}\mathrm{M}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{,}$$

(4)

where ${\overrightarrow{\mathrm{P}}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{s}}^{\mathrm{C}\mathrm{,}\mathrm{I}}({\mathrm{\omega }}_{\mathrm{j}})$ reflects the contribution of the coherent (C) and incoherent (I) processes to the induced polarization of frequency $({\mathrm{\omega }}_{\mathrm{j}})$, while ${\overrightarrow{\mathrm{P}}}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{\mathrm{F}\mathrm{W}\mathrm{M}}({\mathrm{\omega }}_{\mathrm{j}})$ corresponding to the contribution of the coupling process of the beams at frequencies ${\mathrm{\omega }}_1\mathrm{,}{\mathrm{\omega }}_2\mathrm{,}{\mathrm{\omega }}_3$, whose arrangements in the FWM process generate photons at frequency $({\mathrm{\omega }}_{\mathrm{j}})$. Here, we have defined:

$${\overrightarrow{\mathrm{P}}}_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{s}}^{\mathrm{C}\mathrm{,}\mathrm{I}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\mathrm{\chi }}^{\mathrm{S}\mathrm{V}}({\mathrm{\omega }}_{\mathrm{j}})\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{+}{\mathrm{\chi }}_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathrm{\omega }}_{\mathrm{j}})\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_{\mathrm{j}})$$

(5)

with j=1,2,3

where ${\mathrm{\chi }}_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\mathrm{\chi }}_{\mathrm{I}}^{(1)}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{+}{\mathrm{\chi }}_{\mathrm{C}}^{(3)}({\mathrm{\omega }}_{\mathrm{j}})$, and ${\mathrm{\chi }}^{\mathrm{S}\mathrm{V}}({\mathrm{\omega }}_{\mathrm{j}})$is the solvent contribution. On the other hand, we define:

$${\overrightarrow{\mathrm{P}}}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{\mathrm{F}\mathrm{W}\mathrm{M}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_{\mathrm{j}}){\overrightarrow{\mathrm{E}}}_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{,}$$

(6)

where ${\mathrm{\omega }}_{\mathrm{j}}\mathrm{=}{\mathrm{\omega }}_{\mathrm{q}}\mathrm{+}{\mathrm{\omega }}_{\mathrm{p}}\mathrm{-}{\mathrm{\omega }}_{\mathrm{r}}$ and ${\overrightarrow{\mathrm{E}}}_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\overrightarrow{\mathrm{E}}}_{\mathrm{q}}{\overrightarrow{\mathrm{E}}}_{\mathrm{p}}{\overrightarrow{\mathrm{E}}}_{\mathrm{r}}^{\mathrm{*}}$, valid for the following values: $\mathrm{j}\mathrm{=}1$(pump), with $\mathrm{q}\mathrm{=}2\mathrm{,}\mathrm{p}\mathrm{=}3\mathrm{,}\mathrm{r}\mathrm{=}1$; $\mathrm{j}\mathrm{=}2$ (probe), with $\mathrm{q}\mathrm{=}\mathrm{p}\mathrm{=}1\mathrm{,}\mathrm{r}\mathrm{=}3$; $\mathrm{j}\mathrm{=}3$(signal), with $\mathrm{q}\mathrm{=}\mathrm{p}\mathrm{=}1\mathrm{,}\mathrm{r}\mathrm{=}2$. Thus, we have the optical susceptibilities in incoherent processes (those due to absorption), the following:

$${\mathrm{\chi }}_{\mathrm{I}}^{(1)}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\mathrm{a}}_1{\mathrm{\langle }{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_{\mathrm{j}}}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{,}\mathrm{j}\mathrm{=}1\mathrm{,}2\mathrm{,}3\mathrm{.}$$

(7)

For coherent processes (those due to saturation), we have the following optical susceptibilities:

$${\mathrm{\chi }}_{\mathrm{C}}^{(3)}({\mathrm{\omega }}_1)\mathrm{=}\mathrm{-}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\mathrm{a}}_2{\mathrm{\langle }{\displaystyle \frac{{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}}{\mathrm{J}}}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_1}}\mathrm{\rangle }}_{\mathrm{\xi }}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_2)\mathrm{\right|}}^2\mathrm{;}{\mathrm{\chi }}_{\mathrm{C}}^{(3)}({\mathrm{\omega }}_2)\mathrm{=}\mathrm{-}{\mathrm{a}}_2{\mathrm{\langle }{\displaystyle \frac{{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}}{{\mathrm{J}}^{\mathrm{*}}}}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_2}}\mathrm{\rangle }}_{\mathrm{\xi }}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_1)\mathrm{\right|}}^2$$

$${\mathrm{\chi }}_{\mathrm{C}}^{(3)}({\mathrm{\omega }}_3)\mathrm{=}\mathrm{-}{\mathrm{a}}_2{\mathrm{\langle }{\displaystyle \frac{{\mathrm{u}}_{3\mathrm{,}\mathrm{-}1}}{\mathrm{J}}}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_3}}\mathrm{\rangle }}_{\mathrm{\xi }}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_1)\mathrm{\right|}}^2\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\mathrm{a}}_3{\mathrm{\langle }{\displaystyle \frac{{\mathrm{u}}_{1\mathrm{,}\mathrm{-}1}{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_3}}\mathrm{\rangle }}_{\mathrm{\xi }}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_1)\mathrm{\right|}}^2{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_2)\mathrm{\right|}}^2$$

(8)

For the coupling processes at the different frequencies, we obtain:

$${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_1)\mathrm{=}\mathrm{-}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\mathrm{a}}_2{\mathrm{\langle }\mathrm{\left(}{\displaystyle \frac{{\mathrm{u}}_{3\mathrm{,}\mathrm{-}1}}{\mathrm{J}}}\mathrm{+}{\displaystyle \frac{{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}}{{\mathrm{J}}^{\mathrm{*}}}}\mathrm{\right)}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_1}}\mathrm{\rangle }}_{\mathrm{\xi }}$$

$${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_2)\mathrm{=}\mathrm{-}{\mathrm{a}}_2{\mathrm{\langle }{\displaystyle \frac{{\mathrm{u}}_{1\mathrm{,}\mathrm{-}3}}{{\mathrm{J}}^{\mathrm{*}}}}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_2}}\mathrm{\rangle }}_{\mathrm{\xi }}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_1)\mathrm{\right|}}^2\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\mathrm{a}}_3{\mathrm{\langle }{\displaystyle \frac{{\mathrm{u}}_{1\mathrm{,}\mathrm{-}1}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_2}}\mathrm{\rangle }}_{\mathrm{\xi }}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_2)\mathrm{\right|}}^2$$

$${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_3)\mathrm{=}\mathrm{-}{\mathrm{a}}_2{\mathrm{\langle }{\displaystyle \frac{{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}}{\mathrm{J}}}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}}{{\mathrm{D}}_3}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{.}$$

(9)

where ${\mathrm{a}}_1\mathrm{=}\mathrm{i}{\mathrm{\left|}{\overrightarrow{\tilde{\mathrm{\mu }}}}_{\mathrm{B}\mathrm{A}}\mathrm{\right|}}^2\mathrm{N}\mathrm{/}\mathrm{\hslash }$, ${\mathrm{a}}_2\mathrm{=}2\mathrm{i}{\mathrm{\left|}{\overrightarrow{\tilde{\mathrm{\mu }}}}_{\mathrm{B}\mathrm{A}}\mathrm{\right|}}^4\mathrm{N}\mathrm{/}{\mathrm{\hslash }}^3$, ${\mathrm{a}}_3\mathrm{=}4\mathrm{i}{\mathrm{\left|}{\overrightarrow{\tilde{\mathrm{\mu }}}}_{\mathrm{B}\mathrm{A}}\mathrm{\right|}}^6\mathrm{N}\mathrm{/}{\mathrm{\hslash }}^5$. The superscripts ($\mathrm{k}\mathrm{=}1\mathrm{,}3$) are the minimum order required in the process. We note that the case of $\mathrm{m}\mathrm{=}0$ the probe- and FWM-signal beams acquire a similar but not symmetrical coupling susceptibility in their process, and a photon-generating coupling process at the frequency of the pumping beam is not reflected. That is, the strong intensity beam (pump) acts through a single absorption process in the medium, with no FWM processes that allow the couplings to generate photons at the frequency ${\mathrm{\omega }}_1$ (pump) (that is ${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_1)\mathrm{=}0$). However, for the coupling processes the weak beams are differentiated through the pump intensity (eq. 9), a fact that is reflected by the same perturbative condition given to both beams. Similarly, the coherent process becomes strictly symmetric in the lower intensity beams (probe and signal) and the coherent process disappears at the pump beam frequency ${\mathrm{\omega }}_1$, and again the susceptibilities are symmetrized. It is important to draw attention to these optical responses, since considering a probe beam at second order and differentiating it from the emerging signal beam, allows the appearance of a coupling process that now involves the optical response at the pump frequency, and they become very different for the probe and signal beams. This fact is repeated in the same way in the coherent responses. It is necessary to indicate that the incoherent process, independently of the type of beam considered, remains invariant, since it is only attributed to the absorption processes, independent of the treatment given to the test. Making use of eqs. (3), it is possible to show that:

$${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}\mathrm{=}{\displaystyle \frac{{\mathrm{T}}_2^2{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^2{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^2}{{\mathrm{T}}_2^2{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^2{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^2\mathrm{\left(}1\mathrm{-}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\mathrm{T}}_1{\displaystyle {\mathrm{\sum }}_{\mathrm{j}\mathrm{=}1}^3{\mathrm{\varphi }}_{\mathrm{j}}}\mathrm{\right)}\mathrm{+}4{\tilde{\mathrm{S}}}_1\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}4{\tilde{\mathrm{S}}}_2}}{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}\mathrm{.}}\mathrm{,}$$

(10)

${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}$ is the population difference component at zero frequency, and where the saturation parameters are defined as: ${\tilde{\mathrm{S}}}_1\mathrm{=}{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2{\mathrm{T}}_1{\mathrm{T}}_2{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^2$ and ${\tilde{\mathrm{S}}}_2\mathrm{=}{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2{\mathrm{T}}_1{\mathrm{T}}_2{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^2$, for the pump- and probe-beams, respectively, and in presence of the coupling intramolecular. The functions ${\mathrm{\varphi }}_{\mathrm{j}}$ are defined as follows:

$${\mathrm{\varphi }}_1\mathrm{=}{\displaystyle \frac{4{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}\mathrm{\left(}{\mathrm{J}}^{\mathrm{*}}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}^2\mathrm{+}\mathrm{J}{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}^2\mathrm{\right)}\mathrm{;}{\mathrm{\varphi }}_2\mathrm{=}{\displaystyle \frac{4{\tilde{\mathrm{\Omega }}}_1^2{\tilde{\mathrm{\Omega }}}_2^{\mathrm{*}}{\tilde{\mathrm{\Omega }}}_3^{\mathrm{*}}}{{\mathrm{\left|}\mathrm{J}\mathrm{\right|}}^2}}\mathrm{\left(}\mathrm{J}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}3}{\mathrm{u}}_{2\mathrm{,}\mathrm{-}1}\mathrm{+}{\displaystyle \frac{1}{{\mathrm{D}}_1^{\mathrm{*}}}}{\mathrm{J}}^{\mathrm{*}}{\mathrm{u}}_{1\mathrm{,}\mathrm{-}2}\mathrm{\right)}\mathrm{;}{\mathrm{\varphi }}_3\mathrm{=}{\displaystyle \frac{{\tilde{\mathrm{\Omega }}}_1^2}{{({\tilde{\mathrm{\Omega }}}_1^{\mathrm{*}})}^2}}{\mathrm{\varphi }}_2^{\mathrm{*}}\mathrm{.}$$

However, taking a particular case of organic dyes as our system of study, it is possible to have ${\mathrm{T}}_1{\displaystyle {\mathrm{\sum }}_{\mathrm{j}\mathrm{=}1}^3{\mathrm{\varphi }}_{\mathrm{j}}}\mathrm{\ll }1$, and therefore eq. (10) is reduced to:

$${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}\mathrm{=}{\displaystyle \frac{{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}}}{1\mathrm{+}\frac{4{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2{\mathrm{T}}_1}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^2{\mathrm{T}}_2}\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}\frac{4{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2{\mathrm{T}}_1}{{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^2{\mathrm{T}}_2}}}\mathrm{.}$$

(11)

Indicating that in situations of low pump saturation and in conditions of perturbative treatments at first order of probe, the steel-frequency component of the population difference is defined by a Boltzmann distribution in equilibrium regimes. We also note that for cases where $\mathrm{m}\mathrm{=}1$ we have that: $4\mathrm{\left(}{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^2\mathrm{/}{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^2\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2\mathrm{/}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^2\mathrm{\right)}{\mathrm{T}}_1\mathrm{/}{\mathrm{T}}_2\mathrm{\ll }1$, and can be developed in series and finally obtain:

$${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}\mathrm{=}{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}\mathrm{.}}{\displaystyle \mathrm{\sum }_{\mathrm{n}\mathrm{=}0}^{\mathrm{\infty }}{\mathrm{q}}_{\mathrm{n}}}{\displaystyle \mathrm{\sum }_{\mathrm{k}\mathrm{=}0}^{\mathrm{n}}{\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}}{\displaystyle \frac{1}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}}}{\displaystyle \frac{1}{{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}}}\mathrm{.}$$

(12)

Taking as example the incoherent susceptibility of type ${\mathrm{\chi }}_{\mathrm{I}}^{(1)}({\mathrm{\omega }}_1)$ , we have

$${\mathrm{\chi }}_{\mathrm{I}}^{(1)}({\mathrm{\omega }}_1)\mathrm{=}{\displaystyle \mathrm{\sum }_{\mathrm{n}\mathrm{=}0}^{\mathrm{\infty }}{\mathrm{q}}_{\mathrm{n}}}{\displaystyle \mathrm{\sum }_{\mathrm{k}\mathrm{=}0}^{\mathrm{n}}{\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}}{\mathrm{\lambda }}_1{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}$$

(13)

$${\mathrm{q}}_{\mathrm{n}}\mathrm{=}{(\mathrm{-}1)}^{\mathrm{n}}{\mathrm{\left(}{\displaystyle \frac{4{\mathrm{T}}_1}{{\mathrm{T}}_2}}\mathrm{\right)}}^{\mathrm{n}}\mathrm{,}{\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}\mathrm{=}{\displaystyle \frac{\mathrm{n}\mathrm{!}}{\mathrm{k}\mathrm{!}(\mathrm{n}\mathrm{-}\mathrm{k})\mathrm{!}}}{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^{2\mathrm{k}}\mathrm{,}\mathrm{and}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_1^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k}\mathrm{+}1)}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{.}$$

Considering a Gaussian probability distribution function, we have for ${\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}$in the OUP [41]:

$${\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\displaystyle \frac{1}{\sqrt{2\mathrm{\pi }\mathrm{\gamma }\mathrm{\tau }}}}{\displaystyle \underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathrm{\infty }}{\mathrm{\int }}}\frac{{\mathrm{T}}_2^{\mathrm{-}1}\mathrm{-}\mathrm{i}\mathrm{\left(}\mathrm{\xi }\mathrm{-}{\mathrm{\omega }}_1\mathrm{\right)}}{{\mathrm{\left[}{\mathrm{T}}_2^{\mathrm{-}2}\mathrm{+}{\mathrm{\left(}\mathrm{\xi }\mathrm{-}{\mathrm{\omega }}_1\mathrm{\right)}}^2\mathrm{\right]}}^{\mathrm{n}\mathrm{-}\mathrm{k}\mathrm{+}1}{\mathrm{\left[}{\mathrm{T}}_2^{\mathrm{-}2}\mathrm{+}{\mathrm{\left(}\mathrm{\xi }\mathrm{-}{\mathrm{\omega }}_2\mathrm{\right)}}^2\mathrm{\right]}}^{\mathrm{k}}}}{\mathrm{e}}^{\mathrm{-}{\displaystyle \frac{{\mathrm{\left(}\mathrm{\xi }\mathrm{-}{\mathrm{\omega }}_0\mathrm{\right)}}^2}{2\mathrm{\gamma }\mathrm{\tau }}}}\mathrm{d}\mathrm{\xi }\mathrm{.}$$

(14)

From eq. (14) it is possible to derive the following recursive formulas (real and imaginary parts)

$$\begin{gathered} 2(\mathrm{n}\mathrm{-}\mathrm{k})\mathrm{Re}\mathrm{\left[}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right]}\mathrm{+}2\mathrm{k}\mathrm{Re}\mathrm{\left[}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{,}\mathrm{k}\mathrm{+}1)\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right]}\mathrm{-}{\mathrm{T}}_2^2{\displaystyle \frac{\mathrm{\partial }\mathrm{\left(}{\mathrm{T}}_2\mathrm{Re}\mathrm{\left[}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{-}1\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right]}\mathrm{\right)}}{\mathrm{\partial }{\mathrm{T}}_2}}\mathrm{=}0 \\ \mathrm{Im}\mathrm{\left[}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right]}\mathrm{+}{\displaystyle \frac{{\mathrm{T}}_2}{2(\mathrm{n}\mathrm{-}\mathrm{k})}}{\displaystyle \frac{\mathrm{\partial }\mathrm{Re}\mathrm{\left[}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(\mathrm{n}\mathrm{-}1\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right]}}{\mathrm{\partial }{\mathrm{\Delta }}_1}}\mathrm{=}0\mathrm{n}\mathrm{\ne }\mathrm{k} \end{gathered}$$

(15)

To reduce the computational effort required for this type of increasingly complex integrals as we advance in the development of the orders of perturbations and the type of susceptibilities, we have for $\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(0\mathrm{,}0)\mathrm{\rangle }$ and all those terms that cannot be achieved from the recursive relation, the use of Parseval's formula:

$${\displaystyle \underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathrm{\infty }}{\mathrm{\int }}}{\mathrm{x}}_1(\mathrm{t}){\mathrm{x}}_2^{\mathrm{*}}(\mathrm{t})\mathrm{d}\mathrm{t}\mathrm{=}\frac{1}{2\mathrm{\pi }}}{\displaystyle \underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathrm{\infty }}{\mathrm{\int }}}{\mathrm{X}}_1(\mathrm{\omega }){\mathrm{X}}_2^{\mathrm{*}}(\mathrm{\omega })\mathrm{d}\mathrm{\omega }}$$

(16)

${\mathrm{X}}_1(\mathrm{\omega })$ and ${\mathrm{X}}_2^{\mathrm{*}}(\mathrm{\omega })$ are the Fourier transforms of the complex functions ${\mathrm{x}}_1(\mathrm{t})$ and ${\mathrm{x}}_2^{\mathrm{*}}(\mathrm{t})$, respectively. Taking a particular case to illustrate the calculation, we have:

$$\mathrm{Re}\mathrm{\left[}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(1\mathrm{,}1)\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right]}\mathrm{=}{\mathrm{b}}^{\mathrm{\Delta }}{\mathrm{e}}^{{\displaystyle \frac{{\mathrm{T}}_2^{\mathrm{-}2}\mathrm{-}{\mathrm{\Delta }}_1^2}{2\mathrm{\gamma }\mathrm{\tau }}}}\mathrm{[}\mathrm{(}{\mathrm{T}}_2\mathrm{\Delta }\cos \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{-}2\sin \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{)}\mathrm{+}\mathrm{\left(}{\mathrm{T}}_2\mathrm{\Delta }\cos \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{+}2\sin \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_2}\mathrm{\right)}\mathrm{\right)}\mathrm{]}$$

$$\mathrm{Im}\mathrm{\left[}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(1\mathrm{,}1)\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right]}\mathrm{=}\mathrm{-}{\mathrm{b}}^{\mathrm{\Delta }}\mathrm{\left[}{\mathrm{e}}^{{\displaystyle \frac{{\mathrm{T}}_2^{\mathrm{-}2}\mathrm{-}{\mathrm{\Delta }}_1^2}{2\mathrm{\gamma }\mathrm{\tau }}}}\mathrm{(}2\cos \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{+}{\mathrm{T}}_2\mathrm{\Delta }\sin \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{)}\mathrm{-}{\mathrm{e}}^{{\displaystyle \frac{{\mathrm{T}}_2^{\mathrm{-}2}\mathrm{-}{\mathrm{\Delta }}_2^2}{2\mathrm{\gamma }\mathrm{\tau }}}}\mathrm{\left(}\mathrm{\left(}{\mathrm{T}}_2^2{\mathrm{\Delta }}^2\mathrm{+}2\mathrm{\right)}\cos \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_2}\mathrm{\right)}\mathrm{+}{\mathrm{T}}_2\mathrm{\Delta }\sin \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_2}\mathrm{\right)}\mathrm{\right)}\mathrm{\right]}$$

(17)

where ${\mathrm{b}}^{\mathrm{\Delta }}\mathrm{=}\sqrt{{\displaystyle \frac{\mathrm{\pi }}{2\mathrm{\gamma }\mathrm{\tau }}}}\mathrm{\left(}{\displaystyle \frac{{\mathrm{T}}_2}{\mathrm{\Delta }\mathrm{\left(}4\mathrm{+}{\mathrm{T}}_2^2{\mathrm{\Delta }}^2\mathrm{\right)}}}\mathrm{\right)}$; ${\mathrm{b}}^{{\mathrm{\Delta }}_{\mathrm{j}}}\mathrm{=}{\mathrm{\Delta }}_{\mathrm{j}}{\mathrm{T}}_2^{\mathrm{-}1}\mathrm{/}\mathrm{\gamma }\mathrm{\tau }$ , with ${\mathrm{\Delta }}_{\mathrm{j}}\mathrm{=}{\mathrm{\omega }}_{\mathrm{j}}\mathrm{-}{\tilde{\mathrm{\omega }}}_0$ for j = 1,2 .

In the case of $\mathrm{\Delta }\mathrm{=}{\mathrm{\Delta }}_1\mathrm{-}{\mathrm{\Delta }}_2\mathrm{=}0$ associated with the maximum population oscillations condition, (MPO) [42], eqs. (17) have non indeterminate and analytically closed solutions, given by:

$$\underset{\mathrm{\Delta }\mathrm{\to }0}{\lim }\mathrm{\left(}\mathrm{Re}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(1\mathrm{,}1)\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right)}\mathrm{=}\mathrm{-}{\displaystyle \frac{1}{2{\mathrm{\left(}\mathrm{\gamma }\mathrm{\tau }\mathrm{\right)}}^{3\mathrm{/}2}}}\exp \mathrm{\left(}{\displaystyle \frac{1\mathrm{-}{\mathrm{T}}_2^2{\mathrm{\Delta }}_1^2}{2\mathrm{\gamma }\mathrm{\tau }{\mathrm{T}}_2}}\mathrm{\right)}\mathrm{\left[}\mathrm{\left(}1\mathrm{-}\mathrm{\gamma }\mathrm{\tau }{\mathrm{T}}_2^2\mathrm{\right)}\cos \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{-}\mathrm{\left(}{\mathrm{T}}_2{\mathrm{\Delta }}_1\mathrm{\right)}\sin \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{\right]}$$

$$\underset{\mathrm{\Delta }\mathrm{\to }0}{\lim }\mathrm{\left(}\mathrm{Im}{\mathrm{\langle }{\mathrm{\Phi }}^{(1)}(1\mathrm{,}1)\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right)}\mathrm{=}{\displaystyle \frac{1}{2{\mathrm{\left(}\mathrm{\gamma }\mathrm{\tau }\mathrm{\right)}}^{3\mathrm{/}2}}}\exp \mathrm{\left(}{\displaystyle \frac{1\mathrm{-}{\mathrm{T}}_2^2{\mathrm{\Delta }}_1^2}{2\mathrm{\gamma }\mathrm{\tau }{\mathrm{T}}_2}}\mathrm{\right)}\mathrm{\left[}\mathrm{\left(}{\mathrm{T}}_2{\mathrm{\Delta }}_1\mathrm{\right)}\cos \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{+}\sin \mathrm{\left(}{\mathrm{b}}^{{\mathrm{\Delta }}_1}\mathrm{\right)}\mathrm{\right]}$$

(18)

It is important to note that the real and imaginary components of the optical responses observed in frequency space obey the detuning conditions of each of the beams. For our derivative, we have estimated non-degenerate conditions at the frequencies of the incident beams, so that when we try to study the MPO, the values are indeterminate. However, if the treatment starts under these resonance conditions, the result obtained in the optical responses is as indicated by eqs. (18), showing clearly that the magnitude of the signal is only obtained analytically when the calculation is subject to the MPO condition. Our results allow us to explore any region of the frequency space ${\mathrm{\Delta }}_1\mathrm{,}{\mathrm{\Delta }}_2$, subject to RWA, but under simultaneous considerations of stochastic effects caused by the presence of the solvent and intramolecular couplings because of the molecular structure of the active substrate in the radiation-matter interaction. The present study differs substantially from previous works, which not only fails to explore the simultaneity of these effects, but restricts the optical response only to considerations in MPO regimes. It should be noted that these terms correspond to the first order of the series development for the ${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}$ component. In general, it is possible to obtain:

$${\mathrm{\chi }}_{\mathrm{I}}^{(1)}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\displaystyle \mathrm{\sum }_{\mathrm{n}\mathrm{=}0}^{\mathrm{\infty }}{\mathrm{q}}_{\mathrm{n}}}{\displaystyle \mathrm{\sum }_{\mathrm{k}\mathrm{=}0}^{\mathrm{n}}{\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}}{\mathrm{\lambda }}_1{\mathrm{\langle }{\mathrm{\Phi }}^{(\mathrm{j})}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{j}\mathrm{=}1\mathrm{,}2\mathrm{,}3$$

(19)

with ${\mathrm{\langle }{\mathrm{\Phi }}^{(\mathrm{j})}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_{\mathrm{j}}^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}{\mathrm{\left|}{\mathrm{D}}_{\mathrm{j}}\mathrm{\right|}}^2}}\mathrm{\rangle }}_{\mathrm{\xi }}$. In the case of the coherent susceptibilities, we obtain:

$${\mathrm{\chi }}_{\mathrm{C}}^{(3)}({\mathrm{\omega }}_{\mathrm{j}})\mathrm{=}{\displaystyle \mathrm{\sum }_{\mathrm{n}\mathrm{=}0}^{\mathrm{\infty }}{\mathrm{q}}_{\mathrm{n}}}{\displaystyle \mathrm{\sum }_{\mathrm{k}\mathrm{=}0}^{\mathrm{n}}{\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}}\mathrm{[}\mathrm{\left(}{\mathrm{\lambda }}_2{\mathrm{\delta }}_{2\mathrm{,}\mathrm{j}}\mathrm{-}{\mathrm{\lambda }}_2^{\mathrm{*}}{\mathrm{\delta }}_{3\mathrm{,}\mathrm{j}}\mathrm{\right)}\mathrm{\left(}{\mathrm{\langle }{\mathrm{G}}_1^{(\mathrm{j})}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{+}{\mathrm{\langle }{\mathrm{G}}_2^{(\mathrm{j})}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right)}$$

(20)

$$\mathrm{+}\mathrm{\left(}{\mathrm{\lambda }}_3{\mathrm{\delta }}_{3\mathrm{,}\mathrm{j}}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}\mathrm{\right)}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_1)\mathrm{\right|}}^2\mathrm{\left(}{\mathrm{\langle }{\mathrm{G}}_3^{(3)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{+}{\mathrm{\langle }{\mathrm{G}}_4^{(3)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right)}\mathrm{/}\mathrm{J}\mathrm{]}$$

Where

$${\mathrm{\langle }{\mathrm{G}}_1^{(\mathrm{j})}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{\left(}{\mathrm{D}}_{\mathrm{j}}^{\mathrm{*}}\mathrm{\right)}}^2}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}{\mathrm{\left|}{\mathrm{D}}_{\mathrm{j}}\mathrm{\right|}}^4}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{;}{\mathrm{\langle }{\mathrm{G}}_2^{(\mathrm{j})}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_1{\mathrm{D}}_{\mathrm{j}}^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^2{\mathrm{\left|}{\mathrm{D}}_{\mathrm{j}}\mathrm{\right|}}^2}}\mathrm{\rangle }}_{\mathrm{\xi }}$$

$${\mathrm{\langle }{\mathrm{G}}_3^{(3)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_1^{\mathrm{*}}{\mathrm{D}}_3^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^2{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^2{\mathrm{\left|}{\mathrm{D}}_3\mathrm{\right|}}^2}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{;}{\mathrm{\langle }{\mathrm{G}}_4^{(3)}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_1{\mathrm{D}}_3^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^4{\mathrm{\left|}{\mathrm{D}}_3\mathrm{\right|}}^2}}\mathrm{\rangle }}_{\mathrm{\xi }}$$

For the coupling susceptibility at frequency (${\mathrm{\omega }}_2$), we have

$${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_2)\mathrm{=}{\displaystyle \mathrm{\sum }_{\mathrm{n}\mathrm{=}0}^{\mathrm{\infty }}{\mathrm{q}}_{\mathrm{n}}}{\displaystyle \mathrm{\sum }_{\mathrm{k}\mathrm{=}0}^{\mathrm{n}}{\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}}\mathrm{\left\{}{\mathrm{\lambda }}_2\mathrm{\left(}{\mathrm{\langle }{\mathrm{L}}_1(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{+}{\mathrm{\langle }{\mathrm{L}}_2(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right)}\mathrm{+}{\mathrm{\delta }}_{1\mathrm{,}\mathrm{m}}{\mathrm{\lambda }}_3{\displaystyle \frac{2}{{\mathrm{T}}_2}}\mathrm{\left(}{\mathrm{\langle }{\mathrm{L}}_3(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{+}{\mathrm{\langle }{\mathrm{L}}_4(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right)}\mathrm{\right\}}$$

(21)

we have defined:

$${\mathrm{\langle }{\mathrm{L}}_1(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_1^{\mathrm{*}}{\mathrm{D}}_2^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k}\mathrm{+}1)}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2(\mathrm{k}\mathrm{+}1)}}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{;}{\mathrm{\langle }{\mathrm{L}}_2(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_2^{\mathrm{*}}{\mathrm{D}}_3}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2(\mathrm{k}\mathrm{+}1)}{\mathrm{\left|}{\mathrm{D}}_3\mathrm{\right|}}^2}}\mathrm{\rangle }}_{\mathrm{\xi }}$$

$${\mathrm{\langle }{\mathrm{L}}_3(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_1^{\mathrm{*}}{\mathrm{D}}_2^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k}\mathrm{+}2)}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2(\mathrm{k}\mathrm{+}1)}}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{;}{\mathrm{\langle }{\mathrm{L}}_4(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{1}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k}\mathrm{+}1)}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2(\mathrm{k}\mathrm{+}1)}}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{.}$$

For the frequency-coupled component of the FWM signal (${\mathrm{\omega }}_2$), we have

$${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_3)\mathrm{=}{\displaystyle \mathrm{\sum }_{\mathrm{n}\mathrm{=}0}^{\mathrm{\infty }}{\mathrm{q}}_{\mathrm{n}}}{\displaystyle \mathrm{\sum }_{\mathrm{k}\mathrm{=}0}^{\mathrm{n}}{\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}}\mathrm{\left\{}\mathrm{-}{\mathrm{\lambda }}_2^{\mathrm{*}}\mathrm{\left(}{\mathrm{\langle }{\mathrm{M}}_1(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{+}{\mathrm{\langle }{\mathrm{M}}_2(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{\right)}\mathrm{\right\}}$$

(22)

Considering the following functions:

$${\mathrm{\langle }{\mathrm{M}}_1(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_1{\mathrm{D}}_3^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k}\mathrm{+}1)}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2\mathrm{k}}{\mathrm{\left|}{\mathrm{D}}_3\mathrm{\right|}}^2}}\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{;}{\mathrm{\langle }{\mathrm{M}}_2(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}\mathrm{=}{\mathrm{\langle }{\displaystyle \frac{{\mathrm{D}}_2{\mathrm{D}}_3^{\mathrm{*}}}{{\mathrm{\left|}{\mathrm{D}}_1\mathrm{\right|}}^{2(\mathrm{n}\mathrm{-}\mathrm{k})}{\mathrm{\left|}{\mathrm{D}}_2\mathrm{\right|}}^{2(\mathrm{k}\mathrm{+}1)}{\mathrm{\left|}{\mathrm{D}}_3\mathrm{\right|}}^2}}\mathrm{\rangle }}_{\mathrm{\xi }}$$

where the J-function is defined as

$$\mathrm{J}\mathrm{=}{\displaystyle \frac{1}{{\mathrm{T}}_1}}\mathrm{\left(}1\mathrm{+}{\displaystyle \frac{4{\mathrm{S}}_1}{{\mathrm{T}}_2^2{\mathrm{D}}_3{\mathrm{D}}_2^{\mathrm{*}}}}\mathrm{\right)}\mathrm{-}\mathrm{i}\mathrm{\Delta }\mathrm{\left(}1\mathrm{-}{\displaystyle \frac{4{\mathrm{S}}_1}{{\mathrm{T}}_1{\mathrm{T}}_2{\mathrm{D}}_3{\mathrm{D}}_2^{\mathrm{*}}}}\mathrm{\right)}\mathrm{+}2{\displaystyle \frac{{\mathrm{\left|}{\tilde{\mathrm{\Omega }}}_2\mathrm{\right|}}^2}{{\mathrm{D}}_1}}\mathrm{\approx }{\displaystyle \frac{1}{{\mathrm{T}}_1}}\mathrm{-}\mathrm{i}\mathrm{\Delta }$$

(23)

In this case, J can be interpreted as a Lorentzian width regulated by the longitudinal relaxation time, for the local and non-saturated cases of FWM signal. We have defined in the above expressions:

$${\mathrm{\lambda }}_1\mathrm{=}{\mathrm{a}}_1{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}\mathrm{.}}\mathrm{;}{\mathrm{\lambda }}_2\mathrm{=}\mathrm{-}{\displaystyle \frac{{\mathrm{a}}_2{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}\mathrm{.}}}{{\mathrm{J}}^{\mathrm{*}}}}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_1)\mathrm{\right|}}^2\mathrm{;}{\mathrm{\lambda }}_3\mathrm{=}{\displaystyle \frac{{\mathrm{a}}_3{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}\mathrm{.}}}{{\mathrm{J}}^{\mathrm{*}}}}{\mathrm{\left|}\overrightarrow{\mathrm{E}}({\mathrm{\omega }}_2)\mathrm{\right|}}^2\mathrm{.}$$

It is important to note that the different processes (I,C,coupling), ${\mathrm{q}}_{\mathrm{n}}$ strictly contains dissipative processes defining the type of molecular system through the ratio of the relaxation times (except in case of ${\mathrm{q}}_0\mathrm{=}1$); ${\mathrm{q}}_{\mathrm{k}\mathrm{,}\mathrm{n}}$defines the intensity of the pump- and probe-beam dispersion processes (except in case of ${\mathrm{q}}_{0\mathrm{,}0}\mathrm{=}1$); while ${\mathrm{\langle }{\mathrm{\Phi }}^{(\mathrm{j})}(\mathrm{n}\mathrm{,}\mathrm{k})\mathrm{\rangle }}_{\mathrm{\xi }}$, controls the stochastic nature $\mathrm{\xi }(\mathrm{t})\mathrm{=}{\tilde{\mathrm{\omega }}}_0\mathrm{+}\mathrm{\sigma }(\mathrm{t})$ of the process by the presence of the solvent and its strong interaction with the solute. It is important to note that in the functions ${\mathrm{\lambda }}_{\mathrm{n}}$ for $(\mathrm{n}\mathrm{=}1\mathrm{,}2\mathrm{,}3)$ reside the molecular structure effects under consideration of the vibronic coupling through the transition dipole moments of the states in the adiabatic basis ${\overrightarrow{\mathrm{\mu }}}_{\mathrm{B}\mathrm{A}}$; on the other site, the terms ${\mathrm{D}}_{\mathrm{n}}$ contain a Bohr frequency evaluated in the states with vibronic-type energies ${\mathrm{E}}_{\mathrm{A}}\mathrm{,}{\mathrm{E}}_{\mathrm{B}}$. Because of the asymmetric treatment of the pump beam at all orders, unlike the perturbative scheme for the weak beams has that ${\mathrm{\chi }}_{\mathrm{C}}^{(3)}({\mathrm{\omega }}_1)\mathrm{=}{\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}^{(3)}({\mathrm{\omega }}_1)\mathrm{=}0$. On the other hand, it is possible to note that, in the incoherent susceptibilities, independent of the optical frequency of calculation, the ${\mathrm{\lambda }}_1$ parameter has no dependence on the intensities of the beams. So far, the treatment we have given to the FWM signal is associated to the local regime, considering the quantum nature of the structure and the strong presence of the solvent through the function $\mathrm{\xi }(\mathrm{t})$. Taking the nonlinear polarization induced in the material as the source of propagation, and using Maxwell's equation [42], it is possible to derive expressions for the absorption coefficient and nonlinear refractive index, given by: ${\tilde{\mathrm{\alpha }}}_{\mathrm{j}}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})\mathrm{=}{\displaystyle \frac{2\mathrm{\pi }{\mathrm{\omega }}_{\mathrm{j}}}{{\tilde{\mathrm{\eta }}}_{\mathrm{j}}\mathrm{c}}}\mathrm{Im}{\mathrm{\chi }}_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})$ and ${\tilde{\mathrm{\eta }}}_{\mathrm{j}}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})\mathrm{=}{\mathrm{\left[}1\mathrm{+}4\mathrm{\pi }\mathrm{Re}{\mathrm{\chi }}_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})\mathrm{\right]}}^{1\mathrm{/}2}$, respectively, for ${\mathrm{\chi }}_{\mathrm{e}\mathrm{f}\mathrm{f}}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})\mathrm{=}{\mathrm{\chi }}_{\mathrm{I}}^{(1)}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})\mathrm{+}{\mathrm{\chi }}_{\mathrm{C}}^{(3)}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})$. In this case, ${\tilde{\mathrm{\alpha }}}_{\mathrm{j}}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})$ and ${\tilde{\mathrm{\eta }}}_{\mathrm{j}}({\mathrm{\omega }}_{\mathrm{j}}\mathrm{,}\mathrm{z})$ represent the quantities under study, which have dependencies both on the structural parameters that define the vibronic coupling both in the Rabi frequency and in the induced dipole moments, and on those parameters of the FWM signal generation process, mainly associated to the longitudinal and transverse relaxation times. In the present work, we will focus on reviewing the nonlinear optical absorption and scattering properties and the sensitivity of these in relation to different parameters of both the molecular structure and the solvent effect. The study of these nonlinear responses can be complemented with the developments of optical propagation and the effect of coherent and incoherent components on the susceptibility responses, proposed in Ref. [42].

3. Results and Discussion

3.1. Nonlinear Optical Properties

The FWM signal was characterized in terms of the behavior of its nonlinear optical properties and relative intensity, considering the stochasticity due to the presence of the solvent and the intramolecular nature of the coupling of the states. We keep the rotating wave approximation in the strictly near-resonance processes and select as a guide for the choice of some parameters, organic dyes such as green malachite chloride. We raise perturbatively the second order of the probe, which allows us to give presence to new photonic processes and how they can influence the specific signal of higher intensity at frequency ${\mathrm{\omega }}_3\mathrm{=}{\mathrm{\omega }}_1\mathrm{+}\mathrm{\Delta }$. To elucidate the behavior of the signal, we focus our attention on the different nonlinear susceptibilities that are generated, both absorptive (I) and saturative (C), as well as those generated by new coupling processes. We study how sensitive they respond to the optical length z, the intensity of the beams, frequency detuning of the pump and probe beams, time ratios of the longitudinal and transverse relaxation processes, and the parameters that specifically characterize the stochastic nature of the solute-solvent collision, and the vibronic effects resulting from the residual Hamiltonian. The latter is conditioned on the Bohr frequency of the adiabatic states ${\tilde{\mathrm{\omega }}}_0$ and the coupling strength $\mathrm{\nu }$ by the generation of the new transition dipoles (${\tilde{\mathrm{\Omega }}}_{\mathrm{k}}\mathrm{=}{\widetilde{\overrightarrow{\mathrm{\mu }}}}_{\mathrm{B}\mathrm{A}}\mathrm{.}{\overrightarrow{\mathrm{E}}}_{\mathrm{k}}\mathrm{/}\mathrm{\hslash }$). To reduce the complexity of the calculations, we have considered the permanent dipoles of the adiabatic basis to be null. The reason why we cannot determine the values at the maxima of the population oscillations is not because the permanent dipole moments of the adiabatic states are equal or because their energies are close. It is because we need to stay in the RWA and avoid including terms associated with non-resonant processes: $\mathrm{\pm }({\mathrm{\omega }}_1\mathrm{+}{\mathrm{\omega }}_2)$,$\mathrm{n}{\mathrm{\omega }}_1\mathrm{-}(\mathrm{n}\mathrm{-}1){\mathrm{\omega }}_2$ (for $\mathrm{n}\mathrm{\ge }3$ ) or $\mathrm{n}{\mathrm{\omega }}_2\mathrm{-}(\mathrm{n}\mathrm{-}1){\mathrm{\omega }}_1$ (for $\mathrm{n}\mathrm{\ge }2$ ). This assumption makes sense, since the higher-order components of ${\mathrm{\rho }}_{\mathrm{B}\mathrm{A}}({\mathrm{\omega }}_{\mathrm{k}})$must become less significant if the Fourier series in eq. (3) are the converge and the RWA is to remain valid. We selected an internal angle between the pump and probe incident beams, of $\mathrm{\theta }\mathrm{\approx }3\mathrm{°}$(0,0523599 radians) [14], given the very fast decay of the FWM signal intensity as the optical length increases. The parameters corresponding to a typical molecular system, i.e., the Malachite Green organic dye, were used to generate sensible realistic values for the calculations. For this well-known dye, the transition dipole moment in the adiabatic representation is ${\mathrm{\mu }}_{\mathrm{b}\mathrm{a}}\mathrm{=}2\mathrm{,}81\mathrm{\times }{10}^{\mathrm{-}18}\mathrm{e}\mathrm{r}{\mathrm{g}}^{1\mathrm{/}2}\mathrm{c}{\mathrm{m}}^{3\mathrm{/}2}$ and the Bohr frequency is ${\mathrm{\omega }}_0\mathrm{=}3\mathrm{,}06\mathrm{\times }{10}^{15}{\mathrm{s}}^{\mathrm{-}1}$[43]. In the same way, to establish a real comparison, we choose a defined optical length to study the nonlinear optical properties. We have centered this study on both the nonlinear optical properties and the intensity of the emerging signal from the FWM process. The profiles presented further obtained, obey an adiabatic representation of the intramolecular character at the crossing of harmonic curves of the considered states and an interaction with the solvent that shifts the natural Bohr frequency to a time-dependent and stochastic function.

3.1.1. Nonlinear Refraction Index

Figure 1 shows the general scattering profiles in the frequency space ${\mathrm{\Delta }}_1$and ${\mathrm{\Delta }}_2$ at different values of the transverse relaxation time. It can be noted that an increase in this parameter, in addition to changing the shape of the surface, causes the extreme values of this property to become increasingly distant from the refractive index of water (${\mathrm{\eta }}_0\mathrm{=}1.3333$), which is used as a solvent. We have assumed in these cases that all other variables such as saturation, propagation, longitudinal relaxation time, chemical concentration, vibronic coupling, angle of incidence between the beams, and coupling parameters with the thermal bath, remain constant.

Figure 2 shows cross sections of the nonlinear refractive index for different ratios between the frequencies ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$ at different times${\mathrm{T}}_2$, showing that the increase of this molecular parameter, in addition to shifting them vertically, also produces a horizontal shift of the maxima and minima towards the vertical axis, thus narrowing the dispersion curves. It is also worth mentioning that this effect is slightly larger for the case of ${\mathrm{\Delta }}_1\mathrm{=}2{\mathrm{\Delta }}_2$and considerably smaller when $2{\mathrm{\Delta }}_1\mathrm{=}{\mathrm{\Delta }}_2$. To evaluate the maximum contrast in the refractive index, we select the parameter where there is the highest maximum of the positive and negative sides ${\mathrm{\Delta }}_1$ and evaluate the quantity: ${\mathrm{\beta }}^{(\mathrm{\epsilon })}\mathrm{=}\mathrm{\Delta }{\mathrm{\eta }}_{{\mathrm{\Delta }}_1\mathrm{=}2{\mathrm{\Delta }}_2}\mathrm{/}\mathrm{\Delta }{\mathrm{\eta }}_{{\mathrm{\Delta }}_2\mathrm{=}2{\mathrm{\Delta }}_1}$, which $\mathrm{\epsilon }$ indicates the parameter being studied. Selecting ${\mathrm{T}}_2\mathrm{=}2.6\mathrm{p}\mathrm{s}$, the result is ${\mathrm{\beta }}^{({\mathrm{T}}_2)}\mathrm{=}4.7$. That is, the variations in the refractive index for this ${\mathrm{T}}_2$ are 4.7 times in the process ${\mathrm{\Delta }}_1\mathrm{=}2{\mathrm{\Delta }}_2$ than those reached in the same conditions when the process is ${\mathrm{\Delta }}_2\mathrm{=}2{\mathrm{\Delta }}_1$. In the latter case, the frequency of the generated signal beam coincides with the natural frequency of the system, so the absorption is considerably increased and hence the dispersion is reduced. In general, this turns out to be true for all subsequent analyses. In the processes ${\mathrm{\Delta }}_1\mathrm{=}{\mathrm{\Delta }}_2$and ${\mathrm{\Delta }}_1\mathrm{=}2{\mathrm{\Delta }}_2$ we observe rapid decay until the value of the refractive index of water is reached.

In Figure 3 the behavior of the dispersion curves can be corroborated by a continuous variation of the ${\mathrm{T}}_2$time. Since the transverse relaxation time is the time in which the molecular system loses the coherence of the induced electric dipoles; a longer delay in the loss of coherence implies a longer persistence of the induced polarization; which results in a higher scattering capacity of the beam. Figure 4 shows the overall scattering profiles in the frequency space ${\mathrm{\Delta }}_1$and ${\mathrm{\Delta }}_2$at different ratios between the longitudinal relaxation times ${\mathrm{T}}_1$ and transverse relaxation times ${\mathrm{T}}_2$. On the surfaces it can be noted that the most significant change is observed for when ${\mathrm{T}}_1\mathrm{=}10{\mathrm{T}}_2$, there is a noticeable difference along the plane ${\mathrm{\Delta }}_1\mathrm{=}{\mathrm{\Delta }}_2$ and its surroundings

The aforementioned can be corroborated by visualizing in Figure 5 the cross sections of the surface at the considered ratios between ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$: the maximum and minimum values reached by the dispersion differ considerably between the first and second case. In the resonance case, ${\mathrm{T}}_1\mathrm{=}10{\mathrm{T}}_2$ a very accelerated variation in the dispersion of the generated signal beam is observed when the detuning of the pumping and test beam are close to zero. This latter fact of an abrupt change in the refractive index at a minimum change in the pumping beam detuning is remarkable. However, it is possible to observe that it is fulfilled: ${\mathrm{\beta }}^{({\mathrm{T}}_1\mathrm{/}{\mathrm{T}}_2)}\mathrm{\approx }7$ in cases where ${\mathrm{T}}_1\mathrm{=}10{\mathrm{T}}_2$.

Figure 6 plots these changes by continuously varying the ratio between when ${\mathrm{T}}_1$ and ${\mathrm{T}}_2$. A higher ratio between these two relaxation times implies, as has been shown, an increase in signal beam dispersion. This longitudinal component is the time it takes for the molecular system to reach thermodynamic equilibrium once its excitation by the electromagnetic field has stopped, so if this time is longer ${\mathrm{T}}_2$ (as is normally the case) the system will take longer to reach equilibrium and may scatter more light. It is possible to observe that the refractive index shows a higher contrast just at the resonance process ${\mathrm{\omega }}_3$ with the Bohr frequency of the coupled states (${\tilde{\mathrm{\omega }}}_0$). Although the intensity differences are not very high, however, this fact of abrupt change can be used as a sensor generation mechanism, since it changes rapidly with small variations in the pumping detuning.

Figure 7 shows the general surfaces of the nonlinear refractive index in the frequency space ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$ at different coupling factors between the molecular system under study and the thermal reservoir used.

As shown in Figure 7, Figure 8 and Figure 9, a decrease in this coupling factor $\mathrm{\sigma }$ increases the light scattering capacity. The coupling factor between the molecular system and the thermal reservoir is a measure of the interaction between these two entities, so a lower interaction would imply a lower number of collisions and, therefore, a longer time to reach thermodynamic equilibrium. The presence of the thermal bath in its collision with the potentials of the active system, acts by modifying the minima vertically without modifying the nuclear coordinates, leading to a randomness in the process and a lower dispersion, a fact that is reflected mainly in the resonance process of ${\mathrm{\omega }}_3$ with ${\tilde{\mathrm{\omega }}}_0$. It is important to note that for values $\mathrm{\sigma }{\mathrm{T}}_1^{\mathrm{-}1}\mathrm{=}0.7$ it is necessary that ${\mathrm{\beta }}^{(\mathrm{\sigma }{\mathrm{T}}_1)}\mathrm{\approx }2.3$.

Figure 10 shows the general dispersion profiles in the frequency space ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$at different values of the saturation parameter. In these profiles and those shown in Figure 11 and Figure 12, the effect of saturation in the system is to increase the dispersion of the signal beam, keeping practically constant the values at which the maximum and minimum of this property of the system are reached. As saturation increases, the population of the upper-level also increases, so there will be less molecules in the fundamental state, thus limiting absorption and, as less is absorbed, the beam generated will be more dispersed.

The saturation values were selected in experimental correlation and conditioning their magnitudes for the fulfilment of the approximations carried in the expansion of the zero-frequency components of the population difference. We note in this case that due to the saturation effect, we have ${\mathrm{\beta }}^{({\mathrm{S}}_{12})}\mathrm{\approx }2.6$

Figure 13 shows the overall scattering profiles in the frequency space ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$ at different values of the intramolecular coupling parameter ν. The maxima and minima of the surfaces start to be closer to the refractive index of water as the value of this parameter is increased. In the cross sections for different ratios between ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$ at different values of ν in Figure 14, the effect of intramolecular coupling is better visualized, which decreases the optical property in question as it increases. However, in this case ${\mathrm{\beta }}^{(\mathrm{\nu })}\mathrm{\approx }4$ when the lowest coupling value is taken. In cases where the intramolecular effect is considerably increased, $\mathrm{\beta }\mathrm{\approx }1$ is taken.

Initially, the effect is not very noticeable between $\mathrm{\nu }\mathrm{=}0.01$and $\mathrm{\nu }\mathrm{=}0.1$, but then it ends up becoming more visible as the coupling continues to increase, as shown in Figure 15. From the expression for ${\tilde{\mathrm{\mu }}}_{\mathrm{A}\mathrm{B}}$ (and thus for ${\tilde{\mathrm{\mu }}}_{\mathrm{B}\mathrm{A}}$) it can be found that an increase in the intramolecular coupling factor decreases the value of this transition dipole in the adiabatic basis, which decreases the electrical susceptibility of the system and thus its ability to polarize and orient its dipoles, so it ends up scattering less light.

In general, since the detuning of the incident beams is exactly zero, the refractive index turns out to be that of the solvent. In addition, it can be also noted that this property presents an inversion symmetry to an origin centered at the point (0,1.333) of the coordinate plane. In summary, we can indicate that the refractive index variations taken in the ratio of the ${\mathrm{\Delta }}_1\mathrm{=}2{\mathrm{\Delta }}_2$ and ${\mathrm{\Delta }}_2\mathrm{=}2{\mathrm{\Delta }}_1$ processes are most sensitive to the relaxation time ratio ${\mathrm{T}}_1\mathrm{=}10{\mathrm{T}}_2$, with similar effects when varying the effects of coupling with the thermal bath $\mathrm{\sigma }{\mathrm{T}}_1^{\mathrm{-}1}$ and beam saturation ${\mathrm{S}}_{12}$, leaving very little sensitivity to the effects of intramolecular coupling $\mathrm{\nu }$. From the information obtained from the parameter $\mathrm{\epsilon }$ we observe the highest sensitivity with the relaxation time ratio ${\mathrm{T}}_1\mathrm{/}{\mathrm{T}}_2$. This denotes that the optical property raises its intensity contrast symmetrically in relation to detuning, parameterized by the relaxation kinetics at thermal equilibrium with respect to the loss of coherence between states. This fact correlates with that proposed by Yajima et al. [44,45] when analyzing the response of organic dyes in Rayleigh-type optical mixing processes.

3.1.2. Nonlinear Absorption Coefficient

Figure 16, Figure 17 and Figure 18 show different absorption profiles as the transverse relaxation time ${\mathrm{T}}_2$ is varied. It is observed that an increase in the transverse relaxation time implies greater absorption since its increase reflects the dissipation mechanism in which the dipoles take longer to lose their coherence, which implies that a greater matter-field interaction is maintained, and with-it greater absorption.

As can be noticed, compared to the other cases, the absorption ${\mathrm{\omega }}_3\mathrm{=}{\tilde{\mathrm{\omega }}}_0$ is significantly higher than the other cases of the ratio between ${\mathrm{\Delta }}_1$ and ${\mathrm{\Delta }}_2$.

Also, for this case the absorption is never negative, so one would not expect to observe some kind of parametric amplification in this frequency ratio when studying the spatial variation of the generated signal intensity. We also observe that in this resonance, the response presents a greater width for a higher pumping detuning frequency range, which we can also predict from the HWHM associated with the inverse of the transverse relaxation time.

Figure 19, Figure 20 and Figure 21 show the change in the absorption profiles due to the increase in the ratio of longitudinal and transverse relaxation times ${\mathrm{T}}_1\mathrm{/}{\mathrm{T}}_2$. For the case where ${\mathrm{T}}_1\mathrm{=}10{\mathrm{T}}_2$we observe a significant distortion in the absorption profile, with apparent zones of parametric amplification around ${\mathrm{\Delta }}_1{\mathrm{T}}_2\mathrm{\approx }1$.

We can appreciate the resonances of the degenerate processes ${\mathrm{\Delta }}_1\mathrm{=}{\mathrm{\Delta }}_2$and $2{\mathrm{\omega }}_2\mathrm{-}{\mathrm{\omega }}_1\mathrm{=}{\tilde{\mathrm{\omega }}}_0$ that in time ratios of ${\mathrm{T}}_1\mathrm{=}10{\mathrm{T}}_2$ the optical response shows a great broadening in contrast to the process $2{\mathrm{\omega }}_1\mathrm{-}{\mathrm{\omega }}_2\mathrm{=}{\tilde{\mathrm{\omega }}}_0$ where the profile is maintained with a very reduced width. However, the apparent zones of amplification in the first two profiles are fundamentally accentuated at these time ratios, totally different in the cases of resonance of the FWM signal where although there is very high broadening because of the decrease of the time ratio and the HWMH, no parametric amplifications (negative absorptions) are shown throughout the detuning zone of the pumping.

In agreement with the previous result, a longer time to reach equilibrium will lead to an intensification in light absorption, which is observed in the profiles of Figure 22, Figure 23 and Figure 24. In this case, the longer time to reach equilibrium is given by the lower coupling factor between the system of interest and the thermal reservoir through a lower number of intermolecular collisions.

We also see that in the ${\mathrm{\omega }}_1\mathrm{=}{\mathrm{\omega }}_2$ resonance the effect is very similar to the one that occurs when varying the time ratio in this profile. The two processes of resonance degenerate and are tuned to $2{\mathrm{\omega }}_2\mathrm{-}{\mathrm{\omega }}_1\mathrm{=}{\tilde{\mathrm{\omega }}}_0$ show negative absorption zones mostly accentuated in weak couplings, however, in the resonance $2{\mathrm{\omega }}_1\mathrm{-}{\mathrm{\omega }}_2\mathrm{=}{\tilde{\mathrm{\omega }}}_0$, the spectrum tends to be extinguished with the greater increase of the strong coupling with the solvent. It is evident that, for higher coupling values, the randomness increases and the movements of the potential minima of the two harmonic curves make the transition less effective and considerably decrease the absorption of the system, producing extreme broadening and decreased absorption in the MPO zone.

The absorption profiles of the signal beam can be seen in Figure 25, Figure 26 and Figure 27. In general, increasing saturation will mean that the upper level is more populated and there is less room for an absorption transition since the population of the ground state will be smaller, which can be seen in the profiles up to a certain detuning value.

However, as the detuning value is closer to zero, the absorption will begin to grow even when the saturation is also growing, since as the electromagnetic wave has a frequency closer to the transition frequency, even when the population of the basal level is lower, the probability of transition and the absorption increase.

The effect of intramolecular coupling on the absorption profiles is shown in Figure 28, Figure 29 and Figure 30. Following a similar reasoning as in the scattering case, the absorption coefficient will increase as the intramolecular coupling decreases, due to the effect of the latter on the transition dipole moment.

In all cases, the absorption reaches its maximum value when all detuning is equal to zero. It can also be noted that the absorption coefficient presents some symmetry to ${\mathrm{\Delta }}_1$, however, because the refractive index and the exact symmetry is not reached. In summary terms, we can represent both profiles for different points in the frequency space, according to Figure 31.

In general, as would have been appreciated earlier, the optical properties do not exhibit symmetry in terms of the exchange between ${\mathrm{\Delta }}_1$and ${\mathrm{\Delta }}_2$. For the refractive index, it can be noted that the more the ratio ${\mathrm{\Delta }}_1\mathrm{/}{\mathrm{\Delta }}_2$ increases, the more the maxima and minima shift more towards the vertical axis, also reaching a higher and lower value, respectively; it can also be observed that the scattering decays faster as the ratio increases. It follows naturally that the radiation scatters to a greater extent as the resulting frequency moves farther and farther away from the transition frequency of the system. In addition, the detuning of the pumping beam should be smaller and smaller so that the resulting frequency of the signal does not differ so much from that of the system and these can interact. Another important detail that can be noticed is that for the case ${\mathrm{\Delta }}_2\mathrm{=}3{\mathrm{\Delta }}_1$the maximum and minimum of dispersion change position with respect to the vertical axis, this because the detuning ${\mathrm{\Delta }}_3$ turns out to be the negative of ${\mathrm{\Delta }}_1$. As for the absorption coefficient, starting from ${\mathrm{\Delta }}_3\mathrm{=}0$, any increase or decrease in the ratio ${\mathrm{\Delta }}_1\mathrm{/}{\mathrm{\Delta }}_2$ will imply a horizontal contraction in the absorption curve. In these cases, as the signal frequency is farther away from the natural frequency, ${\mathrm{\Delta }}_1$would need to be decreased to decrease the frequency difference so that the beam interacts more with the system and an increase or decrease in the radiation intensity is observed.

3.2. Perturbative Influence on the Development of the Zero-Frequency Component

The participation of the probe beam in the FWM generation process to second order in a perturbative manner (eq. 6) is also evaluated. Considering the incident probe beam at this order, not only generates new dynamical processes of importance included in the nonlinear response optical susceptibilities, but also constitutes a difference in treatment to the emerging signal beam treated to the first order. Treating the probe beam to the second order does not necessarily imply losing or changing its role in the signal generation process as an observer of the perturbing effect of the intense pump beam but allows the test beam a recognition in the process not necessarily linked to its intensity.

On the other hand, our stochastic considerations of the solvent fundamentally reflected in all the averages of susceptibilities in the set of realizations of the random variable, have fundamentally their origin in the behavior of the component at zero frequencies of the difference of populations ${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}$between the adiabatic states (eq. 15). To obtain a useful expression of this component, a series development is performed assuming approximation limits subject to the saturation generated by the two beams, but mainly the pumping beam (eq. 16). To carry out convergence studies of such series is extremely complex and not very useful because of the calculation times involved in raising to high orders of perturbation and making representations of the optical properties in the frequency space parameterized by different study variables. Accordingly, we study the effect that the orders have on both the nonlinear refractive index and the nonlinear absorption coefficient. We observe in Figure 32 that for the degenerate and $2{\mathrm{\omega }}_2\mathrm{-}{\mathrm{\omega }}_1\mathrm{=}{\tilde{\mathrm{\omega }}}_0$ tuning processes, we could reach only the first order or even stay at zero order of expansion for our calculations and obtain acceptable results, unlike the FWM process $2{\mathrm{\omega }}_1\mathrm{-}{\mathrm{\omega }}_2\mathrm{=}{\tilde{\mathrm{\omega }}}_0$, where it is strictly necessary to resort at least to the first order to obtain acceptable results in the expansion of the population difference at zero frequencies.

However, for evaluations in the absorption coefficient, we can observe from Figure 33, that almost invariantly in the degenerate and resonance processes of the test beam with the Bohr frequency in the adiabatic states the zero-order in the perturbative development is sufficient to generate reliable data. This indicates that the system is in an equilibrium situation established according to ${\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{d}\mathrm{c}}\mathrm{\approx }{\mathrm{\rho }}_{\mathrm{D}}^{\mathrm{e}\mathrm{q}}$. However, in the process awarded to the FWM signal for the ${\mathrm{\omega }}_3\mathrm{=}{\tilde{\mathrm{\omega }}}_0$ resonance, it is even necessary to scale up to the second order of expansion in the evaluation of the population difference at zero frequencies. This indicates the possibility that with the saturation effects induced by the orders considered, the system remains in situations far out of equilibrium and that the relaxation times ${\mathrm{T}}_1$ used lead the system to leave equilibrium and probably return to other canonical equilibrium situations, thus losing the memory of the original equilibrium value. We have focused mainly on the behavior of ${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}({\mathrm{\omega }}_3\mathrm{,}\mathrm{z})$ and ${\mathrm{\chi }}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{p}}({\mathrm{\omega }}_2\mathrm{,}\mathrm{z})$ as a function of the longitudinal and transverse relaxation times since they are fundamental parameters associated with two dissipation processes intimately linked in the stochastic optical Bloch equations and subject to vibronic effects, but which clearly distinguish different molecular systems of study, preferably organic dyes.

4. Concluding Remarks

We have devised a semiclassical radiation-matter interaction model. This model considers a two-state system and a molecule model based on coupled harmonic potentials with energy and nuclear coordinate shifts. We make use of the Liouvillian scheme for the reduced density matrix with the use of quantum relaxation theory for dissipation processes referring to molecular populations separated from the equilibrium situation, and those where coherence aspects are measured. Furthermore, we consider for the presence of the solvent in a scheme that shifts the Bohr frequency to a stochastic time-dependent function. To simulate the FWM signal, we align our model with experimental facts to represent the high-intensity pumping beam in saturated form and incorporate a low-intensity probe beam according to a perturbative development, which opens up the possibility of new photonic processes linked to the FWM process under study. We have considered as critical variables in the study those derived from the molecular structure, as well as those derived from the optical process in the analysis of the nonlinear absorptive and dispersive properties. As a result, we have found the wide possibilities of study that are opened in the representations of the optical properties in the frequency space and its dependence on the saturation that induce the beams or the duration in time of the relaxation processes. We analyze the refractive indices with possibilities of designs as tools in optical sensors for the possibility of a high change or contrast according to small variations in the frequency detuning of the beams. Dependence of the optical profile on its reduction with increasing intramolecular parameter or sensitivity with the order of test perturbation considered. This study allows us to understand how important is the inclusion of high-order perturbation from low-intensity fields to warrant the measurement of the refractive index or nonlinear absorption coefficient of the employed organic dyes. Both the refractive index and the nonlinear absorption coefficient are affected to a greater or lesser extent with the variation of the molecular and field parameters considered. These variations cause a greater or lesser effect depending on the ratio between the frequencies of the incident electromagnetic fields, finding that the effect generated is significantly lower in the scattering when the frequency of the resulting signal beam is identical to the Bohr frequency of the molecular system under study. Regarding the absorption coefficient for the frequency ratio $2{\mathrm{\Delta }}_1\mathrm{=}{\mathrm{\Delta }}_2$, as expected, the absorption is maximized compared to other frequency ratios considered, coinciding only the absorptions in all profiles when the detuning of the incident beams is zero. In the developed model it can also be seen that the dispersion increases the greater the difference between the resulting frequency of the signal and that of the system, while the absorption decreases as this difference increases. The consideration of a second order test beam would not represent a significant contribution to the predicted results, since its percentage contribution is negligible within the model development framework, however, this contribution starts to become relevant mainly with the increase of the ratio between longitudinal and transverse relaxation times.