Engineering nonlinear optical phenomena by arbitrarily manipulating the phase relationships among the relevant optical fields

Nonlinear optical processes are intrinsically dominated by the phase relationships among the relevant electromagnetic fields, including the phase of nonlinear polarization produced in them. If one can arbitrarily manipulate these phase relationships at a variety of desired interaction lengths, direct and highly designable manipulations for the nonlinear optical phenomenon could be achieved. Here, we report a proof-of-principle experiment in which a high-order Raman-resonant four-wave-mixing process is used as a representative nonlinear optical process and is tailored to a variety of targets by implementing such arbitrary manipulations of the phase relationships in the nonlinear optical process. We show that the output energy is accumulated to a specific, intentionally selected Raman mode on demand; and at the opposite extreme, we can also distribute the output energy equally over broad high-order Raman modes in the form of a frequency comb. This concept in nonlinear optical processes enables an attractive optical technology: a single-frequency tunable laser broadly covering the vacuum ultraviolet region, which will pave the way to frontiers in atomic-molecular-optical physics in the vacuum ultraviolet region. Adaptive optics permits control of linear and nonlinear optical phenomena in order to achieve the desired output signal. Here, arbitrary manipulation of phase relations are used to engineer nonlinear interactions in a higher-order Raman-resonant four-wave-mixing platform.

S ince the field of nonlinear optics was established by Bloembergen et al. 1,2 , as nonlinear optical processes are intrinsically inefficient, one of the core challenges within this field has been how nonlinear optical phenomena excited at each point in space can be coherently accumulated (phase-matched) over a long interaction length. Previous research has established a variety of phase-matching technologies, including use of the birefringence properties of crystals 3,4 and implementation of periodic structure in a medium (quasi-phasematching, QPM) 1,2,5,6 . These achievements have defined a route for applying nonlinear optics to practical uses, including industrial applications such as laser processing and biomedical lasers.
Over the past two decades, a variety of methods, in which engineered functions are incorporated into nonlinear optical processes, have been investigated to begin this new chapter of nonlinear optics. Such methods include QPM techniques that incorporate a variety of periodic [7][8][9] or non-periodic structures [10][11][12][13] for operating more than two nonlinear optical processes within a single device, and gas-filled hollow-core fibers or photonic crystal fibers in which refractive index dispersions are designed to enhance a specific wavelength region in highharmonic 14,15 or supercontinuum generation 16 . Also, an angledbeam geometry in two dimensions has been numerically investigated to manipulate a broad Raman generation in solid hydrogen 17 . Three-dimensional photonic crystals including a metamaterial as a constituent element have also been investigated to engineer nonlinearly converted radiations into a variety of specific emitted directions in which the reciprocal lattice vectors assist in phase matching 18,19 . Furthermore, other techniques involve zero-index materials that enable phase-matchingfree nonlinear optical processes in all directions, demonstrating simultaneous forward and backward four-wave mixings [20][21][22] . The common conceptual idea among all these methods is that each realizes a specific function by implementing an additional design in a coherent accumulation (phase-matching) process of nonlinear optical phenomena (Fig. 1a).
In contrast, there is another approach to incorporate engineered functions into nonlinear optical processes in which we directly manipulate the scheme of the nonlinear optical process itself by manipulating the phase relationships among the relevant electromagnetic fields (Fig. 1b). This approach has not been extensively studied so far, although it has been described in the expression of the nonlinear optical process itself since the birth of nonlinear optics that nonlinear optical processes are intrinsically affected by the phase relationships among the relevant electromagnetic fields. This approach is more direct and designable in respect of engineering nonlinear optical processes. In fact, on the basis of the concept in Fig. 1b, Zheng and Katsuragawa 23 theoretically and numerically showed that a high-order Ramanresonant four-wave-mixing (Rr-FWM) process was tailored to a variety of targets, leading to a single-frequency tunable laser entirely covering the vacuum ultraviolet (VUV) region from 120 to 200 nm as being one of the attractive applications. Furthermore, Ohae et al. 24 experimentally demonstrated that the Rr-FWM process in gaseous para-hydrogen was tailored on the basis of this concept.
Here, we again focus on the Rr-FWM process as a representative nonlinear optical process in investigating the concept in Fig. 1b. Although Ohae et al. 24 demonstrated the concept, the efficiencies of the Rr-FWM generation were on the order of 10 −3 , thus limiting the Rr-FWM process to first Stokes and first anti-Stokes generation, as the experiment was executed at room temperature (non-optimal) owing to a technical limitation. To overcome this limitation, we develop a Raman cell system (Fig. 1c) that can arbitrarily manipulate the phase relationships among the relevant electromagnetic fields at a variety of interaction lengths in a nonlinear optical medium held at the liquid nitrogen temperature (optimal). By realizing this Raman cell system, we experimentally tailor the Rr-FWM process with efficiencies of tens of percentages to a variety of targets; namely, we demonstrate accumulation of the output energy to a selected single-Raman mode on demand; and at the opposite extreme, an equal distribution of output energy over high-order Raman modes in the form of a frequency comb. These results are also well reproduced in a numerical calculation and are discussed in respect of the physical mechanism.

