Large positive correlation between the effective electron mass and the multipolar fluctuation in the heavy-fermion metal Ce$_{1-x}$La$_x$B$_6$

For the last few decades, researchers have been intrigued by multipolar ordering phenomena while looking for the related quantum criticality in the heavy-fermion Kondo system Ce$_{1-x}$La$_{x}$B$_6$. However, critical phenomena induced by substitution level ($x$), temperature ($T$), and magnetic field ($B$) are poorly understood despite a large collection of experimental results is available. In this work, we present $T$-$B$, $x$-$T$, and $x$-$B$ phase diagrams of Ce$_{1-x}$La$_x$B$_6$ ($\mathbf{B}\parallel[110]$). These are completed by analyzing heat capacity, magnetocaloric effect (MCE), and elastic neutron scattering. A drastic increase of the Sommerfeld coefficient $\gamma_0$, which is estimated from the heat capacity down to 0.05 K, is observed with increasing $x$. The precise $T$-$B$ phase diagram which includes an unforeseen high-entropy region is drawn by analyzing the MCE for the first time in Ce$_{1-x}$La$_x$B$_6$. The $x$-$B$ phase diagram, which supports the existence of a QCP at $x>0.75$, is obtained by the same analysis. A detailed interpretation of phase diagrams strongly indicates positive correlation between the fluctuating multipoles and the effective electron mass.

For the last few decades, researchers have been intrigued by multipolar ordering phenomena while looking for the related quantum criticality in the heavy-fermion Kondo system Ce1−xLaxB6. However, critical phenomena induced by substitution level (x), temperature (T ), and magnetic field (B) are poorly understood despite a large collection of experimental results is available. In this work, we present T -B, x-T , and x-B phase diagrams of Ce1−xLaxB6 (B [110]). These are completed by analyzing heat capacity, magnetocaloric effect (MCE), and elastic neutron scattering. A drastic increase of the Sommerfeld coefficient γ0, which is estimated from the heat capacity down to 0.05 K, is observed with increasing x. The precise T -B phase diagram which includes an unforeseen high-entropy region is drawn by analyzing the MCE for the first time in Ce1−xLaxB6. The x-B phase diagram, which supports the existence of a QCP at x > 0.75, is obtained by the same analysis. A detailed interpretation of phase diagrams strongly indicates positive correlation between the fluctuating multipoles and the effective electron mass.

I. INTRONDUCTION
The sixfold degenerate wave functions of an isolated Ce 3+ ion are split into the Γ 8 quartet and the Γ 7 doublet under the crystal electric field with octahedral symmetry in Ce 1−x La x B 6 . The energy difference between the ground-state (GS) quartet and the excited doublet is about 46 meV (Ref. 1). Hence, low-temperature (T ) and low-field (B) properties are governed by the Γ 8 quartet and there are possibilities to observe multipolar moments including a dipole, a quadrupole, and an octupole (atomic 4f current distribution with zero magnetic moment). What has made this material so intriguing for the last few decades is that the above mentioned multipoles can interact through superexchange [2][3][4] . Indeed, antiferromagnetic (AFM) [5][6][7] , antiferroquadrupolar (AFQ) [8][9][10][11][12] , and antiferrooctupolar (AFO) orderings 8, [13][14][15][16] have been observed, and phase transitions between these ordered phases have also been intensively studied. By convention, the spin/pseudospin (orbital) paramagnetic (PM) phase, the AFQ phase, and the AFM phase are referred to as phases I, II, and III, respectively. The phase IV which appears around x ≃ 0.3 is suspected to represent AFO order 16 .
However, even though largely scattered information about the multipolar phase transitions is gathered 17 , the result is rather difficult to understand. The x-T phase diagram at B = 0 reveals a weakening of the phase IV as T is increased, while the x-B phase diagram at T =0 reveals an enhancement of the same phase as B is increased. Therefore, a discontinuity of the phase IV has been introduced along the x-axis in the up to date x-T -B phase diagram of Ce 1−x La x B 6 , and the criticality at the verge of the phase IV is poorly understood.
In this article, we investigate physical properties of Ce 1−x La x B 6 (x=0, 0.18, 0.23, 0.28, 0.5, 0.75) by measuring heat capacity at ambient pressure, magnetic Bragg intensity, and the magnetocaloric effect (MCE) with B [110]. C p and elastic neutron scattering (NS) measurements are analyzed to construct x-dependent T -B phase diagrams and the x-T phase diagram which show intricate transitions between multipolar phases. From the low T analysis of the C p /T , the Sommerfeld coefficient, γ 0 , is extracted. The B-dependent γ 0 exhibits local maxima at critical points. At low B, the overall magnitude of γ 0 drastically increases for x > 0.28. Based on the unprecedented analysis of the MCE in this system, a novel high-entropy region in the x-B plane is revealed. This region expands as x is increased until it merges with the phase IV. Most importantly, we present an understandable x-B phase diagram, which shows weakening of the phase IV as x is increased and predicts the emergence of a putative quantum critical point (QCP) at the verge of the phase IV. The new observations concerning the variation of γ 0 on the x-B plane are explicable by a large positive correlation between multipolar fluctuations and the effective electron mass, m * , within heavy Fermi-liquid (FL) description.

