Field-driven phase transitions in a quasi-two-dimensional quantum antiferromagnet

We report magnetic susceptibility, specific heat, and neutron scattering measurements as a function of applied magnetic field and temperature to characterize the $S=1/2$ quasi-two-dimensional frustrated magnet piperazinium hexachlorodicuprate (PHCC). The experiments reveal four distinct phases. At low temperatures and fields the material forms a quantum paramagnet with a 1 meV singlet triplet gap and a magnon bandwidth of 1.7 meV. The singlet state involves multiple spin pairs some of which have negative ground state bond energies. Increasing the field at low temperatures induces three dimensional long range antiferromagnetic order at 7.5 Tesla through a continuous phase transition that can be described as magnon Bose-Einstein condensation. The phase transition to a fully polarized ferromagnetic state occurs at 37 Tesla. The ordered antiferromagnetic phase is surrounded by a renormalized classical regime. The crossover to this phase from the quantum paramagnet is marked by a distinct anomaly in the magnetic susceptibility which coincides with closure of the finite temperature singlet-triplet pseudo gap. The phase boundary between the quantum paramagnet and the Bose-Einstein condensate features a finite temperature minimum at $T=0.2$ K, which may be associated with coupling to nuclear spin or lattice degrees of freedom close to quantum criticality.


Introduction
Molecular magnets form intricate quantum many body systems that are amenable to experimental inquiry. Their low energy degrees of freedom are quantum spins associated with transition metal ions such as copper (spin-1/2) or nickel (spin-1). The interactions between these spins can be described by a Heisenberg Hamiltonian of the form Experiments on these materials can test existing theoretical predictions and also reveal entirely new cooperative phenomena [1,2,3,4,5].
If the generic ground state of transition metal oxides is the Néel antiferromagnet (AFM), it is singlet ground state magnetism for molecular magnets. When molecular units have an even number of spin-1/2 degrees of freedom they typically form a molecular singlet ground state. If inter-molecular interactions are sufficiently weak the result can be a quantum paramagnet. This is the state that Landau famously envisioned as a likely outcome of antiferromagnetic exchange interactions. [6] There are also indications that geometrical frustration on lattices containing multiple triangular plaquettes can suppress Néel order and/or produce singlet ground state magnetism as is seen for example in the 3D system Gd 3 Ga 5 O 12 (GGG) [7,8] and possibly in CuHpCl. [9] While magnets with an isolated singlet ground state generally have no zero field phase transitions, there are at least two field induced phase transitions that terminate in zero temperature quantum critical points. The lowest lying excitations in such systems typically have spin S = 1, and thus the lower field critical point, H c1 , occurs as the S z = 1 triplet branch is field driven to degeneracy with the singlet groundstate at a certain critical wave vector. The system becomes magnetized even in the zero temperature limit for fields above H c1 . As the field is increased further, the system is ultimately ferromagnetically polarized at a second QCP characterized by an upper critical field, H c2 where the field energy overwhelms the spin-spin interactions. Above this field, a second spin-gap should open and increase linearly with applied field above H c2 .
What happens as a function of field and temperature between these critical points depends critically on the inter-molecular connectivity. Three-dimensional systems such as Cs 3 Cr 2 Br 9 [10] and TlCuCl 3 [11,12,13,14] exhibit a phase transition to long range antiferromagnetic order both as a function of applied field and applied pressure. For quasi-one-dimensional connectivity the high field phase should be an extended quantum critical phase. However, systems that have been examined thus far, such as Cu(NO 3 ) 2 · 2.5H 2 O [15,16], and Cu(L − aspartato)(H 2 O) 2 [17] exhibit long range Néel order indicating that they may not be as one-dimensional as originally envisioned. Another complication can be multiple molecules per unit cell which lead to staggered fields that break translation symmetry elements and preclude a critical phase. Several new spin ladder systems have recently been identified that may enable access to a critical phase above H c1 . [18,19] Two-dimensional (2D) systems are of particular interest since they are at the lower critical dimension for Néel as the parent compounds to high temperature superconductors. BaCuSi 2 O 6 for example features spin-1/2 pairs coupled on a two dimensional square lattice. [20] An anomaly in the critical exponents and in the phase boundary at low temperatures has been interpreted as the result of a dimensional crossover from three to two dimensions close to quantum criticality. [21] In quasitwo dimensional SrCu 2 (BO 3 ) 2 competing interactions play a central role. Its unique geometry frustrates inter-dimer interactions, yielding effectively decoupled dimers [22]. While the singlet ground state of the material is known exactly, the three finitemagnetization plateaus and the critical phases that surround them are more challenging to understand. The purely organic quantum AFM F 2 PNNO forms a quasi-2D network of spin dimers with a singlet ground state and a 1 meV gap to magnetic excitations. [23] The phase diagram includes a magnetization plateau that extends from µ 0 H = 15 T to µ 0 H = 25 T as well as regions with continuously varying magnetization that may indicate unique 2D quantum critical phases. Neutron scattering is required to disentangle these complex and unfamiliar high field phases. Unfortunately few of these phases can be accessed with current high field neutron scattering facilities, which are limited to fields below 17.5 Tesla.
The material examined here, piperazinium hexachlorodicuprate (PHCC), is a highly two dimensional quantum paramagnet with a lower critical field of 7.5 T that is well within reach of current high field neutron scattering and an upper critical field of 37 T that can be reached in pulsed field magnetization measurements. The large difference between the critical fields indicates strong inter-dimer interactions and a potential for novel, strongly correlated two-dimensional physics. PHCC, (Cu 4 H 12 N 2 )Cu 2 Cl 6 , consists of Cu 2+ spin- 1 2 sites interacting predominantly within the crystallographic ac-planes [24,25], as illustrated in figure 1. Prior zero-field inelastic neutron scattering measurements probed the spin dynamics of the low-temperature quantum paramagnetic (QP) phase, and established the strongly 2D nature of magnetism in PHCC. The spin gap ∆ = 1.0 meV, and the magnetic excitations show dispersion with an 1.8 meV bandwidth in the (h0l) plane. There is less than 0.2 meV dispersion perpendicular to that plane, indicating a highly two dimensional spin system. The spin-spin interactions in the a-c-plane have the connectivity of a frustrated orthogonal bilayer, as shown in figure 1(c). Bonds 6 and 3 form two extended 2D lattices linked by Bond 1, which thus is J ⊥ in the bilayer description. Bonds 2 and 8 provide additional frustrated interlayer interactions.
In the current study, we employed inelastic and elastic neutron scattering as well as specific heat and magnetic susceptibility measurements to characterize the lowtemperature phases of PHCC. Unlike 3D coupled-dimer systems, such as TlCuCl 3 , we find a significant separation between the phase boundary to long range order (LRO) and the crossover line delineating the high-field boundary of the QP regime. This allows access to a disordered renormalized classical (RC-2D) regime near the lower QCP. Specific heat measurements in the disordered low-field QP phase agree well with Crystal structure of piperazinium hexachlorodicuprate (PHCC), (C 4 H 12 N 2 )(Cu 2 Cl 6 ). (a) View along the a axis, showing well-separated Cu-Cl planes. (b) View of a single Cu-Cl plane. Dotted lines between Cl sites show potential halidehalide contacts. (c) Second view of a single Cu-Cl plane, showing only the interacting Cu 2+ spin sites (solid circles). Significant antiferromagnetic spin-spin interactions are depicted as light blue (frustrated ) and black (satisfied ) lines with the same bond numbering used in [24]. Line thickness is proportional to the product of exchange constant and spin-spin correlation function, |J d S 0 ·S d |. Vectors represent the ordered spin structure determined by elastic neutron scattering at T = 1.65 K and µ 0 H = 13.7 T. Dashed lines are the a − c plane unit cell with room temperature lattice constants a = 7.984(4)Å, b = 7.054(4)Å, c = 6.104(3)Å, and α = 111.23 (8) • , β = 99.95(9) • , and γ = 81.26(7) • , space group P1, from [26].
series expansions from a quantum bilayer model. Combining elastic neutron scattering and field dependent specific heat results allows a determination of the critical exponents associated with the transition to LRO. We have previously reported a reentrant phase transition at low temperature near H c1 wherein the magnetic LRO can melt into the QP phase as the temperature is either increased or decreased [27]. We have proposed that coupling between other degrees of freedom is important near the QCP, and thus can have a non-negligible effect on the phase diagram in this regime. We also have recently examined the LRO transition and found good agreement with a representation of the phase boundary away from the reentrant phase transition in terms of magnon Bose-Einstein condensation. In this communication, we provide further experimental and theoretical information regarding these phenomena, including field dependent neutron scattering, magnetization, and specific heat data, and a more detailed theoretical discussion of magnon Bose-Einstein condensation in PHCC.

