Basic formulation and first-principles implementation of nonlinear magneto-optical effects

The conventional formulae for second-order photo-current and electric polarization only consider non-degenerate bands, and therefore are only U(1)-gauge invariant [25, 13]. However, in magnetic materials with $\mathcal{PT}$ symmetry, Kramers-like degeneracy appears at every $k$ point and the choice of eigenstates in the doubly degenerate subspace has a gauge freedom determined by a $2\times 2$ unitary matrix. As physical results are invariant in spite of an arbitrary unitary rotation in the degenerate subspace, the formulae for physical observables under $\mathcal{PT}$ symmetry should be U(2)-gauge invariant. Clearly, certain terms are missing in the conventional formula in the presence of $\mathcal{PT}$ symmetry and in what follows, we will demonstrate that the U(2)-invariant formulae can be retrieved with special treatment at degeneracy [25].

In the length gauge and under long-wavelength approximation, the electromagnetic wave is introduced through an $\mathbf{E}(t)\cdot\mathbf{\hat{x}}$ term to the unperturbed Hamiltonian $H_{0}$ where a monochromatic light at frequency $\omega$ is written as $\mathbf{E}(t)=\mathbf{E}(\omega)e^{-i\omega t}+$ c.c. The Hamiltonian in second quantization form can be found in Appendix A.1. The evaluation of the position operator $\mathbf{\hat{x}}$ is the key to the calculation of polarization and current responses, and the matrix elements of position operator can be written as [25]

$|n_{\mu}\mathbf{k}\rangle=e^{i\mathbf{k}\cdot\hat{\mathbf{x}}}|u_{n_{\mu}}(\mathbf{k})\rangle$ is the eigenstate of $H_{\rm 0}$, and $|u_{n_{\mu}}(\mathbf{k})\rangle$ is the periodic part of the Bloch function. We introduced two subscripts $n$ and $\mu$ for bands, where $n$ denotes bands with different energy, and $\mu$ is used in the subspace of degenerate bands. $\mu$ runs from 1 to 2 under $\mathcal{PT}$ symmetry, while the following derivations hold for an arbitrary $\mu$ value. The diagonal and off-diagonal parts of $\mathbf{\xi}_{n_{\mu}m_{\nu}}(\mathbf{k})=i\langle u_{n_{\mu}}(\mathbf{k})|\nabla_{\mathbf{k}}|u_{m_{\nu}}(\mathbf{k})\rangle$ are the single band and two-band Berry connections $\mathbf{A}_{n_{\mu}}(\mathbf{k})$ and $\mathbf{r}_{n_{\mu}m_{\nu}}(\mathbf{k})$ defined as

In the following formulae, we omit the $\mathbf{k}$ dependence for simplicity. The first (second) term of Eq.(1) is the intraband (interband) position matrix element in an infinite crystal and the position operator matrix element between degenerate bands ($n=m,\mu\neq\nu$) are considered in the interband term. As a result, the polarization can be divided into $\mathbf{{P}}=e\mathbf{{x}}=\mathbf{{P}}_{\rm inter}+\mathbf{{P}}_{\rm intra}$ in which $e=-\absolutevalue{e}$ is the charge of an electron. Similarly, the current density along $a$ direction which is defined as the time derivative of electric polarization is divided into two parts as

where $\hat{H}(t)$ is the total Hamiltonian including the oscillating electric field.

With the density matrix derived with perturbation expansion in Appendices A.1 and A.2, the zeroth and first order responses including the finite frequency anomalous conductivity are derived in Appendix A.3. The second-order electric current and polarization with incoming electric field of frequency $\omega_{\beta}$ and $\omega_{\gamma}$ and the response frequency $\omega_{\Sigma}=\omega_{\beta}+\omega_{\gamma}$ are

and

