## Abstract

We explore theoretically the feasibility of using frequency conversion by sum- or difference-frequency generation, enabled by three-wave-mixing, for selectively multiplexing orthogonal input waveforms that overlap in time and frequency. Such a process would enable a drop device for use in a transparent optical network using temporally orthogonal waveforms to encode different channels. We model the process using coupled-mode equations appropriate for wave mixing in a uniform second-order nonlinear optical medium pumped by a strong laser pulse. We find Green functions describing the process, and employ Schmidt (singular-value) decompositions thereof to quantify its viability in functioning as a coherent waveform discriminator. We define a selectivity figure of merit in terms of the Schmidt coefficients, and use it to compare and contrast various parameter regimes via extensive numerical computations. We identify the most favorable regime (at least in the case of no pump chirp) and derive the complete analytical solution for the same. We bound the maximum achievable selectivity in this parameter space. We show that including a frequency chirp in the pump does not improve selectivity in this optimal regime. We also find an operating regime in which high-efficiency frequency conversion without temporal-shape selectivity can be achieved while preserving the shapes of a wide class of input pulses. The results are applicable to both classical and quantum frequency conversion.

© 2013 Optical Society of America

## Corrections

Dileep V. Reddy, M. G. Raymer, C. J. Mckinstrie, L. Mejling, and K. Rottwitt, "Temporal mode selectivity by frequency conversion in second-order nonlinear optical waveguides: erratum," Opt. Express**25**, 7998-7998 (2017)

https://www.osapublishing.org/oe/abstract.cfm?uri=oe-25-7-7998

## 1. Introduction

Efficient multiplexing of signals in and out of multiple optical channels is central to both quantum and classical optical communication networks. This is accomplished using an add/drop device (or filter). The channels are defined by a set of modes, and ideally, these modes are orthogonal field distributions. For example: a discrete set of frequencies (wavelength-division multiplexing, WDM), or time-bins (time-division multiplexing, TDM), or polarization (polarization-division multiplexing, PDM). True multiplexing, meaning the ability to efficiently route, add, and drop signals between channels, can be accomplished using the above-mentioned schemes. A less powerful form of multiplexing is offered by using schemes that do not permit efficient signal routing, adding, and dropping, but do allow detection of signals in different channels. Examples of this are optical code-division multiple access (OCDMA), quadrature amplitude modulation (QAM), and optical orthogonal frequency-division multiplexing (OFDM) [1, 2]. In those cases the detector has special capabilities that allow the separate detection of signals associated with different channels, but cannot efficiently separate them into distinct spatial channels for routing.

A goal is to accomplish true orthogonal-waveform multiplexing, which would allow signals in different optical channels defined by different orthogonal waveforms to be spatially separated [3]. The waveforms making up the orthogonal basis set will be overlapping in both time and frequency spectra, so standard frequency or time separation techniques do not apply. Because waveforms, or wave packets, have unique signatures in both time and frequency, we will call a system based on such a scheme orthogonal time-frequency-division multiplexing (OTFDM).

Nonlinear optical processes such as three-wave mixing (TWM) have previously been applied to quantum networking in the context of improving infrared single-photon detection efficiency by up-conversion [4–6], and WDM by down-conversion [7]. An important step toward optical OTFDM was made by the group of C. Silberhorn, who proposed an optical pulse gate based on nonlinear sum-frequency generation by TWM [8, 9]. The purpose of the present paper is to explore this TWM frequency conversion scheme theoretically, with the introduction of analytical tools that help clarify the physics. Specifically, we want to check a conjecture that the optimum operating regime for shape-selective and efficient frequency conversion is that in which one of the signals copropagates with the same group velocity as the pump pulse [9]. The results are applicable to classical and quantum frequency conversion. Related results have also been explored for the case of four-wave mixing [10–12].

## 2. Equations of motion and figure of merit

We are concerned with sum/difference frequency generation processes involving three-wave mixing in any *χ*^{(2)}-nonlinear medium. We designate the pulse-carrier frequencies of the three participating field channels as *ω _{s}*,

*ω*, and

_{r}*ω*, where

_{p}*ω*is the strong-pump channel, and

_{p}*ω*,

_{s}*ω*are the weak signal and idler channels (we assume

_{r}*ω*<

_{s}*ω*). Though we account for group velocity mismatch between the channels, we restrict our analysis to sufficiently narrow-band (broad duration) pulses so as to neglect higher order effects such as group velocity dispersion. Starting from the standard three-wave interaction equations for these channels and assuming a strong non-depleting pump pulse, an appropriate choice of channel carriers and polarizations that ensures energy conservation (

_{r}*ω*=

_{r}*ω*+

_{s}*ω*), and phase-matching (

_{p}*k*−

_{r}*k*−

_{s}*k*= 0) yields the following evolution equations in the spatio-temporal domain [13]:

_{p}*j*∈ {

*s*,

*r*,

*p}*,

*β*:=

_{j}*β*

^{(1)}(

*ω*) are the group slownesses of pulses with carrier-frequency

_{j}*ω*(in any arbitrary frame), and

_{j}*γ*is a measure of the mode-coupling strength, which is a product of the effective

*χ*

^{(2)}-nonlinearity coefficient and the pump power. We assume the pump pulse is strong enough that it remains unaltered by the interaction, but not so strong as to affect the group velocities of the signal/idler pulses. The mode-amplitudes

*A*(

_{j}*z*,

*t*) can be interpreted either as the quantum wavefunction amplitudes in the single-photon case [14], or as the pulse-envelope functions in the slow-varying envelope approximation in the classical-pulse limit:

*A*(

_{p}*t*)|

^{2}

*dt*= 1). We denote the length of our uniform-medium with

*L*, and assume the interaction starts at

*z*= 0. The solutions to Eqs. (1a) and (1b) can be represented using the Green function (GF) formalism, thus:

*A*(0,

_{k}*t′*) are the input amplitudes and

*A*(

_{j}*L*,

*t*) are the output amplitudes for

*j*,

*k*∈ {

*r*,

*s*}. The overall GF is unitary, but the block transfer functions (

*G*(...)) by themselves, are not. We can affect the GF of the process by varying the medium length (

_{jk}*L*), pump power (

*γ*), pump pulse-shape (

*A*(

_{p}*t*)) and the group-slownesses (inverse group velocities) of the various channels (

*β*). If the GF is ‘separable’, i.e.

