Geometric representation and the adiabatic geometric phase in four-wave mixing processes

The application of the adiabatic geometric phase (AGP) to nonlinear frequency conversion may help to develop new types of all-optical devices, which leads to all-optical modulation of the phase front of one wave by the intensity of other waves. In this paper, we develop the canonical Hamilton equation and a corresponding geometric representation for two schemes of four-wave mixing (FWM) processes (ω1 +ω2 =ω3 +ω4 and ω1 +ω2 +ω3 =ω4), which can precisely describe and calculate the AGP controlled by the quasi-phase matching technique. The AGPs of the idler (ω1) and signal (ω4) waves for these two schemes of FWM are studied systematically when the two pump waves (ω2 and ω3) are in either the undepleted or in the depleted pump cases, respectively. The analysis reveals that the proposed methods for calculating the AGP are universal in both cases. We expect that the analysis of AGP in FWM processes can be applied to all-optically shaping or encoding of ultrafast light pulse. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
When an eigenstate slowly evolves around a closed trajectory in the parameter space by varying its Hamiltonian, an adiabatic geometric phase (AGP) factor is acquired in addition to the dynamical phase, which is given by the time integral of the eigenenergy [1]. Moreover, the AGP can be interpreted as the flux of a virtual field through the surface enclosed by the closed circuit in the parameter space [2]. Over the past couple of decades, important applications of the AGP have been demonstrated in quantum physics [3], for example, in high-precision quantum measurement [4], quantum information processing [5], quantum computing [6], Bose-Einstein condensates [7,8] and condensed-matter physics [9]. AGPs also have extensive applications in optics, in which they are used for beam shaping [10,11], beam shearing [12], guiding [13], holography [14] and encoding both classical [15] and quantum [16] information.
Recently, the geometric phase for χ (2) nonlinear frequency conversion was studied [17][18][19]. Some manifestations of the geometric phase for light were obtained by changing the polarization or the photon direction in nonlinear metasurfaces [20][21][22][23][24]. A different case is when the AGP is accumulated in a bulk material by nonlinear modulation pattern along the nonlinear crystal through the quasi-phase matching (QPM) technique. The analysis in literatures [18,19] provides an effective geometric representation of the latter case and a method to describing and calculating the AGP in a precise way. This method is universal for the undepleted pump (the intensity of the pump wave is much stronger than the intensities of the weak signal and idler) and depleted pump (the intensity of the pump wave is comparable to the intensities of the signal and idler) cases. Moreover, nonreciprocal transmission and wavefront shaping through AGPs were realized experimentally by circular rotation of the QPM parameters using sum-frequency generation (SFG) in the χ (2) process [25]. These methods provide new ways to manipulate AGP by nonlinear optics.
Adiabatic conversion in nonlinear process is not limited only to three wave mixing processes. Recently, adiabatic four wave mixing was theoretically studied, through either a longitudinally varying silicon waveguide structure or a linearly tapered fluoride step-index fiber [26]. It was assumed that two of the four waves serve as undepleted pump waves, so that the adiabatic exchange of energy occurs between the two remaining waves. Very recently, such an undepleted adiabatic passage in FWM was experimentally observed in a photonic crystal fiber with a tapered core diameter [27]. Here we would like to address the question of how to accumulate geometric phase in a four wave mixing process and how it can be calculated. Unlike Refs. [26,27], our analysis is not limited to the case of undepleted pump waves. In order to analyze the systems dynamics, we develop a canonical Hamiltonian equation and a corresponding geometrical representation of all the system states on a closed surface in three dimensions. This surface is the Bloch sphere in the undepleted case, but for the depleted case the shape becomes elongated and one of the poles attains a sharp conical shape if the system reaches the full depletion condition.
The remainder of this paper is structured as follows. The theory and model for the QPM in the FWM process are illustrated in Sec. II; two schemes of the FWM process (ω 1 + ω 2 = ω 3 + ω 4 and ω 1 + ω 2 + ω 3 = ω 4 ) are considered in this section. The AGPs in schemes I and II for the undepleted, and then depleted pump cases are discussed in Sec. III and IV, respectively. The main results of this paper are summarized in Sec. V.