II. RESULTS
In Fig. 1a, phase transitions in Ce 1−x La x B 6 are signaled by sharp anomalies in C p /T if the Ce concentration is high. However, the anomalies are weakened and the overall magnitude of C p (T )/T is diminished with increasing La concentration. To define the critical temper-  ature, it should be reminded that the mean-field (MF) theory predicts a discontinuous jump in C p /T at the transition temperature. We can extrapolate the discontinuous jump as noted by the red and blue lines drawn in Fig. 1b. Here, the so-called method of entropy balance is applied to define the Néel temperature T N and the AFQ transition temperature T Q : Note that hatched regions with the same color have the same area. On the other hand, a transition from a non-collinear AFM phase (III) to a collinear AFM phase (III ′ ) at T ′ N typically results in a kink in C p /T . Each region of multipolar phases and the corresponding transitions are confirmed from the literature 17,18 . Fig 1c shows that the entropy as a function of T loses its features for the phase transitions (see kinks in the black curve at T N = 2.4 K and T Q =3.3 K) as x is increased. Meanwhile, the overall magnitude of the entropy is increased with x. Due to the remaining short-range interactions between multipoles in phase I, we can observe saturated entropy of the randomly populated quartet only well above 10 K (Ref. 19). The field is applied along [110] crystallographic orientation because a brief x-T -B phase diagram is only available for this direction as a reference 17,18 . Figure 2 shows the field-induced Bragg intensity at the (1/2 1/2 1/2) wave vector in three samples of Ce 1−x La x B 6 measured by NS. Unlike the AFM phase [5][6][7]20 , which exhibits Bragg scattering in zero magnetic field, the hidden-order phase II can only be visualized by the elastic intensity that originates from the field-induced dipolar moments modulated by the underlying AFQ structure. Not only does this offer a method to reveal the presence of AFQ correlations by NS, but it may also provide information about the type of the multipoles involved 21 .
The field dependence in Fig. 2 reveals that the onset of intensity measured at (1/2 1/2 1/2) in all the studied samples occurs in several steps. In Ce 0.82 La 0.18 B 6 ( Fig. 2a), the intensity remains almost zero within the phase III, but starts to increase already at B = 1.5 T inside the phase III ′ for T = 0.07 K. This demonstrates that this distinct phase represents a combination of conventional (dipolar) AFM and AFQ order parameters, in contrast to the purely dipolar phase III. Upon entering the phase II, an even steeper jump in intensity is observed around 2.2 T. With increasing temperature, both transitions are shifted towards lower fields as the AFM phases are suppressed. A similar two-step increase without fieldhistory dependence is also seen in the low-temperature measurement on Ce 0.77 La 0.23 B 6 (Fig. 2b). In our third sample, Ce 0.72 La 0.28 B 6 , the zero-field GS is reportedly taken over by the phase IV (Ref. 18), but the phases III and III ′ can still be stabilized by the application of field. As a result, the corresponding field dependence in Fig. 2c now exhibits a series of three steps at T = 0.1 K: a broadened transition near 1 T from phase IV to III, a kink around 1.5 T corresponding to the III-III ′ transition, and the third transition to phase II at 2.5 T (Ref. 18). We note that for the La substitution level of x = 0.28, some initial increase in field-induced intensity is observed already within phase III, unlike at low La concentrations where no signal within phase III was found. This suggests that the field-stabilized phase III in La-substituted samples is different from the zero-field phase III with regard to its intermixing with the AFQ order parameter. On the other hand, within phase IV, i. e. for fields below 1 T, the (1/2 1/2 1/2) intensity is nearly insensitive to the field. Fig. 3. Phase boundaries and the degree of critical fluctuations are readily recognizable by the contrasting colors in the contour plots of C p (T, B)/T . The low T phase boundaries (yellow solid lines) are determined from the MCE analysis which will be explained with Fig. 4. The notations for the critical fields are straightforward, e.g. the critical field across the phases I and II is labelled by B I-II . In Figs. 3b-3d, the red upward-triangle and the red downward-triangle denote B II-III ′ and B III-III ′ determined by NS shown in Fig. 2. In phase II, the field-induced AFO order stabilizes the AFQ order resulting in positive slope of B I-II (T ) 13 . Phase II is gradually weakened and transformed into phase I with increasing x. For x = 0.23, a certain phase appears in the temperature range of 1 K < T < 1.6 K and fields below 1 T (Fig. 3c). When x = 0.28, this phase becomes the GS (Fig. 3d). Since the primary order parameter of this phase has not been conclusively determined, it is customarily called phase IV. The nature of phase IV will be discussed later. In Fig. 3e, broad peaks are almost insensitive to fields below 2 T. It will be argued that the phase IV survives in Ce 0.5 La 0.5 B 6 with a rather broad heat capacity anomaly. In Fig. 3f, Ce 0.25 La 0.75 B 6 shows field-insensitive peaks only in the low T and low B region. In the rest of the T -B plane, a specific-heat anomaly is broadened and moves to the high-T side as B is increased. This is reminiscent of the Kondo impurity behavior, but the detailed shape of C p (T, B)/T deviates from the resonance-level model 22,23 . Now, we present the quasi-adiabatic MCE observed in Ce 1−x La x B 6 . The ideal MCE is defined as reversible temperature changes of a thermally isolated specimen controlled by an external magnetic field. Suppose the S p,Bi (T ) curve is lying higher than the S p,B f (T ) curve. The sample tends to reduce its entropy by releasing a certain amount of heat as the field is changed from B i to B f . But the heat must be dissipated and the temperature has to be increased to keep the entropy unchanged. In the reverse procedure, cooling of the sample occurs. Likewise, a phase transition from a magnetically disordered state to an ordered state is described by a sharp increase of the sample temperature and vice versa. In practice, however, the sample is in a quasi-adiabatic condition in which the sample temperature, T qad , slowly and monotonically relaxes to T mix , the temperature of the mixing chamber of the dilution fridge to which the sample is mounted. Moreover, a field-induced heating around a critical region is found to be essential to describe the irreversible T qad (B) upon field cycling. We have constructed a mathematical model (see Appendix C for more details) for the above observations, and the basic principle is graphically summarized in Fig. 4. Eddy current heating is negligible in this experiment (Appendix C).
The background of Fig. 4a represents the color-coded entropy S p (T, B) of CeB 6 calculated from C p /T . Black (red) solid lines indicate T qad (B) when the field is increased (decreased). In Fig. 4b, representative magnetocaloric (MC) sweeps with T mix = 0.2 K and the sweep-rate r = 0.1 T min −1 are magnified. To simulate T qad (B), we assume a continuous phase transition, and start with a reversible ansatz for the true MCE. The distribution of the excessive field-induced heat, dQ(B)/dB, is also approximated in terms of an analytic function. Then, these factors are taken into the differential equation for the theoretical quasi-adiabatic temperature and consecutively adjusted until the solution is optimized to T qad (B). Fig. 4c is the amplified view of Fig. 4b at low B. Here, the reversible ansatz is the blue solid line. In Fig. 4d, the envelope of the black hatched-area (red-area) pertains to dQ(B)/dB in sweep-up (sweep-down) mode. The optimized solutions are denoted by dash-dotted lines in Fig. 4c. Note that the major cause of the irreversibility is indeed the field-induced heating accompanied by the maximum MCE. Magnified view of T qad (B) near the III-III ′ boundary is shown in Fig. 4e. An about two times larger change of T qad (B) in down-sweep is well reproduced by the hysteretic field-induced heating (Figs. 4e and 4f). The III-III ′ boundary is also explained by the model (not shown). We tried to simulate the total field-induced heat, Q, as closely as possible for both sweep directions, but no constraint regarding dQ/dB was imposed. These conditions reflect hysteretic behaviors across certain phase boundaries observed in magnetic susceptibility, magnetization, and NS 6,18,20 .
The yellow solid lines in Figs. 4a, 4g, and 4m connect critical points referenced to steepest slopes in the reversible ansatzes (see blue vertical dashed lines in other panels of the Fig. 4). On the other hand, a line which fades away above T = 0.4 K is inserted in Fig. 4a and stresses that the MCE is weakened as the temperature is increased. The shape of the reversible ansatz confirms that the magnetic entropy under this line is larger than the magnetic entropy of the well-defined AFM phase. While this unclosed borderline cannot be attributed to a phase transition in the strict sense, AFM domain selection or motion would also be insufficient to explain the origin of this line as these effects usually continue to T N (Ref. 20).  (Figs. 4i and 4j). The phase II-III ′ boundary is also analyzed (Figs. 4k and 4l). The unclosed borderline with the color gradient is shifted to higher fields compared with the one found in CeB 6 .
In Fig. 4m, phase IV becomes the GS. The model is applied to the representative T qad (B) shown in Fig. 4n. Analyses of III-IV, III-III ′ (not shown), and II-III ′ transitions repeatedly confirm that substantial Q is accompanied by the large MCE (Figs. 4o -4r). The III-IV boundary found by the MCE analysis is distinct, but we cannot draw such a well-defined phase boundary by observing a low T AFM structure since it is strongly field-historydependent 18 . Hence, the sudden onset of the AFQ phase with field-induced Bragg intensity (see Fig. 2c) might explain the sharp III-IV boundary.

