Series expansion studies of the $J_1$-$J_2$-Heisenberg bilayer

TL;DR

This study uses series expansions to analyze the phase transitions and stability of quantum spin liquids in a frustrated bilayer Heisenberg model, revealing potential intermediate phases near critical frustration levels.

Contribution

It applies perturbative Continuous Unitary Transformation to determine the spectrum and phase boundaries, providing new insights into the stability of quantum spin liquids against interlayer coupling.

Findings

Quantum critical lines for magnetic transitions identified.

Well-defined (anti)bound states in the two-triplon sector found.

Evidence for bound state condensation near critical frustration.

Abstract

We study a bilayer of the frustrated $J_1$-$J_2$ Heisenberg-model on the square lattice. Starting from the dimer limit at strong interlayer coupling, we perform series expansions using the perturbative Continuous Unitary Transformation, based on the flow equation method, in order to determine the spectrum up to the two-triplon sector. From the one-triplon dispersion we obtain quantum critical lines for transitions from the dimer phase into either N\'eel or collinear magnetic order. For low to intermediate frustration these transitions are consistent with existing findings, based on the magnetic phases. In the region of strongest frustration, i.e. $J_2/J_1\sim 0.5$, we provide an estimate for the stability of the anticipated single-layer quantum spin-liquids against finite interlayer coupling. In the two-triplon sector we find a set of well defined (anti)bound states, which can be classified according to total spin and in-plane rotational symmetry. For vanishing frustration these states agree with previous series expansion analysis. For $J_2/J_1\gtrsim 0.5$ we provide evidence for a close-by condensation of one-triplon and two-triplon singlet bound states, suggesting that between the dimer and the collinear state, additional phases may intervene.

Equations

$$H$$

$$H_0$$

$$H_I$$

$$|t_{+1}⟩$$

$$|t_0⟩$$

$$|t_{-1}⟩$$

$$H_{\text{eff}}=H_0+\sum _{l,m}^{\infty }{C}_{l,m}{J}_1^l{J}_2^m,$$

$$|{\Psi }_n(\mathbf{K})⟩=\sum _{\mathbf{d}}{C}_{\mathbf{d}}|\mathbf{K},\mathbf{d},S_n^{},S_n^z⟩,$$

$$f(x)=\sum _{n=0}^{\infty }{a}_n{x}^n$$

$$f(x)=\frac{1}{2}\sum _{n=0}^{\infty }\frac{1}{2^n}\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)a_k{x}^k,$$

$$E(J_1,J_2)$$

$$E(J_1,\kappa )$$

$$a_n(\kappa )$$

$$E_{\text{Euler}}(J_1,\kappa )=\frac{1}{2}\sum _{n=0}^{N_{\text{max}}}\frac{1}{2^n}\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)a_k(\kappa )J_1^k.$$

Full text

I Introduction

Phases and excitations of quantum magnets are among the keys to understand correlated electron systems [1]. Stepping beyond conventional long-range magnetic order, exchange frustration is a prime ingredient to achieve novel states of matter in these magnets, displaying, e.g., spin liquid behavior, topological order, and exotic excitations [2, 3, 4, 5, 6, 7, 8]. In this context, the planar antiferromagnetic $J_{1}$-$J_{2}$ Heisenberg-model on the square lattice (J1J2HM) [9] is one of the pillars of frustrated quantum magnetism. While first analysis of this model dates back several decades, basic properties, like parts of the quantum phase diagram still remain open issues. Classically, the ground state is a Néel state for $\kappa=J_{2}/J_{1}<1/2$, and comprises two inter-penetrating Néel states with $\sqrt{2}\times\sqrt{2}$ structure for $\kappa=J_{2}/J_{1}>1/2$. The relative degeneracy of the latter Néel vectors is lifted by ”order-by-disorder” [10, 11, 12], leading to columnar order with an Ising, i.e. $Z_{2}$, symmetric order parameter of the pitch vector at $(0,\pi)$ and $(\pi,0)$. At $\kappa=1/2$ the Luttinger-Tzia method [13] results in macroscopic degeneracy of the classical ground state due to line minima of the energy versus the ordering pitch vector. Moreover, the leading order zero-temperature $1/S$ corrections to the order parameter diverge at $\kappa=1/2$ [9].

