Robust orbital diamagnetism in correlated Dirac fermions

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 π-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, while it is monotonically suppressed in the CDW regime. Around the quantum critical point of the CDW phase transition, we find a scaling behavior of the diamagnetism characteristic of the chiral Ising universality class. Besides, the scaling analysis implies that the robust orbital diamagnetism at weak magnetic fields in a Dirac semimetal regime would hold not only in our model but also in other interacting Dirac fermion systems as long as scaling regions are wide enough. 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.


Introduction
Orbital diamagnetism of conduction electrons is a fundamental 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][2][3][4][5][6][7][8][9][10][11][12][13][14][15] even in a mathematically rigorous manner [16]. 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 diamagnetism 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 long-range 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 suppressed also in other interacting Dirac systems. The orbital diamagnetism in interacting systems is theoretically less explored compared with the spin magnetism partly because of its technical difficulty, and further studies are necessary to develop a deeper understanding of diamagnetism.
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. The robust orbital magnetization at weak magnetic fields is a direct consequence of the scaling behavior and thus similar robustness would be seen in other interacting Dirac fermion models as long as scaling regions are wide enough. 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.

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-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 neighbor sites. The hopping t ij = t (0) ij exp(iA ij ) contains the vector potential in the string gauge as shown in figure 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 section 4. 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 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 M CDW ∼ l −β/ν B ∼ B β/2ν with β 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][52][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 . Details of the extrapolation are discussed in appendix A. All the numerical results in the following discussion are extrapolated ones.

Numerical results
In this section, we discuss the numerical results obtained by iDMRG calculations. We compute the orbital magnetization directly from the ground state energy density in section 3 1, and then discuss its scaling behaviors in section 3 2.

Direct evaluation of diamagnetism
Firstly, we briefly explain 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 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][36][37][38][39][40][41][42][43] and B holds also at the QCP [49]. This scaling will be confirmed 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 figure 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][31][32][33][34][35][36][37][38][39][40][41][42][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, we first observe 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 −1 B . 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 figure 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.50t) 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 can be 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], In a simple Fermi liquid picture around the non-interacting limit, the orbital magnetization of Dirac fermions is expected to be renormalized roughly as . 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 argue that the near constant M(V, B) at small magnetic fields is not specific to the present model equation (1) and similar behaviors would hold in other interacting Dirac fermion models, as will be discussed later based on a scaling analysis. This can be compared with the previous results for long-range Coulomb interacting Dirac fermions where the magnetic catalysis is dominant over the Fermi velocity enhancement at zero temperature and consequently the diamagnetism is suppressed [15]. Suppression of diamagnetism is expected to occur also in on-site interacting Hubbard models where the Fermi velocities are decreased by the interactions [51][52][53][54][55]. In figure 3, as the magnetic field becomes stronger, |M(V = 0.50t)| 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 figure 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 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.

Scaling analysis
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][66][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 y g g, bL −1 ). These two ansatzes are related through the energy density at non-zero l −1 B , L −1 : we have the ansatz equation (4) for l B L → ∞, while the conventional one is obtained for l B → ∞ L. In the previous study on the same model equation (1) [49], we have shown that the scaling ansatz similar to equation (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 equation (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 equation (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 equation (2) which are robust to details of the fitting. Then, the calculated ε collapse onto a single curve as shown in figure 4 with the critical exponent ν = 0.80 and critical interaction V c = 1.30t [49]. Here, the interaction range is relatively wide, 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 figure 4 is the fitting function Φ fit (x) and it well agrees with the data [variance of residuals 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 figure 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  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 figure 3. On the other hand, in the CDW regime, M is strongly suppressed by the interaction. The crossover lines separating the different regimes are roughly given by |gl We point out 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 would not be limited to the simplest t-V model on the π-flux square lattice. Indeed, the mechanism for the robust diamagnetism in the present t-V model is quite simple and it would be basically applicable to other Dirac systems as well; if the scaling function M(x) is valid for all −∞ < x = gl Therefore, there are plateaus in the normalized magnetization at small magnetic fields if the scaling region is wide enough, which may be somewhat similar to the Binder ratio in a spin model. This gives an alternative understanding of the robust diamagnetism which was qualitatively attributed to the non-trivial competition between the Fermi velocity renormalization and magnetic catalysis. Although the scaling region is model dependent, we can expect that the diamagnetism is robust to the dominant interaction even if other small interactions are added such as next nearest neighbor interactions because the scaling region would change smoothly with the perturbations. Besides, our argument based on the scaling analysis would be applicable also to other lattices for Dirac fermions such as the honeycomb lattice. It is an interesting future problem to study effects of interactions other than the present V-term and Dirac fermions on the honeycomb lattice. Extensions of the present study to other models such as spinful Hubbard models where the Fermi velocity decreases with the on-site interaction are also a future issue [51][52][53][54].

Discussion and summary
We have examined the orbital diamagnetism as a finite size effect of the magnetic length. Based on this observation, we discuss an analogy between the diamagnetism and a seemingly unrelated phenomena, the Casimir effect in section 4 1. Finally, a summary is given in section 4 2.

Analogy to casimir effect
As mentioned before, the scaling behavior equation (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 correctionΦ(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 Y Tada Figure 7. The ground state energy density at V = V c without a magnetic field, ε(g = 0, l −1 B = 0, L −1 y ), for system sizes L y = 6, 10, 14 with the periodic boundary condition for y-direction. The qualitative behavior is assumed to be ε − ε 0 ∼ 1/L 3 y , as indicated by the solid line.
have been first discussed for a cylinderical space-time 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 decomposed 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 figure 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 equation (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 fermions; (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 analogue 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 (figure 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 figure 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 figure 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.

Summary
In summary, we have discussed orbital diamagnetism in correlated Dirac fermions 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 can be understood as a consequence of non-trivial competition between the enhanced Fermi velocity and mass generation by the magnetic catalysis. Furthermore, the scaling analysis implies that the robust orbital diamagnetism would hold not only in our model but also in other interacting Dirac systems as long as the corresponding scaling regions are wide enough. 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 fermions other than one-dimensional systems [86][87][88][89][90], and it could provide a basis for further theoretical developments in this field. Figure 9. Extrapolation of the ground state energy density ε for the χ → ∞ limit at V = 0.5t. The blue (red) symbols are for L y = 6(L y = 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 L y = 6, L x = 20 and δB = 2π/(10 × 10) for L y = 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 L y = 6, 10 differ at B = 0 because the corresponding magnetic length is l B = ∞ L y . Figure 10. The scaling plot of the ground state energy density ε(g, l −1 B ) at large gl 1/ν B for L y = 6 (squares) and L y = 10 (circles). The solid line represents the qualitative behavior ∼x −ν with ν = 0.80.

Appendix B. Details of fitting with equation (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 figure 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 equation (5) by using simple known functions and to introduce the fitting function Φ fit (·) in equation (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 equation (5) is indeed correct, which is a priori non-trivial. 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][66][67].