Split operator method

Question

I am new to the field of computational physics and have a couple of questions regarding solving the non-linear Schrödinger equation using Operator splitting.

1) If the hamiltonian is of the form $H=\frac{\partial^{2}}{\partial x^{2}}+\gamma|\psi|^{2}$ then the standard procedure I understand is to exponentiate $-i(\gamma |\psi|^{2})\Delta t/\hbar$ and operate it on the initial value of $\psi$, then take a fourier transform to convert it to momentum space and operate it with exponential of $-ip^{2}\Delta t/\hbar$ and convert the resultant back to position space. We repeat this for each time interval $\Delta t$. Instead, why can't we do everything in the momentum space to begin with? Why this back and forth shifting from position to momentum space?

2) Suppose now I have an additional term of $\frac{\partial^{2}}{\partial x^{2}}|\psi|^{2}$ in the Hamiltonian, then how do I accomodate this term in the scheme of split operator method?

Answer 1

The Strang splitting method goes like this. You start with the PDE $$ u_t = (L+B)u, \qquad L = \partial_x^2, $$ and you notice that when $L$ and $B$ are independent of $x$, the exact solution after time $\delta t$ to this is $$ u = e^{(L+B)\delta t}u_0 \approx e^{\frac12 L\delta t}e^{B\delta t}e^{\frac12 L\delta t}u_0. $$ Because $B$ is just a function of $x$, the operator $e^{B\delta t}$ just multiplies by $e^{B(x)\delta t}$. Because $L = F \hat L F^{-1}$ is diagonalizable by the Fourier transform (assuming the right boundary conditions), $$ e^{L\delta t} = F e^{\hat L \delta t} F^{-1},$$ where $e^{\hat L\delta t}$ is the operator that multiplies each Fourier mode by $e^{-k^2\delta t}$ (this depends on choice of normalization).

This is why Fourier transforms are done at each step: in the Fourier basis, and only in that basis, is $L$ diagonal, which makes it trivial to compute its exponential.

For your equation, to get what $e^{B\delta t}$ would look like, you write out the relevant portion of the equation, with only $B$ present: $$ i\hbar u_t = (\gamma |u|^2 + \alpha (|u|^2)_{xx})u. $$ One thing you could do is to approximate $g(x) \approx (|u_0|^2)_{xx}$, so that $$ u(t,x) \approx u(0,x)\exp\left(\frac{\gamma|u(0,x)|^2 + \alpha g(x)}{i\hbar}\,\delta t\right). $$ Because of the nonlinearity, it might work, but I think there isn't a guarantee that it will—I haven't tried it. But the idea is still the same: split the r.h.s. into two operators, and for each operator solve the corresponding PDE, choosing the operators in a way that makes this step easy.