Proximate deconfined quantum critical point in SrCu₂(BO₃)₂

Single crystals of SrCu₂(BO₃)₂ were grown by the floating-zone method. We used a NiCrAl piston cell for the high-pressure NMR measurements and Daphne 7373 oil as the pressure medium. The pressure was calibrated using the low-temperature Cu₂O NQR frequency Reyes_1992 , and the highest pressure achieved was 2.4 GPa. NMR experiments were performed in two types of cryostat: a variable-temperature insert (VTI) was used for measurements in the temperature range from 1.5 K to 50 K and a dilution refrigerator was used to achieve temperatures ranging from 0.07 K to 2 K.

The ¹¹B spectra were obtained by the standard spin-echo technique. ¹¹B nuclei have spin $I=3/2$, whence the NMR spectra and line widths contain contributions from both local magnetism and electric-field gradients (EFGs). The local Hamiltonian for a nucleus with spin $I$ and quadrupole moment $Q$ can be described by the following form,

where the three terms describe, respectively, the coupling of the nuclear spin to the external field, the hyperfine coupling between the nuclear spin and neighboring electronic spins, and the coupling between the nuclear quadrupole moment and the local EFG. The factor $q$ in the last term is from the EFG tensor (the component along the principal axis), which is produced by the local electronic structure of the ions.

To identify the ordered phases and clarify their nature, NMR spectra were measured over a wide field range from $H=0.2$ T to $15$ T along the crystalline $c$ axis, with a top-tuning circuit allowing for the corresponding wide range of frequencies. The full spectrum was obtained by the frequency-sweep method, which covers one center line and two sets of satellites on each side of it. Figure 3a of the main text shows a typical NMR spectrum at 4 T. It contains one center line, with a frequency of approximately zero relative to ${}^{11}\gamma H$, where ${}^{11}\gamma$ is the Zeeman factor, and two sets of satellites aligned from $\pm 1.0$ MHz to $\pm 1.4$ MHz. We define the Knight shift as ${}^{11}K_{n}=f/{{}^{11}\gamma}H-1$, where $f$ is the position of the center peak in the spectrum.

The NMR magnetic field was oriented primarily along the $c$ axis, with a $8.6^{\circ}$ tilting applied to separate the four ¹¹B sites producing satellites in the spectrum. The angle is calibrated by different satellite frequencies $\nu_{Q}$ of four sets of satellites shown in Fig. 3a, as reported by the earlier NMR study Waki_JPSJ_2007 . The four satellites arise because each ¹¹B site has a different orientation of the principal EFG axis relative to the (tilted) external field. Each satellite line cannot be assigned to a specific ¹¹B site because we did not characterize the exact direction of the field in the sample $ab$ plane.

For a system with local inhomogeneity, the FWHM of the NMR center line has two additive terms, the hyperfine field contribution, which scales linearly with the applied field, and a second-order EFG correction, which is weak and scales inversely with the field. In a field of 0.2 T and at temperatures below 5 K, the FWHM of the center line is of order 2 kHz. By contrast, although the FWHM of the satellites has a similar hyperfine-field contribution, it has a large and nearly field-independent first-order EFG contribution, making the FWHM of each satellite larger than 10 kHz for the same temperature range and field.

The relaxation time ${}^{11}T_{1}$ was measured using the spin-inversion method by applying a $\pi$ inversion pulse (of duration 4 $\mu$s). In all phases displaying a line-splitting, the $1/^{11}T_{1}$ data reported in the main text were measured on the higher-frequency peaks, meaning $f_{R}$ in Fig. 4a, although $1/^{11}T_{1}$ measured on $f_{L}$ was verified to be consistent in every case with that on $f_{R}$.

Note that $1/^{11}T_{1}$ overall increases with the field, as shown Fig. 5 of the main text. In particular, this increase follows the expectations in the quantum-critical fan on the PS side of the transition, analyzed in Fig. 5e. This behavior indicates dominant contributions from the hyperfine field fluctuations originating from the critical spin fluctuations. By contrast, the EFG contributions, which are caused indirectly by bond and plaquette singlet fluctuations that couple to the lattice [causing fluctuations of $q$ in Eq. (S1)], are expected to decrease with the field because the PS fluctuations are suppressed by the field.

The most likely reason for the absence of visible EFG contributions in the $1/^{11}T_{1}$ data is that the singlet fluctuations in the frustrated magnet induce bond-length fluctuations of relatively small amplitude, i.e., the fluctuations of $q$ in Eq. (S1) are small and generate insignificant contributions to $1/T_{1}$ compared to the spin fluctuations mediated directly through the hyperfine coupling.

