Highly Squeezed States in Ring Resonators: Beyond the Undepleted Pump Approximation
Colin Vendromin Yan Liu Zhenshan Yang John E. Sipe


  
APS/123-QED



colin.vendromin@utoronto.ca
Department of Physics, University of Toronto, 60 St. George Street, Toronto, Ontario, Canada M5S 1A7

School of Physics Science and Information Technology, Shandong Key Laboratory of Optical Communication Science and Technology, Liaocheng University, Liaocheng 252059, China

School of Physics Science and Information Technology, Shandong Key Laboratory of Optical Communication Science and Technology, Liaocheng University, Liaocheng 252059, China

Department of Physics, University of Toronto, 60 St. George Street, Toronto, Ontario, Canada M5S 1A7




             

We present a multimode theory of squeezed state generation in resonant systems valid for arbitrary pump power and including pump depletion. The Hamiltonian is written in terms of asymptotic-in and -out fields from scattering theory, capable of describing a general interaction.  As an example we consider the lossy generation of a highly squeezed state by an effective second-order interaction in a silicon nitride ring resonator point-coupled to a waveguide. We calculate the photon number, Schmidt number, and the second-order correlation function  of the generated state in the waveguide. The treatment we present provides a path forward to study the deterministic generation of non-Gaussian states in resonant systems.



                         
Highly Squeezed States in Ring Resonators:

Beyond the Undepleted Pump Approximation
    J. E. Sipe
    April 29, 2024
====================================================================================




§ INTRODUCTION

Continuous-variable (CV) squeezed states of light  exhibit an uncertainty below vacuum noise in one quadrature <cit.> making them an essential resource for implementing CV quantum information processing <cit.>, and  with post-selection strategies they can be used to create
non-Gaussian states of light <cit.>, an important resource for quantum computing. 

Squeezed states are typically generated through nonlinear optical processes such as spontaneous parametric down conversion (SPDC) and spontaneous four-wave mixing (SFWM), the first relying on the presence of a χ^(2) susceptibility, and the second on the χ^(3) susceptibility <cit.>.  In the usual treatment of these processes, the pump field or fields are treated classically, with the nonlinear process being considered sufficiently weak that those fields can be treated as undepleted. 









In this case the resulting squeezed state can be approximately written as <cit.>

    |ψ⟩= exp(β/2∫ dk_1 k_2 ϕ(k_1, k_2) a^†_k_1a^†_k_2 - H.c.)|vac⟩,

where  β is an overall squeezing parameter, and ϕ(k_1, k_2) is the normalized biphoton wavefunction that describes the production of pairs of photons; one with wavenumber k_1 and the other k_2. The ket in Eq. (<ref>) is equivalent to a lowest-order Magnus expansion <cit.> that treats the pump as classical and  undepleted. If the squeezing parameter β is sufficiently weak, then time-ordering corrections that arise in the higher-order terms of the Magnus expansion are expected to be negligible. However, time-ordering effects  become important when the amount of squeezing increases <cit.>, and the approximate ket in Eq. (<ref>) then becomes a poorer and poorer approximation.

Including time-ordering corrections, the ket can obtained
under the assumption that the pump is classical and undepleted,  and that the generated light is Gaussian. As a result, the full Gaussian ket can be written generally as a squeezed vacuum state

    |ψ⟩ = exp(1/2∫ dk_1 dk_2 J(k_1,k_2)a^†_k_1a^†_k_2 -  H.c.)|vac⟩,

where the squeezing matrix J(k_1, k_2) can be obtained from the Heisenberg equations of motion <cit.>. If time-ordering corrections are negligible this ket reduces to  Eq. (<ref>), with  the squeezing matrix just being proportional to the biphoton wavefuction  J(k_1, k_2) ≈βϕ(k_1, k_2).

However, if the nonlinear process is strong enough that pump depletion needs to be considered, the generated light can be expected to be non-Gaussian due to entanglement between the pump and the generated fields <cit.>. In this case a squeezed vacuum state such as in Eq. (<ref>) does not provide a complete description of the features of the generated light, such as pump depletion. To go beyond the undepleted pump approximation  one strategy that has been employed is to replace the vacuum state in Eq. (<ref>) with a non-Gaussian ket, which starts from vacuum initially and evolves such that it remains close to the vacuum. This treatment was recently investigated in waveguides <cit.> and for  a single-mode SPDC Hamiltonian <cit.>, and offers a potential route to the deterministic generation of non-Gaussian states of light. In fact experiments have already demonstrated high conversion efficiencies where the undepleted pump approximation is no longer valid <cit.>, and non-Gaussianity may become important.

In this paper  we extend such a treatment to   resonant systems. 
As an example we study a ring resonator, a system of interest in its own right due to its resonant enhancement of nonlinear processes and the precise control over the frequencies involved, as well as its integrability on photonic chips <cit.>. We develop a multimode Hamiltonian theory of non-Gaussian state generation in such systems, including scattering loss, by employing
the asymptotic-in and -out fields used in scattering theory <cit.>. 

We apply the theory to the particular case of generating a highly squeezed state by an effective χ^(2) interaction in a ring resonator <cit.>. For this interaction one relies on the χ^(3) response of the medium, but there is one pump field in a mode C that is continuous wave (CW) and so strong that it can be treated as classical and undepleted; in the presence of the material χ^(3) this leads to an effective χ^(2) through which a much weaker pump in the mode P can generate two photons in  the mode S (see Fig. <ref>(a)).  There is no issue with phase-matching, since the pumps and the generated photons are close in wavelength. The formalism we present can also be applied to a standard χ^(2) interaction as well.

Although the theory we present is capable of studying non-Gaussian states, our focus in this paper is on highly squeezed states that include the effects of pump depletion. The  treatment of the non-Gaussian corrections to the squeezed state will be handled in a subsequent communication, as will the effects of self- and cross-phase modulation, and the excitation of fields in additional ring resonances. While in this paper we specialize to ring resonator systems, most of the formalism can be applied to other systems as well, such as waveguides.




This paper is organized as follows. In Sec. <ref> we write
the system Hamiltonian, including the SFWM interaction in a ring resonator, in terms of the asymptotic-in and -out fields for our system. In Sec. <ref> we introduce the input and output kets  that are used to describe the generation of light in the ring resonator system. In Sec. <ref> we derive the non-Gaussian unitary operator that connects the input ket to the output ket. In Sec. <ref> we show how to obtain the first-order solution for the ket. In Sec. <ref> we numerically study the generation of a highly squeezed state by an effective χ^(2) interaction in a ring resonator with loss. In Sec. <ref> we conclude.  




§ THE HAMILTONIAN


Our system of interest is a ring resonator point-coupled to an actual waveguide and a phantom waveguide (see Fig. <ref>(b)), where the phantom waveguide is used to describe scattering losses <cit.>. We define the “input and  output channels" of the actual waveguide (shown in red) as the parts of that 
waveguide that are to the left and to the right of the coupling point, respectively; likewise, the input and output channels of the phantom waveguide (shown in white) are the parts of that  waveguide  that are
to the right and to the left of the coupling point, respectively. 








As a basis for expanding the electromagnetic field in the full system of waveguides and ring, we employ the asymptotic-in and -out mode fields <cit.>; they are sketched in Fig. <ref>(c), and are eigenstates of the full linear Hamiltonian of the system. Each asymptotic-in or -out mode field has an amplitude in both waveguides and in the ring. In each diagram, the black arrow corresponds to a freely propagating field that is either incoming (asymptotic-in) or outgoing (asymptotic-out), with its amplitude equal to that of a mode field in an isolated waveguide. The white arrows indicate either outgoing (asymptotic-in) or incoming (asymptotic-out) fields in the other channels, and their amplitudes depend on the coupling parameters between the ring and the waveguides <cit.>. Asymptotic-in mode fields are also defined for the phantom channel, but they will not appear in our calculations.