Results and discussion
Phase-engineered Raman-resonant four-wave-mixing process. First, we describe the theoretical framework by which the highorder Rr-FWM process can be engineered by implementing the relative phase manipulations during its evolution. Figure 1d, e represents a Rr-FWM process. In this process, a high-Raman coherence, ρ ab , is adiabatically driven by precisely controlling the two-photon detuning, δ, where a pair of long-pulse laser fields at the two wavelengths, Ω 0 and Ω −1 are applied [25][26][27][28][29] . The produced molecular ensemble with high-Raman-coherence functions as an ultra-high-frequency phase modulator and thereby deeply modulates a variety of arbitrary incident laser radiations, Ω T 0 28-30 . This optical process simultaneously generates a substantial magnitude of high-order Rr-FWM radiations (high-order Stokes, Ω T Àq , and anti-Stokes, Ω T þq , modes) coaxially along the incident laser beam, Ω T 0 , without being restricted by angle phase matching.
Maxwell-Bloch equation. A standard framework for the Maxwell-Bloch equation well describes this Rr-FWM process (see Methods). Equation 1 represents one of the coupled propagation equations, expressing the generation of the qth Raman mode, To clearly visualize the role of the relative phase, Δϕ T q , in this nonlinear optical process, we used the following expressions: where ϕ T q and ϕ ρ represent the phases of the electric-field amplitude at the qth mode, E T q , and the Raman coherence, ρ ab , respectively. Then, the relative phase, Δϕ T q , is defined as: Other terms n T q and ω T q denote the photon number density and the angular frequency at Ω T q , respectively; d T q is a coupling coefficient between modes Ω T q and Ω T qþ1 ; z is the interaction length; N is the medium density; _ is the reduced Planck constant; ε 0 is the vacuum permittivity; c is the speed of light.
The first term on the right-hand side of Eq. 1 implies energy flows of electromagnetic fields from Ω T qÀ1 to Ω T q and the second term from Ω T qþ1 to Ω T q . According to the signs of the relative phases, Δϕ T q and Δϕ T qþ1 , directions of such energy flow change, and their speeds vary depending on their values over a full dynamic range of −π to +π. Therefore, we can tailor the obtained nonlinear optical phenomenon toward a variety of targets by G la n la s e r p r is m  I ) Fig. 1 Conceptual illustrations of the engineered Raman-resonant four-wave-mixing (Rr-FWM) process. Conceptual comparison between engineerings in nonlinear optical processes: a manipulation of the spatial phases, k(ω, r). b Manipulation of relative phases, ϕ(ω, r). c Schematics: the apparatus is a gas cell filled with gaseous para-hydrogen and cooled to 77 K by liquid nitrogen. Six dispersive plates (fused silica) are inserted in the interaction region, where the relative phases, Δϕ T q , are manipulated by precise control of the effective optical thicknesses of the plates. d Scheme of the adiabatic excitation of the vibrational coherence, ρ ab , at a pure vibrational Raman transition of ν 00 ¼ 1 j i j ν 0 ¼ 0i. e Scheme of the high-order Rr-FWM process initiated by laser radiation, Ω T 0 . The relative phases, Δϕ T q , dominate the photon flow between the adjacent Raman modes, Ω T q and Ω T qÀ1 . f-i Expected high-order Stokes and anti-Stokes generations in the engineered Rr-FWM processes by manipulating the relative phases, Δϕ T q ; f no manipulation applied; g, h output energy accumulations onto a specific Raman mode; and i broad comb-like generation.
providing the freedom of arbitrarily manipulating the directions of these energy flows including their flow speeds at a variety of interaction lengths desired in the evolution of this phenomenon ( Fig. 1f-i).
On-demand phase manipulation among many discrete spectral modes. To implement such phase manipulations in a nonlinear optical medium, it is important to determine how to practically realize such a physical mechanism that can simultaneously manipulate many relevant phases among the high-order Raman modes to arbitrary values. Previous work showed that a simple device with which the optical thickness of a transparent dispersive plate can be precisely tuned over a relatively large plate thickness can manipulate the relative phases among many discrete spectral modes nearly on demand [31][32][33][34] . This method includes the solution of the high-order temporal Talbot method as one of the many solutions 35 . Here, we use this tunable plate thickness technology to arbitrarily manipulate the concerned relative phases, Δϕ T q (Fig. 1c).
Experimental. We developed a Raman cell where six transparent dispersive plates (fused silica with a thickness of 5 mm) were installed in a nonlinear optical medium that was refrigerated down to a liquid nitrogen temperature of 77 K. The effective optical thickness of each plate was electronically controlled by adjusting the incident angle from outside the cell (Fig. 1c) to arbitrarily manipulate the relative phases, Δϕ T q . We used gaseous para-hydrogen (purity > 99.9%) as a nonlinear optical medium with a density controlled at 4.2 × 10 19 cm −3 at a temperature of 77 K (Lamb-Dicke regime; experimentally found), providing the smallest inhomogeneous broadening (<200 MHz) at the pure vibrational Raman transition, ν 00 ¼ 1 ν 0 ¼ 0 (124.7460 THz). We generated a pair of single-frequency nanosecond pulsed-laser radiations at 801.0817 nm (Ω 0 ) and 1201.6375 nm (Ω À1 ) 33 that overlapped in time and space (pulse duration of 6 ns), softly focused them into the medium (waist diameter of 200 μm), and then adiabatically drove a vibrational coherence, ρ ab (δ = 0.5 GHz). We used the second harmonic, Ω T 0 (400.5408 nm, 6 ns), of the driving laser radiation, Ω 0 , as the initiation laser for the Rr-FWM process. The initiation laser radiation, Ω T 0 , was temporally and spatially overlapped with the driving laser radiations, Ω 0 and Ω À1 , and introduced into the medium, where the polarization of Ω T 0 was set orthogonal to those of Ω 0 and Ω À1 . A series of highorder Stokes and anti-Stokes radiations, Ω T À2 (600.8188 nm), Ω T À1 (480.6514 nm), Ω T þ1 (343.3195 nm), Ω T þ2 (300.4038 nm), Ω T þ3 (267.0250 nm), was generated via the vibrational coherence, ρ ab , and was manipulated in a variety of ways by controlling the phase relationships among the high-order Raman modes, Ω T q . Then, such high-order Raman radiations, Ω T q , were picked out at appropriate interaction lengths in the series evolution process and detected by photodiodes with relatively calibrated sensitivities.
Physical mechanism in engineering the Rr-FWM process. A high vibrational coherence, ρ ab , was adiabatically driven, and thereby, a series of high-order Stokes and anti-Stokes radiations, Ω T ± 1 , Ω T ± 2 , …, were generated coaxially without being restricted by angle phase matching.
First, we studied how the physical mechanism of the relative phase manipulation implemented in this nonlinear optical process functioned in reality by comparing the resultant Stokes and anti-Stokes generations with those calculated in the numerical simulations. For this purpose, we focused on the generation of the first Stokes and anti-Stokes modes, where only two dispersive plates (A 1 and B 1 in Fig. 2a) were used. Figure 2b, c shows contour plots of the photon number densities of the generated first Stokes, Ω T À1 (Fig. 2b), and anti-Stokes, Ω T þ1 (Fig. 2c), modes as a function of the angles of the two inserted dispersive plates, A 1 and B 1 , corresponding to manipulations of the effective optical thicknesses (resolution of 0:015 , scanning range of ± 7:5 ). In these plots, 1 represents unity quantum conversion efficiency. Enhancement and suppression of the mode densities at Ω T À1 and Ω T þ1 occurred almost periodically, while on the other hand, their near-maximal densities appeared irregularly. These observed features were well reproduced in the corresponding numerical simulations (Fig. 2d, e), including the appearance rates of the nearmaximal densities. That is, the numerical simulations well reproduced the intrinsic physical nature of this engineered nonlinear optical process, although it did not perfectly reproduce all properties, including their absolute values, which should depend on the initial phases and the initial absolute plate thicknesses.
Although it is difficult to measure the phase relationships among the Raman modes, Ω T q , in the experiment, especially during the evolution of the nonlinear optical process, we can discern such phase relationships in detail from the numerical calculation, as the numerical simulations well reproduced the intrinsic physical nature of the engineered Rr-FWM process (Fig. 2f, g). The periodic behaviors of the observed mode densities were due mainly to the inversion of the signs of the relative phases, Δϕ T 0 or Δϕ T 1 (see Eq. 1), and the irregular appearances of the near-maximal densities were caused by more details in which more than three relative phases, −, Δϕ T À1 , Δϕ T 0 , Δϕ T 1 , − in the case of the Raman mode, Ω T À1 , were non-negligibly associated with different sign-inversion periodicities with different magnitudes which determine the speeds of the relevant photon flows (see Supplementary Fig. S1 for the details including the high-order Raman modes).
In Fig. 2h-j, we show how the relative-phase relationships among the high-order Raman modes, Δϕ T q , were manipulated along the interaction length at the point marked by the white cross in Fig. 2d, which denotes that the maximal mode density at Ω T À1 was found within the explored range. In Fig. 2k-m, we illustrate the flows of the relevant photon number densities provided by Δϕ T q , revealing that such phase relationships were successfully manipulated at interaction lengths Z 1 and Z 2 so that the photons were accumulated to form the first Stokes mode, Ω T À1 . Here, the relative phases, Δϕ T 0 and Δϕ T 1 , were manipulated to have negative signs with large amplitudes; therefore, the photon flows were accelerated from Ω T 0 to Ω T À1 and from Ω T 1 to Ω T 0 , respectively. Although Δϕ T À1 was not optimally manipulated for the purpose of maximizing the mode density at Ω T À1 , it was manipulated to have a negligible photon flow from Ω T À1 to Ω T À2 . The photon flows from Ω T þ2 to Ω T þ3 and manipulations of Δϕ T 2 and Δϕ T 3 were ignored here, as their mode densities were nearly zero.
Up to this point, we confirmed how the physical mechanism of the relative phase manipulations implemented in nonlinear optical processes can work in reality by comparing the fundamental experiment with the detailed numerical simulation.