Turning to the quantum limit, i.e., $S=1/2$, solid evidence has been gathered, that for $0.5\lesssim k\lesssim 0.6$ no magnetic order exists. This has been collected from a plethora of approaches, some of which include exact diagonalizations [14, 15, 16, 17], series expansions [18, 19, 20], coupled-cluster theory [21, 22, 23], variational methods [24, 25], density-matrix renormalization group (DMRG) [26, 27, 28], (infinite) projected entangled-pair state ((i)PEPS) [29, 30, 31], functional-renormalization group (fRG) [32, 33], perturbative analysis [34, 35], and variational Monte Carlo (VMC) calculations [36, 37, 38, 39]. Instead of magnetic order, these approaches have predicted various novel ground states, including plaquette valence-bond crystals (PVBC), columnar valence-bond crystals (VBC), and quantum spin liquids (QSL), with and without a spin gap, however no consensus has been reached.

Materials, which may be proximate to the J1J2HM include Ba2CuWO6 [40], Sr2CuMoO6 [41], Sr2CuWO6 [41, 42], Sr2CuTeO6 [43], and Li2VO(Si,Ge)O4 [44, 45]. These materials cover a wide range of $\kappa$-values, realizing both Néel- and collinear ordered states. Unfortunately, a system in the most frustrated region, $\kappa\sim 0.5$, is still missing.

Apart from investigating the single layer case, it is of relevance, to extend the parameter space of spin systems by introducing further interactions. These can introduce additional well defined quantum phases, the connection of which to the single layer case can provide for more insight. Most popular along this line is the replication of a spin system in terms of antiferromagnetic dimers, forming, e.g., ladders [46], bilayers [47], and three-dimensional networks [48]. For strong dimer exchange $J_{\perp}$, these systems display a near product-state of weakly coupled singlets, the quantum dimer (QDM) phase, which features massive triplet excitations (triplons). While in some cases the reduction of the dimer exchange may lead to condensation of triplons into sought-for quantum phases of some non-dimerized original model, bilayer systems host their own unique set of physics. For this reason, a variety of bilayer systems have previously been studied under various objectives, e.g., uncovering rich phase diagrams [49, 50, 51, 52, 53], examining the crossing to 3D bulk materials [54], investigating effects of disorder [55] and hole doping [56] or analyzing emergent bound states [57] and topological excitations [58]. On the material side, an extensive number of systems exist, which are related to this theme, including Li2VOSiO4 [45], BaCuSi2O6 [59, 60], TlCuCl3 [61], Ba3Mn2O8 [62], and SrCu2(BO3)2[63].

The dimer version of the J1J2HM, forming an AA-stacked bilayer (J1J2BHM), has been considered by modified spin-wave theory [64], which results in magnetic order for all $\kappa$ as $J_{\perp}\rightarrow 0$, by dimer series expansion [65] for the spin-gap and the staggered susceptibility, and recently by application of the high-order coupled cluster method (CCM) [66] for the magnetization. The latter two studies find a paramagnetic region for $J_{\perp}\rightarrow 0$ in the range $0.45\lesssim\kappa\lesssim 0.65$ and $0.43\lesssim\kappa\lesssim 0.61$, respectively. This is consistent with studies of the single layer J1J2HM.

Dynamical properties, e.g., one- and two-triplon excitations in the QDM phase of the J1J2BHM, as well as their condensation into the ordered phases versus $\kappa$ and $J_{\perp}$ remain open issues. This provides the main motivation for our work. We will analyze the spectrum of the J1J2BHM up to the two-triplon sector, starting from the limit of decoupled dimers. We will use the perturbative Continuous Unitary Transformation (pCUT) [67], based on the flow equation method [68], in order to perform a series expansion for the excitation energies directly in the thermodynamic limit. pCUT has been applied successfully to a large variety of dimerized and n-merized quantum spin systems, including, but not limited to ladders [69], tubes [70], planar pyrochlores [71], various SU(2)-invariant Heisenberg bilayers [49, 55], as well as to Kitaev bilayers [51, 57].

The paper is organized as follows: In Sec. II the J1J2BHM is described. Sec. III provides for a general explanation of the pCUT method. Sec. IV details our results, i.e., the one-triplon excitations in Sec. IV.1 and two-triplon excitations (Sec. IV.2). Sec. V concludes our work and lists some speculations. A technical appendix, Sec. A , on specifics of a resummation method we use is included.

II Model

