Mathematica does all its computations in an orthonormal basis. You simply need to specify what coordinate system you're working in. So for you example, you just multiply by {0,0,1}

 e[r_, \[Theta]_, \[Phi]_, t_] := (Sin[\[Theta]]/r)[Cos[r - t] - Sin[r - t]/r] {0, 0, 1}

Apprently this is a pure wave in vacuum, as the divergence is zero:

 Div[e[r, \[Theta], \[CurlyPhi], t], {r, \[Theta], \[CurlyPhi]}, "Spherical"]
 (* 0 *)

Similarly, a pure Coloumb electric field would be 

 col[r_,t_,p_] := {1/r^2,0,0}
 Div[col[r,t,p], {r, t, p}, "Spherical"]
 (* 0 *)

I suggest you look at the the tutorials `tutorial/VectorAnalysis` and 
`tutorial/ChangingCoordinateSystems` and the functions linked therefrom for me.