Application of the engineered Rr-FWM process
On-demand enhancement of a selected Raman mode. The concept of relative-phase manipulation in nonlinear optical processes has generality; we can apply it to a variety of purposes 23 . In Fig. 3a-p, we applied such a physical principle to the targets where the output energy in the Rr-FWM process was accumulated to a specific Raman mode selected from the range between the second Stokes mode, Ω T À2 , and the second anti-Stokes mode, Ω T þ2 . For this purpose, we implemented engineered phase manipulations comprising three layers in the Rr-FWM process, each having a pair of transparent dispersive plates, along the interaction length as illustrated in Fig. 3a. The contour plots in Fig. 3b, c are the same as those in Fig. 2b, c. The white cross in each figure indicates the point at which the photon number density (normalized to that at Ω T 0 ) was maximized at the first Stokes mode, Ω T À1 , in Fig. 3b or the first anti-Stokes mode, Ω T þ1 , in Fig. 3c in the explored range. We fixed the phase manipulations at these conditions as the optimal solutions in the first layer for the targets aimed at as above.
Next, we moved to the second layer ( Fig. 3d-g). Here, also by using two dispersive plates, A 2 and B 2 , we engineered photon flows to extend them to each of the four Raman modes, Ω T À2 to Ω T þ2 , from the optimal solutions in the first layer and maximized each of the four Raman mode densities, Ω T À2 in Fig. 3d, Ω T À1 in Fig. 3e, Ω T þ1 in Fig. 3f, and Ω T þ2 in Fig. 3g. In Fig. 3d, we substantially enlarged the magnitude of the photon flow to the mode Ω T À2 from Ω T À1 (maximized in Fig. 3b). We also simultaneously enhanced the photon flows to Ω T À1 from Ω T 0 and to Ω T 0 from Ω T þ1 . Conversely, in Fig. 3g, the directions of such photon flows were manipulated to be the opposite; that is, we enlarged the photon flows from Ω T þ1 to Ω T þ2 , Ω T 0 to Ω T þ1 , and Ω T