_{j}*G*(

_{rs}*t*,

*t′*) = Ψ(

*t*)

*ϕ*

^{*}(

*t′*), then with sufficient pump power, an incident

*s*-channel signal of temporal shape

*ϕ*(

*t′*) can be 100% converted into the outgoing

*r*-channel packet Ψ(

*t*), and any incoming signal that is orthogonal to

*ϕ*(

*t′*) will be left unconverted. In general however, the GF is not separable. The ability to separate temporal modes becomes easier to quantify if we represent the GF with its singular-value decomposition [10, 11, 15, 16]:

*ψ*(

_{n}*t′*),

*ϕ*(

_{n}*t′*) are the input “Schmidt modes” and Ψ

*(*

_{n}*t*), Φ

*(*

_{n}*t*) are the corresponding output Schmidt modes for the

*r*and

*s*channels respectively. These functions are uniquely determined by the GF, and form orthonormal bases in their relevant channels. The “transmission” and “conversion” Schmidt-coefficients (singular values) {

*ρ*} and {

_{n}*τ*} are constrained by |

_{n}*τ*|

_{n}^{2}+ |

*ρ*|

_{n}^{2}= 1 to preserve unitarity. It is convenient to choose the mode-index ‘n’ in decreasing order of Schmidt mode conversion-efficiency (CE) (|

*ρ*|

_{n}^{2}). The process is deemed perfectly mode-selective if the CE |

*ρ*

_{1}|

^{2}= 1, and |

*ρ*| = 0; for every

_{m}*m*≠ 1 (i.e., the TWM process performs full frequency conversion on one particular input mode and transmits all power from any orthogonal mode in the same input channel). We can quantify the add/drop quality of the GF using the ordered-set of conversion efficiencies to define an add/drop ‘selectivity’:

*ρ*

_{1}|

^{2}) is the CE of the dominant temporal mode. The selectivity characterizes both the degree of separability of the GF and the process efficiency. The equality in Eq. (5) holds for a perfect add/drop device. The unitary nature of the transformation imposes a pairing between the Schmidt modes across the

*r*and

*s*channels [10, 17]. Consider arbitrary input and output fields expressed as discrete sums over corresponding Schmidt modes:

*a*,

_{n}*b*,

_{n}*c*,

_{n}*d*} are pairwise related via a unitary beam-splitter-like transformation [10], which, if we assume real

_{n}*τ*and

_{n}*ρ*, are expressed as: where the

_{n}*n*

^{th}-Schmidt mode CE (|

*ρ*|

_{n}^{2}= 1 − |

*τ*|

_{n}^{2}) is analogous to “reflectance”. All time-domain functions described thus far have corresponding frequency-domain analogs. The form of the GF in frequency domain can also provide meaningful insights. If we define functions Ψ̃

*(*

_{n}*ω*) and

*ϕ̃*(

_{n}*ω*) as the Fourier-transforms of the corresponding time-domain Schmidt modes Ψ

*(*

_{n}*t*) and

*ϕ*(

_{n}*t*), then:

## 3. Low-conversion limit

We can develop an important guide to the different regimes of TWM by solving the problem for small interaction strengths (*γ*) for arbitrary group slownesses and pulse shapes (following the discussion in [12]). We define the coupling coefficient as *κ*(*z*, *t*) = *γA _{p}*(

*t*−

*β*). For this calculation, we could allow the nonlinearity

_{p}z*γ*(

*z*) to be position dependent, which can be used as a design feature if desired [8, 20], but for simplicity we continue to assume that the medium is uniform. By integrating Eqs. (1a) and (1b) with respect to

*z*, we get the exact relations:

*t′*:=

_{r}*t*−

*β*(

_{r}*L*−

*z′*) and

*t′*:=

_{s}*t*−

*β*(

_{s}*L*−

*z′*). Treating the coupling as a perturbation, we get

*t*=

_{r}*t*−

*β*,

_{r}L*t*=

_{s}*t*−

*β*. The ≈ symbols indicate that perturbative approximations render Eqs. (10a) and (10b) weakly non-unitary. By defining

_{s}L*t′*=

*t*−

*β*+

_{r}L*β*, where

_{rs}z′*β*=

_{rs}*β*−

_{r}*β*is the difference in slownesses, one can change the integration variable to time, and rewrite Eqs. (10a) and (10b) using the approximate Green function

_{s}*Ḡ*(

_{jk}*t*,

*t′*) in the low-conversion limit [12]:

*t*=

_{j}*t*−

*β*,

_{j}L*β*=

_{jp}*β*−

_{j}*β*; ∀

_{p}*j*∈ {

*r*,

*s*}, and

*H*(

*x*) is the Heaviside step-function.

Eqs. (12a) and (12b) provide a simple way to understand the FC process for arbitrary relations between group slownesses. A first observation is that the Heaviside step-functions represent that because the medium length is finite, the time of interaction is restricted to the *t′* interval *t′* ∈ (*t* − *β _{r}L*,

*t*−

*β*). In the (

_{s}L*t*,

*t′*) domain, this interval corresponds to a 45°-tilted band with width

*β*. Therefore, if the goal is to have the GF separable in

_{rs}L*t*and

*t′*, the shape of this interval poses a challenge. Whatever the pulse shape of the pump, the low-conversion Green function is proportional to a scaled version of that shape. Note that if (

*β*= 0 or

_{rp}*β*= 0), then the factor

_{sp}*A*(...) in Eqs. (12a) and (12b) depends only on