where the full expression of $\langle\hat{P}_{\rm inter}^{a}(t)\rangle^{(2)}$ can be found in Appendix A.3 and the Einstein summation convention is adopted in all formulae for repeated indices. $\hbar$ is the reduced Planck’s constant and $[d\mathbf{k}]=d\mathbf{k}/(2\pi)^{d}$ is the unit volume in reciprocal space of $d$ dimension. An infinitesimal $1/\tau$ is introduced to avoid divergence at resonance, which can also be regarded as a relaxation rate. $a,b,c$ represent Cartesian directions, and $E_{\beta}$ is the shorthand notation for $E(\omega_{\beta})$ which represents the amplitude of an electric field with frequency $\omega_{\beta}$. $(b\beta\leftrightarrow c\gamma)$ is the interchange of $b\beta$ and $c\gamma$ indices which reflects the fact that physical responses are irrelevant to the sequence of electric fields $E^{b}_{\beta}$ and $E^{c}_{\gamma}$ in the expression. $f_{nm}$, $\omega_{nm}$ and $\Delta_{mn}$ are the occupation, energy and group velocity differences between the $n$-th and the $m$-th bands. In gapped semiconductor at zero temperature, $f=1$ for valence band and $f=0$ for conduction band. $r_{n_{\mu}m_{\nu};a}^{b}=\partial_{k_{a}}r_{n_{\mu}m_{\nu}}^{b}-i(A^{a}_{n_{\mu}}r_{n_{\mu}m_{\nu}}^{b}-r^{b}_{n_{\mu}m_{\nu}}A^{a}_{m_{\nu}})$ is the U(1)-covariant derivative of $r_{n_{\mu}m_{\nu}}^{b}$, where the U(1)-covariant derivative of a non-degenerate eigenstate $|u_{n}(\mathbf{k})\rangle$ is $|u_{n}(\mathbf{k})\rangle_{;a}=\partial_{k_{a}}|u_{n}(\mathbf{k})\rangle+iA^{a}_{n}|u_{n}(\mathbf{k})\rangle=(1-|u_{n}(\mathbf{k})\rangle\langle u_{n}(\mathbf{k})|)\partial_{k_{a}}|u_{n}(\mathbf{k})\rangle\,$, which represents the projection of $\partial_{k_{a}}|u_{n\mathbf{k}}\rangle$ onto all states except itself.

It is convenient to use Fourier transformation to convert results into the frequency domain and define the conductivity and susceptibility in the frequency domain as

In the following, we will write $j_{\rm intra}^{a}(\omega)^{(2)}$ instead of $\langle\hat{j}^{a}_{\rm intra}(\omega)\rangle^{(2)}$ as abbreviation to represent the expectation value.

A changing electric polarization gives rise to a current, and the rectification current $(\omega_{\beta}=-\omega_{\gamma}=\pm\omega,\omega_{\Sigma}=0)$ is of particular interest. If bands are non-degenerate which is equivalent to dropping the ${\mu},{\nu}$ indices and $\omega_{ln}\neq 0$ in Eq. (5), $\mathbf{P}_{\rm inter}(\omega_{\Sigma}=0)^{(2)}$ is finite, with $\omega_{\Sigma}\mathbf{P}_{\rm inter}(\omega_{\Sigma}=0)^{(2)}=0$, and $\mathbf{j}(\omega_{\Sigma}=0)^{(2)}=\mathbf{j}_{\rm intra}(\omega_{\Sigma}=0)^{(2)}$ according to Eq. (3), which can reduce to the well-known U(1)-invariant formulae of rectification current densities [25, 13]. However, if band degeneracy occurs, e.g. $\omega_{ln}=0$, degenerate bands contribute a polarization that is diverging as $\frac{1}{\omega_{\Sigma}}$ when $\omega_{\Sigma}$ approaches the static limit according to Eq. (5).

The divergent polarization at $\omega_{\Sigma}\to 0$ corresponds to a DC current as

which contains the Berry connection between degenerate bands, e.g. $r_{n_{\mu}n_{\lambda}}$. The total rectification current is