Experimental techniques
Single crystal samples of deuterated and hydrogenous PHCC were prepared via previously described methods [24].
Magnetic susceptibility measurements were performed at the National High Magnetic Field Laboratory pulsed field facility at the Los Alamos National Laboratory using a compensated coil susceptometer. The sample environment consisted of a pumped He 3 cryostat located in the bore of a 50 Tesla pulsed field resistive magnet. Measurements were performed with the magnetic field along the b-, a * -and c-axes on hydrogenous samples with masses between m = 0.7 mg and m = 1.4 mg. An empty coil measurement was performed for each individual pulsed field measurement of the sample as a direct measure of the sample-independent background. Specific heat measurements were performed as a function of magnetic field, µ 0 H ≤ 14 T, and temperature, T ≥ 1.8 K, on a deuterated single crystal sample, m = 6.76 mg, with H a * using a commercial calorimeter (Quantum Design PPMS).
Neutron scattering measurements were performed using the FLEX triple axis spectrometer at the cold neutron facility of the Hahn-Meitner Institut in Berlin, Germany. The sample was mounted on an aluminum holder and attached to the cold finger of a dilution refrigerator within a split coil µ 0 H = 14.2 T superconducting magnet. A magnetic field in excess of the superconducting critical field for aluminum was always present to ensure good thermal contact to the sample. The sample consisted of two deuterated single crystals with a total mass of m = 1.75 g coaligned within 0.5 • . The sample was oriented in the (h0l) scattering plane with H b. Measurements were performed using a vertically focusing pyrolytic graphite (PG(002)) monochromator followed by a 60 ′ horizontal collimator. Inelastic measurements used a liquid nitrogen cooled Be filter before the sample and a fixed final energy of E f = 3 meV. A horizontally focusing PG(002) analyzer was used for inelastic measurements. Elastic measurements employed a flat PG(002) analyzer with 60 ′ collimation before and after the analyzer. Measurements characterizing the temperature and magnetic field dependence of Bragg peaks used a liquid nitrogen cooled Be filter before the sample and neutron energy E i = E f = 2.5 meV. Measurements of magnetic and nuclear Bragg peaks for determining the ordered magnetic structure were performed using the same collimation with a room temperature PG filter before the sample and E i = E f = 14.7 meV with a vertically focusing PG(004) monochromator and a flat PG(002) analyzer.