The Hamiltonian of the J1J2BHM reads

[TABLE]

where $\vec{S}=\{S^{\alpha}\}$ with $\alpha=x,y,z$ are spin-$\nicefrac{{1}}{{2}}$ operators, $J_{\perp}$ and $J_{1(2)}$ are the Heisenberg interlayer and (next-) nearest neighbor intralayer exchange, respectively. $L=1,2$ labels the two layers, $\mathbf{r}^{(\prime)}$ the sites of the square lattice and $\langle\mathbf{r},\mathbf{r}^{\prime}\rangle$, $\langle\langle\mathbf{r},\mathbf{r}^{\prime}\rangle\rangle$ denote NN and NNN sites. Fig. 1 shows a depiction of the sites and spin exchanges.

In general the couplings can be any combination of ferro- or antiferromagnetic interactions. In this study we only focus on the pure antiferromagnetic case, i.e., $J_{\perp},J_{1},J_{2}>0$. Further $J_{\perp}\equiv 1$ is chosen from here onwards to fix the energy scale and we will use the parameter $\kappa=\frac{J_{2}}{J_{1}}$ to measure the strength of the NNN-interactions.

III Method

The J1J2BHM under study features at least three limiting quantum phases which are adiabatically disjoint. First for $J_{\perp}\gg J_{1},J_{2}$ the system can be viewed as weakly interacting antiferromagnetic dimers and serves as the starting point for our studies. Second, in the case $J_{1}\gg J_{\perp},J_{2}$, Néel order on each layer is present. In a similar manner, for $J_{2}\gg J_{\perp},J_{1}$, each layer realizes a collinear magnetic order, i.e., a Néel order on each bipartite sublattice, while the relative angle between the spins on both sublattices is fixed through order-by-disorder selection.

We study the J1J2BHM starting from the limit $J_{\perp}\gg J_{1},J_{2}$, following the same evaluation scheme as in earlier work on the Kitaev-Heisenberg-Bilayer [51, 57]. For completeness we reiterate the main points:

The non-degenerate ground state of the model is formed as a product state of singlets on each dimer. The corresponding elementary excitations, i.e., triplets excited on single dimers, can be classified by their $S^{z}$-component. We write

[TABLE]

and refer to them as $\alpha=+1,0,-1$-triplons hereafter. These properties allow us to apply the perturbative Continuous Unitary Transformation (pCUT) technique [67], based on the flow-equation method [68], to the model, requiring the unperturbed Hamiltonian $H_{0}$ to have a non-degenerate ground state and an equidistant spectrum. In the present case each energy level of the spectrum of $H_{0}$ can be assigned a particle number $Q\geq 0$, i.e., the number of excited triplets, which describes the energy of the unperturbed states, i.e., $H_{0}~{}=~{}J_{\perp}Q\,+\,\text{const.}$ and especially $[H_{0},Q]=0$. $Q=0$ refers to the product ground state $\ket{}=\prod_{\mathbf{r}}\ket{s_{\mathbf{r}}}$ of singlets, while the one- and two-triplon states, i.e., $Q=1$ and $Q=2$, are writen as $\ket{\mathbf{r}\alpha}=\ket{t_{\mathbf{r},\alpha}}\otimes\prod_{\mathbf{r^{\prime}}\neq\mathbf{r}}\ket{s_{\mathbf{r}^{\prime}}}$ and $\ket{\mathbf{r}\alpha,\mathbf{r}^{\prime}\beta}=\ket{t_{\mathbf{r},\alpha}}\otimes\ket{t_{\mathbf{r}^{\prime},\beta}}\otimes\prod_{\mathbf{r}^{\prime\prime}\neq\mathbf{r},\mathbf{r}^{\prime}}\ket{s_{\mathbf{r}^{\prime\prime}}}$, respectively.

The perturbation $H_{I}$ of the Hamiltonian mixes different $Q$-sectors through creating or destructing triplons. By virtue of pCUT, the full Hamiltonian $H$ is transformed to an effective Hamiltonian $H_{\text{eff}}=UHU^{\dagger}$, which is $Q$-diagonal and can be expressed by a series in the perturbation parameters $J_{1(2)}$ as

[TABLE]

