Orbital magneto-optical response of periodic insulators from first principles

Abstract

Magneto-optical response, i.e. optical response in the presence of a magnetic field, is commonly used for characterization of materials and in optical communications. However, quantum mechanical description of electric and magnetic fields in crystals is not straightforward as the position operator is ill defined. We present a reformulation of the density matrix perturbation theory for time-dependent electromagnetic fields under periodic boundary conditions, which allows us to treat the orbital magneto-optical response of solids at the ab initio level. The efficiency of the computational scheme proposed is comparable to standard linear-response calculations of absorption spectra and the results of tests for molecules and solids agree with the available experimental data. A clear signature of the valley Zeeman effect is revealed in the continuum magneto-optical spectrum of a single layer of hexagonal boron nitride. The present formalism opens the path towards the study of magneto-optical effects in strongly driven low-dimensional systems.

Introduction

Magneto-optical phenomena originating from the loss of symmetry between left and right circularly polarized light in the presence of a magnetic field are widely used for characterization of different kinds of matter.^1,2 Magnetic circular dichroism (MCD) spectra help to assign overlapping bands and give insight into magnetic properties of the ground and excited states. Faraday rotation of the plane of polarization of linearly polarized light serves as a basic operational principle for functional magneto-optical disks and optical isolators.³ Optical excitations in the presence of a magnetic field allow manipulation of valley pseudospin degrees of freedom in two-dimensional monolayers.^4,5,6,7,8,9 Giant Faraday rotation has been revealed in graphene¹⁰ and metal oxide nanosheets.¹¹ These advances cultivate the growing interest to development of a gauge-invariant and computationally efficient ab initio theory of magneto-optical response.

Although ab initio calculations of MCD spectra in molecules can be performed nowadays in a nearly routine fashion^{12,13,14,15,16} (as implemented in quantum chemistry codes^17,18), the complete response theory for extended systems is still under development. The reason is that external electromagnetic fields break the translational symmetry of such systems, which in the formal way is expressed through unboundness of the position operator. Although according to the modern theory of polarization,^19,20,21 the position operator can be replaced by a derivative with respect to the wave vector in responses to electric fields, the description of magnetic fields is more complicated as it introduces vector coupling to electron dynamics and leads to non-perturbative changes in wavefunctions. Three approaches have been considered in literature to deal with these difficulties: (1) taking a long-wavelength limit of an oscillating perturbation,^22,23 (2) using the Wannier function formalism^24,25,26,27, or (3) treating perturbations of the one-particle Green function or one-particle density matrix,^27,28,29 which are two-point quantities summed up over all occupied bands and having periodic and gauge-invariant counterparts. Although wave functions in the presence of even a very small magnetic field differ drastically from those in the absence of the magnetic field (a plane wave for a free electron and a localized Landau level state for an electron in the magnetic field can be considered as an example), the gauge-invariant counterpart of the density matrix changes perturbatively.^27,28,29 In approach (1), proper sum rules^30,31 should be taken into account to control numerical errors arising upon summing up non-gauge-invariant paramagnetic and diamagnetic terms. In approach (3), such a numerical noise is supressed automatically. Approach (3) also allows us to work under purely periodic boundary conditions as opposed to approach (2), where contributions of open boundaries should be treated carefully.^24,25,26

So far the magnetic field has been considered in the context of static responses.^{22,23,24,25,26,27,28,29} In the present paper we demonstrate that density matrix perturbation theory^27,29,32 can be extended to the case of dynamic non-linear phenomena. We focus on second-order magneto-optical effects, i.e., the change of the optical response in the presence of a magnetic field. Although the approach developed here is general and can be adapted to any first-principles framework, we decide to illustrate it using time-dependent density functional theory (TDDFT).^33,34 This method provides a satisfactory level of accuracy at a moderate computational cost and has been widely employed in literature for magneto-optical response of molecules.^{12,13,14,15,16} The account of excitonic effects in the transverse optical response of solids, however, is not straightforward within TDDFT and is performed here using the approach derived in ref. ³⁵ from time-dependent current density functional theory (TDCDFT).

The procedures for solids implemented for the present paper form a part of the open-source code Octopus.^36,37,38 For the sake of simplicity, we limit our consideration to orbital magneto-optical effects for insulators. Although the spin contribution is trivial, the account of the Fermi surface contribution can be done for metals by analogy with ref. ²³.

In the following, we derive the equations implemented, describe the computational scheme, give the expressions for magneto-optial properties measured experimentally, and finally discuss the results of calculations for molecules and solids.

Results and discussion

One-particle density matrix in electromagnetic fields

Let us consider the response to uniform magnetic and electric fields. We use the temporal gauge, in which both of these fields are described by the vector potential A and are given by B = ∇ × A and E = −c⁻¹∂_tA, respectively, where c is the speed of light (atomic units are used throughout the paper). Though the fields are uniform, the vector potential A entering in the Hamiltonian H is non-periodic. This gives rise to ill-defined expectation values of quantum mechanical operators describing physical properties of the system in the periodic basis. However, it turns out that for any operator ${\cal O} = {\cal O}_{{\bf{r}}_1{\bf{r}}_2}$ defined for two points r₁ and r₂ in real space it is possible to distinguish the periodic and gauge-invariant counterpart $\tilde {\cal O} = \tilde {\cal O}_{{\bf{r}}_{\mathrm{1}}{\bf{r}}_2}$ by factoring out the Aharonov–Bohm-type phase^27,28,29

$$\varphi _{12} = - c^{ - 1}{\int}_{{\bf{r}}_2}^{{\bf{r}}_1} {\kern 1pt} {\bf{A}}({\bf{r}}){\mathrm{d}}{\bf{r}}$$

(1)

so that

$${\cal O}_{{\bf{r}}_1{\bf{r}}_2} = {\tilde{\cal O}}_{{\bf{r}}_1{\bf{r}}_2}{\mathrm{exp}}\left( {i\varphi _{12}} \right).$$

(2)

Here we take ħ = e = 1 and the integral is taken along the straight line between points r₂ and r₁ so that r = r₂ + (r₁ − r₂)ξ, 0 ≤ ξ ≤ 1.

This approach was previously used to derive corrections to the gauge-invariant counterpart $\tilde \rho$ of the one-particle density matrix ρ in the static magnetic field.^27,29 In the present paper, we generalize these derivations to the case of time-dependent electromagnetic fields by rewriting the time-dependent Liouville equation

$$\begin{array}{*{20}{l}} { - i\partial _t\rho + [H,\rho ] = 0} \hfill \end{array}$$

(3)

in terms of $\tilde \rho$. Here and below the commutator of two operators ${\cal O}^{(1)}$ and ${\cal O}^{(2)}$ is introduced as

