Orbital magnetoelectric response of insulating crystals

We calculate the orbital magnetoelectric polarizability (OMP) of a periodic insulator as the linear orbital magnetization response to a homogeneous electric field. We begin by considering the orbital magnetization (OM) in a finite field, and find that it can be written as a sum of three terms, one of which has no counterpart at zero field. The extra contribution is parallel to the electric field and is a multivalued quantity, only defined up to a field-dependent quantum. The full OM expression can be implemented in {\it ab-initio} codes, allowing to calculate the OMP by finite differences. Alternatively, linear-response techniques may be used, and for that purpose we obtain an expression directly for the OMP, by taking analytically the field derivative of the OM formula. The resulting expression has two types of terms: one which is expressed solely in terms of the unperturbed valence wavefunctions, and another coming from the first-order change in the wavefunctions. In normal insulators both are present, while in strong topological insulators only the former type is present. The full expression for the OMP tensor is verified numerically for a 3D tight-binding model of a normal insulator with broken inversion and time-reversal symmetries, by comparing with the finite-field orbital magnetization calculated for both finite and periodic samples.