  We extend the recently-developed theory of bulk orbital magnetization to
finite electric fields, and use it to calculate the orbital magnetoelectric
response of periodic insulators. Working in the independent-particle framework,
we find that the finite-field orbital magnetization can be written as a sum of
three gauge-invariant contributions, one of which has no counterpart at zero
field. The extra contribution is collinear with and explicitly dependent on the
electric field. The expression for the orbital magnetization is suitable for
first-principles implementations, allowing to calculate the magnetoelectric
response coefficients by numerical differentiation. Alternatively,
perturbation-theory techniques may be used, and for that purpose we derive an
expression directly for the linear magnetoelectric tensor by taking the first
field-derivative analytically. Two types of terms are obtained. One, the
`Chern-Simons' term, depends only on the unperturbed occupied orbitals and is
purely isotropic. The other, `Kubo' term, involves the first-order change in
the orbitals and gives isotropic as well as anisotropic contributions to the
response. In ordinary magnetoelectric insulators both terms are generally
present, while in strong Z2 topological insulators only the former type is
allowed, and is quantized. In order to validate the theory we have calculated
under periodic boundary conditions the linear magnetoelectric coupling for a
3-D tight-binding model of an ordinary magnetoelectric insulator, using both
the finite-field and perturbation-theory expressions. The results are in
excellent agreement with calculations on bounded samples.
