Enhanced nonlinearity of four-wave mixing via Rydberg–Rydberg interactions

Atoms excited to high-lying Rydberg states interact via strong and long-range dipole–dipole or van der Waals forces, and such strong interactions manifest as a well-known excitation blockade [1] in which an atom promoted to a Rydberg level shifts the energy levels of nearby atoms and suppresses their excitation. The blockade is of great interest for applications in quantum information [2] and provides the basis for numerous proposals, such as the photon–photon phase gates [3], atomic logical gates [4, 5], single-photon switch [6] and generating nonclassical state of light [7, 8].

Due to such pronounced nonlocal interaction between the atoms, the Rydberg atomic gases are regarded as a promising candidate for many-body correlated systems which act as effective tools to control the interactions between photons [9, 10]. One of such applications is to enhance nonlinearity [11–13] with the help of electromagnetically induced transparency (EIT). For example, the Kerr effect, as one of main topics in nonlinear optics with numerous applications in, such as single-photon switches and transistors [9], controlled quantum gates [14–16], quantum teleportation [17] and entanglement [18], can be greatly enhanced in Rydberg-EIT system [19–21]. Calculations beyond mean-field theory [21–25] attribute such improvement to the two-body or even three-body correlators based on perturbation method.

Motivated by the recent investigations on strengthening the Kerr nonlinearity, in this paper we study the enhancement of another third-order nonlinearity, namely, the four-wave mixing (FWM) process which has been used as an important resource for many quantum applications, such as quantum entanglement and steering [26–30], quantum cloning [31], generating correlated beams and photon pairs [32–35], constructing nonlinear interferometer [36] and quantum networks [37], as well as realizing optical nonreciprocity [38], etc.

The system that we are interested in is shown in figure 1(a) where three applied fields, specifically, a pumping field (Ω_p), coupling field (Ω_c) and a driving field (Ω_d) illuminate the four-level atoms with a Rydberg state as the highest level to generate a signal (Ω_s) with the efficiency determined by the third-order susceptibility. Investigations using the similar excitation schemes are reported, for example, to build a single-photon source by utilizing the Rydberg blockage of confined atoms in an excitation volume similar to, or smaller than that of the blockade sphere [39–41], or using the pumping and coupling fields to construct Rydberg dark-state polaritons for, e.g. the efficient light storage which can be retrieved by applying the driving field after a short time period [42, 43]. We also note the early experimental investigations on FWM in Rydberg gases of which the atomic density is too low to invoke significant Rydberg–Rydberg interaction [44–47]. In this paper, we study instead a (virtual) transition |1⟩ → |2⟩ → |4⟩ → |3⟩ → |1⟩ in an ensemble of cold Rydberg atoms of which the size is larger than the blockade sphere as schematically shown in figure 1(b), to show that the corresponding nonlinear susceptibility is significantly enhanced by the Rydberg–Rydberg interaction.

Zoom In Zoom Out Reset image size

At first glance, the nonlinear susceptibility of the FWM process has a small modulus, meaning that the generation of Ω_s cannot be efficient, as the susceptibility is proportional to the product of the four relevant dipole-moment elements in which the ones associating with the transitions of Ω_c and Ω_d are much smaller than the other two. However, our calculations show that this is only true for the local interaction between the atoms and photons. When the Rydberg–Rydberg interation is considered, we need to evaluate the two-body correlator and take into account the contributions from all Rydberg atoms in the atomic gases. This leads to an increment in nonlinearity which is represented by a large nonlocal susceptibility.

The organization of this article is as follows: in section 2, we present the model and equations for the one-body and two-body correlators. Section 3 discusses the perturbation method we used to solve the equations. The properties of the local and nonlocal susceptibilities obtained from the solutions are analyzed as well. In appendix A, we present the detailed equations for the nonlocal susceptibility.

Let us consider an ensemble of ultra-cold (for example, ⁸⁷Rb) atoms each having a four-level double-ladder configuration as shown in figure 1(a) is loaded in a magneto-optical trap where the Doppler broadening caused by the thermal motion of atoms is significantly suppressed. Assuming that a strong coupling field (at angular frequency ω_c, with Rabi frequency Ω_c) drives the transition |2⟩ ↔ |4⟩ and a weak pumping field (ω_p, Ω_p) drives the transition between |1⟩ and |2⟩, then due to the presence of a strong driving field (ω_d, Ω_d), effective polarization could be built between the transition |1⟩ ↔ |3⟩, and subsequently resulting in the generation of a signal field (ω_s, Ω_s) by means of FWM process. The efficiency of the generation depends on the amplitude of (third-order) FWM nonlinear coefficient χ⁽³⁾ which is the key subject that we are interested in this paper.

The excited state |4⟩ is a Rydberg state and assumed to be |60 S_1/2⟩. The long-range interaction [48] between the Rydberg atom at position r and the one at position r' (see, figure 1(b)) is described by the potential V(Δr) = −C₆/Δr⁶ with Δr = |r − r'|.

We assume that the electric fields can be written as ${E}_{i}{=\mathcal{E}}_{i}\,{\text{e}}^{\text{i}({\mathbf{k}}_{i}\cdot \mathbf{r}-{\omega }_{i}t)}+\mathrm{c}.\mathrm{c}.,\enspace i=\left\{p,c,d,s\right\}$ with c.c. standing for complex conjugate. Then the corresponding Rabi frequencies are ${{\Omega}}_{i}={\mu }_{mn}{\mathcal{E}}_{i}/2\hslash$, where μ_mn with m, n = {1, 2, 3, 4} being the dipole moment of the transition |m⟩ ↔ |n⟩ that Ω_i drives. And the detunings are defined as Δ_p = ω₂₁ − ω_p, Δ_c = ω₄₁ − ω₂₁ − ω_c, Δ_d = ω₄₁ − ω₃₁ − ω_d. Δ_s = ω₃₁ − ω_s.

