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. 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.

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.