III. DISCUSSION
Before we discuss our new findings, it has to be emphasized that the MCE analysis has played essential role to present new perspectives on thermodynamic phenomena in Ce 1−x La x B 6 . Experimentally, the completion of the precise B-T phase diagram is possible because T qad (B) is extremely sensitive to the entropy change, and this property allows us to explore relevant and the most interesting area in the B-T plane where anomalies in C p , NS, ultrasonic attenuation 24 , and thermal expansion 25 are not easily detectable. Theoretically, our phenomenological model implies substantial r-dependent quantum mechanical friction appears due to strongly fluctuating local magnetic moments, and suggest a non-equilibrium formulation for a field-driven phase transition. In terms of a quasiparticle dynamics, it is suspected that inelastic magnon-magnon and magnon-electron scatterings might be prominent. A detailed study about a r-dependent Q and a type of transition is now underway: A first-order phase transition is expected if both Q and dQ/dB are highly irreversible as in the III-IV transition (Fig. 4p), while a second-order transition is characterized by reversible Q and dQ/dB (Fig. 4l and 4q). A continuous III-III ′ transition is described by reversible Q and irreversible dQ/dB because of the domain hysteresis (Fig. 4f).
In CeB 6 , we observe S p = Rln2 exactly at T N (Fig. 3a). This indicates that, below T N , the two lowest-lying energy levels are approximately k B T N apart from each other while, above T N , the level scheme transforms into two separate Kramers doublets suitable for the AFQ ordering. In Ce 0.82 La 0.18 B 6 , the AFQ phase is weakened and the AFM phase is confined by a vertical phase boundary at low B (Fig. 3b). The isentropic line with S p = Rln2 follows this vertical boundary validating the above mentioned fundamental change in the energy level scheme.
The phase IV appears in Ce 0.77 La 0.23 B 6 and it is stabilized in Ce 0.72 La 0.28 B 6 with sharp C p anomalies indicating long-range ordering (Figs. 3c and 3d). The I-IV boundary is almost vertical and the isentropic line with S p = Rln2 is following this phase boundary. Whilst up to date resonant X-ray scattering (RXS) experiment [13][14][15] , NS 16 and MF theory 26 suggest the primary order parameter is of AFO type with the AFQ and the AFM moments as field-induced secondary orders, the vertical I-IV boundary may simply indicate that the lowest-lying levels in phase IV are barely affected by the Zeeman effect. In the preliminary quantum mechanical sense, even if an isolated CEF GS is the Γ 8 quartet, exchange interactions between 15 possible order parameters could further lift the degeneracy of wave functions at each Ce-site. In reality, at T I-IV , the level scheme seems to change into two doublets, and S p = Rln2 along the vertical section of the I-IV boundary originates from the lower doublet. We also suggest that, when major exchange terms are taken into account, an effective g-factor with respect to the lowest lying doublet is very small at low B.
In Ce 0.5 La 0.5 B 6 , sharp C p anomalies are absent, but the above proposed properties of phase IV seem to persist below 2 T (Fig. 3e and Fig. 6). The phase IV is easily suppressed and transformed into phase I in the region where B > 1 T and T > 0.25 K for Ce 0.25 La 0.75 B 6 (see Fig. 3f and Fig. 8). Fig. 5a shows that energy scales of AF-exchange interactions between multipolar moments are so close that most of the multipolar order parameters are competing in the limited region of the x-T space. The temperature difference between the point where the entropy of Rln2 is released and the point where the peak in C p /T is developed (hollow square) provides an estimate for the range of critical fluctuations. As noted by the hatched area, the MF-like C p anomalies are smoothed towards the putative QCP above x = 0.75.
In Fig. 5b, γ 0 (x, B) is extracted from C p (T, B)/T by estimating constant contribution at low T . For x = 0, the FL behavior is obvious and the constant value of C p (T, B)/T below 0.5 K (Fig. 1) is taken as γ 0 . As x is increased, the temperature range for the FL behavior diminishes, but we can easily extract γ 0 (x, B). Even though clear FL behavior is not seen down to the lowest temperature as in phase IV (x > 0.28), we can still extrapolate γ 0 because C p (T, B)/T shows linear T -dependence below 0.2 K. On the other hand, an applied field suppresses this T -linear dependence and a clear FL behavior is revealed. In this case, a part of the nuclear Schottky contribution (∝ 1/T 3 ) is included in C p (T, B)/T , and we simply subtract 1/T 3 term to get γ 0 . It turns out that γ 0x (B) has a local maximum at a critical point as marked by vertical dashed lines in Fig. 5b. Also, the magnitude of γ 0 (x, B) drastically increases with x inside phase IV. These tendencies have never been reported although resistivity measurements 27 support strong enhancement of m * in phase IV. In contrast to the previous report 27 , the low T specific heat capacity does not support the non-Fermi-liquid (NFL) ground state. Symbols in Fig. 5c denote critical fields determined by the MCE analysis with T mix = 0.02 K and the background is the contour plot of γ 0 (x, B). Representative T qad (B) curves with x = 0.5 and x = 0.75 are shown in the Figs. 7 and 8. Compared to the latest x-B phase diagram which is completed by collecting largely scattered experimental results 17 , the current x-B phase diagram is more reliable in three aspects. First, critical fields are well matched with the values of the peak positions in γ 0,x (B) (Fig. 5b). Second, as x is increased above 0.5, phase IV boundary shows negative slope towards the putative QCP. Thus, the discontinuity in the previously reported x-B phase diagram due to the positive slope of phase IV boundary at x > 0. 5 (Ref. 17) is now eliminated. Thirdly, a weak phase transition between phase I and phase II is suggested (see dashed line in upper-right area of Fig. 5c) which is consistent with the observation from Figs. 3e and 3f that phase I will occupy the region in phase space for composition x = 0.75 where phase II existed for x = 0.5 (see Fig. 3, Fig. 8 and Appendix D for more details).
Putting aforementioned observations together, the local maximum of γ 0 at the critical point, the increase of γ 0 with x, and the suppression of γ 0 with B strongly imply that multipolar fluctuations and m * in the heavy FL state of Ce 1−x La x B 6 have a large positive correlation.
As to the newly found high-entropy region below B × in Fig. 5c, we speculate that a thermally unstable (as described by the unclosed lines in Figs. 4a and 4g), but weakly correlated portion of heavy electrons is segregated from the system. Such a phase segregation would allow for more electronic degrees of freedom below the black dash-dotted line in Fig. 5c. Upon the mergence of the high-entropy regime with phase IV, it seems the positive correlation between fluctuating multipoles and m * is strongly enhanced.
In the near future, we hope experimental techniques such as scanning tunneling spectroscopy 29 , nonresonant inelastic X-ray scattering 30 and inelastic NS 31 are accessible for the further characterization of both the heavy itinerant electronic state and the multipolar state in Ce 1−x La x B 6 . Also, high-quality single crystals with 0.75 < x < 1 will allow us to scrutinize the existence of a putative QCP and its possible conjunction with the superconductivity which emerges as x approaches to unity 32 . To conclude, a set of temperature (T ) -external field (B) phase diagrams of Ce 1−x La x B 6 is completed by analyzing specific heat capacity, elastic neutron scattering, and the magnetocaloric effect (MCE). The x-T phase diagram reveals competitions between multipolar energy scales at low-temperatures below 4 K. Local maxima of the Sommerfeld coefficient γ 0 are developed on phase boundaries while overall magnitude of γ 0 increases with x. The analysis of the MCE rectifies the previous x-B phase diagram at T =0, and reveals a high-entropy region in lower-left corner of the modified phase diagram. The systematic change of γ 0 (x, B) on the x-B phase diagram explicitly indicates a large positive correlation between fluctuating multipoles and the effective electron mass in the heavy Fermi-liquid state of Ce 1−x La x B 6 . was determined by integrating the intensity either along a longitudinal scan through the R(1/2 1/2 1/2) point, performed for Fig 2a, or through the integration of a rocking scan through the same point in Figs. 2b and 2c.

