{'fraction_non_alphanumeric': 0.04130119608553824, 'fraction_numerical': 0.016110909749909388, 'mean_word_length': 4.265865063460254, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum magnets with weakly confined spinons: Multiple length scales and quantum impurities

In magnets with strong quantum fluctuations, paramagnetic ground states with or without confinement of spinon excitations can be realized. Here we discuss the physics of the confined phase in the vicinity of a confinement--deconfinement transition: This regime, likely relevant to a multitude of frustrated spin systems, is characterized by multiple length scales: In addition to the magnetic correlation length, a confinement length can be defined, which can be probed, e.g., by local static measurements near non-magnetic impurities. We illustrate the ideas by explicit calculations for dimerized spin chains, but our qualitative results remain valid in higher dimensions as well. In particular, we study crossover from weak to strong confinement visible in the antiferromagnetic polarization cloud around a non-magnetic impurity. We also discuss the effective magnetic interaction between impurities, relevant for impurity-induced magnetic order, and consequences for nuclear magnetic resonance and neutron scattering experiments.


I. INTRODUCTION
The last decades have seen tremendous interest in exotic phases of correlated quantum matter characterized by elementary excitations with fractional quantum numbers. Some examples are Luttinger liquids, 1 spin liquids, 2 and fractional quantum Hall systems. 3 While it has long been known that fractionalization generically occurs in one space dimension (1D), the situation in higher dimensions D is quite different. It has become clear that fractionalized phases in D ≥ 2 can be connected to deconfined phases of certain lattice gauge theories. 4,5 Quantum antiferromagnets have played a central role in investigations of exotic phases. 5 Such systems can exhibit both antiferromagnetic and paramagnetic ground states. In the latter situation, which is generic to 1D due to the Mermin-Wagner theorem, quantum fluctuations prevent long-range order of the spins. In D ≥ 2, quantum paramagnetic states may obtained by the formation of singlet bonds of spins 1/2. For lattices with one spin 1/2 per unit cell, the singlet bonds can form either crystalline order accompanied by lattice symmetry breaking, leading to so-called valence bond solid (VBS) phases, or short-range ordered liquid states without any broken symmetry. 5 Remarkably, such spin liquids can be shown to be fractionalized, with elementary excitations being spin-1/2 spinons.
The search for microscopic realizations of spin liquids in D ≥ 2 remains difficult, because the relevant spin models display geometric frustration: this reduces the applicability of otherwise powerful Quantum Monte Carlo techniques due to the so-called sign problem. The first unambiguous demonstration of a stable phase with fractionalization in two-dimensional (2D) magnet was for the triangular-lattice quantum dimer model, which features both VBS and gapped spin liquid phases, the latter displaying topological order. 6 Models with exact gapped spin-liquid ground states have been constructed, like the toric code model 7 and the spin-1/2 honeycomb model, 8 both due to Kitaev. Heisenberg models with longer-range frustrated interactions are believed to be good candidates for topological spin liquids as well, but with the exception of the triangular-lattice ring-exchange model 9 no conclusive numerical results are available. (Gapless spin liquids with power-law correlation functions do not play a role in the subsequent discussion.) A special form of fractionalization has been argued to be possible at a continuous quantum phase transition between two conventional phases with different broken symmetries. 10 A candidate for this scenario of deconfined quantum criticality is the Néel-VBS transition of the square-lattice Heisenberg antiferromagnet supplemented by a ring-exchange interaction. 11 An interesting prediction of the deconfined critical theory is the existence of two distinct small energy scales on the VBS side of the quantum critical point: primarily, the transition is between a Néel state and a spin liquid, and the VBS is a secondary instability of the spin liquid. Consequently, dimerization and confinement set in at a much smaller energy scale compared to paramagnetic behavior.
A common feature of (i) a confined paramagnet near a confinement-deconfinement transition and (ii) a VBS phase near a deconfined critical point is that, in addition to the magnetic correlation length ξ mag , there exists a much larger length scale associated with confinement, which we dub confinement length ξ conf . In the following, we shall understand ξ conf as the typical spatial extension of a two-spinon bound state. Importantly, both ξ conf and ξ mag can be large compared to microscopic lengths, opening the possibility to find universal phenomena characterized by two length scales, in contrast to the conventional situation with ξ mag being the only large length. Indeed, a recent numerical study 12 of the frustrated regime of the J 1 -J 2 -J 3 Heisenberg antiferromagnet on a square lattice gave evidence for two distinct length scales, with ξ conf > ξ mag . The purpose of this paper is a detailed discussion of physical properties of quantum paramagnets near deconfinement. We believe this regime to be relevant for a rather broad class of frustrated magnets: Even if the elementary excitations are ultimately confined, there will be signatures of deconfinement on shorter length scales, which can be probed experimentally. The presence of two length scales, ξ conf and ξ mag , strongly modifies a number of phenomena, both in statics and dynamics. We present general arguments and support them by explicit calculations for 1D Heisenberg models. Although the mechanisms of fractionalization are quite different in 1D and higher dimensions, we argue that various consequences for observables are similar, and hence our results remain qualitatively valid in D ≥ 2 as well.

A. Observables
In the simplest picture, weak confinement can be understood in terms of a (weak) linear potential, V (x) = α|x|, describing the effective interaction between two spinons. The two-particle problem then has only boundstate solutions, i.e., all excitations are of two-spinon type. The lowest two-spinon triplet state takes the role of the conventional spin-1 (triplon) excitation, its size defining the confinement length ξ conf , with ξ conf → ∞ for α → 0.
A static measurement of ξ conf is possible, e.g., by spatial pinning of one of the two partners of the bound state. This can be achieved by doping with non-magnetic impurities: removing one spin-1/2 from the system breaks a singlet and liberates a spinon which now interacts with the vacancy (holon). Provided that the binding between spinon and holon is identical to that between two spinons, the size of the spinon-holon bound state gives a measure of ξ conf . As we show below, the size of the antiferromagnetic polarization cloud around the vacancy is given by ξ conf in the weakly confined regime (whereas it is determined by ξ mag in the conventional case of strong confinement). This polarization cloud, which can be detected e.g. by local magnetic probes, also determines the effective magnetic interaction between impurity moments and is hence relevant for impurity-induced magnetic order.
Weak confinement shows up as well in the dynamical properties, as it implies the existence of multiple bound states below the two-or three-triplon continuum. Hence, inelastic neutron scattering experiments, measuring the dynamic spin susceptibility, will detect a series of sharp dispersing modes below the threshold of the multiparticle continuum.
With increasing confinement potential, ξ conf decreases. The criterion for the crossover to strong confinement may be defined as ξ conf ∼ ξ mag . This is expected to coincide with the condition that the energy spacing between the lowest two-spinon bound states becomes comparable to the spin gap, and hence the higher bound states are no longer discernible in the dynamical susceptibility.

B. Outline
The remainder of the paper is organized as follows. In Secs. II-V we discuss the physics of weak confinement in one dimension: Sec. II summarizes the well-known properties of the frustrated dimerized spin chain, described by the J 1 -J 2 -δ model, and the methods applied in the following sections. Numerical results for the undoped spin chain are discussed in Sec. III while analytical and numerical calculations for the system doped with a single vacancy are presented in Sec. IV. Sec. V describes experimental implications of weak confinement for nuclear magnetic resonance (NMR) and neutron scattering measurements and for impurity-induced magnetic order. Finally, in Sec. VI we discuss the relevance of our findings for higher dimensions. A brief summary of our results closes the paper.

II. FRUSTRATED DIMERIZED SPIN CHAINS
In one space dimension, spinon excitations are quite generically deconfined. The standard spin-1/2 Heisenberg chain has a gapless excitation spectrum, while frustration can induce spontaneous dimerization leading to gapped (but still deconfined) spinons. To induce confinement, one has to resort to models with explicitly broken translational symmetry.

A. Model and phase diagram
We consider the so-called J 1 -J 2 -δ model, describing a frustrated dimerized antiferromagnetic (AF) spin chain: Here S i represents a spin-1/2 at site i, J 1 and J 2 are the antiferromagnetic first-and second-neighbor exchange couplings, and δ denotes the degree of dimerization, see Fig. 1(a). (We will employ J 1 = 1 in the following.) The ground-state phase diagram of the model (1) in the J 2 -δ plane is well known, 13,14 Fig. 2, and will be summarized in the following. We start with δ = 0. For small J 2 < J 2C , the ground state is unique and the excitation spectrum is gapless. At the critical value J 2C = 0.2411, 15,16 a continuous phase transition of Kosterlitz-Thouless type into a gapped phase takes place. Here, the system breaks an Ising symmetry accompanied by bond order, i.e., spontaneous dimerization. At J 2 = 0.5, the so-called Majumdar-Ghosh (MG) point, the two groundstate wavefunctions are given by a nearest-neighbor sin- The state |0 is illustrated in Fig. 1(b).
In the gapped phase, elementary excitations are topological defects (solitons or spinons) separating the states |0 and |0 . Shastry and Sutherland 18 performed variational calculations for the MG point using a reduced Hilbert space {|n m }, where the states |n m consist of two spins at sites 2n − 1 and 2m separated by nearestneighbor singlets, Fig. 1(c). It was shown that the lowlying excitations are pairs of propagating spinons, i.e., deconfined spin-1/2 (fractionalized) quasiparticles. The upper bound for the spinon dispersion, ǫ k = (1/2)(5/4 + cos 2k), was later confirmed by numerical calculations. 19 A finite dimerization δ opens a gap for all J 2 . The perfectly dimerized states |0 and |0 remain eigenstates of H ch along the so-called "disorder line" 2J 2 + δ = 1, but their degeneracy is lifted by δ > 0. The energy of the ground state |0 is E 0 = (−3/8)(1 + δ)(2N ). It was also shown that the static structure factor S(q) = (1/2N ) j exp(iqR j ) S 0 · S j peaks at momentum q = π and q < π, respectively, for points localized on the left (commensurate phase) and right (incommensurate phase) sides of the disorder line in the phase diagram. 13 The fact that the states |0 and |0 are no longer degenerate implies that the low-lying excitations are now bound states of two spinons, i.e., spin-1 confined objects. 19 Indeed, a finite value of δ introduces an effective attractive potential between two spinons: by converting one particular singlet of the state |0 into a triplet  (1). For δ = 0, the system is gapless for J2 < J2C = 0.2411 and spontaneously dimerized for J2 > J2C , whereas δ > 0 induces a gap for all J2. On the so-called disorder line, 2J2 + δ = 1 (dashed), the ground state is given by a singlet product state on the strong bonds; moreover, the disorder line separates regimes of commensurate (C) and incommensurate (INC) spin correlations. The shaded area is the weakly confined regime of main interest in this paper, with ξ conf ≫ ξmag. and then separating the two spinons as in the state |n m , the energy of the system increases by an amount proportional to |m − n|δ. Such an increase is due to the fact the singlets between the two spinons are now located at the weak bonds [compare Figs. 1.(a) and (c)]. 20 As noted above, the spatial extension of the lowest-energy spinonspinon bound state 21 provides a measure of the confinement length ξ conf , with ξ conf → ∞ for δ → 0.
In the following, we shall be interested in the spingapped regime of weak confinement, with ξ conf ≫ ξ mag . Thus we shall work at δ ≪ 1 and J 2 close to 1/2, Fig. 2. As the data analysis is much simpler for commensurate (instead of incommensurate) spin correlations, the regime to the right of the disorder line will be touched upon only briefly.

B. Methods
Ground-state properties of the spin chain (1) will be calculated numerically for finite-size systems with the aid of the density-matrix renormalization group (DMRG) method, 22 using the libraries provided by the ALPS collaboration. 23 The calculations for 2N = 100 . . . 400 sites employ open boundary conditions, such that strong bonds are located at the chain ends, Fig. 1 (a).
In addition, near the disorder line, variational wavefunctions consisting of nearest-neighbor singlets with additional spinons will be utilized, similar to those in Ref. 18. Details can be found in Appendix A.

III. SPIN CORRELATIONS OF UNDOPED SPIN CHAINS
As a reference, we calculate the bulk spin correlation function S Z i S Z j in the ground state of the Hamiltonian (1) using DMRG. We restrict our attention to small δ and the vicinity of the disorder line.
The typical behavior of S Z i S Z j within the bulk of the chain for J 2 = 0.45 and several values of δ is shown in Fig. 3 on a log-linear plot. One can clearly see an exponential decay of the spin correlations, i.e., S Z i S Z j ∼ exp (|i − j|/ξ mag ), as it is expected for a spin gapped system. Similar results were found for J 2 = 0.48.
In Fig. 4, we show the values of ξ mag , obtained from an exponential fit to the long-distance part to the data, as a function of δ for J 2 = 0.45 and 0.48. (In the following, all length scales are measured in units of the lattice spacing.) For a fixed value of J 2 , ξ mag increases as the dimerization decreases and then approaches a constant value in the limit δ → 0. ξ mag increases as J 2 is reduced for a fixed value of δ. This behavior is in agreement with the discussions of Sec. II. Recall that along the line δ = 0, ξ mag diverges as J 2 → J 2C from above. We note that on the disorder line ξ mag is minimal, but the correlations cannot be reasonably fitted to an exponential decay, as the ground state of the infinite system is given by the product state |0 (2).

IV. SPIN CHAINS WITH A NON-MAGNETIC IMPURITY
Vacancy doping in quantum paramagnets can induce effective magnetic moments. These impurity-induced moments reveal themselves in a Curie-like behavior of the uniform susceptibility, χ ∝ C/T , at intermediate temperatures. Remarkably, in the presence of three-dimensional couplings these induced moments can order at sufficiently low temperatures, thus changing the spin-gapped paramagnetic ground state of the pure compound into a magnetically long-range ordered state upon doping.
The appearance of effective moments upon doping vacancies is a property of systems with confined spinons; for spin 1/2 systems it is best visualized in terms of broken singlet bonds where one spin is replaced by a vacancy. The liberated spin 1/2 is confined to the vacancy at low energies, resulting in an effective spin-1/2 moment. 24,25 In contrast, in host systems with deconfined spinons, no moments are generated by introducing vacancies. This theoretical picture has been supported by various numerical studies, in particular on spin chain 26 and ladder systems. 27,28,29 In the following, we present analytical and DMRG results for spin chains with a single vacancy, with focus on the spin polarization cloud associated with the impurity-induced magnetic moment. It has been established that, in the strongly confined regime, the cloud consists of excited triplets, with the spatial size given by ξ mag . 27,28,29,30,31 In contrast, in the weakly confined regime the cloud size is given by the extension of the spinon-holon bound state, i.e., ξ conf . Indeed, we shall find a cloud of size max(ξ mag , ξ conf ), with a distinct crossover for ξ mag ∼ ξ conf . The shape of the spin polarization cloud allows us to distinguish the weakly from the strongly confined regime.
A remark on limits is in order here: the ground state of the infinite system with a single vacancy is doubly degenerate, corresponding to the residual vacancy-induced moment. Applying a uniform field h u ≪ ∆, where ∆ is the bulk spin gap, fully polarizes this moment and leads to a static staggered spin arrangement near the vacancy. A finite concentration of vacancies introduces another energy scale, namely the effective (maximum) interimpurity interaction J eff , see also Sec. V C. The following discussion of the polarization cloud around a vacancy assumes the limit h u → 0, keeping J eff ≪ h u ≪ ∆.

A. Variational calculations
Near the disorder line the low-energy states of the dimerized spin chain (1) are dominated by short-range singlet coverings: The S = 1 states |m n , Fig. 1(c), have a lower energy than any S = 2 state even for a large separation (m − n) between the two spinons. This is the basis of the variational method of Ref. 18, which we adopt here to calculate the spin correlations in the vicinity of a non-magnetic impurity.
Consider the frustrated dimerized spin chain (1) with strong bonds at the chain ends and introduce a vacancy at an odd site. For finite dimerization, δ > 0, the system has no inversion symmetry at the impurity site i imp even in the infinite-system limit. In this configuration, the liberated spin-1/2 stays mainly in the region i > i imp : recall that the low-energy sector is formed by states made out of short-range singlets. The free spin can only occupy a site i < i imp if a (long) singlet is formed between the spins at sites i imp − 1 and i imp + 1, which corresponds to a high-energy state. To simplify the variational calculations, we shall therefore ignore the region i < i imp , i.e., consider a spin chain with the vacancy placed at the chain end, as illustrated in Fig. 5(a). Also, we focus on parameter values on the disorder line, We restrict the variational Hilbert space to the states |p , which consist of an empty site at i = 1 and a spinon at site 2p in a background of nearest-neighbor singlets as illustrated in Fig. 5(b). The states |p are not orthogonal to each other. In fact, the overlap between two states |m and |p reads A long but straightforward calculation shows that the matrix elements of the Hamiltonian (1) with a vacancy at site i = 1 are given by with the step function θ(x) defined as in Eq. (A8). Instead of considering the non-orthogonal basis {|p }, we can define a orthonormal one 20 via the Gram-Schmidt orthogonalization procedure, namely Diagonalizing the Hamiltonian (5) in the basis {|φ p } gives the ground-state energy E and wavefunction |Ψ in the presence of an impurity: In Fig. 6, we show the energy E as a function of δ calculated for a spin chain with 2N = 600 sites using the procedure described above. The data in the region 10 −5 ≤ δ ≤ 10 −3 , which, as we will see in the next section, is deep in the weakly confined regime, can be fitted by a power law a 0 + a 1 δ 2/3 with a 0 = 0.6252 and a 1 = 1.6241. Such δ dependence of the ground state energy is chosen based on the analytical findings of Uhrig and collaborators, 20 who studied a similar problem via an effective one-dimensional Schrödinger equation. They showed that the solutions ψ(x) of the effective equation can be written as ψ(x) ∼ Ai(x/ξ + z i ), where z i are the zeros of the Airy function Ai(x) 32 and ξ = (3mδ/2) −1/3 is the characteristic length scale with m being the spinon mass, and the ground state energy E ∼ δ 2/3 . Once the coefficients C p are known, the spin polariza- can be calculated with the aid of Eq. (6) and the expression The results for δ = 10 −5 . . . 10 −3 are shown in Fig. 7 on a linear-linear and log-linear scales. One can see that the spin polarization has a staggered profile in the vicinity of i imp as is was already realized in previous calculations. 19,28,33 By decreasing the dimerization, not only the position of the maximum of S Z i moves away from the impurity site but also the distribution gets broader and broader. Such features indicate that the system approaches a deconfined regime as the dimerization decreases. Indeed, at δ = 0, the spin polarization peaks at the center of the chain.
We now discuss possible quantitative measures of the size of the polarization cloud. This is not trivial, as the long-distance decay deviates from a pure exponential, in particular for small δ. The asymptotic behavior of the Airy function Ai(x) suggests a fit using the expression for i > i imp , with parameters b 1...3 and ξ 1 , the latter being the length scale corresponding to the cloud size. This formula, being able to interpolate between a pure exponential and a faster decay, depending on whether b 3 = 0 or b 3 = 0, nicely fits the data as is illustrated in Fig. 7. Alternative definitions for the polarization cloud size (to the right of i imp ) are given by the position of the maximum ξ max of S Z i relative to i imp and by the center-of-mass coordinate of the cloud: 34 The variational results for ξ cm , ξ max , and ξ 1 as a function of δ are plotted in Fig. 8. As expected, the cloud size diverges as δ → 0, which corresponds to the MG point, and all measures are much larger than ξ mag ∼ 1. For instance, the data for ξ cm in the region 10 −5 ≤ δ ≤ 10 −3 can be fitted by a 0 + a 1 δ −1/3 with a 0 = −2.3952 and a 1 = 1.7031. Again, the δ dependence follows from the results derived in Ref. 20.

B. Numerical results
We now use DMRG to investigate spin chains with a single vacancy in a parameter regime away from the  Fig. 9(a). This choice is due to the fact that the spin polarization cloud around i imp is asymmetric (see the discussion below). In the following, we only discuss the properties of lowest state in the S Z tot = −1/2 sector. Some representative results for the spin polarization around the vacancy are plotted in Fig. 10 on a log-linear scale (only S Z i > 0 are shown). For the disorder line, the spin polarization curves are asymmetric with respect to i imp . On the left side of the impurity (i < i imp ) an exponential tail is evident while on the right side (i > i imp ) the shape of the spin polarization resembles the one described in the last section. In fact, the results of DMRG and the variational approach (solid lines in Fig. 10) are in quite good agreement. The behavior for J 2 = 0.48 line is richer: Notice that for the largest value of δ, S Z i is symmetric around the impurity site with exponential tails. As δ decreases, the curves become more and more asymmetric, i.e., they display features qualitatively similar to the ones observed for the disorder line. Indeed for δ < 10 −2 a pure exponential decay is no longer visible in the region i > i imp . Similar behavior (not shown here) is observed for J 2 = 0.45. We have also studied the vacancy-induced spin polarization for J 2 > 1/2: in this regime, the polarization displays incommensurate oscillations with an envelope similar to those in Fig. 10. However, the incommensurate behavior complicates numerical fits such that we refrain from a further analysis.
We proceed as in Sec. IV A and perform a quantitative analysis of the DMRG data for S Z i . Our results are summarized in Fig. 11 which shows the decay lengths ξ 1 , ξ 2 (described below), the center-of-mass position of the cloud, ξ cm (10), and the bulk correlation length ξ mag (Sec. III) as function of δ for three different values of J 2 . We also include the variational ξ cm determined in the previous section. One can see that the variational and DMRG results for ξ cm along the three lines are in quite good agreement. This implies that the restriction of the Hilbert space to the subspace {|p } is a reasonable approximation to describe the (weak) confinement physics. Recall that |p includes only singlets made out of neighboring sites and the ones between the i imp and the spinon site are located at the weak bonds in contrast to the ground state of the undoped chain whose singlet pairs are placed at the strong bonds.
The spatial decay of the spin polarization cloud, S Z i , can be analyzed with the aid of the expression (9). As S Z i is asymmetric due to the broken inversion symmetry with respect to i imp , we use Eq. (9) separately for the regions i > i imp and i < i imp to extract the length scales ξ 1 and ξ 2 , respectively. In both cases, we only include data points into the fit where S Z i is above the noise floor of the numerical data, e.g., i imp < i < 140 in the ξ 1 fit of the δ = 0.01 data in the top panel of Fig. 10. The ξ 1 and ξ 2 resulting from the fit are shown in Fig. 11 as well.
Let us discuss the results collected in Fig. 11. In all panels, one can identify three distinct regimes: (a) For large δ, we have ξ 1 ≈ ξ 2 ≈ ξ mag , i.e., both parts of the spin polarization cloud are determined by magnetic (bulk) correlations. We identify this as the strongly con- fined regime. (b) For small δ, ξ 1 ≈ ξ cm ≫ ξ 2 ≈ ξ mag . Hence, the right part of the polarization cloud is determined by the confinement length, which is moreover much larger than the bulk correlation length. These are characteristics of the weakly confined regime, and we conclude ξ 1 = ξ conf . (c) There is a pronounced crossover regime between (a) and (b), e.g. 10 −3 δ 10 −2 for J 2 = 0.48. As anticipated above, ξ conf diverges as δ → 0 while ξ mag approaches a constant value -this proves the existence of two length scales in the regime of weak confinement.
Qualitatively, the behavior of the spin polarization cloud is easily understood: for strong confinement, i.e. large δ, we have ξ conf ∼ 1 which implies that the spinon is close to i imp . Fig. 9(b) illustrates the most probable configuration of the ground state. Notice that the singlet pairs are at the strong bonds as in the ground state of the undoped chain. Therefore the enhancement of the spin correlations around i imp is mainly due to triplets which are excited by the coupling between the spin adjacent to the impurity (i.e. the spinon) and the environment. 27,28,29,30,31 The triplet physics is dominated by ξ mag , hence ξ 1 ≈ ξ 2 ≈ ξ mag . On the other hand, in a weakly confined regime of small δ, ξ conf ≫ 1. The spinon travels to distant sites i > i imp , Fig. 9(c), where it may also generate triplets. In contrast, configurations with the spinon at sites i < i imp are suppressed as they require a long, i.e., high-energy, singlet. Hence, the impurity-induced spin polarization for i > i imp profits from the travelling spinon and is determined by ξ conf (provided that ξ conf > ξ mag ) whereas the polarization for i < i imp decays on the bulk scale ξ mag .
To summarize, the size and shape of the spin polarization induced by a non-magnetic impurity allow to distinguish the regimes of weak and strong confinement. In the present 1D case, both length scales ξ conf and ξ mag can in principle be extracted from S Z i .

V. WEAK CONFINEMENT IN EXPERIMENTS
In this section, we briefly sketch how weak confinement manifests itself in standard experiments. As before, the discussion is focused on 1D systems, but qualitatively applies to D ≥ 2 as well, see Sec. VI.

A. Dynamic susceptibility: Multiple bound states
We start with the zero-temperature dynamic susceptibility χ( q, ω) of a clean (i.e. impurity-free) gapped confined paramagnet. The minimum of the triplon dispersion is located at some wave vector q = Q, where χ( q, ω = 0) displays a maximum.
In the conventional case of strong confinement, χ ′′ ( Q, ω) consists of a single sharp peak at the spin-gap energy ∆, corresponding to the triplon particle, and a continuum of multi-particle states, which typically starts at 3∆. (In low-symmetry situations, like the J 1 -J 2 -δ chain, a two-particle continuum starting at 2∆ at q = 2Q is present.) In contrast, in a situation of weak confinement, there will be additional sharp peaks between ∆ and the onset of the continuum. This can be easily understood: while ∆ is the energy of the lowest spin-1 two-spinon bound state, weak confinement implies that higher two-spinon bound states are energetically only slightly above ∆. The crossover to strong confinement may be defined by the criterion that ∆ becomes equal to the energetic difference between the two lowest S = 1 two-spinon states.
For the J 1 -J 2 -δ chain considered above, the dynamic structure factor has been discussed e.g. in Ref. 35. The energy spacing of the two-spinon bound states scales as δ −2/3 , i.e. the number of bound states below the continuum diverges as δ → 0. The precise structure of χ ′′ ( Q, ω) is determined by the energies and matrix elements of the two-spinon bound states, the latter giving the weights of the peaks in χ ′′ . The energies depend on the bound state number n as ∆ + n 2/3 . Ref. 35 also provides an estimate of the peak weights, but the employed continuum limit is questionable because the weight depends on the shortdistance behavior of the two-particle wavefunction which is strongly influenced by lattice effects due to the nonanalytic nature of the confinement potential V (x) = α|x|. Numerical results for χ ′′ ( q, ω) in the J 1 -J 2 -δ chain can be found in Refs. 19,36,37,38, which clearly show additional sharp modes below the continuum. In particular, Fig. 4 of Ref. 36 reports energies of triplon bound states on the disorder line for δ = 0.4. As both single-triplon states and two-triplon spin-1 bound states correspond to twospinon spin-1 bound states, their results can be directly cast in our language: there are two (four) two-spinon spin-1 bound states below the two-particle continuum at wave vector 0 (π/2), but already the second is very close to the lower edge of the continuum, consistent with δ = 0.4 being located in the crossover regime between weak and strong confinement, Fig. 11.

B. NMR line shape due to vacancies
We now return to confined paramagnets doped with non-magnetic impurities, which generate a spin polarization cloud as discussed in detail in Sec. IV. This spin polarization can be made static by applying an external field larger than the interimpurity interaction (see Sec. V C below). The method of choice would be spinresolved scanning tunneling microscopy (STM) imaging of the cloud which, however, is not readily available with atomic resolution. Alternatively, nuclear magnetic resonance provides a tool to detect the distribution of local fields. For instance, for doped antiferromagnetic spin chains, 39 a broadening of the NMR spectrum at very low temperatures is observed while, for a doped Haldane chain, the satellite peaks in the NMR spectrum even allow the experimental determination of the magnetic correlation length. 40 The Knight shift K i of the resonance field at the site i is proportional to the local electronic moment (omitting the factor of the electron-nucleus hyperfine coupling); the NMR spectrum simply corresponds to the distribution of the local shifts K i . We are interested in the staggered For the spin chain of Sec. IV B, we calculate the distribution of staggered Knight shifts using where f ǫ (x) is chosen as a Lorentzian with a small width ǫ in order to have a smooth I(x) curve.
In Fig. 12 we show examples of the NMR spectrum calculated using Eq. (12) and the DMRG results derived in Sec. IV B for a chain with a single impurity. The results are for J 2 = 0.48 and δ = 0.0001, 0.004, and 0.04 which are in the regions of weak, intermediate, and strong confinement, respectively. As it was already observed in Ref. 28 the main effect is the broadening of the spectral line.
Weak confinement leads to an asymmetric spectrum, with a pronounced shoulder for one sign of x. This asymmetry is directly related to the spatial asymmetry of the spin polarization cloud and can be understood as follows: for an impurity as in Fig. 9, sites with i > i imp (i < i imp ) have S Z i stag < 0 ( S Z i stag > 0), respectively. Hence, the part of the cloud dominated by weak confinement, i > i imp , contributes to I(x < 0), where the large extension of the cloud results in a prominent tail of I(x < 0). Obviously, for a finite concentration of impurities, evenly distributed on both sublattices, the asymmetry of the line shape will disappear, but the shoulder feature remains.  L (bottom panel). The symbols are the same as in Fig. 1.

C. Effective interaction between vacancy-induced moments
In the experimentally relevant situation of a finite concentration of non-magnetic impurities, the induced moments couple via an indirect exchange interaction mediated by the gapped host. If the underlying lattice is bipartite, then the arising magnetism is typically nonfrustrated, despite the strong disorder inherent to the problem. In D = 3 and at sufficiently low temperatures, this interaction may lead to impurity-induced magnetic order in the paramagnetic host. Experimentally, this phenomenon has been observed in a variety of magnets, e.g., Zn-doped CuGeO 3 , 41,42,43 Zn-doped SrCu 2 O 3 , 44 and Mg-doped TlCuCl 3 . 45 On the theoretical side, vacancy-induced magnetic ordering in such compounds is commonly described by effective models, which exclusively contain the impurityinduced moments of spin 1/2, which are the only lowenergy degrees of freedom below the bulk spin gap (for small impurity concentration). The coupling between these moments is written down as a Heisenberg-type model, with spins S i at the random impurity locations, and a pairwise exchange couplings J ij of the form Here r ij is the distance between two impurities at sites i and j. Note that the common assumption is that the decay length of the J ij is the bulk correlation length ξ mag . In the following, we argue that Eq. (13) does not apply to the regime of weak confinement: there, the decay length of the J ij is given by ξ conf instead. The simplest argument is that J ij is determined by the overlap between the spin polarization clouds associated with the vacancies at sites i and j. As shown above, the cloud size crosses over from ξ mag in the strongly confined regime to ξ conf in the weakly confined regime -the same now applies to the decay length of the J ij .
To make this explicit, let us consider the spin chain described by the Hamiltonian (1) with two vacancies, the first one located at an odd site and the second at an even site, as illustrated in Fig. 13. In this configuration, the two spinons stay mainly in the region of size L between the two impurities. As for the strongly confined regime, we assume that in the ground state the two spinons are in a singlet configuration with energy E 0 S=0 and define the effective exchange coupling as where E 0 S=1 is the energy of the lowest triplet configuration. 34 We proceed along the lines of Refs. 20,35 and study the two-impurity problem near the disorder line via an effective Schrödinger equation in the continuum limit, which captures the physics in the Hilbert space sector without excited triplets. The continuum limit is justified here, as the effective interaction turns out to be dominated by the long-distance tail of the wavefunction. Within this approach, the difference between the singlet and the triplet configurations is given by the boundary conditions imposed on the total wavefunction Ψ(x, x ′ ). For the triplet configuration, the two spinons are not allowed to penetrate through each other while such process is permitted for the singlet one. Moreover, if we assume that the total wavefunction can be written as a product, i.e., Ψ(x, x ′ ) ∼ ψ(x)ψ(x ′ ), then we can reduce the two-particle problem to a single-particle one. Such considerations lead us to describe a spinon of mass m at position x by the Schrödinger equation subject to the following boundary conditions singlet : ψ s (x = 0) = ψ s (x = L) = 0, triplet : ψ t (x = 0) = ψ t (x = L/2) = 0.
Notice that Eq. (16) reflects the fact that the motion of the spinon is restricted to the region between the two impurities while the second term of Eq. (17) simulates the nonpenetration condition for the triplet configuration. The differential equation (15) has an analytical solution given in terms of a linear combination of shifted Airy functions Ai(x) and Bi(x), i.e., 32 with the corresponding eigenvalue Based on the discussions on Sec. IV A, ξ = (3mδ/2) −1/3 is identified with ξ conf . x 0 is determined from the boundary conditions (16) and (17), more precisely, from the solutions of the following equation with l = L/ξ conf (singlet) and l = L/(2ξ conf ) (triplet). We solve Eq. (20) numerically, calculate the corresponding effective exchange coupling as a function of L/ξ conf according to Eq. (14), and plot the result in Fig. (14).
One can see that J eff decays exponentially with the renormalized impurity-impurity distance L/ξ conf . Indeed, the log-linear data in the inset of Fig. (14) can be fitted by ln(J eff ) = a 0 − a1(L/ξ conf ) with a 0 = 3.8681 and a 1 = 0.9491. As the calculation applies to a parameter regime of the spin chain where the bulk correlation length is of order unity, this shows that the effective interaction follows Eq. (13) with the replacement ξ mag → ξ conf .

D. Vacancy-induced order: Staggered magnetization
In a situation where the effective interaction discussed above leads to long-range magnetic order, the staggered magnetization, being the order parameter, will be strongly inhomogeneous. Nevertheless, the Bragg peak seen in the elastic neutron scattering signal at the ordering wave vector Q will be proportional to the total staggered magnetization of the system. For dilute impurities and neglecting fluctuation effects, this is just given by the number of impurities times the staggered magnetization per impurity, M stag . Note that this quantity can be much larger than 1/2, in contrast to the uniform moment associated with each impurity which equals 1/2.
In the following, we calculate the staggered magnetization for a spin chain with a single vacancy from the DMRG data of Sec. IV B. The behavior of M stag as a function of δ for the disorder, J 2 = 0.48 and 0.45 lines is displayed in Fig. 15. One can see that |M stag | increases as we move into the weakly confined regime but does not diverge in the limit δ → 0, despite the diverging cloud size, ξ conf → The finite value of M stag as ξ conf → ∞ is to be contrasted with the situation of an impurity near a magnetic quantum critical point in 2 ≤ D ≤ 4, where ξ mag → ∞ leads to M stag → ∞. 30 Again, the difference can be understood within the picture developed above: on the one hand, for ξ conf → ∞ and finite ξ mag , the polarization cloud is primarily made out of many magnetically neutral singlets, see Fig. 9. Indeed, in this case, it is possible to show analytically 20 that M stag approaches a constant as δ → 0. (This constant, however, depends on ξ mag , see Fig. 15.) On the other hand, for ξ mag → ∞, the cloud consists of abundant triplets, which lead to M stag → ∞.

VI. HIGHER DIMENSIONS
So far, we have discussed one-dimensional spin systems where deconfinement occurs generically in models with spin 1/2 per unit cell, and confinement has to be explicitly induced by dimerization. The situation in D ≥ 2 is quite different, as deconfinement is not generic. Instead, it has been shown that the low-energy physics of many spin models can be written in terms of gauge-field theories. Here, the constituents are fractionalized particles (spinons) being minimally coupled to a gauge field of Z 2 or U(1) symmetry. 2,4,5 The gauge theory may have both deconfining and confining phases. While the former display true fractionalization, the gauge-field dynamics in the latter leads to confinement of objects carrying a gauge charge (e.g. spinons). The deconfined phases display topological order, characterized by the suppression of topological defects of the gauge field and by non-trivial ground state degeneracies. In contrast to 1D systems, non-trivial quantum phase transitions between confined and deconfined phases can be driven by varying microscopic model parameters.
This discussion makes clear that the physical mechanisms for deconfinement in 1D and higher dimensions are quite different. However, in both cases, a system can be in a weakly confined regime, i.e., close to deconfinement, with qualitatively similar observable consequences. In the following, we briefly sketch which of our results can be directly carried over to D ≥ 2.

A. Confined and deconfined paramagnets
In spin models with spin 1/2 per unit cell, quantum paramagnets can come in different varieties: valencebond solid phases spontaneously break the lattice translation symmetry and display gapped spin-1 excitations, i.e., confined spinons. In contrast, spin liquids do not break any symmetry of the Hamiltonian and feature deconfined spin-1/2 spinons, which can be gapless or gapped. In addition, spin liquids display low-lying gaugefield excitations: in the Z 2 case, these are gapped vortices of the Z 2 gauge field, dubbed "visons", whereas the U(1) case displays a linearly dispersing "photon" mode. These excitations carry zero spin, i.e., translate into singlet excitations in terms of the original spin variables.
Weak confinement can for instance be expected (i) near a paramagnetic confinement-deconfinement transition, e.g., from a VBS to a Z 2 spin liquid on the triangular lattice 46,47,48,49 , and (ii) near a deconfined critical point, e.g., a transition from a VBS to a Néel state on the square lattice. In both cases, the energy scale associated to confinement will be small in the VBS phase close to the quantum critical point. Consequently, there will be a large length scale ξ conf , with ξ conf ≫ ξ mag , which characterizes the binding of spinons, which may be detected using nonmagnetic impurities as described above. In a VBS phase, ξ conf also plays the role of a characteristic domain wall thickness. 10 Candidates for a deconfined VBS-Néel transition are the J 1 -J 2 and J-Q square lattice spin-1/2 Heisenberg models, where the antiferromagnetic order is frustrated either by a second-neighbor (J 2 ) or ring-exchange (Q) interaction. Whereas the transition in both cases may be weakly first order, observable finite-temperature properties of the J-Q model have been found to be in agreement with predictions of deconfined criticality. 50 Indications for weak confinement have indeed been found in a numerical study 12 of the frustrated J 1 -J 2 -J 3 Heisenberg model on the square lattice, where the couplings J i connecting i-th neighbor sites are all antiferromagnetic. For J 3 = 0, this model displays a paramagnetic phase for 0.4 < J 2 /J 1 < 0.6, which has been proposed to be of columnar 51 or plaquette 52 VBS type. For J 2 = 0, the paramagnetic phase appears to be larger, 0.3 < J 3 /J 1 < 0.7, 53,54 and again columnar 53 and plaquette 55,56 VBS states have been suggested. Ref. 12 studied the J 1 -J 2 -J 3 model using exact diagonalization, near the line J 2 + J 3 = J 1 /2 where the VBS order is likely most stable. Interestingly, these calculations showed clear indications of the existence of two length scales: while the magnetic bulk correlation length, ξ mag , was rather short, a significantly larger scale could be identified in the decay of the spin polarization near a vacancy -consistent with our analysis, this is to be identified with the confinement length ξ conf . For instance, for J 1,2,3 = (1, 0.3, 0.2), ξ mag ≈ 0.7 and ξ conf ≈ 5.5. 12 These findings (i) suggest that weak confinement is common to many frustrated low-dimensional magnets, and (ii) support that the paramagnetic phase of the J 1 -J 2 -J 3 model is located near a deconfined critical point.
Among the various differences between 1D and higher dimensions is that weak confinement in D ≥ 2 implies low-lying singlet excitations, which are energetically located below the triplet gap; in the gauge-field language, these are, e.g., visons near a confinement-deconfinement transition. In contrast, in the J 1 -J 2 -δ chain there are no low-lying singlets, as any excitation requires the creation of two spinons with an energy cost of at least the triplet gap.

B. Implications of linear confinement
In the language of spinons coupled to gauge fields, a treatment of the weakly confined regime is difficult, as the gauge-field effects leading to confinement are nonperturbative. In the following, we simply assume that the gauge field induces weak linear pairwise confinement between spinons (as well as between spinons and holons) and discuss its consequences. In fact, linear confinement is a plausible result of VBS formation in D ≥ 2. 10 Note, however, that bound states can in principle exist in a situation of asymptotic deconfinement, see Sec. VI C.
The simplest possible analysis employs a linear confinement potential V (r) = αr in the continuum limit. Consider a confined paramagnet with a single vacancy, the latter introducing a spinon into the system. Assuming the spinon to be massive with quadratic dispersion and to be linearly bound to the vacancy, its motion is described by a D-dimensional version of Eq. (15): Here r is the radial distance from i imp . The solution of the above differential equation can be written as where Y m l (θ, φ) is a spherical harmonic. 32 Unfortunately, the resulting radial problem has no exact (analytical) solution. We then proceed, assuming that the wavefunction of the lowest energy vacancy-spinon bound state can be approximated by where ξ = (2mα) −1/3 and the coefficient β is variationally determined. This choice is again motivated by the asymptotic form of the Airy function Ai(x) [see Eq. (9)]. It also satisfies the boundary conditions ψ(0) = 0 and ψ(∞) = 0. By minimizing the total energy E = drΨ * (r)HΨ(r)/ dr|Ψ(r)| 2 , with H given by Eq. (22), it is easy to show that, for D = 2 and l = 0, we have where Γ(x) is the Gamma function. 32 The above expression assumes a minimum value E = 2.1967 J/(2mξ 2 ) for β = 1.1989. Similarly, for D = 3, one finds E = 2.5187 J/(2mξ 2 ) and β = 1.2682. It is worth mentioning that the same variational procedure applied to the one-dimensional case gives β = 1.1428 and E = 2.4429 J/(2mξ 2 ), which is a good approximation for the exact analytical solution of Eq. (22), E = 2.3381 J/(2mξ 2 ), calculated in Ref. 20. We conclude that the vacancy physics in the weakly confined regime in D ≥ 2 is also dominated by the confinement length: Eq. (24) has a maximum for r = ξ, which can be identify with ξ conf . As we move deep into the weakly confined region (smaller α), ξ increases. The spin polarization cloud around i imp , whose characteristic length scale is ξ = ξ conf , does not display a pure exponential decay, but instead has a profile similar to the 1D case. This also suggests that analyzing the numerical data in terms of an exponential decay, as done in Ref. 12, might be not appropriate in the weakly confined regime.
Two-spinon bound states can be addressed similarly. The starting point is the D-dimensional version of Eq. (10) from Ref. 35, which is nothing but Eq. (22) with the replacements: m → µ = m/2 and E → E − ∆. Here µ is the reduced mass, ∆ is the spin gap (the energy of the lowest bound state), and r is the relative radial distance between the two spinons. As a consequence, the characteristic length scale is now ξ = (2µα) −1/3 . The solutions of the effective Schrödinger equation, which describes the two-spinon bound states above the spin-gap ∆, also have the form (23). We assume that Eq. (24) is a good approximation for the radial part of the first boundstate wavefunction ψ 1 (r) and variationally determine the parameter β and the energy E 1 for the case l = 0. Once β is known, variational wavefunctions for higher bound states can be easily constructed just by multiplying ψ 1 (r) by a polynomial, e.g., ψ 2 (r) = [(r/ξ) 2 + c 21 ]ψ 1 (r), ψ 3 (r) = [(r/ξ) 4 + c 32 (r/ξ) 2 + c 31 ]ψ 1 (r), respectively, for the second and third lowest-energy bound state above ∆. The coefficients c ij are determined by requiring that the states are orthogonal to each other. Notice that Eq. (25) satisfy the boundary condition ψ(0) = 0 and they vanish for large r. The ansatz (25) follows from a careful analysis of the analytical solutions of the one-dimensional problem. Indeed, we applied the above procedure for the D = 1 case, compared the results with the analytical solution reported in Ref. 35 and verified that it provides good estimates for the energy of the lowest two-spinon bound states. For D = 2 and l = 0, we find again β = 1.1989 and that the energies of the first (E 1 ), second (E 2 ) and third (E 3 ) lowest bound states are given by E i = e i + ∆ with e 1 = 2.1967, e 2 = 4.0025, and e 3 = 5.4875 J/(2µξ 2 ). According to the discussion in Sec. V A, the strong-weak confinement crossover occurs when e 1 ∼ ∆. Similar considerations hold for D = 3.
The above analysis suggests many parallels between 1D and higher dimensions, but falls short of capturing shortdistance features of the bound states as well as effects of non-trivial statistics of the spinons.

C. Asymptotic confinement vs. short-range bound states
The preceding discussion of confinement has been simplistic in that it assumed that the presence or absence of long-distance (asymptotic) confinement also determines the nature of the low-energy excitations. While this is the case for the one-dimensional J 1 -J 2 -δ spin chain, which features a confinement potential of the type V (x) = α|x| for all length scales, this is obviously not true in general. While the long-distance behavior of the effective interaction, mediated by the gauge field, decides about asymptotic confinement or deconfinement, it is its shortdistance behavior which determines the existence and spatial size of possible bound states at low energies.
In other words, it is conceivable that a phase in D ≥ 2 spatial dimensions with asymptotic deconfinement of spinons and associated topological order displays short-distance spinon bound states. 57 In such a situation, the diagnostics proposed in this paper, i.e., an vacancyinduced magnetic moment with cloud size ξ conf and multiple sharp modes in χ ′′ below the continuum, will signal (weak) confinement, whereas the true nature of the phase is deconfined. Obviously, distinguishing such a deconfined phase with bound states from a confined phase is extremely difficult. (Numerically, topological groundstate degeneracies may be used to detect deconfinement.) To our knowledge, no explicit microscopic examples for fractionalized spin liquids with spinon bound states existing below the continuum are known to date, 58 which is also related to the difficulties in numerically studying the excitation spectra of the relevant microscopic models. Possible candidates are states with Z 2 fractionalization which are believed to generically display a non-zero spinon pairing amplitude. 59 It should be kept in mind that the bound states in question, which are determined by short-range physics, are not necessarily part of effective low-energy gauge-field descriptions.
It is worth noting that a spinless impurity in the VBS phase near a deconfined critical point has been shown 60 to bind a spinon on a length scale much shorter than the VBS (or confinement) length.

VII. SUMMARY
In this paper, we have argued that quantum paramagnets, in particular those with frustration, can display an interesting regime of weak confinement of spinons. Here, the long-distance properties are conventional, with elementary excitations being spin-1 triplons; however, the physics of spinons is visible at shorter length scales. Such weakly confined paramagnets are characterized by two potentially large length scales, namely the common magnetic correlation length ξ mag and a confinement length ξ conf , which measures the typical distance between weakly bound spinons. Weak confinement implies ξ conf ≫ ξ mag ; the crossover to strong confinement occurs at ξ conf ≈ ξ mag .
We have illustrated the physics of the weakly confined regime with explicit calculations for frustrated dimerized spin chains, where weak confinement can be realized for weak explicit dimerization on top of strong spontaneous dimerization. In particular, we have studied the static response of the system to nonmagnetic impurities: the staggered polarization cloud induced by the impurity decays on a length scale given by max(ξ conf , ξ mag ). We have also briefly discussed the structure of the dynamic susceptibility in the clean case and the effective interaction between impurity-induced moments. Qualitatively, our results remain valid in higher dimensions as well: in D ≥ 2, weak confinement is reached near a confinementdeconfinement transition, e.g., between a gapped spin liquid and a valence-bond solid, or near a deconfined quantum critical point between a Néel and a VBS phase.
However, a caveat in D ≥ 2 is that long-distance confinement (deconfinement) of spinons is not necessarily equivalent to the presence (absence) of short-distance spinon bound states: In principle, bound states forming conventional triplon excitations can exist even in deconfined spin liquids. This tantalizing scenario is left for future investigations.
Last not least, upon inclusion of mobile charge carriers, weak confinement can be expected to lead to an almost spin-charge separated metallic state. This may be a conventional Fermi liquid, but with small quasiparticle weights and strong incoherent continua in spectral functions: The carriers will get dressed by a large cloud of singlet excitations, in analogy to the dressing with triplet excitations near a magnetic quantum critical point. 61 The above equation can be written in terms of the operator h(i, j) (A1), namely (J 1 = 1) h(2i, 2i + 1) h(2i, 2i + 2).  Notice that |i, p and |i, p ′ are similar, but for the former we have i < p while for the latter, i > p. The results for the two remaining terms can also be expressed in terms of |p , |p ± 1 , |i, p and |i, p ′ . With the help of (A4) it is possible to show that 2(H −Ē 0 )|p = (4pδ + 1 + 3J 2 − 2δ)|p − J 2 |χ p whereĒ 0 = −(1/4)(1 + 2δ + J 2 )(2N ) and the state |χ p is given by Proceeding in the same way as for Eq. (4), one can show that the overlaps between |m and the states |χ p , |i, p , and |i, p ′ are given by where the function θ(x) is defined as Equation (5) follows from Eqs. (A6) -(A7).