I am trying to understand and replicate figure 3 of this paper.

The idea is that we have our box with a periodic drive $$H(t) = \frac{p^2}{2m}-F_0x\cos(\omega t)$$ The author states

In order to calculate the truncated time-evolution matrix $U(T, 0)$ one takes each of the unperturbed box eigenstates as initial condition, $ψ_n(x, 0) = ϕ_n(x)$, and computes the states $ψ_n(x, T )$ resulting after one period, collecting these state vectors as columns of the monodromy matrix

I am having difficulty figuring this out. I know the unperturbed energies and wavefunctions $|\psi(0)\rangle$. I also know that $|\psi(t)\rangle = U(t,0)|\psi(0)\rangle$ with ($\hbar=1$)

$$U(t,0)=\exp\left[-i\int_0^tH \ dt\right]$$ This is where I am unsure of how to numerically calculate. Usually, I would do just $\psi(t)=\psi(0)\exp[-iE_n t]$, but because of the periodicity, we actually are able to expand $$U(t,0)=\sum_n |n\rangle\langle n| \exp[-i\epsilon_n T]$$ In quasi-energies, which I want to solve for, so I am not sure of the time dependance of H to integrate to solve for $|\psi(x,t)\rangle$. Any suggestions on how to approach this are appreciated.

  • $\begingroup$ Note that your first expression for $U(t,0)$ is incorrect - that needs to be a time-ordered exponential. $\endgroup$ – Emilio Pisanty Mar 7 '17 at 21:04
  • $\begingroup$ Even if I have the OE, I am not sure how to generate a state $|\psi(t)\rangle$ from U(t,0)$|\psi(0)\rangle$ due to the integration of p^2 $\endgroup$ – yankeefan11 Mar 7 '17 at 21:30

Here there's no fancy math that will save you - you just need to explicitly solve the TDSE numerically for a bunch of different initial conditions.

The algorithm looks more or less like this:

  • For each $m=1,\ldots,N_\mathrm{max}$:

    • Formulate the Schrödinger equation for the vector $(\psi_n(t))_{n=1}^{N_\mathrm{max}}$ as a set of coupled ordinary differential equations, $$i\hbar \partial_t \psi_n(t) = E_n\psi_n(t)+\sum_{k=1}^{N_\mathrm{max}} Fx_{nk}\psi_k(t) \quad\text{under}\quad \psi_n(t) = \delta_{mn}$$
    • Solve this set of coupled ODEs numerically up to $T=2\pi/\omega$.
    • Store the results at the end of the period as $U_{mn} = \psi_n(T)$. This is known as the monodromy matrix.
  • Assemble the results of all the different $m$ into one big matrix $U_{mn}$.

  • Diagonalize this matrix. Its eigenvalues are of the form $e^{i \varepsilon_k T/\hbar}$, i.e. the exponentiated quasienergies.
  • Repeat for a different coupling.

This is what the authors mean when they say

then [we] have to deal with a system of $n_\mathrm{max}$ complex coupled ordinary differential equations which can be integrated numerically by standard routines.

I gave this a brief run in Mathematica and it's not that hard (though of course this gets harder and harder as $n_\mathrm{max}$ gets bigger). The results look roughly like this:

I haven't checked the constants for normalization or the convergence with respect to $n_\mathrm{max}$, so take that as a rough guide only; the code is at


  • $\begingroup$ This is awesome. I think I was getting all confused over trying to integrate the p^2, when that is just the En Yn. Any chance you could post the Mathematica code (I am dooing this is matlab) $\endgroup$ – yankeefan11 Mar 7 '17 at 21:51
  • $\begingroup$ @yankeefan11 Sorry for the delay, the uploader was playing up. $\endgroup$ – Emilio Pisanty Mar 7 '17 at 21:52
  • $\begingroup$ What is the significance of the delta function? In your code, you have initialCondition[m_] := Table[ [Psi][n][0] == KroneckerDelta[m, n] , {n, 1, Nmax}] Which is just an Nmax identity matrix? $\endgroup$ – yankeefan11 Mar 8 '17 at 15:10
  • $\begingroup$ That's the initial condition $\psi_n(0)=\delta_{mn}$ above. In essence, you are calculating $U(T,0)|\psi_m\rangle$ for all the different $m$, and then collecting the results together to approximate (the truncation of) $U(T,0)$. $\endgroup$ – Emilio Pisanty Mar 8 '17 at 18:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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