Adopting the similar notations used in reference [24], we use ${\hat{S}}_{mn}(\mathbf{r},t)=\vert n(\mathbf{r},t)\rangle \langle m(\mathbf{r},t)\vert$ to denote the the transition operators for the atom at r, and it satisfies the commutation relation $[{\hat{S}}_{ab}(\mathbf{r},t),{\hat{S}}_{mn}({\mathbf{r}}^{\prime },t)]=[{\delta }_{an}{\hat{S}}_{mb}(\mathbf{r},t)-{\delta }_{mb}{\hat{S}}_{an}({\mathbf{r}}^{\prime },t)]{\delta }_{\mathbf{r}{\mathbf{r}}^{\prime }}$ where δ_ab is Kronecker delta symbol. In the derivation to follow, we use ${\hat{S}}_{mn}(t)$ to represent the transition operator ${\hat{S}}_{mn}(\mathbf{0},t)$ for the atom at origin, then under the electric-dipole and rotating-wave approximations, the Hamiltonian depicting the interaction between the atom (at origin) and the four fields takes the following form:

$$\begin{equation}\begin{aligned}\hfill \hat{H}(t)/\hslash & ={{\Delta}}_{\mathrm{p}}{\hat{S}}_{22}(t)+{{\Delta}}_{\mathrm{s}}{\hat{S}}_{33}(t)+({{\Delta}}_{\mathrm{p}}+{{\Delta}}_{\mathrm{c}}){\hat{S}}_{44}(t)\hfill \\ \hfill & \quad +\left[{{\Omega}}_{\mathrm{p}}{\hat{S}}_{12}(t)+{{\Omega}}_{\mathrm{c}}{\hat{S}}_{24}(t)+{{\Omega}}_{\mathrm{d}}{\hat{S}}_{34}(t)+{{\Omega}}_{\mathrm{s}}{\hat{S}}_{13}(t)\right.\hfill \\ \hfill & \quad \left.+\text{H.c.}\right]+{N}_{0}\int {\mathrm{d}}^{3}r{\hat{S}}_{44}(\mathbf{r},t)V(r){\hat{S}}_{44}(t).\hfill \end{aligned}\end{equation} \tag{ 1 }$$

Here N₀ is the atomic density and the symbol ∫d³ r stands for the integration ${\int }_{0}^{\pi }\mathrm{sin}\,\theta \,\mathrm{d}\theta {\int }_{0}^{2\pi }\mathrm{d}\phi {\int }_{{R}_{\mathrm{b}}}^{\infty }{r}^{2}\,\mathrm{d}r$ in spherical coordinates with R_b being the radius of the blockage sphere (see figure 1(b)). Then the features of the generation, absorption, as well as dispersion of the signal Ω_s all hide in $\langle {\hat{S}}_{31}(t)\rangle ={\rho }_{31}(t)$ which is one of the one-body density-matrix elements governed by $\mathrm{i}\hslash \partial \langle {\hat{S}}_{mn}(t)\rangle /\partial t=\langle [{\hat{S}}_{mn}(t),\hat{H}]\rangle$. For instance, the equation of ρ₃₁(t) reads

$$\begin{equation}\frac{\partial {\rho }_{31}}{\partial t}=-{g}_{31}{\rho }_{31}+\mathrm{i}{{\Omega}}_{\mathrm{s}}({\rho }_{11}-{\rho }_{33})-\mathrm{i}{\rho }_{32}{{\Omega}}_{\mathrm{p}}+\mathrm{i}{\rho }_{41}{{\Omega}}_{\mathrm{d}}^{\ast };\end{equation} \tag{ 2a }$$

it further relates to other elements, e.g. ρ₄₁ whose time dependence is:

$$\begin{equation}\begin{aligned}\hfill \frac{\partial {\rho }_{41}}{\partial t}=& -{g}_{41}{\rho }_{41}+\mathrm{i}{\rho }_{21}{{\Omega}}_{\mathrm{c}}+\mathrm{i}{\rho }_{31}{{\Omega}}_{\mathrm{d}}-\mathrm{i}{\rho }_{42}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill & -\mathrm{i}{\rho }_{43}{{\Omega}}_{\mathrm{s}}-\mathrm{i}{N}_{0}\int {\mathrm{d}}^{3}\mathbf{r}V(r){\rho }_{44,41}(\mathbf{r},t).\hfill \end{aligned}\end{equation} \tag{ 2b }$$

Here g₃₁ = γ₃₁ + iΔ_s and g₄₁ = γ₄₁ + i(Δ_c + iΔ_p). The whole (closed) set of the equations for the one-body density-matrix element is listed in the appendix A. The last term on the right-hand side of equation (2b) comes from the Rydberg–Rydberg interation where ${\rho }_{44,41}=\langle {\hat{S}}_{44}(\mathbf{r},t){\hat{S}}_{41}(t)\rangle$ is one of the two-body correlators, generally defined as ${\rho }_{ab,mn}(\mathbf{r},{\mathbf{r}}^{\prime },t)=\langle {\hat{S}}_{ab}(\mathbf{r},t){\hat{S}}_{mn}({\mathbf{r}}^{\prime },t)\rangle$. Clearly, to solve the equations of the one-body density-matrix elements we need to find the values of the two-body correlators first. Akin to these one-body objects, the two-body elements satisfy a new set of the equations. For example, the equation for ρ_44,41 is

$$\begin{equation}\begin{aligned}\hfill \frac{\partial {\rho }_{44,41}}{\partial t}& =-\left[{g}_{41}+{{\Gamma}}_{42}+{{\Gamma}}_{43}+\mathrm{i}V(r)\right]{\rho }_{44,41}\hfill \\ \hfill & \quad -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }{\rho }_{42,41}-\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }{\rho }_{43,41}+\mathrm{i}{{\Omega}}_{\mathrm{c}}({\rho }_{44,21}+{\rho }_{41,24})\hfill \\ \hfill & \quad +\mathrm{i}{{\Omega}}_{\mathrm{d}}({\rho }_{44,31}+{\rho }_{41,34})-\mathrm{i}{{\Omega}}_{\mathrm{p}}{\rho }_{44,42}-\mathrm{i}{{\Omega}}_{\mathrm{s}}{\rho }_{44,43}\hfill \\ \hfill & \quad -\mathrm{i}{N}_{0}\int {\mathrm{d}}^{3}\mathbf{r}V(r)\langle {\hat{S}}_{44}({\mathbf{r}}^{\prime },t){\hat{S}}_{44}(\mathbf{r},t){\hat{S}}_{41}(t)\rangle .\hfill \end{aligned}\end{equation} \tag{ 3a }$$