The asymptotic-in and -out mode fields can be identified by three parameters. We use a discrete index n to indicate the actual output channel or the phantom output channel, the waveguide wavenumbers
k, and an index J that identifies the range of wavenumbers k in the  
waveguide that correspond to a range of frequencies near a ring resonance, 
where J=P,  S,  C. We assume the ring resonances have
sufficiently high quality factors that the ranges over which k varies for each J can be considered distinct; we leave this implicit in the integrals below.

For the system considered here, the displacement field operator 
can be written in terms of the asymptotic-in mode fields  as
<cit.>

    D( r)    = ∑_n, J∫ dk  b_n J k D^in_n J k( r) +  H.c.,

where b_nJk is the operator that annihilates an asymptotic-in photon associated with input channel n, with wavenumber k in the ring resonance  J, and D^in_n J k( r) is the mode field associated with that photon (see the left sketch in Fig. <ref>(c) for the asymptotic-in mode field associated with the actual input channel). Similarly, the displacement field operator can be written in terms of the asymptotic-out mode fields as 

    D( r)    = ∑_n, J∫ dk a_n J k D^out_n J k( r) +  H.c.,

where a_nJk is the operator that annihilates an asymptotic-out photon associated with output channel n, with wavenumber k in the ring resonance J, and D^out_n J k( r) is the mode field associated with that photon (see the middle or right sketch in Fig. <ref>(c)). The input and output operators satisfy the commutation relations 

    [a_n J k, a^†_n J k'] = [b_n J k, b^†_n J k'] = δ(k-k'),

 with all others zero; the linear Hamiltonian can be written in terms of the input or output operators as

    H_L   = ∑_n,J∫ dk ħω_n J kb^†_n J k b_n J k
       = ∑_n, J∫ dk ħω_n J ka^†_n J k a_n J k,

where ω_nJk is the dispersion relation for the waveguide containing channel n, ranging only over wavenumbers k near ring resonance J; see Appendix <ref> for details.







We now include a third-order nonlinear interaction described by
the Hamiltonian <cit.>

    H_NL = -1/4ϵ_0∫ d r Γ^ijkl_3( r) D^i( r) D^j( r)D^k( r)D^l( r),

where i, j, k, and l are Cartesian components, as usual summed over when repeated, and the tensor Γ^ijkl_3( r) is related to the more familiar third-order tensor χ^ijkl_3( r) by <cit.>
    Γ^ijkl_3( r) = χ^ijkl_3( r)/ϵ^2_0ϵ^4( r),

where we have assumed that the relative dielectric function ϵ( r) does not depend on frequency  over the range of frequencies within each ring resonance, and have neglected any dependence of χ_3^ijlk( r) over the frequencies of interest.

Our process of interest is SFWM, where one
pump photon in resonance P and one in resonance C are annihilated, and two photons in
resonance S are generated  (see Fig. <ref>(a)). The asymptotic-out fields for the actual and phantom waveguides are used for the generated photons, while only the asymptotic-in mode field for the actual waveguide is used for the pump photons; hence we can drop the index n on the asymptotic-in mode fields and operators, understanding the actual input channel to be identified. If we consider only the terms responsible for SFWM in Eq. (<ref>), the nonlinear Hamiltonian becomes <cit.>

    H_NL   = -3/ϵ_0∑_n,n'∫ dk_1dk_2dk_3dk_4 K^nn'(k_1,k_2,k_3,k_4) 
       × a^†_n S k_1a^†_n' S k_2 b_P k_3b_C k_4 +  H.c.,

where the nonlinear coefficient for SFWM is given by

    K^nn'(k_1,k_2,k_3,k_4)    = ∫ d rΓ^ijkl_3( r)  [ D^out,i_n S k_1( r)]^*
       ×[D^out,j_n' S k_2( r)]^* D^in,k_P k_3( r)D^in,l_C k_4( r).

We take the integral in Eq. (<ref>) to range over the ring where the fields are concentrated, and the effect of the nonlinearity is most significant.
The sums over n and n', which range over the actual
and phantom output channels, 
take into account the four combinations of two photons exiting the ring through two output channels. That is,  both photons can leave through the phantom or actual output channel, or one through each of the channels.
The nonlinear coefficient K^nn'(k_1,k_2,k_3,k_4) gives the strength of the interaction corresponding to each combination. It is symmetric with respect to exchanging
the order of both n and n' and k_1 and k_2; that is,

    K^nn'(k_1,k_2,k_3,k_4) = K^n'n(k_2,k_1,k_3,k_4).

But it is not symmetric with respect to only exchanging the order of n and n' or k_1 and  k_2; that is,

    K^nn'(k_1,k_2,k_3,k_4) ≠ K^n'n(k_1,k_2,k_3,k_4).









§ INPUT AND OUTPUT KETS





The  unitary time-evolution operator for the system, 𝒰(t,t'), is a solution to the Schrödinger equation

    iħd𝒰(t,t')/dt = H 𝒰(t,t'),

satisfying 𝒰(t',t')=I for all t', where I is the identity operator. 
We write the  full Hamiltonian H as the sum of the linear and nonlinear contributions,

    H = H_L + H_NL,

where 
H_L  is given by Eq. (<ref>) and H_NL is given by Eq. (<ref>). We assume that times t_on and t_off can be identified such that:  (a) For t<t_on the excitations of the quantum fields incident on the ring are sufficiently far from it that H_NL has no effect on the evolution of the ket, and that (b) for t>t_off all the field excitation has propagated sufficiently far from the ring that, again, H_NL has no effect on the evolution of the ket.

Consider then the evolution of a ket from some early time t_0 < t_on to some later time t_1 >t_off. The evolution operator can be written as

    𝒰(t_1,t_0) =  𝒰(t_1,t_off) 𝒰(t_off,t_on) 𝒰(t_on,t_0).

 The evolution operator from t_0 to t_on and  t_off to t_1 only involves the linear Hamiltonian, that is

    𝒰(t_on,t_0)    =  e^-iH_L(t_on - t_0)/ħ,
      
    𝒰(t_1,t_off)    =  e^-iH_L(  t_1-t_off)/ħ.

We can then write Eq. (<ref>) as

    𝒰(t_1,t_0) =  e^-iH_L(  t_1-t_off)/ħ𝒰(t_off,t_on) e^-iH_L(t_on - t_0)/ħ.




We now suppose that an initial ket |ψ(t_0)⟩ is specified at a time t_0<t_on. We define the “input ket"

    |ψ_in⟩≡ e^-iH_L(0 - t_0)/ħ|ψ(t_0)⟩,

which is the ket that the initial ket would evolve to at t=0 were H_NL = 0. And starting with the ket |ψ(t_1)⟩ characterizing the state at a time t_1>t_off, we define  the “output ket"

    |ψ_out⟩≡ e^iH_L(t_1- 0)/ħ|ψ(t_1)⟩,

which is the ket at t=0 that would evolve to  
|ψ(t_1)⟩ at t=t_1 were H_NL = 0. Since the evolution of the ket from from t_0 to t_1 is given by

    |ψ(t_1)⟩ =  𝒰(t_1,t_0)|ψ(t_0)⟩,

in terms of the input and output kets
this can be written as

    |ψ_out⟩ = U(t_1,t_0) |ψ_in⟩,

where we defined

    U(t_1,t_0) ≡ e^iH_Lt_1/ħ𝒰(t_1,t_0)  e^-iH_Lt_0/ħ.

Then taking the limit

    U(∞, -∞) = lim_t_1 →∞, t_0 → -∞ U(t_1,t_0),

we obtain

    |ψ_out⟩ = U(∞,-∞) |ψ_in⟩.

It is then convenient to define

    U(t) ≡ U(t,-∞),

for all times t. We can then write

    |ψ_out⟩ = |ψ(∞)⟩,

where

    |ψ(t)⟩   = U(t) |ψ_in⟩,

and

    |ψ_in⟩=|ψ(-∞)⟩.

Taking the derivative of Eq. (<ref>) and using Eqs. (<ref>), (<ref>), and (<ref>) we obtain

    i ħd/dt  |ψ(t)⟩ = H_NL(t) |ψ(t)⟩,

where we have defined the interaction picture nonlinear Hamiltonian as

    H_NL(t)≡ e^iH_Lt/ħ H_NL e^-iH_Lt/ħ,

and as t ranges from -∞ to ∞ in Eq. (<ref>) the ket evolves from |ψ_in⟩ to |ψ_out⟩.
 





§ UNITARY EVOLUTION


In Appendix A we simplify the nonlinear Hamiltonian Eq. (<ref>) by treating the strong CW pump classically, and discretize the k involved for the weak pump and generated fields to implement numerical calculations. The result can be written as

    H_NL(t)   =  ħ∑_μ, ν∑_lΛ_μν l(t) a^†_μa^†_ν b_l +  H.c.,

where for convenience we removed the discrete labels P and S on the input and output operators, and they  satisfy the usual commutation relations [b_l, b^†_l'] = δ_ll'  and [a_μ, a^†_ν] = δ_μν. Here Λ_μν l(t) is the nonlinear coefficient for the effective χ^(2) interaction, where l labels a discrete wavenumber k_l for the pump, and μ and ν are discrete indices for the generated photons (see Eq. (<ref>)). Each Greek index μ, for example, refers to two properties of the photon, its discrete wavenumber k_i and the output channel n. The sum over μ in Eq. (<ref>) represents two sums; one over the actual and phantom output channel and the second over the discrete wavenumbers k_i. 

We now turn to the unitary evolution in Eq. (<ref>). Here we do not consider a general ket |ψ_in⟩, but restrict ourselves to
an input ket that is a coherent state in the actual input channel and vacuum in all other channels,

    |ψ_in⟩   = D_b(β_in)|vac⟩,

where  D_b(β(t)) is the unitary displacement operator, with displacement parameter β_in≡β(-∞). It is written in terms of the input operators as

    D_b(β(t)) = exp([β(t)]^ T b^† - [β^*(t)]^ T b),

where β(t) = [β_1(t), β_2(t),…]^ T and β^*(t) = [β^*_1(t), β^*_2(t),…]^ T  are column vectors of displacement parameters, and b = [b_1, b_2, …]^ T and b^† = [b^†_1, b^†_2, …]^ T are column vectors of input operators.  We then seek a solution of (<ref>)
of the form 


    |ψ(t)⟩=U_a(t)U_b(t)|ψ̃(t)⟩,

where |ψ̃(t)⟩ describes the effect of a non-Gaussian unitary operator acting on the vacuum state, as discussed below, and
U_a(t) and U_b(t) are Gaussian unitary operators, given by

    U_a(t)    = S_a( J(t))R_a(ϕ(t)) e^iθ_a(t),
        
    
        U_b(t)    = D_b(β(t)) e^iθ_b(t),

where θ_b(t) and θ_a(t) are real functions of time,  S_a( J(t)) is a multimode squeezing operator, and   R_a(ϕ(t)) is a multimode rotation operator.

The unitary squeezing operator S_a( J(t)) is given by <cit.>
 
    S_a( J(t))    = exp(1/2 a^†^ T J(t)  a^† -  H.c.),

where a = [a_1, a_2, …]^ T and a^† = [a^†_1, a^†_2, …]^ T are column vectors of output operators and J(t) is the symmetric time-dependent squeezing matrix, with J(t) =  J^ T(t). The unitary rotation operator R_a(ϕ(t)) is given by <cit.>
 
    R_a(ϕ(t))    = exp(i  a^†^ Tϕ(t)  a),

where ϕ(t) is
the Hermitian time-dependent rotation matrix,
with ϕ (t) = ϕ^†(t). The parameters β(t), J(t), and ϕ(t) are subject to the following initial conditions β(-∞) = β_in, J(-∞) = 0, and ϕ(-∞) = 0, and also θ_b(-∞) = θ_a(-∞) = 0; further, we 
take

    |ψ̃(-∞)⟩ = |vac⟩.

This guarantees that the initial condition (<ref>) is satisfied.

Returning to Eq. (<ref>), our goal is to define the parameters in the Gaussian unitary operators in Eq. (<ref>) and Eq. (<ref>) so that the ket |ψ̃(t)⟩ isolates the non-Gaussian behavior of the full ket in Eq. (<ref>). Now regardless of how the parameters in those Gaussian unitary operators are chosen to depend on time, we find that for Eq. (<ref>) to satisfy the full Schrödinger equation in Eq. (<ref>) we require

    i ħd/dt |ψ̃(t)⟩ = H_eff(t) |ψ̃(t)⟩,

where we have introduced the effective Hamiltonian

    H_eff(t)   =H_I(t) + U^†_a(t)H_a(t)U_a(t) -  iħ U^†_a(t)dU_a(t)/dt
       +U^†_b(t)H_b(t) U_b(t) -iħ U^†_b(t)dU_b(t)/dt,

with

    H_a(t)    ≡ħΛ_μν l(t)[β(t)]_la^†_μ a^†_ν +  H.c.,
        
    
        H_b(t)   ≡ħΛ^*_μν l(t)  [ V(t) W^ T(t)]_μν(b^†_l - [β^*(t)]_l)+  H.c.,
        
    
    H_I(t)    ≡ħΛ^*_μν l(t) ( [ V(t)]_μμ' [ V(t)]_νν' a_μ' a_ν'
       + ([ V(t)]_μν' [ W(t)]_νμ'  + [ W(t)]_μμ' [ V(t)]_νν') a^†_μ'a_ν'
       + [ W(t)]_μμ' [ W(t)]_νν' a^†_μ'a^†_ν')b^†_l  + H.c..

Here repeated indices are summed over, and the matrices V(t) and W(t) are defined in terms of the parameters appearing in the Gaussian unitary operators in Eq. (<ref>) and Eq. (<ref>); see Appendix B for details.

The Hamiltonians appearing above each capture part of the physics: 
H_a(t) describes the generation of photon pairs by an effective χ^(2) interaction modulated by a time-dependent pump amplitude β(t); 
H_b(t) determines the rate of depletion of the pump amplitude β(t), which is being driven by the generation of photon pairs; 
H_I(t), which  only contains non-Gaussian terms, describes the quantum correlations between the pump and the generated light. 

To achieve our goal of moving all the non-Gaussian behavior to |ψ̃(t)⟩, we must choose the time dependence of the parameters appearing in the Gaussian operators in Eq. (<ref>) and Eq. (<ref>) so that
 
    H_eff(t) =H_I(t).

It is easy to show that this condition is satisfied if we have
 
    iħdU_a(t)/dt   =H_a(t)U_a(t),
     
    
    iħdU_b(t)/dt   = H_b(t)U_b(t).
   
We show in Appendix <ref> that this can be achieved by
appropriately choosing the parameters V(t), W(t), β(t), and θ_b(t) so that

    d V(t)/dt   =-2iζ(t) W^*(t),
        
    d W(t)/dt   =-2iζ(t) V^*(t),
           
    
        id[β(t)]_l/dt   = [γ(t)]_l, 
        
    θ_b(t)   =1/2∫_-∞^tdt'[γ(t')]_l[β^*(t')]_l +  c.c.,

are satisfied, where the matrix ζ(t) drives the production of photon pairs in the squeezed state,

    [ζ(t)]_μν =  Λ_μν l(t)[β(t)]_l,

and the vector γ(t) determines the rate of pump depletion,

    [γ(t)]_l = Λ^*_μν l(t)  [ V(t) W^ T(t)]_μν.

 
 We solve the Eqs. (<ref>), (<ref>), and (<ref>) numerically for the parameters V(t), W(t), and β(t), which are subject to the initial conditions β(-∞) = β_in,  W(-∞) = 0, and V(-∞) =  I, where I is the identity matrix. The solutions are 
 put into Eq. (<ref>), which is then integrated to obtain the phase θ_b(t). In Appendix <ref> we show how to extract the squeezing matrix J(t) and rotation matrix ϕ(t) from the solutions V(t) and W(t), and in Appendix <ref> we give an equation for the phase θ_a(t) in Eq. (<ref>). This procedure provides a solution for the Gaussian parameters of the unitary operators in Eqs. (<ref>) and (<ref>), and is valid for high pump power and nonlinear interaction strength. 

Still necessary for a full solution, of course, is to include the non-Gaussian evolution of |ψ̃(t)⟩. Then the full solution can be written as the squeezed, displaced, and rotated version of non-Gaussian terms acting only on the vacuum state,

    |ψ(t)⟩=U_a(t)U_b(t)|ψ̃(t)⟩,
    
where recall Eq. (<ref>). An attractive feature of this approach is that we can work in the Gaussian limit by setting 
|ψ̃(t)⟩→ |vac⟩, and still study the effects of pump depletion and the dynamics of the joint spectral amplitude in squeezed state generation. 


That is what we do for 
the remainder of this
paper. We put 

    |ψ(t)⟩→ U_a(t) U_b(t) |vac⟩≡ |ψ_G(t)⟩,

and we have

    |ψ_G(t)⟩ =   e^iθ_a(t)S_a( J(t))|vac⟩⊗ e^iθ_b(t)|β(t)⟩.

where  |β(t)⟩=D_b(β(t)) |vac⟩ is a coherent state.
The Gaussian state is a product of a squeezed vacuum state written in terms of the output operators and a coherent state written in terms of the input operators, and from (<ref>) its dynamics is governed by 

    iħd |ψ_G(t)⟩/dt = (H_a(t)+H_b(t)) |ψ_G(t)⟩.

This of course is an approximate result, and there are consequences. For example,
the operator Q =  a^†^ T a /2+  b^†^ T b commutes with H_L and H_NL, and therefore with the full Hamiltonian in Eq. (<ref>), and so it is a conserved quantity. 
It therefore also commutes with H_NL(t), which governs the evolution of the full |ψ(t)⟩ according to Eq. (<ref>). 
However, it does not commute with the Hamiltonian H_a(t)+H_b(t) that governs the evolution of the its approximation |ψ_G(t)⟩. Nonetheless, in Appendix <ref> we show that ⟨ψ_G(t)|Q|ψ_G(t) ⟩ is constant as |ψ_G(t)⟩ evolves according to Eq. (<ref>).







§ FIRST-ORDER SOLUTION


If the energy of the pump is sufficiently low or the nonlinear coefficients Λ_μν l(t) are sufficiently weak, then the elements of the matrix ζ(t) in Eq. (<ref>) will be much less than unity, 
[ζ(t)]_μν≪ 1 for all μ and ν. In this case the solution to Eq. (<ref>) and Eq. (<ref>) can be approximately written as <cit.>

    V(t)    ≈ I,
        
     W(t)    ≈  -2i∫_-∞^t dt' ζ(t').

Consequently,  the pump depletion rate in Eq. (<ref>) is negligible

    d[β(t)]_l/dt   ≈ -2i Λ^*_μν l(t)∫_-∞^t dt'[ζ(t')]_μν≈ 0,

since it depends on the product of the 
nonlinear coefficients, which are assumed to be small. The Hamiltonian  H_b(t) in Eq. (<ref>)  can be set to zero, 
H_b(t) → 0, since the coefficients in H_b(t) are given by d[β(t)]_l/dt. Similarly the phase θ_b(t) in  Eq. (<ref>)  is negligible, θ_b(t) ≈ 0.


With these approximations, the Gaussian ket |ψ_G(t)⟩ in Eq. (<ref>) can be written as

    |ψ_G(t)⟩ = S_a( J(t))|vac⟩⊗|β_in⟩,

which satisfies the Schrödinger equation

    iħd |ψ_G(t)⟩/dt = H_a(t) |ψ_G(t)⟩.

Treating H_a(t) as a small quantity we can construct a first-order solution to  Eq. (<ref>) using the Magnus expansion <cit.>, and keeping only the first term we obtain

    |ψ_G(t)⟩≈exp(-i∫_-∞^tdt' [ζ(t')]_μν a^†_μ a^†_ν  -  H.c.) ⊗|β_in⟩.