Model and theory
Generally, the adiabatic process for nonlinear wave conversion requires to control the wave-vector mismatch parameter. Reference [26] introduces a small periodic change in the waveguide width with a longitudinal chirp to the period to achieve this parameter, as illustrated in Fig. 1(a), where the local period Λ(z) adds a quasi momentum K g (z) = 2π/Λ(z) to the wave-vector difference between the guided modes of the interacting waves. This type of corrugated waveguide can be readily fabricated using planar integrated optics technologies. By introducing a duty cycle and its relative position within the periodic cycle (i.e., the phase), three QPM parameters, namely, the modulation wave vector, modulation phase, and duty cycle, are introduced in this setting. In this paper, two schemes of the FWM process are considered; their energy diagrams are displayed in Figs. 1(b,c). Following the definition in Ref. [26], we also employ waves ω 1 and ω 4 to denote the idler and signal, respectively, and we use waves ω 2 and ω 3 to denote the two pumps, respectively. Moreover, we discuss not only the AGP in the undepleted case, i.e., I 2,3 ≫ I 1,4 (where I j is the intensity of wave ω j ), but also the AGP in the depleted case, i.e., I 2,3 ∼ I 1,4 .

Dimensionless equation for the FWM process of scheme I
The energy diagram for the FWM process of scheme I is shown in Fig. 1(b), which satisfies ω 1 + ω 2 = ω 3 + ω 4 . The coupled amplitude equations for this process with QPM [28] are where n j is the linear refractive index for ω j , c is the speed of light in vacuum, where χ (3) is the third-order susceptibility and A eff is the core area of the waveguide, which varies along z.
If an adiabatic modulation is added to d(z), as mentioned above, periodically changing the width of the waveguide [26] [see Fig. 1(a)], d(z) can be expanded by a Fourier series with slowly varying components: where K g (z) is the modulation wave-vector, ϕ d (z) is the modulation phase, and d m = d m (D(z)) with 0 ≤ D(z) ≤ 1 is the m-th Fourier coefficient, that depends on the duty cycle of the modulation.
Here, the influence of D(z) on d 0 is neglected, and d 0 is assumed to be constant. Equation (3) is substituted into Eq. (1), and the right-hand side is kept with the terms d 0 and d 1 because their oscillations are closest to the phase matching. Furthermore, by applying a transformation from A j to a rotating frame byÃ j , and substituting the above expressions into Eq. (1), the coupled four-wave equations change to the following: It is assumed that the first order QPM term, m = 1, of Eq. (3) is used for coupling the first two equations of Eq. (5), whereas the temp m = −1 is used to couple the last two equations. Therefore, we can distinguish q 1 and q 2 by whether the +1 or −1 order is used for the coupling. For the FWM process of scheme I, it is convenient to assume that we can neglect the small differences between the coupling constants owing to the different angular frequencies and refractive indices, so that σ 0 ≈ d 0 ω/cn and σ 1 ≈ d 1 ω/cn (where ω 1,2,3,4 ≈ ω and n 1,2,3,4 ≈ n).
The following definitions are applied: In these dimensionless parameters, η represents the coupling between the waves, τ -the interaction length, Γ -the phase mismatch and g -the duty cycle and the modulation phase. Thus, Eq. (5) can be changed to a set of dimensionless equations as where q j with j = 1, 2, 3, 4 representing the photon fluxes of the four waves. Obviously, the four waves are conserved by the total norm at any time as Actually, by applying a transformation, some of the cross-phase modulation (XPM) terms in Eq. (7) can be removed, and the equations can be further simplified as Generalizing the concepts derived by Luther et. al. for a three wave mixing process [29], and their implementation for adiabatic processes [30,31] to the case of four wave mixing Eqs. (10) can be presented via a canonical Hamiltonian structure, namely, Moreover, the Hamiltonian gives rise to 6 constant motions, which arẽ These constants are the manifestations of the Manley-Rowe (MR) relations for the scheme I of FWM.

