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.


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:

enter image description here

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:

enter image description here

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[
    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"}]
  • $\begingroup$ 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. $\endgroup$
    – user87932
    Nov 14, 2023 at 4:19
  • $\begingroup$ It's not clear what problem you're trying to solve? $\endgroup$ Nov 14, 2023 at 10:11
  • $\begingroup$ @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! $\endgroup$
    – Saesun Kim
    Nov 14, 2023 at 13:42
  • $\begingroup$ @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 $\endgroup$
    – Saesun Kim
    Nov 14, 2023 at 14:33
  • $\begingroup$ I'll post here. I don't have an account on the other forum, so go ahead and close it, or self-answer there. $\endgroup$
    – user87932
    Nov 14, 2023 at 17:54

1 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:

enter image description here

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

  • $\begingroup$ Thank you very much for your post (+1). Could you add some picture to illustrate typo with total momentum definition? $\endgroup$ Nov 15, 2023 at 8:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.