where ${D^{a}}[r_{n_{\mu}m_{\nu}}^{b}]=\partial_{k_{a}}r_{n_{\mu}m_{\nu}}^{b}-i\sum_{\lambda}(r^{a}_{n_{\mu}n_{\lambda}}r_{n_{\lambda}m_{\nu}}^{b}-r^{b}_{n_{\mu}m_{\lambda}}r^{a}_{m_{\lambda}m_{\nu}})$ is the U($N$)-covariant derivative of $r_{n_{\mu}m_{\nu}}^{b}$ and $N$ is the number of degeneracy determined by the index $\mu$. Similarly, the U($N$)-covariant derivative of degenerate eigenstates just projects $\partial_{k_{a}}|u_{n_{\mu}}(\mathbf{k})\rangle$ out of the $N$ dimensional degenerate subspace such that $D^{a}|u_{n_{\mu}}(\mathbf{k})\rangle=\partial_{k_{a}}|u_{n_{\mu}}(\mathbf{k})\rangle+iA^{a}_{n_{\mu}}|u_{n_{\mu}}(\mathbf{k})\rangle+i\sum_{\nu\neq\mu}r^{a}_{n_{\nu}n_{\mu}}|u_{n_{\nu}}(\mathbf{k})\rangle=(1-\sum_{\nu}|u_{n_{\nu}}(\mathbf{k})\rangle\langle u_{n_{\nu}}(\mathbf{k})|)\partial_{k_{a}}|u_{n_{\mu}}(\mathbf{k})\rangle\,$. Up to now, we have presented a thorough derivation of the second-order electric current and have revealed the microscopic origin of the U($N$)-covariant derivative. In Appendix B, we demonstrated an alternative derivation of the above formulae by including the position operator matrix between degenerate bands into the intraband term in Eq. (1).

In the intermediate region where the double degeneracy protected by $\mathcal{PT}$ symmetry is weakly broken with the energy splitting $E_{\rm deg}$, the fundamental formalisms in Eq. (3 - 5) should be used to calculate the rectification current. As $E_{\rm deg}$ increases from 0 to finite, the rectification current turn to zero which is consistent with the U(1)-invariant result and an AC current peaked at $E_{\rm deg}$ appears.

Using the Sokhotski–Plemelj theorem, the denominator can be simplified as $(\omega_{mn}-\omega-i/\tau)^{-1}={\rm P}(\omega_{mn}-\omega)^{-1}+i\pi\delta(\omega_{mn}-\omega)$ where P indicates principal part. Therefore, the total rectification current can be divided into four parts as normal injection current (NIC), magnetic injection current (MIC), normal shift current (NSC), and magnetic shift current (MSC)

with expressions

The integrand $I_{mn}^{abc}=\sum_{{\mu}{\nu}}r^{b}_{m_{\nu}n_{\mu}}D^{a}[r^{c}_{n_{\mu}m_{\nu}}]$, the two-band Berry curvature $\Omega^{bc}_{mn}=-2\sum_{{\mu}{\nu}}\Im{r_{m_{\nu}n_{\mu}}^{b}r_{n_{\mu}m_{\nu}}^{c}}$, the two-band quantum metric $g^{bc}_{mn}=\sum_{\mu\nu}\Re{r_{m_{\nu}n_{\mu}}^{b}r_{n_{\mu}m_{\nu}}^{c}}$ [7][8][39], and the group velocity difference $\Delta_{mn}$ are all gauge invariant and their symmetry properties are listed in Table 1. The above results are consistent with literature [13, 25], and additionally we generalized the formulae of NSC and MSC to a U($N$)-invariant form in the presence of $N$-fold degeneracy.

The general formula for second-order polarization can be derived through $\mathbf{P}(\omega)^{(2)}=\mathbf{j}_{\rm intra}(\omega)^{(2)}/({-i\omega})+\mathbf{P}_{\rm inter}(\omega)^{(2)}$ in which Berry connection terms containing degenerate bands are already included in $\mathbf{P}_{\rm inter}(\omega)^{(2)}$ in the conventional expression. Therefore, the expressions of polarization and susceptibility are compatible with degenerate bands and no additional modifications are required to fulfill the gauge invariant requirement.

However, two terms in the susceptibility with leading factors $\omega_{\Sigma}^{-1}$ and $\omega_{\Sigma}^{-2}$ originated from $\mathbf{j}_{\rm intra}^{(2)}$ give rise to a divergent problem in the $\omega_{\Sigma}\to 0$ limit with the divergent terms

where the imaginary part $i/\tau$ in the denominator is not written out explicitly. The divergence is unphysical in an insulator with a finite energy gap. Even if $\chi_{\rm dvg}$ take a $0/0$ form and the divergence can be removed eventually, the expression in Eq. (14) is unfriendly to numerical calculations at low frequency. This problem has been noticed in SHG susceptibility ($\omega_{\beta}=\omega_{\gamma}=\omega,\omega_{\Sigma}=2\omega$) by Sipe and Ghahramani [22] and they discovered that the divergence in SHG can be removed by applying $\mathcal{T}$ symmetry. However, the simplification does not work in magnetic materials and as a results, the current SHG expressions might be problematic in studying magnetic materials. Furthermore, applying $\mathcal{T}$ symmetry fails to remove the divergence in the general susceptibility in Eq. (14). Therefore, the second-order susceptibility for arbitrary frequencies still suffers from an unphysical divergent problem even for non-magnetic materials.

