Fractional and composite excitations of antiferromagnetic quantum spin trimer chains

Using quantum Monte Carlo, exact diagonalization and perturbation theory, we study the spectrum of the $S=1/2$ antiferromagnetic Heisenberg trimer chain by varying the ratio $g=J_2/J_1$ of the intertrimer and intratrimer coupling strengths. The doublet ground states of trimers form effective interacting $S=1/2$ degrees of freedom described by a Heisenberg chain. Therefore, the conventional two-spinon continuum of width $\propto J_1$ when $g=1$ evolves into to a similar continuum of width $\propto J_2$ when $g\to 0$. The intermediate-energy and high-energy modes are termed \emph{doublons} and \emph{quartons} which fractionalize with increasing $g$ to form the conventional spinon continuum. In particular, at $g \approx 0.716$, the gap between the low-energy spinon branch and the high-energy band with mixed doublons, quartons, and spinons closes. These features should be observable in inelastic neutron scattering experiments if a quasi-one-dimensional quantum magnet with the linear trimer structure and $J_2


I. INTRODUCTION
Many quasi one-dimensional (1D) magnetic materials with spin S = 1/2 moments harbor exotic phenomena originating from the physics of the Heisenberg antiferromagnetic chain (HAC) and its extensions [1]. The dynamic spin structure factor of KCuF 3 [2,3], measured using inelastic neutron scattering, exhibits the characteristic gapless two-spinon continuum [4] of the uniform HAC. A phase transition to a gapped dimerized state, driven by additional frustration and spin-phonon couplings has been realized in CuGeO 3 [5]. In systems with random couplings, the random-singlet state with infinite dynamic exponent forms [6,7], as originally observed in a class of Bechgaard salts [8,9] and more recently in BaCu 2 SiGeO 7 [10] and BaCu 2 (Si 0.5 Ge 0.5 ) 2 O 7 [11]. Furthermore, the resonant inelastic x-ray scattering (RIXS) technique has now enabled specific detection of multi-spinon excitations [12,13] in the HAC material Sr 2 CuO 3 , and string excitations have been identified by terahertz spectroscopy in the Heisenberg-Ising compounds SrCo 2 V 2 O 8 [14] and BaCo 2 V 2 O 8 [15].
A unit cell of more than one spin, which is the context of our work presented in this paper, can lead to an even richer variety of 1D magnetic properties [16][17][18]. For example, ladder systems with rungs consisting of an odd or even number of S = 1/2 spins have a gapless or gapped spectrum, respectively [16], in a way similar to Haldane's conjecture [19] of chains with half-odd-integer or integer spins. Trimer chains have also been studied experimentally, including A 3 Cu 3 (PO 4 ) 4 (A --Ca, Sr, Pb) [20][21][22][23] and (C 5 H 11 NO 2 ) 2 · 3 CuCl 2 · 2 H 2 O [24], where in all cases the structure is such that two spins in one trimer are coupled to two spins of their neighboring trimers. The linear trimer chain with repeated couplings J 1 -J 1 -J 2 (intratrimer J 1 and intertrimer J 2 ), is realized in FIG. 1. Schematic representation of a trimer spin chain. J 1 and J 2 represent the two different antiferromagnetic nearest-neighbor Heisenberg intratrimer and intertrimer couplings, respectively. We here consider systems with J 1 ≥ J 2 > 0. We will use the letters a, b, c as indicated to refer to the three spins within a given unit cell.
From the theoretical perspective, it is interesting to consider the linear trimer chain illustrated in Fig. 1, where for J 2 ≪ J 1 the excitations can be understood from perturbative calculations starting from the eigenstates of the isolated trimers. Surprisingly, though the case J 2 > J 1 has been studied both experimentally and theoretically, the potential of the J 1 > J 2 system (J 1 , J 2 both antiferromagnetic) to realize a host of interesting excitations and their confinement-deconfinement crossovers has not been recognized in the previous literature. We will study this model system extensively here.
An interesting preliminary aspect of the J 1 -J 1 -J 2 system, where we define g ≡ J 2 /J 1 ∈ (0, 1], is that it reduces to the conventional HAC when g = 1, while for g ≪ 1 the lowenergy excitations can be mapped onto an HAC consisting of one effective S = 1/2 spin in each unit cell. Thus, in both these limits the low-energy excitations should be spinons, but they live in different Brillouin Zones (BZs). By reducing the coupling ratio from g = 1 one can expect an evolution from the conventional two-spinon continuum in the window q ∈ [0, π] of width ∝ J 1 into three continua in the windows q ∈ [0, π/3], [π/3, 2π/3], [2π/3, π] with band width ∝ J 2 . Moreover, at small values of g there should also be weakly dispersive modes arising from the internal excitations of the trimers, and these must eventually fractionalize into spinons as g → 1. The interplay between the two types of spinons and the higher-energy modes, and how these coexisting excitations eventually evolve into just the conventional spinon continuum, is not immediately clear. Our calculations reported here are also in part motivated by recent work on twodimensional systems consisting of weakly coupled multi-spin plaquettes of different shapes [30,31]. In the calculations for coupled 3 × 3 plaquettes, intriguing spectral features were found but were not fully explained [30], and the simpler 1D system considered here can guide additional calculations and provide interpretations.
We will develop a physical understanding of the various observed branches of excitations and their intricate evolution with g by interpreting the numerical spectral functions in the light of perturbative calculations as well as the known properties of the spinon continuum in the uniform HAC. The S = 1 excited states of the uniform S = 1/2 HAC are fractionalized into independently propagating particles, spinons, each carrying S = 1/2 [32][33][34]. The leading contributions of these excitations (two-spinon contributions) with total momentum q to the dynamic structure factor S(q, ω) can be calculated relatively easily [4] by Bethe ansatz (BA) calculations, while four-spinon contributions require sophisticated numerical calculations with the BA states [35]. The utility of the BA is limited, however, when perturbing the HAC beyond the solvable XXZ model [35,36]. In general reliable calculations of dynamical properties are very challenging beyond the small lattices accessible to exact diagonalization (ED) techniques [1].
Currently, the density matrix renormalization group (DMRG) [38,39] and related methods formulated with matrix-product states (MPS) [40] are very powerful for 1D systems and are also applicable for calculations of dynamical structure factors [41,42], primarily for systems with open boundaries, due to inefficiency in the case of periodic chains. Quantum Monte Carlo (QMC) simulations with subsequent numerical analytic continuation of the imaginary-time correlation functions [43][44][45][46][47][48][49] can be applied to large periodic system sizes in any number of dimensions as long as the negative sign problem can be avoided. Some very useful results for S(q, ω) have been obtained for a variety of quantum magnets, e.g., in Refs. [48][49][50][51][52][53].
We here study the trimer chain using both the QMC and ED methods. For the former, we compute imaginary-time correlations with the stochastic series expansion [46] QMC method and employ a variant of the stochastic analytic continuation (SAC) technique [43][44][45][47][48][49]. We will discuss results for S(q, ω) in both the regular and reduced BZs. When g is small, three well separated features are observed. A low-energy continuum extending up to ω ∝ J 2 with a characteristic spinon structure is present due to the fact that each trimer hosts an effective spin-1/2 degree of freedom and an effective HAC with coupling ∝ J 2 forms. At higher energies, ω ∝ J 1 , there are two different weakly dispersing modes corresponding to the internal excitations of the trimers. These modes evolve significantly as g is increased, and eventually for g → 1 the standard spinon continuum of the uniform chain is recovered. In order to further understand the excitation mechanism and the properties of these excitations, we construct perturbatively motivated expressions for propagating internal trimer excitations and find excellent agreement with the numerical results for small to moderate values of g. According to the characters of these excitations, we call the quasiparticle corresponding to the intermediate-energy excitation (ω ≈ J 1 ) the doublon and the one corresponding to the higher energies (ω ≈ 1.5J 1 ) the quarton. These excitations may be helpful for understanding similar dual flat bands observed in the inelastic neutronscattering spectra of A 3 Cu 3 (PO 4 ) 4 (A --Ca, Sr, Pb) [20], where there are next-nearest-neighbor interactions but where trimers are also effectively weakly coupled as in the g ≪ 1 system studied here.
We will also use ED to compute the dynamic structure factor in a truncated Hilbert space, where the full low-energy spinon space is included in addition to one each of the doublon and quarton. For g 0.5, the good agreement of the spectral functions obtained in this manner with those from the other calculations confirm the nature of these excitations, while disagreements for larger g show how the doublons and quartons loose their identity when they begin to fractionalize into the standard HAC spinon continuum that emerges when g → 1.
In the intermediate g regime, we have two different spinon continua coexisting with the doublon and quarton. These calculations demonstrate how two branches of trimer excitations gradually broaden out when g increases, then merge together and evolve into the upper part of the two-spinon continuum in a fractionalization process. Based on these results, we also develop an intuitive picture of the quasi-particles as complex domain walls, generalizing the standard domain-wall description of spinons in the HAC. The high-energy spin excitations (energy of order J) have also been under intense scrutiny in the antiferromagnetic parent compounds of the high-T c superconductors [54][55][56][57][58] and other systems described by the 2D Heisenberg model [48,59]. Our results on the fractionalization mechanism of high-energy spin excitations of the trimer chain may be helpful for further exploring the fractionalization mechanism of high-energy magnons in these systems.

