Recipe to find null coordinates Given a metric $g_{\mu\nu}$, is there a general recipe to define the null coordinates for the system?
For the Schwarzschild metric, I think we can find it by the following method.
$$ds^2 = -f(r)dt^2 + \frac{1}{f(r)}dr^2$$
For a null coordinate $ds^2=0$, which gives
$$\int dt = \pm \int \frac{dr}{f(r)} + const.$$
So, the null coordinates are $t\pm\int \frac{dr}{f(r)}$. Is this method correct?
What about more complicated metrics that have off-diagonal terms like $f(r,t) dt dr$? What do we do in those cases?
 A: I don't know if it is the best way to look for such coordinate (probably not), but for general metric $$ds^2=g_{tt}dt^2+2g_{tr}dtdr+g_{rr}dr^2$$
on given timelike subspace you can look for coordinate transformations $t\rightarrow t'(t,r), r\rightarrow r'(t,r)$ that takes metric into the form 
$$ds^2=2g_{t'r'}dt'dr'+g_{r'r'}dr'^2.$$ The coordinate $t'$ will then be null, since $g(\partial_{t'}, \partial_{t'})$ is trivially zero.
Transforming the metric and collecting $dt'^2$ terms you get the condition:
$$
0=g_{tt}\left(\frac{\partial t}{\partial t'}\right)^2+2g_{tr}\frac{\partial t}{\partial t'}\frac{\partial r}{\partial t'}+g_{rr}\left(\frac{\partial r}{\partial t'}\right)^2
$$
Edit to answer the comment:
The solution of above differential condition is certainly not unique. At the subspace, there are two null directions at every point, but once you pick one direction at one point, you need to (in general) stick with it, because the coordinate vector field $\partial_{t'}$ needs to be smooth. So the direction at every point is (in general) unique, but not so the "length" of vectors (the actual length is zero, but one can still compare two vectors of same direction by the parameter $\alpha$ from equation $v_1=\alpha v_2$) of the coordinate vector field. Obviously any vector field of the form 
$$h(r',t')\partial_{t'}$$ 
is still null. For such vector field to lead to a coordinate, you need to be able to find vector field $\partial_{r''}=u_{t'}(r',t')\partial_{t'}+u_{r'}(r',t')\partial_{r'}$ with vanishing commutator:
$$0=[h(r',t')\partial_{t'},\partial_{r''}]=h \partial_{t'} u_{t'}+ h \partial_{t'} u_{r'}-u_{t'}\partial_{t'} h -u_{r'}\partial_{r'} h.$$
Such vector field you can find by requiring these two pairs to vanish independently:
$$0=h \partial_{t'} u_{t'} - u_{t'}\partial_{t'} h $$
$$0= h \partial_{t'} u_{r'}-u_{r'}\partial_{r'} h $$
These can be rewritten to the form :
$$ \partial_{t'} u_{t'} =H_{t'} u_{t'} $$
$$\partial_{t'} u_{r'}=H_{r'}u_{r'},$$
where $H_{i}=\partial_i h / h$ are known functions.These first order partial differential equations can be always solved, so indeed every vector field of the form $$\partial_{t''}=h(r',t')\partial_{t'}$$ leads to null coordinate also (assuming $h \neq 0$). Which one of all possible null coordinates you want is then up to you and the application.  
A: Null coordinates are, as the name imply, coordinates along which some null curves flow. A simple way to do this is to use the flow-box theorem.
Take two vector fields, $X$ and $Y$. $X$ and $Y$ are null vector fields if $g(X,X) = g(Y,Y) = 0$, and $X, Y \neq 0$. There's a variety of such vector fields you can take, roughly corresponding to the rotation along the null cone at each point. 
If, in addition, $X$ and $Y$ are nowhere proportional (ie $X \neq \alpha Y$), you can define locally some coordinates $(u,v)$ such that $X = \partial_u$, $Y = \partial_v$, by taking $u$ and $v$ to be parameters of the integral curves of $X$ and $Y$.
Example : for 
\begin{equation}
ds^2 = -f(r) dt^2 + \frac{1}{f(r)} dr^2 + r^2 d\Omega^2
\end{equation}
We are looking for two null vector fields, obeying the equation
\begin{equation}
-f(r) X_t^2 + \frac{1}{f(r)} X_r^2 + r^2(X_\theta^2 + \sin^2\theta X_\varphi^2) = 0
\end{equation}
There's a variety we can use, but the simplest case is $X_\theta = X_\varphi = 0$. So we are left with 
\begin{equation}
X_r^2 = f^2(r) X_t^2 
\end{equation}
Two simple choices are $X = (1, f(r), 0, 0)$ and $Y = (1, -f(r), 0, 0)$. It's not too hard to see that they are not proportional, and they are nowhere vanishing, as long as $f(r) \neq 0$.
What is the flow of those vector fields? Take a curve $x(\lambda)$, its flow is
\begin{equation}
\dot{x}(\lambda) = X(x(\lambda))
\end{equation}
For our vector fields, that would be
\begin{eqnarray}
\dot{x}_t(u) &=& 1\\
\dot{y}_t(v) &=& 1\\
\dot{x}_r(u) &=& f(x(u))\\
\dot{y}_r(v) &=& -f(y(v))\\
\end{eqnarray}
Pretty obviously, $x_t = u$, $y_t = v$ (I'm picking a zero constant of integration here), and 
\begin{eqnarray}
u &=& \int_0^r \frac{dy}{f(y)}\\
v &=& -\int_0^r \frac{dy}{f(y)}\\
\end{eqnarray}
Therefore, our new coordinates are defined by
\begin{eqnarray}
u &=& \frac{1}{2} (t + \int_0^r \frac{dy}{f(y)})\\
v &=& \frac{1}{2} (t -\int_0^r \frac{dy}{f(y)})\\
\end{eqnarray}
Probably a bit more effort than this is worth, but on the other hand, it will generalize locally to any metric.