Evidence for different scales of the effective EFG and hyperfine couplings come from the way we have detected the PS and AFM order parameters using the NMR line shape. In the case of PS order, in Fig. 3b of the main text we used the FWHM of the central and satellite lines, because the expected peak splitting is too small to observe. Thus, the lattice deformation induced by the frozen singlets, i.e., the modulation of $q$ in Eq. (S1), is very small. In contrast, the peak splitting is very substantial in the AFM phase, Fig. 4. Given the small modulation of $q$ in the ordered PS state, the fluctuations of $q$ above and close to the PS transition temperature, induced indirectly by quantum fluctuations of bond and plaquette singlets, should also be small and contribute insignificantly to $1/T_{1}$ compared to the spin fluctuations mediated directly by the large coupling $A_{\rm hf}$. We are therefore justified in analyzing $1/T_{1}$ solely in terms of spin fluctuations detected through the hyperfine coupling—this assumption is often made in NMR studies of quantum magnets without the supporting evidence we have here from the weak response in the PS state.

It is also possible, in principle, that $1/^{11}T_{1}$, as an anisotropic probe of EFG and magnetic fluctuations, may not significantly sense the anisotropic EFG fluctuations geometrically by accident Yogi_JPSJ_2011 . Because of our direct evidence of insignificant EFG contributions due to the small amplitude of the lattice fluctuations, there is no need to invoke such a mechanism here, however.

To obtain the fitting parameters for the functional forms of $T_{\rm N}(H)$ and $T_{\rm P}(H)$, with these phase boundaries coming together at a low common $T_{\rm c}$ value at $H=H_{\rm c}$, we followed the Bayesian inference procedure of Ref. Allenspach_2021 . Sets of points spanning the multidimensional parameter space are generated by a Markov-chain Monte Carlo process and the output is a statistical model for the probability that a given parameter set can describe the measured data. The errors inherent to the data are then reflected systematically as the uncertainties in each of the fitting parameters. For further details we refer to Ref. Allenspach_2021 .

The functional forms used to model the transition temperatures are given in the main text and the maximum-probability results for pressures $2.1$ and $2.4$ GPa are shown in Fig. 5f. Probability distributions for the system parameters consistent with the measured data are shown in Sec. S6 of the SI, and these distributions underlie the error bars on the parameters reported in the main text.

The Hamiltonian of the CBJQM, i.e., the $J$-$Q$ model JQ_PRL_2007 defined on the checkerboard lattice Zhao_NP_2019 , subject to an applied magnetic field is given by

where $\mathbf{S}_{i}$ is an $S=1/2$ spin operator at site $i$ and $J>0$ is the nearest-neighbor AFM Heisenberg coupling, with the corresponding operator defined for convenience as a singlet projector,

The second term describes a partial ring-exchange-type interaction (or correlated singlet projection) defined on plaquettes $\square_{s}$ with a checkerboard arrangement on the square lattice. The third term describes the external magnetic field of strength $H$.

For the simulation results shown in Figs. 6b-6f of the main text, we took $J=1$ as the unit of energy to define the dimensionless parameter $g=J/Q$ and the reduced field $h=H/J$. We studied this model by stochastic series expansion (SSE) quantum Monte Carlo (QMC) simulations Sandvik_PRE_2002 on systems of sizes $L$$\times$$L$ up to a maximum of $L=64$ and down to a minimum temperature of $T=J/L$. The SSE method is exact, in the sense that it is based on a construction corresponding to a complete representation of the imaginary-time dimension of a quantum system when mapped into a classical model (on which the Monte Carlo sampling is performed), and thus the simulation results are only affected by statistical errors.

To characterize the PS–AFM transition (Sec. S8.1 and S8.2), we used an order parameter for the PS phase of the form

where $\mathbf{i}$ is the position vector of the checkerboard plaquettes and $i_{x}$ denotes the row index of $\mathbf{i}$. We also calculated the uniform magnetic susceptibility,

where $N=L^{2}$. In the simulations we obtained the susceptibility in the standard way from the magnetization fluctuations,

where $\langle m\rangle=0$ in the absence of external field. We extracted the spin excitation gap in the PS phase by fitting our numerical results to the form $\chi(T)\propto e^{-\Delta/T}$ at low temperatures for each field.

