Diffraction manipulation by four-wave mixing

We suggest a scheme to manipulate paraxial diffraction by utilizing the dependency of a four-wave mixing process on the relative angle between the light fields. A microscopic model for four-wave mixing in a Lambda-type level structure is introduced and compared to recent experimental data. We show that images with feature size as low as 10 micrometers can propagate with very little or even negative diffraction. The mechanism is completely different from that conserving the shape of spatial solitons in nonlinear media, as here diffraction is suppressed for arbitrary spatial profiles. At the same time, the gain inherent to the nonlinear process prevents loss and allows for operating at high optical depths. Our scheme does not rely on atomic motion and is thus applicable to both gaseous and solid media.


I. INTRODUCTION
The diffraction of light during propagation in free space is a fundamental and generally unavoidable physical phenomenon. Because of diffraction, light beams do not maintain their intensity distribution in the plane transverse to the propagation direction, unless belonging to a particular class of non-diffracting (Bessel) beams [1]. In nonuniform media, waveguiding is possible for specific spatial modes [2,3], or equivalently arbitrary images may revive after a certain self-imaging distance [4]. However in such waveguides, the suppression of diffraction for multimode profiles is not trivial, as each transverse mode propagates with a different propagation constant or group velocity, resulting in spatial dispersion of the profile.
Recently, a technique was suggested [5] and experimentally demonstrated [6] to manipulate (eliminate, double, or reverse) the diffraction of arbitrary images imprinted on a light beam for arbitrary propagation distances. The technique is based on electromagnetically induced transparency (EIT) [7] in a thermal atomic gas. Unlike other methods utilizing EIT [2][3][4][8][9][10][11][12][13][14][15][16][17], which rely on spatial non-uniformity, this technique operates and prescribes non-uniformity in k ⊥ space. Here, k ⊥ denotes the transverse wave vectors, i.e., the Fourier components of the envelope of the field in the transverse plane, which is the natural basis for paraxial diffraction. The technique of Refs. [5,6] relies on the diffusion of the atoms in the medium and on the resulting so-called Dicke narrowing [18,19]. Due to Dicke narrowing, the linear susceptibility becomes quadratic in k ⊥ and results in motionalinduced diffraction that can counteract the free-space diffraction. Unfortunately in the currently available experimental conditions, the resolution limit of motionalinduced diffraction is on the order of 100 µm, preventing it from being of much practical use. Higher resolution requires a denser atomic gas, in which strong absorption is unavoidable due to imperfect EIT conditions. Very recently, Zhang and Evers proposed to circumvent the absorption by generalizing the model of motional-induced diffraction to a four-wave mixing (FWM) process in combination with EIT [20]. The FWM process further allows the frequency conversion of the image and increases the available resolution [20].
In this paper, we propose a scheme to manipulate diffraction using FWM [21][22][23] without the need for motional-induced diffraction. The mechanism we study originates from phase matching in k ⊥ space and does not require a gaseous medium; it is therefore directly applicable to solid nonlinear media. For our model to be general and accommodate motional-broadening mechanisms (not important in solids), we here still concentrate on describing atomic gases and validate our model against relevant experiments. The inherent gain of the FWM process allows us to improve the spatial resolution by working with relatively higher gas densities while avoiding loss due to absorption.
In Sec. II, we introduce a microscopic model of FWM in a Λ system, based on Liouville-Maxwell equations, similar to the one used in Ref. [24]. In Sec. III, we compare the model to recent experimental results of FWM in hot vapor [24,25]. We use our model in Sec. IV to show that, with specific choice of frequencies, the k ⊥ dependency of the FWM process can be used to eliminate the diffraction of a propagating light beam. We also present a demonstration of negative diffraction, implementing a paraxial version of a negative-index lens [26], similar to the one in Ref. [6] but with positive gain and higher resolution. Finally, we analyze the resolution limitation of our scheme and propose ways to enhance it. We show that, for cold atoms at high densities (∼10 12 cm −3 ), diffraction-less propagation of an image with a resolution of ∼10 µm can be achieved.

