Add Kerr nonlinearity tutorial

.gitignore

@@ -29,5 +29,6 @@ pyenv/
 /_output/
 *.html
 *.svg
+*.gif
 *.json
 /site_libs/

QuantumToolbox.jl/time_evolution/kerr.qmd

@@ -0,0 +1,201 @@
+---
+title: "Kerr nonlinearities"
+author: Li-Xun Cai
+date: 2025-02-08  # last update (keep this comment as a reminder)
+
+engine: julia
+---
+
+Inspirations taken from [this QuTiP tutorial](https://nbviewer.org/urls/qutip.org/qutip-tutorials/tutorials-v5/lectures/Lecture-14-Kerr-nonlinearities.ipynb) by J. R. Johansson.
+
+This tutorial demonstrates the use of 
+
+- [`plot_wigner`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.plot_wigner)
+- [`wigner`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.wigner)
+
+by exploring Kerr nonlinearities.
+
+## Introduction
+
+Kerr nonlinearities arise from the interaction between the electromagnetic field and a nonlinear medium with a significant third-order susceptibility, $\chi^{(3)}$. Since experiments typically use monochromatic light sources such as lasers, we restrict our discussion to a single electromagnetic mode.
+
+To derive the Hamiltonian, we start by reviewing the classical theory:
+
+1. **Nonlinear polarization density:**  
+   In a nonlinear medium, the $i$-th component of the polarization density $\vec{\mathcal{P}}$ is given by
+   $$
+   \mathcal{P}_i = 
+    \sum_j \chi^{(1)}_{ij} E_j + 
+    \sum_{jk} \chi^{(2)}_{ijk} E_j E_k + 
+    \sum_{jkl} \chi^{(3)}_{ijkl} E_j E_k E_l + \dots  
+   $$
+   for $i,j,k,l \in \{x,y,z\}$, where ${\chi}^{(m)}$ is the $m$-th order nonlinear susceptibility and $E$ is the electric field. By the assumption that only the first and third order terms are non-vanishing for our Kerr medium in the single light mode, we can approximate the polarization density as
+   $$
+   \mathcal{P} \simeq
+    \chi^{(1)} E + 
+    \chi^{(3)} E^3
+   $$
+   with the subscripts for spatial components dropped. 
+2. **Interaction of the Electric Field with the Polarization Density:**  
+   The interaction energy is proportional to $\mathcal{P} \cdot E$. The nonlinear term in the Hamiltonian is therefore 
+   $$
+   H_{\text{NL}} = \chi E^4 
+   $$
+   where $\chi$ is the effective coefficient primarily determined by the nonlinear susceptibility.
+3. **Quantization of the Electromagnetic Mode:**  
+   In quantum optics, the electric field operator for a single mode is proportional to the quardrature operator $\hat{x} = \frac{(\hat{a} + \hat{a}^\dagger)}{2}$ where $\hat{a}$ is the annihilation operator.
+
+We then combine the above properties, expand $(\hat{a}^\dagger + \hat{a})^4$, drop the terms that do not conserve the photon number or irrelevant to this tutorial with [rotating wave approximation (RWA)](https://en.wikipedia.org/wiki/Rotating-wave-approximation), and put the remaining terms in normal order. Finally, we arrive at the effective Hamiltonian
+$$
+H = \frac{\chi}{2} (\hat{a}^\dagger)^2 \hat{a}^2,
+$$
+where $\chi$ again absorbed the coefficients. 
+
+## Code demonstration
+
+```{julia}
+using QuantumToolbox
+using CairoMakie
+```
+
+We begin by defining functions for visualization:
+
+1. `plot_variance` plots the expectation value of an operator `op` and shades the variance. 
+2. `plot_Fock_dist` plots the dynamics of the Fock distribution.
+```{julia}
+function plot_variance(op, tlist, states)
+    e = real.(expect(op, states))
+    v = real.(variance(op, states))
+
+    fig = Figure()
+    ax = Axis(fig[1,1])
+    lines!(ax, tlist, e)
+    band!(ax, tlist, e .- v, e .+ v, alpha = 0.3)
+    return fig, ax
+end
+
+function plot_Fock_dist(tlist, states)
+    fig = Figure()
+    ax = Axis(
+        fig[1,1],
+        xlabel = L"N",
+        ylabel = L"t"
+    )
+
+    n_col = prod(states[1].dims)
+    n_row = length(tlist)
+
+    data = zeros(Float64, n_row, n_col)
+
+    for (idx, state) in enumerate(states)
+        data[idx, :] = real.(diag(state))
+    end
+
+    hm = heatmap!(
+        ax, 
+        0:(n_col-1), 
+        tlist, 
+        data', 
+        colormap = cgrad([:white, :magenta]),
+        colorrange = (0,1)
+    )
+    Colorbar(fig[1,2], hm, label = "Probability")
+
+
+    return fig, ax
+end
+```
+
+Next, we define the system parameters and operators.
+```{julia}
+N = 15 # Dimension of the Hilbert space
+χ = 1  # effective susceptibility
+
+a = destroy(N) # annihilation operator
+n = num(N)     # number operator
+
+# quadrature operators
+x =  a + a'    
+p = -1 * im * (a - a')
+
+# Hamiltonian
+H = 0.5 * χ * a' * a' * a * a
+```
+
+Since we are considering unitary dynamics, i.e., no dissipation, the dynamics are fully captured for $\chi t \in \left[0, 2 \pi \right]$, and the coherent initial state is representative for a laser light source. 
+
+Note that if the keyword argument `e_ops` is not supplied to [`mesolve`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.sesolve), the returned `result` contains the state at each time point at field `states`. Consult the [user guide to `TimeEvolutionSol`](https://qutip.org/QuantumToolbox.jl/stable/users_guide/time_evolution/solution) for more details.
+```{julia}
+ψ0 = coherent(N, 2.0)
+tlist = 0:0.01:(2*π / χ)
+
+result = mesolve(H, ψ0, tlist)
+```
+
+We first check the expectation value dynamics of the number operator `n` with the two visualization functions we defined previously.
+```{julia}
+fig1, ax1 = plot_variance(n, tlist, result.states)
+ax1.title = L"N"
+display(fig1)
+
+fig2, ax2 = plot_Fock_dist(tlist, result.states)
+display(fig2)
+```
+As expected, the photon number is conserved throughout. Either the expectation value or the Fock distribution are conserved in the nonlinear interaction.
+
+We now turn to the quadrature operators. The expectation value dynamics of `x` and `p` are plotted below. 
+```{julia}
+titles = [L"x", L"p"]
+for (idx, op) in enumerate([x, p])
+    _fig, _ax = plot_variance(op, tlist, result.states)
+    _ax.title = titles[idx]
+    display(_fig)
+end
+```
+There are two clear observations from these plots that indicate how the nonlinear interaction has modified the initial coherent state:
+
+1. **Non-oscillatory Behavior of Expectation Values:**  
+   In a quantum harmonic oscillator, the expectation values of the quadrature operators typically exhibit simple sinusoidal oscillations. However, under the Kerr nonlinearity, the expectation values deviate from this periodic behavior. This deviation reflects the phase distortions introduced by the nonlinear term in the Hamiltonian.
+
+2. **Departure from Minimum Uncertainty:**  
+   Coherent states in a harmonic oscillator are known to saturate the uncertainty relation, meaning they maintain a constant uncertainty product, ideally at  
+   $$
+   \Delta x \Delta p = (\frac{\hbar}{2})^2.
+   $$
+   In contrast, the evolving state in the presence of the Kerr interaction shows a time-varying variance product. This variation is indicative of squeezing effects and confirms that the state is no longer a minimum uncertainty state. 
+
+We can extend the investigation to the Wigner function with the built-in function [`plot_wigner`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.plot_wigner). 
+```{julia}
+idx = searchsortedfirst(tlist, π/χ) - 1
+fig3, ax3, hm3 = plot_wigner(result.states[idx])
+ax3.title = "χt = $(tlist[idx])"
+Colorbar(fig3[1,2], hm3)
+display(fig3)
+```
+As the plot revealed, the state at $\chi t = \pi$ is in fact the [**cat state**](https://en.wikipedia.org/wiki/Cat_state#Cat_states_in_single_modes) for a single mode. One of the main characteristics of the cat state is the superposition of two coherent states with opposite phase.
+
+The dynamics of the Wigner function offer a powerful phase-space perspective on the quantum state evolution under the Kerr nonlinearity. We again use `plot_wigner` to setup the initial plot and update frames with the function [`wigner`](https://qutip.org/QuantumToolbox.jl/stable/resources/api#QuantumToolbox.wigner) to record the dynamics as an animated GIF.
+```{julia}
+fig4, ax4, hm4 = plot_wigner(result.states[1])
+
+Colorbar(fig4[1,2], hm4)
+
+record(fig4, "kerr_wigner_dyn.gif", 1:length(tlist); framerate=24) do t
+    wig = wigner(
+        result.states[t], 
+        range(-7.5, 7.5, 200), 
+        range(-7.5, 7.5, 200)
+        )
+    ax4.title = "χt = " * string(round(tlist[t]; digits = 2))
+    hm4[3] = transpose(wig)
+end
+```
+![](kerr_wigner_dyn.gif)
+
+
+As the animation progresses, you can observe how the initial Gaussian distribution, typical of a coherent state, is gradually deformed by the nonlinear interaction. 
+
+## Version Information
+```{julia}
+QuantumToolbox.versioninfo()
+```