In what follows, we present an alternate expression for the divergent term which is given by

where $\rho_{\beta}=\omega_{\beta}/\omega_{\Sigma}$ and $\rho_{\gamma}=\omega_{\gamma}/\omega_{\Sigma}$. As $\omega_{\Sigma}$ is completely removed from the denominator, the above expression does not have divergence and numerical instability in the zero-frequency limit and can be applied to both non-magnetic and magnetic materials. Eq. (15) is derived from Eq. (14) by the following techniques without assuming $\mathcal{T}$ symmetry: exchange of dumpy band indices $n$ and $m$, integration by parts, and the equality $-\omega^{-1}\left[(\omega_{mn}-\omega)^{-1}+(\omega_{nm}-\omega)^{-1}\right]=\omega_{nm}^{-1}\left[(\omega_{mn}-\omega)^{-1}-(\omega_{nm}-\omega)^{-1}\right]$. Although the U(1)-derivative appears in the first term of Eq. (15), the full susceptibility can be rearranged in a way that terms with Berry connection between degenerate bands are grouped with terms with the U(1)-covariant derivative to form the U($N$)-covariant derivative.

As SHG is one of the most commonly used tool for symmetry detection and frequency conversion, we take SHG as an example to discuss properties of the ‘divergent’ terms. The full expression of SHG is given in the Appendix C and the ‘divergent’ terms are

The first term in Eq. (16) survives under $\mathcal{T}$ symmetry and is included in the conventional SHG expression [22, 24, 28, 29, 30, 31, 32, 33, 34, 35, 36]. In addition, the first term is purely real in the static limit. However, the second term in Eq.(16) is present only when $\mathcal{T}$ symmetry is broken and therefore it is a unique term for magnetic materials. We denote the second term by $\chi_{\rm new}(-2\omega;\omega,\omega)$ as it has not been calculated previously. Unexpectedly, we discovered that although $\chi_{\rm new}$ has a low-frequency ‘divergent’ form in Eq. (14), $\chi_{\rm new}(\omega=0)$ does not contribute to static SHG at all, which can be concluded by permutation of the dummy band indices in the summation in Eq. (16). Therefore, the low-frequency ‘divergent’ term $\chi_{\rm new}$ can be observed away from low-frequency regime but negligible in low-frequency regime.

The coefficients of second-order electro-optical responses are third-rank tensors, and their nonzero and independent elements are determined by their magnetic point group [40]. As the electric field, polarization and current are all odd under $\mathcal{P}$, the second-order electro-optical responses are only present in $\mathcal{P}$-breaking materials, according to Eq. (6). Additionaly, special attentions are paid to $\mathcal{T}$ and $\mathcal{PT}$ symmetry and the symmetry properties of quantum metric, two-band Berry curvature, and etc. are summarized in Table 1. As a result, both NSC and NIC are even under $\mathcal{T}$ while MSC and MIC are even under $\mathcal{PT}$ symmetry according to Eqs. (10 - 13). Furthermore, the current coupled to $\Re{E^{b}(\omega)E^{c}(-\omega)}$ can be induced both by a linearly-polarized light (LPL) and a circularly-polarized light (CPL), while the current coupled to $\Im{E^{b}(\omega)E^{c}(-\omega)}$ can only be induced by CPL with the direction of current determined by the helicity of light [13][23].

Up to now, we have focused on the charge photocurrent which is the collective motion of both spin-up and spin-down electrons. As electrons also carry the spin degree of freedom, the spin current can also exists. Although it is still challenging to properly define the spin current in materials with spin-orbit coupling (SOC) as spin is not a good quantum number in this circumstance [41], the conventional definition $j^{ab}=1/2(v^{a}s^{b}+s^{b}v^{a})$ with $v$ and $s$ the velocity and spin operators respectively is still appropriate without SOC [42][43].