A. Model
We consider an ensemble of three-level atoms in a Λ configuration depicted in Fig. 1a. The atomic states are denoted as |u , |l , and |r , for the up, left, and right states, and the corresponding energies are ǫ u , ǫ l , and ǫ r . We introduce the optical transition frequencies ω ul = (ǫ u − ǫ l ) /ℏ and ω ur = (ǫ u − ǫ r ) /ℏ, taken to Typeset by REVT E X be much larger than the ground-state splitting ω lr = ω ul − ω ur = (ǫ l − ǫ r ) /ℏ. To simplify the formalism, we assume the same dipole-moment for the two optical transitions µ = µ ul = µ ur , where µ αα ′ = α| µ· x |α ′ , µ being the dipole-moment operator.
The atom interacts with three external, classical electromagnetic fields, propagating in time t and space r, Here ω i are the frequencies of a weak 'probe' (i = p) and two 'control' fields (i = c); ε j are the polarization vectors (with j = p, cl, cr); k 0 i ≡ ω i /c are the wave vectors in the case of plane waves, otherwise they are carrier wavevectors; and Ω i (r, t) are the slowly varying envelopes of the the Rabi frequencies, satisfying ∂ 2 t Ω i (r, t) ≪ |ω i ∂ t Ω i (r, t)| and ∂ 2 z Ω i (r, t) ≪ k 0 i ∂ z Ω i (r, t) . We shall analyze the case of two identical control fields with the same Rabi frequency Ω c , wave vector k 0 c , and frequency ω c . The strong control and weak probe fields stimulate a weak classical 'Stokes' field (or 'conjugate') at a fre- The resonances are characterized by the one-photon detuning ∆ 1p = ω c − ω ur and the two-photon detuning ∆ 2p = ω p − ω c − ω lr (see Fig. 1a) .The population of the excited level |u decays to the ground levels |l and |r with rates Γ l and Γ r . The atomic coherence between the excited level |u and each of the ground levels |l and |r decays with rates Γ d,l and Γ d,r . For simplicity, we assume Γ l = Γ r = Γ d,l = Γ d,r ≡ Γ. Within the ground state, we consider a population relaxation with a symmetric rates Γ l↔r and decoherence with rate Γ lr .
In a frame rotating with the control frequency ω c , the equations of motion for the local density-matrix ρ (r, t) are better written in terms of the slowly-varying densitymatrix R (r, t) , where R u,j (r, t) = ρ u,j (r, t) e iωct−ik 0 c z for j = l, r and R α,α ′ (r, t) = ρ α,α ′ (r, t) for all other matrix elements, Here γ cj = Γ − i (ω c − ω uj ) are complex one-photon detunings of the control fields (j = l, r), are interference fields. Assuming non-depleted control fields, constant in time and space Ω c (r, t) = Ω c , we complete the description of the atom-field interaction with the propagation equations under the envelope approximation for the probe field (5a) and the Stokes field (5b) where ∇ 2 ⊥ ≡ ∂ 2 /∂x 2 + ∂ 2 /∂y 2 the transverse Laplacian, g = 2πN |µ| 2 q/ℏ the coupling strength proportional to the atomic density N , and q ≡ k 0 c ≈ k 0 p ≈ k 0 s . To obtain Eqs. (5), we neglected the second-order t and z derivatives of the envelopes.