A. Model
The Hamiltonian of the spin-1/2 antiferromagnetic trimer chain with periodic boundary conditions reads where S i,α is the spin-1/2 operator at the αth site in the ith trimer, the intratrimer labels α ∈ {a, b, c} are explained by Fig. 1. The total number of trimers is N and the total length of system is L = 3N. The tuning parameter g is defined as g = J 2 /J 1 and we here limit our study to 0 ≤ g ≤ 1. For simplicity, we set the intratrimer interaction J 1 = 1 as the energy unit, so that intertrimer interaction J 2 = g. Our interest is in the whole range of intermediate coupling ratios g ∈ (0, 1) Dynamic spin structure factor S(q, ω) obtained from QMC-SAC calculations for different g values. The color coding of S(q, ω), illustrated with the bar on the right side, uses a piecewise function where the boundary value U 0 = 5 is indicated by the line on the color bar. Below the boundary, the low-intensity portion is characterized by a linear mapping of the spectral function to the color bar, while above the boundary a logarithmic scale is used, U = U 0 + log 10 [S(q, ω)] − log 10 (U 0 ). Since the SAC method generates spectra with intrinsic broadening, no additional broadening is imposed here. The vertical striped features, noticeable especially between the ω ≈ 1 and ω ≈ 1.5 bands in (a) and (b), are typical in analytic continuation when the statistical precision of the QMC imaginary-time data is barely sufficient for resolving two features, thus not completely resolving them for some momenta while resolving them in other cases (c)-(h).
where the system evolves between the isolated trimers and the isotropic BA solvable HAC. The dynamic structure factor describing the time dependent spin-spin correlations at a given transferred momentum q is defined as where β,γ refer to spin components x, y, z, and S β(γ) q is the Fourier transform of the spins that we discuss later (as it depends on the type of BZ used). In frequency space the dynamic structure factor is where we have explicitly indicated the diagonal z component, which already contains all information in the case of the spinrotational invariant model considered here. From now on we will omit the superscripts and define S(q, ω) ≡ 3S zz (q, ω).

