Introduction

Randomness in interactions caused by imperfection is inevitable in solid-state materials. The presence of randomness is commonly believed to be “harmful”, i.e., it may weaken the symmetry-breaking order and/or destabilize the quantum critical point at low dimensions1,2. Nevertheless, it also gives rise to exotic quantum phases and phenomena that are unique in quenched disorder systems3,4,5,6,7,8,9. For example, the exotic Griffiths phase can emerge in the random systems, which is characterized by strong power-law singularity in statical and dynamical properties8,10,11.

The prototypical and minimal model for understanding the quantum Griffiths phase (QGP) is the random transverse-field Ising chain (RTIC), where spin couplings and transverse fields are not uniform but follow a distribution amongst different sites/bonds8,9. The quantum critical point in uniform quantum Ising chain is strongly modified by the random couplings and fields, into to a quantum Griffiths criticality with divergent dynamical exponent, and the conventional paramagnetic phase turned into a peculiar disordered one with off-critical singularity12,13,14,15. The exotic quantum Griffiths singularity, with no counterpart in clean systems, has been reported in experiments, including the ferromagnetic alloy Ni1−xVx16,17, thin gallium superconducting films18, and has been predicted in the randomly layered superfluids19,20, etc. In these examples, the QGP and related singular behaviors were analyzed phenomenologically, yet the complication of these systems hinders a concise microscopic theory for the presence of QGP. Therefore, finding a RTIC model material would then be of great significance and thus provide an unambiguous example and versatile platform studying the quantum Griffiths physics.

However, unlike the uniform Ising chain that has been materialized in crystalline compounds, including, e.g., CoNb2O621 and Sr(Ba)Co2V2O822, the RTIC model remains a toy model and its experimental realization has been absent for decades. A key challenge for the realization of the RTIC model is to introduce randomness in both transverse fields and couplings. It makes the RTIC model, albeit simple, barely realized in real materials. For example, in the diluted rare-earth Ising magnet, LiHoxY1−xF4, there are random couplings present23,24,25,26, while the transverse magnetic field is always external and virtually uniform. Similarly, there is also no quantum Griffiths physics reported in the Ising-like quasi-spin-chain antiferromagnet Sr(Ba)Co2V2O8 under an applied transverse field27,28,29,30. On the other hand, although the ground state, as well as thermodynamics of the RTIC model, are relatively well understood theoretically8,9,12,13,14,15, the dynamical properties of the quantum Griffiths phase remains largely unexplored7.

In this article, we remove the obstacle by uncovering that the rare-earth compound PrTiNbO6 experimentally realizes the RTIC model, through thermodynamic and dynamical measurements on the single-crystal samples. The intrinsic transverse fields in PrTiNbO6 stem from the energy splittings between the two lowest-lying crystal-electric-field (CEF) levels of Pr3+ ions. As the nonmagnetic Ti4+ and Nb5+ ions randomly take the same 8d Wyckoff position31, the CEF splittings follow a distribution and give rise to both random transverse fields and Ising couplings. We perform comprehensive dynamical measurements, including muon spin relaxation (µSR) and inelastic neutron scattering (INS). We find, besides the power-law thermodynamics, that the quantum Griffiths singularity is also evidenced by the pronounced spin fluctuations and continuous spin excitations at low temperatures. It turns out the RTIC model, with the parameters determined from fitting the thermodynamic data, can naturally explain the dynamical properties of PrTiNbO6. The many-body modeling also enables us to nail down the origin of the continuous spin spectra in the compound as the locally ordered quantum Griffiths rare region.

Results

Rare-earth random quantum Ising magnet

The structure of spin-chain material PrTiNbO6 is shown in Fig. 1a, where the Pr3+ ions contribute effective Ising magnetic moments owing to the two lowest-lying CEF singlets, |E1› and |E2›31. In the subspace of |E1› and |E2›, the matrix element of the magnetic moment is given by mj,j’α = gJµB‹Ej|Jα|Ej’› (α = x, y, z, and j, j’ = 1, 2), where gJ = 4/5 is the Landé g factor. At low temperatures (kBT << ‹E3−E1›, see below), the components of the effective Ising magnetic moments of Pr3+ are naturally obtained,

