Computing flux modulation of the energy spectrum in a DC SQUID I've been reading some work by Y. Chen et al's paper on tunable couplers for transmon qubits (you can find the work at https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.113.220502) and wanted to model the circuit but am running into an issue. As I understand it, the circuit can be approximated as a lumped element parallel LCJJ. I'm interested because this is a general class of LC oscillators with a tunable inductance given by the DC and RF SQUIDs in the lumped element circuit limit. 
So I write down the hamiltonian in terms of relevant energies and junction asymmetry $\alpha := E_{J_{2}}/E_{J_{1}}$, and $\chi := \frac{1-\alpha}{1+\alpha}$ and I get
\begin{align}
H_\text{total} &= 4E_c \frac{d^2}{d\phi^2} + \frac{E_L}{2}(\phi-\phi_{ext})^2 \\
&- (1+\alpha) E_{J_{1}} \cos \left(\pi\frac{\Phi}{\Phi_0} \right) \sqrt{1 + \chi^2 \tan \left(\pi \frac{\Phi}{\Phi_0} \right)} \\
& \cos \left(\phi - \arctan \left(-\chi \, \tan \left(\pi\frac{\Phi}{\Phi_0} \right) \right) \right)
\end{align}
which as I understand it is all fine and well with $E_c$ and $E_L$ being the capacitance and inductive energies respectively, $\Phi_0$ the flux quantum, $\phi_{ext}$ the RF squid external bias, and $\Phi$ the DC squid bias.
As a sanity check, I decided to see if I could reproduce the flux tuning of the 0-1 resonance for a symmetric DC SQUID in the transmon regime. I set $\alpha  = 1$, $E_L = 0$,$E_c = 0.3$ GHz, $E_J = 20$ GHz and have the reduced hamiltonian:
$$H_\text{Transmon} = 4E_c \frac{d^2}{d\phi^2} - 2E_{J} \cos \left(\pi\frac{\Phi}{\Phi_0} \right) \cos(\phi) \, .$$
To solve this I discretize the phase $\phi$ and $\Phi/\Phi_0$ over the interval $[-\pi,\pi]$ and $[-2,2]$ respectively, generate the corresponding matrix equation in matlab, and compute the eigenvalues using eig() function. I expect the first transition energy to be approximately $E_{01} (\Phi=0) \approx \sqrt{8E_cE_J}-E_c = 6.628 GHz$ and to modulate to zero frequency every half flux quantum. See below for what I expect to find for various levels or junction loop asymmetry (note this guy has the label asymmetry for $\chi$ in my notation):

For the $\alpha = 1$ case I get the below plot for $f_{01}$ and $f_{12}$. The peak frequency and modulation are wrong. Notably, it does get the nominal anharmonicity correct at 0 flux bias. I have noticed that the outcome changes somewhat by changing the window and step size in phase/flux to be multiple periods with smaller slices which makes me suspect the periodicity and convergence are issues at play. A friend suggested it might be due to the fact that the junction phase shifts as a function of flux bias and that I might need to solve in phase about $arctan(-\chi tan(\pi\frac{\Phi}{\Phi_0})$. I tried deleting the phase shift as a fast test but it doesn't seem to fix the issue. Any suggestions? I can also supply the code I wrote.

 A: After some more investigation with a friend who works with superconducting circuits I have the solution. The first mistake was I wasn't giving it the values I thought, so the peak frequency error was essentially a result of me setting the individual Josephson energies to $E_{J\Sigma} = (1+\alpha)E_{J1}$.  
The bigger issue with flux biasing is actually somewhat subtle (at least to me, someone who doesn't work with flux biased circuits). Even though the potential for the junction loop is $-E_{J\Sigma}\cos\left(\pi\frac{\Phi}{\Phi_0}\right)\cos(\phi)$ to get the "correct" behavior at multiple flux quanta you need to diagonalize with $-E_{J\Sigma}|\cos\left(\pi\frac{\Phi}{\Phi_0}\right)|\cos(\phi)$ for the flux bias in the transmon regime (or at least when $E_L = 0$). I don't have a clean explanation for why it should be simulated this way. I know when the sign changes the location of the potential minima move, but to me that seems like it should be accounted for already during diagonalization. 
Another thing of note that fundamentally limits you even in the perfect junction symmetry case from achieving zero frequency at half a flux quantum is the fact you fail to have a bound state solution for a given "kinetic" energy of $4E_c$. In my simulations it appears you lose the bound state around $3E_c$. I am sure we could probably derive that exact lower bound for a given $E_c$ but that is an exercise for another time. 

