Mode Structure of a Broadband High Gain Parametric Amplifier

High gain parametric amplifier with a single-pass pulsed pump is known to generate broadband twin photon fields that are entangled in amplitude and phase but have complicated spectral correlation. Fortunately, they can be decomposed into independent temporal modes. But the common treatment of parametric interaction Hamiltonian does not consider the issue of time ordering problem of interaction Hamiltonian and thus leads to incorrect conclusion that the mode structure and the temporal mode functions do not change as the gain increases. In this paper, we use an approach that is usually employed for treating nonlinear interferometers and avoids the time ordering issue. This allows us to derive an evolution equation in differential-integral form. Numerical solutions for high gain situation indicate a gain-dependent mode structure that has its mode distributions changed and mode functions broadened as the gain increases.


I. INTRODUCTION
Parametric processes in the high gain regime are the most common and simple processes for generating a variety of quantum states including twin beam states, squeezed states, EPR-entangled states [1] for the applications in continuous variable quantum information processing, quantum communication, and quantum metrology [2]. For achieving high gain operation, ultra-short pulse pumping in single-pass configuration [3] is usually preferred due to its high concentration of energy and ease of operation. This has been done in optical waveguide structure [4] and optical fiber systems [5] in which spatial modes are well defined. However, due to ultrashort pulses in pumping and dispersion in nonlinear media, the spectral correlation is extremely complicated [6].
Fortunately, the complicated spectral correlation can be decomposed into independent temporal modes [7,8]. This was first pointed out by Law et al. [9] for the low gain case where two-photon events dominate and used in the analysis of the mode structure for the generation of high quality single-and two-photon states [10]. In the high gain regime, the existence of independent pairwise entangled temporal modes was recently confirmed experimentally in a direct measurement of the temporal mode profiles and in subsequent correlation measurement [11]. However, the common treatment of parametric interaction Hamiltonian [8] only works in the limit of low gain or in the regime of spontaneous emission but fails at high gain because it does not consider the issue of time ordering of the Hamiltonian and thus leads to the incorrect result that temporal modes do not change in the high gain limit [12]. Indeed, recent studies with approaches that avoid the time ordering issue showed the spectrum broadening as the gain increases [13]. The experiment that directly measured the temporal mode functions also confirmed the change of mode structure and mode func- * zou@iupui.edu tions as the gain increases [11].
The change of mode structure and mode functions with gain is troublesome in the production and applications of high quality quantum sources with quantum entanglement and noise reduction such as EPR entangled states and squeezed states, which require high gain operation in parametric processes. This is because the measurement on these states relies on the homodyne measurement technique in which the mode match between the local oscillator field (LO) and the quantum field is paramount and any mode mismatch is equivalent to losses and introduces extra vacuum noise. The knowledge of the exact profile of the mode functions will enable us to tailor the shape of LO field to match the quantum field [11]. But the change of the mode functions means that we also need to adjust the shape of LO to accommodate the change. The shape of the mode functions is also important for quantum pulse gates [14][15][16] in temporal mode multiplexing.
But so far there is no analysis about how mode structure and mode functions change with the gain. In this paper, we will investigate the pulse-pumped single-pass parametric processes at arbitrary pumping power. We will use an operator input-output approach that is usually employed for treating multi-stage nonlinear interferometers [17]. This avoids the time ordering issue of the interaction Hamiltonian and allows us to derive a set of coupled operator evolution equations in differentialintegral form. We solve them numerically and analyze the mode structure and mode functions at the final output ports as a function of the pump parameter.
The paper is organized as follows. In Sec.II, we introduce the input-output approach for the single-mode case of a multi-stage nonlinear interferometer involving parametric processes. In Sec.III, we apply the approach to the evolution of the fields in broadband parametric processes pumped by a pulse for an arbitrary length of nonlinear medium. We discuss mode decomposition in Sec.IV and solve numerically the state evolution in Sec.V to find how mode structure and mode functions change with the gain. We conclude with a discussion in Sec.VI. In order to reveal the issue of time ordering in the derivation of evolution operator in parametric processes and find ways to tackle it, we consider an SU(1,1) interferometer, shown in Fig.1, which consists of two parametric amplifiers (PA 1 , PA 2 ) characterized by gain parameters g 1 , g 2 together with a phase shift θ in between. This interferometer has recently been studied extensively [18] for precision phase measurement beyond standard quantum limit [19], quantum imaging with undetected photons [20], and quantum state engineering [21,22].
The two PAs are described by the Hamiltonians: where j = 1, 2. The input-output relations can be derived from evolution operatorsÛ j = e (1/i )Ĥ P A (ξj )t (j = 1, 2) and are respectively given aŝ (2) whereâ 1 =â 1 e iθ/2 ,b 1 =b 1 e iθ/2 , the amplitude gains g j ≡ (ξ j /|ξ j |) sinh |ξ j t|(j = 1, 2), and G j = cosh |ξ j t| for interaction time period of t. Here we assume the phase shifts are the same for both fields: θ a = θ b = θ/2. The outputs of the interferometer are then [17] with This shows that we can treat the whole system as one parametric amplifier with equivalent amplitude gains g T , G T . Furthermore, besides a propagation phase of e iθ/2 for both fields, the extra phase shift e iθ can be absorbed in g 2 by redefining g 2 ≡ g 2 e −iθ or ξ 2 ≡ ξ 2 e −iθ . Here ξ 2 = ξ 2 e −iθ takes the propagation phase shift e iθ into consideration. Notice that only when arg ξ 1 = arg ξ 2 − θ, other than a common phase of e iθ/2 for both fields, the whole system can be described by an equivalent overall HamiltonianĤ T =Ĥ P A (ξ T ) =Ĥ P A (ξ 1 ) +Ĥ P A (ξ 2 ) with ξ T = ξ 1 + ξ 2 . But if arg ξ 1 = arg ξ 2 − θ, thenĤ T =Ĥ P A (ξ 1 ) +Ĥ P A (ξ 2 ). This is because This will have an inconvenient consequence when we extend the interferometer to multiple stages [10,22], as shown in Fig.2. Because the phases in each stage are arbitrary, we therefore cannot write the overall Hamiltonian as the sum of each stage: where the phase shifts at each stage are absorbed in the interaction parameter ξ j = ξ j e −iθj .
To solve this problem, we can proceed by using repeatedly Eq.(4) to add each stage and obtain a recursive relation. Specifically for Eq.(4), we treat all the stage added up to stage k as the first PA with equivalent amplitude gains G T (k), g T (k) and the second PA is the k + 1-th stage to be added: Unfortunately, we cannot find an analytical expression for the final outputs. In order to have some general idea about the outputs, we consider each stage has an infinitesimally small gain and phase shift whose sizes are proportional to an infinitesimal length scale ∆x along the field propagation direction: g k ≈ ζ(x)∆x, θ k ≈ η(x)∆x and we use location x = k∆x to denote the k-th stage. When ∆x → 0, G k+1 = 1 + |g k+1 | 2 ≈ 1 + o(∆x). So, Eq.(6) can be approximated as These are the evolution equations for a parametric amplifier with continuous gain function ζ(x). The phase parameter η(x) usually corresponds to phase mismatching. The initial condition is obviously G T (0) = 1, g T (0) = 0. It is hard to solve analytically the differential equations if ζ(x), η(x) depend on location x. For simplicity, let us assume ζ, η be constant. Then, Eq.(8) can be solved analytically and have the following solution: If ζ ≤ η/2 and η 0 ≡ η 2 − 4|ζ| 2 , we have If ζ ≥ η/2 and ζ 0 ≡ |ζ| 2 − (η/2) 2 , we have The exponential growth of the gain when ζ ≥ η/2 is typical of high gain parametric amplifiers. But the oscillatory low gain behavior when ζ ≤ η/2 is the result of interference as we will see in the following.
In the limit of ζ 1, we have The extra phase factor e iηx/2 extracted out of g T is for the consistency with G T and is due to propagation of the fields through the system. The gain parameters in Eq.(11) are equivalent to an overall Hamiltonian of interaction parameter ξ T ≡ ζe −iηx/2 sinc(ηx/2) with evolution time t replaced by x: which can thought of as the sum of all the stages: Here, dξ = ζe −iηx dx is the infinitesimal gain parameter for each infinitesimal stage. Note that this equivalence is true only when ζ(x) = constant and η(x) = constant. As a matter of fact, when gain parameters |g 1 |, |g 2 | 1, we can add the two Hamiltonian to obtain overall Hamiltonian:Ĥ T =Ĥ P A (ξ T ) =Ĥ P A (ξ 1 ) +Ĥ P A (ξ 2 ). This can be seen from evolution operator Here, we assumed |g 1 | = |ξ 1 |∆x 1, |g 2 | = |ξ 1 |∆x 1 and ξ 1 = ξ 1 e −iθ with phase shift e −iθ included. So, when the amplitude gains |g j | 1, we can simply add the Hamiltonian of each stage: and The above can also be thought of as a result of twophoton interference among the pair of photons generated by each stage [10,22]. This can be confirmed by looking at the output state for vacuum input: with |Ψ k as the two-photon state generated by k-th stage.
Notice that Eq. (15) is true only if the overall amplitude gain g T is much smaller than 1 so that the last step of Eq. (14) stands. For the high gain case, Eq.(15) does not stand and we have to resort to Eq.(6) or Eq.(10). This will pose serious problem in finding solution for a broad band parametric amplifier in the high gain regime.