Dimensionless equation for the FWM process of scheme II
The energy diagram for the FWM process of scheme II is illustrated in Fig. 1(c), which satisfies ω 4 = ω 1 + ω 2 + ω 3 . Therefore, the coupled amplitude equations for this scheme are Following the same disposition as Eq. (3) and a new transformation from A j to a rotating frame byÃ j , Equations (13) change to the following: The first order positive term, m = 1 is used for the coupling in the first two equations and the first order negative term, m = −1 is used for the coupling in the last two equations. Because ω 4 >ω 1,2,3 , we must apply the following new definitions to dimensionless Eq. (15). The new definitions areÃ and where |Ã l0 | 2 is the initial intensity of the l-th wave. Hence, Eqs. (15) change to where q j with j = 1, 2, 3, 4 represent the photon fluxes of the four waves. Obviously, Eq. (16) makes the initial condition of the waves satisfy the normalized condition where I j = |q (0) j | 2 is the initial intensity of q j . Because the creation (or the annihilation) one photon on ω 4 needs annihilation (or creation) 3 photons on ω 1,2,3 , respectively, in scheme II, the total norms are not conserved in the entire dynamical process. Hence, Eq. (20) is valid only at τ = 0, which is different from case I, where the norm is conserved for every τ [see Eq. (8)]. The non-conservation of the total norm in scheme II leads to the XPM terms in Eq. (19) cannot be canceled, which results in a Stark shift, i.e. to intensity dependent phase matching.
Equations (19) can also be presented via a canonical Hamiltonian structure, Similarly, the Hamiltonian in Eq. (21) gives rise to 6 new constant motions, i.e., the MR relation for scheme II of FWM, which are To simplify the calculation, we assume that the following approximations can be made for calculating the nonlinear coupling terms s ij Then, according to Eq. (17), the coefficients of the self-phase modulation (SPM) and XPM terms,

Geometrical representation and calculation of the AGP
To depict the geometric motion of Eqs. (10) and (19) in the state space, we introduce a set of coordinates (X, Y, Z) with the surface φ = 0 for these two schemes, which adopt the similar rule of three-wave mixing in Refs. [19,29,31] (using the real and the imaginary parts of product of the waves in the interacting part of the Hamiltonian as the coordinates X and Y , respectively, and using the intensity of one wave as the coordinate Z): Scheme II : In addition, the state vector is given by W = Xî + Yĵ + Zk. Because the last terms in Eqs. (26) and (28) automatically cancel the sum of first two terms, i.e., X 2 + Y 2 , all the states of the system are mapped onto the closed surface that satisfy φ = 0. The surfaces in these two schemes can be termed the Bloch surfaces for these two schemes of the FWM process.
For the QPM control parameters, we can still use the design of the circular rotation of [18,19]: Here, the symbols "+" and "−" correspond to the clockwise and counterclockwise rotation, respectively, of the QPM vector, expressed as Q = (Ξ cos ϕ d , Ξ sin ϕ d , ∆Γ). Θ is the angle of the normal vectorΘ of the plane containing the circular trajectory of vector Q with respect to the vertical axis [see Fig. 1(d)]. Hence, this circular rotation is fully characterized by the angle Θ. Ω = 2π/T is the angular velocity, which denotes the speed of the rotation (here, T = ηL is the scaled length of the waveguide). For convenience, we fix Γ 0 = 1 in the numerical simulation. Moreover, to ensure the rotation satisfies the adiabatic condition, we select T = 1000, which is a sufficient scale for Ω to perform an adiabatic rotation. Following the framework of three-wave mixing [19], the AGP in the QPM process for wave ω j can be generally expressed by where Φ ± j is the total phase shift of the wave ω j for the clockwise and counter-clockwise rotation, respectively, D j is the dynamical phase, which equals to the circular integral of the eigenenergy along the rotation, γ j is the acquired geometric phase. Hence, the AGP can be calculated by two equivalent methods. The first one is based on the half of the difference of the total phase between the clockwise and the counter-clockwise rotation of the vector Q. The expression of this method is written as where Φ ± j are the total phase shift of wave ω j for the clockwise and counter-clockwise circular rotation, respectively. The second one equals to the difference between the total phase accumulated along the circular trajectory and that along its vertical projection. The expression of this method is written as where Φ c j is the total phase shift of wave ω j for the circular rotation (Here "⃝ =⟳ or ⟲" refers to the clockwise or counter-clockwise circular integral, respectively), Φ rt j is the total phase shift of the vertical projection and "↿⇂" denotes a round trip motion integral, which is governed only by Eqs. (29) and (30). The phase factor of QPM is unchanged in this motion. According to the analysis in Ref. [16], the round-trip motion is a degenerated motion with the clockwise and the counterclockwise circular rotation, they share the same eigenenergy. However, unlike the clockwise and the counterclockwise rotation, the round-trip motion does not accumulate the AGP during its motion because the variation of ϕ d is absent. Hence, in Eq. (32), Φ rt j = ∫ ↿⇂ p ′ j dq ′ j = D j , we can calculate the AGP analytically by the difference of total phase shift between these two types of motion. For the FWM of scheme I, We shall now examine the AGP in the the FWM schemes, for different values of the pump intensities with respect to the signal and idler intensities, starting from the undepleted pumps and up to the case of fully depleted pumps. The different results are summarized in Tables 1 and 2.
at which the intensities of the idler and signal are much smaller than those of the two pumps. If we assume I 4 = 0, the dynamics starts from the bottom of the Bloch sphere. We first discuss the simplest case, I 2 = I 3 , at which the two pumps have a the same intensity. According to the conservative condition in Eq. (8), the limit of this condition isq * 3q 2 =q * 2q 3 = |q 2 | 2 = |q 3 | 2 ≈ 1/2, and the magnitudes of |q 1 | 2 and |q 4 | 2 can be neglected. Hence, Eqs. (10) are simplified to the spin-1/2 Bloch equation as For this limit, the Hamiltonian matrix in Eq. (36) is strictly follows the circular rotation in the parameter space, as described by Eqs. (29)-(31), driven by Eqs. (29)- (31), which can give rise to the AGP as [18] which is true for idler waves (our initial condition is q 4 = 0). To test numerically this analytical limit, we adopt the following initial condition: Obviously, this initial condition satisfies the undepleted condition of I 2,3 ≫ I 1,4 . Then, we substitute this initial condition into the full FWM equations described by Eqs. (10) to calculate the AGP through Eq. (33) or Eq. (34). A comparison between γ l and the numerical γ 1 for different values of Θ is shown in Fig. 2(a), which demonstrates that the numerical results are in accordance with the theoretical prediction. Hence, the two methods for calculating the AGP in three-wave mixing remain valid for the case of FWM. Typical examples for 3 types of trajectories on the Bloch sphere, namely, clockwise, counterclockwise, and round-trip motion, are displayed in Fig. 2(b). The AGP accumulation calculated via Eq. (33) is compared with that calculated via Eq. (34) in Fig. 2(c), which demonstrates that these two methods remain equivalent to each other for the case of FWM. Hence, we will use Eq. (34) to calculate the AGP accumulations in the following discussion.
However, if we assume I 2 ≠ I 3 , a factor, which is the modulus ofq * 2q 3 , is multiplied in front of g(τ). This factor transforms the Bloch sphere into an ellipse and changes trajectory of the state vector. Figure 2(d) shows a typical example of the comparison between the Bloch spheres as well as the trajectories of the state vector at the case of I 2 = I 3 and I 2 ≠ I 3 . These comparison indicates that the elliptical deformation the Bloch sphere and the trajectory of the state vector create a phase shift, i.e., ∆γ 1 = γ 1 − γ l on the prediction of γ l . Figure 2(f) shows that this additional phase shift at the case of I 2 ≈ I 3 /3 as a function of Θ. It shows that the maximum additional phase shift is appears at Θ ≈ π/3.