$$[ {{\cal O}^{(1)},{\cal O}^{(2)}} ]_{{\mathbf{r}}_1{\mathbf{r}}_3} = {\int} {\mathrm{d}}{\mathbf{r}}_2\left( {{\cal O}_{{\mathbf{r}}_1{\mathbf{r}}_2}^{(1)}{\cal O}_{{\mathbf{r}}_2{\mathbf{r}}_3}^{(2)} - {\cal O}_{{\mathbf{r}}_1{\mathbf{r}}_2}^{(2)}{\cal O}_{{\mathbf{r}}_2{\mathbf{r}}_3}^{(1)}} \right).$$

(4)

Using Eq. (2) for the relation between $\tilde \rho$ and ρ in real space, the time-dependent Liouville Eq. (3) gives

$$- ie^{i\varphi _{13}}\left( {\partial _t + i\partial _t\varphi _{13}} \right)\tilde \rho _{{\bf{r}}_1{\bf{r}}_3} = {\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}_2e^{i(\varphi _{12} + \varphi _{23})}\left( {\tilde \rho _{{\bf{r}}_1{\bf{r}}_2}\tilde H_{{\bf{r}}_2{\bf{r}}_3} - \tilde H_{{\bf{r}}_1{\bf{r}}_2}\tilde \rho _{{\bf{r}}_2{\bf{r}}_3}} \right).$$

(5)

It should be noted that $\tilde H = H_0 + \delta \tilde H$, where the difference $\delta \tilde H$ between the gauge-invariant counterpart $\tilde H$ of the Hamiltonian and unperturbed Hamiltonian H₀ is related to the local-field effects coming from changes in the electron density induced by the external fields and corresponds to the variation of Hartree and exchange-correlation potentials in TDDFT (see page 1 of Supplementary information).

Eq. (5) is equivalent to

$$- i\left( {\partial _t + i\partial _t\varphi _{13}} \right)\tilde \rho _{{\bf{r}}_1{\bf{r}}_3} = {\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}_2e^{i\varphi _{123}}\left( {\tilde \rho _{{\bf{r}}_1{\bf{r}}_2}\tilde H_{{\bf{r}}_2{\bf{r}}_3} - \tilde H_{{\bf{r}}_1{\bf{r}}_2}\tilde \rho _{{\bf{r}}_2{\bf{r}}_3}} \right),$$

(6)

where φ₁₂₃ = φ₁₂ + φ₂₃ + φ₃₁.

This phase corresponds to the flux of the magnetic field through the triangle formed by points r₁, r₂ and r₃:

$$\begin{array}{*{20}{l}} {\varphi _{123} = \frac{1}{{2c}}{\bf{B}} \cdot \left( {{\bf{r}}_1 - {\bf{r}}_2} \right) \times \left( {{\bf{r}}_2 - {\bf{r}}_3} \right).} \hfill \end{array}$$

(7)

The time derivative of the phase φ₁₃ on the left-hand side of Eq. (6) introduces the electric field

$$\begin{array}{*{20}{l}} {\partial _t\varphi _{13} = {\bf{E}} \cdot \left( {{\bf{r}}_1 - {\bf{r}}_3} \right).} \hfill \end{array}$$

(8)

Combining Eqs. (6)–(8), we arrive at

$$\begin{array}{*{20}{l}} { - i\left( {\partial _t + i{\bf{E}} \cdot \left( {{\bf{r}}_1 - {\bf{r}}_3} \right)} \right)\tilde \rho _{{\bf{r}}_1{\bf{r}}_3}} \hfill & = \hfill & {{\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}_2e^{i{\bf{B}} \cdot \left( {{\bf{r}}_1 - {\bf{r}}_2} \right) \times \left( {{\bf{r}}_2 - {\bf{r}}_3} \right)/2c}} \hfill \\ {} \hfill & {} \hfill & { \cdot \left( {\tilde \rho _{{\bf{r}}_1{\bf{r}}_2}\tilde H_{{\bf{r}}_2{\bf{r}}_3} - \tilde H_{{\bf{r}}_1{\bf{r}}_2}\tilde \rho _{{\bf{r}}_2{\bf{r}}_3}} \right).} \hfill \end{array}$$

(9)

This expression is gauge-invariant and includes all corrections to time-dependent electric and magnetic fields. Therefore, it can be used to derive expressions for responses of any order to electromagnetic fields.

To describe magneto-optical effects on the basis of Eq. (9), we assume that E corresponds to the oscillating electric field of the electromagnetic wave and B to the static magnetic field applied. The magnetic field of the electromagnetic wave is neglected. We, therefore, consider only the first-order corrections in E, B, and E×B. Keeping only the terms to the first order in the magnetic field is reasonable even for strong magnetic fields $B \ll c{\mathrm{/}}a^2\sim 10^5\,{\mathrm{T}}$, where a = 1 Å is taken as a typical interatomic distance.

Eq. (9) for the density matrix then takes the form

$$\begin{array}{l} - i\partial _t\tilde \rho _{{\bf{r}}_1{\bf{r}}_3} - {\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}_2 \cdot \left( {\tilde \rho _{{\bf{r}}_1{\bf{r}}_2}\tilde H_{{\bf{r}}_2{\bf{r}}_3} - \tilde H_{{\bf{r}}_1{\bf{r}}_2}\tilde \rho _{{\bf{r}}_2{\bf{r}}_3}} \right)\\ = - {\bf{E}} \cdot \left( {{\bf{r}}_1 - {\bf{r}}_3} \right)\tilde \rho _{{\bf{r}}_1{\bf{r}}_3}\\ + \frac{i}{{2c}}{\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}_2{\kern 1pt} {\bf{B}} \cdot \left( {{\bf{r}}_1 - {\bf{r}}_2} \right) \times \left( {{\bf{r}}_2 - {\bf{r}}_3} \right)\\ \cdot \left( {\tilde \rho _{{\bf{r}}_1{\bf{r}}_2}\tilde H_{{\bf{r}}_2{\bf{r}}_3} - \tilde H_{{\bf{r}}_1{\bf{r}}_2}\tilde \rho _{{\bf{r}}_2{\bf{r}}_3}} \right).\end{array}$$

(10)

Using that ${\cal O}_{{\bf{r}}_1{\bf{r}}_2}\left( {{\bf{r}}_1 - {\bf{r}}_2} \right) = [{\bf{r}},{\cal O}]_{{\bf{r}}_1{\bf{r}}_2}$ and introducing notations for the anticommutator of two operators ${\cal O}^{(1)}$ and ${\cal O}^{(2)}$

$$\left\{ {{\cal O}^{(1)},{\cal O}^{(2)}} \right\}_{{\bf{r}}_1{\bf{r}}_3} = {\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}_2\left( {{\cal O}_{{\bf{r}}_1{\bf{r}}_2}^{(1)}{\cal O}_{{\bf{r}}_2{\bf{r}}_3}^{(2)} + {\cal O}_{{\bf{r}}_1{\bf{r}}_2}^{(2)}{\cal O}_{{\bf{r}}_2{\bf{r}}_3}^{(1)}} \right)$$

(11)

and velocity operator V = $- i[{\bf{r}},\tilde H]$ computed with account of all non-local contributions to the Hamiltonian, such as from non-local pseudopotentials, Eq. (10) can be finally rewritten as

