A current challenge in atomistic machine learning is that of efficiently predicting the response of the electron density under electric fields. We address this challenge with symmetry-adapted kernel functions that are specifically derived to account for the rotational symmetry of a three-dimensional vector field. We demonstrate the equivariance of the method on a set of rotated water molecules and show its high efficiency with respect to number of training configurations and features for liquid water and naphthalene crystals. We conclude showcasing applications for relaxed configurations of gold nanoparticles, reproducing the scaling law of the electronic polarizability with size, up to systems with more than 2000 atoms. By deriving a natural extension to equivariant learning models of the electron density, our method provides an accurate and inexpensive strategy to predict the electrostatic response of molecules and materials.