3.2.
Partial depletion with I 1 <I 2 Depletion of the system in scheme I occurs when I 1 and I 2 become comparable. Under this circumstance, the AGP at ω 1 plays a dominant role when I 1 <I 2 . If the initial condition switches to I 1 >I 2 , the AGP at ω 2 plays a dominant role. For convenience, we will select the condition I 1 <I 2 for consideration. Here, the full depletion occurs at I 1 = I 2 , which we will discuss in the next subsection. A typical example for the trajectories of the state vector on the Bloch sphere, i.e., a geometric representation, in this case are displayed in Fig. 3(a), and the corresponding accumulations of the geometric phases forq 1,2,3 are displayed in in Fig. 3(b). These panels indicate that the Bloch sphere remains a smooth and symmetric surface and the magnitude of the geometric phase shift is mainly accumulated on γ 1 . Figure 4 displays the phase shifts of the different waves γ 1,2,3 versus Θ with different ratios of (I 1 , I 2 , I 3 ). Figure 4(a1) shows γ 1 (Θ) upon fixing I 2 and different values of I 1 . At Θ = 0 and π/2, γ 1 = γ l . In the region of 0<Θ<π/2, γ 1 deviates from γ l if the system becomes depleted. Moreover, at a typical angle of Θ, namely, Θ c , the AGP becomes independent of the depletion.  Hence, in the region of 0<Θ<π/2, γ 1 satisfies |γ 1 |<|γ l |, Θ<Θ c , To estimate the relative deviation from the case of the two equally undepleted pumps, we define the deviation parameter (39) Table 3 shows that ∆ increases as I 1 increases, which indicates that the deviation between γ 1 and γ l = −π(1 − cos Θ) (the phase in the undepleted case with two equally pumps) increases as the depletion increases. Here, I 2 is fixed as I 2 = 0.5. (a2,b2,c2) γ 1,2,3 (Θ) for different values of I 2 and I 3 . Here, I 1 is fixed as I 1 = 0.2. Note that I 1 + I 2 + I 3 + I 4 = 1 and we select I 4 = 0 for all panels.  Figure 4(a2) shows γ 1 (Θ) with fixed values of I 1 and different proportions of (I 2 , I 3 ). The results indicate that the ratios of (I 2 , I 3 ) also influence the deviation rate between γ 1 and γ l . Finally, as mentioned above, Fig. 4(b1,b2,c1,c2) shows that the AGPs can also accumulate for pump wavesq 2,3 , but their absolute values, i.e., |γ 2,3 |, are much smaller than |γ 1 |.
3.3. Full depletion with I 1 = I 2 At the limit of I 1 = I 2 = I P , the system is fully depleted. At this limit, γ 1 = γ 2 , and the Bloch sphere becomes elongated to an asymmetric surface, where one of the poles is transferred into a cone. We find a threshold angle of Θ, Θ th , such that the AGP can be successfully accumulated only if Θ<Θ th . Typical examples of the geometric representation and of the AGP accumulation of γ 1 with Θ<Θ th and Θ>Θ th are shown in Figs. 5(a1,a2) and (b1,b2), respectively. When Θ>Θ th , the cusp effect is induced, which is similar to the behavior that was observed in three-wave mixing [15]. Under this circumstance, the AGP fails to accumulate. The numerical results in Fig. 5(c) indicates that Θ th is a function of I P . So far we considered the initial condition I 4 = 0, but if instead we start with I 1 = 0, this would mean that rather than starting from the bottom of the Bloch sphere, we start from the top of the sphere. The accumulated geometric phase at ω 4 will have an opposite sign with respect to what was calculated so far for ω 1 . We note that a similar phenomenon of sign inversion (between the signal wave, vs the idler was) was observed in a circular rotation in three wave mixing process [18].