$$- i{\partial _t}\tilde \rho + \left[ {\tilde H,\tilde \rho } \right] = - \frac{1}{2}\left\{ {{\bf{E}} + \frac{1}{c}{\bf{V}} \times {\bf{B}},\left[ {{\bf{r}},\tilde \rho } \right]} \right\}{.}$$

(12)

This is simply the quantum Bolzmann equation with the Lorentz driving force on the right-hand side. Unlike the singular position operator r, the commutator $[{\bf{r}},\tilde \rho ]$ of the position operator with the periodic function $\tilde \rho$ is well defined here and can be substituted by the derivative with respect to the wave vector, i∂_kρ_k, in reciprocal space.^27,28,29

Moving the term coming from the local-field effects to the right-hand side,

$$- i\partial _t\tilde \rho + \left[ {H_0,\tilde \rho } \right] = - \frac{1}{2}\left\{ {{\bf{E}} + \frac{1}{c}{\mathrm{\bf{V}}} \times {\bf{B}},[{\bf{r}},\tilde \rho ]} \right\} - [\delta \tilde H,\tilde \rho ],$$

(13)

we get all terms dependent on the external fields on the right-hand side of the equation. Differentiating the Liouville Eq. (13), one can evaluate the derivatives of the density matrix $\tilde \rho ^{(P)} = \partial \tilde \rho {\mathrm{/}}\partial P$ with respect to perturbations P of parameters of the Hamiltonian, such as the electric field E or magnetic field B.

Numerical solution of Liouville equation

In the following, we consider solution of the Liouville Eq. (13) within TDDFT, i.e. assuming that ρ is the Kohn–Sham density matrix and H is the Kohn–Sham Hamiltonian. The same Liouville equation, however, describes magneto-optical effects in any other first-principles framework and a similar computational scheme can be used.

From the computational point of view, it is convenient to divide the nth order derivative $\tilde \rho ^{(P)}$ of the density matrix describing the joint response to the perturbations P = P₁P₂...P_n into four blocks within and between the occupied (V) and unoccupied subspaces (C):

$$\begin{array}{*{20}{l}} {\tilde \rho ^{(P)} = \tilde \rho _{{\mathrm{VV}}}^{(P)} + \tilde \rho _{{\mathrm{CC}}}^{(P)} + \tilde \rho _{{\mathrm{VC}}}^{(P)} + \tilde \rho _{{\mathrm{CV}}}^{(P)}.} \hfill \end{array}$$

(14)

These blocks correspond to $\tilde \rho _{{\mathrm{VV}}}^{(P)} = P_v\tilde \rho ^{(P)}P_v$, $\tilde \rho _{{\mathrm{CC}}}^{(P)} = P_c\tilde \rho ^{(P)}P_c$, $\tilde \rho _{{\mathrm{VC}}}^{(P)} = P_v\tilde \rho ^{(P)}P_c$, and $\tilde \rho _{{\mathrm{CV}}}^{(P)} = P_c\tilde \rho ^{(P)}P_v$, where P_v = ρ⁽⁰⁾ and P_c = 1 − P_v are the projectors onto the occupied and unoccupied bands.

Following the density matrix perturbation theory,³² to get the elements of the derivative of the density matrix $\tilde \rho _{{\mathrm{CV}}}^{(P)}$ between the unoccupied and occupied subspaces, we project the Liouville Eq. (13) onto unperturbed Kohn–Sham wavefunctions $\left| {\psi _{v{\bf{k}}}^{(0)}} \right\rangle$ of occupied bands v:

$$L_{v{\bf{k}}}({\mathrm{\Omega }})\left| {\eta _{v{\bf{k}}}^{(P)}} \right\rangle = P_cR^{(P)}\left[ {\tilde \rho ^{(n - 1)}, \ldots ,\rho ^{(0)},n^{(P)}} \right]\left| {\psi _{v{\bf{k}}}^{(0)}} \right\rangle .$$

(15)

Here the operator on the left-hand side is given by $L_{v{\bf{k}}}({\mathrm{\Omega }}) = {\mathrm{\Omega }} + H_0 - \epsilon _{v{\bf{k}}}$, where Ω is frequency considered and $\epsilon _{v{\bf{k}}}$ is the energy of the unperturbed state $\left| {\psi _{v{\bf{k}}}^{(0)}} \right\rangle$. The operator R on the right-hand side includes all terms dependent on the perturbation P coming from the right-hand side of Eq. (13) and is determined by the derivatives of the density matrix of the previous orders (see equations for each type of perturbation on pages 1–3 of Supplementary information). If the local-field effects are taken into account, it also depends on the derivative of the electron density n^(P) to the perturbation P, n^(P)(r₁) = ρ ^(P)(r₁, r₂)δ(r₁ − r₂) (see page 1 of Supplementary information).

The solution of Eq. (15) corresponds to

$$\begin{array}{*{20}{l}} {\left| {\eta _{v{\bf{k}}}^{(P)}({\mathrm{\Omega }})} \right\rangle = P_c\tilde \rho ^{(P)}({\mathrm{\Omega }})\left| {\psi _{v{\bf{k}}}^{(0)}} \right\rangle = \tilde \rho _{{\mathrm{CV}}}^{(P)}({\mathrm{\Omega }})\left| {\psi _{v{\bf{k}}}^{(0)}} \right\rangle } \hfill \end{array}$$

(16)

and once it is known, the elements $\tilde \rho _{{\mathrm{CV}}}^{(P)}$ of the derivative of the density matrix between unoccupied and occupied subspaces can be computed as

$$\begin{array}{*{20}{l}} {\tilde \rho _{{\mathrm{CV}}}^{(P)}({\mathrm{\Omega }}) = {\int}_{{\mathrm{BZ}}} \frac{{{\mathrm{d}}{\bf{k}}}}{{(2\pi )^3}}\mathop {\sum}\limits_v \left| {\eta _{v{\bf{k}}}^{(P)}({\mathrm{\Omega }})} \right\rangle \left\langle {\psi _{v{\bf{k}}}^{(0)}} \right|.} \hfill \end{array}$$

(17)

The elements between the occupied and unoccupied subspaces can be found as $\tilde \rho _{{\mathrm{VC}}}^{(P)}({\mathrm{\Omega }}) = \left( {\tilde \rho _{{\mathrm{CV}}}^{(P)}\left( { - {\mathrm{\Omega }}^ \ast } \right)} \right)^ \ast$ and to obtain them, Eq. (15) should be also solved for the frequency ${ - {\mathrm{\Omega }}^ \ast }$. If the local-field effects are taken into account, Eq. (15) has to be solved self-consistently as the derivative $\tilde \rho ^{(P)}$ of the density matrix determines the derivative of the electron density n^(P), which enters on the right-hand side of Eq. (15).