$${\mathbf{m}}^{\mathbf{x}} = {\mathbf{m}}^{\mathbf{y}} = \left( {\begin{array}{*{20}{c}} 0 & 0 \\ 0 & 0 \end{array}} \right),{\mathbf{m}}^{\mathbf{z}} = \frac{{\mu _{\mathrm{B}}}}{2}\left( {\begin{array}{*{20}{c}} 0 & G \\ {G^ \ast } & 0 \end{array}} \right),$$
(1)

where |G| = g ~ 5 is the effective pseudospin-1/2 g factor along the c (z) axis in PrTiNbO6. The eigenstates of Eq. (1) are the effective pseudospin-1/2 states,

$$\left| {S^z = \pm 1/2} \right\rangle = \frac{1}{{\sqrt 2 }}\left( {\frac{G}{{\left| G \right|}}\left| {E_1} \right\rangle \pm \left| {E_2} \right\rangle } \right),$$
(2)

with the Ising eigen-moments, mx = my = 0 and mz = ±µBg/2. The dipolar magnetic moments of Pr3+ is restricted along the c axis (Ising anisotropy), while the xy-plane moments are mostly multipolar. Therefore, the external magnetic field only couples to the z component of the pseudospin-1/2 moments, as seen in the magnetization measurements31. Furthermore, since the two lowest-lying CEF singlets of |E1› and |E2› form a non-Kramers quasi-doublet, the level splitting hx ≡ E2−E1 constitutes an effective transverse-field term in the pseudospin-1/2 Hamiltonian32,33,34, i.e.,

$${\cal{H}}_{SI} = \frac{{E_2 - E_1}}{2}\left( {\left| {E_2} \right\rangle \left\langle {E_2} \right| - \left| {E_1} \right\rangle \left\langle {E_1} \right|} \right) = - h_j^xS_j^x,$$
(3)

where hjx for site j follows a distribution as evidenced by the INS measurements (see below).

Fig. 1: Crystal structure, CEF splitting, and thermodynamic fittings of PrTiNbO6.
figure 1

a Crystal structure of PrTiNbO6 viewed in the ab plane, where the dashed blue lines indicate the unit cells. b Magnetic INS data (red), where the elastic and lattice INS contribution measured on the nonmagnetic reference compound, LaTiNbO6, has been subtracted. The gray squares present the total intensity of PrTiNbO6. The gray line shows the Gaussian fit to the elastic signal with the fitted resolution of σE = 0.48(1) meV, whereas the red line presents the fit to the magnetic data with the convolution function (see main text). c−e Combined fits to measured thermodynamic data31 using models no.1 (without randomness) and 2–4 (random Ising chain with similar parameters, see main text). f Power-law temperature dependence of the magnetic heat capacity below ∼2 K measured at zero field, as compared to the green line calculated by model no. 4. The lattice heat capacity measured on LaTiNbO6 has been subtracted in (e, f)31. Error bars on the experimental data points show a standard error.

Full size image