It depends on other two-body correlators, we present in the following the equations for only two of them, so as to show the key features and save space in the meantime.

$$\begin{equation}\begin{aligned}\hfill \frac{\partial {\rho }_{44,21}}{\partial t}& =-({g}_{21}+{{\Gamma}}_{42}+{{\Gamma}}_{43}){\rho }_{44,21}\hfill \\ \hfill & \quad +\mathrm{i}{{\Omega}}_{\mathrm{c}}{\rho }_{24,21}-\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }({\rho }_{42,21}-{\rho }_{44,41})+\mathrm{i}{{\Omega}}_{\mathrm{d}}{\rho }_{34,21}\hfill \\ \hfill & \quad -\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }{\rho }_{43,21}+\mathrm{i}{{\Omega}}_{\mathrm{p}}({\rho }_{44,11}-{\rho }_{44,22})-\mathrm{i}{{\Omega}}_{\mathrm{s}}{\rho }_{44,23};\hfill \end{aligned}\end{equation} \tag{ 3b }$$

$$\begin{equation}\begin{aligned}\hfill \frac{\partial {\rho }_{41,24}}{\partial t}& =-({g}_{24}-{g}_{41}){\rho }_{41,24}+\mathrm{i}{{\Omega}}_{\mathrm{c}}{\rho }_{24,21}\hfill \\ \hfill & \quad -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }({\rho }_{41,22}-{\rho }_{44,41})+\mathrm{i}{{\Omega}}_{\mathrm{d}}{\rho }_{31,24}-\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }{\rho }_{41,23}\hfill \\ \hfill & \quad +\mathrm{i}{{\Omega}}_{\mathrm{p}}({\rho }_{41,14}-{\rho }_{42,24})-\mathrm{i}{{\Omega}}_{\mathrm{s}}{\rho }_{43,24}\hfill \\ \hfill & \quad -\mathrm{i}{N}_{0}\int {\mathrm{d}}^{3}\mathbf{r}V(r){\hat{S}}_{44}(\mathbf{r},t){\hat{S}}_{41}(t){\hat{S}}_{24}({\mathbf{r}}^{\prime },t)\hfill \\ \hfill & \quad +\mathrm{i}{N}_{0}\int {\mathrm{d}}^{3}\mathbf{r}V(r){\hat{S}}_{41}({\mathbf{r}}^{\prime },t){\hat{S}}_{44}(\mathbf{r},t){\hat{S}}_{24}(t).\hfill \end{aligned}\end{equation} \tag{ 3c }$$

As we can see that the above equations contains three-body correlators, such as $\langle {\hat{S}}_{44}({\mathbf{r}}^{\prime },t){\hat{S}}_{44}(\mathbf{r},t){\hat{S}}_{41}(t)\rangle$. One can easily foresee that the motion of equations for the three-body correlators depend on the four-body ones and so on.

To achieve an effective FWM interaction, the excitation scheme illustrated in figure 1(a) uses a far-off-resonance pumping field Ω_p to avoid the excitation of the atoms to level |2⟩. In virtual of the strong coupling field, a dark state [49] can be formed between |1⟩ and |4⟩ which further reduces the population on |2⟩. Further considering that the Rydberg–Rydberg interation causes energy shift which leads the blockade effect, then we can conclude that the population on |4⟩ is nearly negligible. Under the driving field Ω_d, the FWM process, that is the virtual transition |1⟩ → |2⟩ → |4⟩ → |3⟩ → |1⟩ is initiated and results in the generation of Ω_s without the other transition processes drawing energy from the inputs.

From this point of view, this double-ladder system is very similar to the double-Λ system that has been used in generating the entangled photon pairs [50, 51]. The difference is that the level |4⟩ is one of the ground levels in the double-Λ configuration, instead of a high excited Rydberg state. And we show in the following that this makes a big difference in the underlying physics.

The arrangement of the applied fields allows us to solve the equations of the density-matrix elements using the perturbation method with respect to the far-detuned pumping field and generated (weak) signal field. And the procedure is shown in figure 2 with panel (a) representing our original system. First, we expand the one-body elements with respect to Ω_s as

$$\begin{equation}{\rho }_{mn}={\rho }_{mn}^{(0)}+{\rho }_{mn}^{(1)}+{\rho }_{mn}^{(2)}+\cdots \,.\end{equation} \tag{ 4a }$$

Zoom In Zoom Out Reset image size

And the (0)-order equations obtained from the perturbation is simply the equations (A1)–(A7) listed in appendix A after setting Ω_s = 0. The interaction corresponding to the (0)-order equations is shown in panel (b1) of figure 2, as we can see the signal field is not included. Under the influences of the pumping, coupling and driving fields, a polarization proportional to ${\rho }_{31}^{(0)}$ is built between |1⟩ and |3⟩, and illustrated vividly in panel (b1) that |1⟩ is connected with |3⟩ by the three applied fields. Naturally, the polarization takes a form of ${\mathcal{P}}_{\mathrm{s}}={\chi }^{(3)}{\mathcal{E}}_{\mathrm{p}}{\mathcal{E}}_{\mathrm{c}}{\mathcal{E}}_{\mathrm{d}}^{\ast }$ and leads to the generation of the signal depicted by the effective Hamiltonian [52, 53] that $\propto {\chi }^{(3)}{\mathcal{E}}_{\mathrm{p}}^{\ast }{\mathcal{E}}_{\mathrm{c}}^{\ast }{\mathcal{E}}_{\mathrm{d}}{\hat{\mathcal{E}}}_{\mathrm{s}}^{-}$ with ${\hat{\mathcal{E}}}_{\mathrm{s}}^{-}$ associated with the creation operator at ω_s. We shall not further discuss in detail about the generation but only focus on coefficient χ⁽³⁾. And the first-order system shown in (b2) is not discussed either in this paper, as it relates to the absorption and dispersion of the signal field after being generated, or of a field that is deliberately applied. The solution of ${\rho }_{31}^{(0)}$ provides the exact result of the nonlinear coefficient. However the result is still quite complex. Fortunately, the ladder-type EIT system (|1⟩ − |2⟩ − |4⟩) employs a far-detuned pumping field which naturally becomes our next perturbation parameter with its zeroth-order system shown in panel (c1), the first-order system in panel (c2) and the second in (c3). The corresponding expansion of the elements is