The higher-order terms in the expansion involve the commutator between the Hamiltonian at two different times, such as [H_a(t), H_a(t')] <cit.>. The ket in Eq. (<ref>) is the solution when we assume that [H_a(t), H_a(t')] = 0 so that time-ordering effects are negligible. It was pointed out earlier for waveguides <cit.> that this assumption becomes invalid when the power of the pump is increased. We  demonstrate 
below that this also holds for a squeezed state generated in a ring resonator. Comparing Eq. (<ref>) and Eq. (<ref>) we immediately obtain the approximate first-order solution for the squeezing matrix

    J(t) ≈ -2i∫_-∞^t dt' ζ(t').

This is the expression for W(t) in Eq. (<ref>), so we can write J(t) ≈ W(t). To obtain the form in Eq. (<ref>) we have put θ_a(t) = 0, which is valid for a weak nonlinearity (see Appendix <ref>).




















§ EXAMPLE: SQUEEZED STATE GENERATION IN A RING RESONATOR


In this section we study the example of generating a highly squeezed state in a ring resonator. We calculate the squeezing matrix, number of generated photons, Schmidt number, and the second-order correlation function in the actual output channel including scattering loss, and compare these results to the first-order solution.
 



We consider a SiN ring resonator point-coupled to an actual waveguide and a
phantom waveguide, the latter modelling scattering loss; see
Fig. <ref>(b).  For a sample calculation we consider a structure with
a ring radius of R = 64 μm, ω_S = 2π× 193  THz, a group velocity of 1.5× 10^8  m/s, a nonlinear parameter of γ_nl≈ 1( Wm)^-1 (see Eq. (<ref>)), and an intrinsic quality factor of 2× 10^6 for all three ring resonances. We take the loaded quality factors for the P, S, and C modes to be 4×10^4 (over-coupled), 2× 10^6 (over-coupled), and 1× 10^6 (critically-coupled), respectively. 

