Materializing rival ground states in the barlowite family of kagome magnets: quantum spin liquid, spin ordered, and valence bond crystal states

The spin-12 kagome antiferromagnet is considered an ideal host for a quantum spin liquid (QSL) ground state. We find that when the bonds of the kagome lattice are modulated with a periodic pattern, new quantum ground states emerge. Newly synthesized crystalline barlowite (Cu4(OH)6FBr) and Zn-substituted barlowite demonstrate the delicate interplay between singlet states and spin order on the spin-12 kagome lattice. Comprehensive structural measurements demonstrate that our new variant of barlowite maintains hexagonal symmetry at low temperatures with an arrangement of distorted and undistorted kagome triangles, for which numerical simulations predict a pinwheel valence bond crystal (VBC) state instead of a QSL. The presence of interlayer spins eventually leads to an interesting pinwheel q = 0 magnetic order. Partially Zn-substituted barlowite (Cu3.44Zn0.56(OH)6FBr) has an ideal kagome lattice and shows QSL behavior, indicating a surprising robustness of the QSL against interlayer impurities. The magnetic susceptibility is similar to that of herbertsmithite, even though the Cu2+ impurities are above the percolation threshold for the interlayer lattice and they couple more strongly to the nearest kagome moment. This system is a unique playground displaying QSL, VBC, and spin order, furthering our understanding of these highly competitive quantum states.

A. Barlowite 1 To resolve open questions in the literature regarding the low-temperature structure of barlowite, we have undertaken a comprehensive investigation of barlowite 1. Previous neutron powder diffraction (NPD) reports assign it to orthorhombic space group Pnma (No. 62) [34,35,41], but others conclude that space group Cmcm (No. 63) is plausible. [40] We employ a combination of high resolution synchrotron powder and single crystal X-ray diffraction (PXRD and SCXRD) as well as neutron powder diffraction (NPD) to definitively determine that our sample of barlowite 1 has a reversible phase transition at T ≈ 265 K to Pnma. Figure 1b in the main text shows the emergence in synchrotron PXRD of superlattice peaks and orthorhombic peak splitting related to this transition, which is confirmed by low-temperature synchrotron SCXRD and Rietveld co-refinements of PXRD and NPD data (see  Tables 2-3, 8).
For several samples similar to barlowite 1 but synthesized using different reagents and temperature profile-which can affect the magnetic properties [33]-the recently proposed lowtemperature Pnma orthorhombic structures based on NPD contain only one (off center) interlayer Cu 2+ site, which is fully occupied. [34,35] However, by combining high-resolution synchrotron SCXRD data with co-refined PXRD and NPD data, we find a different model of the interlayer Cu 2+ . We determine that in barlowite 1 the interlayer Cu 2+ s become disordered over three symmetry inequivalent sites with unequal occupancies (Figure 1a in the main text). The much higher resolution of synchrotron SCXRD and PXRD compared to NPD lends credence to the accuracy of this model; using a model with only one interlayer site worsens the R 1 of our SCXRD model from 2.23% to 14.64%. The occupancies of these three disordered interlayer sites vary slightly in different measurements: they are (52%, 33%, 15%) from the T = 15 K SCXRD refinement and (83%, 11%, 5%) from the low-temperature (T NPD = 2 K and T PXRD = 13 K) co-refinement (Tables 3, 8). This difference can be ascribed to small variance in the crystals. The overarching trend, however, is that one site dominates the occupancy-likely why the previous NPD models only identify one interlayer site.
While we are in broad agreement with the literature as to the Pnma symmetry of orthorhombic barlowite, [34,35,41] it is unsurprising that there may be inconsistencies between structural models measured for different samples given the observed dependence of barlowite's physical properties upon factors such as reactant stoichiometry and synthetic conditions. [33] The proposed Cmcm symmetry of a sample similar to barlowite 1 was determined using several local probe measurements at room temperature. [40] Our room temperature P6 3 /mmc model is consistent with the 19 F magic angle spinning (MAS) NMR data, and our previously published FTIR spectroscopy data is also consistent with the reported IR data. [33,40] However, the 1 H MAS NMR and electron diffraction data are incompatible with P6 3 /mmc. [40] Attempting to consolidate structural models for samples synthesized in distinct ways may not be valid. Local probe measurements such as solid-state 1 H MAS NMR would be required to ensure that each sample has true P6 3 /mmc symmetry.