Susceptibility and magnetization
Pulsed field measurements of the magnetic susceptibility χ(H) up to µ 0 H = 50 T with H b are shown in figure 2. For T ≤ 4 K, the susceptibility has two sharp features. The lower field feature marks closure of the spin gap while the upper feature indicates the field where M(H) saturates, and the system enters the fully polarized (FP) regime. M(H) at T = 0.46 K, obtained by integrating the corresponding χ(H) data, is shown as a solid line in figure 2. At low temperatures, these features in χ(H) indicate the location of the critical fields µ 0 H c1 ≈ 7.5 T and µ 0 H c2 ≈ 37 T. We note that although there is a local minimum in χ(H) at µ 0 H = 23.9 T near H m = 1 2 (H c1 + H c2 ), and hence a point of inflection in the magnetization, there is no indication of fractional magnetization plateaus for T ≥ 0.46 K.
As T increases, the two features in χ(H) broaden, and above T = 9 K, they merge into a single rounded maximum that moves to lower fields as the temperature is further increased. For T > 20 K, χ(H) decreases monotonically with H. Despite the thermal broadening above T = 5 K, several regimes may be identified from χ(H, T ). These are most clearly seen in figure 3, which shows χ(H, T ) at fixed field. The rounded peak and exponential decrease in χ(T ) for low magnetic fields is characteristic of a gapped antiferromagnet and is similar to prior results obtained at H = 0 for PHCC [24,28,26]. As the field increases, the peak in χ(T ) shifts to lower T , and the low-temperature activated behavior vanishes with the spin-gap. As shown for µ 0 H = 10 T in figure  3(a), for fields µ 0 H c1 < µ 0 H < 12 T the rounded maximum that signals the onset of 2D antiferromagnetic spin correlations leading to the gapped state is still visible (near  (a) µ 0 H ≤ 20 T, (b) 23 T ≤ µ 0 H < µ 0 H c2 , (c) H > H c2 . Data were obtained by combining all individual pulsed field measurements and averaging over the range µ 0 (∆H) = 0.5 T for each magnetic field shown. Lines connecting data points are guides to the eye. Curves in (c) for T < 8 K are fits described in the text that yield the field dependence of the gap above H c2 shown in the inset to (c). The solid line in the inset to (c) has the slope gµ B H. T = 12 K for µ 0 H = 10 T), but there is then a sharp rise in χ(T ) as the system enters the strongly magnetizable regime at lower T . As shown in figures 3(a) and 3(b) the boundary of this regime continues to move to higher T for µ 0 H < 25 T, before moving back to lower T as H approaches H c2 . In the FP region above H c2 activated behavior returns, as shown for µ 0 H = 40 T and µ 0 H = 47 T in figure 3(c). This indicates a second spin-gap regime. Note that as the field increases above H c2 , the onset of the activated susceptibility moves to higher T , indicating that this second spin gap increases with field. An overview of χ(H, T ) is given in the color contour plot in figure 11(a) and will be discussed in more detail below.
Low-temperature pulsed field susceptibility measurements were also performed at T = 0.46 K for H c and H a * . These are shown together with the H b data near H c1 and H c2 in figure 4. Although the overall shape of χ(H) is similar in all three orientations, the differences in the onset fields shown in figure 4 indicate that the critical fields are orientation-dependent. Because the ratio of the critical fields is

Specific heat
Measurements of the specific heat as a function of temperature and magnetic field are shown in figure 5. Prior measurements were performed on protonated samples in the same orientation, H a * , but were restricted to T ≤ 2 K and µ 0 H ≤ 9 T. The total specific heat C p is shown for representative fields in figure 5(a). In zero field, C p shows a broad rounded maximum at T ≈ 7 K, with an exponentially activated region at lower T associated with the formation of the spin-gap. Above H c1 the form of the high-temperature specific heat changes, and as shown for µ 0 H = 14 T, a lambdatype anomaly can be seen that is associated with the phase transition to magnetic LRO at T = T N . In figure 5(b) we plot the magnetic portion of the specific heat, C m = C p − C lattice , where C lattice ∝ T 3 is the lattice/phonon contribution to the specific heat determined by comparison of the data to the bilayer model discussed later in the text. This lattice contribution is shown as a dashed line in figure 5(a). Figure 5(b) shows the changes in the low-field activated behavior of C m as H increases and the spin-gap closes, as well as the increase in magnitude of the lambda peaks as H increases above H N . The area under the peak is greater for larger applied fields indicating that the T = 0 ordered moment increases with magnetic field. The field dependence of C p is plotted in figure 5(c) for selected temperatures. We observe a sharp peak at the transition to LRO at H = H N , which has the same field and temperature dependence as that seen in C p (T ) at fixed H.