where $C_{l,m}$ are weighted products of terms in $H_{I}$, each comprising $l+m$ non-local creations(destructions) of triplons which in total conserve the $Q$-number. The weights of the $C_{l,m}$ are integer fractions, which are determined analytically, and independent of the specific model at hand, by recursive differential equations [67]. Due to transforming the Hamiltonian as a whole, pCUT works directly in the thermodynamic limit, and due to its perturbative nature, the expansion is exact up to the order calculated. Thus its results are well controlled for small parameters $J_{1(2)}$.

Using the $Q$-number conservation, we evaluate the spectrum by treating each sector independently. For this we determine the irreducible matrix elements of $H_{\text{eff}}$ for each value of $Q=0,1,2$ and solve the corresponding zero-, one- and two-particle problems, described in the following:

Due to the uniqueness of the ground state, the $Q=0$-case is described by a single matrix element, directly equaling the ground state energy $E_{0}=\bra{}H_{\text{eff}}\ket{}$.

The one-particle case is described by a translational-invariant matrix, leading to the one-particle dispersion $\mathbf{E}_{\mathbf{k},\alpha\beta}=\sum_{\mathbf{r}}\mathrm{e}^{i\mathbf{r}\mathbf{k}}\bra{\mathrm{r}\alpha}H_{\text{eff}}\ket{\mathbf{0}\beta}-\delta_{\mathbf{r},\mathbf{0}}\delta_{\alpha\beta}E_{0}^{\text{cl}}$, where $\alpha,\beta\in\{+1,0,-1\}$, $\mathbf{k}$ is a wavevector and $E_{0}^{\text{cl}}$ is the ground state energy calculated for the same cluster as the corresponding one-particle matrix element [67, 72]. In general this is a $3\times 3$-matrix, however, due to the $SU(2)$-symmetry of the model, the total $z$-component of the spins has to be conserved and thus different triplon flavors do not mix in the one-particle sector. In fact, they must have identical dispersions, we write $\mathbf{E}_{\mathbf{k},\alpha\beta}\equiv E(\mathbf{k})\delta_{\alpha\beta}$ and only refer to the dispersion $E(\mathbf{k})$ in the following. To check at least part of our series coefficients, we can compare the case $J_{2}=0$, i.e., the non-frustrated version of the bilayer, to previous results from the literature [73].

The two-particle problem, i.e., $Q=2$, is more challenging [72, 74]. Here the matrix elements $\bra{\mathbf{r}^{\prime}\alpha^{\prime},\mathbf{r}^{\prime}+\mathbf{d}^{\prime}\beta^{\prime}}H_{\text{eff}}\ket{\mathbf{r}\alpha,\mathbf{r}+\mathbf{d}\beta}$ describe two particles with initial(final) positions $\mathbf{r}^{(\prime)}$ and $\mathbf{r}^{(\prime)}+\mathbf{d}^{(\prime)}$ and triplon flavors $\alpha^{(\prime)}$, $\beta^{(\prime)}$. Similar to the one-particle case, an effective two-triplon Hamiltonian matrix $h_{\mathbf{K}}(\mathbf{d},\mathbf{d}^{\prime},\alpha\beta,\alpha^{\prime}\beta^{\prime})$ with respect to states $\ket{\mathbf{K},\mathbf{d},\alpha\beta}$ can be constructed, where $\mathbf{K}$ is the total momentum and $\mathbf{d}$ labels the two-triplon separation. This matrix $h_{\mathbf{K}}$, comprising of the analytical matrix elements of the SE, has a particular structure, directly representing the underlying physics, i.e., the two-particle scattering problem, that can be used to extract the two-triplon spectrum. Most important, $h_{\mathbf{K}}$ is band-diagonal with respect to $\mathbf{d}$, due the model only involving local spin interactions, and comprises of two different types of matrix elements: First, for $|\mathbf{d}^{(\prime)}|<d_{I}$, i.e., a two-triplon separation smaller then some characteristic length $d_{I}$, both triplons interact through an effective coupling determined by pCUT and the corresponding matrix elements describe those irreducible two-triplon interactions. Second, for larger $\mathbf{d}^{(\prime)}$, the triplons are too far separated to interact, but can still move separately across the lattice. This forms a semi-infinite band in $h_{\mathbf{K}}$, describing the propagation of scattering states. The latter allows us to diagonalize $h_{\mathbf{K}}$ numerically on sufficiently large lattices with periodic boundary conditions without neglecting any two-particle interactions. Thus the resulting spectrum will capture all relevant two-triplon states, especially any (anti-)bound states outside the two-particle continuum. In practice we evaluate the spectrum of $h_{\mathbf{K}}$ on a system with $20\times 20$-dimers (800 spins).