We adopted the conventional definition to analyse the symmetry of spin current. As spin is $\mathcal{T}$ odd and $\mathcal{P}$ even, the symmetry of spin current and charge current are opposite under $\mathcal{T}$ and $\mathcal{PT}$ symmetry. For instance, while the charge current from the NSC and NIC mechanisms are $\mathcal{T}$ even, the spin current from those mechanisms are $\mathcal{T}$ odd. As a result, spin current from NSC and NIC mechanisms only exists in magnetic materials. The symmetry analysis of both charge and spin currents are summarized in Table 2. Additionally, in whatever situation, the spin current and charge current are both present. For instance, the charge NSC and spin MIC are present simultaneously in $\mathcal{T}$-symmetric materials under LPL.

Although the four mechanisms listed in Table 2 can generate both charge and spin currents, the microscopic origins are different. $\mathcal{P}$ breaking is the precondition of second-order charge and spin currents and it usually happens in the crystallographic structure. However, if it is the spin order instead of the geometric structure that breaks $\mathcal{P}$, the presence of SOC is imperative for a non-zero charge current but not necessary for generating spin current.

Despite the intricate expression of second-order susceptibility, certain symmetry conditions are fulfilled. Firstly, second-order susceptibility has the intrinsic permutation symmetry which states that $\chi_{2}^{abc}(-\omega_{\Sigma};\omega_{\beta},\omega_{\gamma})$ is unchanged by the simultaneous interchange of its last two frequency arguments and its last two Cartesian as $\chi_{2}^{abc}(-\omega_{\Sigma};\omega_{\beta},\omega_{\gamma})=\chi_{2}^{acb}(-\omega_{\Sigma};\omega_{\gamma},\omega_{\beta})$. The intrinsic permutation symmetry is introduced in Eqs. (4) and (5) as a convention. Secondly, when all frequencies are detuned from resonance and the small imaginary part in the denominator can be ignored, the nonlinear susceptibility possesses full permutation symmetry which states that all of the frequency arguments of the nonlinear susceptibility can be freely interchanged, as long as the corresponding Cartesian indices are interchanged simultaneously: $\chi_{2}^{abc}(-\omega_{\Sigma};\omega_{\beta},\omega_{\gamma})=\chi_{2}^{bca}(\omega_{\beta};\omega_{\gamma},-\omega_{\Sigma})$. The full permutation symmetry can be deduced from a consideration of the form of the electromagnetic field energy within a lossless medium [44]. Thirdly, in the static limit, the non-zero components of the susceptibility tensor are all equal $\chi_{2}^{abc}(0;0,0)=\chi_{2}^{bca}(0;0,0)=\chi_{2}^{cab}(0;0,0)$, which is an extension of the full permutation symmetry in the static limit, also called the Kleinman’s symmetry. Moreover, the static susceptibility is a purely real quantity with the expression given in Appendix C.

Furthermore, the second-order susceptibility which is odd in $\mathcal{P}$ can be decoupled into two contributions, a non-magnetic response $\chi_{\rm N}$ which is even in $\mathcal{T}$ and a magnetic response $\chi_{\rm M}$ which is odd in $\mathcal{T}$ with their expressions given in Appendix C. $\chi_{\rm N}$ reflects the inversion symmetry breaking from the lattice and charge density while $\chi_{\rm M}$ is sensitive to the inversion symmetry breaking by the magnetic ordering. The principal part and $\delta$-function part of the non-magnetic $\chi_{\rm N}$ are purely real and imaginary respectively, while the real and imaginary parts are opposite in the magnetic $\chi_{\rm M}$ [45, 46]. In other words, $\chi_{\rm N}$ is purely real and $\chi_{\rm M}$ is purely imaginary away from resonance. As the susceptibility in the static limit is purely real, $\chi_{\rm M}$ does not contribute to the static SHG.

Generally, the imaginary part of the susceptibility is responsible for energy dissipation of the electromagnetic field. However, in magnetic materials, $\chi_{\rm M}$ is imaginary even away from resonance, therefore, the meaning of dissipation needs to be scrutinized. The derivation of energy dissipation from the first-order process is in Appendix D. The dissipation rate of the electromagnetic energy due to the second-order polarization $\mathbf{P}^{(2)}$ is given by

For frequencies detuned from resonance, the full-permutation symmetry of $\chi_{2}$ guarantees that the term in the bracket is zero. Therefore, only the $\delta$-function part of susceptibility participates in the dissipation no matter whether it is real or imaginary.

Moreover, $\chi_{\rm N}$ and $\chi_{\rm M}$ also corresponds to different measurable quantities. The sum-frequency generation (SFG) intensity ($I$) measured in experiments can be decoupled into a non-magnetic term, a magnetic term and interference terms as

