  We introduce new methods to numerically construct for the first time
stationary axisymmetric black hole solutions in Einstein-aether theory and
study their properties. The key technical challenge is to impose regularity at
the spin-2, 1, and 0 wave mode horizons. Interestingly we find the metric
horizon, and various wave mode horizons, are not Killing horizons, having null
generators to which no linear combination of Killing vectors is tangent, and
which spiral from pole to equator or vice versa. Existing phenomenological
constraints result in two regions of coupling parameters where the theory is
viable and some couplings are large; region I with a large twist coupling and
region II with also a (somewhat) large expansion coupling. Currently these
constraints do not include tests from strong field dynamics, such as
observations of black holes and their mergers. Given the large aether
coupling(s) one might expect such dynamics to deviate significantly from
general relativity, and hence to further constrain the theory. Here we argue
this is not the case, since for these parameter regions solutions exist where
the aether is "painted" onto a metric background that is very close to that of
general relativity. This painting for region I is approximately independent of
the large twist coupling, and for region II is also approximately independent
of the large expansion coupling and normal to a maximal foliation of the
spacetime. We support this picture analytically for weak fields, and
numerically for rotating black hole solutions, which closely approximate the
Kerr metric.