$$\begin{equation}{\rho }_{mn}^{(0)}={\rho }_{mn}^{\left\langle 0\right\rangle }+{\rho }_{mn}^{\left\langle 1\right\rangle }+{\rho }_{mn}^{\left\langle 2\right\rangle }+\cdots \,.\end{equation} \tag{ 4b }$$

Note that the order of perturbation over Ω_p is emphasized by angle brackets $\left\langle \cdot \right\rangle$ in the superscripts. The $\left\langle 0\right\rangle$-order system in panel (c1) is described by a set of equations with Ω_p = 0 and Ω_s = 0. Then one can easily see that this is a trivial system. Without the pumping and the signal, the atoms stay on the ground state and the corresponding elements of density matrix are zeros except ${\rho }_{11}^{\left\langle 0\right\rangle }=1$. Then the $\left\langle 1\right\rangle$-order system in panel (c2) with the corresponding equations listed in appendix A as equation (A8), can help us to reveal the FWM nonlinearity. Simple calculations lead to

$$\begin{equation}\begin{aligned}\hfill {\rho }_{31}^{\left\langle 1\right\rangle }=& -\frac{\mathrm{i}}{\alpha }{{\Omega}}_{\mathrm{c}}{{\Omega}}_{\mathrm{d}}^{\ast }{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill & +{\mathcal{N}}_{a}{g}_{21}{{\Omega}}_{\mathrm{d}}^{\ast }\frac{4\pi }{\alpha }{\int }_{{R}_{\mathrm{b}}}^{\infty }\mathrm{d}rV(r){r}^{2}{\rho }_{44,41}(\mathbf{r},t).\hfill \end{aligned}\end{equation} \tag{ 5 }$$

Here α = g₃₁(|Ω_c|² + g₂₁ g₄₁) + g₂₁|Ω_d|². The polarization built between |3⟩ and |1⟩ is formed via two different mechanisms represented by the two terms on the right-hand side of equation (5). An atom acting with the applied field as a independent (isolated) object leads to the generation of Ω_s via a local nonlinearity. Considering the definitions of the Rabi frequencies and the polarization ${\mathcal{P}}_{\mathrm{s}}=2{N}_{0}{\mu }_{31}{\rho }_{31}^{\left\langle 1\right\rangle }$, we find the corresponding coefficient is

$$\begin{equation}{\chi }_{a}^{(3)}=-{N}_{0}{U}_{0}\frac{\mathrm{i}{\gamma }_{31}^{3}}{{g}_{31}(\vert {{\Omega}}_{\mathrm{c}}{\vert }^{2}+{g}_{21}{g}_{41})+{g}_{21}\vert {{\Omega}}_{\mathrm{d}}{\vert }^{2}}.\end{equation} \tag{ 6 }$$

With

$$\begin{equation}{U}_{0}=\frac{{\mu }_{21}{\mu }_{42}{\mu }_{43}{\mu }_{31}}{4{{\epsilon}}_{0}{\hslash }^{3}{\gamma }_{31}^{3}}.\end{equation} \tag{ 7 }$$

And N₀ U₀ has a dimension of a third-order nonlinear susceptibility. In other words, if the level |4⟩ is not a Rydberg state, then the atomic gases would provide a nonlinearity whose strength is well modeled by ${\chi }_{a}^{(3)}$. However the level |4⟩ being a Rydberg state with a extreme long lifetime makes the local FWM process less efficient because μ₄₂ and μ₄₃ are much smaller than that of a transition involving only the non-Rydberg states.

The second term in equation (5) comes from the RydbergRydberg interaction and we show in the following that this part of interaction dominants in the FWM process. Note that ρ_44,41(r, t) in equation (5) is a full-order quantity with respect to Ω_p. We follow a similar procedure as in equation (4) to expand the two-body correlators as

$$\begin{equation}{\rho }_{ab,mn}={\rho }_{ab,mn}^{(0)}+{\rho }_{ab,mn}^{(1)}+{\rho }_{ab,mn}^{(2)}+\cdots \,,\end{equation} \tag{ 8a }$$

$$\begin{equation}{\rho }_{ab,mn}^{(0)}={\rho }_{ab,mn}^{\left\langle 0\right\rangle }+{\rho }_{ab,mn}^{\left\langle 1\right\rangle }+{\rho }_{ab,mn}^{\left\langle 2\right\rangle }+\cdots \,.\end{equation} \tag{ 8b }$$

The $\left\langle 0\right\rangle$-order elements ${\rho }_{ab,mn}^{\left\langle 0\right\rangle }$ are zeros since the probability of finding an atom on Rydberg state is zero without pumping and signal [see, panel (c1)]. The $\left\langle 1\right\rangle$-order elements ${\rho }_{ab,mn}^{\left\langle 1\right\rangle }$ are zeros as well, due to the fact that the probability of finding a Rydberg atom located at r is still very small under the large-detuned pumping (at most on the second order of Ω_p) [24].

Calculations on perturbation show that the $\left\langle 2\right\rangle$-order non-trivial elements (corresponding to panel (c3)) belongs to two sets of closed equations which are presented in matrix forms as equations (A10) and (A11) in appendix A. The lowest order of two-body correlator ρ_44,41 which appears in equation (5) is $\left\langle 3\right\rangle$, and based on equation (3a), we have

$$\begin{equation}\begin{aligned}\hfill \left[{g}_{41}+{{\Gamma}}_{42}+{{\Gamma}}_{43}+\mathrm{i}V(r)\right]{\rho }_{44,41}^{\left\langle 3\right\rangle }& =\mathrm{i}{{\Omega}}_{\mathrm{c}}\left({\rho }_{44,21}^{\left\langle 3\right\rangle }+{\rho }_{41,24}^{\left\langle 3\right\rangle }\right)+\mathrm{i}{{\Omega}}_{\mathrm{d}}\left({\rho }_{44,31}^{\left\langle 3\right\rangle }+{\rho }_{41,34}^{\left\langle 3\right\rangle }\right)\hfill \\ \hfill & \quad -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }{\rho }_{42,41}^{\left\langle 3\right\rangle }-\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }{\rho }_{43,41}^{\left\langle 3\right\rangle }.\hfill \end{aligned}\end{equation} \tag{ 9a }$$