Solution of Eq. (15) is performed in the present paper using the efficient Sternheimer approach,^38,39,40,41 which corresponds to the iterative search of the function $| {\eta _{v{\bf{k}}}^{(P)}({\mathrm{\Omega }})} \rangle$ that fits into this equation at each frequency Ω. Other approaches, such as sum over states,¹⁵ methods based on Casida’s equation,^15,16 complex polarization propagator^12,13, and real-time propagation¹⁴ have been used to compute absorption and magneto-optical spectra of molecules. The sum over states, Casida’s equation⁴² and complex polarization propagator,^43,44 however, require inclusion of many well converged unoccupied states. Such calculations are not feasible for large systems, where too many KS states should be computed. They also fail to describe properly high-energy excitations due to poor convergence of the corresponding KS states. Casida’s equation⁴² furthermore relies on the use of real wavefunctions and cannot be straightforwardly extended to solids, where KS states are complex.

Neither Sternheimer approach,^38,39,40,41 nor real-time propagation¹⁴ need calculation of unoccupied states. They also have a favorable scaling of O(N²) with the system size N as compared, for example, to O(N³) for the sum over states (refs ^14,39,40). The advantage of the real-time propagation is that it makes possible calculation of responses for all frequencies at once. However, long propagation times are required to achieve a good resolution. The Sternheimer approach is more appropriate for computing the spectra in a narrow frequency region with a high resolution. The calculations for different frequencies can be performed in parallel. Most importantly, it is ideally suited for implementation of the density matrix perturbation theory considered in the present paper (see Eq. (15)).

A small but finite imaginary frequency δ is added to the frequency Ω₀ of the external perturbation to avoid divergences at resonances^{38,39,40,41,43,44} so that Ω = Ω₀ + iδ. This imaginary frequency δ determines the linewidth in the calculated spectra.

To find the derivatives to the density matrix within the occupied, $\tilde \rho _{{\mathrm{VV}}}^{(P)}$, and unoccupied, $\tilde \rho _{{\mathrm{CC}}}^{(P)}$, subspaces, one can, in principle, also look for solution of the Liouville Eq. (12). However, in the case when the density matrix is idempotent, like the Kohn–Sham density matrix, the solution can be found explicitly from the idempotency condition, ρ = ρρ, and this reduces considerably the computational cost. The idempotency condition in terms of the periodic counterpart $\tilde \rho$ of the density matrix and to the first order in the magnetic field can be written as^27,29

$$\tilde \rho = \tilde \rho \tilde \rho + \frac{i}{{2c}}{\bf{B}} \cdot [{\bf{r}},\tilde \rho ] \times [{\bf{r}},\tilde \rho ].$$

(18)

The commutator $[{\bf{r}},\tilde \rho ]$ corresponding to $i\partial _{\bf{k}}\tilde \rho _{\bf{k}}$ in reciprocal space is determined in the present paper within the k·p theory^38,40,41 (see equations on pages 2 and 3 Supplementary information).

The polarizability α_0νμ in the absence of the magnetic field and the contribution α_νμ,γ to the polarizability in the presence of the magnetic field (α_νμ = α_0νμ + α_νμ,γB_γ) are obtained from the current response as

$$\begin{array}{*{20}{l}} {\alpha _{0\nu \mu }({\mathrm{\Omega }}) = \frac{i}{{\mathrm{\Omega }}}{\mathrm{Tr}}\left[ {V_\nu \tilde \rho ^{(E_\mu )}({\mathrm{\Omega }})} \right]} \hfill \end{array}$$

(19)

and

$$\begin{array}{*{20}{l}} {\alpha _{\nu \mu ,\gamma }({\mathrm{\Omega }}) = \frac{i}{{\mathrm{\Omega }}}{\mathrm{Tr}}\left[ {V_\nu \tilde \rho ^{\left( {E_\mu B_\gamma } \right)}({\mathrm{\Omega }})} \right].} \hfill \end{array}$$

(20)

These polarizabilies can be used to compute the experimentally measurable physical properties as described below.

Experimentally measured properties

The capacity of the system to absorb light is characterized using absorbance A = −log(I/I₀), which is defined through the ratio of intensities of the incident, I₀, and transmitted light, I. The magnitudes of the electric field vectors in the transmitted, E, and incident light, E₀, at frequency Ω₀ are related as E = E₀ exp(−n′Ω₀l/c), where n′ is the imaginary part of the refractive index n′ = Im n and l is the distance passed by the light through the sample studied. Since I ~ E², it can be stated that

$$\begin{array}{*{20}{l}} {A = \frac{{2n^\prime {\mathrm{\Omega }}_0l}}{{c\,{\mathrm{ln}}\,10}}.} \hfill \end{array}$$

(21)

The difference in the absorbance of the left (+) and right (−) circularly polarized light corresponds to the MCD response and is determined by the difference in the refractive indices $n_ + ^\prime - n_ - ^\prime$ for these two light components:

$$\begin{array}{*{20}{l}} {{\mathrm{\Delta }}A = A_ + - A_ - = \frac{{2\left( {n_{+} ^{\prime} -\, n_{-} ^{\prime} } \right){\mathrm{\Omega }}_0l}}{{c\,{\mathrm{ln}}\,10}}.} \hfill \end{array}$$

(22)

The refractive index n is determined by the equation

$$\begin{array}{*{20}{l}} {\epsilon _{\nu \mu }E_\mu = n^2E_\mu ,} \hfill \end{array}$$

(23)

where $\epsilon _{\nu \mu }$ is the dielectric tensor. For crystals, the dielectric tensor is related to the electric susceptibility χ_νμ as

$$\begin{array}{*{20}{l}} {\epsilon _{\nu \mu } = \delta _{\nu \mu } + 4\pi \chi _{\nu \mu }.} \hfill \end{array}$$

(24)

The latter corresponds to the polarizability per unit volume so that χ_νμ = α_νμ/w, where w is the unit cell volume and α_νμ is given by Eqs. (19) and (20).

In the case when the light propagation takes place along the optical axis z and no birefingence is observed, the refractive index in the absence of the magnetic field is equal to $n_0 = \epsilon _{xx}^{1/2} = \epsilon _{yy}^{1/2}$. The magnetic field provides just a small correction to this refractive index and it can be shown from Eqs. (23) and (24) (see pages 5 and 6 of Supplementary information) that

$$\begin{array}{*{20}{l}} {n^ \pm - n_0 \approx \pm i\frac{{2\pi \chi _{xy}}}{{n_0}}.} \hfill \end{array}$$

(25)

Using Eq. (22), the difference in the absorbance of the left and right circularly polarized light can be found as

$$\begin{array}{*{20}{l}} {{\mathrm{\Delta }}A_z = \frac{{4\pi \Omega _0l}}{{c\,{\mathrm{ln}}\,10}}{\mathrm{Re}}\left[ {\frac{{\chi _{xy}({\mathrm{\Omega }}) - \chi _{yx}({\mathrm{\Omega }})}}{{n_0}}} \right].} \hfill \end{array}$$

(26)

