# Imaginary Result on NDSolve+NIntegrate

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"}]
• I tried running your code. NDSolve complains about your boundary conditions (ref/message/NDSolve/bcart), and NIntegrate has trouble with convergence on some of your points. In your NIntegrate, you have {x1,-80,80,0.1}. What does the 0.1 do? My recollection is that at one time, you could specify a contour but why would you do that here from 80 to 0.1? I removed it. No sign of your problem, and if I run Max[Im[list]], I get 4.08585*^-9, which is quite a bit smaller than the value you got and likely due to numerical round off etc. So the main change was removing that 0.1 from NIntegrate.
– user87932
Commented Nov 14, 2023 at 4:19
• It's not clear what problem you're trying to solve? Commented Nov 14, 2023 at 10:11
• @jdp Thank you so much for your comment! It indeed solved a problem. would you please wrap up and answer the question? so I can select and close the post? Thank you! Commented Nov 14, 2023 at 13:42
• @jdp I am planning to close the Physics StackExchange question, but feel free to post the answer there too, I will select your answer. Thanks Commented Nov 14, 2023 at 14:33
• I'll post here. I don't have an account on the other forum, so go ahead and close it, or self-answer there.
– user87932
Commented Nov 14, 2023 at 17:54