Here we briefly introduce the ADT formalism. See Ref. [24] for a detailed discussion. In ADT, in order to have a complete description of the dynamics of a quantum system, a complete operator basis set (COBS)[25, 26], which we denote as $\mathscr{U}=\{\hat{u}^{\alpha}\}$ which satisfies the orthogonality and completeness conditions[26]

	$\displaystyle\text{Tr}((\hat{u}^{\alpha})^{\dagger}\hat{u}^{\beta})=C~{}\delta_{\alpha\beta},$		(1)
	$\displaystyle\hat{u}^{\alpha}\hat{u}^{\beta}=\sum_{\gamma}a^{\alpha\beta}_{\gamma}\hat{u}^{\gamma},$		(2)

where $C$ is a normalization factor which will be taken as 1 for convenience unless noted otherwise and $a_{\alpha\beta}^{\gamma}\in\mathbb{C}$. Satisfying both the complete and the orthogonal conditions, each operator $\hat{O}$ of $\mathscr{H}$ can be expanded into a sum of $\hat{u}^{\alpha}$s,

$\displaystyle\hat{O}=\sum_{\alpha}\text{Tr}(\hat{O}\hat{u}^{\alpha})\hat{u}^{\alpha}.$

(3)

It is convenient and conventional to decompose the algebras into a symmetric (bosonic) sector and an anti-symmetric (fermionic) sector

$\displaystyle[\hat{u}^{\alpha},\hat{u}^{\beta}]=\sum_{\gamma}b^{\alpha\beta}_{\gamma}\hat{u}^{\gamma},\quad\{\hat{u}^{\alpha},\hat{u}^{\beta}\}=\sum_{l}f^{\alpha\beta}_{\gamma}\hat{u}^{\gamma},$

(4)

with $b^{\alpha\beta}_{\gamma}=a^{\alpha\beta}_{\gamma}-a^{\beta\alpha}_{\gamma},\quad f^{\alpha\beta}_{\gamma}=a^{\alpha\beta}_{\gamma}+a^{\beta\alpha}_{\gamma}$.

Given a COBS $\mathscr{U}$, instead of specifying an arbitrary state as superposition of eigenstates or unit vectors of $\mathscr{H}$, the state can be specified by the expectation values of the bases operators $|\{\mbox{$\langle\hat{u}^{\alpha}\rangle$}\}\rangle$ (as a vector), or simply as $\langle\mathscr{U}\rangle$. In particular, according to Eq. (4), the density matrix, pure or mixed, can be written as $\rho=\sum_{\alpha}\text{Tr}(\rho\hat{u}^{\alpha})\hat{u}^{\alpha}=\sum_{\alpha}\mbox{$\langle\hat{u}^{\alpha}\rangle$}\hat{u}^{\alpha}$.

Given a COBS $\mathscr{U}$ and a state $\langle\mathscr{U}\rangle$, ADT further considers the complete set of (two-time) dynamical correlation functions (CSDCF) $\mathscr{G}$, whose elements $G^{\alpha\beta}_{\pm}[\text{\bf{i}},\text{\bf{f}}]$’s are defined as

$\displaystyle iG^{\alpha;\beta}_{\pm}[\text{\bf{i}},\text{\bf{f}}]=\mbox{$\langle\langle\mathcal{T}_{\pm}\Big{(}\hat{u}^{\alpha}(t_{i}),\hat{u}^{\beta}(t_{f})\Big{)}\rangle\rangle$},$

(5)

where $\mathcal{T}_{\pm}$ denotes time-ordering with the $\pm$ sign, $\langle\langle~{}\rangle\rangle$ denotes fully dynamical correlations defined as

$\displaystyle\mbox{$\langle\langle\hat{u}^{\alpha}(t_{i})\hat{u}^{\beta}(t_{f})\rangle\rangle$}=\left\langle(\hat{u}^{\alpha}(t_{i})-\mbox{$\langle\hat{u}^{\alpha}\rangle$})(\hat{u}^{\beta}(t_{f})-\mbox{$\langle\hat{u}^{\beta}\rangle$})\right\rangle.$

(6)

