Numerical study of the temperature dependence of the NMR relaxation rate across the superfluid-Bose glass transition in one dimension

TL;DR

This study uses large-scale simulations to analyze how the NMR relaxation rate varies with temperature and disorder in one-dimensional spin chains, revealing its ability to detect quantum phase transitions and characterize the TLL parameter.

Contribution

It introduces a comprehensive numerical analysis of the NMR relaxation rate across the superfluid-Bose glass transition, linking it to disorder effects and the TLL parameter in 1D systems.

Findings

$1/T_1$ detects the quantum phase transition.

$1/T_1$ probes the TLL parameter $K$.

Distribution analysis aids experimental interpretation.

Abstract

We study the nuclear magnetic resonance (NMR) spin-lattice relaxation rate $1/T_1$ in random one-dimensional spin chains as function of the temperature and disorder strength. In the zero temperature limit, the system displays a disorder-induced quantum phase transition between a critical Tomonaga-Luttinger liquid (TLL) phase and a localized Bose glass phase. The $1/T_1$ is investigated across this transition using large-scale simulations based on matrix product state techniques. We find that this quantity can detect the transition and probe the value of the dimensionless TLL parameter $K$. We also compute the NMR relaxation rate distributions for each temperature and disorder strength considered. In particular we discuss the applicability of the stretched exponential fit to the return-to-equilibrium function in order to extract the $1/T_1$ experimentally. The results presented here should provide valuable insights in regards of future NMR experiments in realistic disordered spin compounds.

Tables

Table 1. Table 1: Tomonaga-Luttinger liquid parameters u 𝑢 u and K 𝐾 K as well as the prefactor A 𝐴 A of the static transverse correlation function ⟨ S r ± ​ ( 0 ) ​ S 0 ∓ ​ ( 0 ) ⟩ delimited-⟨⟩ subscript superscript 𝑆 plus-or-minus 𝑟 0 subscript superscript 𝑆 minus-or-plus 0 0 \langle S^{\pm}_{r}(0)S^{\mp}_{0}(0)\rangle at zero temperature. For the first line, corresponding to the clean case, the parameters are known exactly  Luther and Peschel ( 1975 ); Lukyanov and Terras ( 2003 ); Hikihara and Furusaki ( 2004 ) . In presence of disorder, for h = 0.2 ℎ 0.2 h=0.2 and h = 0.3 ℎ 0.3 h=0.3 , the value of the parameters is determined numerically using density-matrix renormalization group (DMRG) and quantum Monte Carlo (QMC). See the appendix for technical details.

$h$	$A$ (DMRG)	$u/J$ (QMC)	$K$ (QMC, DMRG)
$0.0$ (exact)	$0.159572$	$0.429535$	$2.173408$
$0.2$	0.13084(4)	0.419(2)	2.094(8), 2.0899(6)
$0.3$	0.12448(7)	0.404(2)	1.958(9), 1.953(2)

Equations

$$\mathcal{H}=J\sum_{i}\Biggl{[}\frac{1}{2}\Bigl{(}S^{+}_{i}S^{-}_{i+1}+\mathrm{H.c.}\Bigr{)}+\Delta S^{z}_{i}S^{z}_{i+1}+h_{i}S^{z}_{i}\Biggr{]},$$

$$1-M(t)\propto {\mathrm{e}}^{-t/T_1},$$

$$\frac{1}{T_1}=\frac{\gamma }{2}\sum _{a,b\in [x,y,z]}{A}_{ab}^2{\int }_{-\infty }^{+\infty }\mathrm{d}t{\mathrm{e}}^{it{\omega }_0}{C}_{a,b}(t),$$

$$C_{a,b}(t)=⟨S_i^a(t)S_i^b(0)⟩-⟨S_i^a(t)⟩⟨S_i^b(0)⟩,$$

$$\frac{1}{T_{1}}=\gamma\int^{+\infty}_{0}\mathrm{d}t\,\Biggl{\{}A^{2}_{\pm\mp}\operatorname{Re}\Bigl{[}C_{\pm,\mp}(t)\Bigr{]}+A^{2}_{zz}\operatorname{Re}\Bigl{[}C_{z,z}(t)\Bigr{]}\Biggr{\}},$$

$$\frac{1}{T_1}=\frac{A\cos \left(\frac{\pi }{4K}\right)\mathrm{B}\left(\frac{1}{4K},1-\frac{1}{2K}\right)}{u}{\left(\frac{2\pi T}{u}\right)}^{\frac{1}{2K}-1},$$

$$1-\overline{M(t)}\propto \overline{{\mathrm{e}}^{-t/T_1}}={\int }_0^{+\infty }\mathrm{d}{T}_1^{-1}P\left(T_1^{-1}\right){\mathrm{e}}^{-t/T_1},$$

$$\overline{\mathrm{exp}(-t/T_{1})}\simeq\mathrm{exp}\Bigl{[}-(t/\tau_{\mathrm{str}})^{\theta}\Bigr{]},$$