Note that ellipticity θ = (E₊ − E₋)/(E₊ + E₋) gained by the linearly polarized light is different just by a numerical coefficient θ = ΔA_z (ln 10)/4. The angle of Faraday rotation is determined by a similar expression as θ but with the imaginary part of χ_νμ instead of the real one¹ (see page 5 of Supplementary information). In contrast, in the magneto-optical polar Kerr effect for reflected light, the ellipticity and angle of rotation are determined by Im χ_νμ and Re χ_νμ, respectively.²

For molecules, the measurements are usually performed for a small concentration of randomly oriented molecules immersed into a transparent solvent or in vacuum. In this case, the total dielectric tensor of the medium can be presented as

$$\begin{array}{*{20}{l}} {\epsilon _{\nu \mu } = n_S^2\delta _{\nu \mu } + 4\pi \bar \alpha _{\nu \mu }N,} \hfill \end{array}$$

(27)

where δ_νμ is the Kronecker delta, n_S is the refractive index of the solvent or vacuum, $\bar \alpha _{\nu \mu }$ is the orientationally averaged polarizabiltiy of the molecules and N is their number density. The orientationally averaged polarizability is given by

$$\begin{array}{*{20}{l}} {\bar \alpha _{\nu \mu } = \frac{1}{3}\alpha _{0aa}\delta _{\nu \mu } + \frac{1}{6}Be_{abc}\alpha _{ab,c}e_{\nu \mu },} \hfill \end{array}$$

(28)

where e_νμ and e_abc are the Levi-Civita tensors of the second and third order, respectively, and the polarizabilities α_0aa and α_ab,c are computed from Eqs. (19) and (20) considering internal molecular axes.

For molecules, it is common to use molar extinction coefficients $\epsilon = A{\mathrm{/}}Cl$, i.e., absorbance per unit length and molar concentration. The molar concentration C in this expression is related to the number density as C = N/N_A, where N_A is the Avogadro constant. Taking into account that the concentration of the molecules is small, the refractive index in the absence of the magnetic field becomes approximately n₀ ≈ n_S + 2πNα_0aa/(3n_S) and this gives the molar extinction coefficient

$$\begin{array}{*{20}{l}} {\epsilon = \frac{{4\pi \Omega _0N_{\mathrm{A}}}}{{3n_Sc\,{\mathrm{ln}}\,10}}{\mathrm{Im}}\,\alpha _{0aa}.} \hfill \end{array}$$

(29)

The refractive indices for the left and right circularly polarized light can be correspondingly expressed as

$$\begin{array}{*{20}{l}} {n^ \pm - n_0 \approx \pm i\frac{{\pi NB}}{{3n_S}}e_{abc}\alpha _{ab,c}.} \hfill \end{array}$$

(30)

The difference $\Delta \epsilon$ in the molar extinction coefficients for the left and right circularly polarized light per unit magnetic field can, therefore, be found as

$$\begin{array}{*{20}{l}} {{\mathrm{\Delta }}\epsilon = \frac{{{\mathrm{\Delta }}A}}{{BCl}} = \frac{{4\pi N_{\mathrm{A}}}}{{3n_Sc\,{\mathrm{ln}}\,10}}e_{abc}{\mathrm{Re}}\,\alpha _{ab,c}.} \hfill \end{array}$$

(31)

The formalism for calculation of the magneto-optical response proposed in the present paper and expressions for the physical properties listed above have been implemented in the Octopus code.^36,37,38 The results of the tests for molecules and solids are presented below.

Results of calculations for molecules

First the tests of the developed formalism were performed for molecules (Fig. 1) in a large simulation box with periodic boundary conditions. Traditionally, the MCD response of molecules is divided into ${\cal A}$ and ${\cal B}$ terms (see equations on pages 4 and 5 of Supplementary information). The ${\cal B}$ term^12,14,15 comes from perturbations of molecular states in the magnetic field and is present in all systems. The ${\cal A}$ term^12,14,16 comes from perturbations of energies of excited states with non-zero orbital angular momenta. Such states are present only in molecules with rotational symmetry at least of the third order. As transitions to states with opposite orbital angular momenta are coupled to the light of different polarization, Zeeman splitting leads to an energy shift between absorption peaks for the left and right circularly polarized light. The MCD response in this case is described by the derivative of the spectral density^15,16 and has second-order poles.

Fig. 1

Molar extinction coefficient $\epsilon$ (a, c, in M⁻¹ cm⁻¹) and difference ${\mathrm{\Delta }}\epsilon$ in the molar extinction coefficients for the left and right circularly polarized light per unit magnetic field (b, d, in M⁻¹ cm⁻¹ T⁻¹) for adenine (a, b, $\delta = 0.05\,{\mathrm{eV}}$) and cyclopropane (c and d, $\delta = 0.02\,{\mathrm{eV}}$) as functions of the frequency of light Ω₀ (in eV) calculated using the present solid-state formalism (solid blue lines) and standard finite-system formulation (red dashed lines). The corresponding curves are virtually indistinguishable. The results obtained in the finite-system formulation for linewidths $\delta = 0.1\,{\mathrm{eV}}$ and $\delta = 0.2\,{\mathrm{eV}}$ are shown by magenta dash-dotted lines and green-dotted lines, respectively. The experimental data for adenine⁴⁶ in water and cyclopropane⁴⁷ in the gas phase are represented by circles. To show the results for different linewidths and experimental data on the same scale, the following scaling factors are introduced: 2, 1, and 1/2 for the linewidths of 0.2, 0.1, and 0.05 eV for the absorption and MCD spectra of adenine, 1, 1/2, and 1/10 for the linewidths of 0.2, 0.1, and 0.02 eV for the absorption spectra of cyclopropane and 1, 1/2, and 1/20, respectively, for the MCD spectra of cyclopropane. In the calculations for adenine, the refractive index of water is taken equal to 1.35 (ref. ⁶⁵). The parts of the spectra shown lie below the ionization potential at zero temperature (6.7 and 9.4 eV for adenine and cyclopropane, respectively, according to our calculations). Carbon, hydrogen, and nitrogen atoms in the atomistic structures are colored in gray, white, and blue, respectively. The inset shows the first MCD peak of cyclopropane

To check that both ${\cal A}$ and ${\cal B}$ terms are well described within the developed formalism, we have performed the calculations for adenine and cyclopropane (Fig. 1). Adenine is not symmetric and only the ${\cal B}$ term contributes to the magneto-optical response. Although we use the simplest local-density approximation (LDA)⁴⁵ for the exchange-correlation contribution to the electron energy and adiabatic approximation (ALDA) for the response, we find that the changes in the sign of the MCD signal for adenine are properly described as compared to the experimental data⁴⁶ (Fig. 1b). The magnitudes of the peaks for the simple optical absorption and the ${\cal B}$ term of the magneto-optical response scale inversely proportional to the linewidth, which is an input parameter of our calculations. Using a reasonable linewidth of δ = 0.1 eV, we get the absorption (Fig. 1a) and MCD (Fig. 1b) spectra with the magnitude of the peaks comparable to the experimental ones.