B. Steady-state solution
The evolution of the fields is described by a set of non-linear, coupled differential equations for the density matrix elements R α,α ′ and the weak fields Ω p and Ω s [Eqs. (3)-(5)], which require further approximations to be solved analytically. The solutions and assumptions are detailed in the Appendix, where we find the steady state of the system to first order in the weak fields, α,α ′ being the zero-and first-order steadystate solutions. The most important assumption is the proximity to two-photon resonance, such that δ ω is on the order of the ground-state frequency splitting ω lr and much larger than any detuning, Rabi frequency, or pumping rate in the system. Plugging Eqs. (A5)-(A8) and (6) for R u,r and R u,l into the propagation equations (5) and discarding terms rotating at δ ω and 2δ ω , we obtain the well-known FWM form including paraxial diffraction, where Here α j = g (n l /γ jl + n r /γ jr ) are the linear absorption coefficients of the probe (j = p) or Stokes (j = s) fields, with n i ≡ R i,i the populations of the i = r, l levels. β p = g (n l /γ pl + n r /γ * cr ) and β s = g (n r /γ * sr + n l /γ cl ) are twophoton interaction coefficients, γ jk = Γ − i (ω j − ω uk ) [j = p, c, s; k = l, r] are complex one-photon detunings, and is the complex two-photon detuning. Eqs. (7) are similar to those obtained by Harada et al. [24] but here including the diffraction term ±i∇ 2 ⊥ /(2q), which we require in order to explore the spatial evolution of the FWM process.
We start with the simple case of a weak plane-wave Fig. 1). We assume that the generated Stokes field is also a plane wave Ω s (r) = g (z) e ik s ⊥ ·r ⊥ e i(k s z −k s 0 )z . Substituting into Eqs. (7), the phase-matching condition k s ⊥ = − k p ⊥ is easily obtained, and the resulting equations for f and g are [23] Assuming f (0) = 1 and g (0) = 0, we follow Ref. [23] and find along the medium with the eigenvalues In the limit where |B| and |C| are much smaller than |A| and |D|, the solution is governed by independent EIT for the probe and Stokes fields with little coupling between them. In the opposite limit, the fields experience strong coupling, and the real part of the eigenvalues λ 1,2 can be made positive and result in gain.

III. COMPARISON WITH EXPERIMENTS
To verify our model, we have calculated the probe transmission as a function of the two-photon detuning and compared it to the data published in Refs. [24,25]. The Doppler effect due to the motion of the thermal atoms is taken into account by averaging the FWM coefficients, Q = A, B, C, D in Eq. (8), over the Doppler profile [27]. Assuming nearly collinear beams, the mean coefficients are where v th = k B T /m, T the cell temperature, and m the atomic mass. Fig. 2 shows the transmission spectrum in (a) rubidium vapor and (b) sodium vapor (cells length L ≃ 5 cm). Our model reproduces the experimental spectra, including the Doppler-broadened absorption lines and the gain peaks for both rubidium and sodium experiments. The missing peak in Fig. 2b is due to anti-Stokes generation not included in the model.

IV. DIFFRACTION MANIPULATION BY FWM
We now concentrate on a specific choice of frequency detunings, for which the phase dependency of the FWM process can be used to manipulate the diffraction of the propagating probe and Stokes. To this end, we study the evolution of an arbitrary image F (r) imprinted on the probe beam, with the boundary conditions Ω p (r ⊥ , 0) = F (r) and Ω s (r ⊥ , 0) = 0. Our prime examples shall be the propagation of the image without diffraction or with reverse diffraction, both of which while experiencing gain.

A. Image propagation
We start by solving Eqs. (7) in the transverse Fourier basis, where Ω p/s (k ⊥ , z) = dr 2 ⊥ e −ik ⊥ ·r ⊥ Ω p/s (r ⊥ , z), and the coefficientsĀ,B,C, andD are Doppler averaged according to Eq. (12). We notice that Eqs. (13) are identical to Eqs. (9) with k ∆ = 0 and with the substitutions A →Ā − ik 2 ⊥ /(2q) andD →D + ik 2 ⊥ /(2q). The evolution of the probe and Stokes fields then follows from Eq. (10), where We choose |e λ2z | ≫ |e λ1z | and obtain Ω p (k ⊥ , z) = Ω p (k ⊥ , 0) e Z , where determines the changes in the spatial shape of the probe along its propagation. Re Z is responsible for the k ⊥dependency of the gain/absorption, and Im Z is responsible for the k ⊥ -dependency of the phase accumulation, that is, the diffraction-like evolution.

