Orbital magneto-optical response of periodic insulators from first principles

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 \textit{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 magnetooptical disks and optical isolators 3 . 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 10 and metal oxide nanosheets 11 . These advances cultivate the growing interest to development of a gauge-invariant and computationally efficient ab initio theory of magneto-optical response.
While 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. Though 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. While 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-gaugeinvariant 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. While 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. 35 from timedependent 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. While the spin contribution is trivial, the account of the Fermi surface contribution can be done for metals by analogy with Ref. 23.
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.

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 −1 ∂ t A, 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 O = O r1r2 defined for two points r 1 and r 2 in real space it is possible to distinguish the periodic and gauge-invariant counter-partÕ =Õ r1r2 by factoring out the Aharonov-Bohmtype phase [27][28][29] Here we take = e = 1 and the integral is taken along the straight line between points r 2 and r 1 so that r = r 2 + (r 1 − r 2 )ξ, 0 ≤ ξ ≤ 1. This approach was previously used to derive corrections to the gauge-invariant counterpartρ of the oneparticle 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 in terms ofρ. Here and below the commutator of two operators O (1) and O (2) is introduced as Using Eq. (2) for the relation betweenρ and ρ in real space, the time-dependent Liouville equation (3) gives It should be noted thatH = H 0 + δH, where the difference δH between the gauge-invariant counterpartH of the Hamiltonian and unperturbed Hamiltonian H 0 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 exchangecorrelation potentials in TDDFT (see page 1 of Supplementary information). Eq. (5) is equivalent to where ϕ 123 = ϕ 12 + ϕ 23 + ϕ 31 . This phase corresponds to the flux of the magnetic field through the triangle formed by points r 1 , r 2 and r 3 : The time derivative of the phase ϕ 13 on the left-hand side of Eq. (6) introduces the electric field Combining Eqs. (6)-(8), we arrive at 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 c/a 2 ∼ 10 5 T, where a = 1Å is taken as a typical interatomic distance.
Eq. (9) for the density matrix then takes the form Using that O r1r2 (r 1 − r 2 ) = [r, O] r1r2 and introducing notations for the anticommutator of two operators O (1) and O (2) and velocity operator V = −i[r,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 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 [r,ρ] of the position operator with the periodic functionρ 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, we get all terms dependent on the external fields on the right-hand side of the equation. Differentiating the Liouville equation (13), one can evaluate the derivatives of the density matrixρ (P ) = ∂ρ/∂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 equation (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 firstprinciples framework and a similar computational scheme can be used.
From the computational point of view, it is convenient to divide the n-th order derivativeρ (P ) of the density matrix describing the joint response to the perturbations P = P 1 P 2 ...P n into four blocks within and between the occupied (V) and unoccupied subspaces (C): These blocks correspond toρ and P c = 1 − P v are the projectors onto the occupied and unoccupied bands.
Following the density matrix perturbation theory 32 , to get the elements of the derivative of the density matrix ρ Here the operator on the left-hand side is given by where Ω is frequency considered and vk is the energy of the unperturbed state |ψ vk . 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 1 ) = ρ (P ) (r 1 , r 2 )δ(r 1 − r 2 ) (see page 1 of Supplementary information).
The solution of Eq. (15) corresponds to and once it is known, the elementsρ CV of the derivative of the density matrix between unoccupied and occupied subspaces can be computed as The elements between the occupied and unoccupied subspaces can be found asρ CV (−Ω * )) * and to obtain them, Eq. (15) should be also solved for the frequency −Ω * . If the local-field effects are taken into account, Eq. (15) has to be solved self-consistently as the derivativeρ (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 |η (P ) vk (Ω) that fits into this equation at each frequency Ω. Other approaches, such as sum over states 15 , methods based on Casida's equation 15,16 , complex polarization propagator 12,13 and real-time propagation 14 have been used to compute absorption and magneto-optical spectra of molecules. The sum over states, Casida's equation 42 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 42 furthermore relies on the use of real wavefunctions and cannot be straightforwardly extended to solids, where KS states are complex.
Neither Sternheimer approach 38-41 , nor real-time propagation 14 need calculation of unoccupied states. They also have a favourable scaling of O(N 2 ) with the system size N as compared, for example, to O(N 3 ) for the sum over states (Refs. 14, 39, and 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)).
To find the derivatives to the density matrix within the occupied,ρ (P ) VV , and unoccupied,ρ (P ) CC , subspaces, one can, in principle, also look for solution of the Liouville equation (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ρ of the density matrix and to the first order in the magnetic field can be written as 27,29ρ =ρρ The commutator [r,ρ] corresponding to i∂ kρ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 and 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 0 ), which is defined through the ratio of intensities of the incident, I 0 , and transmitted light, I. The magnitudes of the electric field vectors in the transmitted, E, and incident light, E 0 , at frequency Ω 0 are related as E = E 0 exp(−n Ω 0 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 2 , it can be stated that 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 + −n − for these two light components: The refractive index n is determined by the equation where νµ is the dielectric tensor. For crystals, the dielectric tensor is related to the electric susceptibility χ νµ as 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 = Using Eq. 22, the difference in the absorbance of the left and right circularly polarized light can be found as 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 1 (see page 5 of Supplementary information). By contrast, in the magneto-optical polar Kerr effect for reflected light, the ellipticity and angle of rotation are determined by Im χ νµ and Re χ νµ , respectively 2 .
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 where δ νµ is the Kronecker delta, n S is the refractive index of the solvent or vacuum,ᾱ νµ is the orientationally averaged polarizabiltiy of the molecules and N is their number density. The orientationally averaged polarizability is given bȳ 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 = A/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 0 ≈ n S + 2πN α 0aa /(3n S ) and this gives the molar extinction coefficient The refractive indices for the left and right circularly polarized light can be correspondingly expressed as The difference ∆ in the molar extinction coefficients for the left and right circularly polarized light per unit magnetic field can, therefore, be found as 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 A and B terms (see equations on pages 4 and 5 of Supplementary information). The B term 12,14,15 comes from perturbations of molecular states in the magnetic field and is present in all systems. The 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. Since 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 secondorder poles.
To check that both A and 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 B term contributes to the magneto-optical response. Though we use the simplest local-density approximation (LDA) 48 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 45 (Fig. 1b). The magnitudes of the peaks for the simple optical absorption and the 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 A and B contributions. We find that the A term is clearly dominant for cyclopropane at linewidth δ = 0.02 eV (Fig. 1d), in agreement with previous calculations 12 . However, the A and B terms scale differently with the linewidth. B term is inversely proportional to the linewidth, while the 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 A term relative to the B term. For these linewidths, the shapes of the calculated curves and the magnitudes of the peaks approach the experimental ones 46 (Figs. 1c and 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] .  47). The parts of the spectra shown lie below the ionization potential at zero temperature (6.7 eV and 9.4 eV for adenine and cyclopropane, respectively, according to our calculations). Carbon, hydrogen and nitrogen atoms in the atomistic structures are coloured in gray, white and blue, respectively. The inset shows the first MCD peak of cyclopropane.

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. Since 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 [49][50][51] . 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.
While 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. 35 in the TDCDFT framework. In this approach, non-adiabatic local-field effects are introduced through the exchange-correlation electric field E xc mac (Ω) = iΩ w w dr dr f xc (r, r , Ω)δj(r , Ω), (32) where tensorf xc (r, r , Ω) is the TDCDFT exchangecorrelation 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 E KS mac = E mac + E xc mac . The macroscopic polarization is related to the macroscopic Kohn-Sham electric field E KS mac through the Kohn-Sham electric susceptibility tensorχ KS and to the macroscopic electric field E mac through the net susceptibility tensorχ: Neglecting microscopic current components in Eq. (32), i.e. replacing the induced current density δj( r , Ω) by its unit cell average, and using Eq. (33), the exchangecorrelation electric field is written as In the simplest case,β can be assumed static and isotropic, i.e. β νµ = βδ νµ . Then the longitudial and transverse components of the electric susceptibility tensor are given by and respectively. In these expressions, we neglect the terms of the second order in the transverse components ofχ 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 exchangecorrelation kernel 49,50,53 in TDDFT, which corresponds to f (LRC) xc (q) = −β/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 kernelf xc (r, r , Ω) in the TDCDFT framework.
Let us first discuss the results for bulk silicon (Fig. 2). Fig. 2b shows that the spectra Re/Im xy for the transverse component of the dielectric tensor calculated even without account of excitonic effects follow qualitatively the shapes of the experimental curves 2 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 (Γ 25 and Γ 15 , respectively) 54 , reveals significant contributions that can be attributed to the 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 xy / xx for each of them at the resonance frequency characterizes the relative frequency shift in the magnetic field 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. To model excitonic effects in silicon we use Eq. (37) with β = 0.2. This value fulfils the empirical law β = 4.615/ ∞ − 0.213, where ∞ is the static dielectric constant, derived for a set of semiconductors with continuum excitons 49,50 . The account of the excitonic effects further improves agreement of the calculated spectra for silicon with the experimental data ( Fig. 2a and b).
It should be noted, however, that though the magnitudes of peaks in the longitudinal component xx of the dielectric tensor agree very well with the experimental results 52 , the magnitudes of the peaks in the transverse component xy are about a factor of two smaller than in the magneto-optical measurements 2 . 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 A term and those correspond-ing 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 3 symmetry imposes coupling to only one light component of the left (+) or right (−) circular polarization [56][57][58] . 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 A term is dominant at the band edge of boron nitride.
Clearly such a magneto-optical response is related to the valley Zeeman effect [4][5][6][7][8] . Since the density of states in two-dimensional materials tends to the Heaviside step function in the limit of zero linewidth, the 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 ∆m ± z ≈ ∓1.8µ B . Explicit calculations of the magnetic dipole moments using expressions from Refs. 23 and 55 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 vl = −2∆m + z /µ B = 3.6. Up to now we have neglected excitonic effects in boron nitride. They, however, are known to be very strong 51 . 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 51 (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 and d). The valley g-factor deduced from the ratio Im xy /Im xx at the excitonic peak is about 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 2 (Refs. 4-6) and MoSe 2 (Refs. 6-8) monolayers (see page 7 of Supplementary information for discussion of g-factors observed for these materials). It should be noted that the qualita-tive 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 linearresponse calculations of polarizability in the absence of the magnetic field. When local-field effects are included self-consistently, the calculations of magnetooptical spectra for molecules take the same time as polarizability. For solids, the responses at ±Ω 0 ± iδ are needed for magneto-optics as compared only to ±Ω 0 + 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 59 . For molecules, the density-averaged self-interaction correction 60 is applied to avoid spurious transitions to diffuse excited states. The efficient conjugate-gradients solver 61 is used for the calculation of eigenstates with the tolerance of 10 −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. 62 is added in calculations of magneto-optical spectra of the molecules within the finite-system formulation. The quasi-minimal residual (QMR) method 63 (qmr symmetric and qmr dotp for the molecules and solids, respectively) with the final tolerance of 10 −6 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 selfconsistent 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 64 . For boron nitride, we consider the rectangular unit cell of 4.294Å × 2.479Å × 24.0 A with four atoms. For silicon, the cubic unit cell of 5.38 A 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 65 . 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. 3000 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.