Nonlinear cross-Kerr quasiclassical dynamics

We study the quasiclassical dynamics of the cross-Kerr effect. In this approximation, the typical periodical revivals of the decorrelation between the two polarization modes disappear and they remain entangled. By mapping the dynamics onto the Poincare space, we find simple conditions for polarization squeezing. When dissipation is taken into account, the shape of the states in such a space is not considerably modified, but their size is reduced.

Special mention must be made of the role that this cubic nonlinearity has played in the generation of squeezed light. The first proposals employed schemes involving a nonlinear interferometer [45] or degenerate four-wave mixing [46,47]. But quite soon optical fibers became the paradigm for that purpose [48][49][50][51][52][53]. However, due to the typically small values of the nonlinearity in silica glass [54], Kerr-based fibers need long propagation distances and high powers, which bring other unwanted effects [55,56].
In this paper, we direct out attention to this limit of high intensities, in which one could expect a classical description to be pertinent. Under reasonable assumptions, Maxwell's equations lead to a set of coupled nonlinear Schrödinger equations that has long been a useful tool for depicting the behavior of optical fields in nonlinear dispersive media. It has proved valuable in the description of such diverse phenomena as pulse compression, dark soliton formation, and self-focusing of ultrashort pulses [57]. Yet there remain nonclassical features that cannot be explained in this classical manner. To put it differently, at the most basic level, the propagation of light in a Kerr medium is accompanied by unavoidable quantum effects.
The considerations thus far indicate that the regime we wish to explore can be regarded as a problem at the boundary between classical and quantum worlds. Probably, the transition between both can be best scrutinized by exploiting phase-space methods [58][59][60]. This opens up the possibility of gaining some information about the nonclassical behavior with a quasiclassical description that employs essentially classical trajectories, while the correct quantum initial state is taken into account via, e.g., the Wigner function [61][62][63]. Despite some problems with the interpretation, the Wigner function has enjoyed substantial attention in various domains of physics [64] and has already been applied to some nonlinear problems in quantum optics [65][66][67][68][69].
The intensity dependence of the refractive index, which is the hallmark of the Kerr effect, can manifest itself in two different ways: as a self-phase modulation and as a cross-phase modulation. Self-phase modulation refers to the self-induced phase shift experienced by an optical field during its propagation, whereas cross-phase modulation refers to the nonlinear phase shift of an optical field induced by another one having a different wavelength, direction, or state of polarization.
In this paper we focus on the cross-Kerr effect, for it is especially germane to attain polarization squeezing, a major goal in our laboratory [70]. We capitalize on the quasiclassical approach to re-analyze the light propagation in this case in a very concise way: after neglecting higher-order fluctuations, we get an evolution equation for the Wigner function that can be integrated to an analytical form. This allows us to study the dynamics of mode correlations. Since the resulting state in non-Gaussian, the application of common entanglement criteria [71,72] becomes problematic, so we content ourselves with the study of the purity of the reduced states, which can be carried out in a closed form.
The two-mode Wigner function is appropriately cast in Poincaré space in terms of the phase-space version of the Stokes parameters, which affords an intuitive picture. Finally, as the Kerr dynamics is photon-number preserving, the standard models of dissipation [73] based in coupling the modes to lossy reservoirs, seem inadequate. Instead, we allow for dissipation through pure dephasing processes which turns out to be exactly solvable. The resulting evolution reveals that the shape of the Wigner functions is not considerably modified, but their size is shrunk.