B. Overview of numerical results
In this section, we present the spin excitation spectra of the trimer chain obtained by QMC-SAC calculations (see the technical details in Methods IV). Color plots of S(q, ω) representing the full (q, ω) space are shown in Fig. 2 for a chain with 192 spins. We have also calculated the full spectrum by ED for a chain with 30 spins (see Supplementary Figure 1), which exhibits similar features as the QMC-SAC results but with significant finite-size effects. To better reveal the important spectral features, here and in graphs in later sections, we used color coding of S(q, ω) by a piecewise function which is explained further in the caption of Fig. 2.
Let us first examine results for g = 1, shown in Fig. 2(h), where we observe the well understood characteristic asymmetric spinon continuum of the HAC chain. As q → 0, the width of the continuum vanishes along with the spectral weight, while at the other gapless point q = π the continuum has maximum width and the weight in the thermodynamic limit is ω −1 divergent (with a logarithmic correction). At the lower edge away from the gapless points, the dispersion relation ǫ q = (π/2)| sin(q)| reflects that of a single spinon. The spectral weight is also concentrated at this edge, with a (ω − ǫ q ) 1/2 divergence. The divergent features can of course not be strictly observed in the finite systems, and in the case of QMC-SAC results there is also broadening and some distortions due to the incomplete data used in the analytic continuation. Nevertheless, the lower spectral bound is well reproduced and the observed concentration of spectral weight at the lower edge is a true feature of the spinon continuum [4,35].
Next, we discuss how the spectral function evolves as g is increased from 0 and eventually approaches 1. As apparent in Fig. 2(a)-(c), when g = 0.1 − 0.3 the gapless low-energy spectrum comprises spinon excitations originating from an effective HAC of N effective S = 1/2 degrees of freedom. The reduced BZ corresponds to q ∈ [0, π/3] in the full BZ used in the figures. The same excitations appear also in the windows q ∈ [π/3, 2π/3] and q ∈ [2π/3, π], with different weight distributions due to the different phase factors. In addition to q = 0, π, the spectrum is gapless at q = π/3 and q = 2π/3 when g < 1, which can also be explained by the Lieb-Schultz-Mattis theorem [60,61] along with BZ folding effect since the system is rotationally invariant and translationally invari- ant with unit cell containing three spins.
Many of the observed features follow from the fact that the three low-energy spinon continua must evolve into a single continuum as g → 1. Thus, the spectral weight in the central portions q ∈ [π/3, 2π/3] of the low-energy spinon continuum decreases while the leftmost q ∈ [0, π/3] and rightmost q ∈ [2π/3, π] half arches become more prominent. A very interesting aspect of the evolution is how the intermediateenergy and high-energy modes gradually morph into the highenergy part of the standard HAC spinon continuum. Thus, a fractionalization of the quasiparticles takes place. In the following, we will provide a perturbative analysis to explain the evolution of the spectrum.

C. Perturbative analysis: Effective Heisenberg coupling
In Fig. 3 we observe that the low-energy doublets contain a singlet and an unpaired spin; thus, each trimer contains an effective spin-1/2 degree of freedom. Furthermore, there is a clear separation to the higher-energy states, which according to our results in the previous section survives at least for g ≤ 0.4. Therefore, the low-energy excitations of trimer chain can be well described by an effective Hamiltonian whose excitations only contain the spinons. The trimer chain is translationally invariant with a unit cell of three spins and rotationally invariant, and the effective Hamiltonian must conserve the translation symmetry and SU(2) symmetry. We explicitly derive the effective HAC arising from the doublet trimer ground states by applying the Kadanoff method [2][3][4][5], projecting the Hamiltonian onto the a low-energy subspace constructed from the lowest eigenstates of each trimer (see Supplementary Note 2). The effective Hamiltonian for N trimers which describes an isotropic HAC with an effective coupling strength J eff = 4J 2 /9, and S is the effective spin-1/2 operator. We next define a dynamic structure factor in the reduced BZ as follows where S zz αα (q, ω) is the dynamical structure factor as defined in Eq. (2) but including only the spins at location α ∈ {a, b, c}. The momenta are then of the form q = 6nπ/L for the system of L = 3N spins. The low-energy part is the twospinon continuum, with the predicted lower boundary ω l = πJ eff sin(q) /2 and upper boundary ω u = πJ eff sin(q/2) , is indicated in the reduced BZ in Figs. 4(a)-(d). The predicted lower boundary ω ′ l = πJ eff sin(3q) /2 and upper boundary ω ′ u = πJ eff sin(3q/2) in the full BZ is similarly shown in Figs. 4(e)-(f). We observe that the boundaries are in good agreement with the numerical results, and the spectral weight close to the upper bound is very small at q = π, as is well known in the case of the standard HAC. Thus, we have confirmed that the trimer chain reduces to an effective HAC with exchange interaction J eff = 4J 2 /9 even for g as high as about 0.4.

D. Perturbative analysis: Propagating internal trimer excitations
Looking at the full level spectrum and corresponding eigenvectors of one single trimer in Fig. 3, the ground state is a doublet with energy E 0 = −J 1 , total spin quantum number S = 1/2, and magnetic quantum number M = ±1/2. The first and second excited states of the trimer are a doublet with E 1 = 0, S = 1/2 and a quartet with E 2 = J 1 /2, S = 3/2, respectively. We also show the structure of the eigenstates when written as pairs of spins forming a singlet or a triplet, with one or three left over unpaired spin. The two low-energy doublets both contain only singlets and unpaired spins, while the higher quadruplet excitations contain either a triplet pair and an unpaired spin or three unpaired spins. The magnetic quantum number M = ±1/2 or M = ±3/2 corresponds to the number of unpaired spins in each case.
Taking two trimers as an example, when g is small enough the coupling can be regarded as a perturbation of the product state of the isolated trimers. There are four states forming a singlet ground state and a triplet excitation with a gap of order g. In addition, the lowest internal excitations of the trimers have a gap of order J 1 to the low-energy states. The almost degenerate excitations originating from the two different trimers include S = 1 states, which are of our primary interest when considering the dynamic spin structure factor. Increasing the number of trimers, the analogy of the lowest singlet and triplet of the two-trimer system is the spinon continuum, which for non-interacting spinons come with S = 0 and S = 1. The internal excitations of the trimers will form weakly dispersive bands for small g, and lose their single-trimer identity for larger g. For simplicity, the calculations involving these higher excitations will not be carried out with the spin rotation Dispersion relations of the intermediate-energy and high-energy excitations in the reduced BZ. All results are from the case where g = 0.1. The numbers marked on curves represent the degeneracies. In (a) each dispersion relation is two-fold degenerate. In (b) four colors have been used to distinguish the different degeneracies. symmetry maintained. Being interested in S = 1 excitations probed with the dynamic structure factor, we will still focus on excited states with |∆M| = 1.
The internal excitations of the trimer chain are almost localized when g is small, and cannot be classified as magnons or triplons. In order to understand the nature of these excitations, we propose a scheme to calculate their dispersion relations. The results are shown in Fig. 4 as a number of dispersion relations corresponding to propagation of internal trimer excitations in the background of a chain with defects mimicking the complex ground state of the effective HAC. The overall good agreement with the QMC-SAC results on the location of these excitations and their band widths suggest that our picture of the excitation is correct even though the calcula-tion involves a very rough approximation of the ground state. The bending of, in particular, the intermediate band, and the broadening for g 0.3 visible in the full BZ are not captured by the simple ansatz used here, which will be presented in detail below. The approach nevertheless forms a good starting point before considering other approaches.
By assuming the ground-state wave function of the spin chain is the product states of ground states of each trimer and the excited-state wave function contains one excited trimer, we are able to calculate the dispersion relations corresponding to the intermediate-energy excitations with |∆M| = 1. Considering the massively degenerate states obtained from the possible choices of the doublet ground state on each trimer and total magnetization M = 0, it is found that the dispersion relations are mainly dependent on the excited trimers and their neighbours. Translated to finite size, this means that the correct generic result is already obtained from a system with N = 4 trimers. The details can be found in Supplementary Note 3. The result is six remaining dispersion relations, where each case is two-fold degenerate. All three dispersion relations for g = 0.1 are graphed in Fig. 5(a). Two dispersion relations are independent of q since the excitations are localized. According to the characters of this excitation starting from the g = 0 limit, we refer to it as the doublon.
Calculating the dispersion of the high-energy excitations evolving from the trimer quartet at ω = 3J 1 /2 is more complicated but we follow the same strategy as for the doublon. The intermediate-energy band is formed by reconfiguring the singlet bond in the original ground state of the trimer in Fig. 3 at a cost of J 1 . The higher band is formed out of trimers excited into their S = 3/2 highest state, which can be seen as the singlet of two spins excited into a triplet, at a cost of 3J 1 /2. Accordingly, we refer to this high-energy quasiparticles as the quartons. This rough estimation of the quarton energy matches with the QMC-SAC results for small g values.
Our calculation does not conserve the total spin of the collective many-body state, but we consider the excitations corresponding to ∆S = 1 on one trimer. As a result, 96 dispersion relations fall into 33 different cases, which are displayed in Fig. 5(b) for g = 0.1. The calculation details can be found in Supplementary Note 3. The number marked on each curve shows the degeneracy of every dispersion relation of the high-energy excitation depending on the number of times it appears among the total 96 cases. Among these dispersion relations, some are independent of q, see Eq. (13), and the dispersion relations E 2 − E 0 − J 2 /9 and E 2 − E 0 both have the maximum degeneracy, eight. It can be found that most of the dispersion relations in Eq. (13) also have large degeneracies. The reason is that when g is small, these excitations are localized in the trimers and dominate the whole types of excitation in our perturbative calculation. When g is increased, the dispersion relations dependent of q will become more significant, therefore the deviation of our perturbative calculation on the condition of small g will be more obvious.