_{p}*t*(or

*t′*), making that factor in the GF separable. Further insight is obtained by plotting the GF, as in Fig. 1, for the case of a (normalized) Gaussian pump pulse ${A}_{p}\left(t\right)={\left({\tau}_{p}^{2}\pi \right)}^{-1/4}\text{exp}\left[(-{t}^{2}/\left(2{\tau}_{p}^{2}\right)\right]$ with duration

*τ*.

_{p}For the four GF’s plotted in Fig. 1, the computed CE’s for the first four temporal modes are listed in Table 1, where *γ̄* = *γ*/*β _{rs}* is of order 0.01. The corresponding selectivities are:

*S*=

**(a)**0.646

*γ̄*

^{2},

**(b)**0.676

*γ̄*

^{2},

**(c)**0.646

*γ̄*

^{2},

**(d)**0.610

*γ̄*

^{2}. The Schmidt coefficients were numerically computed by performing a singular value decomposition (SVD) of the GF in Eq. (12a).

The slope of the line along the highest part of the band (lightest color) is given by

In an attempt to create an approximately separable GF, one can choose parameters as in Fig. 2. The computed CE’s for the first four temporal modes in this case are {1.0*γ̄*^{2}, 0.029 *γ̄*^{2}, 0.028*γ̄*^{2}, 0.011*γ̄*^{2}}, where *γ̄* = *γ*/*β _{rs}* is of order 0.01. The selectivity is

*S*= 0.913

*γ̄*

^{2}. While the separability is high, the CE of the first Schmidt mode is of the order of

*γ̄*

^{2}. Improved selectivity can be achieved using the strategy proposed in [9], where one of the signals is matched in slowness to the pump, as in Fig. 1(a), and the pump pulse is made very short. The short pump width helps counter the ill effects of the 45°-sloping step-functions on GF separability by selecting a narrow vertical or horizontal region in (

*t*,

*t′*) space. These choices give the GF’s (in the low-CE limit) in Fig. 3. The numerically computed CE’s for the first four temporal modes in Fig. 3(a) are {1.0

*γ̄*

^{2}, 0.022

*γ̄*

^{2}, 0.006

*γ̄*

^{2}, 0.003

*γ̄*

^{2}}, and the selectivity is

*S*= 0.967

*γ̄*

^{2}. In Fig. 3(b) the CE’s and the selectivities are identical to case 3(a).

*γ̄*is of order 0.01.

The temporal Schmidt modes for the case in Fig. 3(b) are shown in Fig. 4. It is seen that the input modes mimic the projection of the GF onto the *t′*-axis, while the output modes mimic the projection onto the *t*-axis. If the input field occupies only the dominant (*j* = 1) Gaussian-like mode, then it is frequency converted with efficiency |*ρ*_{1}|^{2} ≈ *γ̄*^{2} and generates an output pulse that is much longer and rectangular in shape. Such pulse shaping may or may not be desired, depending on the application.

The *s*-output modes for case 3(a) and the *r*-output modes for case 3(b) will have temporal width *β _{rs}L*, which is the maximum duration of interaction within the medium. Since the pump copropagates with a matched slowness with one of the input channels, and the CE is low enough to prevent input-channel depletion, FC occurs throughout the traversed medium length, stretching the generated output mode in the other channel due to difference in slownesses (

*β*).

_{rs}To demonstrate the ability to choose which temporal mode is selected for FC, Fig. 5 shows the results for a pump pulse with the shape proportional to a first-order Hermite-Gaussian function *HG*_{1}(*x*) ∝ *x*exp[−*x*^{2}/2], which has a zero-crossing at its “midpoint”. The efficiencies are {1.0*γ̄*^{2}, 0.049*γ̄*^{2}, 0.007*γ̄*^{2}, 0.005*γ̄*^{2}} and the selectivity is *S* = 0.936*γ̄*^{2}, where *γ̄* = *γ*/*β _{rs}* is of order 0.01. The dominant mode has a shape similar to the pump pulse.

Alternatively, the preceding analysis can be carried out in the frequency-domain, where the GF takes the form:

*ω̄′*= (

*β*−

_{rp}ω*β*)/(2

_{sp}ω′*β*). Varying the pump duration

_{rs}*τ*changes the bandwidth of factor

_{p}*g*

_{1}(

*ω*−

*ω̄′*), and the choice of slowness (

*β*) and medium length (

_{j}*L*) affects the slope and the phase-matching bandwidth of factor

*g*

_{2}(

*ω̄′*) in (

*ω*,

*ω′*) space. The separability of

*Ḡ*(

_{rs}*t*,

*t′*) is also evident in (

*ω*,

*ω′*)-space. As pointed out in [10], for the case in Fig. 3(a) with

*β*= 0,

_{rp}*g*

_{2}(

*ω̄′*) would be a sinc-function parallel to the

*ω*-axis with a phase-matching bandwidth proportional to 1/(

*β*) [a measure of the vertical separation between the edges of the heaviside-step functions in Fig. 3(a)], and the shortness of the pump will cause

_{sp}L*g*

_{1}(

*ω*−

*ω′*) to have a wide-bandwidth, intersecting

*g*

_{2}(

*ω̄′*) at a 45° inclination. Alternatively, one could choose parameters such that

*β*= −

_{sp}*β*, giving

_{rp}*g*

_{2}(

*ω̄′*) a −45° inclination. In the frequency domain, if the pump bandwidth and medium length are optimized, we can reproduce a roughly separable GF as in Fig. 2, with (

*t*,

*t′*) replaced by (

*ω*,

*ω′*). GF separability suffers in this regime if the pump bandwidth is made much larger than the phase-matching bandwidth.

To summarize, in the low CE limit, we are able to achieve high temporal mode separability when the GF is nearly separable (Fig. 3) and moderate separability when the GF has a modicum of symmetry (Fig. 2). But the low CE diminishes the selectivity. We next address whether high selectivity can also be found in cases with higher conversion efficiencies.

## 4. High-conversion regimes

Assuming energy conservation and perfect phase-matching for the channel carrier frequencies, the choice of waveguide/material dispersion is reflected in our equations via the relative magnitudes of the channel group slownesses. We classify the different regimes of operation as follows:

**Single sideband velocity matched**SSVM:

*β*=_{s}*β*≠_{p}*β*or_{r}*β*=_{r}*β*≠_{p}*β*_{s}**Symmetrically counter-propagating**SCuP:

*β*= −_{rp}*β*_{sp}**Counter-propagating signals**CuP:

*β*< 0,_{sp}β_{rp}*β*≠_{s}*β*_{r}**Co-propagating signals**CoP:

*β*> 0,_{sp}β_{rp}*β*≠_{s}*β*_{r}**Exactly co-propagating**ECoP:

*β*=_{r}*β*_{s}

In this section, we employ numerical techniques similar to those used in [11] to construct the GF for any given set of pump parameters. To accomplish this we numerically propagate a large number of ‘test signals’ through the medium (chosen to be members of a complete, orthonormal set of Hermite-Gaussian functions of appropriate temporal width) to find the effects of the process on an arbitrary input. This method (described in Appendix I) enables a comprehensive study of TWM, even for cases for which analytical solutions are not known.

We first present numerical results for the SSVM regime, which has been favored by C. Silberhorn’s group [8], and has yielded the best results in terms of selectivity.

#### 4.1. (*β*_{sp} = 0, *β*_{rs} ≠ 0*) Single sideband velocity matched regime*

_{sp}

_{rs}

The function of an effective add/drop device is to efficiently discriminate between orthogonal temporal modes. Since any channel input enters and traverses through the waveguide in causal sequence (linearly with the pulse function argument), to achieve discrimination the pump pulse must overlap with all segments of the input pulse for a non-zero amount of time within the waveguide. This ensures that: (a) all the power distributed among all the segments of the first input Schmidt mode has a chance of interacting with the pump and being FC’d into the other channel, and (b) the device “measures” the entire shape of the temporal input mode, which is essential for discriminating between different temporal mode shapes.

Two orthogonal temporal modes (say in channel *r*) can have locally similar shapes in certain segments. When these segments overlap with the pump pulse within the nonlinear medium, the only way for the device to react to them differently is for the local instantaneous mode features in channel *s* to differ (Eq. (1a)), which is determined by all the TWM that has occurred until that time instant. Both of these intuitive requirements are satisfied if one of the channel slownesses is matched to the pump slowness (single sideband group-velocity matched or SSVM), and the temporal pump width is much shorter than the interaction time *β _{rs}L*. The preceding low CE-limit analysis has already deemed this regime favorable for separability, and other groups [9] have predicted significant success at higher CE’s as well. We now present the numerical results for the same. We present the complete exact-analytical solution for the SSVM case in section 5. In the SSVM regime, for a given pump shape, the selectivity is influenced most by the GF aspect ratio (

*τ*/(

_{p}*β*)) and effective interaction strength (

_{rs}L*γ̄*=

*γ*/

*β*).

