Robust orbital diamagnetism of correlated Dirac fermions in chiral Ising universality class

We study orbital diamagnetism at zero temperature in $(2+1)$-dimensional Dirac fermions with a short-range interaction which exhibits a quantum phase transition to a charge density wave (CDW) phase. We introduce orbital magnetic fields into spinless Dirac fermions on the $\pi$-flux square lattice, and analyze them by using infinite density matrix renormalization group. It is found that the diamagnetism remains intact in the Dirac semimetal regime as a result of a non-trivial competition between the enhanced Fermi velocity and magnetic-field-induced mass gap, while it is monotonically suppressed in the CDW regime. Around the quantum critical point (QCP) of the CDW phase transition, we find a scaling behavior of the diamagnetism characteristic of the chiral Ising universality class. This defines a universal behavior of orbital diamagnetism in correlated Dirac fermions around a QCP, and therefore the robust diamagnetism in the semimetal regime is a universal property of Dirac systems whose criticality belongs to the chiral Ising universality class. The scaling behavior may also be regarded as a quantum, magnetic analogue of the critical Casimir effect which has been widely studied for classical phase transitions.


I. INTRODUCTION
Orbital diamagnetism of conduction electrons is a fundamemtal property of a material. Intuitively, it arises through the Lorentz force acting on electrons' kinetic motions and therefore it is susceptible to the band structure of the system considered. Especially in a semimetal with linear dispersions, the Landau level structure is qualitatively different from that in a conventional parabolic band system. This leads to anomalous magnetic responses in Dirac semimetals, and their orbital magnetic moment M shows extremely strong diamagnetism with non-analytic dependence on the magnetic field at zero temperature, M ∼ − √ B in two spatial dimensions. This is much stronger than that in conventional metals, M ∼ −B, for small magnetic fields. Extensive theoretical studies have been done mainly for non-interacting Dirac systems 1-15 even in a mathematically rigorous manner 16 , and various properties of diamagnetism have been theoretically discussed such as finite temperature effects 6 , roles of Berry phase 7 , lattice effects, 8 , effects of an elastic life time 9 , effects of a non-zero gap 10 , disorder effects 11 , and weak interaction effects [13][14][15] . In a realistic finite size sample with surfaces, an edge current will flow along the sample surface and generate orbital diamagnetism, where net edge currents are generally robust to surface conditions 12,[17][18][19][20] . Experimentally, strong diamagnetism has indeed been observed in several systems such as graphene and bismuth, and they are well understood based on free electron models as a direct consequence of the Dirac band structures 6,[21][22][23] . Furthermore, the origin of diamagnetim has been identified as orbital contributions in Sr 3 PbO 24 and Bi 1−x Sb x 25 .
Recently, there have emerged a variety of strongly interacting Dirac electron compounds such as molecular crystal α-(BEDT-TTF) 2 I 3 26 , magnetic layered system EuMnBi 2 27 , perovskite oxides Ca(Sr)IrO 3 28 , and twisted bilayer graphene 29 . Given these experimental developments, it is natural to ask how the orbital diamagnetism behaves in a correlated Dirac system. According to the previous theoretical study for graphene with the longrange Coulomb interaction 13 , the orbital diamagnetization M is enhanced if one takes into account the Fermi velocity (v F ) renormalization since M is proportional to v F in the Dirac semimetal phase. However, it is known that an external magnetic field induces a Dirac mass in presence of an electron interaction, which is known as magnetic catalysis [30][31][32][33][34] . The Dirac mass generally suppresses diamagnetism and it competes with an enhancement of the Fermi velocity. A perturbation study for graphene suggests that orbital magnetization is suppressed at zero temperature as a result of the non-trivial competition of these two opposite effects 15 .
Suppression of diamagnetism may occur also in other related systems where effects of mass generations due to the magnetic catalysis are stronger than those of Fermi velocity renormalizations. There are several kinds of magnetic catalysis corresponding to distinct types of field-induced orders, such as antiferromagnetism, superconductivity, and charge density wave (CDW) order. The critical behaviors around a quantum critical point (QCP) of the semimetal-insulator phase transition have been well established mainly in absence of a magnetic field [35][36][37][38][39][40][41][42][43][44][45][46][47][48] , and quantum criticality of magnetic catalysis can also be understood in a similar manner 49 . The scaling analysis shows that the Fermi velocity v F remains regular around a QCP 50 , but numerical calculations demonstrate that v F decreases to some extent in interacting Dirac fermions which exhibit antiferromagnetic quantum phase transitions [51][52][53][54] and furthermore the reduced v F was observed in the molecular compound α-(BEDT-TTF) 2 I 3 55 . Therefore, it is natural to expect suppression of diamagnetism in these systems similarly to the long-range Coulomb interacting graphene 15 . However, it is not clear whether or not diamagnetism generally gets suppresed also in other interacting Dirac systems.
In this work, we study orbital diamagnetism in a representative model of interacting Dirac fermions exhibiting a CDW order by unbiased numerical calculations with the infinite density matrix renormalization group (iDMRG) [56][57][58][59][60][61] . In this system, the Fermi velocity increases in presence of the interaction 62 , and therefore it is a promising candidate system to realize robust diamagnetism. Indeed, we demonstrate that the orbital diamagnetization remains intact for weak interactions in the Dirac semimetal regime, while it monotonically decreases as the interaction strength is increased in the insulating regime. Furthermore, the orbital magnetization M exhibits a universal scaling behavior near the QCP, and the robust orbital magnetization is characterized as a universal property of Dirac systems whose criticality belongs to the chiral Ising universality class. Besides, the scaling behavior of M is analogous to a seemingly unrelated phenomenon, the critical Casimir effect which has been extensively studied for classical phase transitions. Our study would provide a fundamental understanding of the orbital diamagnetism in correlated Dirac fermions based on the quantum critical scaling.