E. Comparisons of numerical and analytical results
Next, we compare the above doublon and quarton dispersion relations with QMC-SAC results in both the reduced and full BZs. At g = 0.1, Fig. 4(a), we clearly observe three different bands of excitations; in addition to the low-energy spinon continuum we have weakly dispersive intermediate-energy (ω ≈ J 1 = 1) and higher-energy bands (ω ≈ 2J 1 /3 = 1.5) exactly in the regions where the QMC-SAC results exhibit large spectral weight. Fig. 4(e) shows the unfolded dispersion relations and QMC-SAC results in the full BZ. The intermediateenergy and high-energy modes are more clearly visible in the full BZ, as the three q windows have different weighting for the a, b, and c trimer spin operators, thus offering more opportunities for an optimal weighting that makes any of the features visible. The advantages of the full BZ are even more clear at g = 0.2 in Figs. 4(b) and 4(f), where the dynamic structure factor in the reduced BZ does not exhibit two separate modes but they have merged into a single continuum. This merger of the two modes may partially be due to the limitations of the QMC-SAC approach, but most likely it reflects to a large extent the actual weight distribution. As g is increased, these two bands merge into each other also in the full BZ. It is remarkable how flat the bands are even at g as large as 0.4, and that the perturbative calculation at least gives the correct region of dominant spectral weight in the reduced BZ. However, from the full BZ results it is also clear that nonperturbative effects set in, with dispersive modes visible in between the two bands predicted by the variational states with only a single excited doublon or quarton. These dispersive modes that grow out of the doublons and quartons eventually evolve into the upper part of the spinon continuum as g → 1, as seen in Fig. (2).
We can also see in Fig. 2 that the low-energy and highenergy bands merge near g = 0.7. As shown in Figs. 6(a)(b), S(π/6, ω) and S(5π/6, ω) both exhibit two peaks when g = 0.7, while in both cases only one peak is present for g = 0.8. Along with Figs. 2(f) and 2(g), we can conclude a threshold value g t between g = 0.7 and g = 0.8 where the spinon upper bound touches the lower edge of the higher-energy band. This point signifies a hybrid of the doublon and quarton excitations becoming part of the conventional spinon continuum, which should be associated with a fractionalization mechanisms. Beyond the threshold value, the spectrum exhibits a continuum with a single peak, tending to the standard two-spinon continuum of the isotropic HAC when g → 1.
We can derive the threshold value for the merger of the different quasiparticle bands using the perturbative dispersion relations, which amounts to solving the equation πJ eff |sin(q/2)| = E 1 − E 0 , with the result g t = 9/(4π) ≈ 0.716. While we do not expect this value to be very precise, in Figs. 6(a)(b), S(q, ω) at q = π/6 and q = 5π/6 with g t = 9/(4π) are seen to contain a single broad peak. Figs. 6(c) and (d) present S red (q, ω) and S(q, ω) obtained from QMC-SAC calculations for g t = 9/(4π), where the spectrum begins to form a single band. The gap between the upper boundary of the two-spinon continuum and the dispersion relations corresponding to the intermediate-energy branch is closed at q = π. This space contains all the single-trimer eigenstates (see Fig. 3) in addition to the full space of states of the background spin chain, restricted to the total M = 0 sector. There are C N/2 N ground states with M = 0 formed by permutation without repetition of |0 1 · · · 0 1 0 2 · · · 0 2 . For states containing a doublon or a quarton, one excited trimer (represented by red color) is present instead of one of the doublet trimer ground states. The excitations correspond to changes in quantum numbers obeying |∆M| = 0, which in some cases require a flip of a background spin (represented by green color). Here, the excitations |0 2 → |1 1,2 and |0 2 → |2 1,2,3,4 are not included since a duplicate of the state can be found from among the already constructed configurations with |∆M| = 0.
Here, we should emphasize that there is no phase transition near this threshold, but still there is a dramatic change in the nature of the excitations of energy ∝ J 1 .

