Revealing three-dimensional quantum criticality by Sr-substitution in Han Purple

Classical and quantum phase transitions (QPTs), with their accompanying concepts of criticality and universality, are a cornerstone of statistical thermodynamics. An exemplary controlled QPT is the field-induced magnetic ordering of a gapped quantum magnet. Although numerous"quasi-one-dimensional"coupled spin-chain and -ladder materials are known whose ordering transition is three-dimensional (3D), quasi-2D systems are special for several physical reasons. Motivated by the ancient pigment Han Purple (BaCuSi$_{2}$O$_{6}$), a quasi-2D material displaying anomalous critical properties, we present a complete analysis of Ba$_{0.9}$Sr$_{0.1}$CuSi$_{2}$O$_{6}$. We measure the zero-field magnetic excitations by neutron spectroscopy and deduce the magnetic Hamiltonian. We probe the field-induced transition by combining magnetization, specific-heat, torque and magnetocalorimetric measurements with low-temperature nuclear magnetic resonance studies near the QPT. By a Bayesian statistical analysis and large-scale Quantum Monte Carlo simulations, we demonstrate unambiguously that observable 3D quantum critical scaling is restored by the structural simplification arising from light Sr-substitution in Han Purple.


I. INTRODUCTION
At a continuous classical or quantum phase transition (QPT), characteristic energy scales vanish, characteristic ("correlation") lengths diverge and hence the properties of the system are dictated only by global and scaleinvariant quantities [1]. The critical properties at the transition then depend only on factors such as the dimensionality of the system, the symmetry group of the order parameter, and in some cases on topological criteria. These fundamental factors are not all independent, but are all discrete, and as a result phase transitions can be categorized by their "universality class." An instructive example of the effects of dimensionality is found for systems where the U(1) symmetry of the order parameter is broken [2]. For free bosons, the symmetry-broken phase is the Bose-Einstein condensate (BEC) [3,4] and, while most familiar in three dimensions (3D), this transition can in fact be found in systems whose effective dimension is any real number d > 2. However, strictly in 2D the physics is quite different, dependent on the binding of point vortices, and the system displays the Berezinskii-Kosterlitz-Thouless (BKT) transition [5,6]. In real materials, the issue of system dimensionality is complicated by the fact that the coupling of different subsystems (such as clusters, chains, or planes) may occur on different low energy scales that are only revealed as the temperature is lowered. Strict 1D or 2D behavior is avoided because, as the critical point is approached on the low-dimensional subsystem, the grow-ing correlation length causes the subsystems to become entangled and the critical properties are expected to be those associated with a 3D universality class [7,8].
Dimerized quantum spin systems, composed of strongly coupled pairs of spin-1/2 ions, constitute a particularly valuable class of materials for the study of quantum criticality. For sufficiently strong dimerization, the zero-field ground state of local singlets is always gapped, with the excitations taking the form of propagating local triplet quasiparticles ("triplons"). In an applied magnetic field, this system undergoes a QPT at which the gap is closed and the new ground state has field-induced transverse magnetic order, which can be described as a triplon condensate [9][10][11][12]. This is a 3D QPT from a quantum disordered phase with U(1) (or XY) symmetry to a BEC phase, whose emerging long-range magnetic order breaks the U(1) symmetry. This QPT has been studied extensively in materials where the inter-dimer coupling is similar in all three spatial dimensions [11,13]. In gapped quasi-1D materials (dimerized chains, Haldane chains, two-leg ladders), fields exceeding the gap reveal the physics of the Tomonaga-Luttinger liquid at higher temperatures, but 3D scaling behavior is restored in the quantum critical regime around the QPT [14][15][16][17][18][19]. However, 2D systems with continuous symmetry are special because of the Mermin-Wagner theorem [20], and extraspecial because of BKT physics [5,6], which has recently been argued to give rise to behavior at the dimensional crossover that is radically different from the 0D and 1D cases [21]. Thus the question arises as to whether the "conventional" emergence of 3D physics can really be expected at the field-induced QPT in a quasi-2D material.
Unfortunately the list of candidate quasi-2D spindimer compounds suitable for answering this question is rather short.
A further prominent material, which does actually realize 2D square lattices of dimers and a field-induced triplon condensation at µ 0 H c1 23.5 T, is the ancient pigment Han Purple (BaCuSi 2 O 6 ) [43,44]. In this compound the stacked dimer bilayers have a body-centered offset, leading to the expectation of an exact frustration of AF coupling between successive bilayers, and further to the interpretation of thermodynamic measurements showing apparent 2D scaling close to the QPT [45] as a frustration-induced "dimensional reduction" quite different from 3D critical properties. However, it was shown subsequently that the weakly orthorhombic structure of this material below 90 K [46][47][48] contains three types of structurally [49] and magnetically [50][51][52] inequivalent bilayers, and later that the inter-bilayer coupling is not in fact frustrated [53,54]. These properties lead ultimately to a regime of anomalous effective scaling in the experimentally accessible parameter space [54], with the 3D critical behavior remaining unobservable. At room temperature, BaCuSi 2 O 6 has a tetragonal structure (space group I4 1 /acd) in which all dimer bilayers are equivalent. Thus a possible route to the ideal quasi-2D material would be to suppress the orthorhombic transition. In pursuing this program, a number of stoichiometric substitutions have been tested on the nonmagnetic Ba site, and it was reported [55] that even a 5% Sr substitution is sufficient to stabilize the tetragonal structure, represented in Fig. 1(a), down to the lowest temperatures. However, it is a significant challenge to grow single crystals of the substituted material [56], and particularly to grow large single crystals suitable for inelastic neutron scattering (INS) experiments [57]. We have solved these problems, enabling the detailed study of the ground state, excitations, and field-temperature phase diagram of Ba 0.9 Sr 0.1 CuSi 2 O 6 that we present.
The structure of our contribution is as follows. In Sec. II we detail our crystal-growth procedures, experimental methods, and numerical simulations. In Sec. III we report INS experiments at zero magnetic field, which demonstrate that this compound contains only one type of dimer bilayer, and hence one triplon excitation branch. We determine the minimal spin Hamiltonian, showing that it contains only the three Heisenberg interactions of Fig. 1(b), which form a clear hierarchy of intradimer (J), intra-bilayer (J ), and inter-bilayer (J ) coupling strengths. Turning to high magnetic fields, in Sec. IV we perform magnetization, magnetic torque, specific-heat, magnetocaloric effect (MCE), and nuclear magnetic resonance (NMR) measurements to obtain a comprehensive experimental picture of the field-temperature phase boundary to the magnetically ordered state. In Sec. V we combine state-of-the-art statistical data analysis and Quantum Monte Carlo (QMC) simulations based on the INS parameters to interpret our results as demonstrating complete consistency with 3D scaling in the quantum critical region. In Sec. VI we discuss the context of these results, particularly regarding possible magnetic disorder effects, and conclude in Sec. VII that the structural consequences of weak chemical substitution by non-magnetic ions are sufficient to realize the simplicity and elegance of the model originally envisaged for Han Purple, namely (e) (f) Figure 1. Structure and inelastic neutron spectrum of Ba0.9Sr0.1CuSi2O6. (a) Tetragonal crystallographic structure of Ba0.9Sr0.1CuSi2O6 [55]. (b) Representation of the minimal magnetic unit cell, indicating the relevant Heisenberg interaction parameters between dimers on the two bilayers, which are structurally equivalent with body-centered geometry. (â,b,ĉ) and (â ,b ,ĉ ) denote respectively the axes of the crystallographic and minimal magnetic unit cells [54]. Dotted grey lines indicate the correspondence between Cu 2+ ion positions in the two cells. (c)-(d) Neutron scattering spectra measured on CAMEA at 1.5 K for two high-symmetry Q-space directions; Q is indexed in reciprocal lattice units (r.l.u.) of the crystallographic unit cell.
(e)-(f) Spectra modelled using the interaction parameters and linewidths extracted from the data collected on both CAMEA and TASP (Fig. 2).
a quasi-2D spin-dimer material displaying observable 3D quantum critical behavior.