B. Suppression of paraxial diffraction
In general, in order to minimize the distortion of the probe beam, one is required to minimize both the real and imaginary k ⊥ -dependencies of Z. To better understand the behavior of Z, we expand it in orders of k 2 ⊥ . Taking the limit where 2E = [ Ā −D 2 + 4BC] 1/2 , we write and find The k ⊥ -dependency, governed by Z (2) , can be controlled through the FWM coefficientsĀ,B,C, andD given in Eq. (8), by manipulating the frequencies of the probe and control fields (ω p , ω c ), the control amplitude Ω c , and the density N . We demonstrate this procedure in Fig. 3, using for example the experimental conditions of the sodium experiment, detailed in Fig. 2. First, we observe the gain of the probe and the Stokes fields in Fig. 3a and 3b, as a function of the one-(∆ 1p ) and two-(∆ 2p ) photon detunings. The gain is achieved around the two-photon resonance (∆ 2p ≈ 0), either when the probe is at the one-photon resonance (∆ 1p ≈ 0) or the Stokes (∆ 1p ≈ ω lr , here ≈ 2 GHz); the latter exhibits higher gain, since the probe sits outside its own absorption line. The real and imaginary parts of Z (2) are plotted in Figs. 3c and 3d. When Re Z (2) = 0 (dashed line), the gain/absorption is not k ⊥dependent, whereas when Im Z (2) = 0 (solid line), the phase accumulation along the cell is not k ⊥ -dependent. When both happen, Z (2) = 0, and a probe with a spectrum confined within the resolution limit k ⊥ ≪ k 0 propagates without distortion. The exact propagation exponent Z as a function of k ⊥ for the point Z (2) = 0 (∆ 2p ≈ 0.4 MHz, ∆ 1p ≈ 0) are plotted in Fig. 3e. As expected, both real (blue solid line) and imaginary (red dashed line) parts of Z are constant for k ⊥ ≪ k 0 (deviation of 1% within k ⊥ < k 0 /2 and 0.1% within k ⊥ < k 0 /4). In the specific example of Fig. 3, the probe's gain is ∼1.4, the Stokes' gain is ∼ 4, and k 0 ≈ 40 mm −1 .
To illustrate the achievable resolution, we shall employ a conservative definition for a characteristic feature size in the image in area units a = (2π/k ⊥ ) 2 [For example: for a Gaussian beam, a 1/2 shall be twice the waist radius, and, for the field pattern E = 1 + cos(k ⊥ x) cos(k ⊥ y), the pixel area is a 2 . The Rayleigh length is qa 2 /8]. Fig. 4 presents numerical calculations of Eqs. (14) in the conditions found above for a probe beam in the shape of the symbol (R) with features of a ≈ 0.025 mm 2 (corresponding to k ⊥ = k 0 = 40 mm −1 ). The propagation distance is L = 45 mm, equivalent to 2 Rayleigh distances as evident by the substantial free-space diffraction. Indeed when Z (2) = 0, the FWM medium dramatically reduces the distortion of the image due to diffraction. Note that the image spectrum (black dashed-dotted line) lies barely within the resolution limit and that the Stokes distortion due to diffraction is also reduced. Direct numerical solutions of Eqs. (7) give exactly the same results. For the hot sodium system, the required control-field intensity is on the order of 100 mW for beams with a waist radius of a few mm, which is practically a plane wave on the length scale of the image.

C. Negative paraxial diffraction
Another interesting application of diffraction manipulation is imaging by negative diffraction, similar to the one proposed in Ref. [6]. Using the same tools as above, one can find the conditions for the reversal of paraxial diffraction, namely when Re Z (2) vanishes and Im Z (2) = 1 (free space diffraction is equivalent to Z (2) = −i).
At these conditions, as demonstrated in Fig. 5, the FWM medium of length L focuses the radiation from a point source at a distance u < L to a distance v behind the cell, where u + v = L. The mechanism is simple: each k ⊥ component of the probe accumulates outside the cell the phase −ik 2 ⊥ (u + v) /(2q) = −ik 2 ⊥ L/(2q) and inside the cell the the phase ik 2 ⊥ L/(2q), summing up to zero phase accumulation. The probe image thus 'revives', with some additional gain, at the exit face of the cell.