II. MODEL
We consider the t-V model for spinless fermions on a π-flux square lattice (also called staggered fermions) at half-filling under a uniform magnetic field, which is one of the simplest realizations of interacting Dirac fermions similarly to the honeycomb lattice model [35][36][37][38][39][40][41][42][43]49,63 . The model has two Dirac cones in the Brillouin zone corresponding to four component Dirac fermions in total. The magnetization arises only from the electron orbital motion since there is no spin degrees of freedom, which enables us to directly study the orbital magnetism. The Hamiltonian is given by where i, j is a pair of the nearest neibghbor sites. The hopping t ij = t (0) ij exp(iA ij ) contains the vector potential in the string gauge as shown in Fig. 1 corresponding to an applied uniform magnetic field 49,64 , where t (0) j+x,j = t exp(iπy j ) and t (0) j+ŷ,j = t corresponding to the π-flux lattice. In the iDMRG calculation, the system is an infinite cylinder whose size is L x × L y = ∞ × L y with the periodic boundary condition for the y-direction, and we introduce superlattice unit cells with the size L x ×L y . The magnetic field is assumed to be spatially uniform and is an integer multiple of the unit value allowed by the superlattice size, B = n × δB (n = 0, 1, 2, · · · , L x L y ) where δB = 2π/L x L y . We consider two different system sizes L y = 6, 10 to discuss finite size effects, and typically use L x = 20 for L y = 6 and L x = 10 for L y = 10. It turns out that the system can be regarded as a two dimensional system when the magnetic length l B = 1/ √ B is effectively shorter than the system size L y 49 , which enables us to study (2+1)-dimensional physics by iDMRG. We also simulate a system with L y = 14 only at zero magnetic field, which will be touched on at Sec. IV. The system size in the present study is rather limited, but we will demonstrate that our results are consistent with those obtained in the previous studies for larger system sizes at B = 0 [35][36][37][38][39][40][41][42][43] . This consistency supports our discussions on the system in presence of the magnetic field which is less understood. In this study, we use the open source code TeNPy 60,61 . In absence of a magnetic field, the Hamiltonian has sublattice Z 2 symmetry which is related to the chiral symmetry at low energy. This symmetry is spontaneously broken for large interactions above the critical strength, V > V c = 1.30t, and a charge density wave (CDW) state is realized where Dirac fermions acquire a dynamical mass [39][40][41]49 . The CDW order parameter in the ground state has been discussed previously, and it was shown that the CDW state is stabilized for any small V > 0 in presence of the magnetic field B = 0 when the system size L y is large enough 49 , which is called the magnetic catalysis [30][31][32][33][34] . Especially near the QCP, V V c , the CDW order parameter behaves as 0.54, ν 0.80 corresponding to the N = 4 chiral Ising universality class. Note that, if present, the long-range part of the Coulomb interaction would be less important around the QCP and would not affect the criticality at zero magnetic field 42,47,51-53 .
Because the system is gapped at B = 0 for any V > 0 due to the magnetic catalysis, the iDMRG numerical calculations with finite bond dimensions χ are stable and extrapolation χ → ∞ works well. In the present study, we extrapolate the calculated ground state energy density at finite χ ≤ 1600 to χ → ∞ to obtain the true ground state energy density ε for each set of V, B, and L y . De-tails of the extrapolation are discussed in Appendix A. All the numerical results in the following discussion are extrapolated ones.