In mode C there is an undepleted CW pump with power P_C = 30 mW, and in mode P there is a weak pulsed pump with energy U_P and duration τ_P = 0.5 ns. The initial amplitude of the pulse is

    [β_in]_l = √(N_P v_Pτ_P/√(π))exp(-1/2(k_l-K_P)^2v_P^2τ_P^2),

where for simplicity we neglect group velocity dispersion, N_P is the initial photon number, v_p is the group velocity, and K_P is the central wavenumber. 







In Fig. <ref>(a) the magnitude of the squeezing matrix is shown for U_P = 100 pJ (N_P≈7.8× 10^8). The x-axis and y-axis represent the detuning from the center of the ring resonance of the mode S, where v_S, K_S, and 1/Γ_S≈ 330  ps  are its group velocity, central wavenumber, and dwelling time, respectively. For this energy the average generated photon number in the squeezed state is about 19, which corresponds to approximately 19  dB of squeezing, putting it in the high-gain regime.  

The squeezing matrix in Fig. <ref>(a) is calculated using the approximate first-order solution in Eq. (<ref>). In order to compare the first-order solution to the full Gaussian solution, in Fig. <ref>(b) the magnitude of the squeezing matrix is shown for the same parameters as above but 
calculated using the full Gaussian solution in Eq. (<ref>) that is valid for arbitrary energy. The full Gaussian solution in Fig. <ref>(b) exhibits less correlation between the wavenumbers (k_1 and k_2) of the generated photon pairs than the first-order solution in Fig. <ref>(a). 