Magnetocaloric effect
Quasi-adiabatic condition were created by supporting the sample holder with thin nylon wires. The relaxation of the sample temperature to the temperature of the mixing chamber of the dilution-fridge is well described by a heat transfer equation and the relaxation time is about 4 hours from 1 K. Therefore, a sweep-rate has been chosen so that the total period of a field-sweep is shorter than 4 hours. Under the above conditions, the sample temperature is monitored as a function of the external magnetic field.

FIG. 6:
(color) Surface plot of the specific heat capacity divided by temperature obtained from Ce0.5La0.5B6. External magnetic field is applied along [110] direction and regions of phases I, II, and IV are briefly noted.

Appendix B: Representative Surface plot of the specific heat capacity
Specific heat capacity under constant pressure and field, C p,B , of Ce 0.5 La 0.5 B 6 is measured as a function of temperature. The external magnetic fields has been varied from 0 T to 4 T with 0.25 T of interval. Then, the obtained C p,B (T ) curves compose the framework for the surface plot. Field-insensitive peaks around 1 K and below 2.5 T are clearly shown and we suspect these peaks indicate rather broadened phase transition between phase I and phase IV. As the field is further increased, phase IV is suppressed and phase II emerges. This can be recognized by the field-induced (1/2 1/2 1/2) Bragg intensity (not shown). Phase transitions are not very clear in 2D plots but these are more distinguishable in the 3D plot (Fig. 6).