B. Barlowite 2
All structural measurements support the lack of an orthorhombic transition in barlowite 2; superlattice peaks related to hexagonal space group P6 3 /m are observed in SCXRD ( Figure 1) and in neutron scattering ( Figure 1d in the main text and Figure 20). No evidence of orthorhombic splitting is observed in the PXRD data, implying hexagonal symmetry is preserved (Section III). SCXRD and PXRD refinements on the same sample at different temperatures (Tables 3, 9) show slight variation in the relative occupancies of the interlayer sites, ranging from (43%, 33%, 24%) from SCXRD at T = 15 K to (42%, 31%, 27%) and (44%, 30%, and 27%) from PXRD at T = 15 and 99 K, respectively. These deviations are much smaller than those observed in barlowite 1 and further support the lower degree of symmetry lowering observed in this subtle hexagonal phase transition compared to barlowite 1's orthorhombic transition.

C. Zn-substituted barlowite
Both compositions of Zn-substituted barlowite (Zn 0.95 and Zn 0.56 ) crystallize in space group P6 3 /mmc at all temperatures measured. Slight variations in the site occupancies of the two interlayer sites (centered D 3h and the set of three off-center C 2v sites observed in all-Cu barlowite) occur in both samples. The following discussion assumes that the centered site is occupied by Zn 2+ and the off-center sites are occupied by Cu 2+ , consistent with their different Jahn-Teller activities.
SCXRD refinements on the same crystal of Zn 0.56 at different temperatures show slight variation in the relative occupancies of the interlayer sites, ranging from (59% Zn, 41% Cu) at T = 300 K to (50% Zn, 50% Cu) at T = 100 K (Tables 2, 4). Rietveld refinements of PXRD data on the sample show similar deviations, ranging from (62% Zn, 38% Cu) at T = 90 K to (58% Zn, 42% Cu) at T = 295 K (Table 12). The average formula for Zn 0.56 from all diffraction measurements is Cu 3.43 Zn 0.57 (OH) 6 FBr, consistent with the ICP-AES results.
Rietveld co-refinements of NPD and synchrotron PXRD data of Zn 0.95 were performed assuming full occupancy of Zn on the interlayer site. The single crystal measured in SCXRD has (85% Zn, 15% Cu) on the interlayer. Although this is reasonably close to the bulk ICP-AES results (Cu 3.05 Zn 0.95 (OH) 6 FBr), this crystal may not be an accurate representation of the bulk (see Section II B). Low-temperature synchrotron X-ray diffraction measurements were performed on single crystals of barlowite 1, barlowite 2, and Zn-substituted barlowite Zn 0.56 and Zn 0.95 at T = 14 K or T = 100 K; Zn 0.56 was also measured at T = 300 K (Table 2). Barlowite 1 and 2 were previously measured at T = 300 K; the details of these measurements and crystallographic information files (CIFs) can be found in Ref. [33]. Crystallographic information is tabulated in Table 3 for barlowite 1 and 2 and in Table 4 for Zn-substituted barlowite. Selected bond angles and distances are found in Tables 6-7. A. Barlowite 2 Precession images for barlowite 2 at T = 300 K and T = 15 K are displayed in Figure  1. All expected superlattice peak positions arising from space group P6 3 /m are circled in red; the presence of peaks at these positions indicates the presence of this low-temperature phase transition. We do observe a very weak [0 0 l] (l = 2n peak in barlowite 2 at low temperature, although the reflection conditions for both P6 3 /m and P6 3 /mmc are 00l, l = 2n. Their presence suggests that the stacking has subtle imperfections or that there may be very weak doubling of the c-axis: on average, the [0 0 -1] peak is 5500 times weaker than the allowed [0 0 ±2] peaks and 8000 times weaker than the [±1 0 0] peaks. Of the two equivalent reflections, one has an intensity within error of 0 (0.01(1)) and the other's is very close (0.04 (2)).

III. POWDER X-RAY AND NEUTRON DIFFRACTION OF BARLOWITE
Synchrotron powder X-ray diffraction (PXRD) was measured on barlowite 1 and 2, and neutron powder diffraction (NPD) was measured on barlowite 1 at a variety of temperatures. Rietveld co-refinements of the PXRD and NPD data were performed on barlowite 1 (Figures 2-4; Table 8). For barlowite 1 below the phase transition, a mixture of two phases was found in the synchrotron PXRD data-both the low-temperature orthorhombic phase (Pnma) and the high-temperature hexagonal phase (P6 3 /mmc). The hexagonal phase fraction refined to approximately 30%. This was not observed in the NPD data, so only the orthorhombic phase was refined.
Rietveld refinements of the PXRD data of barlowite 2 were performed (Figures 5-6; Table 9). Superlattice peaks due to space group P6 3 /m are not visible ( Figure 1c in the main text and Figures 5-6), suggesting that they are too weak to detect when crystals are ground into a powder. [68] This supports the claim that barlowite 2 is not orthorhombic. In addition, an analysis of orthorhombic distortion was performed on the PXRD data of barlowite 1 and 2 (protonated and deuterated of each) at the lowest temperature measured, T = 13 K. The orthorhombic splitting (b − a)/(b + a) was calculated for three in-plane peaks ([1 1 0], [2 1 0], and [2 2 0], using the notation of the high-temperature space group), where a and b are the positions in q (Å −1 ) of the left-and right-hand peaks, respectively. For barlowite 2, these positions were picked within the peak (which does not split and exhibits very little to no broadening in FWHM from T = 300 K to T = 13 K) and represent an upper limit of possible splitting. The average distortion for barlowite 1 is 0.00137, while the upper limit for the average distortion in barlowite 2 is 20 times smaller (0.000072), demonstrating that barlowite 2 remains hexagonal.

V. PHYSICAL PROPERTIES OF BARLOWITE AND ZN-SUBSTITUTED BARLOWITE
Curie-Weiss parameters are tabulated in Table 15 for bulk samples and Table 16 for aligned single crystalline 2 and Zn 0.56 . Aligned magnetization on a single crystal of barlowite 2 at T = 2 K and T = 20 K is shown in Figure 12. Figure 11 shows the aligned magnetization in the ab plane at T = 2-12 K; the curves below T = 6 K display a different slope than the curves where T ≥ 6 K with a crossover at approximately µ 0 H = 0.75-1 T. This change in slope coincides with the lower-temperature transition (T N2 ) observed in AC susceptibility (Figure 4a in the main text) and supports the hypothesis that different spins are primarily responsible for each magnetic transition. DC susceptibility and inverse susceptibility measurements on bulk samples of barlowite 1, barlowite 2, Zn-substituted barlowite Zn 0.95 and Zn 0.56 , and herbertsmithite at several applied magnetic fields are shown in Figures 13-14 An estimate of the susceptibility of the intrinsic barlowite kagome lattice is also plotted. Aligned single crystalline susceptibility and inverse susceptibility data of barlowite 2 and Zn 0.56 are shown in Figure 15. Molar and magnetic heat capacity (C and C mag , respectively) measurements on barlowite 2 are shown in Figure 16. Dilution refrigerator heat capacity measurements on bulk Zn 0.95 and single crystalline Zn 0.56 are shown in Figure 17. A comparison of the heat capacity of all samples is shown in Figure 18.    Deuterated crystals of barlowite 2 were synthesized for elastic neutron scattering measurements. Figure 19 shows the real part of the AC susceptibility (χ') and DC magnetization (M ) as a function of temperature for one representative batch of crystals. The magnetic behavior is slightly different from that of protonated barlowite 2 (Figure 4a in the main text): the higher-temperature transition (T N1 = 10 K) is more a shoulder instead of a distinct peak. The lower-temperature transition (T N2 = 6 K) does not change significantly.

VI. ELASTIC NEUTRON SCATTERING OF BARLOWITE 2
A. Background subtraction and fitting methods for the superlattice peak The notation of the high-temperature space group P6 3 /mmc is used to index peaks. In reciprocal space, θ-2θ scans and θ scans are along the longitudinal and transverse directions, respectively, of the wavevector Q in the scattering plane. Figure 1d in the main text shows the temperature dependence of longitudinal (θ-2θ) scans of the [0.5 1.5 0] superlattice peak, which are fit well by a Gaussian peak with a sloping background. The background, peak center and peak width are fixed to the best overall value based on all temperatures. Figure 20a shows the temperature dependence of transverse (θ) scans of the [1.5 0.5 0] peak, which is an equivalent position of [0.5 1.5 0]. They are superlattice peaks that are allowed in P6 3 /m but forbidden in P6 3 /mmc. Transverse (θ) scans at base temperature (1.6 K) and the highest temperature measured (283 K) are shown in Figure 20b. Although there is still a bump in the T = 283 K data, its center position deviates significantly from the peak position expected from the lattice parameters, which is denoted by the vertical dashed line. Thus, the T = 283 K data is treated as a temperature independent background and is subtracted from the lower temperature scans. The T = 1.6 K scan has the best statistics and is fit well by a Gaussian peak. The peak center is fixed to the expected peak position, and the peak width is fixed to the best overall value based on all temperatures. The fitted integrated intensities from both data sets are shown in the inset of Figure 1d in the main text after scaling to match each other. The structural phase transition temperature derived from linear fits of both data sets is T = 262(8) K.

B. Peak broadening of superlattice nuclear Bragg peaks and magnetic Bragg peaks
Neither the superlattice peaks nor the magnetic Bragg peaks at half-integer positions are resolution limited; however, the integer magnetic Bragg peaks are limited by the experimental resolution. Figure 21 shows representative scans of different types of magnetic Bragg peaks. Within a finite-size domain model, [42] the peak intensity is given by: These peaks are well-described by a Gaussian line shape, so we can make a Gaussian approximation of the above equation to match the peak maxima and integrated intensities. Assuming a cubic crystal (a 1 = a 2 = a 3 ) and isotropic crystallite (N 1 = N 2 = N 3 ), the domain size is related to the FWHM of the peak by: Therefore, we can extract the domain size from the fitted peak width. It is noted that the experimental resolution should be subtracted quadratically. The calculated domain sizes from  Table 17. For magnetic Bragg peaks, there are two coexisting distinct domain sizes: the half-integer peaks have a short domain size of ∼118Å on average, while the integer peaks have a much larger domain size of at least ∼360Å, which is limited by the experimental resolution.

C. Resolution function calculation and magnetic structure refinement method
For both nuclear and magnetic Bragg scattering, the neutron differential cross section can be written as: [69] ( where F is the static structure factor, v is the unit cell volume and A is a constant that depends on the sample geometry, incident flux, sample volume, and counting time. If these factors are all held fixed in a series of measurements, then A is simply an overall scale factor. For a material with a known crystal structure, one can determine the overall scale factor A by measuring nuclear Bragg peaks. Then a magnetic structure model can be fit to the measured magnetic Bragg peaks. However, the measured intensity is not simply the differential scattering cross section but rather a convolution of the differential cross section with the instrumental resolution: [69] where S( Q) = ( dσ dΩ ) Bragg is the differential cross section, and R( Q − Q 0 ) is the resolution function. The resolution function depends on the specific spectrometer configuration, the sample mosaic spreads, the sample shape, and the wavevector Q. In the Gaussian approximation, the resolution function can be expressed as a 3-dimensional (3D) Gaussian distribution: where M = is a 3 × 3 matrix and ∆ Q = Q − Q 0 = (q x , q y , q z ). The conventional coordinate system is the Q-system, in which x is along Q 0 , z is vertical, and y completes the orthogonal right-handed basis set. The resolution function can be visualized as a 3D ellipsoid centered at Q 0 . In the paraxial approximation, the resolution in the vertical direction is uncoupled from the other two in-plane coordinates. Hence the matrix elements M xz = M zx = 0 in M . For a  nuclear Bragg peak, assuming a perfect crystal, we can write: where S 0 is the amplitude of the scattering cross section at the Bragg peak Q b . Then the measured intensity is just the product of S 0 and the resolution function R( Q b − Q 0 ). For a scan through a Bragg peak Q b , the integrated intensity is given by: Thus, the integrated intensity not only depends on the differential cross section at the Bragg condition but also on the integral of the resolution function along the path of the scan in the reciprocal space. Note that Eqn. 5 can only apply to a perfect crystal. For short-range ordered structures such as the magnetic order in barlowite 2, the delta function approximation (Eqn. 4) no longer holds. However, we can still use Eqn. 5 with a correction factor that is based on three simple assumptions. First, the cross section S( Q) can be modeled by a 3D Gaussian function centered at the Bragg peak Q b : Second, the resolution function R( Q− Q 0 ) can be approximated by R( Q− Q b ), which-similar to S( Q)-is peaked at Q b . We assume that the resolution function varies slowly in the small region of a scan. Third, the off-diagonal elements M xy = M yx . = 0 in the matrix M because the major axis of the resolution ellipsoid is approximately aligned with the y-axis due to the large sample mosaic. Since the convolution of two Gaussian functions is also a Gaussian, we can obtain a formula for integrated intensity similar to Eqn. 5 but with a correction factor that depends on the scan direction. For example, the correction factor for θ scan in the (H K 0) scattering plane is 1/ (M xx σ 2 x + 1)(M zz σ 2 z + 1). The resolution function calculation in this work was performed using the Python module NeutronPy. The nuclear structure factor calculation was performed using VESTA [74] with input crystal structure data from the SCXRD refinement results (Table 3). Since the SCXRD cannot unambiguously determine the positions of H (or D) atoms due to their negligible X-ray scattering lengths, the overall scale factor A was calculated by measuring the nuclear Bragg peak [0 0 2], whose intensity is negligibly affected by these atoms. The magnetic structure refinement was accomplished through a least squares minimization routine written in Python to fit the calculated magnetic Bragg peak intensities to the measured integrated intensities. For simplicity, the Debye-Waller factor was neglected; the Landé splitting factor was fixed to 2; and an isotropic Cu 2+ magnetic form factor was used. The measured intensities at equivalent positions were averaged to increase statistics. The small differences in base temperatures between different measurements was neglected based on the fact that the intensities saturate at T ∼ 3 K (see Figure 5b in the main text). The two magnetic peaks at integer positions ([1 0 0] and [1 1 0]) measured at SPINS were excluded from the refinement since these two peaks were measured at T = 4 K, a temperature at which the intensities are not saturated. The peak broadening correction factors for all half-integer peaks used the same mean domain size of 118 A, while the integer peaks used the same domain size of 360Å, obtained from Table 17.

D. Detailed refinement results of the magnetic models
Possible magnetic space groups were generated using the k-Subgroupsmag application on the Bilbao Crystallographic Server. [70][71][72][73] As shown in Table 18, the highest-symmetry magnetic space group allowing a net in-plane moment (necessitated by the measured bulk magneti-  Figure 4c in the main text as well as Figures 11-12) is P2 1 /m with lattice parameters 2a × 2b × c, which is the same as the low-temperature nuclear unit cell. Its symmetry operations reduce the number of independent Cu 2+ sites in a unit cell to 10 (six kagome and four interlayer sites) and constrain the interlayer Cu 2+ moments to be within the ab plane. Since we measured a limited number of magnetic peaks (17 distinct peaks in total and 12 peaks after taking equivalent positions into account), we further reduce the number of fitted sites to six (four interlayer and two kagome). Table 19 shows the refinement results for the pinwheel q=0 kagome model and the orthorhombic model, which is based upon previous models for polycrystalline, orthorhombic samples similar to barlowite 1 [34,35] and is tested for completeness. Table 21 gives the measured and calculated intensities for each model. Both models give equally good fits, and we cannot unambiguously determine which is the best model from our data. However, the orthorhombic structure is less likely to exist in barlowite 2, which has less distorted kagome planes and should therefore prefer a q=0 type spin structure. Additionally, the orthorhombic model is not compatible with the underlying crystal structure, since the two types of kagome spins do not correspond to the two crystallographically distinct sites.
In the pinwheel q=0 kagome model discussed in the main text, we used 10 fitting parameters. There are two types of kagome spins where each has a distinct moment size and in-plane rotation angle, but they share an out-of-plane rotation angle. There is one interlayer spin moment size and four different in-plane rotation angles of the interlayer spins. Figure 22a shows the moment directions when the rotation angles are at zero. Counterclockwise is the positive direction for all in-plane rotation angles. The direction of the out-of-plane rotation for each kagome spin is illustrated in Figure 5c  fitting parameters (10) is close to the number of distinct peaks (12), and the three unequaloccupancy interlayer sites are along the [1 1 0]-type high-symmetry directions within error, we fix those spin directions in the final fit (Table 19). This choice of spin directions for these three interlayer sites leads to zero net moment from them, which reduces the overall net moment size. However, we cannot exclude the possibility that these interlayer spins deviate from the specified directions and contribute to the net magnetization.
In addition, the possibility that the kagome plane has net magnetization is also tested by introducing an in-plane canting angle for the green kagome spins that have an out-of-plane component. The positive direction for this canting angle is towards the direction of the green kagome spins that lie within the kagome plane. To reduce the number of fitting parameters, the moment sizes are fixed. The fit result (Table 20) shows a small but non-zero canting angle (−7 • ) with a large error bar (20 • ) and a tiny change to the χ 2 . No other rotation angles are significantly affected. Although this canting direction is opposite to the net moment direction from the interlayer spins, the χ 2 only worsens a small amount (from 39.5 to 40.2) if the canting angle is fixed to a positive value (Table 20). This positive canting of kagome moments contributes a net moment of 0.016 µ B per formula unit.
In the orthorhombic model, we used seven fitting parameters. There are two types of kagome spins where each has a distinct moment size and in-plane rotation angle, but they share one outof-plane canting angle. There is one interlayer spin moment size and one in-plane canting angle of the interlayer spins. Figure 22b shows the moment directions when the canting angles are fixed at zero. The positive direction for all the in-plane canting angles is towards the [1 1 0] direction. The direction of the out-of-plane canting for each kagome spin is illustrated in Figure  23b. Since our data is insensitive to the two in-plane canting angles of the kagome spins, we fix these angles to be zero in the final fit (Table 19). Both models have a smaller average ordered kagome moment compared to the proposed models for orthorhombic barlowite. It is markedly reduced from 0.32 µ B per Cu (Ref. [35]) to 0.19 µ B per Cu in both models. This is supported by the considerably weaker magnetic peaks of barlowite 2: the ratio between the magnetic and nuclear contributions to the [1 0 0] Bragg peak is ∼1/90, which is ∼13 times smaller than for orthorhombic barlowite 1. [34] Tables 22 and 25 compare the kagome moments' and interlayer moments' contributions at various magnetic Bragg peaks for the pinwheel q=0 kagome model and the orthorhombic model.  Like in the pinwheel q=0 model, in the orthorhombic model the kagome spins (the red and green spins in Figure 22b) have two moment sizes (0.08(1) µ B and 0.24(1) µ B ). However, this classification is not consistent with the two crystallographically distinct sites, which are denoted Cu1 and Cu2 in the data for barlowite 2 in Tables 3 and 9. The out-of-plane canting of the kagome moments also forms alternating stacking between layers, giving rise to the [1 1 1] magnetic peak. It is necessary to have two distinct kagome sites and cant the interlayer site: it matches the measured peak intensity of the [0.5 0 1] and [0.5 0 3] peaks better than using one kagome site does, and including the canting angle of interlayer Cu 2+ as a refined parameter significantly improves the fit, especially for the integer peaks (Tables 23-24). The interlayer spins have a moment size of 0.14(1) µ B . The in-plane canting of the interlayer spins also results in net FM with a net moment of 0.05 µ B per formula unit. The fit becomes much worse if the net moment size is constrained to 0.11 µ B per formula unit (see Tables 23-24). Based on the pinwheel q=0 kagome model, we can estimate the FM exchange couplings between the kagome and interlayer Cu 2+ s relative to the AFM couplings in the kagome plane. If we consider only NN couplings and assume that the coupling between the NN kagome spins is J (J > 0), then the total FM coupling between an interlayer spin and its six NN kagome spins is −2αJ − 4βJ when this interlayer site is fully occupied (Figure 24, inset). The average FM couplings for different interlayer sites are obtained by taking the site occupancies into account (see Tables 3 and 9).
FIG. 24: Contour map showing the energy difference between the pinwheel and the perfect q=0 kagome configurations (∆E = E pinwheel − E perfect ), as described in the text. The inset shows the exchange couplings considered in this calculation. Red and blue/white spheres represent the kagome and interlayer Cu 2+ ions, respectively. The couplings between the NN kagome spins are assumed to be the same (J, red lines). However, since the interlayer Cu 2+ ion randomly occupies one of the three possible off-center interlayer sites (blue sphere) with some probability, there are two types of exchange paths between this site and its NN kagome spins. These couplings are denoted by −αJ (thick blue lines) and −βJ (thin blue lines), respectively.
The classical-spin-based energy is investigated using the moment sizes from our pinwheel q=0 kagome model (Table 19). The out-of-plane rotation angle of the kagome spins is fixed to the value from our model, but the in-plane rotation angles of both the kagome and interlayer spins are treated as tunable parameters. Because the energy only depends on the relative angles between the spins, we fix the in-plane rotation angle of the kagome spin K1. Therefore, for arbitrary α and β, the best in-plane rotation angles of the other spins can be determined by searching for the minimum energy configuration. This minimum energy is also compared with the energy of the perfect q=0 spin arrangement to determine which state is more energetically favorable. This calculation is done for a grid of α and β in the range of 0 to 2 with a step size of 0.01, and the energy difference between the pinwheel and the perfect q=0 configuration is shown as a contour map in Figure 24. The colored region is where a) the pinwheel q=0 structure has lower energy than the perfect q=0 spin configuration, and b) the best in-plane rotation angles are close to the values from our neutron scattering refinement results (Table 19) with deviations no larger than 5 • (approximately the error bars in the refinement). Within this colored region, the lowest energy state is achieved at α = 1.39 and β = 0.64. Interestingly, the strong FM interactions obtained here match quite well with a previous report based on DFT calculations (α = 1.16 and β = 0.18). [31] The average FM couplings for the three interlayer sites are −0.96J (43% occupied), −0.89J (33% occupied), and −0.82J (24% occupied), respectively.
The robustness of the values of α and β against a change in moment sizes is also tested. When the kagome moment sizes are varied from 0.16 to 0.23 µ B (the values of the kagome moments from the refinement), the best kagome in-plane rotation angle shows a deviation from the neutron scattering refinement result (Table 19) smaller than 17 • .
We can attempt to understand the phase transitions upon cooling by modeling the effect of the interlayer spins L on the kagome spins S as an effective Zeeman field, i.e., J FM L· S ∼ h eff · S, where h eff ∼ J FM L . From the DMRG calculations, the kagome plane forms a VBC ordered ground state with a finite spin gap, e.g., ∆ ∼ 0.15J for J = 0.95J. Our experimental values for the interlayer moment (< 0.1 µ B ) in the intermediate temperature regime T N2 < T < T N1 yield an effective Zeeman field h eff < 0.04J, which is too small to close the VBC spin gap, consistent with the kagome moments' absence of order. However, the interlayer ordered moment (and, therefore, h eff ) continues to increase upon cooling. As a result, it becomes strong enough to destabilize the pinwheel VBC state and drive the system into the pinwheel q=0 state with longrange magnetic order.

VII. FIRST-ORDER PHASE TRANSITION BETWEEN THE PINWHEEL VBC AND QSL
We performed DMRG calculations of the ground state energy of the simplified model for barlowite 2 using the "hysteresis plot" method introduced in Ref. [5] to determine the nature of the phase transition. We find a first-order phase transition point between the pinwheel VBC phase and the QSL. The ground state energy is calculated along four paths of the Hamiltonian H(J ) for a 144-site 6-leg cylinder with two numbers of kept states, m = 1600 and 2000. A small step size 0.001J is taken in each of four adiabatic paths shown in Figure 25. Along the paths, the system shows a clear crossing at J c = 0.998 (1). This value is close to 1, but somewhat smaller, with a deviation larger than the resolution of our simulation.
FIG. 25: Plot of the ground state energy of the simplified model for barlowite 2 for a 144-site 6-leg cylinder for two values of the kept states m. For each of the four paths, the Hamiltonian parameter J is adiabatically tuned for every three DMRG sweeps to resemble the evolution of the system under changing J . A linear subtraction of the energy E with J has been applied.