Cross-Kerr quasiclassical evolution
As heralded in the Introduction, the cross-Kerr configuration corresponds to a situation in which the refractive index of a beam (say a) is modified by the intensity of a second one (say b). These beams are excited in two orthogonal polarization modes that, in a quantum description, are characterized by two complex amplitude operators, denotedâ andb, respectively. These operators obey the standard bosonic commutation relations 1) the superscript † standing for the adjoint.
In terms of these annihilation and creation operators, the Hamiltonian accounting for the cross-Kerr interaction is [74] H =hχâ †âb †b , (2.2) where χ is an effective coupling constant that depends on the third-order nonlinear susceptibility. For any state described by the density operatorρ, the evolution is given by 3) whose solution can be formally written aŝ By expanding this equation in the two-mode Fock basis |n a , n b , the evolution may be, in principle, tracked. Take the example of an initially pure, two-mode coherent state |Ψ(0) = |α 0 , β 0 , where henceforth the subscript 0 indicates the value of the corresponding variable at t = 0. The resulting state is The term exp(−iχtn a n b ) arises because of the coupling between the modes and causes that the state cannot be factorized into single-mode states; i.e., it becomes entangled, as we shall examine in the next section. It is apparent that equation (2.5) is of practical use only for few-photon states. Actually, such an exact solution does not facilitate to extract the classical part of the dynamics in a manifest form. To that end, we proceed to decompose the mode operatorsâ andb aŝ that is, a sum of classical amplitudes and quantum noise operators. The average values of the noise operators are assumed to be much smaller than the corresponding coherent amplitudes (|α| 2 , |β | 2 ≫ 1), so we can restrict the analysis to first-order terms in δâ and δb. If we employ the two-mode Wigner function W (α, β ) and the basic techniques outlined in Appendix A, equation (2.4), with this linearization ansatz, can be recast as Two comments are in order here. First, we are ignoring quantum fluctuations, inasmuch we are disregarding higher-order moments of the noise operators; this seems a plausible approximation for highly-excited fields. Second, we underline that the evolution is specified only by classical trajectories, much in the spirit of the quasiclassical approximation.
To shed light on the physics embodied in (2.7), we resort to action-angle variables (I , ϕ) for each mode [75]. In our context, they can be defined as α = I a exp(iϕ a ) , therefore ϕ is the polar angle in phase space, whereas I is related to the mode intensity (see figure 1). With these variables, equation (2.7) can be rewritten in a simple and elegant form As ∂ /∂ ϕ generates rotations in phase space, equation (2.9) reflects that the amplitudes in each mode experience different rotations, with angles proportional to the intensity components of the other mode [76,77]. The result is schematized in figure 1: roughly speaking, the shaded area indicates the region in phase space occupied by the state. For an initial coherent state this area is a circle; the top of the circle corresponds to higher intensity and therefore is more phase shifted than the bottom, resulting in an elliptical noise distribution. Equation (2.9) can be readily solved: (2.10) If again we assume initially the two-mode coherent state |α 0 , β 0 (α 0 = √ I 0a e iϕ 0a , β 0 = √ I 0b e iϕ 0b ) and using (A.5), equation (2.10) reduces to where we have introduced the dimensionless variable τ = χt/2. Observe that at τ = 0 the Wigner function is made of two independent Gaussians, while as time goes by the induced mode correlations lead to a non-Gaussian state.