F. Quasi-particles in a truncated Hilbert Space
Motivated by the results in the preceding sections, we now more formally construct a truncated Hilbert space in which the number of internal trimer excitations is limited to one doublon or quarton (with no restriction on their internal magnetization). We could in principle carry out low-order perturbation theory in this space, but the spinons are not easily accounted for in this way. Here our goal is instead to demonstrate that the spinons, doublons, and quartons can all be present in the spectrum originating from a severe truncation of the Hilbert space of the trimer states. We therefore carry out full ED calculations and construct the dynamic structure factor based on only the approximation of truncated Hilbert space with at most one internal trimer excitation.
We construct a set of states Ω formed by the eigenvectors of isolated trimers as shown in Fig. 7. In order to capture the spinon continuum, we include all C N/2 N combinations of the trimer doublet ground states for which the total magnetic quantum number M = 0. For the single excited trimer in its excited state there are several options for replacing the doublet ground states, as indicated in Fig. 7. The Ω covers all the M = 0 sates including the spinons, one doublon and one quarton, which is a crucial condition for realizing the full spectrum in this truncated Hilbert space. The states contain and the total number of states φ i ∈ Ω, including those without internal trimer excitations, is The effective Hamiltonian matrix in Ω is where the states are visualized in Fig. 7.
To diagonalize the effective Hamiltonian, we first Fourier transform the spin operator S z q [also see Eq. (33) in Methods IV], with a unit cell index R = 0, · · · , N − 1, The spin operators in the truncated Hilbert space are given by where the intratrimer labels α ∈ {a, b, c}. Then, above spin operator in momentum space is written as where the momenta is still q = 2nπ/L, n = 0, 1, · · · , L − 1.
The dynamic structure factor is given bỹ where |Ψ m is the mth eigenstate with energy E m . Results for L = 24 are already enough to reveal the key features of the excitations and spectral functions. As shown in Fig. 8, when g ≤ 0.5 the spectral shapes and weights coincide well with the results of QMC-SAC (Fig. 2). For g > 0.5, the agreement gradually deteriorates, which confirms that the single propagating trimer excitations are no longer a good quasiparticles of the full system. In particular, the doublon and quarton bands fail to merge into the spinon continuum in the truncated Hilbert space, and this failure also naturally corresponds to the inability of the truncated Hilbert space to capture the fractionalization of the internal trimer excitations.
It is remarkable that the agreement with the full and truncated calculations is good up to g as large as 0.5. The truncated calculation also captures some of the arch feature of the upper bound of the spinon continuum for g > 0.5. This also implies that the fractionalization of the doublons and quartons has not yet set in when the arch forms. The fractionalization likely sets in only at g ≈ 0.72, when the arch touches the lower-energy portion of the continuum, and some of the highenergy excitations may remain confined until g = 1 (similar to the 2D Heisenberg model, where only some of the highenergy excitations exhibit signs of fractionalization [48]). At the same time the spectral weight of the low-energy spinons evolves to form the lower edge of the conventional spinon continuum as g → 1.

G. Simplified pictures of doublons and quartons
Here we develop an intuitive picture of the doublons and quartons, generalizing the standard cartoonish picture of spinons as domain walls in an antiferromagnet. In the conventional simplified description of spinons in the HAC, illustrated in Fig. 9(a), one starts from a staggered spin configuration (mimicking the quasi-ordered antiferromagnetic true ground state). Flipping one spin creates an excitation with |∆M| = 1, which in the actual spin-rotationally invariant system corresponds to an excitation carrying spin S = 1. The left and right misalignments of the flipped spin with respect to its neighbors can be regarded as two domain walls, and these domain walls can move by two lattice spacings at a time by flipping pairs of adjacent spins (conserving the magnetization). These mobile domain walls are the spinons, and once the domain walls have separated the originally flipped spin has lost its identity and is completely fractionalized into two independently propagated spinons.
Let us now extend this cartoon-like picture to the doublon excitation. As shown in Fig. 9(b), to create an excitation of this type with |∆M| = 1, we again start from an antiferromagnetic spin configuration with three parallel spins, but now the effective spin in the middle is of the excited type. We can again think of the misaligned spin configurations as associated with domain walls, and once these domain walls move away from the central doublet site they look just like standard spinons. However, they may not completely free, but bound to the still present central doublet spin. The consequence of the binding is that the central doublet propagates through the system dressed by spinons, and there should be a large number of internal modes of this composite excitations, leading to a band of finite width in energy of these excitations. In addition to the bound spinons forming the cloud, there should also be freely propagating spinons coexisting with the propagating central doublet excitations, since on top of such an excitation a pair of low-energy spinons can be created. However, the observed dynamic spin structure factor, where the initial excitation is created locally and later destroyed locally, should be dominated by states with only dressed central doublet and no free spinons, because of matrix element effects in the same way as two-spinon processes dominate the structure factor of the standard HAC.
Next we turn to the quarton, which offers several possibilities for cartoon states with |∆M| = 1 according to the trimer excitations listed in Fig. 3. For simplicity we consider those based on the S z = 3/2 states, where the effective particles can be fractionalized or not, as illustrated in Figs. 9(c) and 9(d). In Fig. 9(c), we first replace a spin S z i = −1/2 by the S z = 3/2 triplet state. Then we have ∆M = 2, and to bring this down to ∆M = 1 we flip one of the neighbor S z = 1/2 spins down. This creates a domain wall, which can travel away from the excited trimer in the way discussed above. On the other side of the excited trimer, a domain wall can also propagate out. Thus, as in the case of the central doublet, Fig. 9(b), the trimer quarton will create a cloud of bound spinons surrounding it, with the difference that the separation between the domain walls is an even number of trimers, instead of an odd number in the central doublet case.
In Fig. 9(d), we replace an S z i = 1/2 state by the S z = 3/2 triplon state, which gives us ∆M = 1 in a single step. Here we just illustrate how this trimer excitation can propagate with assistance of virtual spinons. Such processes are also possible with the central doublet in Fig. 9(b) and the first quarton in Fig. 9(c). These processes represent the motion of the center of mass of the spinon-dressed internal trimer excitations. If we restore spin-rotation symmetry, the cases depicted in Figs. 9(c) and 9(d) would no longer be separable, and the |S z | = 1/2 trimer excitations would also be involved on equal footing.
Clearly the above considerations only provide rough cartoons of the actual excited states, but in the same way as the simple pictures in Fig. 9(a) have been important in forming useful intuition about spinons, these similar pictures for the collective aspects of the internal trimer excitations are also useful as an illuminating complement to the spectral functions and perturbative dispersion relations. In particular, the fact that these high-energy excitations already contain spinons suggests that they eventually fractionalize into the conventional HAC spinons by unbinding when g → 1. At the same time, the internal trimer excitation becomes gradually more ill-defined, involving a larger number of trimers.

