\section{Action variable for the classical two-mode Bose-Hubbard model with $0 \leq u < 1$}

We consider the single-dimer Hamiltonian $H_r$ in Eq.\ (3), with particle number $n_r=|\alpha_{1r}|^2+|\alpha_{2r}|^2$. We henceforth drop the $r$ subscripts, so that we begin from
\begin{equation}\label{Hsup1}
H = -\frac{1}{2}(\alpha^*_1\alpha_2 + \alpha^*_2\alpha_1) +\frac{U}{4}n^2+\frac{U}{4}(n_1-n_2)^2\;.
\end{equation}

As common in studies of the two-mode Bose-Hubbard model we next transform to the Schwinger representation of angular momentum,
\begin{align}
\frac{\alpha_1^*\alpha_2+\alpha_2^*\alpha_1}{2} &= L_1\nonumber\\
\frac{\alpha_1^*\alpha_2-\alpha_2^*\alpha_1}{2i} &=L_2\nonumber\\
\frac{|\alpha_1|^2-|\alpha_2|^2}{2} &=L_3\;.
\end{align}
Although more commonly introduced quantum mechanically, where it is justified by the fact that the canonical commutation relations of the bosonic creation and destruction operators reproduce the angular momentum algebra, this transformation also works classically, with the classical Poisson brackets replacing the quantum commutators. We further note that this representation implies $L_1^2+L_2^2+L_3^2\equiv (n/2)^2$. With this transformation we render (\ref{Hsup1}) into
\begin{equation}
H = -L_1 +U L_3^2+ \frac{Un^2}{4}\;.
\end{equation}

We then introduce a canonical representation of the angular momentum variables for fixed $n$:
\begin{align}
L_1 &= P\nonumber\\
L_2 &= \sqrt{(n/2)^2-P^2}\cos\phi\nonumber\\
L_3 &= \sqrt{(n/2)^2-P^2}\sin\phi\;,
\end{align}
where $\phi,P$ are a canonically conjugate pair with the finite ranges $-n/2\leq P \leq n/2$ and $-\pi < \phi \leq \pi$, $\phi =\pm\pi$ being identified so that we have a finite cylindrical phase space. (It is easy to confirm that the canonical Poisson brackets of $\phi$ and $P$ reproduce those of the angular momentum variables.) This delivers
\begin{equation}
H = -P + U[(n/2)^2-P^2]\sin^2\!\phi + Un^2/4\;.
\end{equation}
A key feature of this our final transformation is that for $0<Un<1$ all the contours of constant $H=\varepsilon$ sweep through the entire $[0,2\pi]$ range of $\phi$. 

For $0<Un<1$ it can also be shown that $Un^2/4-n/2\leq\varepsilon Un^2/4+n/2$.  We therefore represent the energy in terms of the angle $\eta\in [-\pi/2,\pi/2]$ according to
\begin{equation}
\varepsilon = \frac{Un^2}{4} + \frac{n}{2}\sin\eta
\end{equation}
and express the contour of constant energy $H=\varepsilon$ in the $(\phi,P)$ phase space as
\begin{equation}
U\sin^2\phi\, P^2 + P + \frac{n}{2}\sin\eta - \frac{Un^2}{4}\sin^2\!\phi= 0\;.
\end{equation}
Solving the quadratic equation, and discarding one branch because we must have $P\geq -n/2$, gives
\begin{align}
P(\phi,\eta,n) &= \frac{1}{2U\sin^2\phi}\left(-1 + \sqrt{1 - 2Un\sin\eta\,\sin^2\phi + (Un)^2\sin^4\phi}\right)\nonumber\\
&\equiv\frac{n}{2u\sin^2\phi}\left(-1 + \sqrt{1 - 2u\sin\eta\,\sin^2\phi + u^2\sin^4\phi}\right)\equiv n X(\phi,\eta,u)\;,
\end{align}
for $u=Un$. 

The area enclosed for given $\eta$ between this contour $-n/2\leq P(\eta,\phi)\leq n/2$ and the minimum energy contour $P(-\pi/2,\phi)=n/2$ is then by definition equal to $2\pi$ times our action variable $m(\varepsilon,n)$:
\begin{align}
m(\varepsilon,n) &= \frac{n}{2} - \oint\!\frac{d\phi}{2\pi}\,P\Big(\phi,\eta(\varepsilon,n),n\Big)\nonumber\\
& = n\left(\frac{1}{2}-\oint\!\frac{d\phi}{2\pi}\,X\Big(\phi,\eta(\varepsilon,n),u\Big)\right)\;.
\end{align}