Again, the $SU(2)$-invariance of the model fixes the total $z$-component of a two-triplon state, resulting in a block-diagonal form of $h_{\mathbf{K}}$ under the constraint $\alpha+\beta=\alpha^{\prime}+\beta^{\prime}$. More over, the total spin $S=0,1,2$ of the two-triplon state is conserved as well, providing us with a suitable classification for the spin-structure of the two-triplon states, we write $\ket{\mathbf{K},\mathbf{d},\alpha\beta}\to\ket{\mathbf{K},\mathbf{d},S,S^{z}}$. In practice, this is used to simplify the evaluation of the SE, while also providing an additional check for the resulting states after the numerical diagonalization.

All evaluations of required matrix elements of $H_{\text{eff}}$ can be carried out on suitable chosen linked cluster graphs of the lattice. A detailed description of their construction procedure can be found in Ref. [57].

IV Results

In this section we describe our findings on the low energy spectrum of the J1J2BHM. Sec. IV.1 covers one-particle excitations, while Sec. IV.2 contains the results for the two-particle states. In both cases, we investigate the structure of the spectrum as well as the wave function of the excitations and determine a presumed outline of the phasediagram.

IV.1 One-particle excitations

Fig. 2 shows the one-triplon dispersion $E(\mathbf{k})$ along some symmetry lines of the Brillouin zone (BZ) for different values of $\kappa=\frac{J_{2}}{J_{1}}$ at a fixed value $J_{1}=0.1$, calculated to 7th order in the series expansion using pCUT. From this, the overall tendencies of the model can already be visualized:

First, for $\kappa\ll 0.5$, the dispersion has a clear minimum at $\mathbf{k}=(\pm\pi,\pm\pi)$, which coincides directly with the formation of Néel order in the single layer model for small NNN coupling $J_{2}$. Second, for $\kappa\gg 0.5$, the global minimum can be found at $\mathbf{k}=(\pi,0)$ (and its equivalent points). This fits the expectation from the single layer model, that for sizeable $J_{2}$ a collinear magnetic order ground state is formed. In both cases the minima remain at their respective $\mathbf{k}$-points even for higher $J_{1}$ to the point when the excitation gap vanishes. Last, for $\kappa\approx 0.5$, the dispersion lacks a clear minimum. Instead, almost all wave vectors along the edge of the BZ acquire nearly identical energies. This is consistent with the behavior of the single layer model, where at the maximally frustrated point, a line-degeneracy at zero energy of the dispersion indicates the absence of a well defined ordering vector. In the present case this scenario is accompanied by the dimer gap.

Based on those observations, it is natural to investigate the excitation gap at the most likely critical points $\mathbf{k}=(\pi,\pi)$ and $\mathbf{k}=(\pi,0)$ with respect to the interaction strengths and varying orders of the series expansion. Fig. 3 (a) displays the triplon energy at $\mathbf{k}=(\pi,\pi)$ in the Néel-near parameter space for $\kappa=0.2$ and varying $J_{1}$. It shows a comparison between 3rd to 7th order bare series expansion and equivalent results after an Euler resummation scheme is applied (details can be found in appendix A). In this case, it appears that the convergence of the series is satisfying, even without the application of a resummation scheme for most orders investigated. Only the series at $O(6)$ is an outlier, diverging in the vicinity of the tentative critical point. The resummation corrects this and leads to a monotonic convergence of the critical point. We find $J_{1,c}\approx 0.64$ for $\kappa=0.2$ at our highest expansion order using the resummation.

For larger $\kappa$ and at $\mathbf{k}=(\pi,0)$ the situation differs, see Fig. 3 (b). Here, we observe an alternating behavior between even and odd orders of the bare series expansion, with odd orders providing a gap closure, while even orders diverge. Applying Euler resummation as for the case of $\kappa<0.5$, this issue can be resolved, again leading to a monotonic decrease of the critical coupling with increasing expansion order. Using resummation, we find $J_{1,c}\approx 0.84$ for $\kappa=0.8$ at the largest expansion order.