III. BROADBAND PARAMETRIC AMPLIFIER
When parametric processes are pumped by high power pulses, broadband parametric amplification is achieved. They can be used to produce quantum entangled fields with a wide bandwidth. The traditional treatment of this situation is to start with a multi-mode Hamiltonian of the form [6,8] where subscript "M" denotes multi-mode, χ is some parameter proportional to the nonlinear coefficient of nonlinear medium of length L 0 , A p (ω 3 ) is the spectral amplitude of the pump field, and Ψ(ω 1 , ω 2 , ω 3 ) is obtained from spatial integration: with β ≡ ∆kL 0 /2 and ∆k ≡ k 1 + k 2 − k 3 as the phase mismatching. We then find evolution operator aŝ where time integration gives rise to a delta-function δ(ω 1 + ω 2 − ω 3 ) and the integrated Hamiltonian has the form of with G ≡ χ/C as a dimensionless gain parameter such that is normalized: dω 1 dω 2 |Φ(ω 1 , ω 2 )| 2 = 1. In general, this gives rise to a complicated coupling of different frequency components at the outputs: are the outputs at the end of nonlinear medium andâ(ω 1 ),b(ω 2 ) are those at the start. But we can make a singular value decomposition (SVD) of the joint spectrum function (JSF) Φ(ω 1 , ω 2 ): where {ψ k , φ k } are two sets of ortho-normal functions: and {r k } are non-negative numbers satisfying k r 2 k = 1. Then Eq.(21) becomes Together with Eq. (20), this leads to de-coupling of the different temporal modesÂ k andB k for the two output fields: Unfortunately, it was pointed out [12] that the evolution operator in Eq. (20) is not correct for the Hamiltonian in Eq. (18) The reason is the same as those for Eq.(1). But Eqs. (23,24) are still correct for the Hamiltonian in Eq. (18).
In order to treat this in a correct manner, we can apply the same approach in previous section. Consider a nonlinear medium of length L which is divided into a small segment of size ∆L, as shown in Fig.3. Let us treat the small segment first. The Hamiltonian is given from Eq. (18) for the small segment aŝ (28) where the spatial integration starts at z instead of 0 because the ∆L segment is located at z inside the medium and we assume that there is no pump depletion. For ∆L = dL → 0, we havê The evolution operator for this segment is given by the Dyson series:Û Eq.(30) becomes Eq. (20) if [Ĥ(z, dL, t 1 ),Ĥ(z, dL, t 2 )] = 0 but it is not true forĤ(z, dL) in Eq. (29). On the other hand, for dL → 0, we havê So, we obtain the evolution of the field operators: where the time integral gives rise to a δ-function for ω 1 + ω 2 − ω 3 and ∆k = ∆k |ω3=ω1+ω2 . With dâ( Likewise, we obtain and From Eqs.(36,37) and initial conditions in Eq.(38), we can verify that These relations guarantee the commutation relations

IV. EIGEN-MODES OF HIGH GAIN PARAMETRIC PROCESSES
Although h-functions have very complicated form, we can in general use singular value decomposition method to decompose them as Because of relations in Eq.(39), it can be shown (see Appendix) that ψ k (ω) are four sets of orthonormal mode functions satisfying dωψ (ω) = δ kk and r k 's are normalized mode coefficients satisfying k r 2 k = 1 with G being some parameter depending on χ. So, h-functions are in the form of With these relations and orthonormal relations for ψ (a,b) k (ω), Eqs. (23,24) can be recast aŝ kb (ω) define the annihilation operators for the corresponding output and input temporal modes, similar to those given in Eq. (27). Because of orthonormal relations, they satisfy the Boson commutation relation for annihilation operators.

V. NUMERICAL SOLUTIONS
Let us use a nonlinear fiber as the nonlinear medium. We use the parameters similar to those given in Ref. 24 for a real piece of 300m-long dispersion-shifted nonlinear fiber. We obtain the dimensionless parameters 1/∆ 1 = 0.785, 1/∆ 2 = −0.471. The numerical solution of the amplitude of h 2a (Ω 1 , Ω 2 ) is shown in Fig.4 for four values of pump parameter K: K = 0.01, 2, 4, 10. For small K( 1, Fig.4(a)), h 2a (Ω 1 , Ω 2 ) is exactly the joint  The change of the transfer function h 2a (Ω 1 , Ω 2 ) with pump parameter K will lead to mode structure change. This is reflected in the change of the distribution of the mode coefficients {r k }, whose normalized values to r 1 are plotted in Fig.6 as a function of the pump parameter K. The trend shows the increasing weight of the higher order modes in addition to the broadening of the mode functions as K increases.
The multi-mode nature and the broadening of the mode functions with pump parameter K are due to the values of ∆ 1 , ∆ 2 in ∆kL 0 for a realistic nonlinear fiber case, which gives rise to an asymmetric sinc-function in Φ of Eq. (22) or h 2a (Ω 1 , Ω 2 ) at small K (Fig.4(a)). In principle, we can adjust the dispersion parameters of the fiber to change ∆ 1 , ∆ 2 in ∆kL 0 . It is found that the initial h 2a (Ω 1 , Ω 2 ) when K 1 is nearly round or factorized with parameters 1/∆ 1 = 2.198, 1/∆ 1 = −2.198, as shown in h 2a (Ω 1 , Ω 2 ) in Fig.7(a). The initial r k distribution is indeed close to single mode with high order r k much smaller than 1. This can be seen in Fig.8 for K = 0. However, the trend of mode spreading for large K in Fig.6 persists in Fig.8.
On the other hand, a close look at Fig.7(c,d) for large K shows that function h 2a is still pretty round or nearly single mode. In fact, we find from mode decomposition in Eq.(41) that the coefficient of each mode is sinh 2 (r k G), which becomes 0.25e 2r k G for large r k G. So, the ratio of coefficients of the first mode to higher mode is then e 2(r1−r k )G 1 for large G, that is, the first mode will dominate in the mode decomposition in Eq.(41) at large K (or G) [25], which leads to a nearly round (or factorized) h 2a .
Furthermore, we can look at the mode purity of the output state. It is known that if we discard the other  Fig.4 (orange, (a)) and Fig.7 (blue,  (b)). The dashed line corresponds to M = 1.
field, one of the output fields from a parametric amplifier is in a thermal state. For broadband pulsed pumping, it becomes a multi-mode thermal field. Mode purity for a pulsed multi-mode thermal field was studied in Ref. 26 with the temporal mode format described here. It was found that mode purity can be characterized by the normalized intensity correlation function g (2) defined as where M is defined as the number of modes (M = 1 gives the pure single-mode case) and I k is the intensity of mode k. From Eq.(42), we find that I k = sinh 2 (r k G) for vacuum input to the parametric amplifier. So, Eq.(47) becomes Obviously, for the single mode case with r 1 = 1, r k = 0 (k > 1), we have M = 1. Equation (48) was first derived by Christ et al. [8] for temporal modes and by Sharapova et al. [27] and Dyakonov et al. [28] for spatial modes in high gain parametric processes. For the two cases shown in Fig.4 and Fig.7, we evaluate mode number M in Eq.(48) as a function of pump parameter K and plot the results in Fig.9. It can be seen that mode number M indeed drops first as K increases. This drop is consistent with the dominance of the first mode as K increases. However, the drop stops at around K = 2 when a minimum of M is reached. M starts to slowly increase after K > 2. This is due to the increase of r k for higher order modes as K increases, as shown in Fig.6 and Fig.8.
The minimum M value depends on the initial M at low K( 1). Thus, it is better to have a nearly single mode situation at low pump power.

VI. SUMMARY AND DISCUSSION
We studied the mode structure for a broadband parametric amplifier at different gains by using an inputoutput approach that avoids the crucial issue of timeordering in Hamiltonian. Contrary to previous studies where the time-ordering issue was not treated, the mode structure changes as the gain increases in the sense that both the mode distribution and the mode functions broaden with the increase of the gain. Although the mode number, a quantity that characterizes the total number of modes, drops initially as the gain changes from low to high due to the dominance of the first mode, it reaches a minimum value before slowly increases due to the broadening of the mode distribution.
The mode structure change with the gain will have profound impact on the application of broadband parametric processes in quantum technology with continuous variables, which relies on homodyne detection method as the dominant measurement technique. Even though the output of the broadband parametric processes is of multimode nature, because of the orthogonality of the modes in the processes, we can choose the local oscillators of homodyne measurement to match the mode functions of the output and select specific modes for study [11]. The change of mode structure and functions with the gain means that we will need to measure them constantly at different gains as we change the operation condition. Fortunately, methods are recently developed for the direct measurement of mode functions and mode coefficients of parametric processes [11,30].

ACKNOWLEDGMENTS
This work was supported by US National Science Foundation (Grant No. 1806425).

VII. APPENDIX
In general, we can use singular value decomposition to write where ψ Using the completeness of mode function {ψ , we can rewrite the first relation in Eq.(39) as Treating k as the row index and ω as the column index of a matrix, we can consider mode function ψ (k) (ω) as matrix element of Ψ with {Ψ} k,ω ≡ ψ (k) (ω) (we dropped subscript 1a, 2a for clarity). Then Eq.(51) is equivalent to the following matrix equation: Multiplying both sides with ψ (k1) * 1a (ω)ψ (k2) * 1b (ω ) and integrating ω, ω , with orthonormal relations for ψ where we set r b , and R k2,k1 ≡ dωφ 1a (ω). Switching back notations: k = k 1 , k = k 2 , we have which, in matrix form, is simply R = C b R C −1 a or R = C Now notice that RR † = I = R R † , or using matrix form of Eq.(56), we obtain Rewrite the above, we have Since matrices C a , C b are both diagonalized, the above expressions are true only if R = I = R, or φ k (ω) in Eq.(49), we obtain Eq.(41).