A. Crystal growth
Single-crystal growth of BaCuSi 2 O 6 poses a challenge because the starting CuO component decays before the reagents and compound melt. This decay has been suppressed by increasing the oxygen partial pressure in optical floating-zone growth [43] and by using oxygen-donating LiBO 2 flux [45]. However, the latter method cannot be applied to the substituted system, Ba 0.9 Sr 0.1 CuSi 2 O 6 , because the metaborate flux reacts with the SrO reagent [56]. Thus two methods were employed to grow the single crystals used in the present work. All of the high-field experiments were performed using small single crystals (sizes of order 10 mg) grown from the congruent melt in a tube furnace with an oxygen flow at a pressure under 1 bar [56]. The INS experiments were performed using one rod-shaped single crystal of mass 1.3 g grown by a floating-zone method. While this method was developed [57] for BaCuSi 2 O 6 , it could not be applied directly to Ba 0.9 Sr 0.1 CuSi 2 O 6 because large oxygen bubbles from the reduction of CuO disconnected the two rods and terminated the crystal growth. In this sense, obtaining several gram-sized single crystals of Ba 0.9 Sr 0.1 CuSi 2 O 6 constituted a significant achievement, made possible by using one single-crystalline seed in a mixed argon/oxygen atmosphere to reduce the sizes of the disruptive oxygen bubbles.

B. Inelastic neutron scattering measurements
INS measurements at zero applied magnetic field were performed on the triple-axis spectrometers TASP and on the multiplexing triple-axis spectrometer CAMEA [58,59], both at the SINQ neutron source [60] at the Paul Scherrer Institute. INS probes the dynamical mag- netic structure factor, S( Q, ω), which is the spatial and temporal Fourier transform of the spin-spin correlation function, allowing direct access to the magnetic excitation spectrum as a function of the energy ( ω) and momentum ( Q) transfer of the scattered neutrons [61]. On TASP, a temperature of 1.5 K was used and the final neutron energy was fixed to 3.5 meV. The instrumental parameters of a previous study of BaCuSi 2 O 6 [54] were selected to optimize the energy resolution. On CAMEA, a temperature of 1.5 K was used and sequential measurements at fixed incident neutron energies of 7.5 and 9 meV were combined to obtain a map of the spectrum over multiple Brillouin zones and over an energy range from 2.5 to 5.5 meV. For both measurements, the rod-shaped single crystal was aligned with (h 0 l) in the scattering plane [notation from the crystallographic unit cell, Fig. 1(c)] [55]. The CAMEA data shown in Figs. 1(c)-1(f) were extracted from the raw measured intensities using the software MJOLNIR [62], which corrects these intensities using measurements on a standard vanadium sample. The two datasets shown in Fig. 2

C. High-field measurements
Magnetization. Magnetization measurements up to 50 T were carried out in a capacitor-driven short-pulse magnet at the NHMFL Pulsed Field Facility in Los Alamos using an extraction technique [63]. These measurements are subject to a strong reduction of the sample temperature in the changing magnetic field, which was calibrated by explicit MCE measurements.
Specific heat. The specific heat was measured with a heat-pulse technique inside a motor-generator-driven long-pulse magnet at the International MegaGauss Science Laboratory of the University of Tokyo [64]. The experimental data points were collected during the flat-topped portion of the magnetic-field profiles created by feedback loop control [64,65], at fields up to 37.8 T and temperatures down to 1 K.
Magnetic torque. Measurements of the magnetic torque were conducted in a 31 T resistive magnet at the NHMFL DC Field Facility in Tallahassee using a silicon-membrane cantilever manufactured by NanoSensors TM and a 3 He cryostat to reach temper-atures down to 0.3 K. Torque measurements were made during constant-temperature field sweeps across the phase boundary. The temperature was fixed by maintaining a constant 3 He pressure in the sample space; two field-calibrated sensors, one located close to the cantilever and the other at the sample holder, were used for a precise in situ monitoring of the temperature. The angle between the field and the crystallographic axes of the sample was estimated by rotation tests to be 10.5 • .
Magnetocaloric effect (MCE). MCE measurements were performed in vacuum at temperatures down to 1 K at the International MegaGauss Science Laboratory of the University of Tokyo [66], and down to 0.4 K at the NHMFL in Los Alamos. In both experiments, the temperature was recorded during magnetic-field pulses up to 50 T using the response of a thin-film AuGe thermometer sputtered onto the sample surfaces. Under adiabatic vacuum conditions, the lines defined by T (H) may be understood as isentropes. The sign of the temperature change is given by the Maxwell relation (∂Q/∂H)/T = −(∂M/∂T )| H ; the changes of sign in the gradient of the isentropes correspond to the maxima or minima in M (T ) that occur at the ordering transitions.
Nuclear Magnetic Resonance (NMR). 29 Si NMR measurements were performed in the field range 22-26 T using a resistive magnet at LNCMI. A 13 mg singlecrystal sample (5×1×0.5 mm 3 ) was placed in the mixing chamber of a 3 He-4 He dilution refrigerator with its c-axis oriented parallel to the applied field, H. NMR spectra were taken by the spin-echo pulse sequence. In the quantum disordered phase, the spectra are relatively narrow and were covered fully in a single-frequency recording; inside the magnetically ordered phase, the spectra are broader and their reconstruction required several recordings taken at equally spaced frequency intervals, by summation of the individual Fourier transforms. The frequency (ν) of the 29 Si spectra was measured relative to the reference ν 0 /(µ 0 H) = 29 γ = 8.4577 MHz/T. The magnitude of the applied field was calibrated using a metallic 27 Al reference placed in the same coil as the sample. The same 27 Al reference was also used to obtain a relative temperature calibration at temperatures below 1 K from the Boltzmann-factor dependence of the NMR signal intensity; this was then related to the values measured directly by the field-calibrated RuO 2 temperature sensor at all higher temperatures, where stable and fieldindependent T values were ensured by fixing the pressure of the He bath.