Defining the new integration variable $z=\sin^2\phi$ and then integrating by parts allows us to write $m(\varepsilon,n)$ as
\begin{equation}
m = n\left(\frac{1}{2}- \frac{1}{\pi} \int^{1}_{0} dz \sqrt{\frac{1-z}{z}} \frac{uz-\sin \eta}{\sqrt{(1+iue^{i\eta}z)(1-iue^{i\eta}z)}} \right)\;.
\end{equation}
The integral over $z$ is in the general class of elliptic integrals but it does not reduce immediately to any of the three standard Legendre forms of elliptic integral. Mathematica evaluates it exactly in terms of Appell $F_1$ functions, giving a form which is manifestly real but unfamiliar and also slow to compute numerically. If instead of taking $z=\sin^2\phi$ we change variables to $\xi=\cot\phi$ we can obtain an exact evaluation in terms of incomplete elliptic integrals at infinite argument. Applying elliptic integral identities we can then reduce these to complete elliptic integrals of the three basic Legendre types, and confirm by numerical plotting that the result is exactly equal to the exact result in Appell functions. Most of our analytical results in this paper are based on this evaluation:
\begin{align}
\frac{m}{n} =&\frac{1}{2} - \frac{1}{\pi}\mathrm{Re}\left[ 2\cos\eta\frac{K\!\!\left(\frac{1-i e^{-i \eta } u}{1+ie^{i \eta }u}\right)}{\sqrt{1+i e^{i \eta } u}}+\frac{2i}{u}\sqrt{1+i e^{i \eta } u}\, E\!\!\left(\frac{1-i e^{-i \eta} u}{1+ie^{i \eta }u}\right)
 +\frac{i e^{i \eta } \left(1-i e^{-i \eta } u\right) \Pi\!\!\left(i e^{-i \eta } u|\frac{2 i u \cos  \eta }{1+ie^{i \eta } u}\right)}{\sqrt{1+i e^{i \eta } u}}\right]\;,
\end{align}
where $K(m)$, $E(m)$, and $\Pi(n|m)$ are respectively the complete elliptic integrals of the first, second, and third kinds. 

This final expression is admittedly ponderous but it evaluates efficiently numerically and the complete elliptic integrals obey a number of well-known identities, so that derivatives of $m$ with respect to $\varepsilon$ can also all be evaluated exactly in terms of complete elliptic integrals. In particular the microcanonical entropy, temperature, and specific heat can all be obtained straightforwardly, because through further elliptic integral identities the derivative of $m$ with respect to $\varepsilon$ simplifies dramatically to just
\begin{equation}
\frac{\partial m}{\partial\varepsilon} = \frac{2}{\pi}
\frac{K\!\!\left(
	1-\frac{1-i e^{-i \eta} u} {1+ie^{i\eta}u}
	\right)}
	{\sqrt{1+i e^{i\eta}u}}\;,
\end{equation}
which is real in spite of its complex arguments because of an identity satisfied by complete elliptic integrals of the first kind. 

Since $2\pi m$ is the action variable for the Bose-Hubbard dimer, $2\pi m'(\varepsilon)$ not only gives the quantum density of states in the semiclassical limit, and therefore the microcanonical entropy, but is also the period of the classical orbit at energy $\varepsilon$. $K(0)\equiv\pi/2$, so we can confirm that
\begin{equation}
\lim_{\varepsilon\to\varepsilon_\pm} \frac{\partial m}{\partial\varepsilon} = \lim_{\eta\to\pm\pi/2}\frac{\partial m}{\partial\varepsilon} = \frac{1}{\sqrt{1\mp u}}\;,
\end{equation}
which agree with the inverse Josephson frequencies of small oscillations around the ground and highest excited states.

