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

### Phase-engineered Raman-resonant four-wave-mixing process

First, we describe the theoretical framework by which the high-order 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, \({\rho }_{{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, \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\)^{28,29,30}. This optical process simultaneously generates a substantial magnitude of high-order Rr-FWM radiations (high-order Stokes, \({\Omega }_{-q}^{{{{{{\rm{T}}}}}}}\), and anti-Stokes, \({\Omega }_{+q}^{{{{{{\rm{T}}}}}}}\), modes) coaxially along the incident laser beam, \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\), 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 *q*th Raman mode, \({\Omega }_{q}^{{{{{{\rm{T}}}}}}}\).

$$\frac{\partial \sqrt{{n}_{q}^{{{{{{\rm{T}}}}}}}}}{\partial z}=\, \frac{N{{\hslash }}{{|}}{\rho }_{{ab}}{{|}}}{{\varepsilon }_{0}c}\left[{d}_{q-1}\sqrt{{\omega }_{q-1}^{{{{{{\rm{T}}}}}}}{\omega }_{q}^{{{{{{\rm{T}}}}}}}}\,{{\sin }}\left(\triangle {\phi }_{q}^{{{{{{\rm{T}}}}}}}\right)\sqrt{{n}_{q-1}^{{{{{{\rm{T}}}}}}}}\right.\\ \left. -{d}_{q}^{* }\sqrt{{\omega }_{q}^{{{{{{\rm{T}}}}}}}{\omega }_{q+1}^{{{{{{\rm{T}}}}}}}}\,{{\sin }}\left({\triangle \phi }_{q+1}^{{{{{{\rm{T}}}}}}}\right)\sqrt{{n}_{q+1}^{{{{{{\rm{T}}}}}}}}\right]$$

(1)

To clearly visualize the role of the relative phase, \({\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\), in this nonlinear optical process, we used the following expressions:

$${{{\mbox{E}}}}_{q}^{{{{{{\rm{T}}}}}}}=\frac{\sqrt{{{\hslash }}{\omega }_{q}^{{{{{{\rm{T}}}}}}}}}{{\varepsilon }_{0}}\sqrt{{n}_{q}^{{{{{{\rm{T}}}}}}}}\,{{\exp }}\left(i{\phi }_{q}^{{{{{{\rm{T}}}}}}}\right)$$

(2)

$${\rho }_{{ab}}=\left|{\rho }_{{ab}}\right|{{\exp }}\left(i{\phi }_{\rho }\right)$$

(3)

where \({\phi }_{q}^{{{{{{\rm{T}}}}}}}\) and \({\phi }_{\rho }\) represent the phases of the electric-field amplitude at the *q*th mode, \({E}_{q}^{{{{{{\rm{T}}}}}}}\), and the Raman coherence, \({\rho }_{{ab}}\), respectively. Then, the relative phase, \({\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\), is defined as:

$${\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\equiv {\phi }_{q}^{{{{{{\rm{T}}}}}}}-{\phi }_{q-1}^{{{{{{\rm{T}}}}}}}+{\phi }_{\rho }$$

(4)

Other terms \({n}_{q}^{{{{{{\rm{T}}}}}}}\) and \({\omega }_{q}^{{{{{{\rm{T}}}}}}}\) denote the photon number density and the angular frequency at \({\Omega }_{q}^{{{{{{\rm{T}}}}}}}\), respectively; \({d}_{q}^{{{{{{\rm{T}}}}}}}\) is a coupling coefficient between modes \({\Omega }_{q}^{{{{{{\rm{T}}}}}}}\) and \({\Omega }_{q+1}^{{{{{{\rm{T}}}}}}}\); \(z\) is the interaction length; *N* is the medium density; \(\hslash\) is the reduced Planck constant; \({\varepsilon }_{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 \({\Omega }_{q-1}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{q}^{{{{{{\rm{T}}}}}}}\) and the second term from \({\Omega }_{q+1}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{q}^{{{{{{\rm{T}}}}}}}\). According to the signs of the relative phases, \({\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\) and \({\Delta \phi }_{q+1}^{{{{{{\rm{T}}}}}}}\), 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 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, \({\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\) (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, \({\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\). 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, \({\nu }^{{\prime} {\prime} }=1\leftarrow {\nu }^{{\prime} }=0\) (124.7460 THz). We generated a pair of single-frequency nanosecond pulsed-laser radiations at 801.0817 nm (\({\Omega }_{0}\)) and 1201.6375 nm (\({\Omega }_{-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} (\(\delta\) = 0.5 GHz). We used the second harmonic, \({\Omega }_{0}^{T}\) (400.5408 nm, 6 ns), of the driving laser radiation, \({\Omega }_{0}\), as the initiation laser for the Rr-FWM process. The initiation laser radiation, \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\), was temporally and spatially overlapped with the driving laser radiations, \({\Omega }_{0}\) and \({\Omega }_{-1}\), and introduced into the medium, where the polarization of \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\) was set orthogonal to those of \({\Omega }_{0}\) and \({\Omega }_{-1}\). A series of high-order Stokes and anti-Stokes radiations, \({\Omega }_{-2}^{T}\) (600.8188 nm), \({\Omega }_{-1}^{T}\) (480.6514 nm), \({\Omega }_{+1}^{T}\) (343.3195 nm), \({\Omega }_{+2}^{T}\) (300.4038 nm), \({\Omega }_{+3}^{T}\) (267.0250 nm), was generated via the vibrational coherence, \({\rho }_{ab}\), and was manipulated in a variety of ways by controlling the phase relationships among the high-order Raman modes, \({\Omega }_{q}^{T}\). Then, such high-order Raman radiations, \({\Omega }_{q}^{T}\), 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, \({\Omega }_{\pm 1}^{{{{{{\rm{T}}}}}}}\), \({\Omega }_{\pm 2}^{{{{{{\rm{T}}}}}}}\), …, 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, \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) (Fig. 2b), and anti-Stokes, \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\) (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}^{\circ }\), scanning range of \(\pm {7.5}^{\circ }\)). In these plots, 1 represents unity quantum conversion efficiency.

Enhancement and suppression of the mode densities at \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) and \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\) 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 near-maximal 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, \({\Omega }_{{{{{{\rm{q}}}}}}}^{{{{{{\rm{T}}}}}}}\), 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, \({\Delta \phi }_{0}^{{{{{{\rm{T}}}}}}}\) or \({\Delta \phi }_{1}^{{{{{{\rm{T}}}}}}}\) (see Eq. 1), and the irregular appearances of the near-maximal densities were caused by more details in which more than three relative phases, −, \({\Delta \phi }_{-1}^{{{{{{\rm{T}}}}}}}\), \({\Delta \phi }_{0}^{{{{{{\rm{T}}}}}}}\), \({\Delta \phi }_{1}^{{{{{{\rm{T}}}}}}}\), − in the case of the Raman mode, \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\), 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, \({\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\), 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 \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) was found within the explored range. In Fig. 2k–m, we illustrate the flows of the relevant photon number densities provided by \({\Delta \phi }_{q}^{{{{{{\rm{T}}}}}}}\), 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, \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\). Here, the relative phases, \({\Delta \phi }_{0}^{{{{{{\rm{T}}}}}}}\) and \({\Delta \phi }_{1}^{{{{{{\rm{T}}}}}}}\), were manipulated to have negative signs with large amplitudes; therefore, the photon flows were accelerated from \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) and from \({\Omega }_{1}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\), respectively. Although \({\Delta \phi }_{-1}^{{{{{{\rm{T}}}}}}}\) was not optimally manipulated for the purpose of maximizing the mode density at \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\), it was manipulated to have a negligible photon flow from \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\). The photon flows from \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{+3}^{{{{{{\rm{T}}}}}}}\) and manipulations of \({\Delta \phi }_{2}^{{{{{{\rm{T}}}}}}}\) and \({\Delta \phi }_{3}^{{{{{{\rm{T}}}}}}}\) 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, \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\), and the second anti-Stokes mode, \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\). 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 \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\)) was maximized at the first Stokes mode, \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\), in Fig. 3b or the first anti-Stokes mode, \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\), 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, \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\), from the optimal solutions in the first layer and maximized each of the four Raman mode densities, \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\) in Fig. 3d, \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) in Fig. 3e, \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\) in Fig. 3f, and \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\) in Fig. 3g. In Fig. 3d, we substantially enlarged the magnitude of the photon flow to the mode \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\) from \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) (maximized in Fig. 3b). We also simultaneously enhanced the photon flows to \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) from \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\) and to \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\) from \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\). Conversely, in Fig. 3g, the directions of such photon flows were manipulated to be the opposite; that is, we enlarged the photon flows from \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\), \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\), and \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) to \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\). 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: \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\) in Fig. 3h, \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\) in Fig. 3i, \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\) in Fig. 3j, and \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\) 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: \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\): 0.22\(\pm 0.01\) (no manipulation (nom: 0.047), \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\): 0.50\(\pm 0.03\) (nom: 0.18), \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\): 0.50\(\pm 0.02\) (nom: 0.13), and \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\): 0.29\(\pm 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 no-manipulation 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: \({\Omega }_{-2}^{{{{{{\rm{T}}}}}}}\): 0.19\(\pm 0.02\), \({\Omega }_{-1}^{{{{{{\rm{T}}}}}}}\): 0.19\(\pm 0.02\), \({\Omega }_{0}^{{{{{{\rm{T}}}}}}}\): 0.22\(\pm 0.01\), \({\Omega }_{+1}^{{{{{{\rm{T}}}}}}}\): 0.17\(\pm 0.01\), \({\Omega }_{+2}^{{{{{{\rm{T}}}}}}}\): 0.14\(\pm 0.01\), and \({\Omega }_{+3}^{{{{{{\rm{T}}}}}}}\): 0.10\(\pm 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 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 light-matter 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}.