V. CONCLUSIONS AND DISCUSSION
We have suggested a scheme utilizing the k ⊥dependency of the four-wave mixing process for manipulating the diffraction during the propagation. The inherent gain of the FWM process allows one to take advantage of high optical-depths while avoiding absorption and, by that, achieving higher resolution than with previous EIT-based schemes [5,6]. As oppose to a recent proposal incorporating FWM [20], our scheme does not require atomic motion and is expected to work even more efficiently in its absence. We have introduced a microscopic model for the FWM process, based on Liouville-Maxwell equations and incorporating Doppler broaden-ing, and verified it against recent experimental results. We have delineated the conditions for which, according to the model, the FWM process suppresses the paraxial diffraction. We have also demonstrated the flexibility of the scheme to surpass the regular diffraction and reverse it, yielding an imaging effect while introducing gain. Our proposal was designed with the experimental limitations in mind, and its demonstration should be feasible in many existing setups.
The resolution limit a −1 ∝ k 2 0 of our scheme (and thus the number of 'pixels' S/a for a given beam area S) is proportional to the resonant optical depth . In practice, the latter can be increased either with higher atomic density N or narrower optical transitions. For example, using a density of N = 5 · 10 12 cm −3 , 10 times higher than in the sodium setup of Ref. [24], the limiting feature size would be 250 µm 2 (k 0 ≈ 125 mm −1 ). As long as N L =const , the other parameters required for the suppression of diffraction remain the same. At the same time, avoiding Doppler broadening by utilizing cold atoms or perhaps solids would substantially increase the resolution limit. Assuming cold atoms with practically no Doppler broadening (and ground-state relaxation rate Γ lr = 100 Hz), the same limiting feature of 250 µm 2 can be obtained at a reasonable density of 10 12 cm −3 . Finally, we note that the best conditions for suppression of diffraction are not always achieved by optimizing Z (2) alone (first order in k 2 ⊥ /k 2 0 ); In some cases, one could improve significantly by working with higher orders. As demonstrated in Fig. 6, Combining the aforementioned methods for resolution enhancement with N = 10 12 cm −3 cold atoms, a resolution-limited feature size of down to about 100 µm 2 with unity gain can be achieved. Going beyond this resolution towards the 1 − 10 µm 2 -scale, for applications in microscopy or lithography, would require further work aimed at lifting the paraxial assumption.
The FWM process conserves quantum coherence on the level of single photons, as was previously shown theoretically [28] and experimentally [29] by measuring spatial coherence (correlation) between the outgoing probe and Stokes beams. An intriguing extension of our work would thus be the generalization of the scheme to the single-photon regime. Specifically, the main limitation in the experiment of Ref. [29] was the trade-off between focusing the beams to the smallest spot possible while keeping the 'image' from diffracting throughout the medium. Our scheme could circumvent this trade-off by maintaining the fine features of the image along much larger distances.
In addition, our scheme can be utilized in optical trapping experiments for the production of traps that are long along the axial direction and tight in the transverse direction. An intricate transverse pattern can be engineered, for example, a thin cylindrical shell or a 2D array of narrow wires, that will extend over a large axial distance to allow for high optical depths. This can be further extended by modulating the control fields along the axial direction, such that the probe (and Stokes) diffracts in the absence of the controls and 'anti-diffracts' in their presence. In this arrangement, non-uniform traps in the axial direction can be designed.
Assuming control fields constant in time and space and much stronger than the probe and Stokes fields, the steady-state solution of Eqs. (3)-(5) can be approximated to lowest orders in the weak fields as R α,α ′ ≃ R α,α ′ is the zero-order and R (1) α,α ′ the first-order steady-state solutions. We find R (0) α,α ′ from the zero-order equations of motion by solving (∂/∂t)R (0) α,α ′ = 0. Under the assumption |Ω c /ω lr | ≪ 1, thus R (0) r,l = 0, we obtain for the other elements with the denominator and A l/r = |Ω c | 2 Im[(ω ul/ur − ω c − iΓ d ) −1 ] the optical pumping rates.
To find R (1) α,α ′ , we start from the first-order equations of motion, ∂ ∂t R (1) l,l = −2 Im Ω * p e −iδ k z+iδω t R eliminating the explicit dependency on time. The steadystate solution is obtained from the complete set of linear algebraic equations for the variables P The exact solution of Eqs. (A6) is easily obtained but is unmanageable and bears no physical intuition. Rather,