The above equation depends on other $\left\langle 3\right\rangle$-order two-body elements, such as ${\rho }_{44,21}^{\left\langle 3\right\rangle }$ and ${\rho }_{41,24}^{\left\langle 3\right\rangle }$. Based on equations (3b) and (3c), they satisfy

$$\begin{equation}\begin{aligned}\hfill ({g}_{21}+{{\Gamma}}_{42}+{{\Gamma}}_{43}){\rho }_{44,21}^{\left\langle 3\right\rangle }& =-\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }{\rho }_{43,21}^{\left\langle 3\right\rangle }+\mathrm{i}{{\Omega}}_{\mathrm{p}}{\rho }_{44}^{\left\langle 2\right\rangle }\hfill \\ \hfill & \quad +\mathrm{i}{{\Omega}}_{\mathrm{c}}{\rho }_{24,21}^{\left\langle 3\right\rangle }-\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\left({\rho }_{42,21}^{\left\langle 3\right\rangle }-{\rho }_{44,41}^{\left\langle 3\right\rangle }\right)+\mathrm{i}{{\Omega}}_{\mathrm{d}}{\rho }_{34,21}^{\left\langle 3\right\rangle };\hfill \end{aligned}\end{equation} \tag{ 9b }$$

$$\begin{equation}\begin{aligned}\hfill ({g}_{24}+{g}_{41}){\rho }_{41,24}^{\left\langle 3\right\rangle }& =\mathrm{i}{{\Omega}}_{\mathrm{d}}{\rho }_{31,24}^{\left\langle 3\right\rangle }-\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }{\rho }_{41,23}^{\left\langle 3\right\rangle }\hfill \\ \hfill & \quad +\mathrm{i}{{\Omega}}_{\mathrm{p}}{\rho }_{41,14}^{\left\langle 2\right\rangle }+\mathrm{i}{{\Omega}}_{\mathrm{c}}{\rho }_{24,21}^{\left\langle 3\right\rangle }-\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\left({\rho }_{41,22}^{\left\langle 3\right\rangle }-{\rho }_{44,41}^{\left\langle 3\right\rangle }\right).\hfill \end{aligned}\end{equation} \tag{ 9c }$$

The whole set of the algebraic equations for the $\left\langle 3\right\rangle$-order two-body elements has 27 unknowns including ${\rho }_{44,41}^{\left\langle 3\right\rangle }$, ${\rho }_{44,21}^{\left\langle 3\right\rangle }$ and ${\rho }_{41,24}^{\left\langle 3\right\rangle }$ which are listed in appendix A as an array, see equation (A12). As for their algebraic equations, we find that they are too tedious and not suitable to be presented in this paper even for appendix A, thus we only provide equation (9) here as examples. Note that those $\left\langle 3\right\rangle$-order equations relate to the $\left\langle 2\right\rangle$-order one/two-body elements. And the $\left\langle 2\right\rangle$-order elements depend on $\left\langle 1\right\rangle$-order elements whose equations can be easily solved. Unfortunately, the solution to ${\rho }_{44,41}^{\left\langle 3\right\rangle }$ does not have a compact form, but we managed to write the expression into a power series with respect to Ω_c and Ω_d:

$$\begin{equation}{\rho }_{44,41}^{\left\langle 3\right\rangle }(r)=\sigma (r){{\Omega}}_{\mathrm{c}}{{\Omega}}_{\mathrm{p}}.\end{equation} \tag{ 10 }$$

With

$$\begin{equation}\begin{aligned}\hfill \sigma (r)& =-\vert {{\Omega}}_{\mathrm{p}}{\vert }^{2}\left\{\vert {{\Omega}}_{\mathrm{c}}{\vert }^{2}\left[\frac{1}{2\beta (r){\gamma }_{41}({g}_{24}+{g}_{41})}\left(\frac{{g}_{31}}{\alpha ({g}_{41}+{g}_{12})}+\mathrm{c}.\mathrm{c}.\right)\right.\right.\hfill \\ \hfill & \left.\left.\quad +\frac{1}{\beta (r)({{\Gamma}}_{42}+{{\Gamma}}_{43})({g}_{21}+{{\Gamma}}_{42}+{{\Gamma}}_{43})}\left(\frac{{g}_{31}}{\alpha {g}_{42}}+\mathrm{c}.\mathrm{c}.\right)\right.\right.\hfill \\ \hfill & \left.\left.\quad +\frac{2{g}_{31}}{\alpha \beta (r)({g}_{21}+{g}_{41})\left[2{g}_{41}+\mathrm{i}V(r)\right]\left[{g}_{41}+{g}_{42}+\mathrm{i}V(r)\right]}\right]+\cdots +Q(m,n)\vert {{\Omega}}_{\mathrm{c}}{\vert }^{2m}\vert {{\Omega}}_{\mathrm{d}}{\vert }^{2n}+\cdots \,\right\}.\hfill \end{aligned}\end{equation} \tag{ 11 }$$

With m = 2, 3, ..., and n = 1, 2, .... The parameter α is defined in the previous discussion and β(r) = g₄₁ + Γ₄₂ + Γ₄₃ + iV(r). Q(m, n) represents a general coefficient of the series. The exact solution of ${\rho }_{44,41}^{\left\langle 3\right\rangle }$ should include all terms of the above expansion, and we have to calculate it numerically. From the equation (11) we can see that the resultant nonlinear coefficient depends on the squared modular of the coupling and driving Rabi frequencies, suggesting that nonlocal FWM process is dressed by the strong fields.

Substituting equations (10) and (11) into equation (5), we find nonlocal nonlinear coefficient is