III. DISCUSSION
We have investigated the dynamic spin structure factor and the nature of the visible excitations of the spin-1/2 antiferromagnetic trimer chain by employing the QMC-SAC, ED, and approximate analytical methods. We showed that changes in the intertrimer interactions leads to different types of collective excitations related to the HAC and the internal trimer excitations, and we used a perturbative approach and ED calculation in a truncated Hilbert space to confirm the excitation mechanisms.
When g is small, three well separated excitation branches are present. The low-energy continuum corresponds to the deconfined spinon excitations of an effective HAC with coupling J eff = 4J 2 /9. In the full BZ of the chain, three such continua are present, with different spectral weight distributions. When g = 1, the trimer chain reduces to the conventional HAC with its two-spinon continuum of width ω ∝ J 1 . We have investigated the cross-over behavior from one type of spinon to the other kind, as well as the manifestations of the fractionalization of the higher-energy quasiparticles when they evolve into conventional spinon continuum as g → 1.
When g 0.2, the propagating internal trimer excitations form two weakly dispersive bands, which we named doublons and quartons according to the nature of the excited isolated trimers. For larger g, these two quasiparticle bands merge into each other, and when g → 1 they lose their identity completely as they fractionalize and evolve into the conventional spinon continuum. The coexistence of two kinds of emergent spinon branches and two bands of trimer excitations for intermediate g give rise to interesting spectral signatures. Perturbatively, we identify a threshold value g t = 9/(4π) ≈ 0.716 of the coupling ratio at which these excitations merge into a single continuum, and our numerical results confirm the behavior for g close to this value. Our calculations in a truncated Hilbert space involving all the spin states of the background chain but only one doublon and one quarton. Comparing the results with those of the other calculations confirm the stability of the internal quasi-particle excitations up to g ≈ 0.5, while for larger g the description of the full system requires a larger number of degrees of freedom to describe the fractionalization process.
It would be very interesting to explore the details of the fractionalization mechanism in future calculations using other theoretical and numerical approaches. The trimer chain is perhaps the simplest setting in which such mechanisms can be explored-theoretically as well as experimentally. We note that there are already examples of coupled-trimer quantum magnets as discussed in the introduction, for instance, A 3 Cu 3 (PO 4 ) 4 (A --Ca, Sr, Pb) [20][21][22][23]. From the inelastic neutron-scattering spectra measured at 8K in Pb 3 Cu 3 (PO 4 ) 4 [20], two flat excitations at ω ∼ 9meV and ω ∼ 13.5meV are observed, which are also revealed by the intermediate-energy (at ω ∼ J 1 ) and high-energy (at ω ∼ 1.5J 1 ) excitations in our theoretical results when g = 0.1. Since the intertrimer couplings are small in this material, the trimers are approximately isolated, some features like the energies and relative intensity of two flat bands (see Figs. 2 and 3 in Ref. [20]) can also be compared with our results. However, the trimers in Pb 3 Cu 3 (PO 4 ) 4 do not correspond directly to our linear chain model with only with nearest-neighbor couplings and g < 1.
It is very likely that materials can be synthesized that correspond closely to our model. Once such a quasi-1D material has been identified, our results will be helpful for interpreting inelastic neutron scattering and other experiments probing dynamical properties that beyond the spin waves and conventional spinons, High-energy (∼ J) spin excitations of quantum magnets are less studied than the low-energy modes, but are attract-ing growing interest motivated by the emergence of unusual features in the spectral functions of the high-T c cuprates [54][55][56][57][58], as well as in other antiferromagnets described by the the 2D Heisenberg model [48,59]. These features have been interpreted as partial fractionalization of magnons, which in the presence of additional interactions can evolve into full fractionalization in some parts of the BZ. Our results offer useful insights for further exploring coexisting exotic excitations and fractionalization within a relatively simple theoretical framework.

A. Stochastic analytic continuation of quantum Monte Carlo data
In this section, we outline the SAC of QMC data. The imaginary-time correlation function G q (τ) corresponding to S(q, ω) is given by from which S(q, ω) can in principle be reconstructed by inverting the relationship In reality, the inversion procedure has limited frequency resolution due to the incomplete information available from QMC calculations of G q (τ) on a finite grid of points and with statistical errors. Nevertheless, with the best available analytical continuation tools and small statistical errors achievable with long runs using efficient algorithms, quantitatively useful information can be extracted. In our QMC simulations [46], we obtain unbiased statistical estimates of G q (τ) for a set of imaginary-time points {τ i }. Since the statistical fluctuations for different time points are highly correlated, we also have to compute the covariance matrix. With a number of QMC data bins N B , based on sufficiently long simulation segments to be in practice uncorrelated, we obtain the averagesḠ q (τ i ) = b G b q (τ i )/N B and the covariance matrix .
In practice, we also normalize the correlation functions so that G q (τ = 0) = 1, which automatically removes the covariance corresponding to an overall uniform fluctuation of the correlations. The remaining covariance is still significant and has to be taken into account for the method to be statistically sound.
In the SAC process, S(q, ω) is parameterized with a large number of δ-functions. Instead of the commonly used fixed grid with equally spaced δ-functions with sampled amplitudes, we here use the approach where the δ-functions all have equal amplitude and instead the frequencies are sampled. The mean density of δ-functions, accumulated in a histogram, then represents the normalized spectral function S(q, ω). The normalization factorḠ q (τ = 0) is reintroduced after the analytic continuation so that the correct spectral weight is recovered.
The most complete discussion of the sampling process is currently in Ref. [48]. Here we just briefly review the variant of the method corresponding to Fig. 1(a) of Ref. [48], which was also used in Refs. [51,52] where some additional tests and comparisons with other methods are presented. For a given sampled set of the equal-amplitude δ-functions, the corresponding imaginary-time function G s q (τ i ) on the chosen set of points {τ i } is calculated according to Eq. (28). The goodness of fit, χ 2 , defines the closeness to the corresponding QMC data as: The spectral function is sampled in a Monte Carlo simulation using the probability distribution Here the fictitious temperature Θ is set so that which provides a natural scale of fluctuations at which we do not "fit to the errors" while a good fit is still automatically guaranteed.
We now turn to the Fourier transform of the spin operator, S z q , for the system with three spins per unit cell. Using the full BZ of the spin chain of length L = 3N, we define where q = 2nπ/L, n = 0, 1, . . . , L − 1, and the corresponding imaginary-time correlation function is It is also useful to consider a reduced BZ. Denoting by S z i,α the spin at position α ∈ {a, b, c} of the ith unit cell, we define where q = 2nπ/N, n = 0, 1, . . . , N − 1. The correlation function for the reduced BZ is then assembled as: In principle we could also consider other form factors internally in the unit cells, but for our purposes here it suffices to consider the reduced BZ with the uniform summation, along with results for the full BZ based on Eq. (34). We have performed the QMC calculation with the length of the spin chain up to L = 192 (N = 64). To obtain results representing the low-temperature limit T → 0, we scale the inverse temperature β = J 1 /T = 4L. This low temperature is also necessary in order to resolve the spectral features appearing at very low energies, which are reflected in the imaginary-time correlations at large τ. In the SAC procedure, the statistical noise of the underlying imaginary-time data is a decisive factor governing the frequency resolution. Normalizing G q (τ) by setting G q (0) = 1 as explained above, the statistical errors vanish as τ → 0 and approach roughly a constant value as τ is increased. This almost constant statistical error for large τ is a good measure of the level of the statistical errors [47,48]. We have performed sufficiently long calculations to achieve an error level of approximately 10 −6 for most of the results presented above.

V. DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request. In addition to the QMC-SAC calculations, we have also applied the ED method to obtain reliable results. ED is the most basic numerical method for determining the eigenvalues and eigenstates of a quantum magnet with up to tens of spin-1/2 spins. If all the eigenstates have been obtained for a very small system, Eq. (3) in main text can be computed directly, and somewhat larger systems can be considered within the Lanczos ground-state method combined with the continued fraction expansion [1]. To reach larger sizes, we can further use symmetries to block diagonalize the Hamiltonian and reduce the computation time and memory requirement. Although ED is limited to extracting exact information from relatively small systems, where finite-size effects are still significant, we can gain more reliable insights pertaining to the thermodynamic limit by comparing results from the two different approaches.

VI. ACKNOWLEDGMENTS
In Supplementary Figure 10, we present the spin excitation spectra of the trimer chain obtained by ED calculations for different values of the coupling ratio g. Here the ED results are calculated from a chain with 30 spins. We have applied broadening in the momentum direction to aid in the visualization of spectrum.
In the ED calculations we can clearly see fine structure due to the small number of δ-functions (which are broadened for visualization) in Eq. (3) of main text. The visible bands formed by the individual excitations of the small chain fall within the spinon continuum forming in the thermodynamic limit. The ED and QMC-SAC spectral weight distributions are similar, but the QMC-SAC exhibit the full spinon continua more clearly. The two distinct excitation branches above the spinon continuum are very similar in the two sets of results.

Supplementary Note 2: Effective Heisenberg coupling
To derive the effective HAC arising from the doublet trimer ground states, we formally split the original Hamiltonian into three-spin blocks and apply the generic Kadanoff method [2][3][4][5]. The trimer chain Hamiltonian is thus decomposed into the intratrimer (H t ) and intertrimer (H tt ) parts. The lowest eigenstates of each trimer are then used to construct the basis for the effective Hilbert space. To achieve an effective Hamiltonian with structure similarity to the original one, the full Hamiltonian is projected onto the effective Hilbert space. The intratrimer Hamiltonian, which only contains the J 1 couplings, reads where N represents the number of trimers and we again refer to Fig. 1 of main text for the definition of the intratrimer labels a, b, c. The intertrimer Hamiltonian is given by which only contains the interaction between the c site of jth trimer and the a site of ( j + 1)th trimer. Each trimer has two degenerate ground states (see Fig. 3 in main text); where |↑ and |↓ are the eigenstates of the spin operator σ z . The projection operator of the jth trimer, is built to project the Hamiltonian onto the low-energy subspace, where |⇑ and |⇓ are the renamed base kets in this effective Hilbert space. The effective Hamiltonian up to the first-order correction is where the total projection operator is P = N j=1 P j . The effective Pauli operators are given by where ξ β i is the coefficient required to preserve the SU(2) algebra. Inserting the effective Pauli operators into Eq. (42), we obtain the effective Hamiltonian which describes an isotropic HAC in which the effective coupling strength only depends on the intertrimer coupling, J eff = 4J 2 /9.

Supplementary Note 3: Dispersion relations
As we have discussed in the main text, the internal excitations of the trimer chain are almost localized when g is small, and cannot be classified as magnons or triplons. In order to obtain further insight into the nature of these excitations, we propose a scheme to calculate their dispersion relations. The results are shown in Fig. 4 of main text as a number of dispersion relations corresponding to propagation of internal trimer excitations in the background of a chain with defects mimicking the complex ground state of the effective HAC. The overall good agreement with the QMC-SAC results on the location of these excitations and their band widths suggest that our picture of the excitation is correct even though the calculation involves a very rough approximation of the ground state and a simple ansatz for its excitations.
We begin by assuming that the ground-state wave function of spin chain is a product state of ground states of each trimer, where |0 is one of the two ground states of a single trimer. Thus, we neglect the interactions that cause the effective HAC with its spinon continuum, and instead deal with 2 N degenerate ground states without enforcing any quantum numbers (spin or momentum). Note that we do not apply degenerate perturbation theory here, because we are targeting excitations above the 2 N -fold degenerate ground-state manifold. We also do not carry out a formal perturbation expansion, but construct intuitive variational states including the internal excitations of one trimer. If the rth trimer is excited from its ground state |0 to one of its higher states |1 , then the wave function with this localized excitation reads |ψ e r = |0 1 |0 2 · · · |1 r · · · |0 N .
We can give this excitation a momentum and calculate its dispersion relation. Since the trimer has two types of excitations (see Fig. 3 of main text), we will obtain two bands. Next, we first consider the lower ω = J 1 doublon excitation of the isolated trimer, and then defer the similar but more complicated case of the ω = 3J 1 /2 quarton excitation.