_{rs}In Fig. 6, we plot the numerically determined CE for the first five Schmidt modes for various *γ̄* for a Gaussian pump-pulse with parameters from Fig. 3(b) (*β _{r}* = 2,

*β*=

_{s}*β*= 0,

_{p}*τ*= 0.1,

_{p}*L*= 1). The selectivity values are listed in the inset. A maximum selectivity of 0.81 is found for

*γ̄*= 1.0.

Although the GF displays good mode-separability at low-CE’s, the selectivity is unable to maintain high values beyond a certain *γ̄*. Figure 7 shows the first three input and output Schmidt modes for *G _{rs}*(

*t*,

*t′*) for the same case. Figure 8 shows the first Schmidt modes from both channels for increasing

*γ̄*. Note the strong shape distortion relative to the low-CE case, reflecting the change in the GF shape with increasing

*γ̄*. This illustrates the limits of validity of the approximation used in [9].

Shortening the pump width by a factor of 10 minutely improves the selectivity, while lengthening pump width causes it to decrease. We present the analytical solution for this SSVM regime in section 5, where we show that this case leads to the highest selectivity of all the regimes treated in this study.

#### 4.2. (*β*_{rp} = −*β*_{sp}) *Symmetrically counter-propagating signals regime, shape preserving FC*

_{rp}

_{sp}

We now treat the SCuP regime, in which the signals propagate in opposite directions in the pump frame with the same slowness relative to the pump pulse. Specifically, for this section we work with parameter values: *β _{s}* = 0,

*β*= 2,

_{p}*β*= 4,

_{r}*L*= 1, and Gaussian-shaped pump. For pump width

*τ*= 0.707, the low-CE Green function matches a time-shifted version of the plot in Fig. 2. Increasing

_{p}*γ̄*to higher-CE will cause the selectivity

*S*to rise to a maximum, and then fall back to lower values, just like in the SSVM regime, Fig. 9 plots selectivity vs.

*γ̄*for various pump widths.

We sought to improve this result by varying the low-CE GF aspect ratio (*τ _{p}*/(

*β*)) by changing

_{rs}L*τ*. Starting from narrow pumps and increasing the width, we could optimize the selectivity-maxima to about

_{p}*S*≈ 0.7 for a Gaussian pump width of

*τ*≈ 1.5 and

_{p}*γ̄*≈ 0.75. Further increasing

*τ*stretched the GF shape in the

_{p}*t*=

*t′*direction, reducing its separability/selectivity-maximum, as shown in Fig. 9. The selectivity maximum moves to larger

*γ̄*for increasing

*τ*because longer pump durations correspond to smaller peak intensity.

_{p}Figure 10 shows how the CEs for the first ten Schmidt modes change with *γ̄* for various *τ _{p}*. For large

*γ̄*, higher-order CEs tend to decrease with increasing

*τ*, suggesting mildly improved selectivity. They also appear to oscillate about a decreasing central value in a damped fashion with increasing

_{p}*τ*. For the values plotted, this is most pronounced in the CE of the third Schmidt mode for

_{p}*γ̄*= 1.18.

In this SCuP regime we find the shapes of the output (r) Schmidt modes are essentially identical to those of the input (s) Schmidt modes. Figure 11 shows the dominant *s* input and *r* output Schmidt modes at *γ̄* = 3.36, for select *τ _{p}*. This shape-preserving behavior is related to the GF consisting of the pump shape as a factor sloping parallel to the

*t*= −

*t′*direction, and is independent of

*γ̄*for the values tested. The individual Schmidt mode shapes, however, do change with

*γ̄*. For

*τ*≫ 0.1, the Schmidt mode widths scale linearly with the pump width.

_{p}The time-widths of the Schmidt modes have a lower bound of *β _{rs}L*/2 due to the Heaviside step-function boundaries. Decreasing

*τ*to small values relative to

_{p}*β*(e.g., 0.1) causes the convergence of Schmidt mode shapes to those plotted in Fig. 11(a,b). The dominant CE’s for the short-pump case nearly match each other in values, especially for very low and very high

_{rs}L*γ̄*, making for a non-selective add/drop device. These features cause the short-pump SCuP regime to preserve the shapes of a large family of input pulses during FC (even for CE’s approaching unity). For example, for the

*τ*= 0.1,

_{p}*γ̄*= 3.36 case in Fig. 10, all the first seven Schmidt modes have near unity CE’s. So any input pulse that can be completely constructed by a linear superposition of the first seven input Schmidt modes will FC into the other channel into the exact same superposition of the first seven output Schmidt modes, which also match the corresponding input Schmidt mode shapes (Fig. 11). We hypothesize that this results from the

*t′*= −

*t*direction of the GF being more pronounced for shorter

*τ*(Fig. 12), which maps local time-slices/segments of the input and output pulses in a one-to-one fashion.