D. Quantum Monte Carlo
We performed QMC simulations by applying the stochastic series expansion (SSE) algorithm [67] on 3D lattices of L×L×L/2 sites for L = 8, 10, 12, . . . , 32. To reduce the computational cost we simulated an effective Hamiltonian of hard-core bosons [68], where each site represents a dimer and each site boson a dimer triplon. The dispersion relation of these bosonic quasiparticles was chosen to match the triplon band determined by INS.
For each applied field, the critical temperature, T c , was extracted from the finite-size scaling of the superfluid density, ρ sf , of the hard-core bosons (equivalent to the spin stiffness of an ordered magnet). The quantity L d+z−2 ρ sf (L, T ) has scaling dimension zero, where the dimensionality d = 3 and the dynamical exponent z = 0 for a finite-temperature transition, and hence the Lρ sf (L, T ) curves for all L cross at the same point, which is T = T c [69]. We verified that the values of T c (H), and the associated uncertainties δT c (H), determined in this way were fully consistent with the critical exponent, ν = 0.672, expected for the correlation length [54].

III. MAGNETIC EXCITATIONS AT ZERO FIELD
The first step in characterizing the magnetic excitations of Ba 0.9 Sr 0.1 CuSi 2 O 6 was to succeed in growing gram-sized single crystals by the technique described in Sec. IIA.  Fig. 1(b)], as proposed [53] and observed [54] in BaCuSi 2 O 6 . Similarly, the location of the band minimum in Q l is a consequence of a FM inter-bilayer interaction parameter, J , which is the magnetic coupling between the top ions of dimers in one bilayer and the bottom ions of the dimers in the next [ Fig. 1(b)]. We remind the reader that the FM intra-bilayer correlations ensure that the inter-bilayer interactions are not frustrated for either sign of this latter parameter, excluding [53] the "dimensional reduction" scenario mentioned above. The leading qualitative feature of the intensities is the modulation visible in Figs. 1(c) and 1(e), which is a consequence of the fact that the resolution ellipsoid on CAMEA causes intensity to be spread over a narrower (focusing) or wider (de-focusing) range depending on the value of Q.
Turning to a quantitative analysis of the lowtemperature spectrum, Fig. 2(a) shows the scattered intensity as a function of energy transfer for Ba 0.9 Sr 0.1 CuSi 2 O 6 (blue symbols), and also the comparison with BaCuSi 2 O 6 (red) [54], which we discuss in Sec. IIIB; both data sets were measured on TASP at Q = (0 0 4) using similar experimental conditions. The solid lines are Gaussian fits used to extract the location, linewidth, and intensity of the triplon modes. The results of applying the same procedure at all Q points along the same two high-symmetry directions as in Figs

A. Spin Hamiltonian
We stress that our observation of one triplon mode per minimal magnetic unit cell means that the Ba 0.9 Sr 0.1 CuSi 2 O 6 structure contains only one type of bilayer, in contrast to the complicated situation in BaCuSi 2 O 6 [54]. An accurate determination of the spin Hamiltonian is crucial to verify the quasi-2D characteristics of the material and for a detailed investigation of the resulting critical behavior in an applied magnetic field. We assume that the minimal magnetic Hamiltonian of Ba 0.9 Sr 0.1 CuSi 2 O 6 is that displayed in Fig. 1(b), i.e. that it contains only three interactions, all of purely Heisenberg type, to which, as above, we refer as intradimer (J), intra-bilayer (J ), and inter-bilayer (J ). We note that J is an effective inter-dimer interaction parameter that combines two pairs of ionic (Cu 2+ -Cu 2+ ) interactions within and between the two layers of the bilayer structure [53,54].
By fitting the dispersion of the triplon mode shown in Figs. 2(b) and 2(c) as detailed in App. A, we deduce the optimal parameter set J = 4.28(2) meV, J = −0.52(1) meV, and J = −0.02(1) meV, where negative values denote FM interactions and the error bars represent a conservative estimate of the statistical uncertainties. Figures 1(e), 1(f), and 2(b)-2(e) make clear that these parameters provide an excellent account of all the measured data, leaving no evidence of any requirement for further terms in the spin Hamiltonian. The INS parameters also serve as a benchmark for the magnetic interactions estimated using electronic structure calculations based on measurements of the atomic structure [55], from which one may deduce (using the notation in Table II of that work) that J 1 ≡ J is obtained to very high accuracy, J 3 − J 5 ≡ J is overestimated by approximately 30%, J 4 = 0 is identified correctly, and J 2 ≡ J is surprisingly accurate in both size and sign, given its small value. On this point, we stress that the three terms we have deduced form a hierarchy of interactions whose strengths differ by at least an order of magnitude, and hence that each is responsible for physics at quite different energy scales. Specifically, from the standpoint of triplon excitations, Ba 0.9 Sr 0.1 CuSi 2 O 6 is very much a quasi-2D magnetic material, with an interbilayer coupling 25 times smaller than the energy scale characterizing the intra-bilayer excitation processes.

B. Comparison to BaCuSi2O6
It is instructive to compare the magnetic interactions of the Ba 0.9 Sr 0.1 CuSi 2 O 6 system with those of the parent compound. The transition to a weakly orthorhombic structure below 90 K [46][47][48] in BaCuSi 2 O 6 results in the formation of three structurally inequivalent bilayers of dimers [49,51]. This led to the observation in INS experiments on BaCuSi 2 O 6 of either three separate triplon modes or two separate peaks at Q vectors where the energies of the upper two modes are not distinguishable, which is the case in the data shown in Fig. 2(a) for Q = (0 0 4).
In Ba 0.9 Sr 0.1 CuSi 2 O 6 , the absence of this structural transition [55] means that the system remains tetragonal at low temperatures and only one bilayer type, hosting only one triplon mode, is expected. The neutron spectra presented in Figs. 1(c)-1(f) and 2 confirm this situation. If the tetragonal interaction parameters are compared with the orthorhombic ones, J = 4.28 (2) [54]. This is a consequence of the ABC bilayer stacking in BaCuSi 2 O 6 , by which the three different dimer types cause a dramatic reduction of the interbilayer triplon hopping probability.
In principle, valuable physical information is also contained in the line widths of the triplon modes observed in the two materials. Partial Sr substitution on the Ba sites should result in a degree of bond disorder, meaning a distribution in the values of the interaction parameters, which would act to broaden the triplon in Ba 0.9 Sr 0.1 CuSi 2 O 6 . Comparison with the A triplon in BaCuSi 2 O 6 , which lies at the same energy in Fig. 2(a), reveals that the Full Width at Half Maximum (FWHM) of the Gaussian characterizing the 10%-substituted material is 20% higher at Q = (0 0 4) [0.32(2) meV as opposed to 0.25(2) meV]. At Q = (1 0 4), however, the FWHM is essentially unchanged [0.30(2) meV as opposed to 0.29(4) meV, not shown]. A number of factors affect this comparison, including that the Sr-substituted crys-tal was subject to a grain reorientation during growth [57] and that the different shapes of the two crystals made it necessary to adjust the focusing slits for each case, affecting the experimental resolution. On the basis of these results, we are not able to discern a significant intrinsic broadening, due to bond disorder, from the extrinsic contributions to broadening, and thus we make only the qualitative statement that disorder effects on the magnetic properties arising from the weak structural disorder caused by Sr-substitution are small in Ba 0.9 Sr 0.1 CuSi 2 O 6 .