For materials with $\mathcal{T}$ or $\mathcal{PT}$ symmetry, the SFG intensity is proportional to either the $\chi_{\rm N}$ term or the $\chi_{\rm M}$ term. When both $\mathcal{T}$ and $\mathcal{PT}$ symmetries are broken, e.g. in a polar ferromagnetic, the appearance of the third and the fourth terms reflects the interference between the magnetic and non-magnetic signals and the sign of the interference is determined by the the direction of the magnetic and the polar orders. As a result, domains with opposite magnetic/polar orders could show different SFG intensities which makes SFG an important tool for multiferroics detection.

In the following, we present the implementation of the U(2)-invaraint formulae using both orthogonal and non-orthogonal basis. The general derivative is usually evaluated through a sum-over state method [25, 47] or the Wannier interpolation method [36, 48]. We adopted the Wannier interpolation method proposed in Ref. [36] for orthogonal basis and the scheme proposed in Ref. [49] for non-orthogonal basis, while the Berry connection in both methods is evaluated based on non-degenerate perturbation theory. Although methods with orthogonal and non-orthogonal basis sets have different gauge choices, as the optical responses are gauge invariant, the two methods are equivalent in principle. While the Wannier interpolation method is more computationally efficient, the method with non-orthogonal basis is more suitable for high-throughput calculations [49]. For SHG calculations, the three-band summation in Eq. (36) requires a large number of bands to achieve convergence, therefore, the method with non-orthogonal basis is a more suitable choice to evaluate SHG responses.

We generalized the above mentioned methods by including degenerate perturbation to evaluate the Berry connection and the general derivative, which are unavoidable in the presence of $\mathcal{PT}$ symmetry. For the Wannier interpolation method, we chose a specific gauge which satisfies that $(U^{\dagger}(\mathbf{k}_{0})U(\mathbf{k}_{0}+\delta\mathbf{k}))_{n_{\mu}n_{\nu}}$ is purely real for degenerate bands $n_{\mu}$ and $n_{\nu}$. Therefore, the non-Abelian connection $U^{\dagger}\partial_{k_{a}}U$ between two degenerate states vanishes, while the non-Abelian connections between non-degenerate states can be calculated as usual [36]

where $H^{(W)}$ is the Hamiltonian matrix in the representation of local Wannier functions and $U$ is a matrix of which each column is an eigenvector of $H^{(W)}$. For the method of non-orthogonal basis, we again chose a specific gauge so that the matrix elements between degenerate bands in Ref. [49] are modified as

Here $S$ is the overlap matrix of the non-orthogonal basis, $v$ is a matrix with each column an generalized eigenvector of $H^{(W)}$ which is in the representation of a set of complete but not orthogonal basis. More details can be found in Ref. [49].

First-principles calculations were performed both by Vienna $Ab\ initio$ Simulation Package (VASP) [50] which is a plane-wave basis package and by OpenMX which is a pseudo-atomic basis package [51, 52]. The exchange-correlation functional was parameterized in the PBE form [53]. PAW pseudopotentials [54] and norm-conserving pseudopotentials [55] were used in VASP and OpenMX, respectively. We have included the effect of SOC. For the 3$d$ orbitals in magnetic ions Mn and Cr, the Hubbard $U$ of 4 eV and 3 eV were used, respectively [56]. For layered materials, we used DFT-D3 form van der Waals correction without damping [57]. In VASP, the cut-off energy of plane waves was set to 350 eV and 450 eV for MnBi₂Te₄ and CrI₃, respectively. In OpenMX, Cr6.0-s3p3d2 and I9.0-s3p3d2f1 were chosen as our basis for Cr and I. The convergence criterion for force and electronic calculations were 10 meV/Å and $10^{-6}$ eV. $k$-point samplings of $15\times 15\times 1$ and $13\times 13\times 1$ were used for MnBi₂Te₄ and CrI₃ in VASP, and a $7\times 7\times 1$ $k$-mesh was adopted for CrI₃ in OpenMX.