$$\begin{equation}{\chi }_{R}^{(3)}={N}_{0}^{2}{U}_{0}{\gamma }_{31}^{3}{g}_{21}\frac{4\pi }{\alpha }{\int }_{{R}_{\mathrm{b}}}^{\infty }\mathrm{d}rV(r){r}^{2}\sigma (r).\end{equation} \tag{ 12 }$$

In evaluating (12), we first determine the value of R_b, which under large pumping detuning, is ${R}_{\mathrm{b}}=\sqrt[6]{{C}_{6}/{\delta }_{\text{EIT}}}$ with ${\delta }_{\text{EIT}}={{\Omega}}_{\mathrm{c}}^{2}/{{\Delta}}_{\mathrm{p}}$ [13, 54], then the result is obtained from the integral with simple calculations, and is shown in figure (3).

Zoom In Zoom Out Reset image size

Based on the results of Re ${\chi }_{R}^{(3)}$ in figure 3(a) and Im ${\chi }_{R}^{(3)}$ in (b), we see that the nonlocal nonlinearity reaches its maximal strength at the double-photon resonance Δ_c + Δ_p = 0. This is physically intuitive as, under the same condition, the maximal coherence between |1⟩ and |4⟩ is achived. The center frequency of the generated field must satisfy Δ_p + Δ_c = Δ_d + Δ_s due to the conservation law of the energy. And for the parameters we used in figure 3, Δ_s = 50γ₃₁.

As comparison we plot the local nonlinear susceptibilities in figure 3 as well to show that amplitude of ${\chi }_{R}^{(3)}$ is much larger than that of ${\chi }_{a}^{(3)}$ at the given atomic density. The local nonlinearity comes solely from the interation between the photon and the atoms which can be viewed, as we stated before, a coherence between |1⟩ and |3⟩ built by the three applied fields leading to a polarization that induces the emission of the signal.

On the other hand, the nonlocal nonlinearity has a different mechanism. In virtue of the coupling and pumping field, a few atoms are excited to the Rydberg states and repel each other via the extra large dipole moments. This causes the addition potential energy V(r) which at the same time suggesting the existence of a tendency that the atomic gases would reduce its energy to a low state. Such tendency corresponds to the two-body element in equation (5), that is ${\rho }_{44,41}(\mathbf{r},t)=\langle {\hat{S}}_{44}(\mathbf{r},t){\hat{S}}_{41}(\mathbf{0},t)\rangle$, a correlation between probability of the atom at r staying on Rydberg state |4⟩ and that of anther Rydberg atom (at origin) jumping to the ground state. Since the two operators commutate with each other, the same argument can be made when interchanging the two atoms. Of course, such transition from |4⟩ to |1⟩ is unpractical due to the same parity of two levels. However with the help of driven field (${{\Omega}}_{\mathrm{d}}^{\ast }$ in equation (5)), it manifest as a dominant part of the effective polarization and strengthen the FWM process.

Clearly, higher atomic density leads to shorter distance between atoms, thus significantly increases the two-body correlation and results in a much larger nonlinearity. In figure 4 we show the real part (a) and the imaginary part (b) of ${\chi }_{R}^{(3)}$ as a function of atomic density. Comparing with the local part of nonlinear coefficient, ${\chi }_{R}^{(3)}$ is increased more significantly with the atomic density. Note that the method we used here requires a limitation on the atomic population so that we only need to consider the two-body elements without taking into account of the correlation between three or even more atoms. We can roughly estimate the proper atomic density via relation $1< {N}_{a}{\rho }_{33}^{(2)}\cdot \frac{4}{3}\pi {(3{R}_{\mathrm{b}})}^{3}< 2$ which means that given a blockage sphere and its close-packed one, in the space they occupy which is approximately a sphere having a radius of 3R_b, two atoms can be excited to Rydberg states with considerable opportunity, but never for three. For our parameters that leads to a condition that the proper atomic density should be smaller than 1.0 × 10¹¹ cm⁻³.

Zoom In Zoom Out Reset image size

The possible experimental setup would be very similar to the system used in reference [44]. And since ${\chi }_{a}^{(3)}$ is proportional to N₀, while ${\chi }_{R}^{(3)}$ is proportional to ${N}_{0}^{2}$, see equations (6) and (12). One can verify the mechanism of the generation of Ω_s by examining the relation between the count rate of Ω_s (proportional to squared modulus of the nonlinear susceptibility [50, 55]) and the atomic density. The nonlinearity based on Rydberg–Rydberg interaction would manifest itself as a dependence of the count rate on ${N}_{0}^{4}$.

The FWM nonlinearity plays important roles in lots of applications in quantum information. And in this paper, we have investigated the enhancement of the third-order FWM nonlinearity via the correlation between two Rydberg atoms. Using the one-body and two-body elements of the density matrix, we solved the associated equations using perturbation method and show that the polarization which is induced by the pumping, coupling and driving field, and responsible for generating the signal field is composed of a local part from the interaction between atoms and photons and a nonlocal part due to the Rydberg–Rydberg interaction. Using a Rydberg state with n = 60 as an example, we show via the numerical results that the nonlocal nonlinear susceptibility dominates in the total nonlinear coefficient. For the Rydberg state with a even larger quantum number, beside the above mentioned result, one finds that the Rydberg–Rydberg interaction becomes stronger and that leads to a larger blockade sphere. If we only focus on the two-body correlation, then a relatively low atomic density should be adopt and calculations show that the nonlinearity is actually weaker owning to the relation ${\chi }_{R}^{(3)}\propto {N}_{0}^{2}$. With a higher atomic density, the three-body correlation should be considered and the nonlinearity could be further enhanced.

The work is supported by the National Natural Science Foundation of China (No. 12074061), the Fundamental Research Funds for the Central Universities (No. 2412020FZ028) and the Scientific and Technological Research Program of Jilin Education Department (No. JJKH20211280KJ).

No new data were created or analysed in this study.