Next, we use such findings to determine the phase diagram for the J1J2BHM. Because our approach provides the complete triplon dispersion as power series in $J_{1}$ and $\kappa$, we can easily scan for the critical coupling $J_{1,c}(\kappa)$ at which the triplon energy vanishes for a given wavevector for any value of $\kappa$. I.e., similar to the preceding paragraphs, we now consider the gap closings at $\mathbf{k}=(\pi,\pi)$ (Néel-type) and $\mathbf{k}=(\pi,0)$ (collinear-type) versus $\kappa$ and plot the resulting phase boundaries for various orders in Fig. 4. Panel (a) shows the results of the bare series expansion, while panel (b) shows those for the resummed series. For the two pitch vectors considered, in both cases only those boundaries are shown for which the excitation gap closes first.

Several remarks are in order: First, the phase diagram is displayed in terms of the parameters $(\kappa,\frac{1}{J_{1}})$. This allows for direct comparison of our findings to those from CCM in Ref. [66]. The latter are marked by open dots and squares in the figure. In terms of this parameter space, the decoupled dimer limit corresponds to regions of large $\frac{1}{J_{1}}$, while the single layer model resides on the lower x-axis at $\frac{1}{J_{1}}=0$. Second, contrasting panel (a) against (b) the influence of the series resummation is apparent. While in (a), and in particular for even orders, there are large windows of $\kappa$ with no gap-closing, for (b) these windows shrink. Moreover, the evolution versus expansion order in (b) suggest a well behaved convergence, with $O(7)$ critical lines from resummation not too far from the infinite order limit. In contrast to the CCM [66], our results from the one-triplon gap predict a slightly more extended dimer phase for small $\kappa$ based on the bare series and for all $\kappa$ based on the resummed series. This variance likely stems from a bias imposed on the two methods by virtue of their opposite ”starting phases”, i.e., dimer (LRO) phase for the pCUT (CCM).

For all $\kappa$ and in particular in the region of maximal frustration near $\kappa\approx 0.5$ the critical lines of neither the bare nor the resummed series show a tendency to approach $\frac{1}{J_{1}}\to 0$ upon increasing the expansion order. Rather, as one can see from Fig. 4 (b), there is a clear tendency near $\kappa\approx 0.5$, to stabilize the critical line for single triplon gap-closure at some minimum finite value of $\frac{1}{J_{1}}\sim 0.65$. This is remarkably close to the termination of reentrant behavior of the Néel and collinear phases, observed by CCM [66]. From the latter, the non-magnetic and potentially spin-liquid regime of the single-layer model, for $0.45\lesssim\kappa\lesssim 0.59$, extends upwards, forming an ”hourglass” shaped region at finite $J_{\perp}$, visible in Fig. 4 (b). Combining this with the critical line from the resummed SE, it is very tempting to speculate that the single-layer QSL, anticipated on the line $\frac{1}{J_{1}}=0$, is confined to the lower part of this hourglass and terminates within its constriction. This is very reminiscent of a somewhat similar situation of a QSL surrounded by a QDM and two reentrant LRO phases in the frustrated honeycomb bilayer Heisenberg model [50].

IV.2 Two-particle excitations

In this section we focus on the two-triplon excitations. Fig. 5 shows their spectrum versus the total momentum $\mathbf{K}$ along selected high-symmetry paths in the BZ for different values of $J_{1}$ and $\kappa$. A few notes are in order. First, for each total momentum $\mathbf{K}$, the spectrum comprises of two parts, i.e., a continuum of states and potentially several discrete (anti-)bound states. The continuum is formed from all combinations of two one-triplon states with energies $E(\mathbf{k}_{1})$ and $E(\mathbf{k}_{2})$ and total momentum $\mathbf{K}=\mathbf{k}_{1}+\mathbf{k}_{2}$. The (anti-)bound states can in principle occur at any $\mathbf{K}$. Their formation and size of splitting from the continuum however, i.e., (anti)binding energy, depends strongly on the specifics of the two-triplon scattering potential and the total wavevector. In particular, for small $\kappa$, we find (anti-)bound states primarily in the vicinity of $\mathbf{K}=(\pm\pi,\pm\pi)$, where the two-triplon continuum is narrowest, see Fig. 5 (a). Satisfyingly, these states are consistent with previous results from a different type of series expansion performed at $\kappa=0$ only, i.e., the $J_{1}$-Heisenberg-Bilayer [75].