The microcanonical temperature can be expressed in terms of complete elliptic integrals as well:
\begin{equation}
\frac{1}{k_BT} = \frac{u-\sin\eta}{2n\cos^2\eta (1-i e^{-i\eta}u)}\frac{E\!\!\left(
	1-\frac{1-i e^{-i \eta} u} {1+ie^{i\eta}u}\right)}{K\!\!\left(
	1-\frac{1-i e^{-i \eta} u} {1+ie^{i\eta}u}\right)}+\frac{\sin\eta}{2n\cos^2\eta}\;,
	\end{equation}
which is again identically real in spite of its complex arguments. The temperature is also finite at minimum and maximum energy in spite of the $1/\cos^2\eta$ factors, and finite for $n\to0$ as long as $U>0$. The limits are
\begin{equation}
k_B T_\pm = k_B T(\varepsilon_\mp) = \frac{(1\pm u)^2}{2U(1-\pm u/4)}\;.
\end{equation}


\section{Critical $Un$ for adiabatic invariance of $M$ from Bogoliubov stability analysis of the maximally excited array}

The calculation we perform to fix $u_c$ for the onset of $M$ decay is actually linear stability analysis of the maximally excited array state $\alpha_{1,r}=-\alpha_{2,r}=\sqrt{n/2}$. Since this maximally excited state has maximum $M$ value, its instability necessarily involves decrease of $M$. The reason that the stability of this very particular homogeneous state can be used to diagnose the invariance of $M$ for generic states is that the instability which we will find for the homogeneous state will turn out to be at the \emph{shortest} possible wavelengths. It is thus actually a local instability which can appear whenever even a few neighboring dimers are maximally excited. We have seen in our main text that the cases in which $M$ is adiabatically invariant, and in which it matters that $M$ is adiabatically invariant, are cases with bimodal dimer ensembles in which maximally excited dimers are not uncommon. The calculation of $u_c$ from the maximally excited state of the whole array thus shows when energy bimodality due to adiabatic invariance of $M$ becomes self-inconsistent, because above this $u_c$ threshold invariance of $M$ will lead to clusters of maximally excited dimers in which $M$ will spontaneously decay.

To perform the linear stability analysis we begin with the Hamiltonian equations of motion according to (1) from our main text, which are the two-component discrete Gross-Pitaevskii (GP) nonlinear Schr\"odinger equation
\begin{align}
i\dot{\alpha}_{1,jk} &= -\frac{J}{2} \sum_{\langle r,r^{\prime}\rangle} \alpha_{1,r^{\prime}} - \frac{\alpha_{2,r}}{2}  + U|\alpha_{1,r}|^2 \alpha_{1,r} \\
i\dot{\alpha}_{2,r} &= -\frac{J}{2} \sum_{\langle r,r^{\prime}\rangle} \alpha_{2,r^{\prime}} - \frac{\alpha_{1,r}}{2}  + U|\alpha_{2,r}|^2 \alpha_{2,r} \;,
\end{align}
where $r$ is short for the site indices $j,k$ in the 2D array.

We then linearize the above equations around the array's homogeneous maximally excited state $\alpha_{1,r}=-\alpha_{2,r}=\sqrt{n}$, which has $i\dot{\alpha}_{\sigma,r}=\nu\alpha_{\sigma,r}$ for
\begin{equation}
\nu = \frac{Un+1-J}{2}\;.
\end{equation}
Since the background about which we linearize has discrete translation symmetry in the $L\times L$ torus, we assume a discrete plane wave Ansatz
\begin{equation}
\alpha_{\sigma,jk}(t) = e^{-i\nu t}\left((-1)^\sigma\sqrt{\frac{n}{2}} + a_\sigma e^{\frac{2\pi i}{L}(p_1 j + p_2 k)}e^{-i\omega t} + b^*_\sigma e^{-\frac{2\pi i}{L}(p_1 j + p_2 k)}e^{-i\omega t}\right)
\end{equation}
for integers (discrete two-dimensional wave numbers) $\mathbf{p} = (p_1,p_2)$, linear excitation frequency $\omega(\mathbf{p})$, and infinitesimal two-component amplitudes $a_\sigma$, $b_\sigma$. Writing $u=Un$ as in the main text, this yields the time-independent Bogoliubov-de Gennes (BdG) equations
\begin{equation}
\omega \left(\begin{matrix} a_1 \\ a_2 \\ b_1 \\ b_2\end{matrix}\right) = \left(\begin{matrix} u + \mathcal{J}_\mathbf{p}-2 -\nu & \frac{1}{2} & \frac{u}{2} & 0 \\ \frac{1}{2} & u + \mathcal{J}_\mathbf{p}-2 -\nu & 0 \frac{u}{2}  \\ -\frac{u}{2} & 0 & -(u + \mathcal{J}_\mathbf{p}-2 -\nu) & -\frac{1}{2} \\ 0 & -\frac{u}{2} & -\frac{1}{2} & -(u + \mathcal{J}_\mathbf{p} -2-\nu)\end{matrix}\right)\left(\begin{matrix} a_1 \\ a_2 \\ b_1 \\ b_2\end{matrix}\right)
\end{equation}
where $\mathcal{J}_\mathbf{p} = 2J[\sin^2(\pi p_1/L)+\sin^2(\pi p_2/L)]$. 