In Fig. 1b, the first CEF excitation is probed by INS at 5 K and is found much broader than the instrumental resolution of σE = 0.48(1) meV. Fitting the broad peak with the convolution function32,35, \(I_{\mathrm{m}}\left( E \right) = {\int} {\frac{{2I_0}}{{\pi \sigma _E}}\frac{{\omega _{\mathrm{L}}}}{{4(E\prime - \left\langle {E_2 - E_1} \right\rangle )^2 + \omega _{\mathrm{L}}^2}}\sqrt {\frac{{4{\mathrm{ln}}2}}{\pi }} {\mathrm{exp}}\left( { - \frac{{4{\mathrm{ln}}2(E - E\prime )^2}}{{\sigma _E^2}}} \right)dE\prime }\), we find a Lorentzian distribution of hjx which has a mean value ‹E2−E1› ~ 1.0 meV and a full width at half maxima (FWHM) ωL ~ 1.4 meV, where I0 is the integral intensity. This is in good consistency with our previous estimate based on the analysis of the heat capacity of Pr0.08La0.92TiNbO631. The next CEF excited energy is much larger, i.e., ‹E3−E1› ~ 11.9 meV ~ 140 K31. Henceforth, we focus on the lowest quasi-doublet below ~20 K, and safely ignore higher CEF levels. Besides, there exists antiferromagnetic (AF) Ising coupling ‹J› ~ 4 K between the pseudospin-1/2 moments via the “Pr(f)-O(p)-Pr(f)” superexchange path31,36,37. As the bond coupling Jj depends on the wavefunctions of the lowest-lying quasi-doublets at adjacent sites36,37, it is also randomly distributed in PrTiNbO6.

Determination of the RTIC model parameters

As the Ising moments in PrTiNbO6 are coupled through random couplings and under intrinsic random transverse fields, the microscopic RTIC Hamiltonian for the compound reads,

$${\cal{H}} = \mathop {\sum}\limits_{j = 1}^l {J_jS_j^zS_{j + 1}^z} - \mathop {\sum}\limits_{j = 1}^l {h_j^xS_j^x} - \mu _0H^z\mu _B\mathop {\sum}\limits_{j = 1}^l {g_jS_j^z} ,$$
(4)

where Jj represents Ising coupling strength between Pr3+ ions (Sjz moments), and Hz is the external longitudinal magnetic field. Other symmetry-allowed terms like (Sj+Sj+1++h.c.) are expected to be small (<1 K)31 and not included in constructing our minimal RTIC model for PrTiNbO6.

To accurately determine the parameters in pseudospin-1/2 Hamiltonian (Eq. (4)) of PrTiNbO6, we fit the thermodynamic data using four different models (Fig. 1c–e). In model no. 1, we neglect the randomness, and find an optimal parameter set J = 3.5 K, hx = 12 K, and g = 5.0 within the uniform spin chain model, which, however, fails to fit the thermodynamic data well (see Fig. 1), quantified by a rather large residual Rp = 35 (see Eq. (8) for the definition of Rp). It thus confirms the indispensable ingredient of randomness in PrTiNbO6.

Therefore, we introduce independent Lorentzian distribution of J, hx, and g (as the lineshape of the first CEF excitation indicates in Fig. 1b) in the models no. 2−3, and find their optimal values by fitting the experimental thermodynamic data. In model no. 2, tensor renormalization group (TRG) calculations are performed on rather long spin chains with length l = 1000, averaged over 20 sample chains, which leads to a very accurate estimate of parameters as ‹J› = 3.5 K, ∆J = 3.1 K (∆ represents the FWHM of distribution, same below), ‹hx› = 12 K, ∆hx = 11 K, ‹g› = 5.0, and ∆g = 0.05. The fitting with model no.2 is excellent above 1.5 K, with significantly reduced least-Rp ∼ 10. In no. 3 model, we perform exact diagonalization (ED) calculations of l = 12 chains averaged over 48 samples, and get ‹J› = 3.8 K, ∆J = 2.8 K, ‹hx› = 11 K, ∆hx = 10 K, ‹g› = 5.0, and ∆g = 0.06 with the least-Rp = 8.6. Moreover, to further confirm the microscopic RTIC model for the compound, in model no. 4 we take a fixed distribution of {hjx} measured by LET at 9 K with a high resolution (Supplementary Note 2) that is slightly asymmetric. We then can fit the thermodynamic data down to 0.4 K, by fine tuning other parameters, and obtain ‹J› = 3.5 K, ∆J = 11 K, ‹g› = 4.9, and ∆g = 0.05, with Rp = 14.

