Imaginary Result on NDSolve+NIntegrate

Question

I posted the same question in the Physics StackExchange, wondering if it is physics problem. I haven't received any answers, so I am posting in Mathematica StackExchange to see if it is programming issue. Please let me know if it is against any policy.

I am trying to solve the Schrödinger equation in position space.

$$i\hbar\frac{i\hbar\partial\Psi\left(x,t\right)}{\partial{t}}=-\frac{\hbar^2}{2m}\nabla^2\Psi\left(x,t\right)+V\left(x,t\right)\Psi\left(x,t\right)$$

Let's assume it is 1D in wave function with two level system atom, $\Psi(x,t)=C_g(x,t)|{g}\rangle+C_e(x,t)\exp[I(kx-wt)]|{e}\rangle$.

$V\left(x,t\right)$ is atom-light interacting potential, it is $\vec{d}\cdot\vec{E}$ potential.

Solving two equations gives:

$$i\hbar\frac{\partial C_g\left(x,t\right)}{\partial{t}}=-\frac{\hbar^2}{2m}\frac{\partial^2 C_g\left(x,t\right)}{\partial{x}^2}+\Omega \exp[-i(kx)] C_e\left(x,t\right)$$ $$i\hbar\frac{\partial C_e\left(x,t\right)}{\partial{t}}=-\frac{\hbar^2}{2m}(\frac{\partial^2 C_e\left(x,t\right)}{\partial{x}^2}+2 i k \frac{\partial C_e\left(x,t\right)}{\partial{x}})+\Omega^* \exp[i(kx)] C_g\left(x,t\right)$$

where $ \langle e| \vec{d}\cdot \vec{E} |g\rangle=\Omega \exp[- i w t]$.

I wish to solve this analytically, but it seems very difficult due to the term $\exp[i(kx)]$, so if I solve them numerically, by giving initial condition of $ C_g(x,0)=\sqrt[4]{\frac{2}{\pi }} \sqrt{\frac{1}{\text{$\sigma $x}}} \exp \left(-\frac{x^2}{\text{$\sigma $x}^2}\right),C_e(x,0)=0$, and setting constants $\hbar,\Omega,m$ to 1, $k=0.1$, and $\sigma x=8$, we got:

where the red is ground state and blue is excited state. The problem is that when I try to see the expectation value of the momentum, I am getting imaginary number. It might be either computational problem/incorrect boundary condition/some mistake of calculation.

Momentum is calculated with simple method: $\Psi^*(x,t) \frac{\hbar}{I} \frac{\partial}{\partial x} \Psi(x,t)$, then we got:

It shows the imaginary and real part of the momentum is much below what I am expecting momentum ($k=0.1$). However, momentum cannot be imaginary, so there must be a problem. Anyone has any clue what is the problem? Thanks!

Here is Mathematica code:

sx = 8.;
k := 0.1;
sol1 = NDSolve[{
   I D[a[x, t], t] == -1/2 D[a[x, t], {x, 2}] +  Exp[-I k x] b[x, t],
   I D[b[x, t], 
      t] == -1/2 (D[b[x, t], {x, 2}] + 2. k I D[b[x, t], x]) + 
     Exp[I k x] a[x, t],
   a[x, 0] == (2./Pi)^(1/4) Sqrt[1/sx] Exp[-x^2/(sx^2)],
   b[x, 0] == 0.}
  , {a[x, t], b[x, t]},
  {x, -80., 80.}, {t, 0, 10.}, MaxStepSize -> 0.1]

Plot3D[{Evaluate[Abs[a[x, t] /. sol1]^2], 
  Evaluate[Abs[b[x, t] /. sol1]^2]}, {x, -20, 20}, {t, 0, 10}, 
 PlotRange -> All, PlotPoints -> 50, PlotStyle -> {Red, Blue}]


momentum[cg_, ce_] := -I (Conjugate[cg] D[cg, x] + Conjugate[ce] D[ce, x])

axis = Table[t1, {t1, 0, 10, 0.1}];
list = Table[
   NIntegrate[
    momentum[sol1[[1, 2, 2]], sol1[[1, 1, 2]]] /. {t -> t1, 
      x -> x1}, {x1, -80, 80, 0.1}], {t1, 0, 10, 0.1}];
ListPlot[{Transpose[{axis, list // Re}], 
  Transpose[{axis, list // Im}]}, Joined -> True, 
 PlotLegends -> {"Re", "Im"}]

Answer 1

@Alex Trounev Solved in the comments, but I'll post an answer so this doesn't sit in the unanswered queue.

It looks like the question was edited, so the problem isn't visible now, but in the original posting, the iterator for NIntegrate was {x,-80.,80.,0.1}, which was likely a typo, or confusion due to the code being imbedded inside Table. NIntegrate doesn't complain because this is a valid iterator syntax for this function, used to specify contours for the complex plane, or to avoid singular points, etc. Here, the integral went from -80 to 80 and backtracked to 0.1, so effectively the integral was performed from -80 to 0.1. Somehow this asymmetry in the integration region resulted in the effect the OP saw. After correcting this, the problem went away and the complex terms are now small and consistent with numerical round off. Use of Chop would likely zero them out completely.

EDIT: @Alex Trounev requested I post a picture of the original definition. It can be found by selecting "edit" for the original question. Here's a screenshot of the original definition:

Note that the NIntegrate is over the range {x1,-80,80,0.1} instead of {x1,-80,80}.