$$K=\pi \sqrt{{\rho }_{\mathrm{s}}\chi }\mathrm{and}u=\sqrt{{\rho }_{\mathrm{s}}/\chi }$$

$$F(i,j,L)=2A\Biggl{[}\frac{f(i+j,2L)f(i-j,2L)}{2\sqrt{f(2i,2L)f(2j,2L)}}\Biggr{]}^{-\frac{1}{2K}}$$

Full text

I Introduction

Understanding phase transitions is one of the cornerstones of condensed matter physics. Among these, disorder-induced quantum phase transitions can lead to fascinating phenomena with exciting new phases of matter. A famous example is the Anderson localization in the absence of interaction Anderson (1958); Abrahams et al. (1979), where the electronic wave function is spatially confined due to impurities and resulting destructive quantum interferences. Here, the random environment — disorder — can completely block the transport and drive a metal-to-insulator phase transition. The presence of interactions and its interplay with disorder can qualitatively change the picture and lead to the many-body localization (MBL) phenomenon. Whereas MBL commonly refers to high energy properties of interacting disordered systems Gornyi et al. (2005); Basko et al. (2006); Altman and Vosk (2015); Nandkishore and Huse (2015); Abanin and Papić (2017); Alet and Laflorencie (2018), low-energy physics is also concerned with remaining open questions.

At zero temperature, this elusive quantum many-body localized phase of matter is known as Bose glass (BG), and first appeared in the context of 4He in porous media Fisher and Fisher (1988); Fisher et al. (1989). While no simple microscopic picture clearly emerges, the set of properties defining this quantum phase of matter, describing a bosonic fluid lacking superfluid coherence in a random environment, is pretty well-established: It has a finite compressibility, its low-energy spectrum is gapless, correlations are short-ranged, and there is no global phase coherence. Besides helium-4 in random media, the phase has been reported in various experimental systems such as amorphous indium oxide films with a transition from a superconducting to insulator phase Shahar and Ovadyahu (1992); Sacépé et al. (2008, 2011). It has also been observed in an array of quasi-one-dimensional samples of 39K cold atoms, subject to a quasiperiodic optical lattice D’Errico et al. (2014). Another type of systems in which the Bose glass phase has been investigated are antiferromagnetic Mott insulators Nohadani et al. (2005); Hida (2006); Roscilde and Haas (2006); Zheludev and Roscilde (2013) such as the spin-$1/2$ ladder compound $($CH${}_{3})_{2}$CHNH3Cu$($ClxBr${}_{1-x})_{3}$ Hong et al. (2010) and the Br doped spin-$1$ system Ni(Cl1-xBrx)2-4SC(NH2)2 Yu et al. (2010, 2012a, 2012b); Wulf et al. (2013); Orlova et al. (2017); Dupont et al. (2017a, b); Orlova et al. (2018) to cite but a few Yamada et al. (2011); Kamieniarz et al. (2016).

In one dimension, the microscopic mechanisms driving the superfluid-to-Bose-glass transition are still controversial Altman et al. (2004, 2008, 2010); Ristivojevic et al. (2012a); Pollet et al. (2013); Pielawa and Altman (2013); Pollet et al. (2014); Ristivojevic et al. (2014); Yao et al. (2016). Even in absence of disorder, in one dimensional systems, quantum fluctuations prevent the emergence of a global phase coherence that would result from the spontaneous breaking of a continuous symmetry. At best, one can expect a critical Tomonaga Luttinger Liquid (TLL) phase with a finite superfluid density and quasi-long-range order characterized by power-law decaying transverse correlations $\propto r^{-1/2K}$ at large distance $r$, where $K$ is the so-called TLL parameter Giamarchi (2004). In the presence of disorder, the TLL phase is expected to be unstable towards a Bose glass phase, unless $K>3/2$. In the latter case, a critical disorder strength is required to drive the system from a superfluid to a many-body-localized phase. Giamarchi and Schulz have shown in their seminal work that this transition belongs to the Berezinskii-Kosterlitz-Thouless (BKT) universality class Berezinskii (1971); Kosterlitz and Thouless (1973); Kosterlitz (1974) with $K\equiv K_{c}$ taking a finite value at criticality Giamarchi and Schulz (1987, 1988). The TLL parameter might take the universal value $K_{c}=3/2$ at weak disorder (compared to the bandwidth) Ristivojevic et al. (2012a, 2014) but the extension of this weak disorder regime remains unclear and other approaches Altman et al. (2004, 2008, 2010); Pollet et al. (2013); Pielawa and Altman (2013); Pollet et al. (2014); Yao et al. (2016) based on strong disorder analysis suggest a non-universal value. Only recently a numerical study attempted to clarify the different scenarios proposed and to precisely define the weak versus strong disorder regimes Doggen et al. (2017).