Cyclopropane has a rotational symmetry of the third order and its magneto-optical response has both ${\cal A}$ and ${\cal B}$ contributions. We find that the ${\cal A}$ term is clearly dominant for cyclopropane at linewidth δ = 0.02 eV (Fig. 1d), in agreement with previous calculations¹². However, the ${\cal A}$ and ${\cal B}$ terms scale differently with the linewidth. ${\cal B}$ term is inversely proportional to the linewidth, while the ${\cal A}$ is inversely proportional to square of the linewidth. Therefore, raising the linewidth to the experimental values of δ = 0.1–0.2 eV decreases the ${\cal A}$ term relative to the ${\cal B}$ term. For these linewidths, the shapes of the calculated curves and the magnitudes of the peaks approach the experimental ones⁴⁷ (Fig. 1c, d).

The calculations for the molecules (Fig. 1) demonstrate that the present formalism gives the results indistinguishable from the formulation using the position operator r (see page 4 of Supplementary information), which is commonly applied in literature for finite systems.^{12,13,14,15,16}

Results of calculations for solids

To test the developed formalism for solids we have applied it to bulk silicon and a monolayer of hexagonal boron nitride. For these periodic systems, we set the linewidth at δ = 0.1 eV, which is sufficient to resolve the important features of the spectra. As we use LDA for our test calculations, the excitation energies are systematically underestimated. To adjust the position of the peaks we apply the scissor operator, i.e., rigidly shift the spectra, to include the correction to the band gap known from GW calculations.^48,49,50 It should be, nevertheless, emphasized that the same code can be used with more advanced functionals like hybrid ones, which provide an improved description of the excitation energies. The approach can be also straightforwardly translated into the many-body framework.

Although account of local-field effects through Eq. (13) even within the simplest ALDA approximation is very important for molecules, for silicon and boron nitride, such adiabatic effects provide a minor correction to the spectra (see Fig. 2 of Supplementary information). The account of long-range exchange and correlation interactions in solids is, on the other hand, crucial for description of excitons. To take them into account we follow the approach proposed in ref. ³⁵ in the TDCDFT framework. In this approach, non-adiabatic local-field effects are introduced through the exchange-correlation electric field

$$\begin{array}{*{20}{l}} {{\bf{E}}_{{\mathrm{mac}}}^{{\mathrm{xc}}}({\mathrm{\Omega }}) = \frac{{i{\mathrm{\Omega }}}}{w}{\int}_w {\kern 1pt} {\mathrm{d}}{\bf{r}}{\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}^\prime \hat f_{{\mathrm{xc}}}({\bf{r}},{\bf{r}}^\prime ,{\mathrm{\Omega }})\delta {\bf{j}}({\bf{r}}^\prime ,{\mathrm{\Omega }}),} \hfill \end{array}$$

(32)

where tensor $\hat f_{{\mathrm{xc}}}({\bf{r}},{\bf{r}}^\prime ,{\mathrm{\Omega }})$ is the TDCDFT exchange-correlation kernel and δj(r′, Ω) is the induced current density. This field together with the macroscopic electric field E_mac gives the macroscopic Kohn–Sham electric field ${\bf{E}}_{{\mathrm{mac}}}^{{\mathrm{KS}}} = {\bf{E}}_{{\mathrm{mac}}} + {\bf{E}}_{{\mathrm{mac}}}^{{\mathrm{xc}}}$.

The macroscopic polarization

$$\begin{array}{*{20}{l}} {{\bf{P}}_{{\mathrm{mac}}}({\mathrm{\Omega }}) = \frac{{ - i}}{{{\mathrm{\Omega }}w}}{\int}_w {\kern 1pt} {\mathrm{d}}{\bf{r}}\delta {\bf{j}}(\vec r,{\mathrm{\Omega }}),} \hfill \end{array}$$

(33)

is related to the macroscopic Kohn–Sham electric field ${\bf{E}}_{{\mathrm{mac}}}^{{\mathrm{KS}}}$ through the Kohn–Sham electric susceptibility tensor $\hat \chi ^{{\mathrm{KS}}}$ and to the macroscopic electric field E_mac through the net susceptibility tensor $\hat \chi$:

$$\begin{array}{*{20}{l}} {{\bf{P}}_{{\mathrm{mac}}}({\mathrm{\Omega }})} \hfill & = \hfill & {\hat \chi ^{{\mathrm{KS}}}({\mathrm{\Omega }})\left( {{\bf{E}}_{{\mathrm{mac}}}({\mathrm{\Omega }}) + {\bf{E}}_{{\mathrm{mac}}}^{{\mathrm{xc}}}({\mathrm{\Omega }})} \right)} \hfill \\ {} \hfill & = \hfill & {\hat \chi ({\mathrm{\Omega }}){\bf{E}}_{{\mathrm{mac}}}({\mathrm{\Omega }})} \hfill \end{array}$$

(34)

Neglecting microscopic current components in Eq. (32), i.e., replacing the induced current density $\delta {\bf{j}}(\vec r\prime ,{\mathrm{\Omega }})$ by its unit cell average, and using Eq. (33), the exchange-correlation electric field is written as

$$\begin{array}{*{20}{l}} {{\bf{E}}_{{\mathrm{mac}}}^{{\mathrm{xc}}}({\mathrm{\Omega }}) = \hat \beta ({\mathrm{\Omega }}){\bf{P}}_{{\mathrm{mac}}}({\mathrm{\Omega }}),} \hfill \end{array}$$

(35)

where

$$\begin{array}{*{20}{l}} {\hat \beta ({\mathrm{\Omega }}) = - \frac{{{\mathrm{\Omega }}^2}}{w}{\int}_w {\kern 1pt} {\mathrm{d}}{\bf{r}}{\int} {\kern 1pt} {\mathrm{d}}{\bf{r}}^\prime \hat f_{{\mathrm{xc}}}({\bf{r}},{\bf{r}}^\prime ,{\mathrm{\Omega }}).} \hfill \end{array}$$

(36)

Substitution of Eq. (35) into Eq. (34) gives

$$\begin{array}{*{20}{l}} {\frac{1}{{\hat \chi ({\mathrm{\Omega }})}} = \frac{1}{{\hat \chi ^{{\mathrm{KS}}}({\mathrm{\Omega }})}} - \hat \beta ({\mathrm{\Omega }}).} \hfill \end{array}$$

(37)

In the simplest case, $\hat \beta$ can be assumed static and isotropic, i.e. β_νμ = βδ_νμ. Then the longitudial and transverse components of the electric susceptibility tensor are given by

$$\begin{array}{*{20}{l}} {\chi _{xx}({\mathrm{\Omega }}) = \frac{{\chi _{xx}^{{\mathrm{KS}}}}}{{1 - \beta \chi _{xx}^{{\mathrm{KS}}}({\mathrm{\Omega }})}}} \hfill \end{array}$$