Mode correlation dynamics
Two-mode Gaussian states constitute the simplest example of a continuous-variable bipartite system, the workhorses of quantum information. Accordingly, the theoretical aspects of these states have been extensively worked out and a variety of quantitative characterizations are available for them [78][79][80][81][82].
The unique feature of these Gaussian states is that they are fully specified (up to local displacements) by the covariance matrix γ, with elements γ i j = Tr[ρ{R i ,R j }/2], where {, } denotes the anticommutator andR = (x a ,p a ,x b ,p b ) is the vector of phase-space operators. This covariance matrix can be jotted down as Here, A and B are the covariance matrices associated to the reduced state of the modes a and b, while C describes the correlation between these modes. The symplectic eigenvalues of γ are with ∆ = detA + det B + 2 detC. These symplectic eigenvalues encode all the essential information and provide powerful, simple ways to express fundamental properties. For example, a Gaussian state is entangled if and only if  .11), as a measure for the entanglement between the two modes. We have taken both modes with the same intensity I 0a = I 0b = 10 6 . Entanglement is proven forν − < 0.5, a region which can be observed only in the inset.
where the smallest symplectic eigenvalueν − of the covariance matrix corresponding to the partially transposed state is obtained from ν − by replacing det C with − det C; i.e., by time reversal in the second system and thus a flip of its canonical momentum.
In figure 2 we have plotted the time evolution ofν − for the state (2.11), showing a rapid increase (in the inset, we observe a fluctuating behavior that is smeared out in a larger scale). The main caveat with this approach is that, as mentioned before, our state rapidly becomes non-Gaussian, and criterion (3.3) gives then only a sufficient condition. Consequently, we can certify entanglement just in the short-time window displayed in the inset. This actually holds for any available criterion [71,72]: if the state is entangled, a given test may or not detect its entanglement; in turn, if a particular test does not detect entanglement, we can not conclude separability of the state.
Genuine non-Gaussian entanglement can only be revealed by measures involving higherorder moments. In this vein, Shchukin and Vogel [83] (see also [84,85]) have introduced a general hierarchy of necessary and sufficient conditions for any state to be entangled. Nevertheless, the application of this technique to our problem turns out to be very arduous for it involves checking non-trivial inequalities, which can be performed only numerically. Moreover, the method involves the determination of moments that are extremely oscillatory and noisy [86].
In view of these difficulties, we content ourselves with assessing the purity of the reduced state of both modes. This is related to the linear entropy and intimately connected to the intermodal correlations [87]. These local purities are being the reduced density matrices of modes a and b, respectively. If we employ now the two-mode Wigner function (2.11), the purity, say P a (τ), can be written as For a bipartite system, both purities in (3.4) coincide for pure states [88,89]. In general, these quantities are different for mixed states. In our case, after a long but otherwise straightforward calculation (which, for completeness, is sketched in Appendix B), P a (τ) can be displayed as where I n (z) are the modified Bessel functions of first kind and the last equality has been carefully checked by numerical experiments. This surprising symmetry can be ascribed to the way in which the modes enter the Kerr Hamiltonian (2.2). Accordingly, we drop the mode subscripts in the purities.
As we are dealing with highly-excited fields (I 0a ≫ 1), we can make use of the asymptotic expansion [90] In addition, as τ ≪ 1 and the functions in (3.6) do not oscillate, we can replace the summation by an integral; the final result being In figure 3 we plot the time evolution of this P(τ) in the same scale as in figure 2. At τ = 0 the reduced purity is unity, in agreement with the fact that initially the state consists of two uncorrelated Gaussians. As time evolves, the purity smoothly decreases (much in a similar way as the symplectic eigenvalueν − decreases), which indicates the presence of mode correlations. This is supported by the following asymptotic estimate of (3.8) valid for I 0b τ 1. It is clear that this form of P(τ) is invariant under mode permutations. Finally, P(τ) tends to its stationary value.
One might wonder how quantum fluctuations, neglected thus far, could modify this quasiclassical picture. For the particular case of initial coherent states we are treating here, we can analytically compute the purity for the exact quantum solution. Indeed, from (2.5) we havê wherefrom one easily derive the exact expression for the purity: Using the properties of the Bessel functions, we redraft this as which explicitly exhibits the aforementioned symmetry. In fact, taking into account (3.7), P exact (τ) appears as a bidimensional Jacobi theta function [90], which is periodic. However, in the time scales we are considering here, such a periodicity is unnoticeable and we can replace again the sum by an integral, getting precisely equation (3.9).
In the inset of figure 3 we have plotted the difference between the exact solution (3.12) and the quasiclassical one (3.8). As we can see, both solutions coincide for any practical purpose. This means that the correlations examined before are of quantum nature, but higherorder correlations play no relevant role here.