In addition to the above observables, which can be obtained using diagonal estimators in the $S^{z}$ basis used, we also study cross-correlations between the AFM and PS order parameters. In the case of $h=0$, we can use the diagonal correlator $\langle m_{z}^{2}m_{p}^{2}\rangle$, where $m_{z}$ is $z$-component of the staggered magnetization, but when $h>0$ the magnetization operator in $\langle m_{x}^{2}m_{p}^{2}\rangle=\langle m_{y}^{2}m_{p}^{2}\rangle$ is off-diagonal and a more complicated estimator based on the directed-loop SSE updates is required. Essentially, in the directed-loop updates used to evolve the configurations Sandvik_PRE_2002 , the probability distribution of the two open ends of a uncompleted loop (i.e., a string) at the construction stage is related to the off-diagonal operator $m_{x}^{2}+m_{y}^{2}$ (here in the case where the string ends are located at the same “time slice” in the SSE operator space). We follow the implementation discussed by Dorneich and Troyer Dorneich_PRE_2001 , with the modification that each contribution is also multiplied by the diagonal operator $m_{p}^{2}$ when accumulating the cross-correlation function $\langle(m_{x}^{2}+m_{y}^{2})m_{p}^{2}\rangle$.

To complement the simulations of the CBJQM, we also apply SSE QMC simulations to study an anisotropic XXZ Heisenberg $S=1/2$ spin model given by the Hamiltonian

where $\lambda$ expresses the spin anisotropy as a deviation from the Heisenberg model.

The XXZ model clearly has an exact O($3$) symmetry at $\lambda=0$, while perturbations with $\lambda<0$ imply an O($2$) and with $\lambda>0$ a Z₂ (Ising) order parameter. For $\lambda>0$, the system undergoes a first-order “spin-flop” transition versus the magnetic field; for small $h$ remaining in the Ising phase and transitioning through a level crossing into a canted XY (planar) AFM. We use this transition as a well understood benchmark case for comparing and contrasting with the CBJQM results in Sec. S8 of the SI.

In Fig. 4a and 4b of the main text we showed the splitting of the NMR center line as the system is driven from the PS to the AFM phase by increasing the applied field. The two field-induced peaks are located symmetrically around the zero of frequency measured relative to $f_{0}={{}^{11}\gamma}H$. This symmetrical line splitting indicates the onset of equal negative and positive hyperfine fields and constitutes direct evidence for AFM ordering. To complement the results for the field-driven PS–AFM transition, here we discuss the line splitting observed at constant field as a function of the temperature. The center line of the NMR spectrum at 2.1 GPa and 8.5 T is shown in Fig. S4a over a range of low temperatures. A single line is observed at and above 0.7 K, but below 0.65 K the spectrum shows two peaks, labeled as $f_{L}$ and $f_{R}$.

We again use the frequency difference $f_{R}-f_{L}$ as a proxy for the AFM order parameter, which we show as a function of temperature in Fig. S4b. We observe that $f_{R}-f_{L}$ rises rapidly on cooling below $T_{\rm N}$, and note that for this field, $1/^{11}T_{1}$ displays a clear peak at the same temperature, as shown in Fig. 5b of the main text.

Because no obvious three-peak signature of phase separation can be observed in the NMR spectra in Fig. S4a, our results in this case establish a continuous or very weakly first-order AFM transition on cooling in a constant field 8.5 T, which is significantly away from the common transition field $H_{\rm c}\approx 6.2$ T. Close to the first-order triple point that we have established in the main text at $(T_{\rm c},H_{\rm c})$ and $2.1$ GPa, the transition should remain first-order, while sufficiently far away from the triple point one should expect continuous transitions (as has been observed in model studies Sun_CPB_2021 ).

The curve fitted to the data points in Fig. S4b is of the form $f_{R}-f_{L}=F(T_{\rm N}-T)(T_{\rm N}-T)^{\beta}$, where $F(x)$ is a second-order polynomial in $x=T_{\rm N}-T$ to account for the cross-over from the asymptotic critical form to almost temperature independent at lower $T$. We do not have enough data for such a fit to produce quantitatively conclusive results, e.g., to rigorously test the exponent $\beta\approx 0.35$ corresponding to the expected 3D XY universality class for the $T>0$ XY AFM transition. Nevertheless the behavior is consistent with a continuous or a very weakly first-order transition, and the exponent $\beta\approx 0.30$ obtained from the fit is reasonably close to the expected value.

We do not have enough $f_{R}-f_{L}$ data at other fields to systematically study the change from first-order to continuous AFM transitions as the field strength is varied. Similar to the situation at 2.1 GPa, splitting of the center NMR line is also observed at 1.95 and 2.4 GPa in the same range of temperature.