_{p}Note that the pump in the GF is constant along the −45° direction. So the interaction of each time segment of the signal with the corresponding segment in the idler is driven by the same pump profile. The resolution of such one-to-one segment mappings is determined by the pump width. For broad pumps, any given time segment of the signal would then influence a larger portion of the idler pulse, and vice versa. The inability of the global shape of an input pulse to influence its CE results in a poor add/drop device, but this feature, which we call shape-preserving FC, has potential applications in multi-color quantum interference [10, 21].

#### 4.3. (β_{rp}β_{sp} < 0) counter-propagating signals regime

The previous section dealt with the parameter set *β _{s}* = 0,

*β*= 4,

_{r}*β*= 2. Holding the

_{p}*β*and

_{s}*β*slownesses at these values, we now vary the pump slowness (

_{r}*β*) within the range [0, 4] and chart the properties of the GF. At values 0 and 4, the results matched those of the SSVM regime. The range

_{p}*β*∈ [0, 2] showed a one-to-one symmetrically-mapped correspondence with the range [2, 4]. That is, for every Δ in the range [0, 2], the GF had the same selectivities for

_{p}*β*= (2 − Δ) as well as (2+Δ). Even the Schmidt modes were identical but interchanged between the signal channels.

_{p}As the low-CE GF plots in Fig. 1 show, the pump-shape factor in *G _{rs}*(

*t*,

*t′*) has slope

*β*/

_{sp}*β*, defined in Eq. (13). For fixed

_{rp}*L*= 1, changing this slope, particularly for small pump widths, will change the projected width of the GF on the

*t*and

*t′*axes, which changes the widths of the Schmidt modes. This is also true for arbitrary CE’s, as is shown in Fig. 13. Bringing

*β*closer to

_{p}*β*will tend to align

_{r}*G*(

_{rs}*t*,

*t′*) with the vertical

*t′*-direction. This increases the

*s*-channel Schmidt mode widths, and decreases the

*r*-channel Schmidt mode widths.

Figure 14 shows the plots of selectivity vs. *γ̄* for different *β _{p}* values and pump widths. While in the SCuP regime, the selectivity-maximum was highest for

*τ*≈ 1.5. As

_{p}*β*drew closer to

_{p}*β*= 4, the optimum-selectivity-pump width was seen to decrease. This is consistent with our finding for the SSVM regime, which shows larger selectivity-maxima for shorter pumps. The selectivity-maximum also increased as we approached the SSVM regime.

_{r}The selectivities for short-pumps were hyper-sensitive to changes in *β _{p}* since the shape of the GF is affected the most (due to pump-factor slope defined in Eq. (13)) for shorter pumps. This implies that the closer we are to SSVM regimes (but not in it), the shorter our pump needs to be for the FC to still be shape-preserving. Selectitivies for wider-pumps did not show the same sensitivity to changes in

*β*.

_{p}#### 4.4. (β_{rp}β_{sp} > 0) co-propagating signals regime

In this section, we explore the regime in which the slope of the pump factor in the low-CE GF, i.e. the quantity *β _{sp}*/

*β*is positive. We do this by fixing

_{rp}*β*= 4,

_{p}*β*= 0,

_{s}*L*= 1, and varying

*β*within the range [0, 4]. Selectivity behavior for negative values of

_{r}*β*mapped bijectively to the corresponding positive

_{r}*β*that resulted in an inversion in pump-factor slope, while the Schmidt modes swapped across the

_{r}*r*and

*s*channels.

Figure 15 consists of selectivity vs. *γ̄* plots for various pump widths and *β _{r}*. The selecivity maximum for any given

*τ*, apart from decreasing in magnitude with decreasing

_{p}*β*, also migrates to higher

_{r}*γ̄*values. This effect is more pronounced for shorter pumps. The optimum pump width (with the highest selectivity maximum) also increases with decreasing

*β*.

_{r}As *β _{r}* →

*β*, the pump-factor slope approaches unity. This allows for shape-preserving FC behavior when using short pumps, through a mechanism analogous to that illustrated in Fig. 12, except here the idler pulse convects through the pump in the same direction as the signal pulse. CE’s for the first ten Schmidt modes for small

_{s}*β*tended to match each other, confirming non-shape-descriminatory GF. This “rotation” of the GF pump-factor causes the Schmidt mode widths to track the GF projection on the (

_{r}*t*,

*t′*)-axes. The difference is most noticeable for short pumps (Fig. 16). The dominant Schmidt modes in both channels converge to matching shapes, as expected for shape-preserving FC.

The pump-factor slope can also be made to approach unity by keeping *β _{r}* and

*β*fixed and increasing

_{s}*β*to very high magnitudes. This approach would maintain the spacing between the Heaviside step-functions and prevent the selectivity maximum from migrating to higher

_{p}*γ̄*values. Numerical constraints restrain us from covering the entire range of the pump-factor slope using this method.

#### 4.5. (β_{r} = β_{s}) Exactly co-propagating signals regime

The ECoP regime is special in that we cannot plot the low-CE GF as we did for all the other regimes. As *β _{rs}* → 0, the separation between the Heaviside step-functions also converges to zero. We can however, explicitly write down the complete analytical solution for real pump-functions. If

*β*=

_{s}*β*= 0 and

_{r}*A*(

_{p}*x*) ∈ ℝ, then:

_{βp→0}

*P*(

*z*) =

*γA*(

_{p}*t*)

*z*. The GF are

*δ*-functions in

*t′*, and do not lend themselves to numerical Schmidt decomposition. This regime is beyond the scope of our simulation methodology (detailed in Appendix I).

The absence of walk-off between the two signal channels implies that the evolution of *A _{r}*(

*z*,

*t*) for a given local time index ‘

*t*’ is insensitive to the global shapes of the input wavepackets (

*A*(0,

_{r}*t′*),

*A*(0,

_{s}*t′*)). This results in a poor add/drop-device. Different time slices of arbitrary input pulses will undergo the same transformation as they sweep across the pump, allowing for distortionless conversion.

The SSVM regime (*β _{s}* =

*β*), under the

_{p}*β*→ 0 limit also converges to Eqs. (16a) and (16b) (as verified in Appendix II). The exact solution in this regime for arbitrary complex-valued (chirped) pump functions will not be dealt with in this publication.

_{rs}## 5. Analytical solution for single sideband velocity matched regime

The SSVM regime, where *β _{s}* =