In Fig. <ref>(a) the average generated photon number in the squeezed state is shown for pump energy (U_P) between 1  pJ and 100  pJ. The calculation is done with the full Gaussian solution (blue line) and the first-order solution (pink dashed). The photon number is calculated with tr( W^*(t)  W^ T(t)), where the matrix W(t) is obtained using the full Gaussian solution or the first-order solution. The first-order solution gives a photon number that is approximately 20% less than the full Gaussian solution for pump energies at or below 60  pJ. This corresponds to an underestimation of the squeezing by about 1  dB at most. But at 80  pJ and 100  pJ the first-order solution underestimates the  squeezing by 1.6  dB and 2.3  dB, respectively. Despite the high level of squeezing, the level of pump depletion is negligible. As a result the neglect of the non-Gaussian corrections to the vacuum, which are captured by the evolution Eq. (<ref>) of |ψ̃(t) ⟩ – and which we expect to only become important when there is significant pump depletion <cit.> – is justified. 

In Fig. <ref>(b) the Schmidt number <cit.> of the squeezing matrix is shown for the same parameters as above. The first-order solution consistently overestimates the Schmidt number, and thus
the squeezing matrix is more separable than what is predicted by the first-order solution.   













Each Schmidt mode typically involves light exiting into both the actual and the phantom waveguide, the latter modelling the scattering of light off the chip. To focus on the properties of the light exiting the actual waveguide, we
consider a measurement  of the second-order correlation function G^(2)(t_1,t_2) of the squeezed light exiting the actual output channel. It is convenient to introduce a channel operator a̅(x,t) as the Fourier transform of the output operators for the actual waveguide a_k as <cit.>

    a̅(x,t) = ∫dk/√(2π) a_k e^i(k-K_S)(x-v_St),

where x is the direction of propagation in the waveguide. The operator a̅(x,t) is associated with an electric field at the position x and time t.  Just beyond the coupling point between the ring and waveguide (x=0^+), the second-order correlation function is given by <cit.>

    G^(2)(t_1,t_2) = v_S⟨a̅^†(0^+,t_1)a̅^†(0^+,t_2) a̅(0^+,t_2)a̅(0^+,t_1)⟩,

where the expectation value is done with the Gaussian ket in Eq. (<ref>). This can be used to predict the probability of detecting coincidence counts at the two times t_1 and t_2. We normalize it using Φ^2, which is defined as the ratio of the number of generated photons per pump pulse duration, given by 

    Φ≡⟨ n_S⟩/τ_P,

where ⟨ n_S⟩ is the number of generated photons.  

In Figs. <ref>(a) and <ref>(b) the second-order correlation function for the actual output channel is shown using the first-order solution, as well as in Figs. <ref>(c) and <ref>(d) using the full Gaussian solution. When the energy of the pump is 10  pJ the correlation functions from the first-order solution and full Gaussian solution appear to agree. For the high energy pump at 100  pJ the first-order solution begins to show correlations that are not present in the full Gaussian solution, and the maximum of the first-order solution occurs at an early time (t_1 = t_2 = 0.9/Γ_S) than the full Gaussian solution (t_1 = t_1 = 1.5/Γ_S). 









§ CONCLUSION


In conclusion, we presented a multimode theory of Gaussian state generation for an effective χ^(2) interaction, in terms of asymptotic-in and -out fields, that includes scattering loss and pump depletion. We showed that the full ket for our system can be written as a Gaussian unitary acting on non-Gaussian
ket that satisfies a Schrödinger equation with an effective Hamiltonian that describes the quantum correlations between the pump and the generated light. In order for this to be valid, we required that the squeezing, rotation, and displacement parameters of the Gaussian unitary were solutions to a set of coupled differential equations.  As an example we studied the lossy generation of a highly squeezed state using an effective χ^(2) interaction in a silicon nitride ring resonator, and we compared this full Gaussian solution to a lowest-order Magnus expansion for the ket. In future work, we will study the evolution of the  non-Gaussian  ket, the effects of self- and cross-phase modulation, and the application to more general resonant systems. This work opens up a path to study the deterministic generation of non-Gaussian states in resonant systems. 




§ ACKNOWLEDGEMENTS

We thank Marco Liscidini for valuable discussions. J. E. S. and C. V. acknowledge the Horizon-Europe research and innovation program under grant agreement ID: 101070168 (HYPERSPACE) and Natural  Sciences  and  Engineering  Research Council of Canada (NSERC) for funding.



 










§ EFFECTIVE Χ^(2) INTERACTION IN A RING RESONATOR


In this Appendix we derive the nonlinear Hamiltonian for an effective χ^(2) interaction in a ring resonator. 

We begin with the nonlinear Hamiltonian for SFWM in Eq. (<ref>). In the interaction picture it takes the form

    H_NL(t)     = -3/ϵ_0∑_n,n'∫ dk_1dk_2dk_3dk_4 
       × K^nn'(k_1,k_2,k_3,k_4)  e^-iΩ_nn'(k_1,k_2,k_3,k_4)t
       × a^†_n S k_1a^†_n' S k_2 b_P k_3b_C k_4 +  H.c.,

where we have used the definition in Eq. (<ref>), and 
defined

    Ω_nn'(k_1,k_2,k_3,k_4) ≡ω_Ck_4+ ω_Pk_3 - ω_n'Sk_2-ω_nSk_1.

In order to build the effective χ^(2) interaction, the pump in the mode C is treated classically with a cw field that is  undepleted. 
The  operators b_Ck_4 are replaced with the amplitudes <cit.>

    b_Ck_4→√(2π P_C/ħω_C v_C)δ(k_4 - K_C),

where P_C is the cw power and v_C is the group velocity of the cw pump in the actual input channel.
Under this approximation, Eq. (<ref>)  becomes

    H_NL(t)     = ħ∑_n,n'∫ dk_1 dk_2 dk_3Λ^nn'(k_1,k_2,k_3,t)  
       × a^†_n S k_1a^†_n' S k_2b_Pk_3+  H.c.,

 where the effective χ^(2)  coefficient  Λ^nn'(k_1,k_2,k_3,t) is given by

    Λ^nn'(k_1,k_2,k_3,t) =    -3/ħϵ_0 √(2π P_C/ħω_C v_C)K^nn'(k_1,k_2,k_3,K_C)
       × e^-iΩ_nn'(k_1,k_2,k_3,K_C)t .