Technically, because the ¹¹B sites are located slightly above the Cu plane and have dipolar coupling to the Cu moment Waki_JPSJ_2007 , establishing in-plane AFM moments on the Cu sites ought to produce a dipolar hyperfine field along the $c$ axis at the ¹¹B sites, which due to the AFM order is antiparallel on different sites. Our results are then consistent with the expected planar rather than $c$-axis AFM order, because the latter would produce a much larger uniform $c$-axis magnetization, which, as we discuss below in Sec. S4.4, is not observed at such low fields.

The low-temperature ¹¹B NMR spectra at 2.1 GPa and 2.4 GPa and at 0.07 K, presented in Fig. 4a and 4b of the main text, show that the field-induced AFM ordering is accompanied by a narrow regime of magnetic fields exhibiting phase coexistence around the AFM transition detected at our lowest temperature of $0.07$ K.

As shown in Fig. 4a, at $2.1$ GPa, a significant line splitting is detected all the way down to the field $H_{\rm c}$ determined from the fits illustrated in Fig. 5f, demonstrating that the transition is rather strongly first-order in this case. In contrast, at $2.4$ GPa, the line cannot be reliably separated into three different peaks below $5.8$ T, which is still some distance away from the transition field $H_{\rm c}=5.72$ T. Thus, while the transition is also first-order at this higher pressure, the discontinuity has weakened significantly relative to that at $2.1$ GPa, indicating that the system is approaching a continuous QCP (which we argue is a DQCP) located at only slightly higher pressure. The transition field also clearly moves down with increasing pressure, and the trends observed here suggest that $H_{\rm c}$ at the DQCP should be below $5.7$ T.

Because the hyperfine coupling between the ¹¹B nuclear spins and the Cu moments in SrCu₂(BO₃)₂ is dominated by dipolar interactions for the off-diagonal tensor elements of the hyperfine coupling Kodama_JPCM_2005 , the hyperfine field at each ¹¹B nucleus can be estimated in the AFM state for different ordered configurations of the local moments. In general, the hyperfine field, ${\bf H_{\rm in}}$ = ($H_{\rm in}^{x}$, $H_{\rm in}^{y}$, $H_{\rm in}^{z}$), on a ¹¹B nucleus due to one neighboring Cu spin can be written as

where ${\bf m}=(m^{x},m^{y},m^{z})$ is the local moment of the Cu ion and the hyperfine coupling tensor has the form

All of the tensor elements can be estimated by using the general form of the dipolar field, with the relative positions of the ¹¹B nucleus and the Cu ion as input.

In the AFM phase induced by a field applied along the $c$ axis, the Cu moments, which have strong Heisenberg interactions, should orient in the crystalline $ab$ plane. The AFM order is collinear type, with parallel moments on each Cu dimer being antiparallel to those on all four neighboring dimers (Fig. 1c of the main text). If one assumes that the Cu moments are oriented along the $[100]$ direction, the hyperfine fields at the four ¹¹B nuclei shown in Fig. S1 can be calculated by summing the contributions from all of the neighboring Cu sites, giving the result shown in Table S1(a). We note that the real hyperfine field in a material is usually enhanced by various factors beyond the simple dipolar-field calculation, but this approximation is sufficient for the qualitative conclusions we will draw. If one assumes that the Cu moments are ordered along the $[010]$ or $[110]$ directions, the resulting hyperfine fields at the four ¹¹B sites are those shown respectively in Tables S1(b) and S1(c).

These calculations show that $H_{\rm in}^{z}$ has one pair of negative and one pair of positive values for the four B sites when the AFM moments are oriented along the $[100]$ or $[010]$ directions. Such a pattern of hyperfine field components create a double NMR line splitting, which is consistent with our experimental observations (Fig. S4). By contrast, moments aligned in the [110] direction will split the NMR line into three, which is not consistent with our observations. Similarly, moments oriented in other directions, which are necessarily of lower symmetry, will also create three or more split NMR lines. Thus we conclude that the ordered moment in the field-induced AFM phase is aligned in the $[100]$ or the $[010]$ direction. Note that the small non-zero angle of the sample alignment in the field was not taken account in the calculations.

Because the hyperfine field components $H_{\rm in}^{x,y}$ are not zero for planar AFM moment orientations (Table S1), the spin-lattice relaxation rate $1/^{11}T_{1}$, measured with $H\parallel c$, should pick up transverse fluctuations from the AFM phase. This is demonstrated by the peaks appearing in $1/^{11}T_{1}$ at the transition temperature, $T_{\rm N}$, in Figs. 5b and 5d of the main text.