Increasing $\kappa$, the width of the continuum in the vicinity of the $K$-point increases and the dominant (anti)bound states start to move. For intermediate $\kappa$ at the lower edge of the region of strongest frustration, i.e., $\kappa\sim 0.5$ and Fig. 5 (b), well split-off, low-energy bound states can be found all along $\Gamma$-$M$. For even larger $\kappa$, i.e., at the upper edge of the most frustrated region, Fig. 5 (c), bound states at the $\Gamma$-point are lowest in energy. Most remarkably, here, we find these latter energies to be lower than all of those from the complete one- and two-triplon spectrum. This renders the two-particle bound states the low-energy elementary excitations of the system. We will elaborate on this later.

Two additional features are visible in Fig. 5 (c). First, for larger $\kappa$, small ripples at the boundary of the continuum occur. This is not due to the series being of insufficiently high order, rather it is a finite size effect of particular nature. I.e., for any finite size of the one-triplon momentum space, the actual total momentum of the lower edge of the two-triplon continuum may be off from the available momenta $\mathbf{k}_{1}+\mathbf{k}_{2}$, leading to the ripples visible. Second, and interestingly, a contour plot of the continuum two-triplon density of states (DOS) in Fig. 5 (c) suggests, that a qualitative understanding of the positions of the bound states can be drawn from a T-matrix type of argument. I.e., in regions of high two-particle DOS close to the lower edge of the continuum, along $\Gamma$-$M$, bound states are ’pushed’ away from the continuum, while closer to the $K$-point, large two-particle DOS is observed only further into the continuum, leaving the bound states very close to the continuum boundary.

Now we turn to the internal structure of the bound states uncovered. As described in Sec. III, the spin-component of the two-triplon eigenstates of $H_{\text{eff}}$ can be classified according to total spin $S=0,1,2$ and its $z$-component $S^{z}=-S,\cdots,+S$. We find that all bound states we obtain are $S=0$ and $S=1$ states, while the antibound states satisfy $S=2$. Apart from the spin-quantum number, the two-triplon wave function can be classified according to lattice harmonics, i.e., the ’angular’ momentum loosely speaking. For that purpose we consider the $\mathbf{d}$-dependence of the two-triplon wave function

[TABLE]

where the $C_{\mathbf{d}}$ are wave function amplitudes, which in principle are complex. However, at $\mathbf{K}=(0,0)$, and because all exponentials of type $\sim\exp(i{\bf K}(\dotsc))$ are unity, the effective Hamiltonian in each $S,S^{z},Q{=}2$-subspace is real and symmetric. In turn, the $C_{\mathbf{d}}$ are real numbers at that point. By construction, $C_{-\mathbf{d}}=C_{\mathbf{d}}$ due to the triplons being indistinguishable. This defines a basic symmetry for all spatial wave functions considered. Morover, we can now identify additional symmetries of the (anti)bound states, from the wave functions exemplified in Fig. 6 (b) for $\kappa=0.6$ at $\mathbf{K}=(0,0)$. Namely first, there are states with the largest $C_{\mathbf{d}}$ for $\mathbf{d}$ along NN $J_{1}$-bonds. These are the lowest energy bound states. At $\mathbf{K}=(0,0)$ these states are sign reversed along the two orthogonal lattice directions, e.g., $C_{(1,0)}=-C_{(0,1)}$, i.e., they are odd under rotation by $\pi/2$ (states ’1’ to ’6’ and ’B’, Fig. 6 (b)). Second, there are states with largest amplitude along the NNN $J_{2}$-bonds (state ’A’, Fig. 6 (b)). These states are invariant under rotation by $\pi/2$. Finally, states of type ’C’ in Fig. 6 (b) show a quasi-cylindrical symmetry. This seems unnatural for the underlying square lattice.

Fig. 6 (a) highlights the evolution of the energies of the various bound states versus $J_{1}$. The figure details, that for all $J_{1}<0.5$, the lowest lying bound states display ’B’-symmetry, i.e., they are odd under $\pi/2$ rotation and are tightly bound with a dominant NN-amplitude. Only for $J_{1}\gtrsim 0.5$ lowest-energy bound states of type ’C’ (see lower right corner of Fig. 6 (b)) emerges. We argue that the latter states are unphysical artifacts of our SE. First, as can be seen in Fig. 5 (c), their dispersion is unusually steep, with only very few $\mathbf{K}$-points outside of the continuum. Second, their wave function is rather extended with no clear-cut lattice symmetry present. In turn, we discard ’C’-type states.