Spin compounds are one of the best candidates to address these open issues experimentally. For instance, the nuclear magnetic resonance (NMR) spin-lattice relaxation rate $1/T_{1}$ has proven to be a formidable probe for one-dimensional physics in Mott insulators Berthier et al. (2017) and might be very well-suited for the purpose. The $1/T_{1}$ in a TLL phase is expected to diverge algebraically $\propto T^{1/2K-1}$ at low temperature Sachdev et al. (1994); Chitra and Giamarchi (1997); Barzykin (2001); Klanjšek et al. (2008). This prediction has been perfectly checked against numerics in paradigmatic (clean) $S=1/2$ XXZ spin chains Coira et al. (2016); Dupont et al. (2016) for $T/J\lesssim 10$, where $J$ is the antiferromagnetic exchange coupling. Moreover, the parameter $K$ has been reliably extracted by fitting experimental $1/T_{1}$ measurements versus $T$ in various quasi-one-dimensional spin compounds, such as the spin$-1/2$ Heisenberg antiferromagnetic ladder $($C7H10N$)_{2}$CuBr4 (DIMPY) Jeong et al. (2013, 2016); Möller et al. (2017). The effect of inherent (albeit weak) three-dimensional couplings $J_{\mathrm{3D}}$ on the $1/T_{1}$ in spin materials has also been recently studied in Ref. Dupont et al., 2018 where the authors could properly define a temperature window $J_{\mathrm{3D}}/100\lesssim T\lesssim J/10$ for observing genuine one-dimensional physics. Therefore, this quantity might provide a final experimental answer to the disputed issue of $K_{c}$. Indeed, strong-leg ladders subject to an external magnetic field can realize a TLL with $K>3/2$ Giamarchi and Tsvelik (1999); Hikihara and Furusaki (2001); Giamarchi (2004) and numerous examples of observation of the BG phase have been given above, making us confident that, eventually, an ideal spin compound combining a strong enough attractive regime and controllable chemical disorder will be experimentally available in the future.

In this work, using large-scale simulations based on matrix product state, we study the $1/T_{1}$ across the superfluid—Bose glass transition in paradigmatic one-dimensional XXZ spin$-1/2$ chains with a negative Ising anisotropy and subject to a random magnetic field, as a function of the temperature and disorder strength. The results presented here should provide valuable insights in regards of future NMR experiments.

The rest of the paper is organized as follows. In Sec. II, we provide an overview of the studied model, the definition of the NMR spin-lattice relaxation rate $1/T_{1}$ and briefly present the numerical techniques. Results are then discussed in Sec. III. Finally, we summarize our conclusions in Sec. IV.

II Model, definitions and methods

II.1 Model

We consider the spin-$1/2$ XXZ chain in a random magnetic field described by the following Hamiltonian:

[TABLE]

with $i$ labeling the lattice sites, $J$ the overall energy coupling, set to unity in the following, and $\Delta\in(-1,1]$ is the Ising anisotropy. The random variables $h_{i}$ are drawn independently from a uniform distribution $\in[-h,h]$ where $h$ characterizes the disorder strength.

In the clean ($h=0$) case, it is well-known that in the low-energy limit, the model (1) can be described as a Tomonaga-Luttinger liquid Giamarchi (2004) with algebraically decaying correlations and a finite superfluid density (spin stiffness). This description only relies on two phenomenological parameters: $u$, the propagation velocity of the excitations in the system and $K$, the dimensionless TLL parameter governing for instance the decay of correlations. Here, $u$ and $K$ can be related to the microscopic parameter $\Delta$ from Bethe ansatz equations Luther and Peschel (1975). Through a Matsubara-Matsuda transformation Matsubara and Matsuda (1956), the spin Hamiltonian (1) can be mapped exactly to a hard-core bosons (HCB) model with $\Delta$ controlling the nearest-neighbor HCB interaction. Precisely, $\Delta<0$ ($K>1$) corresponds to an attractive interaction between the bosons and this part of the phase diagram is often referred to as the attractive regime. In presence of disorder ($h>0$), the critical TLL phase is expected to be unstable towards a localized Bose glass phase with exponentially decaying correlations and a vanishing superfluid density. However, for strong enough attractive interaction, $K>3/2$ ($\Delta<-0.5$), a finite disorder strength $h_{c}$ is necessary to drive the system from a TLL to a Bose glass phase, with the TLL parameter $K_{c}(\Delta,h_{c})$ taking a finite value at the transition Giamarchi and Schulz (1987, 1988). A schematic phase diagram of model (1) is shown in Fig. 1.

In this work, we focus on $\Delta=-0.75$ and vary the disorder strength $h$ from zero to $1/2$, covering the TLL and Bose glass phases with the transition happening at $h_{c}\simeq 0.37$ Doggen et al. (2017). Importantly, the temperature will be a parameter to be as close as possible to realistic experimental setups. Although the model (1) is unrealistic to describe spin materials because of the random fields with no reasonable origin, the essential ingredients are: (i) in the low-energy limit and in the clean case, a TLL description of the spin compound with a TLL parameter $K>3/2$, (ii) disorder whose microscopic origin is likely compound-dependent but which drives the system to a localized phase if sufficiently strong.

II.2 The NMR spin-lattice relaxation rate $1/T_{1}$