Overall, from the fittings with random Ising models no. 2–4, we obtain a reliable and consistent estimate of random couplings and fields in PrTiNbO6. In particular, we find ‹hx›/‹J› (or ∆hx/∆J) ∼ 3 is much larger than the critical point of 1/2 (We use the spin operators, S = σ/2, in the RTIC Hamiltonian (Eq. (4)), and thus there exists a pre-factor of 1/2 in front of the critical value reported in, e.g., ref. 12, where Pauli matrices σ were taken.), located deeply in the disordered QGP regime. Nevertheless, due to interaction randomness, there exists quantum Griffiths singularity in the low-temperature scalings of thermodynamics. The susceptibility rapidly increases as T lowers, and the magnetic specific heat exhibits an algebraic behavior below 2 K. The latter indicates a nonzero dynamical exponent z ∼ 0.45 characterizing the quantum Griffiths singularity, which otherwise should be zero if the transverse fields were uniform.

Quantum spin fluctuations

Zero-field (ZF) µSR is an extremely sensitive probe detecting the fluctuations of internal local fields in magnetic materials. Figure 2a, b shows general Gaussian dynamic Kubo-Toyabe (DKT) relaxations38 (see the Supplementary Note 3), due to the relatively large distribution of the local magnetic fields induced by the strong hyperfine couplings of 141Pr3+ 39,40, 㵋|Bloc|› ∼ ν, where γµ = 135.5 MHz/T is the µ+ gyromagnetic ratio, Bloc is the local magnetic field around the implanted µ+, and ν is the fluctuation rate. As shown in Fig. 2c, the local fields, perpendicular (‹|B⊥|›) and parallel (‹|B|||›) to the incident µ+ beam, are found to have distribution widths in the order of ∼1 mT, similar to other quantum disordered Pr-based magnets39,40. We further fit to the temperature dependence of the average absolute value of the local field using ‹|Bloc|› ∝ tanh[‹hx›/(2kBT)] by ignoring the distribution of hx and spin–spin interactions39, and find ‹hx› ≡ ‹E2−E1› ∼ 1.2−1.8 meV (see Fig. 2c). The average value of the intrinsic random transverse fields obtained from µSR measurements is consistent with the above INS results (Fig. 1b).

Fig. 2: µSR measured on the single crystal of PrTiNbO6.
figure 2

a, b Selected ZF-µSR spectra with the incident beam perpendicular to the c axis in the dilution refrigerator and 4He cryostat, respectively. The colored lines represent the Gaussian DKT fits. c Temperature dependence of the average absolute value of the local-field components perpendicular (‹|B⊥|›) and parallel (‹|B|||›) to the incident beam, with the colored lines showing the fits (see main text). d Temperature dependence of the fluctuation rate of the local field around the implanted µ+. Error bars on the data points show a standard error.

Full size image

The zero-temperature limit of the fluctuation rate of the local field reflects the quantum fluctuations of spins, which are pronounced in quantum spin liquid candidates. For instance, a very strong quantum fluctuation rate of ∼9 MHz has been reported in YbMgGaO4 at ∼0.07 K41,42. As shown in Fig. 2d, the zero-temperature fluctuation rate in PrTiNbO6 is ν0 ∼ 0.8 MHz, indicating the existence of significant quantum fluctuations in the compound. Both the magnetically disordered bulk (medium) and the putative fluctuating rare regions under the tunneling fields (see below) can contribute to the dynamic µSR signal. Notably, broad spin fluctuations have also been reported in the random system Ni1−xVx with quantum Griffiths singularity observed17.

Spin dynamics and excitation continuum

Neutron scattering measurements provide valuable information on the dynamical spin correlations and excitations in quantum magnets43. For PrTiNbO6, despite the AF Ising coupling of ‹J› ∼ 4 K, no magnetic freezing is evidenced, as there is neither splitting of zero-field-cooling and field-cooling susceptibilities down to 1.8 K, nor sharp “λ” heat capacity peak down to ∼0.1 K. Moreover, no frequency dependence of the ac susceptibilities was observed31. This is confirmed by the µSR and neutron scattering measurements in this work (see Figs. 2 and 3). As indicated by the columnar-like shape in Fig. 3a, we find the correlations are mainly along the b axis (i.e., spin-chain direction) but not along the orthogonal a axis, suggesting a negligible interchain correlation32,33, due to relatively large interchain distance (≥ 5.4 Å)31.