(38)

and

$$\begin{array}{*{20}{l}} {\chi _{xy}({\mathrm{\Omega }}) \approx \frac{{\chi _{xy}^{{\mathrm{KS}}}}}{{\left( {1 - \beta \chi _{xx}^{{\mathrm{KS}}}({\mathrm{\Omega }})} \right)\left( {1 - \beta \chi _{yy}^{{\mathrm{KS}}}({\mathrm{\Omega }})} \right)}},} \hfill \end{array}$$

(39)

respectively. In these expressions, we neglect the terms of the second order in the transverse components of $\hat \chi ^{{\mathrm{KS}}}$.

It should be noted that Eq. (38) for the longitudinal response is equivalent to the head term of the long-range contribution (LRC) to the exchange-correlation kernel^48,49,51 in TDDFT, which corresponds to $f_{xc}^{{\mathrm{(LRC)}}}({\bf{q}}) = - \beta {\mathrm{/}}q^2$ in reciprocal space. However, the latter model does not describe properly the transverse response. Eq. (39) gives an adequate expression for the transverse response thanks to the tensorial nature of the exchange-correlation kernel $\hat f_{{\mathrm{xc}}}({\bf{r}},{\bf{r}}^\prime ,{\mathrm{\Omega }})$ in the TDCDFT framework.

Let us first discuss the results for bulk silicon (Fig. 2). Figure 2b shows that the spectra ${\mathrm{Re/Im}}\,\epsilon _{xy}$ for the transverse component of the dielectric tensor calculated even without account of excitonic effects follow qualitatively the shapes of the experimental curves² at the direct absorption edge. The analysis of optical transitions at the Γ point of the Brillouin zone, where the highest valence and lowest conduction bands are formed by triply degenerate p-like states (${\mathrm{\Gamma }}_{25}^\prime$ and Γ₁₅, respectively),⁵² reveals significant contributions that can be attributed to the ${\cal A}$ term (Fig. 2d). Two inequivalent contributions come from excitations with the change in the magnetic quantum number l_z from 0 to ±1 and vice versa. The ratio $\epsilon _{xy}{\mathrm{/}}\epsilon _{xx}$ for each of them at the resonance frequency characterizes the relative frequency shift in the magnetic field

$$\begin{array}{*{20}{l}} {\frac{{\epsilon _{xy}}}{{\epsilon _{xx}}}\sim \frac{{{\mathrm{\Delta }}m_zB_z{\mathrm{\Delta }}l_z}}{\delta },} \hfill \end{array}$$

(40)

where Δm_z is the change in the orbital magnetic dipole moment and Δl_z is the change of the magnetic quantum number (see explanation on page 5 of Supplementary information). Correspondingly, we can estimate the effective g-factors g = −Δm_z/μ_BΔl_z, where μ_B is the Bohr magneton, and they are found to be g = 3.5 in ${\mathrm{\Gamma }}_{25}^\prime \to {\mathrm{\Gamma }}_{15}$ transitons with l_z = 0 → ±1 and g = −0.40 for l_z = ±1 → 0. Note that nearly the same values are obtained using explicit expressions for the band magnetic dipole moments from refs ^23,53 (see page 4 of Supplementary information). Thus, unlike absorption, transitions l_z = 0 → ±1 prevail in the magneto-optical response at the band edge. The domination of the ${\cal A}$ term is consistent with the experiments, where ${\mathrm{Re/Im}}\,\epsilon _{xy}$ (Fig. 2b) look similar to derivatives of ${\mathrm{Im/Re}}\,\epsilon _{xx}$ (Fig. 2a).

Fig. 2

Calculated components $\epsilon _{xx}$ a and $\epsilon _{xy}$ b of the dielectric tensor of silicon as functions of the frequency of light Ω₀ (in eV) for the magnetic field of 1 T along the z axis. The real and imaginary parts are shown by solid and dashed lines, respectively. The results obtained with and without account of excitonic effects correspond to red and black lines, respectively. The calculated data are blue-shifted in energy by 0.7 eV to take into account the GW correction to the band gap.^48,49 The experimental data from refs ^2,54 for $\epsilon _{xx}$ and $\epsilon _{xy}$ are shown by symbols. The experimental data for $\epsilon _{xy}$ are scaled by a factor of 1/2. Squares correspond to the real parts and circles to the imaginary ones. The transitions at the Γ point of the Brillouin zone are indicated by the vertical gray line. Calculated contributions to ${\mathrm{Im}}\,\epsilon _{xx}$ c and ${\mathrm{Im}}\,\epsilon _{xy}$ d from the Γ point: total contribution (triangles), contribution from all transitions ${\mathrm{\Gamma }}_{25}^\prime \to {\mathrm{\Gamma }}_{15}$ to the ${\cal A}$ term (blue dashed lines) and contributions from transitions ${\mathrm{\Gamma }}_{25}^\prime \to {\mathrm{\Gamma }}_{15}$ with the magnetic quantum number $l_z = 0 \to \pm 1$ (green solid lines) and $\pm 1 \to 0$ (black dash-dotted lines) to the ${\cal A}$ term. Total ${\mathrm{Re}}\,\epsilon _{xy}$ is shown by diamonds

To model excitonic effects in silicon, we use Eq. (37) with β = 0.2. This value fulfils the empirical law $\beta = 4.615{\mathrm{/}}\epsilon _\infty - 0.213$, where $\epsilon _\infty$ is the static dielectric constant, derived for a set of semiconductors with continuum excitons.^48,49 The account of the excitonic effects further improves agreement of the calculated spectra for silicon with the experimental data (Fig. 2a, b).

It should be noted, however, that although the magnitudes of peaks in the longitudinal component $\epsilon _{xx}$ of the dielectric tensor agree very well with the experimental results,⁵⁴ the magnitudes of the peaks in the transverse component $\epsilon _{xy}$ are about a factor of two smaller than in the magneto-optical measurements.² As discussed above for molecules, the magnitudes of peaks in magneto-optical calculations are strongly dependent on the linewidth assumed. The ratio of the magnitudes of peaks coming from the ${\cal A}$ term and those corresponding to the simple absorption scale inversely proportional to the linewidth (see Eq. (40)). Therefore, agreement with the experimental magneto-optical spectra should be improved once the linewidth in the calculations is reduced. Fine-tuning of the linewidth is, however, beyond the scope of the present paper.

In boron nitride (Fig. 3), the magneto-optical response of continuum states starts from a prominent peak at the band edge (Fig. 3b). In this material, the first optical transitions take place at the K^± points in the corners of the hexagonal Brillouin zone, where phase winding of wavefunctions related to the C₃ symmetry imposes coupling to only one light component of the left (+) or right (−) circular polarization.^55,56,57 Accordingly, contributions to the magneto-optical spectra from the K^± points can be described by a second-order pole (Fig. 3c). The map of contributions from different k-points (Fig. 3d) shows that the response is mostly provided by narrow regions in reciprocal space and the sign of the response is opposite in two such regions. Therefore, it can be concluded that the ${\cal A}$ term is dominant at the band edge of boron nitride.