Undepleted limit
Similarly, an undepleted case occurs at the initial condition of I 1,4 ≪ I 2,3 . Here, we consider first the undepleted case with two equally strong pumps. According to the normalization in Eq. (20), the limits of this condition are I 2 = I 3 ≈ 1/2 and I 1 = I 4 ≈ 0. Under this circumstance, q 2,3 (τ) ≈ √︁ I 2,3 e iφ 2,3 (τ) and |q 1,4 | 2 can be neglected. Hence, Eqs. (19) can be simplified as where ϕ P (τ) = ϕ 2 (τ) + ϕ 3 (τ). The following transformation is applied: Then, Eq. (40) can be changed to a spin-1/2 Bloch equation as where is a Stark shift [26], which means here that the phase matching also depends on the pump intensities I 2 and I 3 (in analogy to the shift of atomic energy levels by electric field). It is interesting to note that intensity dependent phase matching (Stark shift) is not limited to FWM processes, and is also observed in simultaneous three wave mixing processes, under the conditions of adiabatic elimination [32]. In this case there is also additional phase shift of the pumps Here, we let ϕ 2,3 (0) = 0. According to the above descriptions, in the case of ∆ s = ϕ P = 0, Eq. (42) driven by Eqs. (29)- (31) can give rise to the AGP as γ l = ∓π(1 − cos Θ) for idler and signal waves (i.e.,q 1 andq 4 ) if we adopt I 4 = 0 and I 1 = 0, respectively.
In comparison to to these the two theoretical results, γ l = ∓(1 − cos Θ), we also select the following two initial conditions for the numerical calculation: Obviously, these two initial conditions satisfy the undepleted condition. In the numerical simulation, we use Eq. (24) to define all the Kerr coefficients in Eqs. (19). Typical examples of the geometric representation and the accumulations of the geometric phase for the dynamics being initialed from these two conditions are displayed in Figs. 6(a1,a2) and (b1,b2), respectively. They show that the geometric phases calculated via Eqs. (33) and (34) are perfectly in accordance with each other, which demonstrates the validity of the calculation. Hence, we plot the results of the AGPs for the idler and signal (i.e., γ 1,4 ) versus Θ for these two conditions in Fig. 6(c). Unlike the AGP in the FWM process of scheme I, which can perfectly match the theoretical prediction at the undepleted limit with equal pump waves, the AGPs in the FWM process of scheme II show little deviation from the theoretical prediction. If we eliminate all the Kerr coefficients, i.e., if we set s jk = 0 (j, k = 1, 2, 3, 4), such deviations can be removed and γ 1,4 can perfectly match the theoretical prediction, γ l = ∓π(1 + cos Θ), again. It demonstrates that these deviations originate from the nonzero values of ∆ s and ϕ P , which are owing to the Stark shift.. First, we will consider the accumulation of the AGP of the idler rotated from the bottom of the Bloch sphere. This process is termed SFG because the initial signal is zero (i.e., I 4 = 0). and the three waves at ω 1,2,3 sum up to generate the higher frequency wave at ω 4 . When the magnitudes of the intensities of the other three inputs, i.e., I 1,2,3 , are comparable to each other, the system becomes depleted. For simplification, we only study the following two cases:

