NDSolve Time-dependent Schrodinger Equation

Question

I am currently trying to solve a time dependent Schrodinger equation (TDSE) for a rectangular potential barrier and a varying magnetic potential with a time dependent magnetic potential is oriented on the x-Axis and has a sinusoidal form.

The equation I'm trying to solve is: $$ i\omega\partial_{\xi}u(x,\xi) = -\frac{\hbar^2 c^2}{2\mu c^2}\partial_{x,x}u(x, \xi) + V(x)u(x,\xi) + i \frac{\hbar c}{\mu c^2}q* A0*cos(\xi)\partial_{x}u(x,\xi) $$

$\xi = \omega t$.

The initial state of the particle is a stationary state with energy Q > 0 found by:

hbarc = 197.326;
mu = 4000.;
S0 = 3.0;
omega = 0.0001;
ec = 0.0854;
q = 2.*ec;
Q = 5.0;
V0 = 0.;
V1 = 20.;
a = 10.;
l = 5.;
x0 = 0.;
xN = a + l + 10.;
A0 = mu omega S0/(hbarc q);
k = Sqrt[2 mu Q]/hbarc;
potential[x_] = Piecewise[{{V1,a<=x<=a+l},{V0,x<a||x>a+l}};

eqSchUnpert={Q*u0[x]==-hbarc*hbarc/(2*mu)*D[u0[x],{x,2}]+potential[x]*u0[x]};
bc={u0[xN]==beta*Exp[I*k*xN],(D[u0[x],x]/.x->xN)==Exp[I*k*xN]*beta*I*k};
sol=ParametricNDSolve[Join[eqSchUnpert,bc],u0,{x, x0, xN},beta]

and, since it is not a bound state, any beta will do so I choose 10^-4:

psi0=u0[10.^-4]/.sol

Now, when adding the magnetic potential the code takes forever and doesn't bring any result:

eqSchPert = {I*omega*
 D[u[t, x],t]==-hbarc*hbarc/(2 mu)D[u[t, x],{x, 2}]+
 potential[x] u[t, x] + I*S0*omega*Cos[t]*D[u[t, x], x]};

with the initial and boundary contitions:

icPert = {u[0.,x]== psi0[x]};
bcPert = {u[t,x0+S0]== psi0[x0+S0 + S0*Sin[t]]*Exp[-I*Q*t/omega], u[t,xN-S0]== psi0[xN-S0 + S0*Sin[t]]*Exp[-I*Q*t/omega]};

NDSolve never finishes. If anyone knows what's wrong, please tell me. What I'm doing now is:

eqSyst=Join[eqSchPert, icPert, bcPert];

NDSolve[eqSyst,u,{t,0.,2.},{x,x0+S0,xN-S0},Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","DifferenceOrder"->"Pseudospectral","MinPoints"->300,"MaxPoints"-> 300}}]

EDIT 1

After searching a little more I gave FiniteElement a try and I think I'm getting close, but I'm not there yet. If I set $\omega$ = $10^{4}$ and solve the equation with:

omega = 10.^4;
icPert = {u[r, 0.] == psi0[r]};
bcPert = {u[k*(x0 + S0), t] == 
psi0[k (x0 + S0 + S0*Sin[t])]*Exp[-I*Q*t/omega],
    u[k*(xN - S0), t] == 
psi0[k (xN - S0 + S0*Sin[t])]*Exp[-I*Q*t/omega]};
schEq = {I*omega*D[u[r, t], t] == -Q (D[u[r, t], {r, 2}]) + I*k*S0*omega*Cos[t] D[u[r, t], r] + potential[r]*u[r, t]};
eqSyst = Join[schEq, icPert, bcPert];

and I use a mesh made like this:

Needs["NDSolve`FEM`"];
mesh = ToElementMesh[Rectangle[{k (x0 + S0), 0.}, {k (xN - S0), tFinal}], "MaxBoundaryCellMeasure" -> 0.05, "MeshElementType" -> QuadElement];

the solution becomes:

solution = NDSolve[eqSyst, u, {r, t} \[Element] mesh];
u1 = u /. First[solution];
Plot3D[Abs[Evaluate[u1[r, t]]]^2, {r, k (x0 + S0), k (xN - S0)}, {t, 0., 2*\[Pi]}, PlotRange -> All, Mesh -> All]

I obtain:

which is exactly what I expected for $\omega$ very small. If I try this with the previous $\omega$ I get this:

which is absurd. In both cases, I receive a warning message:

NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

Please note that in the equations I passed to the dimensionless variable $x->\rho = k x$, with $k$ defined in the beggining of this post.

Answer 1

Minimum working code, in which all equations have solutions. I made a replacement t/omega->t to reduce computation time.

hbarc = 197.326;
mu = 4000.;
S0 = 3.0;
omega = 0.0001;
ec = 0.0854;
q = 2.*ec;
Q = 5.0;
V0 = 0.;
V1 = 20.;
a = 10.;
l = 5.;
x0 = 0.;
xN = a + l + 10.;
A0 = mu omega S0/(hbarc q);
k = Sqrt[2 mu Q]/hbarc;
potential[x_] = 
  Piecewise[{{V1, a <= x <= a + l}, {V0, x < a || x > a + l}}];
eqSchUnpert = {Q*u0[x] == -hbarc*hbarc/(2*mu)*D[u0[x], {x, 2}] + 
     potential[x]*u0[x]};
bc = {u0[xN] == beta*Exp[I*k*xN], (D[u0[x], x] /. x -> xN) == 
    Exp[I*k*xN]*beta*I*k};
psi0 = ParametricNDSolveValue[{eqSchUnpert, bc}, 
   u0, {x, x0, xN}, {beta}];
Plot[Abs[psi0[10^-4][x]], {x, x0, xN}]
eqSchPert = {I*
    D[u[t, x], t] == -hbarc*hbarc/(2 mu) D[u[t, x], {x, 2}] + 
    potential[x] u[t, x] + 
    I*S0*omega*Cos[t*omega]*D[u[t, x], x]}; icPert = {u[0., x] == 
   psi0[10^-4][x]};
bcPert = {u[t, x0 + S0] == psi0[10^-4][x0 + S0], 
   u[t, xN - S0] == psi0[10^-4][xN - S0]};

sol1 = NDSolveValue[{eqSchPert, icPert, bcPert}, 
   u[t, x], {t, 0., 2}, {x, x0 + S0, xN - S0}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 300, "MaxPoints" -> 1000, 
       "DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6];
Plot3D[Abs[sol1], {t, 0, 2}, {x, x0 + S0, xN - S0}, Mesh -> None, 
 ColorFunction -> Hue, AxesLabel -> {"t", "x"}]