Appendix C: Analysis of the quasi-adiabatic MCE
Here, we derive the differential equation for the theoretical quasi-adiabatic temperature, T 0 qad (B). First, we consider a term for the genuine (reversible) MCE, and then add a term which reflects an excessive field-induced heating. Given that the specimen with a virtual temperature T 0 qad is linked to the mixing chamber at the temperature of T mix , there should also be a term representing the natural relaxation of T 0 qad to T mix . Taking pieces together, we can write the general differential equation with respect to dT 0 qad /dB as where the double signs are in the same order and the equation with the upper (lower) signs is for the increasing (decreasing) field. dT rev is an infinitesimal temperature change following an isentropic line. dT ± irr is an infinitesimal temperature change due to a field-induced heating other than an eddy current heating. It is quantitatively discussed in the last paragraph of the current section that the eddy current heating is negligible. The thermal conduction coefficient K(T, B) between the sample and the mixing chamber is extracted by analyzing the relaxation behavior of T qad,B (t) to the T mix , and we found K(T, B) = αT δ with negligible B-dependence below 5 T. Here, α ≃ 100 nW K −1 and δ=1.2. r in the last term is the field sweep-rate.
By experience, we already know that many of analytic expressions for thermodynamic functions are written by various combinations of hyperbolic functions. Hence, T rev ∝ {tanh(B − B c ) + 1} provides a good starting point to construct a reversible ansatz for the MCE since this function well depicts stepwise temperature change upon the phase transition. Here, dT 0 qad /dB is the largest at an arbitrary critical field B c . We further elaborate this idea to imitate various shapes of isentropic lines calculated in the main text with the following expression, where C rev 1 determines the size of the stepwise MCE in a narrow range of field close to B c , and w rev 1 determines sharpness of the transition. Shapes of isentropic lines in a wider range of field are determined by adjusting C rev 2 , w rev 2 , and w rev 3 . All the blue curves for the reversible ansatz in Fig. 4 of the main text are generated by eq. (2). B i and B f are the initial and the final fields for the eq. (1), respectively.