To further explore the nature of the AFM phase, we calculated the uniform magnetization, $M(H)$, at different applied fields from the NMR Knight shift. ${}^{11}K_{n}$ for the AFM phase is calculated from the average frequency, $f$, of the center line of each spectrum and the uniform magnetization is obtained from its magnetic part, $K_{s}$, using the expression $M={{}^{11}K}_{s}H/^{11}A_{\rm hf}$. For a pressure of 2.1 GPa, ${}^{11}K_{n}$ at our base temperature of 0.07 K is shown as a function of field in Fig. S5a. The resulting $M(H)$ is shown normalized to its saturation value, $M_{s}$, in Fig. S5b, for which we used a $g$-factor of $g_{c}=2.28$ kageyama_JPSJ_1998 .

With ${\vec{H}}\parallel{\hat{c}}$, we find that the uniform magnetization never exceeds $2\%$ of its saturation value at $15$ T. Such a high field is about $9$ T above $H_{\rm c}$, and therefore offers an energy scale of $14$ K to partially polarize the system. This energy scale is not very small compared to the microscopic AFM interactions ($J^{\prime}\approx 57$ K and $J\approx 39$ K at $2.1$ GPa Guo_PRL_2020 ), and because the couplings are frustrated the effective interaction scale may be even smaller. Experimentally, however, the $M$ values are still very low as shown above, indicating that the effective AFM interactions still lock the spins strongly to the XY plane with very little canting

An important aspect of the magnetization is the absence of significant (not clearly detectable) discontinuity at the transition field; $H_{\rm c}\approx 6.2$ T at $2.1$ GPa in Fig. S5b. We do expect some discontinuity, similar to spin-flop transitions in uniaxially anisotropic quantum magnets, but, as we discuss further in Sec. S8, an anomalously small discontinuity observed here at the PS–AFM transition is likely a consequence of the emergent O($3$) symmetry associated with the nearby DQCP.

The small induced magnetization also makes it very unlikely that a magnetization plateau could exist up to the highest field strength reached here. Thus, our data do not support a supersolid phase at these fields and pressures. Different studies suggest that a much larger uniform magnetization is required to stabilize such a phase Haravifard_NC_2016 ; Shi_arxiv_2021 .

It is tempting on the basis of Fig. 5f of the main text to fit our data for $T_{\rm P,N}$ directly to the form expected at a QCP, i.e., $T_{\rm P,N}\propto|H-H_{\rm c}|^{\nu z}$. At a DQCP, the dynamical exponent is $z=1$ and, for the most well studied case of the transition from an O($3$) AFM to a four-fold degenerate dimerized state, the emergent symmetry is SO($5$) and the correlation-length exponent $\nu$ has been estimated by QMC simulations to $\nu\simeq 0.46$ Zhao_PRL_2020 ; Nahum_PRX_2015 ; Sandvik_CPL_2020 (which applies to both the PS and the AFM order parameter). At a PS–AFM DQCP in zero field the emergent symmetry should instead be O($4$) on account of the PS order parameter being a scalar instead of the two-component dimer order parameter. In this case a similar value of $\nu$ as above was obtained Qin_PRX_2017 .

In the presence of an external magnetic field the symmetry is further reduced to O($3$) at a putative DQCP separating phases with scalar PS and O($2$) AFM order parameters. For this DQCP, $\nu$ has not been determined, but, in analogy with the slowly evolving exponents of the conventional O($N$) transitions (where there is only one order parameter), one can expect that $\nu$ would remain close to its SO($5$) and O($4$) values. It is also very possible that the O($3$) DQCP does not strictly exist but is a weakly first-order triple point, at which scaling properties may still be governed by the O($4$) DQCP as long as the external field is not very large, i.e., the zero-field system is close to the DQCP. This closeness to O($4$) symmetry is supported by the very small field-induced magnetization at $H_{c}$ in Fig. S5, i.e., the AFM order parameter is still O($3$) here for all practical purposes and can combine with the scalar PS order to form an effectively O($4$) emergent symmetry.

An argument against the above identical treatment of the PS and AFM transition temperatures is that $T_{\rm N}=0$ in an ideal 2D system, i.e., the system orders only exactly at $T=0$ on account of the continuous order parameter symmetry (which precludes $T>0$ ordering by the Mermin-Wagner theorem). In our case of O($3$) AFM symmetry reduced to O($2$) by a magnetic field there is a Kosterlitz-Thouless (KT) transition into a critical phase existing down to $T=0^{+}$, which we do not consider further.