Fig. 3: INS spectra of the PrTiNbO6 single crystal measured with Ei = 7.5 meV.
figure 3

a Wave-vector dependence of the spin excitation continuum measured at 0.4 ≤ E ≤ 0.8 meV at 0.1 K. b Energy dependence of the continuum along the b axis. c, d Wave-vector dependence of the dynamic spin-correlation function measured at selected longitudinal fields at 0.1 K and measured at various selected temperatures under zero field, respectively. The colored lines are guide for the eye and the dashed gray lines mark the high-symmetry reciprocal-space points (Γ and M points). In c the FWHM of the hump (ωQ) is marked, and the inset shows the data integrated over different energy ranges. e Energy dependence of the measured intensities at the M (solid scatters) and Γ (open scatters) points, respectively. The 8.8 T data provides the elastic background below 2 meV, which has been subtracted in the INS data shown in (a−e). Error bars on INS data indicate one standard error propagated from neutron counts.

Full size image

In Fig. 3, we can observe a continuum of excitations broadly distributed in both momentum and energy space even at the lowest temperature of 0.1 K, where the magnetic entropy is estimated to be almost zero31. As shown in Fig. 3e, the applied longitudinal magnetic field gradually shifts the spectral weight toward higher energies, thus indicating the magnetic origin of these excitations. Owing to the dominant single-ion physics, ‹E2−E1› ≡ ‹hx› >> ‹J›, the wave-vector dependence of the spectra is negligible at E >> ‹J›, which can be evidenced in Fig. 3e, while below ∼1 meV ∼3‹J› the wave-vector dependence of the dynamic spin correlation Sm is significant. Therefore, in the analysis below, we focus on the energy range between 0.4 and 0.8 meV, and the energy integration of the entire continuum slightly weakens the wave-vector dependence in the integrated spin structure Sm (see the inset of Fig. 3c).

There exist two magnetic Pr3+ ions per unit cell along the chain (see Fig. 1a), and the broad humps appear around the M points, i.e., K = −3, −1, etc. (see Fig. 3c, d), suggesting the formation of AF correlations in the system. By applying longitudinal field (e.g., ∼2 T) or increasing temperature (e.g., ∼4.5 K), the wave-vector dependence of Sm as well as the overall intensity are clearly suppressed in PrTiNbO6 (see Fig. 3c, d). These observations indicate that the AF correlations, driven by the coupling J of ∼4 K, are mostly restricted within the nearest neighbors, i.e., with a short correlation length ξb ∼ b/(πωQ) ∼ 2.4 Å (see Fig. 3c for ωQ)44,45,46. Even for a low-energy transfer down to ∼0.05 meV, the peak of the INS spectrum (around the M point) remains broad in PrTiNbO6 (Supplementary Note 1). This is even quite different from other disordered magnets, e.g., YbMgGaO442,46, where the correlation length increases at lower transfer energies. We ascribe the difference to the intrinsic random effects of the non-Kramers magnet PrTiNbO6 (see below).

RTIC dynamics and quantum Griffiths rare region

With the RTIC models (no. 1−4) determined from fitting thermodynamic properties, in Fig. 4 we compute the dynamical spin structure factor and spectral weight of the RTIC model, and compare directly to experimental measurements. Firstly, model no. 1 without randomness again fails to capture the E and K dependence of dynamical spectra. The uniform chain exhibits a narrow magnon “band” in the E dependence (see Fig. 4i), and sharp peaks in spin structure factor at K = −3, −1 in Fig. 4g, clearly different from the experimental measurement at low temperatures. On the other hand, models no. 2−4, with similar random parameters (see above), which can fit the thermodynamics very well (see Fig. 1) can also semi-quantitatively produce the spin spectra in Fig. 4h, j. The anomalous excitation continuum observed in experiments, including the large density of states at small energies and broad peaks in Sm vs. K, are well produced by our RTIC models no.2−4, without any further adjustment of the interaction parameters (see Fig. 4h, j).