The equal-time-limit of $\mathscr{G}$ is a complete mapping $\mbox{$\langle\mathscr{U}\times\mathscr{U}\rangle$}\mapsto\mbox{$\langle\mathscr{U}\rangle$}$

	$\displaystyle iG^{\alpha;\beta}_{+}[\text{\bf{i}},\text{\bf{i}}]=\mbox{$\langle\{\hat{u}^{\alpha},\hat{u}^{\beta}\}\rangle$}/2-\mbox{$\langle\hat{u}^{\alpha}\rangle$}\mbox{$\langle\hat{u}^{\beta}\rangle$},$		(7)
	$\displaystyle iG^{\alpha;\beta}_{-}[\text{\bf{i}},\text{\bf{i}}]=\mbox{$\langle[\hat{u}^{\alpha},\hat{u}^{\beta}]\rangle$}/2,$		(8)

as $\mathscr{G}$ covers the full algebra table and a convention of $\theta(0)=1/2$ is adopted.

Next, we take a time-derivative on all the elements of CSDCF and apply the HEOM to obtain a many-body analogue of Schwinger-Dyson-equations-of-motion (SDEOM):

$\displaystyle\begin{split}i\partial_{t_{i}}G_{\pm}^{\alpha;\beta}[\text{\bf{i}},\text{\bf{f}}]&=\delta^{\alpha\beta}_{\pm}[\text{\bf{i}},\text{\bf{f}}]+i\mbox{$\langle\langle\mathcal{T}_{\pm}\big{(}[\hat{u}^{\alpha}(t_{i}),\mathcal{H}],\hat{u}^{\beta}(t_{f})\big{)}\rangle\rangle$}\\ &=\delta^{\alpha\beta}_{\pm}[\text{\bf{i}},\text{\bf{f}}]+\sum_{\eta\gamma}h_{\eta}b^{\alpha\eta}_{\gamma}G_{\pm}^{\gamma;\beta}[\text{\bf{i}},\text{\bf{f}}],\end{split}$

(9)

where $\delta^{\alpha\beta}_{\pm}[\text{\bf{i}},\text{\bf{f}}]=2\delta(t_{i}-t_{f})iG^{\alpha;\beta}_{\mp}[\text{\bf{i}},\text{\bf{i}}]$ and the Hamiltonian is given as $H=\sum_{\alpha}h^{\alpha}\hat{u}^{\alpha}$. Note that the equal-time-limit of SDEOM is also a complete mapping $i\partial_{t}\mbox{$\langle\mathscr{U}\times\mathscr{U}\rangle$}\mapsto\mbox{$\langle\mathscr{U}\rangle$}$.

For lattice problems, the COBS $\{\hat{u}^{\alpha\beta\cdots}_{ij\cdots}\}$ is typically constructed from local COBS $\{\hat{u}^{\alpha}_{i}\}$ as Cartesian products or their superposition $\hat{u}^{\alpha\beta\dots}_{ij\dots}=\hat{u}^{\alpha}_{i}\hat{u}^{\beta}_{j}\cdots$. When we represent a state with such many-body COBS, it is more physical to use their cumulant $\mbox{$\langle\hat{u}^{\alpha}_{i}\hat{u}^{\beta}_{j}\dots\rangle$}_{c}$ instead of the plain expectation values.

Finally, we can arrange the CSDCF into a vector form, and rewrite the complete set of SDEOM of the CSDCF in a matrix form:

$\displaystyle i\partial_{t_{i}}[{\mathbf{G}}_{\pm}]=[\mathbf{\Delta}]_{\pm}[\text{\bf{i}},\text{\bf{f}}]+[[\mathbf{L}]]\cdot[{\mathbf{G}}_{\pm}],$

(10)

where we use $[\mathbf{\Delta}]_{\pm}$ as the vector form of $\{\delta^{\alpha\beta}_{\pm}[\text{\bf{i}},\text{\bf{f}}]\}$ (with ${\alpha}$ as the component index) and $[[\mathbf{L}]]$ as the matrix form of $\{\sum_{\eta}h_{\eta}b^{\alpha\eta}_{\gamma}\}_{\alpha\gamma}$ (with different $\alpha$ and $\gamma$ as the element indexes). For time-translational invariant problems, we can transform Eq. (11) into frequency space as