In experiments, the NMR spin-lattice relaxation rate $1/T_{1}$ is the inverse characteristic time that the targeted nuclear spins, polarized through a static external magnetic field, take to relax back to their thermodynamic equilibrium after a perturbation by an electromagnetic pulse of frequency $\omega_{0}$, changing their polarization Abragam and Carr (1961); Horvatić and Berthier (2002); Slichter (2013). The return-to-equilibrium process can be described as function of time $t$ as

[TABLE]

where $M(t)$ is the component of the nuclear spins along the static magnetic field. One can show that $1/T_{1}$ provides information on the electronic spins $\mathbf{S}_{i}=(S^{x}_{i},S^{y}_{i},S^{z}_{i})$, i.e., those described by the Hamiltonian (1), through the following expression,

[TABLE]

with $\gamma$ the gyromagnetic ratio, $A_{ab}$ the hyperfine coupling constant describing the interaction between nuclear and electronic spins and $C_{a,b}(t)$ the dynamical correlation function

[TABLE]

which is a local quantity because measured on the same site $i$ for both operators. Here $\langle\rangle$ indicates the thermal average and $S^{a}_{i}(t)=\mathrm{e}^{i\mathcal{H}t}S^{a}_{i}\mathrm{e}^{-i\mathcal{H}t}$ is the time-dependent spin operator in the Heisenberg picture. It is theoretically justified to take the limit $\omega_{0}=0$ since the NMR frequency is usually the smallest energy scale of the problem 111Taking this limit assumes some smoothness in the Fourier transform of $C_{a,b}(t)$ (time $t$ to frequency $\omega$), with no sharp contribution at $\omega\rightarrow 0$ that would not be captured by the NMR measurements because $\omega_{0}$ is actually finite. This is the case for the model studied here.. Moreover, due to the $\mathrm{U}(1)$ symmetry of the Hamiltonian (1), the $1/T_{1}$ can be reduced to the components $(a,b)=(\pm,\mp)$ and $(z,z)$, where the $x$ and $y$ spin components have been expressed using the raising and lowering operators 222In the $\omega_{0}\rightarrow 0$ limit, the contribution of $(+,-)$ becomes equivalent to $(-,+)$ and one can choose to work with the most convenient to compute their contribution to the NMR spin-lattice relaxation rate $1/T_{1}$.. The definition (3) now reads

[TABLE]

where the first and second terms of Eq. (5) can be labeled as transverse and longitudinal contributions respectively. In the following, we only focus on the transverse contribution which dominates over the longitudinal one from intermediate to low temperatures and set the experiment-dependent prefactor $\gamma A^{2}_{\pm\mp}$ to unity. For a TLL, the transverse contribution to the NMR relaxation rate has been found for $\omega_{0}/T\ll 1$ in the form Sachdev et al. (1994); Chitra and Giamarchi (1997); Barzykin (2001); Klanjšek et al. (2008)

[TABLE]

with $\mathrm{B}(x,y)$ the Euler beta function, $u$ and $K$ the TLL parameters and $A$ the prefactor of the static correlation function $\langle S^{\pm}_{r}(0)S^{\mp}_{0}(0)\rangle$ at zero temperature. The prediction (6) perfectly checks against numerics at low temperature, $T/J\lesssim 10$, without any adjustable parameter Coira et al. (2016); Dupont et al. (2016).

II.3 Numerical methods

II.3.1 Finite temperature

To compute the local dynamical correlation function (4) at finite temperature, we use the Matrix Product State (MPS) formalism Schollwöck (2011). Whereas a MPS represents a pure state, it can also be used for mixed state through the purification method Verstraete et al. (2004). The basic idea is to write the density matrix as a pure state in an enlarged Hilbert space with half physical and half auxiliary degrees of freedom (they can be taken as a copy of the physical ones). From a practical point of view, the corresponding infinite temperature pure state can be written down exactly as a MPS of bond dimension $m=1$: A product state of maximally entangled pairs of physical and auxiliary degrees of freedom. One can show that the pure state at inverse temperature $\beta=1/T$ is obtained by time-evolving the infinite temperature one with $\mathrm{exp}{(-\beta\mathcal{H}/2)}$, where $\mathcal{H}$ only acts on physical degrees of freedom. We perform the imaginary-time evolution using the time-evolving block decimation (TEBD) algorithm Vidal (2004) along with a fourth order Trotter decomposition Hatano and Suzuki (2005) and time-step $\delta_{\beta}=0.1$.

II.3.2 Real-time evolution