The eigenfrequencies $\omega(\mathbf{p})$ of the four-by-four matrix have the two branches $\pm \sqrt{\mathcal{J}_\mathbf{p}(\mathcal{J}_\mathbf{p} + u)}$, which is always real, and
\begin{equation}
\omega =\pm \sqrt{(1 - \mathcal{J}_\mathbf{p})(1 - u - \mathcal{J}_\mathbf{p})}
\end{equation} 
which is imaginary, indicating dynamical instability, when $1-u < \mathcal{J}_\mathbf{p} < 1$. Since $\mathcal{J}_\mathbf{p}$ can always be made smaller by choosing smaller $\mathbf{p}$, the requirement for $\omega$ to be real for \emph{all} $\mathbf{p}$ is the requirement that the largest possible value $\mathcal{J}_\mathbf{p}$ is still smaller than $1-u$. The largest possible value of $\mathcal{J}_\mathbf{p}$ occurs for maximum $p_1 = p_2 = L$, giving $\mathcal{J}_\mathbf{p}\to 4J$. Thus the onset of instability, which occurs at the shortest possible wavelengths in the array, is $4J > 1- u$, so that the instability threshold is $u_c = 1-4J$ as stated in our main text. For the numerical evolutions shown in our main text we have used $J = 0.05$, giving $u_c = 0.8$.

We can confirm by numerical evolution of the array that the total dimer action $M=\sum_r m_r$ really does remain adiabatically invariant as long as the average $u<u_c$ but ceases to be invariant for average $u>u_c$. We do this by taking the initial state in which all the dimers of the array are at maximum energy ($\alpha_{2r}=-\alpha_{1r} = \sqrt{n/2}$), then adding a small random perturbation, and evolving over a long time. In Fig.\ \ref{joson_stability}a we have plotted the evolution of $M$ over time for two cases of average $u$, one just below $u_c$ and one just above it. The late-time value of $M$ is then plotted as a function of $u$ in Fig.\ \ref{joson_stability}b. The total $M$ is quite accurately conserved for $u$ below critical $u_c$, but for $u$ just above $u_c$ we see that $M$ drops substantially. (It does not necessarily decay to zero; total energy remains conserved and the single-dimer energy held in $m_r$ attains an equilibrium with the energy stored in long-wavelength sound waves in the Bose-Hubbard array.)  As explained above, the breakdown of adiabatic invariance of $M$ is due to short-wavelength instability in the array and so any small cluster of maximally excited dimers, such as tends to form if $M$ is invariant, will actually allow $M$ to decrease if $u> u_c$ in the cluster. 
\begin{figure}[ht]
\centering
\includegraphics[clip=true,width=0.5\columnwidth]{FigS2.pdf}
\caption{(a) Dynamics of total dimer action variable $M$ for two representative values of $u$ below and above the critical $u_c$ for $J=0.05$. The corresponding $u_c$ is $0.8$ as marked by the vertical dotted line in (b), where the long time value $M_{\infty}$ is plotted against $u$. $M_0$ denotes the the total $M$ at $t=0$. (c-d) Final distribution of initially prepared maximally excited array of dimers in the $\varepsilon, n$ plane for $J=0.05$. Below $u_c$ the nearly uniform state remains stable, as seen in (c), while above $u_c$ it evolves into a dispersed ensemble of many different single-dimer configurations, scattered over the array.}
\label{joson_stability}
\end{figure}

    
\end{document}


