Numerical simulation of Kerr nonlinear systems; analyzing non-classical dynamics

We simulate coherent driven free dissipative Kerr nonlinear system numerically using Euler’s method by solving Heisenberg equation of motion and time evolving block decimation (TEBD) algorithm, and demonstrate how the numerical results are analogous to classical bistability. The comparison with analytics show that the TEBD numerics follow the quantum mechanical exact solution obtained by mapping the equation of motion of the density matrix of the system to a Fokker-Plank equation . Comparing between two different numerical techniques, we see that the semi-classical Euler’s method gives the dynamics of the system field of one among two coherent branches, whereas TEBD numerics generate the superposition of both of them. Therefore, the time dynamics determined by TEBD numerical method undergoes through a non-classical state which is also shown by determining second order correlation function.


Introduction
The Kerr effect was discovered by John Kerr in 1875 [1], which exhibits quadratic electro-optic (QEO) effect, is seen in almost all materials, but certain materials display more strongly than others, for example organic molecules and polymers [2], Se-based chalcogenide glasses [3] and silicon photonic devices [4]. The non-linear phenomenon introduced by Kerr effect has been observed experimentally and it has broad range application in many optical and magnetic devices. For example, the optical Kerr effects have been useful for nonlinear signal processing which has shown several applications including NRZ-to-RZ conversion [5], multi-casting, demultiplexing, regeneration, monitoring, multiple-wavelength source [6,7], and many more. The magnetooptical Kerr effect (MOKE) has potential application in ultrathin magnetic devices, e.g. films [8], multilayers [9], and magnetic superlattices [10], and, the surface magneto-optic Kerr effect (SMOKE) has remain a powerful tool for in situ characterization [11] . All the setups have shown bistability as their predominant characteristics which causes due to nonlinear susceptibility.
The multistability in the steady state solution of Kerr nonlinearity has been encountered theoretically in two different ways: semiclassicaly where the state of the system is approximated to the nearest coherent state, and quantum mechanically where we estimate the exact solution from the master equation formalism of the density matrix of the system. The semicalssical solution of both dispersive and absorptive bistability has been derived by using the quantum Langevin equation [12,13], and the theory of quantum mechanical solution for the absorptive case [14] and the dissipative case [15] has been obtained by mapping the master equation to the Fokker plank equation. Both the techniques have been used widely to study the dynamics of open quantum systems which is considered as one of the most fundamental problems in quantum mechanics.
The theoretical techniques for open quantum system have been developed over decades and applied successfully for detecting and preparing quantum states of matter and radiation [16], sensing electromagnetic fields [17], quantum communication [18] and detection of gravitational waves [19]. The recent development of nanoscale fabrication techniques, in general, circuit quantum electrodynamics (QED) setups exhibit the technological application of quantum mechanics, particularly in superconducting qubits and nanomechanical Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. resonators ( [20,21]). Within this framework, recently the theory has been used to study various nonlinear systems, e.g. two state systems [22,23], microwave quantum optomechanics [24] and impurities in solid state systems [25]. In all cases, including Kerr nonlinear systems, linearized approximated theory has been implemented to study the dynamical behavior of the system, where one transforms the nonlinear Hamiltonian to a linear one by accounting the quantum fluctuation over nonlinear steady state field amplitude. Appreciating its simplicity, the technique, however, cannot provide a satisfactory platform. The limitation of the analytics provokes us for the numerical simulation of the time evolution of nonlinear systems, explicitly. In order to implement the numerical model, we transform the environmental degrees of freedom to a one dimensional many body system with nearest neighbor interaction and simulate the whole chain using time-adaptive density matrix renormalisation group (t-DMRG) method. The computational method consists of numerical diagonalization and renormalization process.
The t-DMRG technique is considered as one of the most powerful numerical schemes in optical, atomic and condensed matter physics to be applied on strongly-correlated many-body quantum systems. The technique have already been used for some of the renowned models of quantum mechanics, e.g. Hubbard model [26][27][28][29], Bose-Hubbard model [30][31][32] and Ising model [33][34][35], especially aiming to study the quenching dynamics, magnetization and phase transition properties.
In this paper, we introduce a TEBD numerical model for the simulation of open quantum system, and use the model for the first time to study the time dynamics of the Kerr nonlinear system. The article is composed by starting with the theoretical model that explains how the continuous modes of the bath together with the nonlinear system is transformed into a one dimensional discrete chain. Hereafter, we use two different numerical methods: time propagation of the system field by solving the Heisenberg equation of motion using Euler's method and TEBD numerical method to determine the time dynamics and steady state behavior of the system. We also compare between the numerics and with analytics, discussing in detail explaining the physical significance of our result.