*β*, and all other parameters are arbitrary, was shown above to be the optimal regime for the drop/add process. Fortunately, in this same regime the problem can be solved analytically, following [22](detailed in Appendix III). The exact GF is found to be

_{p}*τ*=

*t*−

*β*,

_{s}L*τ′*=

*t′*,

*ζ*=

*β*−

_{r}L*t*,

*ζ′*= −

*t′*,

*ξ*=

*ζ*−

*ζ′*,

*γ̄*=

*γ*/

*β*, and $\eta ={\int}_{{\tau}^{\prime}}^{\tau}{\left|{A}_{p}\left(x\right)\right|}^{2}dx$.

_{rs}*J*(...) is the Bessel function of order

_{n}*n*, and

*H*(

_{H}*τ*,

*τ′*,

*ζ*,

*ζ′*) =

*H*(

*τ*−

*τ′*)

*H*(

*ζ*−

*ζ′*),

*H*(

*x*) being the Heaviside step-function.

For an analysis of selectivity/separability, we need only consider the structure of *G _{rs}*(

*t*,

*t′*), which has two non-separable factors in (

*t*,

*t′*): the Bessel function ${J}_{0}\left\{2\overline{\gamma}\sqrt{\eta \xi}\right\}$, and the step-functions

*H*(

_{H}*τ*,

*τ′*,

*ζ*,

*ζ′*). Decreasing the pump width relative to the effective interaction time (

*β*) can diminish the ill-effects of the step-functions on GF-separability, but the effect of the Bessel function worsens at higher

_{rs}L*γ̄*. A numerical singular value decomposition of this analytical GF in Eq. (17c) for high

*γ̄*plotted in Fig. 17 confirms our numerical results from section 4.1.

Increasing *γ̄* improves the CE of the first Schmidt mode by scaling the peak of the GF, but via the Bessel function, decreases the separability (Fig. 18). Hence, selectivity, being a product of the two, attains a maximum value at around *γ̄* ≈ 1.15. While decreasing pump width (*τ _{p}*) improves selectivity, the maximum asymptotically approaches a limiting value of approximately 0.85 (Fig. 19). This Bessel function induced distortion in GF shape is reflected in the shape of the Schmidt modes (section 4.1).

One might suspect that including a frequency chirp in the pump field could improve the selectivity. We prove here that for the SSVM regime this is not the case. Note that the pump-squared integral *η*(*t*, *t′*), and consequently the Bessel function, is independent of any pump-chirp. To demostrate this, we rewrite Eqs. (1a) and (1b) in the pump’s moving frame:

Replacing the pump envelope-function by its real-amplitude and phase (*A _{p}*(

*t*) :=

*P*(

*t*)exp[

*iθ*(

*t*)]) and setting

*β*= 0 for the SSVM constraint, we get:

_{sp}*s*-channel envelope function as

*Ā*(

_{s}*z*,

*t*) =

*A*(

_{s}*z*,

*t*) exp[

*iθ*(

*t*)], we can recover Eqs. (1a) and (1b) with a real-pump envelope in the SSVM regime. Any time dependent complex phase in the pump gets absorbed into the Schmidt modes, without affecting the CE’s or GF selectivity. Nevertheless, for any given

*γ̄*, the shape of the pump gives us some control over the shapes of the Schmidt modes, and this may be used to tune the add/drop device to accept easy-to-produce pulse shapes as input Schmidt modes.

The parameter (*β _{rs}L*) is responsible for the Schmidt mode width for the channel with velocity mismatched with that of the pump. This parameter has units of time, and is a measure of the duration of “interaction”. Increasing

*β*in the SSVM regime will make higher CE’s attainable at higher pump powers but lower

_{rs}L*γ̄*. The selectivity maximum also follows a similar trend until

*β*becomes comparable to pump width

_{rs}L*τ*(at which point the slope of the Heaviside step-functions reduces overall GF separability), as shown in Fig. 20.

_{p}The non-separability arising from the Bessel function in the GF can be traced to the oscillations shown in Fig. 18. These are similar to those in Burnham-Chiao ringing [23] seen in fluorescence induced by short-pulse excitations by the propagation of short, weak pulses through a resonant atomic medium. To model an analogy, the phase- and energy-matched wave-mixing process may be represented by a 2-level pseudo-atomic-medium with a ground-state energy at *ω _{p}* and an excited state at

*ω*. Any finite-width input pulse in

_{s}*r*channel with energy resonant with the atomic-medium (

*ω*=

_{r}*ω*−

_{s}*ω*) will have a non-zero bandwidth in the frequency domain. As it interacts with the medium, its spectral-components detuned above resonance will acquire a different phase shift than the spectral components below resonance. These two spectral components will beat to produce the ringing effect, resulting in the oscillations seen in Fig. 18. The Bessel function factor is a fundamental barrier that restricts selectivity in the SSVM regime (which is thought to be the optimal one). Some groups [8, 20] have sought to make the waveguide properties non-uniform (

_{p}*γ*→

*γ*(

*z*)) in an attempt to overcome this limitation, with limited success. A full analysis of that regime is beyond the scope of this publication.

## 6. Summary and concluding remarks

In this paper, we modelled sum/difference frequency generation processes involving three-wave mixing with a strong pump in a uniform but finite *χ*^{(2)}-nonlinear medium. We used the notion of Green function separability for such processes and used it to study the feasibility of using these processes for implementing an orthogonal time-frequency-division multiplexer (OTFDM). Such a device would selectively discriminate between orthogonal weak classical pulses or single-photon wavepackets that overlap both in temporal and frequency domains. We employed singular-value-decomposition (or Schmidt decomposition) of said Green functions to define selectivity: a figure of merit that quantifies the process’s viability for application as an OTFDM add/drop device. The decomposition also produced parameter-dependent input and output Schmidt modes, which functioned as a natural orthogonal basis set of channel waveforms for the device to multiplex.

Under perfect phase-matching conditions, and ignoring higher-order dispersion (valid for long-enough pulses) as well as self- and cross-phase modulation effects, we identified the group slownesses of the participating optical channels, the medium length, the pump pulse-shape and pump power as all the parameters that determine the Green function. We were able to perturbatively approximate the Green functions for low pump powers and contrast the separabilities associated with various configurations and parameter regimes, subject to the above mentioned constraints. We then undertook an exhaustive numerical computation of selectivities [defined in Eq. (5)] in all regimes for real chirp-free Gaussian pump-pulses of arbitrary power and width within the above mentioned constraints.