When the desired finite temperature state is obtained, a real-time evolution with $\mathrm{exp}{(-i\mathcal{H}t)}$ is carried out using the same TEBD algorithm as for the imaginary-time evolution (fourth order and time-step $\delta_{t}=0.1$). This is the most limiting part since the real time evolution of a quantum state produces a rapid growth of entanglement entropy Laflorencie (2016) while the efficiency of the MPS representation relies on low entangled states. In practice, this limits the maximum time $t$ one can reach in the simulation. Some workarounds have been developed in order to push the limits further like the linear prediction Barthel (2013) or evolving the auxiliary degrees of freedom with $-\mathcal{H}$ in real time Karrasch et al. (2013). We have used the latter, which can be seen as a local disentangling operation, of which the real time evolution by $-\mathcal{H}$ is just one possibility Hauschild et al. (2018). In this work, the maximum bond dimension of the MPS was set to $m=500$. Furthermore, the NMR relaxation rate being a local quantity, the dynamical correlation (4) is computed in the middle of the chain to avoid open boundary and finite size effects. Indeed, at long time and at finite temperature, the correlation decays exponentially over time with a characteristic time $\tau(T,h)$ which makes the choice of a sufficiently large system $L(T,h)$ enough to consider the results in the thermodynamic limit, i.e., $u(T,h)\tau(T,h)\ll L(T,h)$ with $u(T,h)$ some velocity accounting for the spreading of the excitation $S_{i}^{\pm}$ in the chain. This is especially true in regards of the definition (5) for the NMR relaxation rate as an integral over time of this exponentially decaying correlation function. This integration is performed numerically using the standard Simpson’s rule.

II.3.3 Disorder sampling

In the following, the local dynamical correlation (4) is computed over $N_{\mathrm{s}}\approx 500$ independent samples for each temperature and disorder strength considered. Therefore, it is the most demanding part numerically since the simulation of each independent sample is already quite demanding itself.

III Results

III.1 Time dependence of the local spin-spin correlation function

We first look at the real part of the transverse dynamical correlation function $C_{\pm,\mp}(t)$ defined in Eq. (4) and whose integral over time gives the NMR spin-lattice relaxation rate $1/T_{1}$ according to Eq. (5).

It is instructive to consider the disorder-free $h=0$ and $\Delta=0$ case, which in the low-energy limit, also belongs to the TLL phase. At this specific point of the phase diagram, an expression for $C_{\pm,\mp}(t)$ can be derived exactly and expressed as a Pfaffian Stolze et al. (1992, 1995). In some cases, it can be brought into a more explicit form, and one finds three distinct regions in time. After a very short time of order $\sim\mathcal{O}(1)$ (region A), the correlation shows a power-law decay (region B) before displaying an exponential decay $\propto\mathrm{exp}[-t/\tau(T)]$ at longer time (region C), where the decay time diverges algebraically with temperature $\tau(T)\propto 1/T$ for $T\ll J$. It is also well-known that the thermal correlation length $\xi(T)$ diverges at low temperature in the same way, $\propto u/T$ with $u$ the TLL velocity Pereira (2012), so that in the end, one can relate those characteristic thermal time and space quantities through $u\tau(T)\sim\xi(T)$. The intermediate power-law region grows larger and larger as the temperature is reduced to eventually take over the long time exponential decay at zero temperature, with no characteristic length scale in this limit, the system being critical. We note $\tau_{*}(T)$ the crossover time between regions B and C.

Away from the $\Delta=0$ point, no exact expression is available, but no major qualitative difference is expected at low temperature, with a finite thermal correlation length $\xi(T)\propto u/T$ and the existence of a characteristic thermal time scale $\tau(T)$: It is still expected to follow $\propto 1/T$ Pereira (2012), which we verified numerically at $h=0$ for $\beta J\gtrsim 5$ (data not shown). At high temperature however, and as pointed out in Ref. Stolze et al., 1995, it might not be as simple to relate the two quantities, displaying different temperature dependences: For instance at infinite temperature $\xi(\infty)$ is zero while $\tau(\infty)$ remains finite. In any case there exists, for the transverse dynamical correlation function in the TLL regime, characteristic thermal length and time scales in the temperature window considered in this work. More generally, they also depend on the disorder strength, i.e., $\xi(T,h)$, $\tau(T,h)$, $\tau_{*}(T,h)$, and correspond in that case to the disorder averaged quantity.

Plugging in disorder, the position of the crossover time scale $\tau_{*}(T,h)$ becomes more manifest with the decay time $\tau(T,h)$ being shorter as the disorder strength and temperature are increased. This results in a steepest slope of the long-time exponential decay, as visible in all three panels of Fig. 2. By increasing these parameters, one expects quantum coherence to be generally weaker over time, hence the shorter characteristic and crossover times $\tau(T,h)$ and $\tau_{*}(T,h)$. In the localized regime shown in Fig. 2 (c) for $h=0.5$, the width of the intermediate power-law region should be bounded, not only by thermal, but also by localization effects. Indeed, taken separately, both effects induce a finite length scale in the system: A thermal length $\xi_{\mathrm{th}}(T,h)$ or a localization length $\xi_{\mathrm{loc}}(T=0,h\neq 0)$. The dominant effect will be associated to the smallest of the two length scales. Deep in the localized regime, the disorder averaged localization length follows $\xi_{\mathrm{loc}}(T=0,h\neq 0)\propto h^{2/(2K-3)}$ Giamarchi and Schulz (1987, 1988) and it diverges exponentially when approaching the transition from the localized phase $h\rightarrow h_{c}$, expected behavior for a BKT transition Berezinskii (1971); Kosterlitz and Thouless (1973); Kosterlitz (1974).