Model: theory
The Hamiltonian of a system that describes the Kerr effect is, where ω S is the frequency of the cavity mode of oscillation, c¢¢ is the anharmonicity parameter which is dependent on the real part of the third order nonlinear susceptibility tensor ( [15]), and ( ) † a a are the creation (annihilation) operators of the system. E is the amplitude of an external driving field with an oscillation frequency ω L , expressing as In order to make the Hamiltonian time independent, we switch to the frame of the driving field. Eventually, the detuned cavity frequency becomes Δ=ω S −ω L . Considering the system is coupled to a thermal reservoir, the total Hamiltonian is given by, represents the Hamiltonian of a multimode bosonic reservoir which is at zero temperature, and x are the creation (annihilation) operators, and g(x) and h(x) are the frequency of oscillation and the coupling strength between the system and environment, respectively, for the environmental mode x. The properties of bath can be characterized by a uniquely defined spectral density function J(ω). Considering the linear dispersion relation (g(x)=g. x, where g is the inverse of density of states), and implying wide band limit approximation , we get the spectral density function [37] w gq w w q w 2 0 2 is the decay rate of the system and θ is the Heaviside step function. With this choice of hard cutoff, we fix the frequency limit ω c =g. x m to run over the entire spectrum of bath.

Model: tebd numerics
To simulate an open quantum system numerically, we transform the Hamiltonian of the system/bath coupling model to a semi-infinite chain model, by mapping the bath operators to the operators of lattice chain, using a unitary transformation: . In a case where the spectral density is defined by the equation (3), the normalized shifted Legendre polynomial is a natural choice as the unitary operator: and satisfies the orthogonality condition [38]. The transformed Hamiltonian of the semi infinite chain is å å . The schematic diagram of the transformation is given in figure 1(a). Similar mapping is introduced recently in [38] to simulate open quantum systems aiming to be applied to spin-boson models [39] and biomolecular aggregates [40].
Hereafter, we express the state of the chain as a matrix product state (MPS) to do the numerical simulation using TEBD. The MPS state is expressed by a a a a a a a a a = = + - The Γ and λ tensors are obtained through the Schmidt decomposition of the pure state of N sites where χ is the Schmidt number and M is the dimension of local Hilbert space. Figure 1(b) shows the method of numerical simulation for the real time evolution diagrammatically, where we choose 2nd order Suzuki Trotter (ST) expansion ( [42]) which expresses the unitary evolution operator as . The ST expansion minimizes the error in 3rd order of the time step by evolving the pairs of alternate sites.
The simulation parameters are estimated by minimizing errors which appear in two ways: during the modeling of the S/B formalism to a 1D chain and the simulation of each step of the real time evolution. We discuss the errors extensively in appendix with an estimation of simulation parameters.

Model: time propagation of semi-classical equation
Since the bath is at zero temperature, the time dynamics of the system field is obtained by the Heisenberg equation of motion: which is a first order differential equation. Therefore, using Euler's method ( [43]), we do the time propagation of the system field numerically to obtain the time dynamics and steady state.

Results: steady state situation
The exact analytical expression of the moment calculating generalized functions of the system field operators are derived by mapping the master equation into the Fokker-Plank equation [15], which determines the steady-state field amplitude and second order correlation functions: , , , is the F 0 2 hypergeometric function. Furthermore, the relation between the input drive and the semi-classical stationary value of the system field, determined from the Heisenberg equation of motion given in equation (7), is given by where α is the steady state system field.
Since the TEBD numerical model for the Kerr nonlinear system is introduced for the first time, we justify its applicability by comparing with the analytically estimated results. Here, we plot the steady state system field and the second order correlation function in figure 2 which presents both the numerical results determined through TEBD and the time propagation of the system field using Euler's method, along with the analytically determined semi-classical and quantum mechanical solution, which shows how the TEBD numerical result is analogous to classical bistability. For the numerical simulation of the time dynamics, the initial state of the system is chosen to be in a ground state, and therefore, no photon existed initially. We see the stationary value of the system field loses its linear nature when the driving field is increased far. It is also to be noted that the semi classical solution exhibits bistability whereas the exact quantum mechanical solution does not. The peak in the plot of g 2 (0) indicates the increase of the quantum fluctuations near the transition point, which happens due to the superposition of two coherent states in the quantum mechanical solution; and, as the coherent states are not mutually orthogonal, the state loses its classical nature. The TEBD determined numerical result matches to the quantum mechanical exact analytical solution, whereas the numerical time propagation of the system field using Euler's method follows the analytically determined semi-classical solution, and the branch shift occurs within the boundary of analytically predicted transition region. The extent to which bistability is observed depends on the fluctuations of the input driving field, which in turn determine the time for random switching from one branch to the other. The time scale of the change of driving field must be larger than time intervals of the random switching between branches.

Results: time dynamics
The time dynamics of the Kerr nonlinear system is generally estimated analytically by linearizing the quantum fluctuation over nonlinear steady state field amplitude. Anticipating its simplicity, however, this does not provide accuracy when the impact of nonlinearity is predominant which is observed, especially when the system is driven by a stronger pump. Even though an effort to study the classical dynamics of Kerr system has been made in [44], but this does not provide sufficient information regarding quantum dynamics; which provokes us to opt numerical methods. We plot the time dynamics of system field in figure 3 which shows that the field stabilizes after suffering initial oscillation . However, the plots exhibit difference for two different methods, which comes from the fact that the semi-classical Euler's method determines the coherent field of the system which lies on one among two branches, whereas TEBD determined result gives the superposition of both the branches, and as a consequence, the difference enhances around the transition region (figures 3(c) and (d)) . The interesting phenomenon is also noticed when we plot the trajectory of the time evolution of the system field in phase space. From the plots given in insets, we see that two different branches reach different steady states following a completely opposite trajectory. However, in case of the TEBD numerical result, we don't find any particular pattern in the transition region, for the trajectory of the system field due to the dominant superposition of both the coherent states.