À1
to Ω T 0 . In Fig. 3e, f, we again implemented the same manipulations as in the first layer. The white crosses in Fig. 3d-g indicate the conditions that maximized the photon number densities for each of the four Raman modes in the explored range, which were then fixed as the optimal solutions at the second layer.
Finally, in the third layer, we repeated the same conceptual phase manipulations as those implemented in the second layer, with the two dispersive plates, A 3 and B 3 , then, we completed this output energy accumulation onto each of the four Raman modes in the Rr-FWM process. The white crosses in Fig. 3h-k indicate the optimal solutions, which indicated the conditions for obtaining maximal photon number densities in the final layer: Ω T À2 in Fig. 3h, Ω T À1 in Fig. 3i, Ω T þ1 in Fig. 3j, and Ω T þ2 in Fig. 3k. Figure 3m-p shows the photon number density distributions among the Raman modes that were ultimately achieved through this phase engineering process. When no manipulation was applied ( Fig. 3l; dispersive plates not inserted), this nonlinear optical process intrinsically evolved to broaden both the positive and negative high-order Raman modes. Here, by implementing the three relative-phase manipulation layers in the Rr-FWM process, we achieved significant enhancements of specific single-Raman-mode densities: Ω T À2 : 0.22 ± 0:01 (no manipulation (nom: 0.047), Ω T À1 : 0.50 ± 0:03 (nom: 0.18), Ω T þ1 : 0.50 ± 0:02 (nom: 0.13), and Ω T þ2 : 0.29 ± 0:02 (nom: 0.038). The light red bars overlaid in Fig. 3m-p show the mode densities at the targets, simulated numerically with the optimal phase manipulations.
Generation of broad comb-like Raman modes. As already described, the physical concept of the engineered nonlinear optical process described here can be applied for a variety of purposes. In Fig. 3q-t, we examined another target that is regarded as an opposite extreme against those executed in Fig. 3m-p, showing the generation of an equal photon number density distribution over broad high-order Raman modes. Compared with the nomanipulation scenario in Fig. 3l, a very flat photon number density distribution in the form of a frequency comb was realized via the same three layers of relative-phase manipulations: Ω T À2 : 0.19 ± 0:02, Ω T À1 : 0.19 ± 0:02, Ω T 0 : 0.22 ± 0:01, Ω T þ1 : 0.17 ± 0:01, Ω T þ2 : 0.14 ± 0:01, and Ω T þ3 : 0.10 ± 0:01. The produced spectrum was phase coherent in time and space 33 , having an ability to generate ultrafast pulses with a temporal duration of 1.2 fs at an ultrafast repetition rate of 125 THz (inset in Fig. 3t).
Discussion on the conceptual differences from related studies. Last, we briefly comment on the conceptual differences between to Ω T þ2 , observed after each of the three relative-phase manipulation layers (1st, b, c; 2nd, d-g; 3rd, h-k). The crosses in the contour plots indicate the optimal conditions explored in each of the three layers in each of the targets, which maximize a specific Raman mode density of Ω T À2 to Ω T þ2 . l-p Photon number density distributions among the Raman modes, Ω T q (q = −2 to 3; Ω T À2 : 600.8188 nm, Ω T À1 : 480.6514 nm, Ω T 0 : 400.5408 nm, Ω T þ1 : 343.3195 nm, Ω T þ2 : 300.4038 nm, Ω T þ3 : 267.0250 nm), which were finally achieved at the apparatus exit; for l, no manipulation was applied; for m-p, the plates were engineered towards each of the targets to maximize a specific Raman mode density (red bars) from Ω T À2 to Ω T þ2 . The overlaid light red bars show the photon number densities calculated by the numerical simulations. q-t Generation of an equal photon number density distribution over broad high-order Raman modes, Ω T À2 to Ω T þ2 , was examined. The measurements were conducted at four interaction lengths of Z 1 (q), Z 3 (r), Z 5 (s), and Z 7 (t) at the exit of the cell. This target was regarded as an opposite extreme to those executed in m-p. Each of the error bars in l-t implies a standard deviation in one hundred measurements executed for each of the Raman mode densities.
previous studies and the work we present here. In the terminology of adaptive optics, a large number of studies involving the control of nonlinear optical phenomena have been reported. The key concept in such works is the implementation of an artificial design in the distributions of amplitude, phase, or polarization of the incident light to achieve the optimal solution for a specific target, where the feedback control is applied to the light at the incidence and the lightmatter interaction is, in general, considered a black box [36][37][38][39] . Recently, Tzang et al. demonstrated the manipulation of a multimode stimulated Raman scattering cascade in a multimode fiber by applying adaptive-optic control to the wavefront of the incident light, where the angle phase-matching conditions were substantially shaped 40 . In addition, the conceptual idea of metamaterial: designing optical susceptibility tensors including their phases by creating artificial structures at the nanoscale, has been extended to nonlinear optical processes. Many works on nonlinear beam shaping including the nonlinear optical hologram have been reported [41][42][43] .

Conclusion
We conducted a proof-of-principle experiment to show that nonlinear processes can be tailored in a variety of ways by manipulating the relative phases among the relevant electromagnetic fields during the evolutions of nonlinear processes. A high-order Rr-FWM process was used as a representative nonlinear process and was tailored to a variety of targets by installing practical phase manipulation devices (thickness-controlled dispersive plates) consisting of three layers in gaseous para-hydrogen controlled in the Lamb-Dicke regime (77 K, 4.2 × 10 19 cm −3 ). We showed that a specific, intentionally selected Raman mode was enhanced on demand with a maximum photon number density of 0.5. We also demonstrated that at the opposite extreme, the photon number densities were equally distributed over the highorder Raman modes in the form of a frequency comb.
The physical concept described here is simple and has generality, paving the way to avenues for incorporating engineered functions in nonlinear optical processes. An attractive potential application of this technology could be the realization of a single-frequency tunable laser entirely covering the VUV region from 120 to 200 nm, as suggested in a previous study 23 where Ω T 0 was set at 210 nm. As the mature single-frequency tuneable solid-state laser technology in the near-infrared region established a major trend in atomic-molecularoptical (AMO) physics leading to the realization of Bose-Einstein condensation from the initiation of laser cooling, if an equivalent laser technology can be established in the VUV region, frontiers in AMO physics, such as the laser cooling of anti-matter (Lyman α: 121.6 nm) 44 or optical frequency standards locked at nuclear transitions (Th, 149 nm) 45,46 , may be deeply explored.

Methods
Lasers and operating conditions. The Ω 0 laser radiation was generated by an injection-locked nanosecond pulsed Ti:sapphire laser (repetition rate: 10 Hz), with a homemade external-cavity-controlled continuous-wave diode laser (800 nm) as the seed laser. The Ω À1 laser radiation was generated by an injection-locked optical parametric oscillator (OPO) followed by an optical parametric amplifier (OPA), with the Ω 0 laser radiation used as the pump laser radiations for both the OPO and the OPA. A homemade external-cavity-controlled continuous-wave diode laser (1200 nm) was used as the seed laser for the OPO. The third laser radiation, Ω T 0 , was generated by taking the second harmonic of the driving laser radiation, Ω 0 . Temporal overlap of the three laser radiations (Ω 0 ; Ω À1 ; and Ω T 0 ) was very stable (pulse duration: 6 ns), as all of them were provided by the single Ti:sapphire laser system. We also coaxially overlapped the three laser beams and softly focused them into the Raman cell. The beam diameters of the three laser radiations at the waist were set to 220 μm at 1 e 2 for Ω À1 and Ω 0 and to 90 μm at 1 e 2 for Ω T 0 . We set wavelengths of the two driving laser radiations to 1201.6375 nm for Ω À1 and 801.0817 nm for Ω 0 , and adiabatically drove the molecular coherence, ρ ab 25, [27][28][29][30]33 , where the frequency difference of the two driving laser radiations was set so that the two-photon detuning, δ; to the vibrational Raman transition ( b j i : ν 00 ¼ 1 j i a j i : jν 0 ¼ 0i) was optimal (0.5 GHz). The pulse energies of the two driving laser radiations were adjusted to 5.3 mJ for Ω À1 and 5.0 mJ for Ω 0 .
Theoretical framework: Maxwell-Bloch equation. The high-order Rr-FWM process can be well described by the standard framework of the Maxwell-Bloch equations. In the far-off resonant Λ-scheme, the entire system can be reduced to a two-level system with an effective Hamiltonian as: where Ω aa and Ω bb are ac Stark shifts for the ground state, jai, and the excited state, jbi, respectively, and Ω ab indicates the effective two-photon Rabi frequency. Equation 6 represents the equations of motion in the reduced two-level system with a density matrix formula: here, ρ aa and ρ bb are the populations of the ground state, jai, and the excited state, jbi, respectively, and ρ ab is coherence associated with the Raman transition between states jai and jbi. The coefficients γ a , γ b , and γ c are the decay rates of the populations ρ aa and ρ bb , and the coherence, ρ ab , respectively. The coupled propagation equation for the complex electric-field amplitude, E q (qth Raman mode), propagating in the z direction, is expressed with the slowly varying envelope approximation as: where N is the molecular density of para-hydrogen, ω q is the angular frequency at the qth Raman mode, _ is the reduced Planck constant, and ε 0 is the vacuum permittivity. Parameters, a q and b q determine the dispersions of para-hydrogen, and d q determines the coupling strength between neighboring Raman modes. To show the role of the relative phase more explicitly, we transformed Eq. 7 to Eqs. 8 and 9, where the complex electric-field amplitude and molecular coherence are transformed as E T q ¼ E T q e iϕ T q and ρ ab ¼ ρ ab e iϕ ρ , respectively, by using the expressions of photon number density, n T q / E T q 2 _ω q , and the relative phase, Δϕ T q . The relative phases, Δϕ T q and Δϕ T qþ1 are defined as Δϕ T q ϕ T q À ϕ T qÀ1 þ ϕ ρ ; ð10Þ Data availability The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
The codes that support the findings of this study are available from the corresponding author upon request.