Finally, in Fig. 7 we discuss the options for bound-state criticality, i.e., a gap-closure of the QDM phase comprising two-triplon bound states. To appreciate this, we first realize, that for $\kappa\lesssim 0.5$ where the QDM phase is expected to condense into Néel LRO, the ordering pitch vector is at the $K$-point. There, and as can be seen from Fig. 5 (a), the two-triplon bound states close to the $K$-point are at high energies. Comparing the latter with those at the $K$-point in Fig. 2 (a) it is clear that the breakdown of the QDM phase is driven by condensation of one-triplon states only. This is very different for $\kappa\gtrsim 0.5$, where the pitch vector for collinear LRO is at the $M$-point. Here, by comparing the gap-closures obtained from the bare SE for the one-triplon states with the lowest two-triplon bound states, i.e., ’B’-type, we are faced with the remarkable fact, that down to $\kappa\sim 0.67$ one-triplon, as well as $S=0$ bound-states condense almost simultaneously within a reasonable level of accuracy and with the bound-state transition slightly above the magnetic one. Below $\kappa\sim 0.67$, a proper separation of bound and continuum states from the $Q{=}2$-sector turns infeasible. However it seems very likely that the black dashed line in Fig. 7 remains close to the red line, until the critical pitch switches from collinear to Néel. Since ’B’-states have $S=0$, it is now very tempting to speculate that these results indicate an additional, intervening non-magnetic quantum phase between the QDM and collinear region, along its upper edge in the quantum phase diagram of Fig. 4 (b). Even more surprising, such a phase, separating QDM from magnetic spiral states has been found also in the frustrated honeycomb bilayer Heisenberg model [50], where the intervening phase displayed nematic character. Consolidating our present findings by resummation or higher-order SE remains an open question beyond this work.

V Conclusion

To conclude, using pCUT series expansion, we have investigated the elementary excitations of the antiferromagnetic frustrated $J_{1}$-$J_{2}$ bilayer spin-$\nicefrac{{1}}{{2}}$ Heisenberg model on the square lattice in its quantum dimer regime, with a particular focus on the stability of the dimer phase as well as on the spectrum of (anti)bound two-triplon states. We have determined the quantum critical lines for a condensation of the one-triplon excitation into magnetic phases. The location of these transitions is in good agreement with other published analysis, starting from the magnetic phases. For a sizeable pocket in the $J_{1}/J_{2}-J_{\perp}$ plane, situated above the region of the potential quantum spin-liquid of the single layer, we have provided an approximate upper bound for its stability against the formation of a dimer gap. Regarding the two-triplon spectrum we have uncovered a rich structure of collective (anti)bound states which we have classified according to their total spin and rotational symmetry. We found, that frustration impacts the continuum density of states such, that not only the binding energy, but in particular the location of the lowest lying bound states can be shifted strongly within the Brillouin zone. For large frustration, this leads to a scenario in which the one-triplon and the $S=0$ bound states condense very close to each other, with a slight preference for the latter. This may suggests that the transition from the dimer into the collinear magnetic state could involve an additional intermediate nonmagnetic phase. This perspective calls for additional analysis. For an experimental probe of our result in the two-triplon sector, future investigations of optical probes, including phonon-assisted magnetic absorption, as well as magnetic Raman scattering should be of interest.

Acknowledgements.

This work has been supported in part by the DFG through Project A02 of SFB 1143 (Project-Id 247310070). Work of W.B. has been supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. W.B. also acknowledges kind hospitality of the PSM, Dresden.

Appendix A Euler resummation

Here we briefly describe the resummation of a power series using a variation of the Euler resummation. An infinite power series in $x$ of the form

[TABLE]

can be rewritten as

[TABLE]

which is called the Euler resummation. We are working with finite power series, depending on two parameters $J_{1}$ and $J_{2}$, which can be written as

[TABLE]

Now $E(J_{1},\kappa)$ is a polynomial in $J_{1}$ and the Euler resummation from Eq. (6) can be applied by limiting its $n$-sum to the same order $N_{\text{max}}$ as the bare series, yielding a resummed series

[TABLE]

In practice this limits the impact of the highest order of the series expansion, counteracting an alternating behavior between orders (see Fig. 3 (b)). This can improve the convergence.