Criticality-Enhanced Magnetocaloric Effect in Quantum Spin Chain Material Copper Nitrate

In this work, a systematic study of Cu(NO3)2·2.5 H2O (copper nitrate hemipentahydrate, CN), an alternating Heisenberg antiferromagnetic chain model material, is performed with multi-technique approach including thermal tensor network (TTN) simulations, first-principles calculations, as well as magnetization measurements. Employing a cutting-edge TTN method developed in the present work, we verify the couplings J = 5.13 K, α = 0.23(1) and Landé factors g∥= 2.31, g⊥ = 2.14 in CN, with which the magnetothermal properties have been fitted strikingly well. Based on first-principles calculations, we reveal explicitly the spin chain scenario in CN by displaying the calculated electron density distributions, from which the distinct superexchange paths are visualized. On top of that, we investigated the magnetocaloric effect (MCE) in CN by calculating its isentropes and magnetic Grüneisen parameter. Prominent quantum criticality-enhanced MCE was uncovered near both critical fields of intermediate strengths as 2.87 and 4.08 T, respectively. We propose that CN is potentially a very promising quantum critical coolant.


I. INTRODUCTION
Heisenberg spin chains and nets, owing to their strong quantum fluctuations and correlation effects, can accommodate plentiful interesting quantum phases like topological spin liquids [1,2], unconventional excitations like anyontype quasi particles [3], and inspiring behaviors like Bose-Einstain condensation in magnets [4], which continues stimulating both condensed matter theorists and experimentalists.What is more, these low-dimensional systems, which at a first glance are of purely academic interest, can actually have their experimental realizations.People have successfully discovered and/or synthesized plenty of spin materials which are very well described by the low-dimensional Heisenberg-type spin models.The long list includes, to name only a few, the diamond spin chain material azurite [5], the kagome spin liquid herbertsmithite [6], and copper nitrate as an alternating Heisenberg antiferromagnetic chain (AHAFC) .Therefore, low-dimensional quantum magnets arouses long-lasting research interest both theoretically and experimentally.
Among many other interesting properties of lowdimensional quantum magnets, we emphasize the enhanced magnetocaloric effect (MCE) in quantum critical regime.MCE is an intrinsic property of magnetic materials which exploits the reversible entropy changes caused by varying magnetic fields.MCE has a long history of study [29][30][31], and in the past decades, developing novel MCE materials which have prominent MCE properties, like the Gadolinium alloys with giant MCE [32,33], has raised great research interest.This is due to that MCE has appealing applications in eco-friendly refrigeration near room temperature [34,35], providing a good substitute to conventional vapor compression refrigeration, and also in space technology [36,37].In addition, MCE materials, in particular adiabatic dimagnetization refrigerant (ADR), serve as efficient coolants for ultra low temperature (sub-Kelvin) regime [38][39][40][41].People pursues MCE refrigerant which have higher isothermal entropy change (∆S), larger adiabatic temperature difference (T ad ), and also lower hysteresis dissipation [33].
Recently, quantum spin chain materials are shown to exhibit enhanced MCE even at ultra low temperatures, and thus raised great research interest [40][41][42][43][44][45][46][47][48][49].On one hand, through exploring low-T MCE properties of spin chain model materials [41,42] which shows divergent Grüneisen parameter near field-induced quantum critical points (QCP), people are able to directly detect and study quantum criticality [43,44].On the other hand, one can inversely utilize this low-temperature thermodynamic anomaly to realize enhanced cooling effects near quantum critical points (QCPs) [46,47].Very recently, Sharples et al. realized temperatures as low as ∼ 200 mK using the enhanced MCE of a molecular quantum magnet [40], and Lang et al. experimentally studied a spin-1/2 Heisenberg antiferromagnetic chain material [Cu(µ-C 2 O 4 )(4-aminopyridine) 2 (H 2 O)] n (CuP, for short) [48], and demonstrated this quantum critical coolant is a perfect alternative to standard ADR salts, due to its wider operating temperature range, longer hold time and high efficiency.[49] In order to study the thermodynamic information including the appealing MCE property of these strongly correlated arXiv:1607.04238v1[cond-mat.str-el]14 Jul 2016 spin systems, accurate thermal algorithms are of crucial significance, which is indispensable in establishing links between theoretical spin models and experimental measurements at finite temperatures.In one spatial dimension (1D), the transfer matrix renormalization group (TMRG) method [50] has been long accepted as the method of reference, owing to its high accuracy and versatility.In Ref. 52, Li et al. proposed an alternative approach for calculating thermodynamics of lowdimensional quantum lattice models called linearized tensor renormalization group (LTRG) method, which also adopts the Trotter-Suzuki decomposition [53] to express the partition function as a thermal tensor network (TTN) and linearly contract the resulting d + 1 dimensional (d = 1, 2 for 1D and 2D lattice, respectively) TTN along Trotter direction via renormalization group (RG) techniques.
In this work, combining three different methods, i.e., thermal quantum manybody computation, first-principles calculations, and experimental measurements of magnetization, we performed a comprehensive investigation of an alternating quantum spin chain material Cu(NO 3 ) 2 • 2.5H 2 O (copper nitrate hemipentahydrate, hereinafter referred to as "CN").CN is one of the earliest spin chain material ever studied experimentally 54], while continues intriguing people for its abundant physics including triplon wave excitation [23] and precise Tomanaga-Lutting liquid behavior [28].We notice that, despite many efforts, discrepancy in coupling constants still exists: the exact diagonalization (ED) fittings (J = 5.16 K, α = 0.27) to thermodynamic measurements measurably deviates from those obtained from inelastic neutron scattering (INS) experiments (J = 5.14 K, α = 0.24) [23].
In order to resolve this discrepancy in couplings, we generalize the LTRG approach to a bilayer formulation (dubbed as LTRG++) which further improves the accuracy of calculations.With this cutting-edge TTN method at hand, we revisit the previous experimental data in Ref. 22 including specific heat curves (at various fields) and magnetization curves, augmented with magnetization measurements done in this work.The couplings are determined precisely to be J = 5.13 K, α = 0.23(1), g = 2.31, and g ⊥ = 2.14, consistent with that from INS experiments.In addition, first-principles calculations present electron density distributions and therefore visualized superexchange paths, thus providing direct and indubitable proof on the spin-chain alignment in material CN.Furthermore, through TTN simulations, we show that CN has large entropy change and pronounced peaks (and dips) in Grüneisen parameter around QCPs at low temperatures, and the calculated adiabatic temperature changes can fit strikingly well to previously measured isentropes, revealing that CN may be a remarkably ideal quantum critical refrigerant.
The rest of the article is arranged as followings: Section II presents a brief overview of the previous experimental and theoretical studies on material CN.An ab initio study is also carried out which reveal intra-and inter-chain couplings directly and explicitly; Section III is devoted to a review of the TTN algorithms and the development of bilayer LTRG; In Section IV we show the TTN simulations and its fitting to experiments data, including both taken from previous measurements and done in the present work; Section V discusses the criticality-enhanced MCE in material CN and its potential usage as low-temperature coolant; Lastly, Section VI presents the conclusion and discusses some possible extensions of the present work.

