# Numerical solution of driven infinite well with Floquet matrix

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.

• Note that your first expression for $U(t,0)$ is incorrect - that needs to be a time-ordered exponential. – Emilio Pisanty Mar 7 '17 at 21:04
• 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 – 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

Import["http://halirutan.github.io/Mathematica-SE-Tools/decode.m"]["http://i.stack.imgur.com/rPjH2.png"]

• 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) – yankeefan11 Mar 7 '17 at 21:51
• @yankeefan11 Sorry for the delay, the uploader was playing up. – Emilio Pisanty Mar 7 '17 at 21:52
• 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? – yankeefan11 Mar 8 '17 at 15:10
• 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)$. – Emilio Pisanty Mar 8 '17 at 18:13