Finally, the NMR spin-lattice relaxation rate is the integral over time of the real part of the dynamical correlation function displayed in Fig. 2. Because it decays exponentially with time, reaching a finite maximum time $t_{\mathrm{max}}$ is actually sufficient numerically to get an accurate estimate of the $1/T_{1}$.

III.2 NMR spin-lattice relaxation rate distributions

Computing the dynamical correlation (4) for independent disorder samples $(N_{\mathrm{s}}\approx 500)$ allows us to establish the corresponding distribution fo the NMR relaxation rate $P(1/T_{1})$ for each temperature and disorder strength, covering both TLL and BG phases as shown in Fig. 3. Overall, it is at high temperature that the distributions seem to be the narrower, independently of the disorder strength. With a very short thermal correlation length $\xi(T/J\gg 1,h)\ll 1$, one expects the effect of temperature to be dominant over disorder surrounding the site $i$ on which the local dynamical correlation function is being computed.

Interestingly, in the BG phase, see Fig. 3 (d), the $1/T_{1}$ distributions for $\beta J\lesssim 10$ seem to be double-peaked. This can be understood as a competition between thermal and localization effects. At fixed temperature and for a given disordered sample, the relevant characteristic length $\xi(h,T)$ surrounding the local site from which the $1/T_{1}$ is computed, is either going to be $\xi_{\mathrm{th}}$ or $\xi_{\mathrm{loc}}$, whichever is the shortest one. The distribution of $\xi(h,T)$ over many samples is going to be a weighted combination of the respective distributions of thermal and localization lengths, leading to a bimodal structure of the overall distribution. As discussed previously, a finite length scale induces a finite time scale, either through the decay time $\tau(T,h)$ or the crossover time $\tau_{*}(T,h)$, in the local dynamical correlation function, which will be reflected on the $1/T_{1}$ distributions by also displaying a double-peaks structure. In Fig. 3 (d) for $\beta J=16$, we do not observe two modes anymore because localization is dominant at that temperature. For $h=0.4$, see Fig. 3 (c), no double-peaks structure clearly emerges despite being in the localized BG phase. This is because we are very close to the transition, $h_{c}\simeq 0.37$ Doggen et al. (2017), with the localization length diverging exponentially, making it way larger than the thermal length.

Experimentally, one does not have access to the $1/T_{1}$ distributions but, for a disordered system to the average value of the return-to-equilibrium function $M(t)$ of Eq. (2),

[TABLE]

from which the average NMR spin-lattice relaxation rate $\overline{1/T_{1}}$ cannot be readily obtained. In order to access it, one usually makes the assumption that the sum of $\mathrm{exp}(-t/T_{1})$ originating from the disorder averaging in Eq. (7) can be approximated by a stretched exponential Johnston (2006),

[TABLE]

with $\theta$ and $\tau_{\mathrm{str}}$ two parameters fitted against the experimental data. One then considers that $\tau_{\mathrm{str}}\equiv T_{1}$, which is exact in the disorder-free case where $\theta=1$. Under the stretched exponential assumption, the distribution $P(1/T_{1})$ corresponds to the stretched exponential one Johnston (2006). It is clear that the distributions displayed in Fig. 3 do not correspond to pure stretched exponential distributions defined on the semi-infinite interval $[0,+\infty]$. For the disorder strengths considered, we do not observe rare events in the $1/T_{1}$ value leading to very broad distributions, over orders of magnitudes, as it was for instance observed in spin chains with random exchange couplings Shu et al. (2018). Here, the NMR spin-lattice relaxation rate seems to be bounded by its clean value at $h=0$. Moreover the stretched exponential distribution would not capture the double peaks at high temperature in the BG phase. Yet, surprisingly, the value of $1/\tau_{\mathrm{str}}$ extracted from the approximation (8) is very close to the average $\overline{1/T_{1}}$ value, as discussed in the following.

III.3 Temperature dependence of the NMR spin-lattice relaxation rate $1/T_{1}$

In the TLL regime, the NMR spin-lattice relaxation rate $1/T_{1}$ takes the form (6) in the low temperature limit, with no free parameter. This expression has been perfectly checked against numerics in clean XXZ chains for various $\Delta$ anisotropies and for $T/J\lesssim 10$ Coira et al. (2016); Dupont et al. (2016). For completeness, we display the disorder-free case at $\Delta=-0.75$ in Fig. 4 (a), where the crossover towards genuine the TLL physics with $1/T_{1}\propto T^{1/2K-1}$ is indeed recovered at low temperature. In presence of disorder, but below the critical disorder strength $h_{c}\simeq 0.37$ Doggen et al. (2017), we still expect a TLL phase and the algebraic temperature dependent expression in Eq. (6) should still apply. However for $h\neq 0$, there is no exact expression for the TLL parameters $K$, $u$ or the prefactor $A$ of the static transverse correlation function at zero temperature, and they must be computed numerically. These values are reported in Tab. 1 for $h=0$, $0.2$ and $0.3$. The technical details for their determination are explained in the appendix.