III. NUMERICAL RESULTS
Firstly, we breifly exlpain qualitative behaviors of the ground state energy density ε in simple limiting cases before discussing numerical results of the iDMRG calculations. In the free Dirac fermions with a linear dispersion, single-particle energies are ∼ l −1 B with degeneracy ∼ l −2 B , which leads to the ground state energy These qualitative behaviors should hold not only deep inside each phase but also in general interaction strength in the phases. Besides, the low energy Lorentz symmetry of the Dirac semimetal phase is kept up to V = V c

35-43
B holds also at the QCP 49 . This scaling will be confimed later. Now we show the ground state energy density ε(V, B) as a function of the magnetic field B calculated by iDMRG with extrapolation χ → ∞ in Fig. 2, where ε(V, B = 0) has been shifted for the eyes. (The results ε(V = 0) have been simply obtained by direct diagonalization of the non-interacting Hamiltonian for sufficiently long cylinder geometry.) We see that the results for two different system sizes L y = 6, 10 coincide for relatively large magnetic fields B 0.02B 0 where the magnetic length l B is effectively shorter than L y , although there are some deviations for small magnetic fields B 0.02B 0 with longer l B . Therefore, finite system size effects are negligible as long as the magnetic length is effectively shorter than the system size L y as previously mentioned. This means that our system is essentially two-dimensional with the size ∼ l B × l B . In this scheme, we focus only on the magnetic length effectively shorter than L y = 6. It may seem difficult to discuss (2+1)-dimensional physics because l B < L y = 6 is too small, and in general, L y = 6, 10 would not be sufficient to discuss the true two dimensionality. However, as was demonstrated in the previous study for the magnetic catalysis 49 , it is indeed possible to investigate (2+1)-dimensional physics with this range of magnetic fields and the resulting physical quantities are consistent with those obtained by the previous studies for larger system sizes at B = 0 30-43 . This consistency supports our argument based on the small system sizes for the present model. The valid range of l B (< L y ) would be changed and correspondingly results could be improved if we include numerical data for larger system sizes, although it is computationally expensive. It has also been shown that a scaling analysis at zero magnetic field with an even shorter length scale works well for a related model 43 .
By numerically fitting the discrete data for B 0.02B 0 in our system, one can obtain continuum curves which smoothly connect them. To this end, as discussed above, in the CDW phase. Then, we introduce the following fitting functions so that their leading functional forms are consistent with these behaviors, where a j are fitting parameters. We have also included the higher order terms. We note that the zero-field energy a 0 is robust to the fitting even when we include further higher order terms in l Given the extrapolated ground state energy density ε, we can now evaluate the orbital magnetization, Although the magnetic field has to be a continuum variable in this formula, it is discrete B = n × δB in our calculations and we find that numerical differentiation δε/δB is not so reliable as will be seen in the following. Therefore, we mainly focus on the fitting function ε fit (B) and differentiate it analytically to obtain the magnetic moment M . In Fig. 3. we show the results obtained from the fitting function ε fit (B) (solid curves), and also the direct forward differentiation of the calculated discrete data with symbols for a comparison. For weak interactions, we clearly see that M (V = 0) and M (V = 0.5t) are very close each other, and think that the small difference is not so physically relevant as will be revisited later. The robustness of M (V ) for small V is understood as a result of non-trivial cancellations of two opposite effects. Firstly, at B = 0, the Fermi velocity v F (V ) is increased by V according to the recent numerical studies of the t-V model 62 and it remains regular even at the critical point according to the scaling analysis 50 , v F ∼ (V c − V ) ν(z−1) with the dynamical critical exponent z = 1. In a simple Fermi liquid picture around the non-interacting limit, the orbital magnetization of Dirac fermions is expected to be renormalized roughly as M (V )/M (0) ∝ v F (V )/v F (0). At the same time, however, the magnetic catalysis generating fermion mass ∝ B at weak interactions [30][31][32][33][34]49 will suppress the magnetization M . As a result of the non-trivial cancellation between these two effects, M can remain almost unchanged for small V . Furthermore, we aruge that the near constant M (V, B) in presence of the weak interaction V is not specific to the present model Eq. (1) and it is a general property of Dirac systems whose criticality belongs to the chiral Ising universality class, as will be discussed later based on a scaling analysis. This can be compared with the previous results for long-range Coulomb interacting Dirac electrons where the magnetic catalysis is dominant over the Fermi velocity enhancement at zero temperature and consequently the diamagnetism is suppressed 15 . In addition, suppression of diamagnetism is expected to occur in other correlated Dirac fermions where the Fermi velocities are decreased by the interactions 51-55 . In Fig. 3, as the magnetic field becomes stronger, |M (V = 0.5t)| becomes even larger than |M (V = 0)| within the present model calculation. This behavior is related to the subleading terms in ε fit (B) and it seems to be a non-universal, model-dependent property at least in Fig. 3. This point will also be revisited later.
When the interaction becomes stronger V V c , the B-dependence of the energy density ε gets weaker, which means that the orbital magnetic moment M is simply suppressed by the interaction V , as seen in Fig. 3. By increasing the interaction, the magnetization M decreases monotonically with the qaulitative change from We see that M ∼ B indeed holds in the direct numerical differentiation of the discrete data and they agree well with the fitting result. As the interaction increases further, V → ∞, the Dirac mass becomes larger and finally the orbital magnetization approaches zero, M → 0. To elucidate universal aspects of the diamagnetism in the present Dirac system whose criticality belongs to the chiral Ising universality class, we now introduce a scaling ansatz for the singular part of the ground state energy density in the thermodynamic limit L y → ∞, where g is the reduced interaction g = (V − V c )/V c 49,65-67 . The scaling dimension of g is y g = 1/ν with the correlation length exponent ξ ∼ |g| −ν , and the dimensionality is D = 2 + z = 3 with the dynamical critical exponent z = 1 for the present Lorentz symmetric criticality. This scaling ansatz can describe the critical behaviors around the QCP as a function of B, which belongs to the chiral Ising universality class with four component Dirac fermions in the present case. The proposed scaling ansatz is formally similar to the conventional finite size scaling ansatz for the isotropic system size L in absence of the magnetic field, ε sing (g, L −1 ) = b −D ε sing (b yg g, bL −1 ). These two ansatzes are related through the energy density at non-zero l −1 B , L −1 : We have the ansatz Eq.(4) for l B L → ∞, while the conventional one is obtained for l B → ∞ L. In the previous study on the same model Eq. (1) 49 , we have shown that the scaling ansatz similar to Eq. (4) indeed holds and obtained the critical exponents ν = 0.80(2), β = 0.54 (3), and the critical interaction strength V c = 1.30(2)t, where β is the CDW order parameter exponent M CDW ∼ g β for g ≥ 0. These values are consistent with those obtained in other studies at zero magnetic field [36][37][38][39][40][41][42][43] . In the present study, we simply use these previous results and examine quantum criticality of the orbital magnetization.
From the scaling ansatz Eq.(4) with b = l B , the total ground state energy density is regarded as a function of the single variable gl We have included a correction to scaling to improve the scaling description, and similar corrections with respect to the system size L have been often used in numerical calculations 38 , although physical origin of the introduced corrections may not be so clear in general. In this study, we regard the correction to scaling as a working ansatz to evaluate large l B behaviors in a systematic way. The scaling function Φ(x) is universal in the sense that it is determined only by the universality class, and it should be independent of boundary conditions of the system since Eq. (5) is the energy density in the thermodynamic limit. This is in sharp contrast to the conventional finite size corrections which depends on boundary conditions. Note that the universal function Φ(x) behaves as Φ(x −1) ∼ const. corresponding to ε − ε 0 ∼ l −3 B in the Dirac semimetal phase for g < 0 , while Φ(x 1) ∼ x −ν corresponding to ε − ε 0 ∼ l −4 B in the CDW phase for g > 0. Around the QCP, Φ(x) should be analytic in x since there would be no phase transitions for any nonzero l −1 B , and Φ(0) at g = 0 may contain some useful information about the criticality as will be discussed later.
To show a scaling plot of ε(g, l −1 B ), we use ε 0 (g) = a 0 (g) from Eq. (2) which are robust to details of the fitting. Then, the calculated ε collapse onto a single curve as shown in Fig. 4 with the critical exponent ν = 0.80 and critical interaction V c = 1.30t 49 . Here, the interaction range is relatively wide, V = 0.50t ∼ 2.0t, and the magnetic length is l B 3.2l B0 ∼ 7.1l B0 measured in unit of l B0 = 1/ √ 2π. The overall behavior of Φ data (gl is consistent with the above mentioned general expectation. Based on this observation, we introduce the following fitting function as a working ansatz, where α 0 = −α 1 and α j are parameters to be determined from numerical fitting with the calculated data (see Appendix B for details). The solid curve in Fig. 4 is the fitting function Φ fit (x) and it well agrees with the data (variance of residuals χ 2 = O(10 −3 )). Once the scaling function has been obtained, we can find the universal scaling of the orbital magnetization M = −∂ε/∂B = (1/2l 3 B )∂(Φl −3 B )/∂l B with suppressing the non-universal correction term near the QCP for sufficiently large l B , This equation clearly means that the orbital magnetization in the form M(gl B) is a universal function of gl 1/ν B characteristic of the associated quantum criticality, namely, the N = 4 chiral Ising universality class in D = (2 + 1)-dimensions. We show M obtained from Φ fit in Fig. 5. For a comparison, we also show the results calculated with forward differentiation of the numerical data. Although the numerical differentiation of our data is less accurate due to its discreteness, overall behaviors are in agreement with the one obtained from the analytic differentiation of Φ fit (x). As explained above, the scaling function behaves as Φ(x −1) ∼ const and therefore we have , which means M ∼ −B in the CDW phase as expected. At the QCP, the magne- where the amplitude is expected to be universal as will be discussed in the next section.   Fig. 6 (a), the diamagnetism is highly robust to the interaction. For example at a small magnetic field, B/B 0 = 0.01, the orbital magnetization is robust up to the interaction V /t 1.0 and then sharply drops in the quantum critical regime around the QCP, V c = 1.30t, as seen in Fig. 6 (b). Such a behavior is commonly seen for other small values of B and M (V ) as a function of V gets smeared for larger values of B. In the CDW regime, M is strongly suppressed by the interaction. The crossover lines separat-ing the different regimes are roughly given by |gl In addition, by looking at the magnetization M closely, one can see that M is slightly enhanced by the interaction V at small magnetic fields, and it is free from non-universal (model dependent) finite size corrections in contrast to Fig. 3. We stress that the present argument is based on the scaling ansatz, and therefore the resulting robust and slightly enhanced orbital magnetization in presence of the interaction V should be a common property in general Dirac systems whose criticality belongs to the chiral Ising universal class as long as the ansatz holds. As already mentioned, the robust M in presence of the interaction V is a consequnce of the non-trivial competition between the renomalized Fermi velocity and magnetic catalysis. The present results imply that the former effect is dominant over the latter when V is weak in the chiral Ising universality class. This is contrasted with the diamagnetism in the Hubbard models and also in systems with the long-range Coulomb interaction, where the Fermi velocity renormalization cannot cancel out effects of the magnetic catalysis 15,51-54 .

IV. DISCUSSION AND SUMMARY
As menioned before, the scaling behavior Eq. (5) is seemingly similar to the conventional finite system size scaling at zero magnetic field, ε(g, L −1 ) = ε 0 (g) + Φ(gL 1/ν )/L D + · · · . The leading finite size correctioñ Φ(0)/L D is called the Casimir energy density in field theories and contains universal information of the criticality. In D = 1+1 dimensions, the Casimir amplitude is written asΦ(0) ∼ cv with a boundary condition dependent coefficient, where v is the speed of light (velocity of excitations) and c is the central charge of the underlying conformal field theory [68][69][70] . Generalizations to higher dimensional systems have been first discussed for a cylinderical spacetime geometry 68 and also recently examined in torus and infinite systems 62,71,72 . It was proposed that the Casimir amplitude in a (2 + 1)-dimensional torus system is de-composed asΦ torus (0) = C torus v, where C torus contains some universal information of the underlying field theory. For a comparison, we also calculate the Casimir energy density in our model as shown in Fig. 7, where the ground state energy density is assumed to behave as ε(g = 0, l −1 B = 0, L −1 y ) = ε 0 (0) +Φ iDMRG (0)/L D y + · · · , as in other Lorentz symmetric critical systems 62,71,72 . The amplitude is found to be negativeΦ iDMRG (0) < 0 in contrast to Φ(0) > 0 in Eq. (5). Within infinite projected entangled pair states (iPEPS) calculations, the correlation length ξ D due to a finite bond dimension can be a new length cut-off scale in a thermodynamically large system and will play a similar role to that of the system size L, leading to ε = ε 0 + C iPEPS v/ξ 3 D + · · · at g = 0 72 . We have demonstrated in this study that, in presence of a magnetic field, the leading term Φ(0)/l 3 B in a thermodynamically large system can be regarded as a magnetic analogue of the conventional Casimir energy in a finite size system, and may be called magnetic Casimir energy. Note that similarity between conventional Casimir energy and magnetic Casimir energy is already implied in single-particle spectra of the non-interacting Dirac electrons; (L) ∝ v F /L in a finite size system without a magnetic field, and (l B ) ∝ v F /l B in an infinite system with a magnetic field. Structures of the single-particle spectrum are governed by the Lorentz symmetry and the characteristic length scale is either L or l B . These properties are also common to a correlated Dirac system around a QCP, leading to similar functional forms of the critical Casimir energies. However, there is an essential difference between them; the conventional Casimir energy is geometry (boundary condition) dependent, while the magnetic Casimir energy is independent of boundary conditions since it is the energy in the thermodynamic limit. Besides, the magnetic Casimir energy can be controlled by an external magnetic field, which is a difference from the Casimir energy in iPEPS for an infinite system that is a purely theoretical quantity.
It is known that the Casimir energy L 2Φ (gL 1/ν )/L D leads to the critical Casimir force f L = −∂E/∂L, E = L 2 ε, and related physics has been extensively studied in various systems which exhibit finite temperature classical phase transitions [73][74][75][76][77][78][79][80][81][82][83] . The universal nature of the classical critical Casimir force has been experimentally observed for example in a binary liquid mixture in thin film geometry [74][75][76] . The quantum critical orbital diamagnetism discussed in this study can be regarded as a magnetic, quantum analouge of the classical critical Casimir force. Indeed, the magnetization is rewritten as and the "effective force" f B is repulsive for diamagnetism. This is contrasting to the attractive real space Casimir force in our model (Fig. 7), corresponding to the sign difference between Φ(0) andΦ iDMRG (0). The repulsiveness of the effective force f B could be compared with general Casimir forces, where they are usually attractive and repulsive forces are rarely realized 84,85 .
Interestingly, the critical orbital magnetization M or equivalently the effective force f B could be measured in experiments by carefully controlling experimental parameters, which may be advantageous over the formidable challenge for a direct observation of the real space Casimir force in a solid crystal. In a bulk magnetization measurement, observing the magnetization M is just equivalent to measuring f B , and one can understand the above analogy in a visible manner. For example, the effective force f B could be measured by the conventional magnetization measurement with the Faraday balance, where a mechanical force f z (z) under a macroscopically non-uniform magnetic field B(z) is observed in the typical setup shown in Fig. 8 (a). The effective force is related with the mechanical force simply as f z = −∂E/∂z = f B · ∂l B /∂z. Besides, note that the direction of f B is determined by the macroscopic configuration of the magnetic field, since the energy density ε(l B ) of a diamagnetic Dirac system favors a larger l B (a smaller B). If the z-axis magnetic field varies in the xy-plane in a macroscopic length scale as shown in Fig. 8  (b), the effective force f B will be parallel to the plane. Such an in-plane force would be more analogous to the critical Casimir force f L which also acts within the plane, but the system may exhibit very complex behaviors in an inhomogeneous setup.
In summary, we have discussed orbital diamagnetism in correlated Dirac electrons with use of iDMRG for the t-V lattice model which exhibits the CDW quantum phase transition. The orbital diamagnetism is robust to the short-range interaction V in the Dirac semimetal regime, while it is suppressed for a strong V in the CDW regime. The robustness of the diamagnetism to the interaction is understood as a consequence of non-trivial competition between the enhanced Fermi velocity and mass generation by the magnetic catalysis. Furthermore, it is concluded that the robust orbital diamagnetism is a universal property of Dirac systems in the chiral Ising universality class based on the scaling analysis in terms of mag- If B(z) is decreasing (lB(z) is increasing) with z, the effective force is repulsive fB > 0 (or equivalently, the real force is fz = −∂E/∂z > 0). Similary, for B(y) decreasing in the y-direction, the effective force will be fB +ŷ (or fy = −∂E/∂y > 0). netic length. The analogy between the quantum critical diamagnetism and critical Casimir effect was discussed. To our best knowledge, this is a first unbiased numerical calculation of the orbital diamagnetism in strongly interacting electrons other than one-dimensional systems [86][87][88][89][90] , and it could provide a basis for further theoretical developments in this field. consider only the extrapolated energy density ε(χ → ∞) in our discussion. FIG. 9. Extrapolation of the ground state energy density ε for the χ → ∞ limit at V = 0.5t. The blue (red) symbols are for Ly = 6(Ly = 10), and the curves are polynomial fittings. The corresponding magnetic fields for each curve are B/δB = 0, 1, 2, 3, 4, 5, 6 from the bottom to the top, where δB = 2π/(6 × 20) for Ly = 6, L x = 20 and δB = 2π/(10 × 10) for Ly = 10, L x = 10. As marked by a broken circle, at B = 2π/20 = 2π × 6/(6 × 20) = 2π × 5/(10 × 10), the extrapolated ε(χ → ∞) for two system sizes coincide within the error (which is smaller than the symbols). On the other hand, ε(χ → ∞) for Ly = 6, 10 differ at B = 0 because the corresponding magnetic length is lB = ∞ Ly.
Appendix B: Details of fitting with Eq. (6) The numerically obtained scaling function Φ data (·) is consistent with the limiting behaviors of the ground state energy density at x → ±∞. For example, Φ data (x 1) ∼ x −ν can indeed be confirmed in Fig. 10. We also find that Φ data (x −1) is roughly Φ data (x) − Φ data (−∞) ∼ |x| −ν (not shown), which leads to a natural behavior, ε = ε 0 + const × l −3 B + const × l −4 B + · · · in the Dirac semimetal phase. These observations enable us to evaluate the universal scaling function Φ(·) in Eq. (5) by using simple known functions and to introduce the fitting function Φ fit (·) in Eq. (6). We note that it is not trivial to have a successful scaling plot over a wide range of x = gl 1/ν B where Φ(x) does not have a simple Taylor expansion with a small order in x. However, such a non-trivial scaling has been often examined in classical statistical models for the Casimir effect [73][74][75][76][77][78][79][80][81][82][83] . The successful description of scaling behaviors of the ground state energy density consequently suggests that the scaling ansatz Eq. (5) is indeed correct, which is a priori nontrivial. Together with the previous scaling argument for the CDW order parameter 49 , the present study supports the magnetic length scaling ansatz for which there are only few studies 65-67 . for Ly = 6 (squares) and Ly = 10 (circles). The solid line represents the qualitative behavior ∼ x −ν with ν = 0.80.