For simplicity, we use complex decay rates g_mn in the following equations, and they are defined as the composed of decoherence rates and detunings. g₂₁ = γ₂₁ + iΔ_p, g₃₁ = γ₃₁ + iΔ_s, g₄₁ = γ₄₁ + i(Δ_c + iΔ_p), g₃₂ = γ₃₂ + i(Δ_s − iΔ_p), g₄₂ = γ₄₂ + iΔ_c, g₄₃ = γ₄₃ + iΔ_d. The equations for one-body equations for

$$\begin{equation}{\partial }_{t}{\rho }_{11}={{\Gamma}}_{21}{\rho }_{22}+{{\Gamma}}_{31}{\rho }_{33}-(\mathrm{i}{\rho }_{12}{{\Omega}}_{\mathrm{p}}+\mathrm{i}{\rho }_{13}{{\Omega}}_{\mathrm{s}}+\mathrm{c}.\mathrm{c});\end{equation} \tag{ A1 }$$

$$\begin{equation}{\partial }_{t}{\rho }_{22}={{\Gamma}}_{42}{\rho }_{44}-{{\Gamma}}_{21}{\rho }_{22}+(\mathrm{i}{\rho }_{12}{{\Omega}}_{\mathrm{p}}-\mathrm{i}{\rho }_{24}{{\Omega}}_{\mathrm{c}}+\mathrm{c}.\mathrm{c}.);\end{equation} \tag{ A2 }$$

$$\begin{equation}{\partial }_{t}{\rho }_{33}={{\Gamma}}_{43}{\rho }_{44}-{{\Gamma}}_{31}{\rho }_{33}+(\mathrm{i}{\rho }_{13}{{\Omega}}_{\mathrm{s}}-\mathrm{i}{\rho }_{34}{{\Omega}}_{\mathrm{d}}+\mathrm{c}.\mathrm{c});\end{equation} \tag{ A3 }$$

$$\begin{equation}{\partial }_{t}{\rho }_{21}=-{g}_{21}{\rho }_{21}+\mathrm{i}{{\Omega}}_{\mathrm{p}}({\rho }_{11}-{\rho }_{22})-\mathrm{i}{\rho }_{23}{{\Omega}}_{\mathrm{s}}+\mathrm{i}{\rho }_{41}{{\Omega}}_{\mathrm{c}}^{\ast };\end{equation} \tag{ A4 }$$

$$\begin{equation}{\partial }_{t}{\rho }_{32}=-{g}_{32}{\rho }_{32}+\mathrm{i}{\rho }_{12}{{\Omega}}_{\mathrm{s}}-\mathrm{i}{\rho }_{31}{{\Omega}}_{\mathrm{p}}^{\ast }-\mathrm{i}{\rho }_{34}{{\Omega}}_{\mathrm{c}}+\mathrm{i}{\rho }_{42}{{\Omega}}_{\mathrm{d}}^{\ast };\end{equation} \tag{ A5 }$$

$$\begin{equation}{\partial }_{t}{\rho }_{42}=-{g}_{42}{\rho }_{42}+\mathrm{i}{{\Omega}}_{\mathrm{c}}({\rho }_{22}-{\rho }_{44})+\mathrm{i}{\rho }_{32}{{\Omega}}_{\mathrm{d}}-\mathrm{i}{\rho }_{41}{{\Omega}}_{\mathrm{p}}^{\ast }-\mathrm{i}{N}_{a}\int {\mathrm{d}}^{3}\mathbf{r}V(r){\rho }_{44,42}(\mathbf{r},t);\end{equation} \tag{ A6 }$$

$$\begin{equation}{\partial }_{t}{\rho }_{43}=-{g}_{43}{\rho }_{43}+\mathrm{i}{\rho }_{23}{{\Omega}}_{\mathrm{c}}+\mathrm{i}{{\Omega}}_{\mathrm{d}}({\rho }_{33}-{\rho }_{44})-\mathrm{i}{\rho }_{41}{{\Omega}}_{\mathrm{s}}^{\ast }-\mathrm{i}{N}_{a}\int {\mathrm{d}}^{3}\mathbf{r}V(r){\rho }_{44,43}(\mathbf{r},t).\end{equation} \tag{ A7 }$$

Equations for ρ₃₁, ρ₃₂ and ρ₄₁ are presented in the main text as equations (2a) and (2b), respectively. In the above equations, the phase-matching condition Δ_s = Δ_c + Δ_p − Δ_d is already assumed. Under the double perturbation we find that the $\left\langle 1\right\rangle$-order equations for one-body elements are

$$\begin{equation}\left(\begin{matrix}\hfill -{g}_{21}\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill \\ \hfill 0\hfill & \hfill -{g}_{31}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill \\ \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill -{g}_{41}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\rho }_{21}^{\left\langle 1\right\rangle }\hfill \\ \hfill {\rho }_{31}^{\left\langle 1\right\rangle }\hfill \\ \hfill {\rho }_{41}^{\left\langle 1\right\rangle }\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill -\mathrm{i}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill 0\hfill \\ \hfill \mathrm{i}{{\Delta}}_{44,41}\hfill \end{matrix}\right).\end{equation} \tag{ A8 }$$

Here ${{\Delta}}_{44,ij}={\mathcal{N}}_{a}\int {\mathrm{d}}^{3}{r}^{\prime }V({\mathbf{r}}^{\prime }-\mathbf{r}){\rho }_{44,ij}({\mathbf{r}}^{\prime },\mathbf{r},t)$. In order to solve the $\left\langle 3\right\rangle$-order equations for the two-body elements we also need the value of $\left\langle 2\right\rangle$-order one-body element, and they reads

$$\begin{equation}\left(\begin{matrix}\hfill -{{\Gamma}}_{21}\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill & \hfill {{\Gamma}}_{42}\hfill \\ \hfill 0\hfill & \hfill -{g}_{23}\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill \\ \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill -{g}_{24}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{g}_{32}\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{{\Gamma}}_{31}\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill {{\Gamma}}_{43}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill -{g}_{34}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill \\ \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{g}_{42}\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill \\ \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -{g}_{43}\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill -{{\Gamma}}_{42}-{{\Gamma}}_{43}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\rho }_{22}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{23}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{24}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{32}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{33}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{34}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{42}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{43}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{44}^{\left\langle 2\right\rangle }\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill 2\,\mathrm{I}\mathrm{m}[{\rho }_{12}^{(1)}{{\Omega}}_{\mathrm{p}}]\hfill \\ \hfill -\mathrm{i}{\rho }_{13}^{(1)}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill -\mathrm{i}{\rho }_{14}^{(1)}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill \mathrm{i}{\rho }_{31}^{(1)}{{\Omega}}_{\mathrm{p}}^{\ast }\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill \mathrm{i}{\rho }_{41}^{(1)}{{\Omega}}_{\mathrm{p}}^{\ast }\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ A9 }$$