We found that the best selectivity (∼ 0.83) is obtained when the group slowness of one of the signal channels is matched with that of the pump, a regime first proposed in [9], and discussed in sections 4.1 and 5. This SSVM regime resulted in Schmidt mode widths that equaled the pump width for the group-slowness-matched signal channel, and equaled the effective inter-pulse-interaction time (*β _{rs}L*) for the other channel. We then presented the complete analytical solution for this regime, which sheds light on certain parameter-scale invariances, whilst imposing a strict upper bound on the selectivity and proving its independence of pump-chirping.

We also found that symmetrically mismatching the group slownesses of the signal channels, making one of them slower than the pump and the other one faster, also yielded reasonable selectivity (up to 0.7) for moderate pump widths. This regime also had Schmidt modes of equal time widths in both channels, and for short pumps resulted in pulse shape-preserving frequency conversion with high efficiency.

We conjecture that further boosts to selectivity can be acheived only by the introduction of non-uniformity in the nonlinear medium, introduction of higher-order dispersion, or reliance on higher order processes such as four-wave-mixing.

**Note added in proof:** A very recent related publication proposed that approximately single-mode QFC can be obtained independently of the phase-matching regime by choosing a pump pulse duration comparable to and/or shorter than the inverse of the phase-matching bandwidth [24]. In the present study we did not find evidence that high selectivity, as defined by our figure of merit Eq. (5), can be obtained generally in this manner. We find that the SSVM regime, proposed by Eckstein [9], results in significantly higher selectivity than all other regimes.

The work of DR, MR and CM was supported by the National Science Foundation through ECCS and GOALI. We thank Dr. Craig Rasmussen for his help with parallelization of the simulation code. All numerical computations presented in this paper were performed on the ACISS cluster at the University of Oregon.

## Appendix I Procedure for numerical derivation of Green function

In order to derive the numerical GF for TWM, we first implemented a coupled mode equation solver that accepts arbitrary input functions (*A _{r}*(0,

*t′*),

*A*(0,

_{s}*t′*)) as arguments, and computes the resultant output functions (

*A*(

_{r}*L*,

*t*),

*A*(

_{s}*L*,

*t*)) for Eqs. (1a) and (1b). This is achieved using a Runge-Kutta based method. The solver iterates over differential steps in pulse-propagation (Δ

*z*) from

*z*= 0 to

*L*(medium length). Every iteration consists of an upwinded

*z*-propagation scheme for all three pulses (signal, pump and idler) by a step (Δ

*z*), followed by a Runge-Kutta implementation of the coupled nonlinear interaction, all in spacetime domain. Next, we compute the GF by computing the outputs for an orthogonal set of input ‘test signals’. To ellaborate, consider the GF submatrix and its singular-value-decomposition (SVD):

*ρ*, Ψ

_{j}*(*

_{j}*t*),

*ϕ*(

_{n}*t′*)) on the right-hand-side. We first pick a spanning-set of basis functions {

*B*} and {

_{r,k}*B*} and re-express the Schmidt modes:

_{s,l}*A*(0,

_{r}*t*) = 0,

*A*(0,

_{s}*t′*) =

*B*(

_{s,l}*t′*)) as inputs for the solver, and decomposing the resulting

*r*-channel outputs

*A*(

_{r}*L*,

*t*) in the {

*B*} basis will yield the entire

_{r,k}*l*-column of the complex matrix

^{th}*Ḡ*. Once this matrix is determined, its SVD will directly reveal {

_{rs}*U*}, {

_{jk}*V*}, and {

_{jl}*ρ*}, and through them, the Schmidt modes. For the results presented in this publication, we chose Hermite-Gaussian functions for the spanning-set of basis functions for both input and output Schmidt modes, since the low-CE Schmidt modes for Gaussian pump shapes are nearly Hermite-Gaussian.

_{j}## Appendix II SSVM solution converging to ECoP in the *β*_{rs} → 0 limit

_{rs}

Taking the *β _{rs}* → 0 limit whilst still enforcing SSVM conditions will cause

*β*→ 0 as well. We now verify that the exact SSVM analytical solution consistently reduces to the expected sinusoidal form expressed in section 4.5. Consider

_{p}*G*(

_{rs}*t*,

*t′*) from Eq. (17c) for the input condition

*A*(0,

_{r}*t′*) = 0:

*t*″ is being integrated from 0 to

*β*,

_{rs}L*β*→ 0 ⇒

_{rs}*t*″ → 0. Then the integral inside the bessel argument reduces to (

*β*−

_{rs}L*t*″)|

*A*(

_{p}*t*−

*β*)|

_{r}L^{2}.

We then use:

*β*= 0 and

_{r}*A*(0,

_{r}*t*) = 0.

## Appendix III Analytical derivation of the TWM Green function in the single sideband group-velocity matched regime

Consider the equations of motion (Eqs. (1a) and (1b)) for frequency conversion (FC) by three-wave mixing (TWM):

*β*=

_{s}*β*, it is convenient to define the retarded time variable

_{p}*τ*=

*t*−

*β*and the normalized distance variable

_{s}z*ζ*=

*β*−

_{r}z*t*. By using these variables, one can rewrite the equations of motion in the simplified forms

*γ̄*=

*γ*/

*β*and the differential slowness (walk-off parameter)

_{rs}*β*=

_{rs}*β*−

_{r}*β*> 0. The equations are to be solved for −∞ <

_{s}*τ*< ∞ and 0 ≤

*ζ*< ∞. An easy way to do this is by Laplace transformation in space (

*ζ*→

*s*) [25, 26]. If an impulse is applied to the signal at the input boundary (

*z*= 0), the transformed equations are

*τ′*is the source time. If

*H*(

*τ*−

*τ′*) is the Heaviside step function and the effective time variable $\eta \left(\tau ,{\tau}^{\prime}\right)={\int}_{{\tau}^{\prime}}^{\tau}{\left|{A}_{p}\left(x\right)\right|}^{2}dx$, the solutions are:

*J*is a Bessel function of order

_{n}*n*,

*ξ*=

*ζ*−

*ζ′*,and

*H*(

_{H}*τ*,

*τ′*,

*ζ*,

*ζ′*) =

*H*(

*τ*−

*τ′*)

*H*(

*ζ*−

*ζ′*).

In the text, these solutions are referred to as the Green functions *G _{rs}* and