Fig. 4: Calculated ground state and dynamical properties of PrTiNbO6.
figure 4

a, c Transverse field distribution, and effective g factors of each Pr3+ ion, together with the distribution of the Ising couplings, for the sample chain A. b Calculated local magnetic moments and related spin-up and -down states of QGRR in the chain A. d−f Corresponding plots of the sample chain B. g, h Dynamical structure factors are calculated using the uniform chain model no. 1 and the random models no. 2−4, respectively. i, j The measured energy dependence of the INS intensities (Ei = 3.7 meV), with the colored lines showing the calculated data of models no. 1−4. Error bars on INS data indicate one standard error propagated from neutron counts. A more detailed comparison between experimental and numerical data under applied fields can be found in Supplementary Note 2. k Illustration of the energy spectrum where the two-lowest levels represent a QGRR.

Full size image

Given the accurate RTIC modeling of PrTiNbO6, we can now analyze the origin of quantum Griffiths singularity in the compound. In Fig. 4, we show the ED calculations of magnetic moments along the spin chain in two typical samples, with the site-dependent parameters Jj, hjx, gj generated from the Lorentzian distribution of model no. 3. In each sample chain, we find two lowest-lying states, |E1tot› and | E2tot› that are nearly degenerate and with a small energy splitting (inner gap) proportional to the minimal absolute value amongst all random fields, min(|hjx|). This collective quasi-doublet is well separated from higher-energy states by a relatively large gap of ≥ ‹J› ∼4 K (Supplementary Note 2). Similar to Eq. (2), at low temperatures (T << ‹J›) we can reconstruct the spin-up and spin-down states out of |E1tot› and | E2tot› through a coherent superposition,

$$\left| {S_{{\mathrm{tot}}}^z = \pm 1/2} \right\rangle = \frac{1}{{\sqrt 2 }}\left( {\left| {E_1^{{\mathrm{tot}}}} \right\rangle \pm \left| {E_2^{{\mathrm{tot}}}} \right\rangle } \right).$$
(5)

The two orthonormal states, |Stotz = ±1/2›, as illustrated in Fig. 4k, have a nonzero net moment restricted within the quantum Griffiths rare region (QGRR)5,6,7 well localized in a RTIC and thus can be regarded as a collective “qubit”.

By resetting (E1tot+E2tot)/2 = 0, we write down an effective Hamiltonian for such a qubit under the external field Hz,

$${\cal{H}}_{{\mathrm{QGRR}}} = - H^xS_{{\mathrm{tot}}}^x - \mu _0H^z\mu _Bg_{{\mathrm{tot}}}S_{{\mathrm{tot}}}^z,$$
(6)

where gtot = 2|Σgj‹Stotz = ±1/2|Sjz|Stotz = ±1/2›| is the overall effective g factor, and Hx ≡ E2tot−E1tot enables quantum tunneling between spin up and down, and Stotx ≡ (|Stotz = 1/2›‹Stotz = −1/2|+|Stotz = −1/2›‹Stotz = 1/2|)/2 is highly similar to Sjx in Eq. (3). The resulting distribution of local magnetic moments, ‹mjz›± = gj‹Stotz = ±1/2|Sjz|Stotz = ±1/2›, are shown in Fig. 4b for the sample chain A and in Fig. 4e for the sample chain B, from which we find that each QGRR has nonzero moments centered at hjx ∼ 0 K (indicated by the green dashed lines). The AF correlation extends to only few neighbors, while rest spins of the chain remain disordered, i.e., ‹mjz›± ~ 0, due to the relatively large transverse field ‹hx›.

