# Solving Schrödinger equation for Dirac comb potential (kicked rotor)

I need to solve the Schrödinger equation for a Dirac delta potential. I could not find the correct way to write the time-dependent potential and how to solve the time-dependent equation for it.

The time-dependent Hamiltonian is given as

$$\mathcal H = \frac{p^2}{2} + k_0V(q)\sum_n \delta(t-nT), \quad V(q)=\cos(q) + i\gamma q,$$

where $$L = π, k_0 = 5, γ = 0$$ and $$T=1$$.

L = π; k0 = 5; γ = 0; T=1
U[q_] := Cos[q] + I γ q;
V[q_, T_] := k0 U[q] Sum[DiscreteDelta[t - n T], {n, Infinity}]


There is a problem with writing V[q,t]. How can I correct it and find the eigenvalues and eigenstates?

{ℒ, ℬ} = {-1/2*Laplacian[u[q, t], {q, t}] + V[q, t]*u[q, t],
DirichletCondition[u[q, t] == 0, True]};

{ev, ef} = NDEigensystem[{ℒ, ℬ}, u[q, t], {q, -L, L}, {t, 0, 1 }, 5]


Is this correct?

• If you're going to use a numerical routine to solve the equation, you'll run into problems since DiscreteDelta is non-zero only for exact integral values. The numerical routine will vary parameters, treating them as reals, and the likelihood of getting an exact integer is zero, due to numerical precision. You'll probably have to define individual delta functions as narrow peaks, and vary the widths to suit your needs.
– user87932
Aug 30, 2023 at 18:18
• I understood your point. but how should I do that in Mathematica? Could you share a hint or a piece of code Aug 30, 2023 at 18:35
• Are you sure in γ = 0? Aug 31, 2023 at 18:58
• @user64494, Yes, it is. For reference, you can check this DOI: 10.1103/PhysRevLett.104.054102 2nd page Sep 1, 2023 at 3:58
• Eigenvalues are for time-independent Hamiltonian: H psi = E psi. Are you looking for eigenvalues between the delta function pulses? Or for a time-dependent solution, you could break into steps, integrating the equation between each delta pulse and adjust for the instantaneous effect of the pulse to go to the next step.
Jan 20 at 15:58

One problem with your approach is that DiscreteDelta is designed to accept exact integers. A numerical equation solver works with real numbers, and the sampling won't select an exact integer, so DiscreteDelta won't trigger. You need to write your own version of DiscreteDelta, which accepts real numbers as input. For example, something like this:

comb[width_][t_] := Piecewise[{{1/width, Abs[t] <= width/2}}, 0];


The comb function replaces each Dirac delta function with a square wave of a specified width. You can experiment with the value and vary it to get a suitable setting. The 1/width is output when the t value is inside the pulse; the delta function integrates to 1, so this is a normalization factor to get the correct area.

Here's a plot showing the comb function's behavior:

Plot[Table[comb[1/10][t - n], {n, -2, 5}], {t, -1, 4},
Exclusions -> None]


Here's your code with comb replacing DiscreteDelta. I added width and nmax to the list of parameters; summing to Infinity won't work well for a numerical function, so I added a finite cutoff. You can also adjust this to whatever you want. I chose width - 0.1 and nmax = 5 to test this out.

L = Pi; k0 = 5; gamma = 0; T = 1; width = .01; nmax = 5;
U[q_] := Cos[q] + I gamma q;
V[q_, T_] := k0 U[q] Sum[comb[width][t - n T], {n, nmax}]

{lap, bc} = {-1/2*Laplacian[u[q, t], {q, t}] + V[q, T]*u[q, t],
DirichletCondition[u[q, t] == 0, True]};

{ev, ef} =
NDEigensystem[{lap, bc}, u[q, t], {q, -Pi, Pi}, {t, 0, 1}, 5]


There were 5 solutions as you requested. Here are the eigenvalues; I left the eigenfunctions out because pasting them here converts them to OutputForm and they're huge.

{{5.05996, 5.43496, 6.05997, 6.93497, 8.06}, <left out IFuncs>}


Here are plots of the 5 eigenfunctions:

GraphicsGrid[{Table[
Plot3D[Evaluate[ef[[n]]], {q, -Pi, Pi}, {t, 0, 1}], {n, 5}]}]


You can experiment with width and nmax and try to determine how sensitive your results are to the settings.

• Never heard before of eigenfunctions with a time dependent potential. A constant potential switched on for a short time does not change anything. It will reduce the kinetic energy for that short time interval, in case of a DiracDelta to zero. As a mechanical system it does not fit resonably boundary conditions of a hyperbolic system. Aug 31, 2023 at 4:21