A. Intermediate-energy doublon excitation
For simplicity, we first consider the approximate uniform ground-state wave function with all the trimers in the same ground state of the isolated trimer, |0 2 = (|↓↑↑ − 2 |↑↓↑ + |↑↑↓ )/ √ 6. This state clearly has momentum q = 0. In the lowest excitation involving the internal trimer excitation, one of the trimers is put in its first excited state |1 1 = (|↓↓↑ − |↑↓↓ )/ √ 2 with ∆M = −1. Fourier transforming such a state, the one-trimer excited state in momentum space is given by The total Hamiltonian with periodic boundary conditions is the summation of intratrimer Hamiltonian H m and the intertrimer Hamiltonian H m,m+1 Next, we calculate the expectation values of this Hamiltonian in the ground state and the first excited trimer momentum state. For the ground state of trimer chain, the expectation value of H trivially evaluates to For the one-trimer excited state, we consider separately the two terms of the expectation value of H in momentum space: ψ q e N m=1 H m ψ q e and ψ q e N m=1 H m,m+1 ψ q e . The first term only contributes the gap E 1 − E 0 . To evaluate the second term, we let H m,m+1 act on the c spin of mth trimer and the a spin of the (m + 1)th trimer. The excited trimers labeled by l and r in the bra ψ q e and ket ψ q e should take all the values of [1, N]. Practically, only the cases l = r and l = r ± 1 contribute, however, because ψ q e and ψ q e are orthogonal for |l − r| ≥ 2. Therefore, the expectation value of H for the one-trimer excited state is H m,m+1 0 2 · · · 1 1 r · · · 0 2 FIG. 11. Graphical representation of the calculation of the dispersion relation ǫ(q) = H e − H g of the intermediate-energy excitation within our approximation. The doubly-degenerate trimer eigenstates are shown as lighter and darker blue (ground states) and orange or red (first excited states). With these states, the intermediate-energy excitations with |∆M| = 1 originate from a trimer ground state |0 1 r on the trimer located at r when excited to |1 2 r or |0 2 r excited to |1 1 r (and in the corresponding bra states we use the site index l instead of r). The surrounding trimers, which can be in either of their doubly-degenerate ground states, are not affected. The excitations are given momentum q, and the matrix elements contributing to the dispersion relation are indicated. In addition to the case r = l + 1 shown, r = l − 1 contributes equally and r = l gives a different contribution. No other relationships between r and l contribute.
Thus, the dispersion relation is of what we refer to as the intermediate-energy excitation in reduced BZ is independently of the length of the spin chain. Unfolding this dispersion relation by replacing q with 3q, we can obtain the dispersion relation in full BZ The first two terms describe the gap of the localized excited trimer, and the last two terms reveal the dynamical features dominated by the intertrimer interaction. By symmetry, for the excitation from ground state |0 1 to the other first-excited state |1 2 with ∆M = 1, the same dispersion relation, Eq. (51), is also obtained in the reduced BZ.
To better mimic the true many-body ground state, we can consider random arrangements of the states |0 1 i and |0 2 i on each site. Exciting from all possible state in this degenerate manifold will lead to a band of excitations that approximate the true bands arising from combinations of an internal trimer excitation and spinons. We always start from momentum 0 in the ground state, and in the excited state we consider the same configuration of |0 1 i and |0 2 i but with one of these trimer ground states at site r excited to the first trimer excitation |1 1 r or |1 2 r such that |∆M| = 1, and this state is given momentum q. Because of the form of the Hamiltonian, and examining the energy calculation as illustrated in Supplementary Figure 11, we can conclude that the dispersion relations are mainly dependent on the excited trimers and their neighbours. Translated to finite size, this means that the correct generic result is already obtained from a system with N = 4 trimers. Thus, we can generate 16 different dispersion relations for the intermediate-energy excitation, and because of symmetries there are only four different ones. In addition to the one already given in Eq. (51), the other three are the Eqs. (6)- (8) in main text. Although all possibilities have been considered above, the true ground state for g > 0 is a totalspin singlet satisfying M = i M i = 0, and the numbers of |0 1 and |0 2 used when forming ground states should therefore be equal. After considering this restriction, the result of Eq. (51) is eliminated since its approximation corresponds to a maximal-M uniform ground state. Only 6 dispersion relations described by Eqs. (6)- (8) in main text are left for N = 4 trimers. Each equation corresponds to two degenerate dispersion relations. All three dispersion relations are graphed in Fig. 5(a) of main text for g = 0.1. Two of the above dispersion relations are independent of q since the excitations are localized. Only the l = r terms in Supplementary Figure 11 contribute to these excitations, while the l = r ± 1 terms lead to propagation of the internal trimer excitations. According to the characters of this excitation starting from the g = 0 limit, we refer to it as the doublon.

B. Dispersion relations of the high-energy excitation
Calculating the dispersion of the high-energy excitations evolving from the trimer quartet at ω = 3J 1 /2 is more complicated but we follow the same strategy as for the doublon. In the latter case, the intermediate-energy band is formed by reconfiguring the singlet bond in the original ground state of the trimer in Fig. 3 of main text at a cost of J 1 . The higherenergy quarton band is formed out of trimers excited into their S = 3/2 highest state, which can be seen as the singlet of two spins excited into a triplet, at a cost of 3J 1 /2. This rough estimation of the excited energy matches very well the numerical ED and QMC-SAC results for small g values.
Again our calculation does not conserve the total spin of the collective many-body state, but we consider the excitations corresponding to ∆S = 1 on one trimer. In this calculation, we not only need to prepare the ground-state wave functions, but also select the excited states from the quadruplet. Here, N = 4 trimers are also adequate to do the calculation at our level of linear-in-g approximation. With each of the lth and rth excited trimers having four possible states, at least 2 4 ground states and 4 2 higher-lying excited states are needed, and there are 2 4 × 4 2 = 256 cases in total. After implementing the restriction of the ground-state wave function having M = 0, 96 cases are left. As a result, the dispersion relations of the high-energy excitation in the reduced BZ fall into 33 different cases and are displayed in Eqs. (9)- (17) in the main text.