The QGRR, a correlated local cluster, contributes a nonzero net magnetic moment, i.e., gtot > 0, to the random quantum Ising chain. The QGRR is responsible for the observed Griffiths singularity7 evidenced by the low-T diverging susceptibility (Fig. 1c), as well as reflected in the broad continuous spin excitations (Fig. 4) in PrTiNbO6. The latter contributes a virtually uniform background beneath the cosine-like dispersion due to short-range AF correlations. In contrast, the model no. 1 without randomness, and thus without QGRR, cannot reproduce the singular thermodynamic behaviors (Fig. 1c, e) and the spin excitation continuum (Fig. 4g, i).

Discussion

Overall, it is very remarkable that all the measured equilibrium thermodynamic properties of PrTiNbO6 can be captured accurately by this simple, yet profound, quantum RTIC model, with a proper set of random parameters. Moreover, the model also captures the main features of the spin excitation continuum. Given the accurate microscopic spin model, we can well understand the anomalous spin excitations measured in the compound and ascribe them to the existence of QGRR.

Nevertheless, the real material is always more complicated than the effective model (see Eq. (4)), and there may exist other symmetry-allowed interaction terms and additional microscopic details in PrTiNbO6 that we did not fully take into account. For example, the asymmetric distribution of random transverse fields in model no. 4 can reproduce very well the low-temperature properties below 1 K (see Fig. 1f), suggesting some deviation from symmetric Lorenzian distribution of randomness in model no. 2 and 3 may further improve their fittings to PrTiNbO6 measurements. Besides, the inclusion of small non-RTIC interactions in PrTiNbO6 may also take effects below 1 K. To understand their effects, we have computed the RTIC model under a small random XY perturbation (Sj+Sj+1−+h.c.) and observed that the quantum Griffiths singularity is intact (see Supplementary Note 4). Therefore, we believe that the RTIC model well captures the essential quantum Griffiths physics in the strongly correlated quantum material PrTiNbO6.

Lastly, as a unique phenomenon in the random spin systems that do not have clean counterparts, the quantum Griffiths singularity found in PrTiNbO6 can also exist in other non-Kramers rare-earth random compounds with correlated ground-state CEF quasi-doublets. One such example is the triangular-lattice Ising magnet, TmMgGaO432,33,47, where similar site-mixing disorder of Mg2+/Ga3+ was reported and the randomness of the transverse field was also evidenced. Therefore, our present discovery in PrTiNbO6 sheds light on similar Pr- and Tm-based compounds with site mixing, which may constitute a family of non-Kramers quantum Griffiths magnets to be explored.

Methods

Crystal synthesization and experimental measurements

Large single crystals (∼1 cm) of PrTiNbO6 were grown by the optical floating zone technique reported previously31. The µSR data were collected on the MuSR spectrometer at the ISIS pulsed muon facility, Rutherford Appleton Laboratory, United Kingdom, with the incident muon beam along the b axis (perpendicular to the Ising direction) between 0.09 and 4.2 K in a dilution refrigerator48. The additional data between 2 and 40 K were collected by transferring the sample to a 4He cryostat. Moderate-high-energy INS data were collected on the MERLIN spectrometer at the ISIS pulsed thermal neutron facility, Rutherford Appleton Laboratory, U.K., at 5 K using the powder samples of PrTiNbO6 (10.91 g) and LaTiNbO6 (10.85 g). Two incident energies, Ei = 12 and 120 meV were used with the instrumental resolutions of ∼0.5 and 5 meV, respectively49.

Low-energy-transfer INS measurements were carried out on a cold neutron multi-chopper spectrometer LET at the ISIS pulsed neutron and muon source. Incident energies of 22.5, 7.5, 3.7, and 2.2 meV were used with the energy resolutions of 1, 0.22, 0.11, and 0.04 meV, respectively50. The sample temperature down to 0.1 K was achieved using a dilution refrigerator. Magnetic fields of 0, 2, 5, and 8.8 T were applied along the c axis. All datasets in Fig. 3 were integrated over the momentum space, −0.3 ≤ L ≤ 0.3 in [0, 0, L], the data in Fig. 3b–e were further integrated over the momentum space, −0.3 ≤ H ≤ 0.3 in [H, 0, 0]. Moreover, the data at Γ and M points in Fig. 3e were further integrated over −2.1 ≤ K ≤ −1.9 and −1.1 ≤ K ≤ −0.9, respectively, in [0, K, 0]. All measured momenta were integrated over in Fig. 4i, j. The measured dynamic spin-correlation function is obtained from the raw intensity as43,