In SrCu₂(BO₃)₂, the weak 3D couplings allow for order also at $T>0$, as we have found experimentally with $T_{\rm N}$ overall about half of $T_{\rm P}$ at comparable distances from $H_{\rm c}$ (Fig. 2d in the main paper). The 3D couplings also imply $T_{\rm c}>0$, and, therefore, the standard critical forms of $T_{\rm P}$ and $T_{\rm N}$ have to be modified. Phenomenologically, the form

where we have defined $\phi=z\nu$, can be used.

We first fit $T_{\rm P}$ and $T_{\rm N}$ independently, each with four parameters $T_{\rm c}$, $H_{\rm c}$, $a$ and $\phi$, following the Bayesian inference procedure Allenspach_2021 briefly described in the Methods section, with data errors included. Data in about 1 T range from the critical fields are used to obtain reasonable power-law fitting. The quality of the fits are depicted in Fig. S10, with $H_{\rm c}$ and $T_{\rm c}$ as below. With statistics within the 95% credible interval (equivalent to 2$\sigma$), at 2.1 GPa, the fit to $T_{\rm P}$ gives $H_{\rm c}=6.189\pm 0.017$ T, $T_{\rm c}=0.074\pm 0.026$ K, $\phi=0.542\pm 0.053$, and the fit to $T_{\rm N}$ gives $H_{\rm c}=6.185\pm 0.017$ T, $T_{\rm c}=0.070\pm 0.025$ K, $\phi=0.609\pm 0.107$; at 2.4 GPa, the fit to $T_{\rm P}$ gives $H_{\rm c}=5.731\pm 0.016$ T, $T_{\rm c}=0.065\pm 0.025$ K, $\phi=0.533\pm 0.060$, and the fit to $T_{\rm N}$ gives $H_{\rm c}=5.700\pm 0.017$ T, $T_{\rm c}=0.067\pm 0.024$ K, $\phi=0.522\pm 0.103$.

As shown above, the values of $H_{\rm c}$, $T_{\rm c}$, and $\phi$ obtained from $T_{\rm P}$ and $T_{\rm N}$ at each pressure are consistent within the errors. Therefore, it is reasonable to assume the same value of $\nu$ for both order parameters, which again is motivated in light of the expected duality of the PS and AFM phases Wang_PRX_2017 ; Qin_PRX_2017 . Considering that the two transitions meet at a bi-critical point $(H_{\rm c},T_{\rm c})$, a general fitting form

can be applied. A poor fit to the above form would indicate a more generic first-order transition.

To use Eq. (S12), it is important to establish the width of the critical regime, i.e., where $H$ is sufficiently close to $H_{\rm c}$ for the scaling form to apply, while, in the present case where $T_{\rm c}>0$, the transition temperatures also exceed $T_{\rm c}$. To obtain the five fitting parameters in Eq. (S12) describing the field-dependent transition temperatures measured at each of the two pressure values, the Bayesian inference procedure Allenspach_2021 is again applied.

The output of the process of sampling the parameters is a multidimensional probability distribution, and we first visualize this in Figs. S11a-S11d by showing projections of the distribution that illustrate its dependence on $\phi$, $H_{\rm c}$, and $T_{\rm c}$, for the two measurement pressures. It is clear that the distributions specify a single well defined maximum in the space of each parameter, with narrow intervals of uncertainty in the parameters that reflect the error bars inherent to the data.

Before discussing the optimal parameter values determined by the fitting procedure, it is necessary to establish the unknown width of the critical scaling regime. For this we performed the analysis using $9\leq N\leq 25$ data points, ordered in magnetic field by their separation from $H_{\rm c}$, and in Fig. S12 we show the evolution of the optimal fitting parameters $\phi$, $H_{\rm c}$ and $T_{\rm c}$ with $N$. One may anticipate a deteriorating fit at large $N$, where the points no longer obey critical scaling, and at small $N$, where there are simply too few points to extract a reliable functional dependence. Indeed we observe that as $N$ is decreased below 20, all the fitting parameters for both pressures converge towards values that are nearly constant over the range $12\leq N\leq 16$. However, for $N<11$ the fitting becomes less reliable and we neglect these estimates. Thus we take the parameters to follow for $N=12$, where the data at both pressures show good convergence to a reliable fit.