Polarization squeezing
Since the polarization modes a and b have the same frequency and are orthogonal, their superposition results in a general elliptical polarization. This means that one needs only three independent quantities: the amplitudes of each mode and the relative phase between them. To describe this at the quantum level, it is advantageous to use the Stokes operators [91] complemented with the total numberN =â †â +b †b . On account of (2.1), the operators (4.1) satisfy the commutation relations of an angular momentum where the Latin indices run over {x, y, z} and ε kℓm is the Levi-Civita fully antisymmetric tensor. This noncommutability precludes the simultaneous exact measurement of the physical quantities they represent and leads immediately to the Heisenberg inequalities [92][93][94][95] where ∆ 2Â = Â 2 − Â 2 indicates the variance. It is always possible to find pairs of maximally conjugate operators for this uncertainty relation. This is equivalent to establishing a basis in which only one of the operators (4.1) has a nonzero expectation value, say Ŝ k = Ŝ ℓ = 0 and Ŝ m = 0. The only nontrivial Heisenberg inequality reads thus Polarization squeezing can then be sensibly defined by the condition [96-100] Note that squeezed states according to (4.5) are not, in general, minimum uncertainty states. The choice of the conjugate operators {Ŝ k ,Ŝ ℓ } is by no means unique: there exists an infinite set {Ŝ ⊥ (ϑ ),Ŝ ⊥ (ϑ + π/2)} that are perpendicular to the classical excitation Ŝ m , for which Ŝ ⊥ (ϑ ) = 0 for all ϑ . All these pairs exist in the S k -S ℓ plane, which is called the dark plane. A genericŜ ⊥ (ϑ ) can be written aŝ S ⊥ (ϑ ) =Ŝ k cosϑ +Ŝ ℓ sin ϑ , (4.6) ϑ being an angle defined relative toŜ k . Condition (4.5) is then equivalent to whereŜ ⊥ (θ sq ) is the maximally squeezed operator andŜ ⊥ (θ sq + π/2) the antisqueezed one. In many experiments both modes have the same amplitude but are phase shifted by π/2: â = i b . This light is circularly polarized and fulfills Ŝ x = Ŝ z = 0, Ŝ y = 0, so (4.7) directly applies.
The time evolution of the variables involved in those definitions can be evaluated using the Wigner-distribution approach: Here, WŜ ⊥ (ϑ ) (α, β ) refers to the phase-space function corresponding to the operatorŜ ⊥ (ϑ ) (commonly called its symbol). From (4.1) and (4.6) it is clear that the symbol ofŜ ⊥ (ϑ ) can be directly constructed in terms of the symbols of the basic mode amplitudesâ andb, which, from Appendix A, we know are given by Wâ(α) = α/π and Wb(β ) = β /π. Therefore, we get where I 0 = Tr[ρ(0)N] is the initial average number of photons of the state. The second-order moments are calculated much in the same way; the final result being (4.10) A major advantage of this formalism is that we can specify the time evolution of polarization squeezing. In particular, for sufficiently short times τ ≪ 1, we can expand equations (4.9) and (4.10) up to second order, so that Ŝ y (τ) τ , (4.11) so that the optimal squeezing angle is roughly given by i.e., it starts at ϑ sq = π/4 and slowly moves towards 0 as τ goes by.

Mapping the dynamics on the sphere
It is possible to turn the action the Stokes operators discussed in the previous section into a very simple phase-space picture. To this end we introduce the parametrization [101] α = √ I e iϕ a cos(θ /2) , β = √ I e iϕ a e −iφ sin(θ /2) , where ϕ a appears now as a global phase and the pertinent relative phase is φ = ϕ a − ϕ b . The radial variable represents the total intensity. The parameters θ and φ can be interpreted as the polar and azimuthal angles, respectively, on the Poincaré sphere: θ describes the relative amount of intensity carried by each mode and φ is the relative phase between them. In term of these new variables, equation (2.9) becomes In (5.1), ϕ a appears as an irrelevant global phase over which we can integrate without loosing relevant information; the result is whose solution in terms of the adimensional variable τ reads The three numbers (I , θ , φ ) are the spherical coordinates in the Poincaré space: In terms of the Cartesian counterpart equation (5.5) can be compactly expressed as where σ (θ , φ , τ) = 2 [I I 0 + S z S 0z + cos(2S z τ)(S x S 0x + S y S 0y ) + sin(2S z τ)(S y S 0x − S x S 0y )] . (5.8) For the aforementioned case of circularly polarized light, with S 0x = S 0z = 0, S 0y = I 0 , this reduces to In the x-p quadrature phase space, the usual way of representing states is by an uncertainty region which is just a contour of the Wigner function W (x, p) for that state. Much in the same way, for each fixed time, the equation W (S x , S y , S z |τ) = constant defines an isocontour surface in the Poincarś space of axes (S x , S y , S z ), which gives complete information about the fluctuations of the state. In the supplementary material of this paper, we include a movie portraying the time evolution of the Wigner function (5.7) for the particular instance in (5.9). As it can be appreciated, the state gets elongated along the direction of maximal squeezing.
In figure 4 we present three snapshots of the movie, corresponding to different times. . The orthogonal axis are S x , S y and S z , the box is centered at S x = S z = 0,S y = 10 6 and the axis ticks are measured in unit of the (spherical) isocontour at τ = 0, which corresponds to the shot-noise limit.