$$S_{\mathrm{m}}\left( {{\mathbf{Q}},E} \right)\sim \frac{{k_{\mathrm{i}}}}{{k_{\mathrm{f}}\left| {f\left( Q \right)} \right|^2}}I\left( {{\mathbf{Q}},E} \right),$$
(7)

where ki and kf are the incident and final neutron wave vectors, and f(Q) is the magnetic form factor of Pr3+ in the dipole approximation. For the sake of clarity, the raw data (Fig. 3a, b, e) have been symmetrized using the space group symmetry of PrTiNbO6 in the reciprocal lattice space (see Fig. 3c, d).

Quantum many-body calculations

We employ TRG approach for the accurate computation of equilibrium thermodynamics and spectral function of the RTIC model. Thermal TRG approach34,51,52 is exploited to obtain the matrix product operator (MPO) representation of the thermal density matrix \(\rho \left( \beta \right) = e^{ - \beta {\cal{H}}}\), from which the thermodynamic properties, e.g., Cm and χ, can be evaluated. We introduce the distribution of couplings and fields along the chain (with length up to l = 1000), and the observables are evaluated over 20 independent chain samples that have led to fully converged results of RTIC (Xjcal). In both thermodynamic and dynamic calculations, up to D = 200 bond states are kept in MPO, which corresponds to negligible (machine-precision) truncation errors throughout and thus assuring a very high accuracy. We fit the reported thermodynamic data31 by minimizing the following loss function,

$$R_p = \sqrt {\frac{1}{{N_0}}\mathop {\sum}\limits_{j = 1}^{N_0} {\left( {\frac{{X_j^{{\mathrm{obs}}} - X_j^{{\mathrm{cal}}}}}{{\sigma _j^{{\mathrm{obs}}}}}} \right)^2} } ,$$
(8)

where N0, Xjobs, and σjobs are the sample number, measured value, and standard deviation of the experimental data, respectively.

Given the thermal density matrix ρ(β), we perform real-time evolution and compute the dynamical spin structure factor53. We employ time-dependent algorithm at finite temperature to calculate the time-dependent correlation

$$\left\langle {S_j^z(t)S_0^z} \right\rangle _\beta = {\mathrm{Tr}}\left[ {g_jS_j^z \cdot \left( {e^{ - it{\cal{H}}}g_0S_0^ze^{ - \beta {\cal{H}}}e^{it{\cal{H}}}} \right)} \right],$$
(9)

between Sjz (at time t) and S0z at time 0 and center of the chain. The last four terms are bracketed (computed first), so as to eliminate unnecessary entanglement accumulation in the course of time evolution. Given that, the dynamical spin structure factor can be computed via the Fourier transformation

$$S\left( {K,\omega } \right) = \mathop {\sum}\limits_{j = 1}^l {e^{ - i\pi Kj}} \mathop {\int}\nolimits_{\! - \infty }^\infty {\left\langle {S_j^z(t)S_0^z} \right\rangle _\beta e^{i\omega t}dt} .$$
(10)

In practice, we perform dynamical computations on spin chains with length = 48, over 40 samples. The real-time evolution is up to time t/J = 10, guaranteeing a sufficient energy resolution below the instrumental resolution. The simulations are performed on spin chains with open boundary conditions, which, rendering no essential finite-size/boundary effects in the calculated spectral data of RTIC, as a short correlation length ξb << lb/2 is observed up to t/J ≤ 10.

Besides, we have also carried out the ED calculations of 12 sites under periodic boundary conditions (see model no.3), which also shows good consistency with large-scale many-body calculations. The international system of units is used throughout this paper, and ‹› represents thermal and/or sample average.