Results: second order correlation function
We also plot the time evolution of the second order correlation function using TEBD algorithm, in figure 4 which shows that the correlation function does not deviate much from unit value when the system relaxes closer to a stable classical branch. However, as anticipated due to the superposition of two coherent states, around the transition region, the time evolution of the correlation function differs much from the unit, which indicates that the evolution of the system goes through nonclassical states.

Conclusion
We have used TEBD numerical technique and Euler's method successfully for the time propagation of the system field of a Kerr nonlinear system, and studied how the numerical results are analogous to classical bistability. Analyzing the steady state behavior of the system, we see that the TEBD numerical result follows the quantum mechanical exact solution obtained by mapping the equation of motion of the density matrix of the system to a Fokker-Plank equation, whereas the time propagation of the system field obtained using Euler's method follows the semi-classical solution of the Heisenberg equation of motion. The time dynamics determined by two different numerical techniques show that the semi-classical Euler's method determines the coherent field of the system which lies one among two branches, whereas TEBD determined numerical result keeps the superposition of both of them. As a result, there comes a difference in the system field for two different methods, which enhances around the classical transition region. The TEBD determined second order correlation function does not evolve as unit valued, especially around the transition region which is an indication of the generation of non-classical state due to the superposition of the two coherent states. The importance of our work, the analysis of the dynamical behavior of the externally driven Kerr nonlinear system, has been visualized in recent experiments. For example, studding the influence of different magnetic fields on electrical conductivity in a nonlinear media has drawn attention for exhibiting interesting quantum effects [45,46]. Analyzing the performance of the numerical techniques, we conclude by saying that the techniques chosen here are quite promising to work with for the analysis of nonlinear systems, and could be useful for the investigation of nonlinear dynamics reported in [24,25].

Acknowledgments
This work was supported by the Academy of Finland under contract no. 275245.

Appendix: Errors in tebd simulation and estimation of parameters
We investigate the error here introduced due to the finite values of the simulation parameters, and show how they modify with the variation of those parameters. This study will help us to estimate and optimize the parameters in order to do the simulation efficiently.

Modeling error
The modeling error is contributed by the canonical transformation of S/B coupling to 1D chain. In practice, we choose a model where the chain has a finite length considering the fact that the number of modes of the bath is finite, which causes the recurrence of the particle from the end of the chain. Here, we discuss how the recurrence time changes when those parameters change.
• Length of the chain: The recurrence time is dependent on the group velocity which is defined by where ω is the frequency, and k N is the wavenumber determined by the number of lattice sites (k N ∝N). In figure A1(a), we see the recurrence time increases with the increment of the length of the chain, which happens due to the fact that the increment of the number of sites reduces the group velocity for the particle to travel.
• Cut off frequency: The increment of the cutoff frequency increases the group velocity, forcing the particle to travel faster in the lattice, causing the reduction of the recurrence time which is seen in the figure A1(b).

Numerical error
Apart from modeling error, there are two other major sources of simulation error appears due to the finite sizes of the time step and the truncation of the Hilbert space. • Suzuki Trotter error:In case of real time evolution, the Suzuki-Trotter error which is introduced due to the finite size of the time step, tends to concentrate in the overall phase ( [42]). In figure A2(a), we show how the accuracy of the simulation improved with the reduction time step. As time step decreases the curves approach each other, and beyond the time step 10 −2 g −1 , we do not see any substantial improvement in the plot of cavity field amplitude.
• Truncation error:The TEBD numerical method is involved with the truncation of Hilbert space in every step of time evolution ( [41]). The reasonable size of Hilbert space is such a size that has the ability to express the cavity field with negligible error. The coherent field generated in the system has a Poisson probability distribution  2 . The cumulative probability distribution is shown in figure A2(b), which estimates the reasonable size of local Hilbert space, anticipating accuracy upto a significant extent. In order to estimate a reasonable size of the Schmidt number, we plot the von Neumann entropy associated with the entanglement between the system and the first site in figure A2(c). As the evolution of the state starts from a product state, the von Neumann entropy was zero at the beginning. It increases initially and reaches to the maximum value, and then, reduces with time. We estimate the Schmidt number when the entanglement is maximized. The Schmidt number should be chosen in such a way that the eigenvalues after truncation contribute so less that they can easily be neglected. In figure A2(d) we plot the eigenvalues of the reduced density matrices when the entanglement is maximum between the system and the first site of the chain. The truncation line shows that the choice of Schmidt number is reasonable in this case.