After getting the converged electronic structures, we either generated the maximally localized Wannier functions using Wannier90 [58] or used the non-orthogonal atomic basis from OpenMX to build the tight-binding (TB) Hamiltonian and calculate the optical responses [36, 49]. We obtained 100 maximally localized orbitals for MnBi₂Te₄ and 112 maximally localized orbitals for CrI₃. In the calculation of MIC and MSC, we adopted a $400\times 400\times 1$ kmesh. For SHG susceptibility, the relaxation rate $\hbar/\tau$ was set to be 0.05 eV, and the $k$-mesh was set to be $50\times 50\times 1$ which leads to identical results as $100\times 100\times 1$ $k$-point sampling. The degenerate perturbation was applied when the energy difference between two bands are smaller than $E_{\rm deg}$ which was set to be 0.5 meV in our calculations.

The flowchart of our computations including the electronic structure and optical response calculations is illustrated in Fig. 1. Three different routes were adopted to calculate optical responses. In route I, the electronic structure is calculated by VASP and then maximally localized Wannier functions are constructed to obtain the TB Hamiltonian in an orthogonal basis set. The difference between route II and I is that the electronic structure is calculated by OpenMX rather than VASP. In route III, after the self-consistent calculation by OpenMX, we directly get the TB Hamiltonian in non-orthogonal basis. The purpose of using multiple packages (VASP and OpenMX) and different basis sets (orthogonal and non-orthogonal) for TB Hamiltonian is to examine the influence of different electronic structures and gauge choices for optical responses separately.

We considered two prototypical materials: bilayer MnBi₂Te₄ [59] and bilayer CrI₃ [60], both of which have interlayer A-type AFM magnetic orders. Both materials exhibit $\mathcal{PT}$ symmetry, which implies that the spatial inversion symmetry is broken only through the AFM spin order between two layers while the geometric structure is centrosymmetric as illustrated in Fig. 2. In this case, only magnetic responses including MIC, MSC and magnetization-sensitive SHG are present while normal responses NIC, NSC and non-magnetic SHG are forbidden, which makes them a simple platform to investigate the nonlinear magneto-optical responses. Additionally, both materials are feasible in experiments [21, 60] and their MIC has been studied theoretically in literature[12, 13, 12, 61], which provides valuable results for comparison. For bilayer MnBi₂Te₄, we considered AB-type stacking (same stacking pattern as in bulk) in which each layer laterally shifted by [1/3, 1/3, 0]. For bilayer CrI₃, multiple stacking orders have been investigated previously [60], and we focused on the one with AB^′ stacking in which each layer laterally shifted by [1/3, 0, 0] as both MIC and MSC are allowed. The stacking order, magnetic order, magnetic symmetry and independent tensor components are summarized in Table 3.

The MIC of bilayer MnBi₂Te₄ and bilayer CrI₃ have been calculated recently [13, 61]. Although previous works haven’t considered the U(2)-gauge problem under $\mathcal{PT}$ symmetry, fortunately, MIC is not subjected to modifications in $\mathcal{PT}$ symmetric case as demonstrated in Eq. (13). In the following, we performed calculations and compared with previous works to validate our method and results.

In bilayer MnBi₂Te₄, its interlayer coupling is AFM and the intralayer coupling is FM. The energy band is highly dispersive near $\Gamma$ point with a small gap of 0.076 eV at $\Gamma$ point as shown in Fig. 3(a) . In addition, the comparison between our band structure (red solid line) and the one from Ref. [13] (blue dotted line) along high symmetry line showed excellent agreement near band edge. Fig. 3(b) shows the MIC conductivity of bilayer MnBi₂Te₄ along $xxx$ direction. Our result again shows an excellent agreement with the Ref. [13].

For bilayer CrI₃, in order to compare the band structure and MIC, we adopted the same Hubbard $U=1$ eV and the same Cartesian coordinates shown by dashed line in Fig. 2(d) as Ref. [61]. Fig. 3(c, d) shows the calculated band structure and MIC conductivity. The band gap of bilayer CrI₃ is 0.89 eV at $\Gamma$ point in our calculation and 0.78 eV in the reference, therefore we applied a 0.11 eV scissor operation to results in the reference to align the band gap in Fig. 3(c) to 0.89 eV. While the shape of bands is similar near band edge, noticeable differences exist in many bands. In addition, the magnetic moment of Cr atom along $z$ axis is 3.12 $\mu_{B}$ in our calculation, while 3.21 $\mu_{B}$ in the reference. As we adopted the same parameters (Hubbard $U$, exchange-correlation functional, pseudopotentials, Van der Waals correction) in the electronic structure calculation, the above discrepancies might result from small differences in atomic structures. For the MIC shown in Fig. 3(d), the main characters of our results (red solid line) are similar to the reference (blue dashed line) including the location and height of the first peak. However, due to discrepancies in band structures away from the band edge, the location and height of other peaks are shifted and modified. Therefore, combining the results of bilayer MnBi₂Te₄ and CrI₃, we can conclude that we are able to reproduce MIC results in literature while the detailed features of MIC are sensitive to geometric and electronic structures.

