Let $(M,g)$ be a Riemannian manifold and let $k$ be a Killing vector of $g$, that is a vector field such that $\mathcal{L}_kg = 0$.
I would like a general method for finding coordinates 'adapted' to this Killing vector, that is, coordinates $(X^{\mu})=(X^i,\theta)$, such that $$k=\frac{\partial}{\partial \theta}$$
This is frequently referenced in physics papers, but I have no idea how to do this in principle or in practice. One often sees statements along the lines of "...by choosing coordinates adapted to the isometry...", but no description of how to find such coordinates.
A concrete example would be helpful, so consider the following Riemannian metric, given in terms of coordinates $(x,y,z)$ as: $$ds ^2 = dx^2 + (dy - xdz)^2 + dz^2$$ If I am interested in the Killing vector $k=\partial_x + z\partial_y$, how do I find coordinates adapted to this Killing vector? Note that the given coordinates are already adapted to the Killing vectors $\partial_y$ and $\partial_z$, but that's not what I'm asking.
A basic recipe: Choose a hypersurface $\Sigma$ transverse to $k$ (which will end up being $\{\theta = 0\}$) and a coordinate system $X^i$ for this hypersurface, and then transport $\Sigma$ along with $X^i$ using the flow of $k$.
In your example, an easy choice of $\Sigma$ is simply the coordinate surface $\{ x = 0\}$, which comes an obvious choice of coordinates $X^1 = y, X^2 = z$. Now we want to flow this along $k$; i.e. we want a map $(x,y,z) = \phi(X^1, X^2, \theta)$ such that $\phi(X^1, X^2, 0) = (0,X^1, X^2)$ and $$ \frac{\partial \phi}{\partial \theta} = \partial_x + z\partial_y.$$
Expanding this in components gives us $$\frac{\partial x}{\partial \theta} = 1, \; \frac{\partial y}{\partial \theta} = z, \; \frac{\partial z}{\partial \theta} = 0.$$
These equations along with the initial conditions $x = 0, y = X^1, z= X^2$ are easily integrated to yield $x=\theta, z = X^2, y = X^1 + \theta X^2$. Inverting these gives (an example of) the adapted coordinates you are looking for:
$$\theta = x, X^1 = y - xz, X^2 = z.$$