$\displaystyle\left(\omega I-[[\mathbf{L}]]\right)\cdot[{\mathbf{G}}_{\pm}]=[\mathbf{\Delta}]_{\pm}[\text{\bf{i}},\text{\bf{f}}],$

(11)

which can solved formally. In this form, perturbation methods and theories of linear systems [27, 28] can be introduced. The LHS of Eq. (12) is determined by the Hamiltonian, disregarding the state of concern. Therefore, the RHS determines the pivoting for Eq. (12) as a linear system. In other words, the nonzero terms of $[\mathbf{\Delta}]_{\pm}[\text{\bf{i}},\text{\bf{f}}]$, i.e. the state. While exact solutions are generally difficult, if not impossible, it is reasonable to expect that a tractable of number of components could be sufficient to describe a system statistical-mechanically.

For a single spin-1/2, COBS can be chosen as either $\{\hat{I},\hat{s}^{x},\hat{s}^{y},\hat{s}^{z}\}$ or $\{\hat{I},\hat{s}^{+},\hat{s}^{-},\hat{s}^{z}\}$. The latter is used more often for spin dynamics. We shall follow such convention, using the former to for wavefunctions or states and the latter for dynamics. Its algebra table can be decomposed into a conventional form as:

$\displaystyle[\hat{s}^{\alpha},\hat{s}^{\beta}]=i\epsilon^{\alpha\beta\gamma}\hat{s}^{\gamma},\quad\{\hat{s}^{+},\hat{s}^{-}\}=1,\{\hat{s}^{\pm},\hat{s}^{z}\}=0.$

(12)

An arbitrary spin state, conventionally denoted by a wavefunction $\mbox{$|\phi\rangle$}=(a+ib)\mbox{$|\uparrow\rangle$}+(c+id)\mbox{$|\downarrow\rangle$}$, is equivalently represent by COBS as $\{\mbox{$\langle\hat{s}^{x}\rangle$}=ac+bd,\mbox{$\langle\hat{s}^{y}\rangle$}=ad-bc,\mbox{$\langle\hat{s}^{z}\rangle$}=(a^{2}+b^{2}-c^{2}-d^{2})/2\}$. CSDCF of a single spin-1/2 $\{G^{\alpha;\beta}_{\pm}[\text{\bf{i}},\text{\bf{f}}]\}$ contains the usual single particle correlations studied in literature.

To compare with experiments, the dynamical structure factors (DSFs) are often required, which are the Fourier transform of the unordered two-time correlation functions $S^{\alpha\beta}(\bm{x},t)=\mbox{$\langle\hat{s}^{\alpha}_{\bm{x}}(t)\hat{s}^{\beta}_{0}(0)\rangle$}/(2\pi)$. We can express DSFs via $\mathscr{G}$ as

$\displaystyle\begin{split}&S^{\alpha\beta}(t)=-\frac{i}{4\pi}(G^{\alpha\beta;R}_{+}(t)-G^{\alpha\beta;A}_{+}(t)\\ &+G^{\alpha\beta;R}_{-}(t)-G^{\alpha\beta;A}_{-}(t)).\end{split}$

(13)

With the above expression, the space-time DSF can be obtained by integrating the correlation functions in frequency space. For the rest of this work, we shall work with the causal Green’s functions only, and simply inspect $S^{\alpha\beta}[\omega]\sim-i(G^{\alpha\beta}_{+}[\omega]+G^{\alpha\beta}_{-}[\omega])/(4\pi)$ without carrying out the subtraction $G^{R}-G^{A}$.

Although the two-spin problems can be solved easily by diagonalization, the dynamical correlation functions have not been studied in detail. The two-spin problems can play the role of a Hubbard atom in Hubbard models, which is still being studied, generating useful insights[29] for lattice model studies.

For quantum spin models, this “atom” can describe states including i) a single bond correlation in RVB VBS etc, ii) ground states of TFIMs for which only two-spin cumulant correlations are present and iii) exactly solvable models with fully dimerized ground states such as the Shastry-Sutherland lattice model[30], spin- Su-Schrieffer-Heeger (SSH) models[31, 32], etc.. Therefore, detailed studies of of such atomic limits are necessary ground work for dynamical perturbation studies on top of these states and models. Moreover, analysis of the exact solutions can help us understand the analytic properties of the dynamical correlation functions in different states and contribute to improving the approximations one can make for similar states on a lattice.