The correction to ensure U(2)-gauge invariant is essential in the presence of $\mathcal{PT}$ symmetry which hasn’t been considered in the previous computational works. Here we calculated the MSC conductivity of AB^′ stacking bilayer CrI₃ through three different routes shown in Fig. 1 and alternating between U(1)-covariant derivative and U(2)-covariant derivative to demonstrate the influence of the gauge choice. As the results are very sensitive to details in electronic structure (see Fig. 6 in the Appendix), we stay with the electronic structure obtained from OpenMX for comparison.

Figure 4(a) shows the results from orthogonal basis method (blue dashed line) versus non-orthogonal basis method (red solid line) (route II vs route III in Fig. 1) with the same U(2)-gauge-invariant formula. Since the two methods have different gauge choices, the excellently agreed results in Fig. 4(a) confirm that our formula of MSC is truly gauge independent.

In contrast, Fig. 4(b) shows the results of orthogonal basis method (blue dashed line) versus non-orthogonal basis method (red solid line) (route II versus route III) both using U(1)-form formula adopted in previous works. Although the same electronic structure guarantees that the peak positions resulted from the delta function in Eq. (11) are exactly the same, differences in peak height are prominent, especially at low-frequency region. For instance, they differ by more than three times at 1.2 eV and there is a shoulder in the red solid line at 0.8 eV while absent in the blue dashed line. More importantly, the differences demonstrate that the previous formula with U(1)-covariant derivative is gauge dependent and the results are not reproducible due to different gauge choices. Therefore, the previous MSC formula with U(1)-covariant derivative cannot be used to describe physical observables in $\mathcal{PT}$-symmetric materials.

Furthermore, Fig. 4(c) shows results from U(1)-invariant formula (blue dashed line) versus U(2)-invariant formula (red solid line) using non-orthogonal basis (both are in route III). Noticeable differences are observed at peak values due to the improper implementation of the covariant derivative in $\mathcal{PT}$-symmetric materials, implying the necessity of using U(2)-invariant formula.

In terms of the SHG response in magnetic materials, there are three sets of formulas that are available to use: the ‘DVG’ expression which is the original expression in Ref. [22] with low-frequency divergent problem shown in Eq. (14), the ‘$\mathcal{T}$-SMP’ expression in which the low-frequency divergent terms are simplified by $\mathcal{T}$ symmetry and have been implemented in many calculation packages [30, 31, 32, 33, 34, 35, 36], and the full-frequency convergent expression ‘CVG’ proposed in this paper as demonstrated in Eq. (15). As the gigantic SHG susceptibility in AB^′ CrI₃ has been measured [21] and calculated using the $\mathcal{T}$-SMP expression recently [38], we compare the SHG response of bilayer CrI₃ along $xxy$ direction calculated using the DVG, $\mathcal{T}$-SMP and CVG expressions.

Figure 5(a) shows the real and imaginary parts of SHG calculated using the DVG and the CVG expressions. As the CVG expression is derived from the DVG expression, results from both expressions are exactly the same away from the low-frequency region ($>$ 0.1 eV). However, in the low-frequency region ($0-0.1$ eV), both the real and imaginary parts calculated from the DVG expression exhibit divergent problem which is absent in the CVG expression. Additionally, in the static limit, the CVG expression converges to exactly zero which is consistent with our symmetry analysis in Sec. III.

Figure 5(b) shows the real and imaginary parts of SHG calculated using the $\mathcal{T}$-SMP and the CVG expressions. As the SHG response in bilayer CrI₃ is gigantic, the difference between the two expressions is hardly observed in Fig. 5(b). Detailed analyses in Fig. 5(c) demonstrates that the differences between the two expressions, which is $\chi_{\rm new}$, are roughly 10 pm/V. The materials and frequency range in which $\chi_{\rm new}$ is considerable is subject to further investigation with more magnetic materials.