We consider a ring resonator that is point-coupled to an actual waveguide that couples the light out of the ring and a phantom waveguide that handles the scattering loss; recall Fig. <ref>(b). The circumference of the ring is given by ℒ = 2π R, where R is its radius. The resonant wavenumbers for the ring are given by κ, such that κℒ  =2π m, where m is a positive integer. We assume that around each ring resonance we can neglect group velocity dispersion, and that the dispersion relation is the same for the actual and phantom waveguides. Under these assumptions the dispersion relation for either channel is given by

    ω_J(k) = ω_J + v_J(k - K_J),

where ω_J is the center frequency of the resonance J, v_J is the group velocity in either channel, and K_J is the wavenumber for the light in either channel with frequency ω_J. The linewidth of the ring resonances are given by 2Γ_J, which includes coupling losses to the actual waveguide and scattering losses to the phantom waveguide; the loaded quality factor is then Q^load_J = ω_J/(2Γ_J).

The asymptotic -in and -out displacement fields are defined over all space, but we assume that the nonlinear interaction takes place mostly inside the ring resonator.  Thus to obtain the nonlinear coefficient we only require expressions for the field that hold at positions in the ring. To simplify the integration of that region of space, we introduce the coordinates of the ring-frame as ( r_⊥, ζ), where r_⊥ = (ρ, z) and ζ = ϕ R, and (ρ, ϕ, z) are cylindrical coordinates. 
For positions inside the ring, the asymptotic-in pump fields can be written as <cit.>

    D^in_P  k( r)    = -√(ħω_P/4π) d_P( r_⊥;ζ) e^i ζκ_PF_P-(k) ,
    
     D^in_C  k( r)    = -√(ħω_C/4π) d_C( r_⊥;ζ) e^i ζκ_CF_C-(k) ,

and the asymptotic-out fields can be written as

    D^out_n S k( r)    = -√(ħω_S/4π) d_S( r_⊥;ζ) e^i ζκ_SF^(n)_S+(k).

Here we introduced the field enhancement factors for the actual input channel, F_J -(k), given by <cit.>

    F_J -(k) =1/√(ℒ)( (γ_J)^*/v_J(K_J - k)- iΓ_J)  ,

where J=P,C and γ_J is a coupling constant between a discrete ring mode and the continuous waveguide mode of the actual input channel at the frequency band J. The field enhancement factors for the output channels, F^(n)_S +(k), are given by

    F^(n)_S+(k) =1/√(ℒ)( (γ^(n)_S)^*/v^(n)_S(K^(n)_S - k)+ iΓ_S)  ,

where now the coupling constant, group velocity, and center wavenumber depend on the output channel n. The coupling constant for a given channel is related to the intensity decay rate into that channel, as introduced earlier <cit.>. For simplicity we assume that the coupling constants are real. 