*G*, respectively.

_{ss}If an impulse is applied to the idler at the input boundary, the transformed equations of motion and solutions are

These solutions are referred to as the Green functions *G _{rr}* and

*G*, respectively. These Green functions are stable analogs of the Green functions for stimulated Brillouin scattering (SBS) and stimulated Raman scattering (SRS) [25, 26], and are equivalent to the Green functions for anti-Stokes SRS [28].

_{sr}## References and links

**1. **T. H. Lotz, W. Sauer-Greff, and R. Urbansky, “Spectral efficient coding schemes in optical communications,” Int. J. Optoelectronic Eng. **2**, 18–25 (2012).

**2. **W. Shieh, X. Yi, Y. Ma, and Q. Yang, “Coherent optical OFDM: Has its time come? [Invited],” J. Opt. Netw. **7**, 234–255 (2008) [CrossRef] .

**3. **M. G. Raymer and K. Srinivasan, “Manipulating the color and shape of single photons,” Phys. Today **65**, 32–37 (2012) [CrossRef] .

**4. **A. P. Vandevender and P. G. Kwiat, “High efficiency single photon detection via frequency up-conversion,” J. Mod. Opt. **51**, 1433–1445 (2004).

**5. **M. A. Albota and F. C. Wong, “Efficient single-photon counting at 1.55 *μ*m by means of frequency upconversion,” Opt. Lett. **29**, 1449–1451 (2004) [CrossRef] [PubMed] .

**6. **R. V. Roussev, C. Langrock, J. R. Kurz, and M. M. Fejer, “Periodically poled lithium niobate waveguide sum-frequency generator for efficient single-photon detection at communication wavelengths,” Opt. Lett. **29**, 1518–1520 (2004) [CrossRef] [PubMed] .

**7. **Y. Ding and Z. Y. Ou, “Frequency downconversion for a quantum network,” Opt. Lett. **35**, 2491–2593 (2010) [CrossRef] .

**8. **B. Brecht, A. Eckstein, A. Christ, H. Suche, and C. Silberhorn, “From quantum pulse gate to quantum pulse shaper - engineered frequency conversion in nonlinear optical waveguides,” New J. Phys. **13**, 065029 (2011) [CrossRef] .

**9. **A. Eckstein, B. Brecht, and C. Silberhorn, “A quantum pulse gate based on spectrally engineered sum frequency generation,” Opt. Express **19**, 13770–13778 (2011) [CrossRef] [PubMed] .

**10. **M. G. Raymer, S. J. van Enk, C. J. McKinstrie, and H. J. McGuinness, “Interference of two photons of different color,” Opt. Commun. **283**, 747–752 (2010) [CrossRef] .

**11. **H. J. McGuinness, M. G. Raymer, C. J. McKinstrie, and S. Radic, “Quantum frequency translation of single-photon states in a photonic crystal fiber,” Phys. Rev. Lett. **105**, 093604 (2010) [CrossRef] [PubMed] .

**12. **L. Mejling, C. J. McKinstrie, M. G. Raymer, and K. Rottwitt, “Quantum frequency translation by four-wave mixing in a fiber: low-conversion regime,” Opt. Express **20**, 8367–8396 (2012) [CrossRef] [PubMed] .

**13. **G. Cerullo and S. De Silvestri, “Ultrafast optical parametric amplifiers,” Rev. Sci. Instrum. **74**, 1 (2003) [CrossRef] .

**14. **B. J. Smith and M. G. Raymer, “Photon wave functions, wave-packet quantization of light, and coherence theory,” New J. Phys. **9**, 414 (2007) [CrossRef] .

**15. **G. Strang, *Introduction to Linear Algebra*, 3rd ed. (Wellesley-Cambridge, 1998).

**16. **G. J. Gbur, *Mathematical Methods for Optical Physics and Engineering* (Cambridge, 2011).

**17. **S. L. Braunstein, “Squeezing as an irreducible resource,” Phys. Rev. A **71**, 055801 (2005) [CrossRef] .

**18. **M. T. Rakher, L. Ma, O. Slattery, X. Tang, and K. Srinivasan, “Quantum transduction of telecommunications-band single photons from a quantum dot by frequency upconversion,” Nat. Photonics **4**, 786–791 (2010) [CrossRef] .

**19. **R. Ikuta, Y. Kusaka, T. Kitano, H. Kato, T. Yamamoto, M. Koashi, and N. Imoto, “Wide-band quantum interface for visible-to-telecommunication wavelength conversion,” Nature Commun. **2**, 1544 (2011) [CrossRef] .

**20. **A. M. Brańczyk, A. Fedrizzi, T. M. Stace, T. C. Ralph, and A. G. White, “Engineered optical nonlinearity for quantum light sources,” Opt. Express **19**, 55–65 (2011) [CrossRef] .

**21. **J. Huang and P. Kumar, “Observation of quantum frequency conversion,” Phys. Rev. Lett. **68**, 2153–2156 (1992) [CrossRef] [PubMed] .

**22. **C. J. McKinstrie, L. Mejling, M. G. Raymer, and K. Rottwitt, “Quantum-state-preserving optical frequency conversion and pulse reshaping by four-wave mixing,” Phys. Rev. A **85**, 053829 (2012) [CrossRef] .

**23. **D. C. Burnham and R. Y. Chiao, “Coherent resonance flourescence excited by short light pulses,” Phys. Rev. **188**, 667–675 (1969) [CrossRef] .

**24. **Y. Huang and P. Kumar, “Mode-resolved photon counting via cascaded quantum frequency conversion,” Opt. Lett. **38**, 4 (2013) [CrossRef] .

**25. **M. G. Raymer and J. Mostowski, “Stimulated Raman scattering: Unified treatment of spontaneous initiation and spatial propagation,” Phys. Rev. A **24**, 1980–1993 (1981) [CrossRef] .

**26. **R. E. Giacone, C. J. McKinstrie, and R. Betti, “Angular dependence of stimulated Brillouin scattering in homogeneous plasma,” Phys. Plasmas **2**, 4596–4605 (1995) [CrossRef] .

**27. **R. V. Churchill, *Operational Mathematics*, 3rd ed. (McGraw-Hill, 1971).

**28. **M. G. Raymer, “Quantum state entanglement and readout of collective atomic-ensemble modes and optical wave packets by stimulated Raman scattering,” J. Mod. Opt. **51**, 1739–1759 (2004).