# Split operator method

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?

• Are you sure it's not $\gamma|\psi|^2$, which is real, unlike $\gamma\psi|\psi|^2$? You're describing a particular operator splitting method, there are others. It might help to describe it explicitly using formulas and linear algebra, instead of just words. In particular, you have to see that the Fourier transform diagonalizes $\partial/\partial x$—something I find impossible to describe in words. Feb 1, 2017 at 21:55
• @Kirill My bad. Edited the mistake. Also, I see that the fourier transform diagonalizes ∂/∂x, but my question is not that. What I am asking is that why can't we work only in the momentum space where the partial derivative operator, as you mentioned, can be written as p, i.e. ∂/∂x→p? I am talking about the split-step method in the following wikipedia page, to be specific: en.wikipedia.org/wiki/Split-step_method Feb 2, 2017 at 6:48
• Because then $\gamma|\psi|^2$ isn't diagonal? The whole point of solving $x'=Ax$ with $e^{At}x_0$ is that it's really easy when $A$ is diagonal, and hard otherwise. It's not really clear to me what you're asking. Feb 2, 2017 at 18:14
• @Kirill I guess I can put forth my difficulty more clearly when you address point 2 of my question. To re-state it, if I put an additional term of $\partial^2|\psi|^2/\partial x^2$ , then how do you evolve the system using the split step method? Feb 3, 2017 at 17:29
• Can you write out the exact expression for the whole PDE? That term looks so odd in the context of NLSE, that it would really help to be explicit. Feb 3, 2017 at 17:36

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.
• Ok, so you mean to say that I numerically evaluate $\nabla^{2}|\psi|^{2}$ and use it while calculating $exp([|\psi|^{2}+\nabla^{2}|\psi|^{2}]\delta t)$ at each time step?(I have omitted the constants). Feb 5, 2017 at 7:18