II. ALTERNATING HEISENBERG ANTIFERROMAGNETIC SPIN CHAIN MATERIAL COPPER NITRATE
Copper nitrate hemipentahydrate [Cu(NO 3 ) 2 • 2.5H 2 O] was the first inorganic S=1/2 spin chain material studied experimentally [7,8].As one of the common copper salts, CN possesses some special thermodynamic properties at the low temperatures, including the zero-magnetization plateau [11,20], 1D Luttinger liquid behavior under magnetic fields [28,55], and 3D magnetic transition at ultra low temperature (150 ∼ 160 mK) [16,19,28], etc, which has aroused people's research interest for more than half a century, significantly promoting developments of the research on low-dimensional quantum magnets.
Before diving into detailed discussions on CN, we point that there exist other copper nitrate hydrates.Among others, the Cu(NO 3 ) 2 • 3H 2 O (copper nitrate trihydrate) has similar formulas and exactly the same X-ray diffraction pattern [56], and thus is easy to be confused with CN hemipentahydrate studied here.We paid special attention to discriminate between these two similar hydrates, and according to thermogravimetric analysis [24,26], it is concluded that Cu(NO A. Crystal Structure Figure 1 depicts the crystallographic structure of CN, which is monoclinic with space group I12/c1 [9].The corresponding lattice constants are found to be a = 16.1 Å, b = 4.9 Å, c =15.8 Å, and β = 92.9• at low temperatures (∼ 3 K) [23,25].As shown in Fig. 1

B. Alternating Heisenberg Antiferromagnetic Chain
Based on early experimental research on CN, including the measurements of magnetic susceptibility [7] and specific heat [8], a binary cluster model for describing its magnetic properties was proposed: Cu 2+ ions having spin S = 1/2 are coupled in pair (with coupling strength 1  2 J/k B = 2.56 K) [7], and the system thus comprises of independent spin binary clusters.In addition, people also perform proton magnetic resonance (PMR) experiments and confirm that the short-range antiferromagnetic order are within spin-pairs, thus supporting the dimer model.[54] Although the binary cluster model can capture some of the main features of spin-spin correlation in CN, discrepancy between this simple model and experimental measurements still remains.A weak interdimer exchange interaction (J 2 ) was then brought into the binary cluster model in Ref. 11, which substantially improves the fitting to the isothermal magnetization curve.In addition, van Tol, et al., studied the magnetic phase transitions via PMR and uncovered a crossover to shortrange magnetic order at 350 mK, as well as a critical transition to long-range order at 160 mK, under a magnetic field of 3.6 T, which can be understood only by introducing inter-dimer couplings in the model [13].Besides the PMR data, this is also evidenced in the specific heat measurement, where the low temperature hump(peak) signals the crossover(transition) into short-(long-) range ordered state [15,18], also suggesting the existence of weak inter-dimer couplings.
Early attempts to introduce inter-dimer interactions include two possible model structures: a ladder model (with dimers on the rungs) and an alternating chain model [57].These two models both fit thermodynamic measurements of CN well since their corresponding thermal predictions are essentially equal [18,22].However, the discrimination between these two models was later done by the angular dependent PMR [19], according to which one identifies the possible superexchange paths and thus rules out the ladder model.In addition, neutron-diffraction results also support the spin chain scenario [21].Therefore, it has been concluded that an AHAFC model can very well describes the magnetic properties of CN (in the temperature regime above ∼ 160 mK), which reads where S = {S x , S y , S z } is the vector spin operator containing three spin operators in different directions; J = J 1 is the strongest superexchange coupling; α = J 2 /J 1 is the relative strength of dominant inter-dimer interaction, whose precise value was measurably different in various experiments and left undetermined between 0.24 and 0.27 [22,23].Also note that in the magnetic-field coupling (Zeeman) term, the Landé factors are different (g = g ⊥ ) on the direction along b axis and that perpendicular to it.This magnetic anisotropy has been observed experimentally in the magnetic susceptibility measurements for a period of time [7], while little understanding has been achieved yet.
Figures 1(b-e) depict the spin chain structure and the spinspin interaction paths.The distances between one Cu 2+ ion to three neighbors are 5.33 Å , 6.22 Å , and 6.32 Å [18], which leads to, through Cu-O-H-O-N-O-Cu bridges, three distinct couplings J 1 , J 2 , and J 3 , respectively [Fig.1(b)].Therefore two possible inter-dimer superexchange paths J 2 and J 3 are shown in Fig. 1(d), and the spin chain thus could have had two possible routes on any (10 1) plane [see Fig1(d)].Until recently, inelastic neutron scattering (INS) determines that J 3 =-0.01meV (of magnitude about 1/10 of J 2 ) [27], so that J 2 is the dominant inter-dimer interaction, which connects dimers to form a tilted alternating chain, as shown n Figs.1(c-e).
Moreover, from Fig. 1(e), we can see that there exist four inequivalent types of ( 10

C. First-Principles Calculation and Electron Density Distributions
In Subsection II B, we scrape together quite a number of experimental observations, including various thermodynamic measurements and INS results, and arrive at the conclusion of an AHAFC model description for CN.However, a thorough study of electronic structures in CN via ab initio calculations is indispensable, which may provide an direct check for the existence of spin-chain type magnetic interactions in CN and offer insight into exchange path other than intra-chain couplings.
In this work, we employ a self-consistent field calculation, based on the all-electron projector augmented wave (PAW) method [58] implemented in VASP [59], to investigate the electron density distributions in CN.We adopt the generalized gradient approximation of Perdew, Burke, and Ernzerhof exchange-correlation functional [60].The cutoff energy for the plane wave expansion is chosen as 1000 eV, and the kpoint mesh is 2×3×2.In practical calculations, little changes both in the cell shape and atomic positions have been obseved after structure relaxation, hence the experimental lattice parameters shown in Fig. 1(a) are used, and two unit cells which comprise 264 atoms (including 16 copper atoms) are selected.
In Fig. 2, we show the simulated results of electron density distributions.Remarkably, in Figs.2(a,b) the spin chain alignment in (10 1) plane is clearly demonstrated, where the electrons tend to reside along the chain directions and thus leads to larger exchange integrals J 1 and J 2 [see Fig. 1(a)].Note that from the calculated results, we can discriminate J 2 from J 3 without any ambiguity, where the Fig. 2 show that the electron densities (hence also the coupling strength) are different in J 2 and J 3 bonds for orders of magnitudes.This difference between J 2 and J 3 , as well as the fact that tilted chains are along difference directions between I, III and II, IV planes [Figs.2(a,b)], agree with the INS observations in Refs.23, 27.Moreover, in Fig. 2(c) we show the electron densities in (010) plane, where the J 1 dimers are highlighted, from which we can see that there exists a weak dimer-dimer exchange coupling J m between every pair of dimers along [001] direction, this again has been observed experimentally [23].

D. Inter-chain Couplings and 3D Heisenberg spin model
The AHAFC model can explain well most experimental observations of CN in a very wide range of temperatures (above a few hundreds of milikelvins).Nevertheless, ultra low temperature measurements of PMR [13], adiabatic susceptibility [17], and heat capacity [22], uncovered a 3D phase transition at about 150 ∼ 160 mK, under magnetic fields between B c and B s .The existence of such a finite-temperature transition from magnetic disordered phase to 3D long-range order phase [28] is clearly beyond the 1D spin model given in Eq. ( 1).
In Subsection II B, we mentioned that the ferromagnetic coupling J 3 whose magnitude is about 1/10 of J 2 [Figs.1(b,d)], and thus plays the role of a weak inter-chain interaction.In addition, INS experiments reveal that there exist some other inter-chain interactions, namely, J m = 0.018 (2) meV, between dimers and along [001] directions [Fig 1(e)] [23].These important facts have been confirmed by our firstprinciples calculations in Subsection II C, as Fig. 2 illustrates.
Put all these intra-and inter-chain interactions together, we arrive at a 3D Heisenberg spin model with Hamiltonian H 3D = H (10 1) + H CTC , which reads: and Equation ( 2) is a 2D honeycomb lattice model with three different coupling constants along î, ĵ, and k directions [see Fig. 1(d)].p is the coordination vector in {A} sublattice [Eq.2] or running over all lattice sites [Eq.3] of honeycomb lattice on (10 1) planes.One typical site in {A} sublattice has been marked in Fig. 1(d).Note that J j = J 2 and J k = J 3 on planes I and III, while J k = J 2 and J j = J 3 on planes II and IV.
Equation ( 3) represent the inter-layer couplings, J m is the dimer coupling, p is coordination vector of every localized spin site, and c has magnitude of lattice constant c = 15.8Å and is along [001] direction (there must be a spin site at p + c/2 according to CN structure shown in Fig. 1).Some additional remarks on 3D inter-chain couplings in Eq. ( 3) are in order: Since there exist four kinds of 2D honeycomb planes [I to IV in Figs.1(c,e)] which are different from each other by lattice shifts and reflection operations about (010) plane (see related descriptions in Section II B), and the chains are arranged along different paths ([111] or [11 1]) on different planes, Eq. 3 thus describes a 3D coupled tilted chains (CTC) model.In this model, two spins are coupled via J 1 (a dimer) when both involved sites are of equal height in b axis and have a shift in distance of c/2 7.9 Å in [001] direction in ac-plane.Therefore, the seemingly simple H CTC in Eq. ( 3) actually represents a quite peculiar 3D inter-chain coupling model which looks bizarre while are actually feasible in the material CN.This 3D model has not been reported before as far as we know, neither has its properties been explored.In Ref. 23, INS experiments also show there exists interchain interaction between nearest dimers along [100] direction.However, we find that by shifting dimers along [1/2 0 0] as indicated by the authors in Ref. 23, there locates no dimer in the supposed position (see Fig. 1).This is also verified in our Abinitio calculations, where Fig. 2(c) shows clearly that there is no visible dimer-dimer coupling between a dimer and its nearest neighbor along [100] direction.Therefore, we include only the inter-dimer coupling J m along [001] direction in Eq. ( 3), while leaving it as an open problem about the possibility of adding more inter-chain coupling terms to Eq. (3).Note that the inter-chain interactions are rather weak and does not alter the physical properties except at ultra low temperatures.In the followings, except for the discussions on criticality-enhance MCE in CN, this 3D model Eqs.(2,3) will not be involved, and we focus on the AHAFC model description in Eq. ( 1) exclusively.

III. THERMAL TENSOR NETWORK APPROACH
High-precision thermal quantum manybody calculations are indispensable for relating the spin models discussed in Sec.II to the thermodynamical pmeasurements of CN.In this section, we firstly review the LTRG method proposed by some of the authors in Ref. 52.Then, we promote LTRG to a double-layer form with significantly improved accuracy, and then discuss its intrinsic relation to the well-established TMRG approach.

A. Thermal Tensor Networks and Linearized Tensor Renormalization Group Method
Employing the Trotter-Suzuki decomposition [53], we obtain a 2D TTN for the 1D Heisenber chain, which consists of rank-four tensors In order to calculate the thermodynamic properties one needs to accurately contract the TTN, which, however, is a NP-hard problem and thus can not be solved exactly.Therefore, people has to resort to approximate methods for efficient contractions of TTN.Among others, renormalization group algorithms constitute an important class of approaches which are developed to accurately contract the TTN and calculate interested quantities including free energy per site, specific heat, magnetic susceptibility, entropy, and others.
As shown in Fig. 3, the interconnected local tensors constitute a 2D checkerboard-style TTN, which subsequentlyp can be regarded as repeated 1D vertical (T 1 , T 2 ) or horizontal (V 1 , V 2 ) stripes, i.e., transfer matrices, as shown in Figs. 3 (b,c).The full contraction of TTN and consequently the calculations of thermal properties can be accomplished with the help of these transfer matrices.For instance, the vertical stripes T 1 , T 2 transfers the spin indices {σ i } between different lattice sites and the partition function Z thus reads where L denotes the total length of the chain.In the thermodynamic limit the dominant eigenvalue λ max of transfer matrix T = T 1 T 2 determines the free energy per site f = 1 2β ln λ max and also other thermodynamic quantities.
In order to calculate the extreme eigenvalue (and corresponding eigenvector), in Ref. 50, Xiang and Wang utilized the density-matrix renormalization group (DMRG) method, originally developed for Hamiltonian system, to solve the transfer matrix problem.They perform RG process along the Trotter direction and truncate the accumulated {σ i } indices/states into a fixed number (M ) of renormalized states.TMRG can determine the thermodynamic properties with high precision, and has been established as the method of reference in calculating thermodynamics of 1D strongly correlated quantum lattice systems [61].
Alternatively, the efficient contraction of the 1+1D TTN can also be performed by making use of the horizontal transfer matrices V 1 , V 2 [see Fig. 3 (c)].Li et al. proposed a TTN algorithm dubbed as LTRG [52], which project continually the transfer matrix V 1(2) to the density matrix of the system (in a form of matrix product operator, MPO).LTRG method can be used to accurately calculate the thermodynamics in 1D chains [52,62] and applies also to higher dimensional lattices [52,63].At a first glance, it seems that these two RG methods, TMRG and LTRG, are quite different, since they are dealing with transfer matrices along two distinct axes, i.e., the horizontal and vertical directions, which are clearly inequivalent: the spatial direction is infinite while the thermometric axis is finite and subject to a periodic boundary condition.However, a closer look into it below reveals that they are actually intimately related to each other, on the ground that both methods manage to accurately contract the TTN, while such globally optimized contraction schemes should treat both directions (implicitly) in equal footing.

B. LTRG and LTRG++ Algorithms
In this subsection, we firstly briefly review the single-layer LTRG proposed in Ref. 52, and then propose its bilayer form (Fig. 4), which improves the accuracy prominently and is dubbed as LTRG++.For more technical details on LTRG, we refer to Appendix A.
In the LTRG algorithm, only a single MPO [upper or lower one in Fig. 4(a)] is involved in the process of contractions: Starting from the single MPO at initial inverse temperature τ , and by continually projecting ν tensors onto upper (or lower) MPO, we cool down the system from very high temperatures (1/τ ) to various lower temperatures, β = 1/((n+1)τ ) at n-th step and thus the thermodynamics at temperature 1/β can be calculated (see more details of single-layer LTRG in Ref. 52).LTRG can produce quite accurate results for thermodynamic properties even at quite low temperatures, and the precision is comparable to that of TMRG.[52] Nevertheless, besides the single-layer algorithm, we devise here a double-layer LTRG++ algorithm, as shown in Fig. 4, for contracting the TTN.The main idea is as following: we contract the TTN into two (instead of one) MPOs, this is what "++" means.Figure 4(a) exploits a symmetric construction, where the upper and lower layers are contracted into two MPOs in exactly the same manner, saving one half of the projection time.Due to the checkerboard structure of TTN, each projection step comprises of two substeps, called as the odd and even substeps.After n steps, one reaches the inverse temperature β = (2n + 2)τ .The free energy per site f (β) = lim N →∞ − 1 N log(Z) can now be calculated from the series of renormalization factors κ a and κ b , extracted at odd and even substeps, respectively [64].In addition, we also need to calculate the dominant eigenvalue λ max of transfer matrix consisted of T a , T b and their conjugates in Fig. 4(b), with corresponding left (right) eigenvectors E l (E r ).With these results, the free energy per site f at β = 2(n + 1)τ can be computed via The advantages of LTRG++ over the previous single-layer algorithm can be summarized in mainly two aspects: Firstly, LTRG++ further improves the accuracy of LTRG, and is now practically of the same precision as TMRG; Secondly, LTRG++ can save half of the projection time.Related discussions and numerical evidence can be found in Appendix A. Note that in the LTRG++ algorithm one can also adopt an asymmetric construction such that two MPOs are (approximately) conjugate to each other, and recovers then the purification scheme [65] used in finite-temperature DMRG [66].Therefore, TTN approaches are very flexible, and can be designed to rewrite/recover the well-established TMRG, purification, and potentially other thermal RG methods.

IV. THERMAL TENSOR NETWORK SIMULATION AND FITTINGS TO CN EXPERIMENTAL RESULTS
Fitting the numerical simulated results to experimentally obtained magnetothermodynamic data is an important way to determine the parameters in Hamiltonian Eq. ( 1) describing the CN chain.In particular, in the present study, the purposes of performing numerical fittings are two-fold.
Firstly, we notice that fittings to the thermal properties includes the magnetic susceptibility, magnetization curves, and variousp specific heat curves have been performed in Ref. 22, where part of the data can be nicely fitted based on ED calculations on system length up to L = 12 (and their extrapolations).The authors get a set of fitting parameters, among which the coupling ratio between weak and strong bond is α = 0.27, and J/k B = 5.16(4)K) [see Hamiltonian Eq. ( 1)].However, the specific heat curves at large magnetic fields (2.82 and 3.57 T) were poorly fitted by the ED calculations.Actually, since at such strong field (2.8 T < B < 4.4 T) the ground state of the system is in a quantum critical region (Luttinger liquid phase [28]) and is supposed to have rather long correlation lengths at low temperatures.Therefore, ED calculations with such a small system size (up to L = 12) might be insufficient for an accurate estimation in that circumstance, and a large-scale thermodynamic algorithm is essentially needed to clarify this ambiguity.
Secondly, and more importantly, it is also noticed that the coupling ration α = 0.27 determined in Ref. 22 is "measurably different" from that from INS experiments, where it is determined that α 0.24 [with strong bond J = 0.442(2) meV = 5.13(2) K]. [23] The authors in Ref. 23 ascribed this discrepancy to the influence of interchain interactions in thermodynamic fittings.However, this argument can not be fully plausible since the interchain interaction is so weak (∼ 0.06K) [18,19] and could hardly induce any sizable effects on the thermodynamic properties in the temperature range of 2∼20 K so as to shift the fitting from α = 0.24 to 0.27.
Therefore, we perform state-of-the-art TTN simulations developed in Sec.III B and fit the magnetothermodynamic data both taken from Refs. 7, 8, 22 and measured in the present work.In our lab, we prepared various CN single-crystal specimens with sizes up to one or two centimeters (see Appendix .B for details of sample preparation) and measured their magnetization in high-precision SQUID devices.In practical calculations, Trotter slice is set as τ = 0.025, the lowest temperature reached is T /J = 1/150 (i.e., inverse temperature β = 150), and χ = 400 ∼ 600 bond states are retained, with truncation error smaller than 10 −13 .Numerical convergence versus χ of various concerned quantities including free energy, specific heat, magnetization curve, etc, has always been checked.

A. Specific Heat at Various Magnetic Fields
We start from the specific heat curves at various magnetic fields B = 0, 0.87, 2.82, and 3.57 T, as shown in Fig. 5. Ex-perimental data (symbols) are taken from Ref. 22, and the coupling constant J, chosen to be 5.13 K, is within the fiducial range of 5.16(4) K from previous thermal fitting and 5.13(2) K from scattering fitting.Actually we find out that small change of J (say, ±0.3 K) does not cause significant changes calculated results of quantities so as to affect other fitting parameters.
In Figs.5(a), 5(b), we plot specific heat curves of low magnetic fields (B = 0, 0.78 T), and both fittings with α = 0.23 (solid lines) and 0.27 (broken lines) are displayed as comparisons.It is seen clearly that the calculated curves of both α values can fit the magnetic specific heat curves almost equally well, for either B = 0 or B = 0.78 T case.Therefore, it is difficult making a preferable choice amongst these two α values, as well as potentially many other values in between.
In Fig. 5(c) the measured specific heat curve C p shows double peak structure, and α = 0.24 and 0.27 curves start to show some qualitatively different behaviors: While the α = 0.27 curve only presents a shoulder below 1 K, the calculated α = 0.23 curve correctly captures the double-peak structure, making the latter fitting noticeably better than the former.Moreover, difference between two fittings becomes more strikingly in Fig. 5(d) [67], where C p in the region 0.3 ∼ 1.5 K is quite sensitive to the change of α, and α = 0.23 is obviously superior than 0.27 in this case.
Therefore, from the direct comparisons in Fig. 5, we conclude that α = 0.23 is an overall better parameter than α = 0.27 in the fittings of specific heat curves at various fields.We would like to stress that the discrimination between α = 0.23 and 0.27 can be done only if an accurately calculation is possible for the low-temperature thermodynamic property of CN at high fields (2.82 and 3.57 T) where the ground state is a critical Luttinger liquid.The authors in Ref. 22 performed fittings to these data based ED simulations on rather limited system size (up to 12 sites) and get poor agreements, thus did not manage to discriminate α = 0.24 from α = 0.27 [68].

B. Magnetic Susceptibility and Magnetization Curve
In Sec.IV A, we find the parameter α = 0.23 prevails α = 0.27 for fitting specific heat curves under magnetic fields.In this subsection, we check whether the preferred parameter α = 0.23 can also fit other magnetothermodynamic quantities such as zero-field magnetic susceptibility measurements χ and the magnetization curves at various temperatures.
Figure 6 illustrates the fittings to magnetic susceptibility reuslts, which comprises data measured in the present work and those taken from Ref. 22.In particular, the present magnetic susceptibility measurements are performed in order to fill up the gap in the temperature range 5 K < T < 15 K where the old susceptibility data are absent.It is seen that in Fig. 6 the TTN calculations can fit the experimental results very well.Note that the magnetic susceptibilities are measured both along and perpendicular to the crystal b axis, Fig. 6 reveals that there exists quite prominent annisotropy in the spin chain material.It turns out, through the fittings of both susceptibilities with Hamiltonian Eq. 1, that this annisotropy can be attributed to different Laudé factors in the directions parallel (g = 2.31) and perpendicular (g ⊥ = 2.14) to the b axis.
Besides the zero-field χ, we also fitted the magnetization curves at various temperatures (517 mK, 2, 2.03, and 5 K).In Fig. 7 the measured magnetization curves with magnetic fields perpendicular to b axis, parallel magnetization curves taken from Ref. 22, as well as 517 mK data from Ref. 28, are quantitatively fitted with the set of parameters J = 5.13 K, α = 0.23, g = 2.31, and g ⊥ = 2.14.
To summarize, through above high-precision fittings, we conclude that the coupling ratio α determined is close to that (α = 0.24) obtained from INS experiments [23], while "measurably" different from α = 0.27 in previous thermodynamic fittings [22].This finding reveals that the thermal and scattering experiments are actually consistent with each other, and the previously supposed discrepancy may be due to limited simulations in fitting low-T thermal data of gapless Luttinger liquid phase.

V. CRITICALITY-ENHANCED MAGNETOCALORIC EFFECT
As early as 1968, magnetic refrigeration in spin chain material Cu(NO 3 ) 2 • 2.5H 2 O has been experimentally observed, where a temperature decrease ranging from 1 to 58 mdeg was measured by increasing fields form 0.46 to 0.96 T [8].Later on, van Tol et al. [15] explored the isentropic lines with fields in the range from 2 to 5 T, and observed two prominent dips with temperature as low as ∼ 100 mK (near two transition fields B c and B s ).This phenomenon was later attributed to large entropies near two fields, through ED calculations and theoretical analysis [20].However, it was unclear then that this enhanced MCE is due to quantum criticality at two transition fields, and a close comparison to experimental measured isentropes was absent, perhaps due to the lack of powerful manybody computation tools.
In this section, we employ the TTN simulations to explore the isentropes and magnetic Grüneisen parameters, and revisit the early isentropic data in Ref. 20.We reveal that there exists criticality-enhanced MCE near two field-induced QCPs.In Fig. 7 the calculated magnetization at T =40 mK [69], where two QCPs, i.e., the plateau-closing field B c = 2.87 T and saturation fields B s = 4.08 T, are clearly shown.Between these two QCPs, there exists a continuous critical Luttinger liquid phase which hosts gapless magnetic excitations.
In Fig. 8(a), we plot the isentropic curves of various magnetic entropies (from S/R = 0.05 to 0.5).For curves with relatively large entropies (0.2 ≤ S/R ≤ 0.48), the lowest temperature appears at around B 3.5 T, roughly located in the center of gapless region.However, with further lowering temperatures, we see that the broad dip eventually splits into two sharper dips in the isentropic curves, signalling two QCPs.Fig. 8(a) therefore manifests that in the both vicinities of two QCPs and in the quantum critical region, the thermal entropies are relatively large, which in turn results in criticalityenhanced MCEp.
Along each isentropic curve, one can read out the adiabatic temperature changes with varying magnetic fields.A quite distinct future of Fig. 8(a) for CN chain is that on both small and large field sides, one experiences large temperature changes by varying fields.This is in contrast to uniform Heisenberg model (see, for instance, Fig. 2 in Ref. 49 for spin chain material CuP), where only on large field (right) side of saturation QCP.For CuP chains, one can observe large MCE by decreasing from very large field to saturation point, while little temperature change is seen by increasing fields from 0 to saturation due to the presence of Luttinger liquid all along the magnetization curve.On the contrary, for the CN chain, the situation is different due to the existence of dimerization, which opens up a gap at low fields < B c .This fact enables us to realize criticality-enhanced MCE for relatively small fields (< 4 T), and one could even properly design a thermal cycling to make use enhance MCE around both low and high critical fields in one complete cooling process, with presumably much higher refrigeration efficiency than, say, uniform CuP chain [49].
In Fig. temperature changes evidences that CN indeed has criticalityenhanced MCE characterized by large temperature change even for moderate fields (say, from 0 to 3.5 T).Another important quantity measuring MCE property is the magnetic Grüneisen parameter Γ B = 1 T ( ∂T ∂B ) s , which is a differential characterization on the temperature change ∆T over small magnetic field variation ∆B in an adiabatic process.In the vicinity of QCPs, Γ B diverges as T tends to zero, whose scaling behavior is intimately related to the quantum criticality [43,44].In Fig. 9, we show the calculated Γ B of CN, and also the measured Γ B of uniform spin-1/2 Heisenberg chain material CuP as a comparison (taken from Ref. 49), from which it is seen that the CN chain has much larger Γ B around either one of its two QCPs, twice as large as that of CuP around the saturation field.The latter has been proposed as a perfect alternative for ordinary demagnetization refrigerant due to its wide operating range, large cooling power, and high efficiency [49].Our simulations here shows that the dimerized spin chain CN studied in the present work has even more promising potential as quantum critical coolant, not only because it has two sharp dips at suitable fields (Fig. 8), one at B c = 2.87 and one at B s = 4.08 T, but also due to large temperature change in response to the same field variation as revealed by calculated Γ B shown in Fig. 9.

VI. CONCLUSION AND OUTLOOK
In this paper, we generalize the linearized tensor renormalization group (LTRG) method to a bilayer form which is essentially equivalent to the well-established TMRG method, and employ this cutting-edge TTN method to accurately study the thermodynamics of a 1D dimerized spin chain material copper nitrate.We calculate and fit the experimental data of specific heat, magnetic susceptibility, and magnetization curves, some of which are measured experimentally in the present work.The previous discrepancy in coupling constants determined from different experiments are resolved, and we conclude that the set of parameters J = 5.13 K, α = 0.23(1), g = 2.31, and g ⊥ = 2.14 yielded from fitting thermal properties, is actually in remarkable consistency with those determined from INS experiments.In addition, based on electron density distribution pattern, we have for the first time visualized the spin-chain exchange path in CN, through ab initio calculations.On the grounds of these calculations, and also according to INS experiments, we obtain a 3D spin model describing magnetic properties of CN below about 160 mK.
With this set of parameters characterizing the magnetic properties of spin chain model material CN, we uncover, though accurate TTN simulations, that there exists criticalityenhanced large MCE near two quantum phase transition points, even at very low temperatures.Judged from the quantum anomaly in low-temperature isentropic curves and the extraordinarily good agreements to experimentally measured isentropes, as well as the large peaks/dips in magnetic Grüneisen parameters, we propose that CN is a very promising quantum critical coolant with significant temperature changes in response to magnetic field variations of moderate values.
There are still a number of interesting questions deserving further discussions, on both experimental and theoretical sides.To name a few, the direct experimental measurement of adiabatic temperature change for larger field ranges, instead of the rather limited field range between 2 to 4 T in previous experiments, is important to verify our prediction of CN as a promising quantum critical coolant.In addition, the performance characteristics such as operation temperature range, cooling power, and efficiency, are also in due to be investigated.Another important ingredient missing in the present work is the effect of inter chain couplings, whose Hamiltonian has been given in Eqs.(2,3).The inter-chain couplings could be of importance since the coolant is supposed to work in a circumstance with lowest temperature T < 100 mK, an energy scale comparable to that of inter-chain exchange constant in CN.In this appendix, we provide more technical details of linearized tensor renormalization group (LTRG) and LTRG++ algorithms and relate them to the well-established transfermatrix renormalization group (TMRG) approach.In Fig. A1(a), we firstly specify the local two-site Hamiltonian as h i,j = J i,j S i • S j between spin-1/2's.Such local spin-spin coupling term can be seen in the AHAFC model Eq.(1) of main text.Given the local tensor ν σ1,σ2,σ 1 ,σ 2 = e −τ hi,j as in Figs.4(b), one decomposes ν into P a(b) and then recombine them in a way [Figs.4(d,e)] to form the rank-four tensor U (2)   U (1)   U (2)   Λ scheme shown in Fig. A1(f).In either LTRG or LTRG++, we absorb the evolving operator ν to local tensors T a , T b in the MPO, as well as two bond weights Λ, and arrive at the M tensor of rank six.Then one reshape M into a matrix and decompose it to obtain two updated tensors Ta , Tb , and new bond weight Λ (find more details in Ref. 52).Note that since the so obtained MPO is quasi-canonical due to that every single evolving operator ν is close to identity.In order to calculate the free energy per site at any specific temperature, following Eq.( 5) in the main text, we need to collect the normalization factors κ a(b) i (obtained by normalizing Λ after each projection process) and the dominant eigenvalue θ max of ladder transfer matrix shown in Fig. 4(f) of main text.Notice that a proper combination of κ i 's and λ max , i.e., θ max = (Π n i=1 κ a i κ b i ) • λ max , is nothing but an accurate estimate of dominant eigenvalue of vertical transfer matrix [see Fig. 3(b) of main text] one pursues in the TMRG algorithm.
To be specific, we compare in Fig. A2 the LTRG++ algorithm with traditional TMRG.As shown in Fig. A2(a), TMRG exploits the "s-µ-e-σ" construction, adding two blocks (µ and σ) per iteration, and truncate enlarged system (environment) block s-µ (e-σ) with the aid of dominant eigenvectors of transfer matrix.Note that the transfer matrix is non-Hermitian and have two sets (left and right) of eigenvectors.Take the dominant left eigenvector as an example: one can construct the reduced density matrix of s-µ (or e-σ), and retain the largest m states (White's rule) to reduce the enlarged direct-product space s-µ (e-σ), executing the truncation of left indices of the transfer matrix.The procedure is similar for the truncation of right indices, with the help of dominant right eigenvector.
At first sight, this algorithm, firstly proposed in the language of traditional DMRG [50], is quite different from the bilayer LTRG++ algorithm introduced in section III B of main text.The latter is developed in the language of tensor networks, and it contracts the tensor networks along Trotter direction, seemingly different from the former.However, in fact, these two algorithms are, despite some minor technical details, essentially equivalent.In Figs.A2(b,c), we show that there exists a hidden matrix product state (MPS) in LTRG++ algorithm, which is nothing but the dominant eigenvector of vertical transfer matrix in TMRG.In Fig. A2(b), take the left eigenvector E l as an example, it is shown that E l U Λ 2 U † where Λ is a diagonal matrix storing the Schmidt coefficients, obtained from singular value decomposition on the specific bond [52].Given that, we can further expand E l in terms of the series of U (i) and (U (i) ) † , the truncation matrices generated in each projection step (not stored during the process of computations, though), and arrive at Fig. A2(c) which reveals explicitly the hidden MPS representation of dominant eigenvector of transfer matrix.
To confirm the above argument, in Fig. A3 we compare numerically the accuracy of single-layer LTRG, TMRG and bilayer LTRG++ results for the exactly soluble XY spin chain model.Note that there exists a point where the truncation error happens to cancel the Trotterp error completely, leading to a very steep dip in Fig. A3: For small β (high temperatures regime) Trotter error is dominant, while for large β (low temperature regime) truncation error takes the place.Since the Trotter error accumulates linearly with β and is independent of the specific truncation schemes, the position of the cancellation point thus reveals how fast truncation errors accumulate and plays the role as an indicator of "goodness" of the truncation scheme.In principle, the further the cancellation point dwells on β axis, the better the truncation scheme is.Therefore, according to Fig. A3, bilayer LTRG++ has significant improvement compared to the single-layer LTRG, and bears practically the same accuracy as TMRG.This observation strongly supports our argument above that LTRG++ is in essence equivalent to TMRG.
To summarize, although it seems that LTRG++ and TMRG are following reversed orders in contracting the TTN (Trotter or spatial direction first), a more careful analysis, though, reveals that the truncation matrices in LTRG++ constitute a hidden MPS and these two methods are essentially equiva-lent.This can be attributed to the fact that a full contraction of TTN must consider both directions in equal footing and the superficial order of contractions does not really matter.In fact, these two renormalization group (RG) algorithms, and potentially other RG based thermal algorithms, can be unified in the framework of TTN simulations, more details on this topic will be published elsewhere.
Figure1depicts the crystallographic structure of CN, which is monoclinic with space group I12/c1[9].The corresponding lattice constants are found to be a = 16.1 Å, b = 4.9 Å, c =15.8 Å, and β = 92.9• at low temperatures (∼ 3 K)[23,25].As shown in Fig.1(a), a conventional unit cell comprises of eight formula units.Each Cu 2+ ion is surrounded by five nearest-neighboring oxygen atoms, which constitutes a distorted pyramidal polyhedron [Fig 1(b)].Four oxygen atoms resides at the vertices of the basal plaquette of the polyhedron, roughly on the same plane: two of these oxygen atoms belong to H 2 O molecules and the other two are from NO − 3 groups; the rest apical oxygen atom belongs to a third NO − 3 group.The pyramidal polyhedrons of opposite orientations are arranged alternatively along a line [Fig.1(b)].

FIG. 1 :
FIG. 1: (Color online) Crystal structure and magnetic exchange couplings in Cu(NO3)2 • 2.5H2O.(a) The unit cell of Cu(NO3)2 • 2.5H2O, where the coordinate axes coincide with the crystal axes.The lattice constants are shown in the figure, indicating that CN belongs to the monoclinic system.(b) The expanded view of the quare pyramidal polyhedrons comprise of oxygen atoms and Cu 2+ ions.(c) Superexchange paths between spins along chains in four inequivalent (10 1) planes which are adjacent to each other.(d) Structure in a typical (10 1) plane, where the Cu 2+ ions are highlighted while other atoms left transparent.The heavy solid lines are the intradimer interactions J1, the interdimer J2 and interchain J3 couplings are plotted differently (in black solid and red dashed lines, respectively).A labels one out of two sublattices of honeycomb structure in (10 1) plane, and î, ĵ, k are vectors connecting one site (in A sublattice) with its three nearest neighbors.(e) Projected view of the crystal structure in (010) plane, where the alternating solid lines represents the J1-J2 chain.We denote the four existing (10 1) planes as I II III and IV, respectively, where the chains have different paths in each plane.Arrows indicate the directions along [100] and [001], which represent interchain exchange path JL, Jm.

1 )
planes where the spin chains are arranged in different ways, namely, planes I to IV shown in Fig. 1(e).In I and III planes, the AHAFC stretch along [111] direction [from left top to right bottom, see Fig. 1(c)], with a shift of b/2 2.45 Å between chains in I and III planes; while in planes II and IV, the chains go from left bottom to right top ([1 11] direction), with the same shift of distance (b/2) between spin chains in both planes [Fig.1(c)].

FIG. 2 :
FIG. 2: (Color online) The projected electron densities on (a) III-type (10 1), (b) IV-type (10 1), and (c) (010) planes.a0 0.53 Å is the Bohr radius, the projection range of electron density is of thickness [-0.5, 0.5], respect to [10 1] unit vector for (a,b) and to [010] vector (i.e., primitive vector b) for (c) (refer to Fig. 1 for the specific crystal directions).The positions of copper ions are marked by solid balls.In (a,b) the tilted chain structures are clearly shown by high electron densities along the chain direction [111] for (a) and [1 11] for (b).In (c) the dimers with different hight on b axis are labeled in different colors, from which it is clear that there exist weak inter-dimer interactions (denoted as Jm) along [001] direction, while there exists no visible exchange path between two nearest neighboring dimers connected by 1/2 × [111] or [1 11] vector.

2 FIG. 3 :
FIG. 3: (Color online) (a) The 2D TTN represents the partition function Z of a 1D quantum lattice model, which exhibits a checkerboard pattern.(b) and (c) are the transfer matrices along spatial and Trotter directions, respectively.(d) depicts the rank-four local tensor ν, the elementary unit in the TTN.

FIG. 4 :
FIG. 4: (Color online) LTRG++ algorithm adopts a double-layer scheme which contracts simultaneously upper and lower layers into two MPOs.(b) T a(b) and its conjugate counterpart constitute the transfer matrix, with the left(right) dominant eigenvector E l(r) and the corresponding eigenvalue λmax.

FIG. 5 :
FIG. 5: Fitting to experimental data of specific heat curves under various magnetic fields, (a) B = 0, (b) 0.78, (c) 2.82, and (d) 3.57 T. The experimental data (symbols) are taken from Refs. 8, 22, and the dashed fitting lines are calculated with α = 0.27, while the solid lines are fittings with α = 0.23.The B = 0.78 T curve in (b) was measured with powder samples [22], thus is fitted using average Landé factor gav 2.2; and the B = 2.82 and 3.57 T curves in (c,d) are measured along crystal b axis, with Landé factor g = 2.31.

31 ⊥FIG. 7 :
FIG. 6: (Color online) Fittings to measured magnetic susceptibility χ taken from previous experiments (Ref.22), as well as those obtained in the present work (squares and circles).The latter is measured under a small magnetic field (B = 0.6 T) to mimic the zero-field susceptibility.χ has clear annisotropic g factors along the crystal b axis and the direction perpendicular to it.

FIG. 8 :
FIG. 8: Numerically simulated and experimentally measured isentropes of CN.(a) The contour lines represent entropy per site 0.08 ≤ S/R ≤ 0.48 (bottom to top) with interval ∆S/R = 0.04, where R 8.314 J • K −1 • mol −1 is the gas constant.(b) Comparisons to measured adiabatic isentropes of CN around the critical field Bs 2.87 T, the experimental data are taken from Refs. 15, 20.

FIG. 9 :
FIG.9:(Color online) Magnetic Grüeisen parameter ΓB (definition in main text), which characterizes differentially the temperature change over an unit magnetic field change.The lines plotted, with different hight of peaks, correspond to different ΓB at various temperatures, which decrease from 640 mK to 320 mK (top to bottom).

FIG
FIG. A1: (Color online) (a) A local coupling term Si • Sj, where S = {S x , S y , S z } is a vector spin operator.(b) Elementary tensor ν = exp (−τ h) can be decomposed into two rank-three tensors Pa and P b .(c) Matrix product operator representation of density matrix of spin chain.(d,e) Rank-four tensors T a(b) are obtained by recombining P b and Pa, which interconnects via the geometric (horizontal) bonds and form an MPO.(f) Local projection and truncation scheme.
T a,b 's are interconnected with each other via geometric bonds α, β and constitute an matrix product operator (MPO) of initial bond dimension D = 4 [Fig.A1(c)].After the initialization, we need to evolve the MPO along Trotter direction with the local projection and truncation FIG. A3: (Color online) Relative errors of free energy per site in LTRG, LTRG++, and TMRG, for the exactly soluble model 1D free fermion chain.M and D are the number of retained states during the process of projection and truncation in LTRG/LTRG++ and TMRG, respectively.
FIG. A4: (Color online) Examples of synthesized single-crystal CN specimens, which are transparent and with needle-like shapes.The needle axis is parallel to crystal b axis, and the lengths of samples range from a few mm to more than 1 cm.
U † is the dominant left eigenvector of the transfer matrix, where Λ is a diagonal matrix obtained in LTRG++.(c) Hidden matrix product state representation of E l in the LTRG++, revealed explicitly by expanding E l with {U (i) } series [and also (U (i) ) † ].