FWM with depletion
In the first case [cf. Equation (47)], we can focus mainly on γ 1 because the AGP accumulates primarily at the idler, and the contributions from the two pumps (i.e., γ 2 and γ 3 ) can be neglected.
In the second case [cf. Equation (48)], γ 1 = γ 2 , while γ 3 can be neglected. For these two cases, the limits of full depletion both appear at I 1 = I 2 = I 3 = 1/3, and the domain under consideration is 0<I 1 ≤ 1/3. If I 1 >1/3, the AGP accumulates mainly on one of the pump wave, which owns a smaller initial intensity. In these two cases, the depletion condition makes the curves of γ 1 (Θ) deviate from their counterparts at the undepleted limit. An example for the geometric representation for the case in Eq. (47) is displayed in Fig. 7(a1). It shows that the Bloch sphere is elongated to an asymmetric egg-shape surface when the system becomes depleted. Figure 7(a2) displays some examples for γ 1 (Θ) with different values of I 1 for the case of I 1 <I 2 = I 3 and shows that the degree of deviation increases as the level of depletion increases. Note that the level of depletion depends on the value of I 1 (here, 0 ≤ I 1 ≤ 1/3, where I 1 → 0 is the undepleted limit and I 1 = 1/3 is the fully depleted limit). For the case of I 1 <I 2 = I 3 , when I 1 >0.2, a threshold value of Θ appears, i.e., Θ th . When Θ>Θ th , the cusp effect is induced, as the AGP fails to accumulate under these circumstances. A typical example of a geometric representation for such a cusp effect is displayed in Fig. 7(b1). It shows that the spiral lines appear after the state vector passes the cusp, which indicates that the state vector cannot follow the motion of the vector of the QPM parameters after they pass through the cusp of the sphere. Figure 7(b) displays Θ th versus I 1 in the domain of I 1 ∈ (0.2, 1/3]. In contrast, for the case of (48)], Θ th exists for the entire domain of I 1 (i.e., I 1 ∈ (0, 1/3]). The cusp effect is induced if Θ>Θ th . A typical example of the geometric representation for the cusp effect in this case is shown in Fig. 7(c1) and Θ th versus I 1 is displayed in Fig. 7(c2). In this subsection, we consider the motion of the state vector from the top of the Bloch sphere. In this case, we need to set one value of I 1,2,3 to zero (here, for convenience, we set the initial intensity of the idler to zero, i.e., I 1 = 0). Under this circumstance, the system is in the process of difference-frequency generation, since a low frequency signal at ω 1 = ω 4 − (ω 3 + ω 2 ) is generated. The depletion for this process is induced when the initial intensity of the signal becomes comparable to the intensities of the pumps. Similar to the case of three-wave mixing, (b2) Θ th , beyond which the cusp effect may be induced, versus I 1 for the case of I 1 <I 2 = I 3 . This panel indicates that when I 1 >0.2, the cusp effect may be induced when Θ>Θ th . (c2) Θ th versus I 1 for the case of I 1 = I 2 <I 3 . Points A, B, and C in panels (a2), (b2) and (c2) are the position of the cases in panels (a1), (b1) and (c1), respectively. In panels (b1,c1), cusp effects are induced. In this case, spiral lines appear after the state vector passes the cusp, which indicates that the state vector cannot follow the motion of the vector of the QPM parameters after they pass through the cusp of the sphere. once the initial intensity is set at the top of the Bloch sphere, the AGP accumulates only at the signal wave.
Similarly, we also consider two cases in this subsection: This study finds that the AGP of the signal can be acquired in the domain of 0<I 4 ≤ 1/3. Figures 8(a) and (b) display the curves for γ 4 (Θ) with different values of I 4 for these two cases. Similar to the case study in SFG, the depletion creates a deviation of the AGP from the undepleted limit. Here, we find that this depletion causes γ 4 to acquire a larger phase magnitude than the undepleted counterpart when Θ< ∼ 0.44π. Typical examples of the 3 types of trajectories on the Bloch sphere are shown in Fig. 9(a1), and the calculations to this γ 4 by Eqs. (33) and (34) are shown in Fig. 9(a2). These two panels demonstrate the validation of the related calculations. However, when I 4 >1/3, the state vector W cannot immediately follow the motion of the QPM vector Q, therefore an angle is created between them. This angle induces a precession of W around Q, which can be approximately described by Obviously, if W can immediately follow the motion of Q, the angle ∼ 0, which yields W × Q ∼ 0, and the precession of W around Q is suppressed. This is what we observed at the case of  Here, for all the panels, we select Θ = π/3. I 4 <1/3. On the other hand, vector Q itself performs a circular rotation (or round-trip motion) and drags the vector W follow with it. Hence, when I 4 >1/3, the combination of these two types of motion on W (rotating around Q and following the motion of Q) draws a spiral trajectory on the Bloch sphere. A typical example of this phenomenon is displayed in Fig. 9(b1). Moreover, this spiral motion creates strong oscillations in the accumulation of γ 4 , as illustrated in Fig. 9(b2).
In this panel, it shows that the amplitude of the oscillation of γ 4 gradually increases along τ, which indicates that the angle between W and Q is also gradually increase with τ. Finally, if I 4 continues to increase, the magnitude of the angle between W and Q is also increased. The spiral motion of W becomes very strong and eventually destroys all 3 types of motion on the Bloch sphere. Hence, the accumulation of γ 4 is also completely destroyed under this circumstance [see an example in Figs. 9(c1,c2)].