We first focus on the average $\overline{1/T_{1}}$ value, which for $\beta J\gtrsim 16$ agrees with the TLL prediction combining Eq. (6) and the parameters of Tab. 1. The crossover towards the TLL regime in presence of disorder happens at slightly smaller temperatures than the clean case. This can be explained by the behavior of $\tau_{*}(T,h\neq 0)$ and $\tau(T,h\neq 0)$: They both decrease with disorder strength, making the dynamical correlation drop much faster than in the clean case, see Fig. 2 (a). Nonetheless, the parameter-free TLL prediction makes it possible to precisely define the crossover temperature, even in the disordered TLL regime, and this should be a key element to be computed in a more realistically experimental model in order to interpret carefully experimental measurements.

As we cross the superfluid to BG transition at $h_{c}\simeq 0.37$ Doggen et al. (2017), we observe that on the log-log scale of Fig. 4 (a), the gradient of the average value of the NMR spin lattice relaxation rate starts to decrease at low temperature as the temperature is lowered, while it increases in the TLL phase. The gradient is plotted in Fig. 4 (c) where this effects is clearly visible and we believe is a signature of the transition at low temperature. In the TLL phase, the gradient should saturate to $1-1/2K$ as visible for the clean case with $1-1/2K(\Delta=-0.75,h=0)\simeq 0.77$. One would need to access lower temperatures in order to observe the saturation at $h=0.2$ and $0.3$.

In all cases, the average $\overline{1/T_{1}}$, as previously discussed, is not directly accessible in experiments where a stretched exponential fit to the return-to-equilibrium function $M(t)$, as in Eq. (8), is usually used and the parameter $1/\tau_{\mathrm{str}}$ interpreted as the “relevant” $1/T_{1}$. Numerically, for each disordered sample, we can compute $M(t)$ and fit its disorder averaged value $\overline{M(t)}$ against a stretched exponential, as it would be done in experiments. The parameters $1/\tau_{\mathrm{str}}$ and $\theta$ are reported in Fig. 4 (a) and (b) respectively. Although it is clear that the distributions $P(1/T_{1})$ shown in Fig. 3 are not those corresponding to a stretched exponential $\overline{M(t)}$, the value of $1/\tau_{\mathrm{str}}$ is in very good agreement with the exact average $\overline{1/T_{1}}$ value, in both the TLL and BG phases. Yet, we can note that a discrepancy seem to appear as we go deeper in the BG phase. For comparison, this is in sharp contrast with results in spin$-1/2$ Heisenberg chains with random exchange couplings, which realize in the low-energy limit the so-called random singlet phase Shu et al. (2018): At low temperature, it was for instance found that the average $\overline{1/T_{1}}$ diverges while the stretched exponential estimate was found to go to zero as the temperature was decreased. In that case, the authors found that the stretched exponential did not capture the average value but rather the typical value of the NMR spin-lattice relaxation rate, characterized by the median value of the distribution. We report the median value $\mathrm{med}(1/T_{1})$ in Fig. 4 (a) which behaves similarly to the average and the stretched exponential estimates. As thoroughly discussed in Ref. Johnston, 2006, the physical interpretations of the stretching exponent 􏰐$\theta$ are not straightforward, although it is commonly related to the width of the $1/T_{1}$ distribution. The value of $\theta$ displayed in Fig. 4 (c) decreases with the disorder strength and seem to decrease with the temperature in the BG phase while it remains roughly constant with temperature in the TLL phase.

IV Summary and conclusions

Using large-scale simulations based on matrix product state techniques, we have computed the NMR spin-lattice relaxation rate $1/T_{1}$ in random spin chains displaying a disorder-induced phase transition in the low temperature limit between a critical Tomonaga-Luttinger Liquid phase and a many-body localized phase, known as Bose glass. We have provided numerical evidences that this quantity versus temperature detect the transition and that it should be able to address the still controversial value of the TLL parameter $K$ at criticality, which might be universal Ristivojevic et al. (2012a, 2014). One indeed expects at low temperature that $1/T_{1}\propto T^{1/2K-1}$, where a clear identification of the crossover temperature below which this algebraic dependence becomes valid was possible.

We were also able to access quantities which are not experimentally, such as the $1/T_{1}$ distributions for different temperatures and disorder strengths in both phases. From this, we discussed the applicability of approximating the disorder averaged return-to-equilibrium function $\overline{M(t)}$ (the only quantity accessible in such NMR experiments) as a stretched exponential $\propto\mathrm{exp}[-(t/\tau_{\mathrm{str}})^{\theta}]$. Through this approximation, $1/\tau_{\mathrm{str}}$ is usually referred as the “relevant” NMR spin-lattice relaxation rate but there is no guarantee that it is equal (or simply related) to the actual disorder averaged value of $\overline{1/T_{1}}$. To emphasize this difference between $\overline{1/T_{1}}$ and $1/\tau_{\mathrm{str}}$, it was found that for spin$-1/2$ Heisenberg chains with random exchange couplings Shu et al. (2018), $\overline{1/T_{1}}$ diverges at low temperature while the stretched exponential estimate goes to zero. Here, we show that they both behave in a very similar way: A decisive point for experiments to reliably extract $K$ from $1/\tau_{\mathrm{str}}(T)$.