Another means of demonstrating the quality of the near-critical fit and the width in field of the critical scaling regime is to show the dependence of $|T_{\rm P,N}-T_{\rm c}|$ on $|H-H_{\rm c}|$ on logarithmic axis, as we do in Figs. S13a and S13b. In this form, data in the critical regime fall on straight lines whose gradient is $\phi$, and scaling of $T_{\rm P}$ to $T_{\rm N}$ is ensured by the ratio of the prefactors, $a_{\rm N}$ and $a_{\rm P}$, obtained from the Bayesian procedure. Only very close to $H_{c}$ do some points deviate from the expected fit, which nevertheless lies well within their error bars (expanded by the logarithmic axes), because of the difficulty we have in identifying $T_{\rm P,N}$ in this regime. Far from $H_{c}$, we find that many more than $N=12$ data points are well described by the critical scaling form with a single exponent $\phi$, and the fact that $T_{\rm P}$ and the $T_{\rm N}$ continue to fall on the same crossover curve out of the scaling regime underlines the duality between the PS and AFM ordered phases. From Figs. S13a and S13b we conclude not only that the critical fit is fully consistent but also that the width of the critical regime extends well beyond 1 T in both directions ($N=17$ at 2.1 GPa and $N=20$ at 2.4 GPa).

Returning now to the projected probability distributions as functions shown in Figs. S11a-S11d, the optimized fitting parameters, with uncertainties determined from the the 2$\sigma$ level, are $H_{\rm c}=6.184\pm 0.010$ T, $T_{\rm c}=0.066\pm 0.020$ K, $\phi=0.569\pm 0.040$, and $a_{\rm N}/a_{\rm P}=0.378\pm 0.034$ for $P=2.1$ GPa and $H_{\rm c}=5.719\pm 0.012$ T, $T_{\rm c}=0.062\pm 0.023$ K, $\phi=0.498\pm 0.070$ and $a_{\rm N}/a_{\rm P}=0.367\pm 0.030$ for $P=2.4$ GPa. The errors presented in the main text are at the $1\sigma$ level ($68\%$ credible interval). The most important single property of these fits is the remarkably low value of $T_{\rm c}$ at both pressures. Critical values of order 0.06 K are more than one order of magnitude below the typical (non-critical) values of both $T_{\rm P}$ and $T_{\rm N}$, which is not a feature of a generic bicritical or triple point and strongly suggests that the finite $T_{\rm c}$ is a weak, or residual, effect, for example one arising from 3D perturbations to a system controlled by 2D physics.

Turning to the critical exponent, $\phi$, we observe at the $2\sigma$ level that our $2.1$ GPa result, $\phi=0.569\pm 0.040$, and our $2.4$ GPa result (which has a somewhat larger error bar), $\phi=0.498\pm 0.062$, are mutually compatible. Certainly the exponent at 2.4 GPa, a pressure we have argued appears close to the DQCP of SrCu₂(BO₃)₂, is compatible with the estimates obtained for the SO($5$) DQCP, $\phi=z\nu$ with generically $z=1$ and $\nu\approx 0.46$ Sandvik_CPL_2020 ; Nahum_PRX_2015 . A very similar value has been obtained for the O(4) case Qin_PRX_2017 , which is the more likely type of DQCP controlling the behavior here.

In the O(4) DQCP case, the model system studied in Ref. Qin_PRX_2017 has anisotropic, planar spin interactions, i.e., O($2$) AFM order parameter, and four-fold degenerate dimerized state. However, it has been argued Lee_PRX_2019 that the DQCP is the same as that obtaining when a scalar PS order parameter combines with an O($3$) AFM order parameter. It is not clear whether an O($3$) DQCP exists (to our knowledge this case has not yet been considered theoretically), but given the extremely small uniform magnetization at the PS-AFM transition in SrCu₂(BO₃)₂ (Sec. S4), it is possible that the relevant DQCP in this case is still the O($4$) one. In other words, the uniaxial deformation of the O($4$) symmetry is insignificant at the temperatures $T\gg T_{\rm c}$ for which we study critical scaling.

Given that at least the $2.1$ GPa system has a rather strongly discontinuous behavior of the AFM order parameter (Fig. 4c of the main text), the applicability of the near-critical form Eq. (S12) of the transition temperatures can be questioned. The very low $T_{\rm c}$ values compared to the transition temperatures when $H$ is not close to $H_{\rm c}$ is then unusual and calls for a mechanism not requiring very close proximity to a QCP. As we have discussed in the main text, emergent symmetry induced by a 2D DQCP far into a first-order line is such a mechanism.

In an ideal 2D system with an O($N$) order parameter and $N>2$, a perturbation making one component of the interactions larger pushes the system to the “Ising side,” where it has an excitation gap and undergoes a finite-temperature phase transition with a scalar order parameter. In contrast, a perturbation making one component smaller still has a continuous order parameter symmetry, O$(N-1)$, and there can be long-range order only at $T=0$ (with a quasi-ordered KT phase at $T>0$ in the special case $N=3$).