To provide a ADT solutions, we first specify a grade-2 COBS for the problem. Since there are only two site involved, grade-2 is also the largest possible grade for the problem. Consider two spin-1/2’s $\hat{s}_{a}$ and $\hat{s}_{b}$. We use the following grade-2 COBS construction:

$\displaystyle\begin{split}&\hat{S}^{\alpha}_{ab}=\hat{s}_{a}^{\alpha}+\hat{s}_{b}^{\alpha},~{}\hat{\eta}_{ab}^{\alpha}=\hat{s}_{a}^{\alpha}-\hat{s}_{b}^{\alpha},\\ &\hat{B}_{A,ab}^{\gamma}=2\varepsilon_{\alpha\beta\gamma}\hat{s}_{a}^{\alpha}\hat{s}_{b}^{\beta},\\ &\hat{B}_{S,ab}^{\gamma}=2\varepsilon_{\alpha\beta\gamma}^{2}\hat{s}_{a}^{\alpha}\hat{s}_{b}^{\beta},\\ &\hat{D}_{ab}^{\alpha}=2\hat{s}_{a}^{\alpha}\hat{s}_{b}^{\alpha},\end{split}$

(15)

where the Einstein summation notation over dummy indices is assumed.

Although a direct product construction is equally valid and more convenient to implement in programming, this construction provides more physical insight which we shall discuss later. This particular construction follows from the geometric algebra[33]. But for the discussion of dynamical correlation functions, we stick to the Cartesian COBS. Since $B^{\gamma}_{S(A)ba}=\pm B^{\gamma}_{S(A)ab}$, we drop the _ab index and denote $B_{S(A)}^{\gamma}=B^{\gamma}_{S(A)ab}$. So the grade-2 COBS is written as

$\displaystyle\mathscr{U}_{g2}=\{\hat{I},\hat{S}^{\alpha},\hat{\eta}^{\alpha},\hat{B}_{S}^{\alpha},\hat{B}_{A}^{\alpha},\hat{D}^{\alpha}\}.$

(16)

Any state should be described by $\langle\mathscr{U}_{g2}\rangle$, which has 16 real numbers, exceeding the 8 real parameters allowed by considering the 4-dimensional $\mathscr{H}_{ab}$. This reflects the algebraic constraints on the operators’ expectation values. While formally we choose $\mathscr{U}_{g2}$ as Eq. (15) for convenience of algebras, we shall use $\mathscr{U}^{\prime}_{g2}=\{\hat{I},\hat{s}^{\alpha}_{i},\hat{B}_{S}^{\alpha},\hat{B}_{A}^{\alpha},\hat{D}^{\alpha}\}$, i.e. use local single particle operators $\hat{s}^{\alpha}$ instead of the Fourier components, when discussing Green’s functions.

To accommodate Eq. (17) and to simplify notations, we switch to notations of Green’s functions in the form as $G_{\pm}[\hat{O}_{1}^{\alpha};\hat{O}_{2}^{\beta}](t\text{ or }\omega)$, where $\hat{O}^{\alpha}$ can be an element of Eq. (17) or a product operator such as $\hat{s}^{\alpha}_{i}\hat{s}^{\beta}_{j}$. We also omit $(t\text{ or }\omega)$ when there is no ambiguity. Since it is too cumbersome to display the complete form of Eq. (12), we shall only put down and discuss necessary terms, i.e. the pivoting equations, for a given state.

Conventionally, the single particle Green’s functions, i.e. $G_{\pm}[\hat{s}^{\alpha}_{i};\hat{s}^{\beta}_{j}]$, can be solved from SDEOM, if and only the vertex functions due to interactions can be solved. For example, the typical single particle GFs of concern are those of $\hat{s}_{i}^{\pm}$. With only $\mathcal{H}_{0}$, the HEOM reads

$\displaystyle i\partial_{t}\hat{s}^{\alpha}_{i}=\alpha(-J_{z}\hat{s}^{z}_{\bar{i}}\hat{s}^{\alpha}_{i}-h^{z}\hat{s}^{\alpha}_{i}),$

(18)