Conclusion
In this paper, we study the accumulation of the adiabatic geometric phase through quasi-phase matching (QPM) for two schemes of a nonlinear four-wave mixing process. A geometric representation and a canonical Hamiltonian equation are built to describe the dynamics of these two processes. The parameters of the Hamiltonian are presented by three QPM parameters, namely, the period, phase, and duty cycle of the nonlinear modulation pattern. Furthermore, the geometric representation maps all the states of the system onto a close surface. The circular rotation of these QPM parameters is characterized by an angle, Θ. Under adiabatic conditions, the state vector follows the rotation of the QPM vector and draws a circular trajectory on the state surface. The AGPs of these waves are equal to the half of the difference of the total phase between the clockwise and the counter-clockwise rotation of the vector Q. An identical result is obtained by the difference between the total phase accumulated along the circular trajectory and that accumulated along its vertical projection. We note that, in this respect, the calculation of the geometric phase follows the same rules as the case of three-wave mixing [19]. We studied two different schemes of four wave mixing, and the main results are summarized in Tables 1 and 2. Under the undepleted condition, where the two pump waves are much stronger than the idler and signal waves, the system of scheme I (ω 1 + ω 2 = ω 3 + ω 4 ) is completely in accordance with the spin-1/2 dynamics, and the AGPs for the idler and signal can be predicted analytically. However, for the system of scheme II (ω 1 + ω 2 + ω 3 = ω 4 ), the existence of the Stark shift and the phase shift of the pump waves induce a deviation in the AGP from the theoretical result of the spin-1/2 system. In the depleted case, where the idler and signal waves become comparable to the pump waves, the magnitudes of the AGPs for the idler and signal deviate from their counterparts at undepleted limits. In addition, the deviation ratio increases as the level of depletion increases. In the case of full depletion, the threshold angle of Θ is found, beyond which the system becomes unstable and the accumulation of the AGP fails. The ability to control the geometric phase can find interesting applications in photonic devices. One possible application is that of an all-optical phased array device for fast steering and focusing, where the array elements are a set of parallel waveguides, and the phase of of each waveguide can be all-optically controlled by varying the pump intensities.