In the weak disorder regime (not considered in this work, corresponding to the TLL phase close to $\Delta\sim-1/2$), there might be a multiplicative logarithm correction to the $1/T_{1}$ expression (6) in the TLL phase since the static correlations have also been found to display such a correction in presence of disorder Ristivojevic et al. (2012b). This would be qualitatively similar to the one appearing in the better-known isotropic $\Delta=1$ (and clean, $h=0$) case Barzykin (2000, 2001); Dupont et al. (2016). However, its effect in presence of disorder would probably be way smaller since it arises with a log at the exponent $2/9$ in the static correlations, very difficult to capture numerically. For instance, there is no mention of including this correction to fit the static correlations in Ref. Doggen et al., 2017.

It would be interesting to consider next the dynamical spin structure factor $S(q,\omega)$, measured in inelastic neutron scattering experiments. Though it would be much more challenging numerically because involving dynamical correlations at all distances in order to compute the Fourier transform from real to momentum space. In comparison, the study carried out in Ref. Shu et al., 2018 on the dynamical properties of the random singlet phase, was performed in the canonical ensemble (imposing $S^{z}_{\mathrm{tot}}=0$), more amenable with the MPS techniques Nocera and Alvarez (2016). One possibility would be to use quantum Monte Carlo with analytic continuation, where the disorder average could be done before doing the continuation. For the Bose glass phase, the equivalent to the dynamical structure factor in the context of cold atom experiments and Bragg spectroscopy has been investigated at zero temperature in Ref. Roux et al., 2013, but the nature of the excitations in this phase remains an open question.

Acknowledgements.

We are grateful to S. Capponi, M. Horvatić, N. Laflorencie and J. E. Moore for interesting comments and discussions. This work was funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division under Contract No. DE-AC02-05-CH11231 through the Scientific Discovery through Advanced Computing (SciDAC) program (KC23DAC Topological and Correlated Matter via Tensor Networks and Quantum Monte Carlo). This research used the Lawrencium computational cluster resource provided by the IT Division at the Lawrence Berkeley National Laboratory (Supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231). This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. The code for calculations is based on the ITensor library 333ITensor library, http://itensor.org..

Appendix A Determining the Tomonaga-Luttinger Liquid parameters numerically

The TLL parameters $u$ and $K$, as well as the prefactor $A$ of the static correlation function $\langle S^{\pm}_{r}(0)S^{\mp}_{0}(0)\rangle$ at zero temperature have no exact expression at finite disorder strength and have to be determined numerically. The parameters $u$ and $K$ can be related through the hydrodynamic relations to the uniform susceptibility $\chi$ and the spin stiffness $\rho_{\mathrm{s}}$ Giamarchi (2004),

[TABLE]

which are both readily computed in quantum Monte Carlo (QMC) through the stochastic series expansion algorithm Syljuåsen and Sandvik (2002); Bauer et al. (2011). We compute these quantities for $N_{\mathrm{s}}\approx 10^{3}$ independent random systems of size $L=256$ with periodic boundary conditions at inverse temperature $\beta J=2^{11}$ using the $\beta$-doubling scheme Sandvik (2002). The temperature is sufficiently low to consider that the QMC algorithm is only probing the ground state. Moreover, these quantities do not show strong finite-size effects for the disorder strengths considered here, hence $L=256$ can be taken as the thermodynamic limit result. Similar simulation parameters have been used in Ref. Doggen et al., 2017, providing very satisfactory results. The average of $\rho_{\mathrm{s}}$ and $\chi$ over the $N_{\mathrm{s}}$ samples is then computed to extract the values of $K$ and $u$ according to Eq. (9). The values for $h=0.2$ and $h=0.3$ are reported in Tab. 1 in the main text.

Independently, simulations were also carried out on systems of size $L=256$ with open boundary conditions using the variational density-matrix renormalization group (DMRG) algorithm White (1992, 1993) at exactly zero temperature, and enforcing $S^{z}_{\mathrm{tot}}=0$. Here, the susceptibility and the spin stiffness are not as easily computed but one can compute the transverse static correlation function $\langle S^{\pm}_{i}(0)S^{\mp}_{j}(0)\rangle$ versus the distance $r=|i-j|$. It has been computed for $N_{\mathrm{s}}\approx 10^{3}$ independent samples and the average value has been fitted to the form,

[TABLE]

with $f(i,L)=L|\sin(\pi i/L)|/\pi$. This corresponds to the expected form of $|\langle S^{\pm}_{i}(0)S^{\mp}_{j}(0)\rangle|$ in a finite system with open boundary conditions, with $i$ and $j$ far enough from the boundaries, and with $|i-j|$ large enough Cazalilla (2004). In practice, we only consider $|i-L|>L/4$, $|j-L|<3L/4$ with $i<j$ and $|i-j|>10$. The fitted values of $K$ and $A$ are reported in Tab. 1 in the main text for $h=0.2$ and $h=0.3$. The values of $K$ extracted from the QMC and the DMRG agree well with each others.