Dissipative effects
As light propages through the Kerr medium, it experiences a decorrelation of the relative phase between both basic polarization modes. A sensible approach to deal with this decorrelation is through the notion of decoherence, by which we loosely understand the appearance of irreversible and uncontrollable quantum correlations when a system interacts with its environment [102]. Usually, decoherence is accompanied by dissipation, i.e., a net exchange of energy with the environment. However, giving the nature of the Kerr nonlinearity, we are interested in the case of pure decoherence (also known as dephasing), for which the processes of energy dissipation are negligible. Models in which the number of photons do not change, while the coherences are strongly decaying, are at hand [103][104][105][106][107]. Surprisingly enough, however, they have not been applied in the context of the phase-number preserving Kerr dynamics. In consequence, we model such a dephasing by the master equation where Lâ[ρ] is the Linblad superoperator Using again the variables (5.6) and integrating over the irrelevant overall phase ϕ a , this equation turns out to be with γ = γ a + γ b . Its general solution can be represented by with Θ(φ |t) = ∑ k exp(ikφ − tk 2 ). In the limit I 0 ≫ 1, this exact result simplifies to The snapshots of the evolution of this Wigner function can be again appreciated in figure  4, with γ = 2γ a = 2γ b = 0.5χ. While at the beginning one can only observe a very gentle difference with the non-dissipative case, this difference gets more visible as time goes by. The shrinking of the isolevels of the Wigner function for the dissipative evolution means that it gets "smeared out" over the phase space due to the dephasing. Note that the shape and the direction of the ellipsoids are not changed; only their size is different, indicating a lower degree of polarization.
We can also investigate the impact of dephasing on squeezing. To this end we need to calculate the corresponding quantities as in equations (4.9) and (4.10). One finally gets For short times τ ≪ 1 one can show that the optimal squeezing angle is approximately given by Note that in contradistinction with equation (4.12), for a given time τ the optimal squeezing angle is closer to 0 in the presence of dephasing. In a certain sense, dephasing makes the isocontour ellipsoid rotate faster (yet also making it smaller). Finally, the optimal squeezing amount turns out to be In figure 4 we plot this optimal squeezing for several values of the ratio γ/χ. The degradation of this quantity with γ/χ can be clearly observed.

Concluding remarks
In summary, we have presented a quasiclassical approximation to the light propagation in a cross-Kerr medium. Even if the states considered are bright and we neglect quantum correlations, we still observe nonclassical effects such as entanglement or squeezing. Interestingly, in the quasiclassical limit the correlations remain in the system once induced, as opposed to the periodical decorrelation observed in the exact evolution. We have also constructed a model for dephasing processes in these media, demonstrating that dissipation does visibly affect the degree of polarization, but not so much its vectorial direction.

Appendix A. Two-mode Wigner function
In this appendix a brief review of the Wigner distribution is given for the problem at hand. For a single mode a, the Wigner function for a state given by the density matrixρ is defined as where the kernelŵ(α) readŝ so it appears as the Fourier transform of the displacement operatorD(α), witĥ Note that the standard coherent states |α are generated by the action ofD(α) on the vacuum, i.e. |α =D(α)|0 .
For a coherent state |α 0 , the Wigner function is In a more general context, the Wigner function can be interpreted as the phase-space symbol of the density matrixρ. This notion can be extended to any operatorÔ in such a way that its symbol is given by (A. 6) In particular, for the basic mode operatorâ we have In terms of W (α), we can map any operator evolution into a differential equation using the following rules [105] aρ → α + 1 2 and after performing the decomposition (2.6), this reads The rest of the properties needed in the paper can be extended to this two-mode case in a direct way.

(B.2)
In the second line, we can expand the exponential in power series in e iϕ a and e −iϕ a , considering the rest of variables as fixed coefficients. Then, the integration over ϕ a can be explicitly carried out, with the result Together with the term e −4I a , this can be immediately integrated over dI a , yielding π