The nonzero two-body elements on $\left\langle 2\right\rangle$ order are the solutions to two set of closed equations. In matrix form, they are

$$\begin{equation}\left(\begin{matrix}\hfill {M}_{1}\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {M}_{2}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill \\ \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill {M}_{3}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {M}_{4}\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {M}_{5}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill {M}_{6}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill \\ \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {M}_{7}\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill \\ \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {M}_{8}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill {M}_{9}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\rho }_{21,12}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{31,12}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{41,12}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{21,13}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{31,13}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{41,13}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{21,14}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{31,14}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{41,14}^{\left\langle 2\right\rangle }\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill 2\,\mathrm{I}\mathrm{m}[{\rho }_{12}^{(1)}{{\Omega}}_{\mathrm{p}}]\hfill \\ \hfill \mathrm{i}{\rho }_{31}^{(1)}{{\Omega}}_{\mathrm{p}}^{\ast }\hfill \\ \hfill \mathrm{i}{\rho }_{41}^{(1)}{{\Omega}}_{\mathrm{p}}^{\ast }\hfill \\ \hfill -\mathrm{i}{\rho }_{13}^{(1)}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill -\mathrm{i}{\rho }_{14}^{(1)}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ A10 }$$

$$\begin{equation}\left(\begin{matrix}\hfill {M}_{10}\hfill & \hfill 0\hfill & \hfill 2\mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {M}_{11}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill & \hfill 0\hfill \\ \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill {M}_{12}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}^{\ast }\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {M}_{13}\hfill & \hfill 2\mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill {M}_{14}\hfill & \hfill \mathrm{i}{{\Omega}}_{\mathrm{d}}^{\ast }\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 2\mathrm{i}{{\Omega}}_{\mathrm{c}}\hfill & \hfill 0\hfill & \hfill 2\mathrm{i}{{\Omega}}_{\mathrm{d}}\hfill & \hfill {M}_{15}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\rho }_{21,21}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{31,21}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{41,21}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{31,31}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{41,31}^{\left\langle 2\right\rangle }\hfill \\ \hfill {\rho }_{41,41}^{\left\langle 2\right\rangle }\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill -2\mathrm{i}{\rho }_{21}^{(1)}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill -\mathrm{i}{\rho }_{31}^{(1)}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill -\mathrm{i}{\rho }_{41}^{(1)}{{\Omega}}_{\mathrm{p}}\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ A11 }$$

With M₁ = −g₁₂ − g₂₁, M₂ = −g₁₂ − g₃₁, M₃ = −g₁₂ − g₄₁, M₄ = −g₁₃ − g₂₁, M₅ = −g₁₃ − g₃₁, M₆ = −g₁₃ − g₄₁, M₇ = −g₁₄ − g₂₁, M₈ = −g₁₄ − g₃₁, M₉ = −g₁₄ − g₄₁. M₁₀ = −2g₂₁, M₁₁ = −g₂₁ − g₃₁, M₁₂ = −g₂₁ − g₄₁, M₁₃ = −2g₃₁, M₁₄ = −g₃₁ − g₄₁, M₁₅ = −2g₄₁ − iV.

The $\left\langle 3\right\rangle$-order equations depends on the above $\left\langle 2\right\rangle$-order elements. To obtain ${\rho }_{44,41}^{\left\langle 3\right\rangle }$ we need to solution a set of equations with 27 unknowns, and they are:

$$\begin{equation}\begin{aligned}\hfill & \left({\rho }_{22,21}^{\left\langle 3\right\rangle },{\rho }_{23,21}^{\left\langle 3\right\rangle },{\rho }_{24,21}^{\left\langle 3\right\rangle },{\rho }_{32,21}^{\left\langle 3\right\rangle },{\rho }_{33,21}^{\left\langle 3\right\rangle },{\rho }_{34,21}^{\left\langle 3\right\rangle },{\rho }_{42,21}^{\left\langle 3\right\rangle },{\rho }_{43,21}^{\left\langle 3\right\rangle },{\rho }_{44,21}^{\left\langle 3\right\rangle },{\rho }_{31,22}^{\left\langle 3\right\rangle },{\rho }_{41,22}^{\left\langle 3\right\rangle },{\rho }_{31,23}^{\left\langle 3\right\rangle },{\rho }_{41,23}^{\left\langle 3\right\rangle },\right.\hfill \\ \hfill & {\left.{\rho }_{31,24}^{\left\langle 3\right\rangle },{\rho }_{41,24}^{\left\langle 3\right\rangle },{\rho }_{32,31}^{\left\langle 3\right\rangle },{\rho }_{33,31}^{\left\langle 3\right\rangle },{\rho }_{34,31}^{\left\langle 3\right\rangle },{\rho }_{42,31}^{\left\langle 3\right\rangle },{\rho }_{43,31}^{\left\langle 3\right\rangle },{\rho }_{44,31}^{\left\langle 3\right\rangle },{\rho }_{41,32}^{\left\langle 3\right\rangle },{\rho }_{41,33}^{\left\langle 3\right\rangle },{\rho }_{41,34}^{\left\langle 3\right\rangle },{\rho }_{42,41}^{\left\langle 3\right\rangle },{\rho }_{43,41}^{\left\langle 3\right\rangle },{\rho }_{44,41}^{\left\langle 3\right\rangle }\right)}^{\mathrm{T}}.\hfill \end{aligned}\end{equation} \tag{ A12 }$$

They are obtained following the same procedure of finding the matrices of $\left\langle 2\right\rangle$-order equations. The matrix equation are not shown here since the size of the matrix is too large.