Fig. 3

Calculated components $\epsilon _{xx}$ a and $\epsilon _{xy}$ b of the dielectric tensor of boron nitride monolayer as functions of the frequency of light Ω₀ (in eV) for the magnetic field of 1 T along the z axis directed out of the plane. The real and imaginary parts are represented by solid and dashed lines, respectively. The results obtained with and without account of excitonic effects correspond to red and black lines, respectively. The data for $\epsilon _{xx}$ and $\epsilon _{xy}$ obtained with account of excitonic effects are multiplied by 1/10 and 1/50, respectively, to show all the results on the same scale. The calculated data are blue-shifted in energy by 2.6 eV to take into account the GW correction to the band gap⁵⁰. The transitions at the K and M points of the Brillouin zone are indicated by vertical gray lines. Boron and nitrogen atoms in the atomistic structure are colored in magenta and blue, respectively. c Calculated contributions to $\epsilon _{xx}$ (black lines) and $\epsilon _{xy}$ (blue lines) from the K points of the Brillouin zone. d Calculated contributions to ${\mathrm{Re}}\,\epsilon _{xy} \times 10^3$ from different points (k_x, k_y, 0) (in Å⁻¹) of the Brillouin zone of the 4-atom cell for Ω₀ = 7.8 eV

Clearly such a magneto-optical response is related to the valley Zeeman effect.^4,5,6,7,8 As the density of states in two-dimensional materials tends to the Heaviside step function in the limit of zero linewidth, the ${\cal A}$ term related to its derivative approaches a delta peak. Thus, discrete peaks in continuum magneto-optical spectra of two-dimensional materials are indicators of the Zeeman splitting.

From the comparison of magneto-optical and optical spectra for boron nitride, we estimate that the change of the magnetic dipole moment upon the excitation at the K^± points is ${\mathrm{\Delta }}m_z^ \pm \approx \mp 1.8\mu _{\mathrm{B}}$. Explicit calculations of the magnetic dipole moments using expressions from refs ^23,53 give ∓0.95μ_B and ∓2.8μ_B for the valence and conduction bands, respectively, which agrees very well with our estimate. The valley g-factor for the edge of the continuum spectrum according to our calculations is, therefore, $g^{{\mathrm{vl}}} = - 2{\mathrm{\Delta }}m_z^ + {\mathrm{/}}\mu _{\mathrm{B}} = 3.6$.

Up to now we have neglected excitonic effects in boron nitride. They, however, are known to be very strong.⁵⁰ To describe the first bound exciton in boron nitride we set the parameter β in Eq. (37) at β = 17.5 to reproduce the binding energy of 1.4 eV that follows from the Bethe–Salpeter calculations⁵⁰ (Fig. 3a). The absorption (Fig. 3a) and magneto-optical (Fig. 3b) spectra computed using this parameter are very similar to those of symmetric molecules like cyclopropane (Fig. 1c, d). The valley g-factor deduced from the ratio ${\mathrm{Im}}\,\epsilon _{xy}{\mathrm{/Im}}\,\epsilon _{xx}$ at the excitonic peak is ~1.8. It is, therefore, reduced twice compared to the result for the edge of the continuum spectrum. To confirm our estimate, a photoluminescence experiment for boron nitride could be performed by analogy with the measurements for WSe₂ (refs ^4,5,6) and MoSe₂ (refs ^6,7,8) monolayers (see page 7 of Supplementary information for discussion of g-factors observed for these materials). It should be noted that the qualitative shapes of the spectra computed with account of the excitonic effects do not depend on the parameter β used (see Fig. 3 of Supplementary information) and the valley g-factor changes only by 30% in the interval of β from 10 to 20.

To summarize, in spite of simplifications made in the present paper for the test calculations, the developed formalism gives realistic results for the magneto-optical response. It provides a unified description of finite and periodic systems and automatically takes into account gauge invariance. Furthermore, it can be straightforwardly extended to the case of higher-order responses to arbitrary electromagnetic fields.

The efficiency of the implemented procedures for magneto-optics is comparable to standard linear-response calculations of polarizability in the absence of the magnetic field. When local-field effects are included self-consistently, the calculations of magneto-optical spectra for molecules take the same time as polarizability. For solids, the responses at ±Ω₀ ± iδ are needed for magneto-optics as compared only to ±Ω₀ + iδ for simple optics (see the detailed explanation on pages 7 and 8 of Supplementary information) and, therefore, the calculations of magneto-optical spectra take twice as long as those of polarizability.

Methods

The interaction of valence electrons with atomic cores is described using Troullier-Martins norm-conserving pseudopotentials.⁵⁸ For molecules, the density-averaged self-interaction correction⁵⁹ is applied to avoid spurious transitions to diffuse excited states. The efficient conjugate-gradients solver⁶⁰ is used for the calculation of eigenstates with the tolerance of 10⁻¹⁰ and mixing parameter for the Kohn–Sham potential of 0.2 for molecules and 0.1 for solids. The semiconducting smearing is applied. The magnetic gauge correction from ref. ⁶¹ is added in calculations of magneto-optical spectra of the molecules within the finite-system formulation. The quasi-minimal residual (QMR) method⁶² (qmr_symmetric and qmr_dotp for the molecules and solids, respectively) with the final tolerance of 10⁻⁶ is used to solve linear equations for projections of derivatives of the density matrix onto unperturbed wavefunctions (Eq. (15)). The local-field effects in the ALDA approximation are taken into account through a self-consistent iteration scheme similar to the ground-state DFT.

For molecules, the size of the simulation box of 24 Å and the spacing of the real-space grid of 0.14 Å are sufficient for convergence of the magneto-optical spectra. Only the Γ point is used in this case. The geometry of the molecules is optimized till the maximal residual force of 0.01 eV/Å using the fast inertial relaxation engine (FIRE) algorithm.⁶³ For boron nitride, we consider the rectangular unit cell of 4.294 Å × 2.479 Å × 24.0 Å with four atoms. For silicon, the cubic unit cell of 5.38 Å size with 8 atoms is studied and the grid spacing is increased to 0.25 Å. Integration over the Brillouin zone is performed according to the Monkhorst-Pack method.⁶⁴ Time-reversal and crystal symmetries are taken into account to reduce the number of k-points considered. To take into account time-reversal symmetry, the average of the polarizabilities at frequencies Ω and −Ω is computed for irreducible k-points. Three thousand irreducible k-points are needed for convergence of the magneto-optical spectra for boron nitride and 6600 for silicon and these are achieved using shifted k-point grids (see the results of calculations using different k-point grids in Figs. 1 and 2 of Supplementary information).

Code availability

Our implementation is available through the development version of the Octopus code at https://gitlab.com/octopus-code/octopus.git and will be available in future releases at https://octopus-code.org. The code is provided under the GNU General Public License. The manual and tutorials can be found at https://octopus-code.org.