with $\alpha=\pm$, which leads to the following SDEOM

$\displaystyle\begin{split}&\alpha\omega G_{+}[\hat{s}^{\alpha}_{i};\hat{s}^{\beta}_{j}]=\delta_{ij}\delta_{\alpha\beta}\mbox{$\langle 2\hat{s}^{z}_{i}\rangle$}\\ &-h^{z}G_{+}[\hat{s}^{\alpha}_{i};\hat{s}^{\beta}_{j}]-J_{z}G_{+}[\hat{s}^{z}_{\bar{i}}\hat{s}^{\alpha}_{i};\hat{s}^{\beta}_{j}],\end{split}$

(19)

for $\alpha,~{}\beta,~{}\gamma\in\{+,-\}$ and $\alpha$ in the equations is interpreted as a sign $\pm$.

Conventional approaches rely on the decoupling the vertex function $\mbox{$\langle\hat{s}^{z}_{\bar{i}}\hat{s}^{\alpha}_{i};\hat{s}^{\beta}_{j}\rangle$}\rightarrow\mbox{$\langle\hat{s}^{z}_{\bar{i}}\rangle$}\mbox{$\langle\hat{s}^{\alpha}_{i};\hat{s}^{\beta}_{j}\rangle$}$, which is only exact for not-time-ordered correlation functions of a product state. Otherwise, the decoupling requires extra free parameters, which is termed “random phase approximation”(RPA) or Tyablikov-decoupling [19, 20]. As we shall find next, for entangled states and time-ordered correlation functions, the vertex functions are independent dynamical degrees of freedom. RPA and its parameters are of a dynamical origin.

Consider a state specified in Table 1, the single particle GFs and the associated vertex functions are

$\displaystyle\begin{split}&G_{+}[\hat{s}^{\alpha}_{a};\hat{s}^{\beta}_{a}]=\frac{\delta_{\alpha\bar{\beta}}\alpha(2\mbox{$\langle\hat{s}^{z}_{a}\rangle$}\omega-J_{z}/2)}{(\omega+J_{z}/2)(\omega-J_{z}/2)},\\ &G_{+}[\hat{s}^{\alpha}_{a};\hat{s}^{\beta}_{b}]=\frac{\delta_{\alpha\bar{\beta}}J_{z}\mbox{$\langle\{\hat{s}_{a}^{\alpha},\hat{s}_{b}^{\beta}\}\rangle$}}{(\omega+J_{z}/2)(\omega-J_{z}/2)},\\ &G_{+}[\hat{s}^{z}_{b}\hat{s}^{\alpha}_{a};\hat{s}^{\beta}_{a}]=\frac{\delta_{\alpha\bar{\beta}}(\mbox{$\langle\hat{s}^{z}_{a}\rangle$}J_{z}-\omega)/2}{(\omega+J_{z}/2)(\omega-J_{z}/2)},\\ &G_{-}[\hat{s}^{\alpha}_{a};\hat{s}^{\beta}_{a}]=\frac{\delta_{\alpha\bar{\beta}}(\omega-\alpha J_{z}\mbox{$\langle\hat{s}^{z}_{i}\rangle$})}{(\omega+J_{z}/2)(\omega-J_{z}/2)},\\ &G_{-}[\hat{s}^{\alpha}_{a};\hat{s}^{\beta}_{b}]=\frac{\delta_{\alpha\bar{\beta}}\mbox{$\langle\{\hat{s}_{a}^{\alpha},\hat{s}_{b}^{\beta}\}\rangle$}\omega}{(\omega+J_{z}/2)(\omega-J_{z}/2)},\\ &G_{-}[\hat{s}^{z}_{b}\hat{s}^{\alpha}_{a};\hat{s}^{\beta}_{a}]=\frac{\delta_{\alpha\bar{\beta}}\omega\mbox{$\langle\hat{s}^{z}_{a}\rangle$}}{(\omega+J_{z}/2)(\omega-J_{z}/2)},\\ \end{split}$

(20)

where $\bar{\pm}=\mp$. First, we verify the solutions by taking the equal-time limit, which should reduce to Eq. (8) and (8). While the above expressions are valid for an arbitrary state of Table 1, we shall discuss two of the most important ones: the fully polarized states and the fully entangled states.