Elastic neutron scattering
Neutron diffraction was used to explore the nature of the LRO transition indicated by heat capacity data. Field and temperature dependent antiferromagnetic Bragg peaks were observed at wave vector transfer Q = ( (2n+1) 2 , 0, (2m+1) 2 ) for integer n and m throughout the (h0l) plane, corresponding to a magnetic structure with a doubling of the unit cell along the a-and c-axes. Figure 6 depicts a wave-vector scan through the ( 1 2 01 2 ) Bragg peak as a function of (a) magnetic field and (b) temperature. The scattering intensity at T = 1.7 K increases as a function of applied field. At µ 0 H = 14.2 T (figure 6(b)) there is only a flat background at T = 10.1 K; however, for T = 3.7 K the broad maximum indicates short  range AFM correlations. At T = 1.7 K the resolution limited peak indicates long range antiferromagnetic order. Scans along the (h01 2 ) and (ζ0ζ) directions were also resolution limited. Magnetic Bragg peaks were found at locations that were indistinguishable from the nuclear half-wavelength, λ/2, Bragg peaks, which were observed upon removal of the Be filter to admit λ/2 neutrons from the PG(004) monochromator Bragg reflection. No indication of incommensurate order was found along the (h01 2 ) and ( 1 2 0l) directions through the entire Brillouin zone at T = 1.9 K and µ 0 H = 11.8 T.
To characterize the onset of antiferromagnetic order, the ( 1 2 01 2 ) Bragg peak was measured as a function of magnetic field at sixteen temperatures for 0.065 ≤ T ≤ 2.85 K. Data were obtained by continuously sweeping the magnetic field (µ 0 dH/dt = 0.01 or The intensity increases continuously upon cooling below a distinct critical temperature as expected for a second order phase transition. Figures 9(a) and 9(c) show low-temperature magnetic field sweeps of the ( 1 2 01 2 ) magnetic Bragg intensity plotted as a contour plot and as intensity versus field in the vicinity of the critical field for the transition to Néel order, which we denote by H N . The data indicate that the ordered phase is thermally reentrant for µ 0 H ≈ 7.5 Tesla. Upon heating from 0.065 K, H N initially decreases with increasing temperature. However, there is an apparent minimum in the critical field curve at T = 0.20 K, above which temperature H N increases on heating. The same data are plotted as a function of temperature at fixed fields in figure 9(b). These data show a peak in the scattering intensity at T = 0.20 K for fields between µ 0 H ≈ 7.45 T and µ 0 H ≈ 7.9 T. The peak location remains at T ≈ 0.2 K until vanishing for fields above µ 0 H ≈ 8.1 T.

Inelastic neutron scattering
Inelastic neutron scattering was used to probe spin-dynamics as a function of magnetic field and temperature. In particular, constant wave vector scans at Q = ( 1 2 01 2 ) were performed as a function of magnetic field at T = 0.06 K and T = 4.5 K. A portion of these data is shown in figure 10. The low-temperature scan at H = 0 T features a resolution limited excitation athω ≈ 1 meV, consistent with the previous zero field   measurements [24]. In an applied field this excitation splits into three components indicating a zero field triplet 1 meV above a singlet ground state. The intensity distribution is consistent with the upper and lower peaks resulting from transitions to S z = ±1 states while the central highest intensity peak corresponds to a transition to an S z = 0 state. The lower energy, S z = +1, peak is observable until µ 0 H ≈ 6 T, where it is obscured by elastic nuclear incoherent scattering. Both remaining gapped excitations increase in energy continuously as the field is increased above H c1 . The higher energy excitation increases in energy and decreases in intensity at a faster rate than the excitation which evolves out of the S z = 0 component of the H < H c1 triplet. There are also subtle differences in the temperature and field dependence of the width of the excitations, which we shall return to later.

Susceptibility and magnetization
A second view of the magnetic susceptibility data shown in figures 2 and 3 is provided by the color contour map in figure 11(a). Since the sharp features in χ(H) at the lowest-T herald the onset of finite magnetization at H c1 and the entry into the fully polarized (FP) regime at H c2 , we identify the continuation of these features for T ≤ 5 K as experimental signatures of the closing of the low-field spin gap and the opening of the high-field spin gap with increasing field in the vicinity of the critical fields H c1 and H c2 , respectively. These features provide a starting point in constructing an H-T phase diagram for PHCC, and are shown as open circles in figure 11(b). By extending these points through inclusion of the sharp features visible in χ(T ) for T > 5 K, they are seen to demarcate the entire large-susceptibility region in figure 11(a). At non-zero T , these points delineate a crossover, rather than a phase transition. Near the lower QCP  We also identify two additional crossovers in the H-T plane of PHCC that are defined by rounded maxima in χ(T ) for H < H c1 and H > H c2 . These are demarcated by open diamonds in figure 11(b). At low field, the resulting boundary indicates where antiferromagnetic correlations begin to develop as the system is cooled toward the spin gap regime. This process has been well-studied in other spin gap systems, such as Haldane spin-1 chains [29]. The intersection of this crossover with the zero-field axis in figure 11(b) marks the onset of a high-temperature semi-classical paramagnetic regime (P) where the system has not yet attained significant antiferromagnetic correlations and the magnetic susceptibility is that of a disordered ensemble of exchange-coupled spins with one energy scale, i.e. Curie-Weiss behavior.
The high field crossover regime (open diamonds above H = H c2 in figure 11(b)) is not as well characterized experimentally, but should mark the onset of appreciable ferromagnetic correlations as the system is cooled at high field toward the FP state. The excitations in the FP state are expected to be ferromagnetic spin waves. As for a conventional ferromagnet in a field, these excitations will show a field-dependent gap, although here the energy scale for the gap should be set by ∆ F P = gµ B (H − H c2 ). The evolution of the magnetic susceptibility as a function of magnetic field and temperature has been examined numerically for the two-leg spin-ladder [30]. Calculations depict a low-temperature field dependent magnetization, which is symmetric about H m = 1 2 (H c1 + H c2 ) with square root singularities near the upper and lower critical fields, consistent with results for other 1D spin-gap systems [31]. Its phase diagram is thus similar to our findings for PHCC. The calculations also report similar susceptibility anomalies: a double peaked function at low temperatures, a single rounded peak that moves to lower fields as the temperature is increased, and a monotonically decreasing field dependence as the temperature is increased beyond a certain temperature, T cl . A crossover temperature marking the boundary between the classical and quantum regimes, was identified as T cl = 1.109∆/k B for the two-leg spin ladder. For PHCC, this corresponds to T cl = 1.109∆/k B = 13.1 K which agrees well with the zero-field lower temperature bound of the crossover regime defined in figure 11 where we find T cl ≈ 15.4 K based upon susceptibility measurements.
Close to the critical fields the low-temperature field dependent magnetization should follow power laws of the form M ∝ |H − H c | 1/δ , where H c is a critical field [32]. In terms of magnetic susceptibility, this power law becomes immediately above (χ + ) and below (χ − ) the respective lower and upper critical fields. For PHCC Figures 2 and 4 indicate that δ 2 > δ 1 ≈ 1 to account for the observed curvature near the critical fields. Because δ 1 ≈ 1, we fit M(H) rather than χ(H) to determine δ. This is done simultaneously for all T = 0.46 K single crystal orientations with a common value of δ 1 over a range of one Tesla above H c1 . The single crystal low-temperature, T = 0.46 K data were also simultaneously fit over a range of one Tesla below H c2 using Eqn. (2b). The simultaneous fits shared a common δ parameter; whereas critical fields and overall amplitudes were allowed to vary between the different orientations. The resulting values, δ 1 = 0.91(1) and δ 2 = 1.8 (2), and their corresponding power-laws are shown in figures 2 and 4. For comparison, 1D singlet ground state systems are predicted to have symmetric susceptibility and magnetization curves with δ 1 = δ 2 = 2 [32], and the 2D square lattice Heisenberg model has an analogous upper critical field with predictions ranging from δ 2 = 1 to δ 2 = 1 with logarithmic corrections [33,34]. An interpretation of H c1 as magnon Bose-Einstein condensation yields the prediction of δ 1 = 1 (see section 4.5).

Specific heat
We compare our specific heat measurements to the frustrated bilayer model of Gu and Shen [35]. Their model contains a strong interlayer bond J 0 , a weaker intralayer bond J 1 , and a frustrating interlayer bond J 2 . This exchange topology is similar to that in PHCC. We include a magnetic field dependent term in their thermodynamic expansion via a linear Zeeman splitting of the triplet excited state as a function of applied magnetic field. We carried out a global fit of all our specific heat data with µ 0 H ≤ 9 T and T ≤ 30 K using only the free energy terms up to powers of 1/E −2 in Gu and Shen's expansion, one independent prefactor for each magnetic field, and global terms consisting of the three exchange constants, a g-factor and a lattice contribution ∝ T 3 . The results of this fit are shown as solid lines in figure 5 and account very well for the temperature and field dependence with J 0 = 2.08(1) meV, J 1 = 0.34(3) meV and J 2 = −0.13(3) meV The exchange parameters for the frustrated bilayer model determined from the specific heat data for PHCC can be used to evaluate the corresponding dispersion relation, equation 4 of Ref. [35]. Figure 12 depicts the calculated dispersion for PHCC along several high-symmetry directions in reciprocal space, compared to the zero-field triplon dispersion measured by inelastic neutron scattering [24,25]. The values derived from the fits to C m for the spin gap ∆ = 1.11(1) meV and bandwidth of 1.92(2) meV agree quite well with the measured spin gap, ∆ = 0.99 meV and bandwidth, 1.73 meV. Furthermore, the frustrated bilayer model reproduces the shape of the dispersion remarkably well given that specific heat data provided the only experimental input.
The lambda-type peaks in the magnetic heat capacity yield an additional measure of the boundary on the magnetic field vs. temperature phase diagram for PHCC, figure 11(b), between the quasi-2D RC-2D disordered phase and the LRO phase. However, C p was measured with H a * , while all the other data included in figure 11(b) were measured with H b. As noted above, differences in the critical fields measured for χ(H) for different field orientations are accounted for by g-factor anisotropy. After  figure 13. The two different power-laws observed for H < H N , as shown by the fits for 0.02 ≤ h ≤ 0.2 and the asymptotic change in power-law as one approaches H N is a potential indication of an additional disordered portion of the phase diagram for H < H N . This is in fact directly measured as a gapless disordered region through comparisons of the magnetic susceptibility and neutron scattering measurements as we discuss later in the text and as shown in figures 11(b) and (c).
We fit the H > H N data in figure 13 to a power law of the form where α is the specific heat exponent. A global fit for T = 2 K and 3 K results in α = 0.097(2) while individual fits yield α(T = 2K) = 0.104(3) and α(T = 3K) = 0.072 (2). Individual fits to Eq. 3 for H < H N over the range 0.02 ≤ h ≤ 0.2 yield α(T = 2K) = 0.175(2) and α(T = 3K) = 0.2626 (7). These values are larger than those obtained for H > H N , which indicates that the fits include data beyond the critical regime. Table 1 summarizes theoretical and experimental results for the 3D magnets.

Elastic neutron scattering
4.3.1. Ordered magnetic structure To determine the field induced magnetic structure we measured fourteen magnetic Bragg peaks in the (h0l) plane. The measurements were carried out at T = 1.65 K and µ 0 H = 13.7 T which is within the ordered phase as established by heat capacity measurements. Table 2 shows raw experimental results in the form of wave vector integrated intensities, I m (τ m ), for transverse (rocking) scans.
I m (τ m ) is related to the magnetic structure of the material as follows: Here C is a normalization factor that depends on aspects of the instrument such as the incident beam flux that are not considered in the elastic resolution function R τ (q) [40] which is sharply peaked around q = 0. N m and v m are the number of magnetic unit cells and the magnetic unit cell volume respectively. The magnetic Bragg scattering intensity is given by where r 0 = 5.38 fm, [41] is the magnetic structure factor, g d is the Lande g-factor for the magnetic ion on site d within the magnetic unit cell, and F d (Q) is the magnetic form factor of the Cu 2+ ion [42,41]. The magnetic structure enters through the expectation values for the spin vector, S d in units ofh. Table 2 compares the measured scattering intensity I m (τ m ) to the calculated intensity for the uniaxial magnetic structure that best fits the data and is shown in Values of (T N , H N ) extracted from the analysis appear as filled squares on Figure 11. There is a clear distinction between the onset of a state of high differential susceptibility as determined from χ(H, T ) and the boundaries of the long range ordered phase (LRO) extracted from neutron diffraction data. We denote the state of the system between these boundaries as 2D renormalized classical. Spin correlations in this phase were examined by quasi-elastic neutron scattering measurements along the ( 1 2 0l) and (h01 2 ) directions through the ( 1 2 01 2 ) wave-vector at T = 3 K and µ 0 H = 9 T. No magnetic Bragg peaks were encountered. Instead the broad peak centered at Q = ( 1 2 01 2 ) indicates Table 2. Measured and calculated relative magnetic Bragg peak scattering intensities at the indicated magnetic reciprocal lattice points τ m = ha * + lc * . The parentheses indicate statistical errors only. The minimum in the crystallographic reliability index R = |I obs m − I calc m |/ I obs m = 0.3 is found for the magnetic structure depicted in figure 1 and described  short range correlations with the same local structure as the long range order that sets in at slightly higher fields.

Critical exponents
The order parameter critical exponent, β, determined for each of the fits to the magnetic field and temperature dependent scans is shown versus temperature in the inset to figure 7. For T ≥ 0.5 K, the average value is β = 0.34 (2) which is indistinguishable from the critical exponent for the 3D Heisenberg [37], Ising (β = 0.326), and XY (β = 0.345) models [36]. However, the positive specific heat critical exponent, α, uniquely identifies the phase transition as being in the 3D Ising universality class. This is consistent with expectations for a phase transition to long range magnetic order in space group P1. We also calculate the remaining critical exponents based upon the Rushbrooke , α + 2β + γ = 2, Fisher, γ = ν(2 − η), Josephson, 2 − α = νd, and Widom, δ = 1 + γ β , scaling laws, assuming d = 3 dimensions [36,43]. The values are again consistent only with those of 3D Ising transitions as summarized in table 1. Recalling the exponents δ 1 = 1.022(2) and δ 2 = 1.8(1) in (2a) and (2b) determined from pulsed field magnetic susceptibility data, note that these values are different from δ derived from the scaling relations. This is because δ 1 and δ 2 are associated with the uniform susceptibility while δ in table table 1 characterizes the response to a staggered field conjugate to the Néel order parameter.
The apparent order parameter critical exponent increases as temperature is reduced, such that β = 0.40(1) for T ≤ 0.43 K (see inset to figure 7). A possible interpretation of this behavior is that the system enters a crossover regime associated with the approach to quantum criticality. Indeed, in the T = 0 limit the system should be above its upper critical dimension such that mean field behavior with a critical exponent β mf = 0.5 should be expected.
The field sweeps of the Q = ( 1 2 01 2 ) magnetic Bragg peak intensity shown in figure 9(a)-(c), indicate a reentrant phase transition. The H − T phase boundary obtained by fitting these data ( figure 11(b)) show a local minimum in the phase boundary between the LRO and QP phases. A possible origin of this anomaly is coupling to non-magnetic degrees of freedom such as nuclear spin or lattice distortions that are unimportant at higher temperatures but become important sufficiently close to the quantum critical point [27]. Figure 10 shows the magnetic field and temperature dependence of the excitation spectrum at the critical wave vector Q = ( 1 2 01 2 ) for 0 ≤ H < H c2 . We fit each spectrum to a sum of Gaussian peaks to account for the magnetic excitations with a constant background and a field independent Gaussian+Lorentzian lineshape centered at zero energy to account for incoherent elastic nuclear scattering. The field dependence of the peak locations are plotted for the T = 0.06 K and T = 4.5 K data in figure 10(b). The data demonstrate the triplet nature of the excited state through its Zeeman splitting. Above H c1 , the rate of change of the peak position as a function of field for the higher energy peak is twice that of the lower energy excitation, indicating a magnetized ground state. We fit the low T field dependent peak locations to the spectrum of an antiferromagnetically coupled spin-1/2 pair with singlet-triplet splitting ∆ subject to a magnetic field described by H Z = −gµ B µ 0 H. Results of this two-parameter fit to the low-temperature peak locations are plotted as solid lines in figure 10(b), and yield values of ∆ = 1.021(2) meV, g = 2.33(1), and µ 0 H c1 = ∆/(gµ B ) = 7.58(3) T, which are in agreement with zero-field inelastic neutron scattering and bulk magnetization data.

Gapped Excitations
The average resolution-corrected half width at half maximum, Γ, of the three Gaussians that were fitted to each spectrum for T = 0.065 K and T = 4.5 K are shown in figure 14. While Γ, considering systematic uncertainties in the energy resolution, is indistinguishable from zero at low T , it increases with applied field at elevated temperature. Also increasing with field and temperature is the magnon density. Hence the result is consistent with an increased rate of magnon collisions as their density increases. Self consistent mean field theories have been developed to account for these effects. [44] The inset to figure 14 shows the field dependent ratio of the integrated intensities of the central Gaussian peak to that of the satellite peaks. Based upon the polarization factor and the Wigner-Eckart theorem as applied to a Zeeman split triplet, this ratio should be 0.5 for H < H c1 [45]. Both the low-and high-temperature constant wavevector scans are consistent with this result. Above H c1 , Bose condensation modifies the ground state and the experiment shows a decrease of intensity for the highest energy mode. We use inelastic neutron scattering to examine the crossover between the low-field gapped (QP) and gapless (RC-2D) portions of the phase diagram. The T = 0.06 K lower critical field determined by comparing the peak locations in low-temperature constant energy scans to Eq. ?? is plotted in figure 11(a) and (b) (solid diamond at T = 0.06 K). With wave vector transfer at Q = ( 1 2 01 2 ) and energy transfer fixed athω = 0.3 meV the data in figure 15 show the field dependent intensity for various fixed temperatures. At the lowest temperature a sharp peak indicates the field at which the S z = 1 mode passes through the 0.3 meV energy window of the detection system. With increasing temperature and hence increasing magnon density not only does the peak broaden as expected from an increased relaxation rate, it also moves to higher fields indicating a finite slope for the crossover line in the H − T plane. Just as for the finite magnon lifetime, the origin of this effect is magnon repulsion. The increasing background is likely due to population of low-energy gapless excitations which develop for magnetic fields above the lower critical field. By extrapolation to thehω = 0 line in Fig. 15(b), these data were used to determine the temperature dependence of the field required to close the spin gap and the corresponding data points are included in figure 11(a) and (b) (solid diamonds). The coincidence between the crossover field determined in this manner and the fields associated with differential susceptibility anomalies establish a not unexpected association between effectively closing the gap at the critical wave vector and maximizing the differential uniform susceptibility.

Magnetic Excitations in the high field paramagnetic state
The previous results indicate a well defined crossover in PHCC from a low field and low temperature state dominated by a gap in the excitation spectrum to a high field state with an excitation spectrum that has finite spectral weight at zero energy. The crossover clearly occurs at lower fields than required to induce long range magnetic order. There are possible analogies to certain S = 1/2 1D Heisenberg antiferromagnetic systems including the linear chain, the next nearest neighbor frustrated chain, and the alternating chain where a low-energy continuum of excitations is predicted to exist upon closure of the spin-gap [45,46,47]. Calculations examining both gapped and gapless S = 1/2 spin-orbital chains describe an incommensurate gapless continuum for large portions of exchange and orbital interaction parameter space [48]. A characteristic of these continua are soft incommensurate modes moving through the Brillouin zone from the antiferromagnetic zone center to the ferromagnetic point as a function of applied magnetic field [31,49,50,51]. Incommensurate magnetic field dependent soft modes are also found in numerical calculations for the 2D antiferromagnetic square lattice with next nearest neighbor interactions, a system with similar antiferromagnetic interactions as PHCC [34].
An important distinction between PHCC and the systems mentioned above lies in the fact that for PHCC the putative gapless state so far was only observed at finite temperatures. Ascribing this to the inevitable effects of inter-plane coupling, the analogue of which occurs via interchain coupling in quasi-one-dimensional systems, there exists a potential for field induced incommensurate correlations in PHCC. Here we present preliminary experiments at the highest fields presently accessible to neutron scattering. The field was µ 0 H = 14.2 Tesla and the temperature chosen was T =  Fig. 16 (a) are difference data indicating excess magnetic spectral weight at the antiferromagnetic zone center ( 1 2 01 2 ) as compared to ( 1 2 01). The solid line shows a fit to the following resolution convolved Lorentzian relaxation spectrum, which obeys detailed balance and features a relaxation rate,hΓ = 0.025 ± 0.2 meV. The spectrum extends to and is approximately symmetric about zero energy, which is consistent with a 2D renormalized classical phase close to the phase transition to commensurate magnetic order. Figure 16 (b) shows two constant energy scans athω = 0.36 meV andhω = 0.55 meV that indicate short range commensurate magnetic correlations. The fitted Lorentzian peaks, which are much broader than the experimental resolution (solid bar) correspond to a dynamic correlation length scale of 3.2 ± 2.5Å and 4.2 ± 2.1Å for thehω = 0.55 meV andhω = 0.36 meV measurement respectively. The magnetization of PHCC is only 20% of the saturation magnetization in these measurements so any incommensurate shift, which should be a similar fraction of the distance to the Brillouin zone center is not clearly visible in the present data.

Bose condensation of magnons
In this section the phase transition to long range magnetic order in PHCC is analyzed from the perspective of magnon condensation [52,53,54,55]. Following Nikuni et al. [54] and Misguich and Oshikawa [56] we treat magnons as bosons whose density n is proportional to the deviation of the longitudinal magnetization from its plateau value M z0 : gµ B n = |M z − M z0 |. The low-density magnon gas is described by the Hamiltonian where µ = gµ B (H − H c1 ) is the boson chemical potential and v 0 is a short-range repulsion. The energy dispersion ǫ k is assumed to have the global minimum at k 0 = [ 1 2 0 1 2 ]. For convenience, excitation energies are measured from the bottom of the band, so that ǫ k 0 = 0.
At zero temperature the system acquires a condensate a 0 = α √ V when µ > 0. The energy density of the ground state in the presence of a condensate (neglecting quantum corrections) is v 0 |α| 4 − µ|α| 2 , which is minimized when |α| 2 = µ/2v 0 . In this approximation the boson density is n is related to the longitudinal magnetization At a finite temperature we obtain an effective Hamiltonian where n is the density of thermally excited bosons. The bosons condense when the renormalized chemical potentialμ ≡ µ − 2v 0 n vanishes. This yields the critical density The critical field is then The boson scattering v 0 can be extracted from low-temperature magnetization curves. The magnon density per unit cell n is equal to the ratio M/M sat (there is exactly one magnon per unit cell when M = M sat ). Based on the magnetization curve at the lowest available temperature T = 0.46 K we obtain v 0 = 1.91 meV × the unit-cell volume (u.c.v.). In other words, if we model the magnons as bosons on a cubic lattice with anisotropic hopping and on-site interactions U, the interaction strength would be U = 2 meV. The linear behavior extends over a range of some 10 T, or 13 K in temperature units.
Using this value of the boson coupling we can test the predicted relation between the critical field and critical magnon density (13) at finite temperatures from the data for n = M/M sat in the temperature range 0.46 K ≤ T ≤ 6 K. From these we extracted pairs (M, H) whose H(T ) values were sufficiently close to the line of phase transitions H c (T ). The resulting graph is shown in Fig. 17. The zero-temperature critical field The line H c (T ) is obtained from Eqs. (12) and (13). The magnon dispersion is taken from Stone et al. [24] (the 6-term version). That dispersion relation ǫ 2D k is strictly twodimensional, which precludes Bose condensation at a finite temperature. We therefore added weak tunnelling between the planes: The experimental data [57] give an upper bound for the total bandwidth in the [010] direction: 4τ < 0.5 meV. When T ≪ τ , only the bottom of the magnon band has thermal excitations and one may replace the exact band structure with a parabolic dispersion to obtain where m 3D = (m a m b m c ) 1/3 is the 3D effective mass. In the limit of a very weak interplane tunnelling we can also access a regime where τ ≪ T ≪ W , where W ≈ 2 meV is the magnon bandwidth. In this regime the in-plane dispersion can still be treated as parabolic. This yields the critical density for a quasi-2D Bose gas [58] n c (T ) = m 2D T where m 2D = (m a m c ) 1/2 and b is the interplane distance. In the intermediate regime T = O(τ ) the critical density first rises faster than T 3/2 [56] before crossing over to the T log T behavior. Curves for n c (T ) were obtained by integrating Eq. (12) numerically for several values of the interplane tunnelling τ . Since v 0 is already fixed, there are no further adjustable parameters. The data below 2.5 K (H c < 9 T, n c < 0.05 per unit cell) are best described by the curves for τ = 0.03(1) meV. The results are shown in Fig. 18.
The theory based on the Hamiltonian (11) assumes that the boson interaction v 0 is weak. However, it is also applicable to strongly interacting bosons, as long as the interaction is short-range and the gas is sufficiently dilute. In that case v 0 is not the bare interaction strength (which can be large) but rather the scattering amplitude (the t-matrix) for a pair of low-energy bosons. For hard-core bosons on a cubic lattice with a nearest-neighbor hopping amplitude t, one obtains v 0 ≈ 8t [59]. In our case t ≈ W/8 = 0.25 meV giving v 0 ≈ 2 meV. Another estimate along these lines, using the "exact" magnon dispersion yields, via yields a similar value of v 0 = 2.5 meV. These numbers demonstrate that the inferred interaction strength of 1.91 meV u.c.v. is quite reasonable. At higher temperatures the agreement worsens. At T = 3.7 K-the highest temperature where the condensation is observed-the mean-field theory predicts H c = 11 T, well below the observed value of 14.2 T. The observed magnon density at the critical field is 0.19 per unit cell (the theoretical estimate is 0.10). It is possible that the mean-field theory of magnon condensation, appropriate for a dilute Bose gas of magnons, becomes unreliable due to overcrowding. (At a technical level, the assumption of pairwise interactions may break down: two particles no longer collide in a vacuum.) Another possible reason for this failure is an increased strength of (thermal) critical fluctuations.

Kosterlitz-Thouless scenario
If the planes in PHCC were fully decoupled from each other the Bose condensation would change into the Kosterlitz-Thouless (KT) vortex-unbinding transition. It does not seem likely that the data can be reconciled with this scenario: in d = 2 spatial dimensions the crossover exponent φ = 1, so that the line of KT transitions H c (T ) = H c (0) + CT φ is expected to be linear.
The quantum criticality of the two-dimensional Bose gas was recently examined by Sachdev and Dunkel [60] who used a combination of scaling arguments and numerics to determine H c (T ) for a weakly interacting Bose gas (11). (Analyzing the RG flow alone does not solve the problem: the KT fixed point lies in a strongly coupled limit.) The RG analysis focuses on two variables: the (renormalized) chemical potentialμ and the two-particle scattering amplitude v. Unlike in d = 3 dimensions, the scattering amplitude vanishes in the limit of zero energy. Therefore the role of the coupling constant is played by the t-matrix for particles with energies O(T ): where Λ is a high-energy cutoff on the order of the magnon bandwidth. While bothμ and v flow under RG transformations, a certain combination of these variables, namely T v/μ, remains unchanged. The value of this constant determines whether the system flows to the disordered phase or to the phase with long-range correlations. The KT transition corresponds to 2mT v h 2μ = −K ≈ −34.
At the KT transition, the renormalized chemical potential is We see that |μ| is a bit below T (particularly if T ≪ Λ), which means that there are many thermally excited particles. The magnetic field is connected to the bare chemical potential: The density is computed for an ideal gas: This yields the final result for H c (T ): It depends weakly on the ultraviolet cutoff Λ. For Λ = 2.5 meV, we have µ ≈ 1.5T in the temperature range of interest. The curve is labelled 'KT1' in Fig. 18.
The result (23) can also be obtained directly from Eqs. (3) and (4) of Sachdev and Dunkel [60] if one takes the classical limit µ ≪ T . The value µ ≈ 1.5T suggests that we are (just) outside the classical regime and one needs to solve the Sachdev-Dunkel equations numerically. The resulting curve (for Λ = 2.5 meV) is shown in Fig. 18 as 'KT2'. It is also approximately linear, with a slope dµ/dT ≈ 1.
Even with the refined approach, the KT theory does not seem to agree with the data. As anticipated, it produces a linear dependence H c (T ) (up to logarithmic corrections) instead of a higher power law. The rise of the data above the KT line at high temperatures may signal the inadequacy of the continuum approximation at high magnon concentrations.

Conclusions
Originally thought to be a weakly alternating spin chain, our measurements establish PHCC as a quasi-two-dimensional spin-dimer system. The gap is sufficiently small to allow access to the lower critical field and the bandwidth sufficiently large to produce several distinct phases.
The long range Néel ordered phase is a commensurate antiferromagnet with moments perpendicular to the applied field direction. The critical exponents for the phase transition are similar to those of a 3D Ising model. The transition into this phase can be described as magnon Bose-Einstein condensation. Mean-field (Hartree-Fock) theory does an adequate job describing the magnetization M(H) at the lowest temperatures (T = 0.46 K) and the shape of the critical line H c (T ) for T < 3 K. The theory predicts a small but nonzero magnon dispersion of 0.08 ≥ 4τ ≤ 0.16 meV perpendicular to the 2D planes. The higher power of the H c (T ) dependence is caused by deviations of the magnon density of states from the E 1/2 behavior of the continuum theory. The extra DOS (due to a van Hove singularity) is positive, hence an increase over the T 3/2 law. The continuum KT theory gives an approximately linear dependence H c (T ) and thus does not seem to do a good job.
Close to the lower quantum critical point there is an anomaly in the phase boundary that is not accounted for by the BEC theory (see figure 11(b) and figure 9(a)-(c).) Low T anomalies near the zero temperature termination of phase boundaries are also seen experimentally in the quantum liquids 3 He[61] and 4 He [62,63], the frustrated spin system GGG [7,8], and in numerical and theoretical results for the BCC Ising antiferromagnet [64], the Ising antiferromagnet on a cubic lattice [65], and S = 1/2 XXZ 1D- [66] and 2D- [67] antiferromagnets. In 4 He the anomaly is associated with there being more low energy degrees of freedom in the solid than the liquid and a low−T anomaly in the phase boundary of LiHoF 4 has been associated with nuclear spin degrees of freedom. [68] For PHCC we speculate that features beyond a spin model such as the relatively soft lattice or nuclear spins become relevant sufficiently close to the quantum critical point. Such a mechanism which appears to be rather general merits further exploration through experiments that probe such auxiliary degrees of freedom directly.
The high field paramagnetic phase is perhaps the most intriguing aspect of PHCC. The crossover from the quantum paramagnet to this phase remains remarkably sharp upon heating and the experiments show that a finite temperature pseudo spin gap closes there. Further analysis is required to understand why this crossover remains so distinct. A proper examination of the gapless high field phase requires inelastic neutron scattering magnetic fields of order 30 Tesla where the system is more heavily magnetized. This work must await development of a higher field facility for inelastic neutron scattering.