A renormalization-group analysis of the Heisenberg model ($N=3$) with this type of anisotropy Irkhin_PRB_1998 deduced a logarithmic form of the transition temperature on the Ising side; the same analysis should also apply to $N>3$. In the CBJQM at zero external field, a first-order PS–AFM transition was found with a coexistence state whose order-parameter distribution indicated the emergence of O($4$) symmetry up to the largest system sizes studied, $L\simeq 100$ Zhao_NP_2019 . Then, moving into the PS side corresponds to a uniaxial deformation of the order parameter and the logarithmic form is expected; specifically Irkhin_PRB_1998

where $g=J/Q$. This form was also confirmed in Ref. Zhao_NP_2019 for $g$ close to $g_{c}$. In the most likely scenario, the emergent symmetry is violated at some long length scale $\Lambda$, corresponding to some small energy scale $\epsilon$ (which can be taken as the gap of the Goldstone mode that corresponds to the Ising direction). The above form should then only be valid down to $T_{\rm P}$ of order $\epsilon$, but this energy scale was not reached in Ref. Zhao_NP_2019 .

In our experiments on SrCu₂(BO₃)₂, we study the phase transitions at fixed pressure versus the magnetic field, and the emergent symmetry is reduced to O($3$). Also in this case must there be a smooth approach of $T_{\rm P}$ to $0$ (assuming now momentarily that the symmetry is asymptotically exact) and Eq. (S13) should apply with $g_{c}-g$ replaced by $(H_{\rm c}-H)^{a}$, where the exponent $a$ accounts for the different forms of gap closing versus $g$ and $h$; $(g_{c}-g)^{1/2}$ Irkhin_PRB_1998 and $H_{\rm c}-H$, respectively. However, the value of the exponent can be absorbed into the factor of proportionality and we do not need to consider it further.

Experimentally, we expect violations of the O($3$) symmetry not only from the fundamental 2D effect related to the distance to the DQCP, but also because of 3D effects (and potentially other sources, like impurities, that we have not discussed). Thus, the logarithmic form should break down below some temperature. As in the case of the near-critical form Eq. (S12), we may again add $T_{\rm c}$ as an offset optimized in the fitting procedure. However, our fits to are actually optimal with $T_{\rm c}=0$, i.e., Eq. (S13) without any added constant, though statistically acceptable fits are also obtained with $T_{\rm c}$ values up to those obtained in the previous section. The fit shown in Fig. 5f of the main paper are for zero offset, and the further analysis presented below also was done with $T_{\rm c}=0$.

It is not immediately clear whether the AFM transition temperature can also be incorporated into a common fitting form in this case. Interestingly, however, weak inter-layer (3D) couplings $J_{\perp}$ between 2D Heisenberg AFM layers also leads to a logarithmic form; $T_{\rm N}(J_{\perp})=\ln^{-1}(C/J_{\perp})$ Irkhin_PRB_1998 . In the experiments the microscopic coupling $J_{\perp}$ is fixed, but it is still possible that an effective, “renormalized” 3D coupling $J_{\perp}(H)$ can be defined, in light of the fact that 2D AFM order is present only above $H_{\rm c}$ and increases with $H$. Thus, we assume here that both transition temperatures take the common form

though we are less confident in the form of $T_{\rm N}$ than $T_{\rm P}$ and only show results for the latter in Fig. 5f of the main paper. Here we will discuss both transitions.

At both pressures, $2.1$ and $2.4$ GPa, the optimal $H_{\rm c}$ value is statistically indistinguishable from the value obtained using the near-critical form of Eq. (S12). To gauge the quality of the logarithmic fit, in Figs. S13c and S13d we show the inverse transition temperatures on semi-log axes. The logarithmic fit remains valid over a significantly larger range of fields at $2.4$ GPa (where it includes 15 data points) than at $2.1$ GPa, which would indeed be expected if $2.4$ GPa lies closer to the DQCP, because the emergent symmetry would then be better established (i.e., manifested on longer length scales).

Although the near-critical forms in Figs. S13a,b remain valid somewhat further away from $H_{\rm c}$, corrections are expected here in both cases, and the non-universal width of the scaling regime is not a criterion for selecting the best critical form. As mentioned in the main text, each of the fitted forms has its range of validity in principle, but the experiments may be in a cross-over region where there are corrections to both forms but they still work reasonably well.