## 2 The Model
The modified \(\rho_{nm}^{(2)}\) at presence of weak probe field is  
\[\rho_{nm}^{(2)^{\prime}}=\frac{-i}{\hbar}\int_{-\infty}^{t}[\hat{V}_{c}(t^{ \prime})+\hat{V}_{p}(t^{\prime}),\rho^{(1)}]_{nm}e^{(i\omega_{nm}+\gamma_{nm})(t ^{\prime}-t)}dt^{\prime}, \tag{8}\]  
the difference \(\delta\rho_{nm}^{(2)}=\rho_{nm}^{(2)^{\prime}}-\rho_{nm}^{(2)}\) brought up by the insertion of weak probe field \(\hat{V}_{p}\) is  
\[\delta\rho_{nm}^{(2)}=\frac{-i}{\hbar}\int_{-\infty}^{t}[\hat{V}_{p}(t^{\prime }),\rho^{(1)}]_{nm}e^{(i\omega_{nm}+\gamma_{nm})(t^{\prime}-t)}dt^{\prime}, \tag{9}\]  
which is the \(\lambda^{3}\) order correction to \(\rho_{nm}^{(2)}\). This corresponds to the susceptibility induced by weak probe field in the presence of strong coupling fields. To figure out the expression of \(\delta\rho_{nm}^{(2)}\), we first calculate the first order correction \(\rho_{nm}^{(1)}\). For the strong coupling fields, we suppose \(\hat{V}_{c}(t^{\prime})\) is in the form  
\[\hat{V}_{c}(t^{\prime})=\sum_{nm}\hat{V}_{c(nm)}(t^{\prime})=-\sum_{nm}\!\! \mathbf{d}_{nm}\cdot\mathbf{E}_{c(nm)}(t^{\prime}), \tag{10}\]  
with  
\[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]0  
being the field at frequency \(\omega_{f(nm)}\) coupling atomic states \(|n\rangle\) and \(|m\rangle\) with relevant dipole moment \(\mathbf{d}_{nm}\). Under the assumption that the unperturbed density matrix element is given by \(\rho_{nm}^{(0)}=0\) for \(n\
eq m\), we have  
\[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]1  
and we therefore can write down  
\[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]2  
We then assume the weak probe field is applied between the atomic states \(|k\rangle\) and \(|l\rangle\), so that  
\[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]3  
where \(\mathbf{E}_{probe}\) is the amplitude of the probe field, \(\mathbf{d}_{kl}\) is the relevant dipole moment between the atom's energy levels \(|k\rangle\) and \(|l\rangle\), and \(\omega_{f(kl)}\) is the frequency of the probe field. Then  
\[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]4 \[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]5 \[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]6  
Since the above expression is analogous to the one in Equation (12), we can easily deduce from analogy between Equations (5) and (6) and comparison to Equation (13) that  
\[\rho_{nm}=\rho_{nm}^{(0)}+\lambda\rho_{nm}^{(1)}+\lambda^{2}\rho_{nm}^{(2)}. \tag{2}\]7  
and this is the \(\delta\rho_{nm}^{(2)}\) in Equation (9) we are looking for. From the above results, we can get the expression of probe susceptibility. Normally, the expectation value of the probe field-related induced dipole moment \(\langle\widetilde{\mathbf{d}}(t)\rangle=tr(\hat{\rho}^{(2)}_{probe}\hat{\mathbf{d}})\) should be decomposed into sum of different frequency components as \(\langle\widetilde{\mathbf{d}}(t)\rangle=\sum\limits_{p}\langle\mathbf{d}( \omega_{p})\rangle e^{-i\omega_{p}t}\). Since in our case there is only one frequency \(\omega_{f(kl)}\) existing in the probe field, we have  
\[ \begin{array}{rcl}\langle\widetilde{\mathbf{d}}(t)\rangle&=&tr(\hat{\rho}^{ (2)}_{probe}\hat{\mathbf{d}})\\ &=&\rho^{(2)}_{kl}\mathbf{d}_{lk}\\ &=&\frac{i}{\hbar}(\rho^{(1)}_{ll}-\rho^{(1)}_{kk})\frac{\mathbf{d}_{lk}[ \mathbf{d}_{kl}\cdot\mathbf{E}_{probe}(\omega_{f(kl)})]e^{-i\omega_{f(kl)}t}}{ i(\omega_{kl}-\omega_{f(kl)})+\gamma_{kl}}.\end{array}  \tag{17}\]  
On the other hand the probe susceptibility can be defined in equation  
\[P(\omega_{p})=\mathcal{N}\langle\mathbf{d}(\omega_{p})\rangle=\mathcal{N} \langle\widetilde{\mathbf{d}}(t)\rangle=\epsilon_{0}\chi^{(2)}_{kl}(\omega_{ p})\cdot\mathbf{E}_{probe}(\omega_{p}), \tag{18}\]  
where \(\omega_{p}=\omega_{f(kl)}\) and \(\mathcal{N}\) is the number density of atoms.