IV. FIELD-TEMPERATURE PHASE BOUNDARY
The parameters we have determined for the gapped, zero-field ground state have essential consequences for the field-temperature phase diagram. The onset of fieldinduced magnetic order, when the field is strong enough to close the triplon gap, is generally referred to as a Bose-Einstein condensation of magnons; while the quadratic magnon dispersion at the transition ensures the Bose-Einstein universality class, the effective dimensionality of the system remains an open issue, as explained in Sec. I. Given the questions raised by the parent compound BaCuSi 2 O 6 , and the more general questions raised about possible fingerprints of BKT physics in coupled quasi-2D systems [21], we have performed the most comprehensive study within our capabilities of the QPT and the surrounding quantum critical regime.

A. High-field thermodynamic measurements
For this we applied a combination of thermodynamic measurements up to the highest available pulsed and static magnetic fields (Sec. IIC) to map the fieldtemperature phase diagram of Ba 0.9 Sr 0.1 CuSi 2 O 6 . Magnetization data obtained in pulsed fields are shown in Fig. 3(a). While zero magnetization indicates a spin gap, an abrupt jump and increasing M (H) values are a signature of gap closure and magnetic polarization. The phase transitions both into and out of the ordered state are identified as sharp features in the second field derivative [inset, Fig. 3(a)]. Although field pulses up to 50 T make it possible to reach saturation even at the lowest temperatures, the M (H) curves taken in short-pulse magnets are adiabatic, rather than isothermal. BaCuSi 2 O 6 and Ba 0.9 Sr 0.1 CuSi 2 O 6 have a strong MCE, meaning that the rapidly changing applied field drives a dramatic reduction of the sample temperature in the absence of heat transfer. We performed explicit MCE measurements (below) in order to obtain an accurate determination of the real measurement temperatures shown for both phase transitions in Fig. 3(a).
We have measured the specific heat in pulsed fields up to 37.8 T and display the results in Fig. 3(b). Be-cause this field value is closer to the upper critical field of Ba 0.9 Sr 0.1 CuSi 2 O 6 than to the lower, these measurements make it possible to access the top of the phaseboundary "dome" shown in Fig. 3(d). The λ-shaped anomalies in C p (T ) are characteristic of second-order phase transitions and we used an equal-area method to estimate the critical temperatures for each magnetic field. The λ shapes have also been associated directly with the 3D XY universality class in a number of systems showing magnon BEC [40,41,43,44,[70][71][72], while below the transition C p (T ) is consistent with the T 3 form expected from AF magnons in 3D.
In a static external field, we have measured the magnetic torque up to 29 T, as shown in Fig. 3(c). In these experiments the sample temperature is known exactly but, because a torque acts only when the sample is tilted away from its principal axes, the tilt angle provides a degree of uncertainty that affects the absolute field. By contrast, the relative values of the critical fields at each temperature are extremely accurate, allowing one to probe the shape of much of the left side of the phaseboundary dome. Again the second field-derivative [inset Fig. 3(c)] gives the most accurate determination of the critical field, H c , for each measurement temperature, T .
Here we have corrected these fields to the NMR values of H c (T ) (below) in the regime of temperature where the two methods overlap, thereby obtaining the phaseboundary points shown up to 3.03 K in Fig. 3(d).
The MCE, meaning the change in sample temperature due to an applied field, is related directly to the first temperature derivative of the magnetization (Sec. IIC). MCE measurements in pulsed fields up to 50 T were used to obtain the lines of constant entropy shown in Fig. 3(d), whose local minima mark the location of the ordering transition in both field and temperature. We carried out two separate MCE experiments, at ISSP Tokyo and NHMFL Los Alamos, and used the latter data to deduce the true sample temperature for the magnetization results shown in Fig. 3(a). The MCE shows a discernible asymmetry [ Fig. 3(d)], in that the phase-boundary temperature for the same isentrope is higher on the left side of the dome than on the right [42], indicating that the isothermal entropy, S(H), is higher on the right than on the left. This result implies a larger specific heat on the right side of the dome and has been explained as a magnon mass anisotropy [71] that is connected to the ratio of the upper and lower critical fields. In Fig. 3(d) we show also the magnetic Grüneisen parameter, whose sign change becomes increasingly abrupt as the temperature is lowered, thereby locating the phase boundary with increasing accuracy.
We also performed NMR measurements at very low temperatures to determine the phase boundary close to the lower critical field, but defer a detailed discussion of these experiments to the next subsection. Figure 3 is rather symmetrical in shape between a lower critical field µ 0 H c1 22.5 T and an upper critical field µ 0 H c2 46.5 T, with a maximum height of T max c 3.8 K around 35 T. In Fig. 3(d) we show also the phase-boundary dome determined [43] for the parent compound, BaCuSi 2 O 6 , which is altered only slightly by the Sr substitution. Specifically, the small reduction of H c1 is a consequence of the smaller gap in Ba 0.9 Sr 0.1 CuSi 2 O 6 , which arises due to the combination of increased bandwidth (controlled by J ) and rather constant band center (controlled by J).
The reduction of the upper critical field, H c2 , is a consequence of the fact that the stronger intra-dimer coupling constants of the structurally inequivalent bilayers in BaCuSi 2 O 6 are absent in Ba 0.9 Sr 0.1 CuSi 2 O 6 .

B. Nuclear Magnetic Resonance
For a detailed investigation of the quantum critical regime, we performed NMR measurements at high magnetic fields and very low temperatures using the facilities at the LNCMI (Sec. IIC). In the magnetically ordered phase above the critical field, which at any temperature we denote by H c (T ), the presence of a staggered magnetization, m ⊥ (T ), transverse to the direction of the external field (which is applied along theĉ-axis) lifts the degeneracy between sites that were equivalent in the gapped dimer phase [73,74]. Their NMR spectra therefore show a splitting into two separate peaks. Figure 4(a) shows 29 Si NMR spectra obtained using a 13 mg single crystal of Ba 0.9 Sr 0.1 CuSi 2 O 6 at T = 120 mK over a range of static magnetic fields crossing the value of H c for this temperature. Whereas in BaCuSi 2 O 6 it was possible by 29 Si NMR to detect only a broad distribution due to the presence of structurally different bilayer types [48,50,52], whose local spin polarization could be determined quantitatively only by 63 Cu NMR [75], in the 29 Si spectra of the Sr-substituted variant we observe a second peak that splits from the first with increasing applied field, becoming a separate and clearly resolved entity in the upper panels of Fig. 4(a).
Because of the low temperatures of these measurements, the field value at which the NMR peak-splitting appears constitutes our most accurate estimate of the critical field, H c (T ), for each temperature, T , in the regime most likely to approach quantum critical scaling. The splitting of the peaks in frequency is directly proportional to the magnetic order parameter [73], allowing us to obtain the field-induced evolution of the ordered moment up to a single, non-universal constant of proportionality, which is the hyperfine coupling constant of the NMR experiment. To determine the fielddependence of the peak-splitting, and hence of the order parameter, we measured NMR spectra up to fields significantly higher than H c (T ) for one example temperature, 120 mK [ Fig. 4(a)], and the results are shown in the inset of Fig. 4(d). To determine the phase boundary for a study of its quantum critical properties, we have measured NMR spectra for a number of fixed temperatures between 49 and 921 mK, over a range of fields close to H c (T ) in each case, and a further set of spectra is shown for T = 412 mK in Fig. 4(b).
To optimize our extraction of the peak-splitting from the spectral data, which contain multiple additional features and complications, we exploit its relation to the order parameter to model it as a power-law function of the field difference from the phase transition, H − H c (T ), and perform a simultaneous fit of one peak (blue) or two split peaks [red and green in Figs. 4(a) and 4(b)] to all NMR spectra measured at the same T . Details of this procedure, which is valid over a field-temperature regime near the phase boundary, are provided in App. B. Only at fields well above H c (T ) was it possible to obtain reliable fits to the two clearly separated spectral peaks of each individual spectrum, and examples of these fits are shown in the upper panels of Fig. 4(a).
The result of the simultaneous fit procedure is the set of fitting functions shown by the lines in Fig. 4(c), where the crosses represent the fields at which spectra were collected, and these lines yield a set of estimated critical points, H c,j , with uncertainties δH c,j . The corresponding temperature values, with their uncertainties, were determined relative to a field-calibrated RuO 2 thermometer (Sec. IIC). The data pairs (H c,j , δH c,j ) and (T c,j , δT c,j ) together constitute the NMR phase boundary near the quantum critical point, which we show as the blue and open circles in Fig. 4(d). The blue line and shading represent the results of our statistical analysis of the quantum critical properties of the NMR phase boundary, which we discuss below.
To assist with the interpretation of the measured phase boundary, we performed large-scale QMC calculations of a model with only one type of bilayer, following the established procedure [54] summarized in Sec. IID. For this we mapped the spin Hamiltonian of Fig. 1(b) to a hard-coreboson model [68] for triplons with the low-energy triplon dispersion determined from our zero-field INS measurements. The QMC phase boundary in the regime close to the quantum critical point (QCP) is shown by the orange and open squares in Fig. 4(b), while the corresponding lines and shading are the results of the statistical analysis we present next (Sec. V). We remark that the precise agreement in location of the QCP obtained by using the zero-field parameters implies that magnetostriction effects are negligible, a result also found in BaCuSi 2 O 6 . We also performed QMC simulations to determine the order parameter at low temperatures for fields well inside the magnetically ordered regime, and in the inset of Fig. 4(d) these are shown for a temperature corresponding to 120 mK. By matching the field-induced ordered moment in Bohr magnetons to the peak-splitting measured in MHz over this field range, we obtain the hyperfine coupling constant as noted above.

V. CRITICAL EXPONENT
From the statistical mechanics of classical phase transitions, the field-temperature phase boundary should be described by a power-law form in the critical regime close to the QCP, and hence can be expressed as Here H c1 is the critical field of the QCP, i.e. at zero temperature, α is a nonuniversal constant of proportionality and φ is the critical exponent, which is universal and  29 Si NMR spectra measured for a range of field values around the corresponding value of Hc(T ) at 120 mK (a) and 412 mK (b). The spectra show a splitting from a single, sharp peak (blue) in the quantum disordered phase at H < Hc(T ) to two generally broader peaks (red and green) in the magnetically ordered phase at H > Hc(T ). The measured intensities were modelled as the sum (grey) of the single or double peak(s) and a weak, constant contribution. The double peaks were modelled by a simultaneous fitting procedure in all panels other than the highest four fields shown at 120 mK (a), where individual peak fits were used. (c) Splitting of peaks in the NMR spectra. Solid lines show splitting functions deduced from simultaneous fits to all spectra measured at the same temperature, and crosses mark the field values at which data were taken. (d) Low-field phase boundary of Ba0.9Sr0.1CuSi2O6 extracted from the NMR peak-splitting data (blue and white circles). Orange and white squares show the phase boundary obtained from QMC simulations. Blue and orange lines show the optimal fit to the phase boundary based on each set of points, obtained from the mean of the posterior distribution determined by a Bayesian inference procedure based on Eq. (1). The blue and orange shaded areas represent the estimated uncertainties in the form of the 68% (dark) and 95% (light shading) credible intervals (CIs). NMR and QMC data points marked by white circles and squares are those excluded from the statistical analysis of the quantum critical properties. Inset: NMR peak-splitting measured up to 25 T at T = 120 mK (blue circles) juxtaposed with the magnetic order parameter, m ⊥ (T ), obtained from QMC simulations.
should reflect only the dimensionality of space and the symmetry of the order parameter. In a true BEC one has φ = 2/d [45,[76][77][78], where the numerator is the dynamical exponent for free bosons (ω ∝ k 2 ) and d the dimensionality of space. The form of Eq. (1) means that the parameters (α, H c1 , φ) required to model the experimental data are highly correlated; specifically, the best estimate for φ depends strongly on the value of H c1 . Further, the width of the quantum critical regime, the ranges of T c and H − H c1 over which Eq. (1) remains valid, is a non-universal quantity about which little specific information is available. Attempts to deal with both the correlation problem [45] and the critical-regime problem [79,80] on the basis of limited available data have led in the past to significant confusion, which has been cleared up for certain models [76]. In order to obtain an accurate estimate of φ and to quantify its uncertainty reliably, we have applied a thorough statistical analysis using Bayesian inference. Bayesian inference is a general method for the statistical analysis of all forms of experimental data, and thus its applications can be found throughout physics [81]. The quantity of interest is the posterior probability distribution, of the model parameters, θ, conditional on the observed data, D. In this expression of Bayes' theorem, p(D|θ) denotes the likelihood of making the observations D for a given θ and p(θ) is a prior probability distribution summarizing any information previously available about the model parameters. In the problem of describing the phase boundary, θ = (α, H c1 , φ) is the set of parameters in the model of Eq. (1) and D = {H c,j , T c,j , δH c,j , δT c,j } is the set of phase-boundary points with their associated uncertainties; clearly the total number of points (j) in the input information will have a leading role in determining the accuracy of the output. The posterior distribution, p(θ|D), is a high-dimensional density function that contains all of the information available about the parameters in θ, including their uncertainties and (non-linear) correlations. To obtain individual parameter estimates and error bars, the posterior distribution is compressed in various ways. We use the mean of the marginal distribution of each parameter, which provides the prediction with the smallest error given the input data [81]. To quantify these errors, we construct the highest posteriordensity credible intervals (CIs), which are the shortest intervals containing a fixed fraction of the posterior density [82], and show the commonly chosen 68% and 95% CIs. Figure 4(d) shows the mean of the posterior distributions of the model function of Eq. (1) as a blue line for the NMR dataset, D NMR , and an orange line for our QMC data, D QMC , while the shaded regions indicate the 68% and 95% CIs. Referring back to the issue of the width of the quantum critical regime, the two posterior distributions were constructed on the basis of the 6 NMR phase-boundary data points falling below the temperature T NMR max = 412 mK (blue circles) and the 17 QMC points below T QMC max = 426 mK (orange squares). The criterion by which these maximum values were determined is discussed below, and an extension of the fit to higher temperatures is shown for reference.
To explain these results, Figs. 5(a)-5(c) show the marginal distributions, p(θ i |D NMR ), of each of the three model parameters in θ. These are normalized density distributions, obtained from Monte Carlo sampling of the posterior by integrating over the remaining two parameters in each case, which allow the full shape of each marginal to be investigated. These are clearly skewed, with the distributions of α [ Fig. 5(a)] and φ [ Fig. 5(b)] having tails extending towards larger values whereas that of H c1 [Fig. 5(c)] extends towards smaller values. Qualitatively, the origin of these distorted shapes lies in the interplay of the uncertainties in the field and temperature measurements, whose simultaneous consideration is a key strength of the Bayesian approach; details of its quantitative implementation are presented in App. C. By contrast, a regression analysis such as a least-squares approach would yield biased estimates under these circumstances. Figure 5(d) shows a projection of the posterior distribution as a pair plot depending on φ and H c1 , which highlights the strong negative correlation between these two parameters. In the absence of accurate knowledge of H c1 , the NMR phase boundary data of Fig. 4(d) can often be described by decreasing H c1 while simultaneously reducing the height of the curve, i.e. increasing φ, or vice versa. In view of this situation, one may consider the precise goal of the present analysis as being to localize the most probable (H c1 , φ) pair accurately on the axes of Fig. 5(d).
Using the mean and 68% CI boundaries of the marginal posterior distributions, from our NMR data for T ≤ T NMR max we obtain the model parameter estimates α = 1.35 +0. 15 −0.16 K, µ 0 H c1 = 22.550 +0.007 −0.003 T, and φ = 0.666 +0.053 −0.074 . The non-universal parameter α has no influence on the critical properties. By comparing H c1 , for a field applied along theĉ-axis, with the value of the spin gap extracted from INS, ∆ = 3.00(1) meV, we deduce a g-factor of g c = 2.298 +0.008 −0.008 . This result is in agreement with both the value 2.32(2) determined for the Sr-doped material [55] and the value g c = 2.306(3) obtained for the parent compound [83]. Clearly the mean value of φ we obtain agrees spectacularly well with the 3D exponent of 2/3, although we stress that the CIs are broad (as one may expect from the number of data points in the analysis), and to interpret these results we turn to QMC modelling.
Because our QMC simulations were performed with an effective quasiparticle Hamiltonian matching the INS dispersion, H c1 is fixed by the lower band edge (using the same g-factor) and the set of points {H c,j } has no errors ({δH c,j } = 0), which simplifies the statistical analysis (App. C). We deduce the values α = 1.46 +0.03 −0.03 K and φ = 0.663 +0.008 −0.009 , based on our data below T QMC max = 426 mK. Figure 5(e) shows the putative critical regime on logarithmic axes, a format in which the gradient of each line is the critical exponent, φ, and in which both our NMR and QMC results are indistinguishable from the line φ = 2/3 expected of a 3D BEC.
Returning to the width of the critical regime, this is a non-universal quantity in studies of critical behavior and hence is an unknown in any fitting procedure. Under these circumstances, the windowing fit was introduced [76] to discern the trend of the data towards a critical value at the QCP when the width of the critical regime may be arbitrarily small, and has been applied [45,54] for the analysis of BaCuSi 2 O 6 data. To investigate this issue, we performed an effective windowing analysis by applying the Bayesian inference procedure using 10,9,8,7,6, and 5 NMR data points, and the results for the fitted exponents are shown as the blue symbols in Fig. 5(f). The horizontal lines show the mean and the heights of the shaded boxes show the 68% and 95% CIs of the marginal posterior distribution of φ for each case. We observe that the estimated mean value of φ is almost exactly 2/3 for the three narrowest windows, but then dips on inclusion of the 8th data point and recovers slightly towards 2/3 when the 9th and 10th points are considered. However, by a criterion that the mean φ should include 2/3 within the 68% CI, all of these windows are consistent with 3D critical scaling. For further insight as to whether it is realistic to ascribe this behavior in its entirety to quantum critical physics, we perform the same type of analysis using our QMC data. While QMC studies of 3D models have revealed quantum critical regimes with relative widths on the order of 10-20% of the control parameter [76,84], simulations of quasi-1D models have struggled to reach such a regime. Our QMC results, which have better data coverage and much narrower CIs, can be expected to reflect the physics of quasi-2D bosons at the energy scales of their interlayer coupling. The data in Figs. 5(f) show a systematic and near-monotonic evolution best described as the estimated φ converging to 2/3 from below as the number of data points is reduced. If one defines T max as the highest temperature at which φ has converged to 2/3 within the 68% CI, one obtains T QMC max = 426 mK, which corresponds to the use of 17 phase-boundary points. On this basis we deduce that effects beyond quantum critical fluctuations become discernible above approximately 0.4 K, and hence that this is an appropriate width to adopt for the scaling regime.
This criterion terminates our NMR dataset at T NMR max = 412 mK, corresponding to the 6 phase-boundary points used for the analysis shown in Figs. 4(d) and 5(a)-5(d). However, we note that the criterion is far from exact, and in the present case an almost identical estimate with lower uncertainties, φ = 0.665 +0.040 −0.054 , is obtained by including the 7th data point. The evolution of the estimated φ on inclusion of the 9th and 10th data points indicates both the onset of complex and manifestly nonmonotonic crossover behavior beyond the critical regime and the accuracy of the analysis beyond features visible in Figs. 4(d) or Figs. 5(e). The discrepancy between the forms of behavior exhibited by our NMR and QMC data could lie in interaction terms neglected in our minimal magnetic Hamiltonian, or in possible disorder effects arising due to Sr substitution, which we discuss further in Sec. VI. However, we stress that these differences in critical behavior are extremely subtle, and are certainly not discernible in the NMR and QMC phase boundaries shown in Fig. 4(d).
To summarize, our NMR measurements and statistical analysis, backed up by our QMC simulations, present a quite unambiguous demonstration of 3D critical physics in a material with a strongly quasi-2D structure and magnetic properties. With reference to the prediction that unconventional magnetic ordering bearing the hallmarks of BKT physics should appear uniquely in 3D coupled quasi-2D systems [21], our results [Figs. 3(c) and 4(d)] demonstrate that the field-induced onset and evolution of the ordered phase are quite conventional. From this we conclude that the inter-bilayer coupling in Ba 0.9 Sr 0.1 CuSi 2 O 6 is still significantly too large to access BKT physics, despite being only 1/100 of the magnon bandwidth (1/25 of the intra-bilayer interaction, a value that by any other measure appears worthy of the designation "quasi-2D"). Nevertheless, the restoration of observable 3D quantum criticality from an apparently minor chemical substitution constitutes a somewhat remarkable consequence of structural "disorder."

VI. DISCUSSION
We have found that the interaction parameters of . It is notable that these interaction parameters are not changed, within the error bars of our zero-and highfield experiments, by applied fields in excess of 20 T. More noteworthy still is the hierarchy of energy scales contained within the minimal spin Hamiltonian, which appear to define a quasi-2D material. The zero-field spectrum measured in Figs. 1 and 2 probes largely the intrabilayer energy scale, which gives the triplons a bandwidth of 2 meV. However, it is clear in Fig. 5(f) that the interbilayer energy of fractions of a Kelvin is the scale on which the 3D physics emerges.
We remark again that the FM effective intra-bilayer interaction removes any frustration of the inter-bilayer interactions, reinforcing the expectation of a simple scaling form around the QCP. In a windowing analysis of the critical scaling in BaCuSi 2 O 6 , an abrupt change in behavior towards an exponent φ = 1 was reported below 0.75 K, which was interpreted as a frustration-induced dimensional reduction [45]. Subsequent investigation has revealed the absence of frustration [53], the presence of structurally inequivalent bilayers [49,51], and the problem of the strong negative correlation between φ and H c1 [shown clearly in our Fig. 5(d)]. The experimental confirmation of the absence of frustration was accompanied by detailed QMC simulations of the quantum critical regime, from which it was possible to conclude that the true behavior of the system includes an intermediate fieldtemperature regime of non-universal scaling governed by an anomalous exponent whose value turned out to be φ > ∼ 0.72 over much of the effective scaling range [54]. This result raises the question of whether our present analysis, which obviously relies on a limited number of experimental data points, can be regarded as sufficiently precise to claim that Ba 0.9 Sr 0.1 CuSi 2 O 6 shows true 3D scaling. For this we appeal first to the fact we demonstrated by INS, that the Sr-substituted material has only one bilayer type and one species of triplon ( Figs. 1 and 2), whereas all the complications leading to anomalous scaling in BaCuSi 2 O 6 result from a new effective energy scale dictated by the separation of the inequivalent triplon modes. Thus in the Sr-substituted case one expects both a wider critical regime and a lack of alternatives other than the universal 2D or 3D forms, and our results are certainly not consistent with 2D scaling (φ = 1). Second, we note from the error bars in the QMC data for the layered model of BaCuSi 2 O 6 [54] that the anomalous effective exponent is in no way compatible with 2/3; by contrast, our NMR data for Ba 0.9 Sr 0.1 CuSi 2 O 6 favor strongly an exponent of 2/3 in Fig. 5(f) rather than any anomalous value, as do our QMC simulations for the model with only one type of bilayer.
The crux of our analysis is the materials-synthesis result that Sr substitution improves the regularity of the BaCuSi 2 O 6 structure, maintaining tetragonal symmetry and eliminating the orthorhomic transition that also pro-duces the structurally inequivalent dimer bilayers. However, the introduction of the smaller Sr ion in the interbilayer position of the Ba ions also creates a local strain most likely to act along theĉ-axis within the bilayers. The most important statement concerning any magnetic disorder produced by this structural disorder is that it can be at most a bond disorder, and not a site disorder, because the Cu 2+ sites are unaffected. Because of the inter-bilayer location of the Ba 2+ and Sr 2+ ions, this bond disorder is most likely to affect the local value of J , a parameter so small that our results are not very sensitive to the possibility of a distribution of values. The most important observation in the context of possible magnetic disorder is the result of our INS measurements that any resulting triplon broadening (lifetime) effects are too small to be distinguished from extrinsic factors (Sec. IIIB). From this we may conclude that local magnetic effects arising from this disorder are minimal, and that it has no discernible consequences for the magnetic Hamiltonian. Thus the disorder can be called "ideal" in the sense that only its qualitative and quantitative effects on the average structure are relevant. Regarding its density, while the nominal substitution was x = 0.1 and a compatible value of x = 0.08(2) was determined by energy-dispersive x-ray spectroscopy (EDX), a somewhat lower value of x = 0.05(1) was obtained from Rietveld refinement of synchrotron diffraction data for crushed crystals [55].
Further disorder effects within our experiment are the following. The density of magnetic impurities was estimated by fitting the Curie tail in the magnetic susceptibility, which arises at the lowest temperatures due to free spin-1/2 levels, and was found to be 3% in the best crystals selected for our experiments [55,57]. A typical nonmagnetic impurity in BaCuSi 2 O 6 is Cu 2 O, which can be created by the decay of CuO during crystal growth; it is visible on surfaces and in thinner crystals, allowing its presence to be minimized by careful selection of crystals. In our INS experiments, the floating-zone-grown single crystals have a reorientation in growth direction of up to 15 • as a consequence of the growth process [57]. However, it is important to stress that none of these possible sources and consequences of disorder are significant when compared to the complex and incommensurate crystal structure and magnetic Hamiltonian realized in the parent material [46,49,51,52,54].

VII. CONCLUSION
In summary, we have performed a comprehensive experimental and numerical analysis of the quasi-2D dimerized antiferromagnet Ba 0.9 Sr 0.1 CuSi 2 O 6 . The nonmagnetic site disorder due to this weak Sr-substitution of Han Purple affects only the lattice, acting to stabilize the tetragonal room-temperature structure of the parent material down to the lowest temperatures. We have measured the spin excitation spectrum at 1.5 K and zero ap-plied magnetic field to verify that indeed there is only one triplon branch, and hence that all dimer bilayers are structurally and magnetically equivalent in this compound. This triplon has a spin gap of 3.00 meV and is described by only three Heisenberg superexchange interactions, all on very different energy scales.
To examine the consequences of these interactions for field-induced magnetism and universal critical properties at the magnetic transition, we have gathered extensive thermodynamic data at high fields and low temperatures using a range of experimental techniques, including a detailed NMR study of the quantum critical phase boundary. For an optimal analysis of these latter data, we have performed Bayesian inference to obtain the critical exponent φ = 0.67 +0.05 −0.07 , determined by assuming a quantum critical regime of approximate width 0.4 K. QMC simulations of a 3D model with a single type of bilayer yielding the triplon band determined from the zero-field INS measurements provide excellent agreement with the quantum critical phase boundary, return a critical exponent φ = 0.66 +0.01 −0.01 , and indicate both a convergence to φ = 2/3 and the actual width of the regime over which this quantum critical scaling is obeyed. We conclude that the structural alteration caused by the small chemical modification of Han Purple to its 10% Sr-substituted variant fully restores the originally envisaged properties of 3D quantum critical scaling at the field-induced magnetic ordering transition of a structurally and magnetically quasi-2D system. Centre (SDSC) through project BISTOM C17-12, and of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through Grants No. SFB/TR49, UH90/13-1, and UH90/14-1.

Appendix A: Neutron spectrum calculations
The single-particle excitations in an AF spin-dimer material with a global singlet ground state are coherently propagating triplets, commonly referred to as triplons. To determine the interaction parameters of Ba 0.9 Sr 0.1 CuSi 2 O 6 , we have modelled the magnetic dynamical structure factor, S( Q, ω), by applying a triplonbased RPA approach [54,85]. In this framework, the zero-temperature dispersion of the single triplon mode is given by The INS intensities, I( Q, ω), are given by where A is an overall prefactor. The magnetic form factor, F mag ( Q), is known to be spatially anisotropic [54] and hence was modelled using the minimal parametrization Theoretically, the dynamical structure factor at zero temperature due to the triplon is given by where N is the number of dimers in the crystal and R the vector connecting the two ions of a single dimer. In experiment, N is absorbed in A and the δ-function is replaced by a Gaussian of unit integrated intensity and width Γ( Q), which includes both intrinsic dynamical and extrinsic instrumental resolution effects. In the modelled spectra shown in Figs. 1 and 2, Γ( Q) was fitted separately at every Q point. A and F 2 mag ( Q) were fitted using the intensity data shown in Figs. 2(d) and 2(e) while keeping the interaction parameters fixed to their optimal values determined using the dispersion data of Figs. 2(b) and 2(c).

Appendix B: Extraction of NMR phase boundary
To extract the phase boundary of Ba 0.9 Sr 0.1 CuSi 2 O 6 from the NMR spectra, we modelled the opening and increase of the peak-splitting visible in Figs. 4(a) and 4(b) by exploiting its direct proportionality to the order parameter [inset, Fig. 4(d)]. In the critical regime we therefore assumed that the splitting of peak frequencies is given by for each measurement temperature, T . This allowed us to perform simultaneous fits of all spectra at each T , and hence to optimize the accuracy of our estimate for the corresponding H c (T ).
To model the shapes of the NMR peaks [Figs. 4(a) and 4(b)] we used the Student t-distribution. A single peak was used below H c (T ) and two peaks above, where H c (T ) is the key parameter to be fitted. Implementing this behavior required a number of additional assumptions. The widths of the peaks below H c (T ) [blue in Figs. 4(a) and 4(b)] and of the lower-frequency peaks above H c (T ) (red) were treated as individual parameters, while the widths of the higher peaks above H c (T ) (green) were parametrized using a low-order polynomial in order to avoid spurious contributions from background features. The shape parameters of the t-distributions were modelled using a function linear in H − H c (T ), which allowed us to confirm that the shape remains almost constant for different fields. The relative intensity ratio of the two peaks, and the size of a constant background parameter, were taken to be the same for all spectra measured at the same temperature. To assist the fit convergence, and because the important parameter is the position or positions of the NMR peak(s), all integrated intensities were normalized to unity before processing.

Appendix C: Bayesian analysis of phase-boundary data
To analyze the NMR and QMC data close to the QCP within the framework of Eq. (1), it is necessary to specify the functions entering Eq. (2). Beginning with QMC, the critical fields, H c1 and {H c,j }, are input values for the simulation and hence are known exactly ({δH c,j } = 0). The corresponding critical temperatures, {T c,j }, are estimates based on extrapolation from repeated simulations with increasing grid sizes and we assumed them to obey a normal distribution around the true critical fields with standard deviations {δT c,j }. The likelihood function was therefore taken to be p(D QMC |θ QMC ) = jmax j=1 N (T c,j |α[H c,j − H c1 ] φ , δT c,j ), (C1) where N (x|µ, σ) denotes the probability density function of the normal distribution with mean µ and standard deviation σ, evaluated at x. In this expression θ QMC ≡ (α, φ) contains only two unknowns and j max , the index of the data point with temperature T c,j = T max , introduces the finite width of the quantum critical regime.
The statistical relationship between the phase-boundary parameters and the observations in this case is therefore that of a regression model.
In the phase-boundary data obtained by NMR, both the critical fields and the corresponding temperatures are estimated quantities. Each temperature, T c,j , is subject to both measurement uncertainties and possible systematic variations, and was treated as a noisy measurement of a single value. Each critical field, H c,j , was obtained from the NMR peak-splitting with an uncertainty that arises from background, noise, and the finite number of measurements, as a result of which it is essentially uncorrelated with δT c,j . The presence of uncertainty in two measurement dimensions means that the true point on the phase boundary at each j is unknown. Its location can be parametrized by either a temperature or a field coordinate, and here we chose the latter, introducing the parameters {H * c,j } to represent the true, unknown critical fields of which {H c,j } are uncertain estimates. Taking {H c,j } and {T c,j } to be independent and to have a normal distribution around these true phase-boundary points, with respective standard deviations {δH c,j } and {δT c,j }, the likelihood function was taken as ×N (H c,j |H * c,j , δH c,j ), (C2) where θ NMR denotes the extended parameter set (α, H c1 , φ, H * c,j ). In the Bayesian methodology, once the posterior is obtained its dependence on the introduced parameters {H * c,j } can be removed in an optimal manner by averaging over all possible values weighted by their posterior probability.
Turning to the prior probability distributions in Eq. (2), their selection should guide the inference process efficiently without introducing undue bias in the final result. α is known from studies of BaCuSi 2 O 6 [45,54] only to be of order unity (for measurements in T and K), and hence we assigned it an uninformative normal distribution on a logarithmic scale. For µ 0 H c1 , we used the gap ∆ = 3.00(1) meV obtained from the INS data and the g-factor g c = 2.32(2) [55] to obtain an initial estimate of 22.34 ± 0.24 T. Because φ is positive and is not expected to exceed unity, we used as its prior a half-normal distribution, N half (x|σ), with the scale parameter cho-sen to ensure that it is approximately flat in the region where φ is physically plausible, but does permit larger values. These considerations gave the prior probability distributions p(log 10 α) = N (log 10 α|0, 3), p(H c1 ) = N (µ 0 H c1 |22.34 T, 0.24 T), for the three model parameters. For the true critical fields, {H * c,j }, required to analyze the NMR phase boundary, these are unknown before measurement and thus we assigned a flat prior probability distribution, where U(x|a, b) denotes a uniform distribution with limits a and b. We made the minimal assumption that the range of possible values does not exceed 0-50 T, which includes the entirety of the ordered phase. Because all of these probabilities are independent, the joint prior probability distributions, p(θ QMC ) and p(θ NMR ), entering Eq. (2) are given by their product.
To investigate the unknown width of the quantum critical regime, we repeated the analysis of our NMR data multiple times with 339 mK ≤ T max ≤ 921 mK (5 ≤ j max ≤ 10 data points) and of our QMC data with 234 mK ≤ T max ≤ 939 mK (8 ≤ j max ≤ 30 data points). In each case we recovered the posterior by drawing random samples using emcee [86], a Markov-chain Monte Carlo implementation of an affine invariant ensemble sampler [87]. An ensemble of 200 chains was sampled for 40 000 steps to ensure equilibration of the Markov process before a further 10 000 steps were recorded, yielding a total of 2 000 000 parameter samples for each j max . These samples form point clouds in the parameter space of θ, whose density is proportional to the posterior probability, and one of these (obtained using the NMR data with j max = 6) is visualized using different projections in Figs. 5(a)-5(d). As a heuristic confirmation of convergence we note that the maximal autocorrelation time [87] was 117 steps. The posterior standard deviations are at least 20 times smaller than the prior ones in all cases, confirming that, after convergence, the results depend only on the data and not on our choice of weakly informative prior.