Putting Eqs. (<ref>), (<ref>), and (<ref>) into Eq. (<ref>), the nonlinear coefficient for SFWM can be written as
    
    K^nn'(k_1,k_2,k_3,k_4)    =ħ^2 ϵ_0 v_Pv_C/12π^2γ_nlω_Sℒ F^(n)*_S+(k_1)  
       × F^(n')*_S+(k_2) F_P-(k_3)F_C-(k_4),

where we have introduced the nonlinear parameter <cit.>

    γ_nl   =  3√(ω_Pω_C)/4 ϵ_0 v_Pv_Cℒ∫_ ring d r_⊥ dζΓ^ijkl_(3)( r)
       × d^*i_S( r_⊥;ζ) d^*j_S( r_⊥;ζ)  d^k_P( r_⊥;ζ)  d^l_C( r_⊥;ζ)  e^i Δκζ,

which has units of ( Wm)^-1, and  Δ k = κ_P + κ_C - 2κ_S is the phase matching function. 

Now to obtain an effective second-order nonlinear coefficient in the ring, we put Eq. (<ref>) into Eq. (<ref>), and find

    Λ^nn'(k_1,k_2,k_3,t)   =-ħ v_Pv_Cγ_nlω_S ℒ/4π^2 √(2π P_C/ħω_C v_C)
       × F^(n)*_S+(k_1) F^(n')*_S+(k_2)F_P-(k_3)  
       ×F_C-(K_C) e^-iΩ_nn'(k_1,k_2,k_3,K_C)t.

Then putting Eq. (<ref>) into Eq. (<ref>) we obtain the nonlinear Hamiltonian for an effective χ^(2) interaction in a ring resonator.  In order to work with this Hamiltonian numerically the wavenumbers in all channels have to be discretized.   Each waveguide mode has a discrete label j, such that its wavenumber is given by k_j = 2π m_j/L, where m_j is an integer mode number, and L is the quantization length <cit.>. Using this discrete notation, the nonlinear Hamiltonian in Eq. (<ref>)  can be written as

    H_NL(t)   =  ħ∑_n,n'∑_i,j,lΛ^nn'_ijl(t)a^†_nSia^†_n'Sj b_Pl +  H.c.,

where we defined

    Λ^nn'_ijl(t) ≡(2π/L)^3/2Λ^nn'(k_1i,k_2j,k_3l,t),

where Λ^nn'(k_1i,k_2j,k_3l,t) is the effective second-order nonlinear coefficient for the ring given in Eq. (<ref>),
and the operators satisfy the commutation relations

    [a_n J i, a^†_n J j]  = [b_P i, b^†_ P j] = δ_ij.

The Hamiltonian in Eq. (<ref>)  can be written equivalently as

    H_NL(t)   =  ħ∑_(i,n),(j,n')∑_lΛ_(i,n)(j,n')l(t)
       × a^†_S(i,n)a^†_S(j,n') b_Pl +  H.c..

We then introduce the discrete indices μ and ν that are constructed by grouping the labels i and n and j and n' together as

    μ   ≡ (i,n),
        
    ν   ≡ (j,n'),

where  μ refers to a photon in the output channel n with wavenumber k_i. With this labelling we obtain the Hamiltonian in Eq. (<ref>).












§ DERIVING EQUATIONS FOR THE GAUSSIAN PARAMETERS


In this Appendix we show how the Eq. (<ref>) for U_a(t) leads to the Eqs. (<ref>) and (<ref>) for V(t) and W(t), and how the Eq. (<ref>) for U_b(t) leads to the Eqs. (<ref>) and (<ref>) for β(t) and θ_b(t).

Beginning with Eq. (<ref>) for U_a(t), it was pointed out earlier
<cit.> that, with an H_a(t) of the form (<ref>),     its unique solution can indeed be written as 
(<ref>). The task here is just to identify how the parameters in U_a(t) depend on time for a specified H_a(t). The strategy that we employ relies on introducing time-dependent output operators a(t) defined as

    a(t)    = U^†_a(t)  a U_a(t) 
       =  V(t)  a +  W(t)  a^†,

where V(t) and W(t) are found to be given by

    V(t)    ≡cosh( u(t)) e^iϕ(t) ,
    
     W(t)   ≡sinh( u(t)) e^iα(t) e^-iϕ^ T(t),

and where the Hermitian matrices u(t) and α(t) that appear here come from the polar decomposition of the squeezing matrix, 

    J(t) =  u(t) exp(iα(t)).
 
The matrices V(t) and W(t)  in Eqs. (<ref>) and (<ref>) satisfy the constraints

    V(t)  V^†(t) -  W(t) W^†(t)    = I,
    
     V(t)  W^ T(t) -  W(t)  V^ T(t)   =0,

which are necessary so that the operators in Eq. (<ref>) satisfy the usual equal time bosonic commutation relations.
Taking the time-derivative of the first line in Eq. (<ref>) we obtain the  equation of motion

    d  a(t)/dt = i/ħU_a^†(t)[H_a(t) ,  a]U_a(t),

where we used Eq. (<ref>). Now using  Eqs. (<ref>) and (<ref>), Eq. (<ref>) becomes

    d  a(t)/dt   =-2iζ(t) W^*(t)  a - 2iζ(t) V^*(t)  a^†,

where  ζ(t) is given by Eq. (<ref>). On the other hand, taking the time derivative of the second line in Eq. (<ref>), we have equivalently

    d  a(t)/dt   =d V(t)/dt a+d W(t)/dt a^†.

Since the expressions for d a(t)/dt in Eq. (<ref>) and Eq. (<ref>) are equivalent, we can equate the elements of the matrix coefficients multiplying the a and a^† operators in each expression, leading to the coupled equations Eqs. (<ref>) and (<ref>). 



 
 

Next we turn to deriving Eq. (<ref>) for β(t) and Eq. (<ref>) for θ_b(t). We start by multiplying  Eq. (<ref>) by U^†_b(t) from the left

    iħ U^†_bdU_b(t)/dt   = ħ [γ(t)]_l b^†_l +  H.c.,

where γ(t) is given by Eq. (<ref>). Putting Eq. (<ref>) into the left hand side of Eq. (<ref>) we obtain

    iħ D^†_bdD_b(t)/dt - ħdθ_b(t)/dt   = ħ [γ(t)]_l b^†_l +  H.c..

The derivative of the displacement operator is obtained by disentangling it and using the chain rule. Doing this we obtain

    iħ D^†_bdD_b(t)/dt   =iħd[β(t)]_l/dtb^†_l    +i/2ħd[β(t)]_l/dt[β^*(t)]_l+ H.c..

Putting Eq. (<ref>) into Eq. (<ref>) and collecting terms, we obtain

    0   =(id[β(t)]_l/dt  - [γ(t)]_l )b^†_l    - (id[β^*(t)]_l/dt  + [γ^*(t)]_l )b_l 
       +i/2(d[β(t)]_l/dt[β^*(t)]_l - d[β^*(t)]_l/dt[β(t)]_l )- dθ_b(t)/dt,

which leads to
Eqs. (<ref>) and (<ref>). 





















§ EXTRACTING J(T) AND Φ(T) FROM V(T) AND W(T)


We begin by doing a numerical polar decomposition of V(t). Note that Eq. (<ref>) is already in the form of a polar decomposition, where  cosh( u(t)) is a Hermitian positive semi-definite (PSD) matrix and  e^iϕ(t) is a unitary matrix. And since the Hermitian PSD matrix in a polar decomposition is always unique, the diagonal elements of cosh( u(t)), which we identify by the set {s_i(t) }, will be unique.  Further, since the matrices cosh( u(t)) and u(t) commute they can be diagonalized by the same unitary transformation. So we can write 

    cosh ( u(t)) =  F(t)  diag(s_1(t), …)  F^†(t).

and

    u(t) = F(t) 
        diag(r_1(t),…)  F^†(t),

where the unitary matrix F(t) is the same in both decompositions, and the set { r_i(t) } are the diagonal elements of u(t). Since Eqs. (<ref>) and (<ref>) are diagonalizations of PSD Hermitian matrices, we have 

    s_i(t) ≥ 0

and

    r_i(t) ≥ 0.

But a condition stronger than Eq. (<ref>) holds. Writing
Eq. (<ref>)  as 

    u(t) =  F(t)  r(t)  F^†(t),

where r(t) is a diagonal matrix we have 

    cosh( u(t))    = cosh( F(t)  r(t)  F^†(t)) 
       =  F(t) cosh( r(t))  F^†(t),

or 

    cosh ( u(t)) =  F(t)  diag(cosh r_1(t),…)  F^†(t).

Comparing with Eq. (<ref>) we see that 

    s_i(t) = cosh(r_i(t)),

and so, given Eq. (<ref>), we have 

    s_i(t) ≥ 1.


A first consequence of this, which follows from using it in Eq. (<ref>),
is that cosh ( u(t)) is 
invertible and then so is V(t); thus from the numerical
polar decomposition of V(t)
we not only extract a unique Hermitian matrix cosh ( u(t)) but also a
unique unitary matrix e^i ϕ(t). A second consequence is that from Eq. (<ref>) we can write Eq. (<ref>) as

    u(t)=  F(t)  diag
        (  arccosh (s_1(t)), …)  F^†(t)

where by “arccosh" we mean the non-negative arccosh, by virtue of Eq. (<ref>). 
In summary, then, from the numerical polar decomposition of V(t) we immediately obtain the unique matrices e^i ϕ(t) and cosh ( u(t)), from the diagonalization Eq. (<ref>) we determine the unitary matrix F(t) and the { s_i(t) }, and then from Eq. (<ref>) we construct u(t).

With u(t) in hand we can then form 

    sinh ( u(t))    = sinh ( F(t)  r(t)  F^†(t)) 
       =  F(t) sinh ( r(t))  F^†(t),

or equivalently 

    sinh( u(t)) =  F(t)  diag(sinh(r_1(t)),…)  F^†(t),

which can be written as

    sinh ( u(t))    =  F(t)  diag(sinh (r_1(t))/r_1(t),…) 
       × diag(r_1(t),…)  F^†(t),

or

    sinh ( u(t)) =  K(t)  u(t),

where we have used Eq. (<ref>) and let

    K(t) =  F(t)  diag(sinh (r_1(t))/r_1(t),…)  F^†(t),

with 

    K_i(t) = sinh (r_i(t))/r_i(t)

Forming K_i(t) seems problematic if r_i(t)=0. Of course, in such an instance the natural choice would be to take K_i(t) to be the limit of Eq. (<ref>) as r_i(t) → 0, which would give K_i(t) = 1. But one could worry about issues of uniqueness.  In fact, we show below that if r_i(t)=0 we can take K_i(t) to be any nonzero number, and the final result of the following calculation of J(t) will be the same.


The matrix K(t) is clearly invertible, since we can identify 

    K^-1(t) =  F(t)  diag(1/K_1(t),…) F^†(t),

and furthermore we can numerically construct K^-1(t) since we know F(t) and the r_i(t). From Eq. (<ref>) we then have 

    K^-1(t) sinh( u(t)) =  u(t),

and from Eq. (<ref>)
we have

    K^-1(t)  W(t)    =   u(t)  e^iα(t) e^-iϕ^ T(t) =  J(t)  e^-iϕ^ T(t),

so

    J(t) =  K^-1(t)  W(t)  e^iϕ^ T(t).

All the quantities on the right-hand-side are known, and so we can determine J(t). 


To investigate what happens to J(t) when some of the eigenvalues r_i(t) are zero, consider

    K^-1(t)  W(t)     =  F(t)  diag(sinh (r_1(t))/K_1(t) ,…) 
       × F^†(t) e^iα(t) e^-iϕ^ T(t),

where we used Eqs. (<ref>) and (<ref>).
Now we can see that for any i for which r_i(t)=0 we will have sinh(r_i(t))/K_i(t)=0, regardless of which nonzero value K_i(t) is set. 
So the eigenvalues r_i(t) that are zero do not contribute to the matrix J(t).

At this point we have identified J(t) and e^iϕ(t); we still must extract ϕ(t). The ϕ(t) appearing in R_a(ϕ(t)) of (<ref>) is a Hermitian matrix, and so can be diagonalized,

    ϕ(t)=  M(t)  diag (χ_1(t),…)  M^†(t),

where M(t) is a unitary matrix. Thus we have 

    e^i ϕ(t)   = e^i  M(t)  diag (χ_1(t),…)  M^†(t) 
       =  M(t) e^i  diag (χ_1(t),…) M^†(t).

And so numerically diagonalizing the identified e^i ϕ(t) would naturally lead us to identify ϕ(t) as given by Eq. (<ref>). But although the eigenvalues of e^i ϕ(t) are well-defined, the values of the phases appearing in those eigenvalues are not. For we could just as easily write 
    e^i ϕ(t)   =  M(t) e^i  diag (χ_1(t)+2π m_1(t),…) M^†(t) 
       =  e^i  M(t)  diag (χ_1(t)+2 π m_1(t),…)  M^†(t) ,

where m_1(t), etc., are any integers.  So instead of Eq. (<ref>) we could just as well claim that ϕ(t)
is given by 

    ϕ(t)   =  M(t)  diag (χ_1(t)+ 2 π m_1(t),...)  M^†(t).

However, this ambiguity is benign, because both Eq. (<ref>) and Eq. (<ref>) lead to the same effective R_a(ϕ(t)).

To see this, note that if we posit Eq. (<ref>) we have 

    [ϕ(t)]_kl=  [ M(t)]_kp(χ_p(t)+2 π m_p(t)) [ M^†(t)]_pl,

and so 

    R_a(ϕ(t)) =  e^i  a^†_k [ M(t)]_kp (χ_p(t) +2 π m_p(t)) [ M^†(t)]_pl a_l.

Now put 

    c_p(t) = [ M^†(t)]_pl a_l = [ M^*(t)]_lp a_l,

so

    c^†_q(t)= [ M(t)]_jq a^†_j.

Thus

    [c_p(t),c^†_q(t)] = δ_pq,

and the new operators c_p(t) and c^†_q(t) satisfy equal time commutation relations of raising and lowering harmonic oscillator operators.  In terms of them, we can write (<ref>) as

    R_a(ϕ(t))    =  e^i  c^†_p(t)(χ_p(t)+2 π m_p(t)) c_p(t)
       =  e^i c^†_p(t) χ_p(t) c_p(t) e^ 2 π i  m_p(t) c^†_p(t) c_p(t).

As basis kets we can choose the eigenkets | {n},t ⟩ of all the number operators c^†_q(t)c_q(t),
where 

    c^†_q(t) c_q(t) | { n },t ⟩  = n_q| {n },t ⟩.
 Then

    e^ 2 π i m_p(t) c^†_p(t) c_p(t) | { n } ,t ⟩ =  e^2 π i  m_p(t) n_p  | { n } ,t ⟩ = | { n } ,t ⟩,

and since any ket | φ(t) ⟩ can be written as a superposition of the | { n } ,t ⟩, from (<ref>) and (<ref>) we have 

    R_a(ϕ(t)) | φ(t) ⟩   =  e^i  c^†_p(t)(χ_p(t)+2 π m_p(t)) c_p(t) | φ(t) ⟩
       =  e^i c^†_p(t) χ_p(t) c_p(t) | φ(t) ⟩.

Thus the result we find for R_a(ϕ(t)) will be the same whether we use  Eq. (<ref>) or Eq. (<ref>) in identifying the logarithm of e^i ϕ(t) to construct R_a(ϕ(t)); any logarithm of e^i ϕ(t) will do.











§ DETERMINING THE PHASE Θ_A(T)


The phase θ_a(t) appears in the unitary operator U_a(t) in Eq. (<ref>). In this section we determine this phase following the approach of Ma and Rhodes <cit.>. Recall that U_a(t) is the solution to the Schrödinger equation in Eq. (<ref>) with the Hamiltonian H_a(t) defined in Eq. (<ref>). Given the form of H_a(t) in Eq. (<ref>), it follows <cit.> that the phase θ_a(t) is given by

    θ_a(t) = -1/2∫_-∞^t dt'  tr(ζ^*(t') tanh( u(t')) e^iα(t')) + c.c.,

where ζ(t) is defined in Eq. (<ref>) and is related to the coefficients in H_a(t); recall Eq. (<ref>) for the polar decomposition of the squeezing matrix J(t) =  u(t) exp(iα(t)).

The first-order solution that is discussed in Sec. <ref> is valid when the amount of squeezing is small or the nonlinearity is weak.  In this subsection we show that the phase θ_a(t) in Eq. (<ref>) is approximately zero in the first-order solution. In this case because the entries of the squeezing matrix J(t) are small  we replace tanh( u(t'))exp(iα(t')) in Eq. (<ref>) with u(t')exp(iα(t')) =  J(t') (recall Eq. (<ref>)). Doing this we can write Eq. (<ref>) approximately as

    θ_a(t) ≈ -1/2∫_-∞^t dt'  tr(ζ^*(t')  J(t')) + c.c..

But from Eq. (<ref>) we have that J(t') = -2i∫_-∞^t' dt”ζ(t”). Putting this into Eq. (<ref>), the term inside the trace is ζ^*(t')ζ(t”), which is much smaller than unity. Thus θ_a(t) is approximately zero. 











§ A CONSERVED QUANTITY


We define the operator Q̂ 

    Q̂≡1/2 a^†^ T a +  b^†^ T b,

which is a conserved quantity with respect to the full Hamiltonian in Eq. (<ref>), that is 

    [Q̂,H_NL(t)] = 0.

But Q̂ is not a conserved quantity with respect to the effective Hamiltonian in Eq. (<ref>). In this Appendix we show that the average value of Q̂ in the Gaussian limit is nonetheless a conserved quantity. 


We begin with expectation value of Q̂ using the ket in the Gaussian limit given by Eq. (<ref>) 

    ⟨Q̂⟩   ≡⟨ψ_G(t)|Q̂|ψ_G(t)⟩
    
    
    
    
    
    
    
    
    
    
    
    
    
       = 1/2 tr( W^†(t)  W(t))   + [β^*(t)]^ Tβ(t).
 
In what follows we drop the time-dependence of W(t) and β(t) for convenience. This quantity is conserved if its time derivative is zero

    d ⟨Q̂⟩/dt = 0.

Putting Eq. (<ref>) into Eq. (<ref>) we obtain

    1/2 trA + B=0

where we have defined

    A    ≡d  W^†/dt W+ W^†d W/dt,
    
    
    B   ≡∑_l( dβ^*_l/dtβ_l + β^*_ldβ_l/dt).

Now we use the differential equations for β_l and W in Eqs. (<ref>) and (<ref>) to simplify A and B. First, to simplify  B we put Eq. (<ref>) into Eq. (<ref>)

    B   = i  ∑_i,k,jγ_kjiβ_i [ V^*  W^†]_kj -  i∑_i,k,jγ^*_kjiβ^*_i [ V  W^ T]_kj
       =i∑_k,j[ζ]_kj[ V^*  W^†]_kj  -i∑_k,j[ζ^*]_kj[ V W^ T]_kj
       =i tr(ζ W^*  V^†)- i tr(ζ^*  W  V^ T).
 
Next we simplify A by putting Eq. (<ref>) into Eq. (<ref>), and
obtain

    A    =2i V^ Tζ^* W-2i W^†ζ V^*.

Putting Eqs. (<ref>) and (<ref>) into Eq. (<ref>)

    0   =1/2 tr[  2i V^ Tζ^*W-2i W^†ζ V^*  ]  +i tr(ζ W^*  V^†)
       - i tr(ζ^*  W  V^ T) 
       = i tr(ζ^* W  V^ T)-i tr(ζ V^* W^†)+i tr(ζ W^*  V^†)
       - i tr(ζ^*  W  V^ T),

and since W^* V^† =  V^* W^† the right hand side is zero.








