Can a gradient vector field with no equilibria point in every direction?
A quick solution for $n=2$, and an explanation of where it came from: $$\nabla (e^x \sin y) = (e^x \sin y, e^x \cos y).$$ The right hand side is never zero, but does assume every nonzero value in $\mathbb{R}^2$.
Motivation: If $f: \mathbb{C} \to \mathbb{C}$ is holomorphic, then $\nabla(\mathrm{Re}(f)) = (\mathrm{Re}(f'), -\mathrm{Im}(f'))$. I looked for an $f'$ (namely $e^z$) which takes $\mathbb{C}$ onto $\mathbb{C}_{\neq 0}$, and then integrated it to find $f$.
Inspired by this, let's try $$\nabla(e^{x_0} \cos(x_1^2+x_2^2+\cdots + x_n^2)) =$$ $$e^{x_0} (\cos(x_1^2+\cdots + x_n^2), - x_1 \sin(x_1^2+\cdots + x_n^2), \cdots, -x_n \sin(x_1^2+\cdots + x_n^2)).$$
First, we note that the gradient is not zero. The first coordinate only vanishes if $x_1^2+ \cdots + x_n^2$ is of the form $(2k+1) \pi$. But, in this case, at least one of $x_1$, $x_2$, ..., $x_n$ are nonzero; say $x_j$. Then $- x_j \sin(x_1^2+\cdots + x_n^2) = \pm x_j \neq 0$.
Now, let's take a nonzero vector $(v_0, \ldots, v_n)$. We need to go through cases.
If $v_0 = 0$, choose $(x_1, \ldots, x_n)$ proportional to $(v_1, \ldots, v_n)$ and such that $\sin(x_1^2+\cdots+x_n^2)$ has the right sign.
If $v_1 = \cdots = v_n = 0$ and $(-1)^k v_0>0$, take $x_1^2+\cdots+x_n^2 = k \pi$.
If neither of those cases holds, we'll take $(x_1, \ldots, x_n)$ of the form $(t v_1, \ldots, t v_n)$ for some $t$ to be determined soon. Set $s = x_1^2+ \cdots + x_n^2$. Then our vector points in direction $\pm (- \cot(t^2 s), v_1, v_2, \ldots, v_n)$. Since cotangent is surjective, we can choose $t$ such that $\cot (t^2 s) = v_0$, and we can get the sign right.
It is possible to have a vector field without equilibrium point and takes every possible directions in $\mathbb{R^n}$ for any $n \ge 2$.
The proof for the 2-dimension case is given below. From that, vector fields for higher dimension cases can be constructed.
Let $\rho : (-\epsilon,\infty) \to \mathbb{R}$ be any $C^\infty$ function such that
- $\rho(t)$ is even over $(-\epsilon,\epsilon)$ and $\rho(0) = \rho'(0) = 0$.
- $|r\rho'(r)| < 1$ for all $r \ge 0$.
- for some $r_c > 0$, $\rho(r_c) = 2\pi$.
- for some $r_m > 0$, $N \in \mathbb{Z}$, $\rho(r) = 2\pi N$ for all $r \ge r_m$.
For any point $(x,y) \in \mathbb{R}^2$, let $(r,\theta)$ be the corresponding
polar coordinates.
Let $\hat{e}_x, \hat{e}_y, \hat{e}_r, \hat{e}_\theta$ be unit
vectors along the $x$, $y$, $r$ and $\theta$ directions respectively.
Consider following function $$U(x,y) = r\cos(\theta - \rho(r)) = x\cos\rho(r) + y\sin\rho(r)$$ Its gradient is given by:
$$\begin{align} \vec{\nabla} U &= \hat{e}_x \left( \cos\rho + \frac{x\rho'}{r}( -x \sin\rho + y\cos\rho )\right) + \hat{e}_y \left( \sin\rho + \frac{y\rho'}{r}( -x \sin\rho + y\cos\rho )\right)\\ &= \hat{e}_r\left(\cos(\theta - \rho) + r\rho'\sin(\theta-\rho)\right) -\hat{e}_\theta\left(\sin(\theta-\rho)\right) \end{align} $$ From these two expressions, it is easy to deduce
- $\vec{\nabla} U$ is well defined at $(0,0)$ because $\rho'(0) = 0$.
- $\vec{\nabla} U(0,0) = \hat{e}_x$ because $\rho(0) = \rho'(0) = 0$.
- $\vec{\nabla} U \ne 0$ for $r > 0$ because $$\begin{align} |\vec{\nabla} U| &= |\hat{e}_r\left(\cos(\theta - \rho) + r\rho'\sin(\theta-\rho)\right) -\hat{e}_\theta\left(\sin(\theta-\rho)\right)|\\ &\ge |\hat{e}_r\left(\cos(\theta - \rho)\right) -\hat{e}_\theta\left(\sin(\theta-\rho)\right)| - |\hat{e}_e \left(r\rho'\sin(\theta-\rho)\right)|\\ &\ge 1 - |r\rho'| > 0 \end{align} $$
Combine these, we find $\nabla U \ne 0$ over $\mathbb{R}^2$. This allow us to define following unit vector field globally on $\mathbb{R}^2$.
$$\hat{u}(x,y) = \frac{\vec{\nabla}U(x,y)}{\left|\vec{\nabla}U(x,y)\right|}$$
Over the curve $\displaystyle\;\gamma : [0,r_c] \ni t \quad\mapsto\quad (t\cos\rho(t),t\sin\rho(t)) \in \mathbb{R}^2\;$, we have $$\hat{u}(t) = \hat{e}_r = \hat{e}_x \cos\rho(t) + \hat{e}_y\sin\rho(t)$$ As we move from $\gamma(0)$ to $\gamma(r_c)$ along $\gamma$, $\hat{u}(\gamma(t))$ will cover the whole $S^1$ at least once. The map $\hat{u}$ is surjective on $\gamma$ and hence on $B(0,r_c)$ and $\mathbb{R}^2$. This means $U$ is a counter-example of $V$ for the 2-dimenisonal case.
As an illustration, following is a picture of what $U(x,y)$ will typically look like under this construction.
$\hspace0.75in$
This particular $U(x,y)$ is computed using
$$\rho(t) = 2\pi + (f(t)-2\pi)g(t) \quad\text{ where }\quad \begin{align} f(t) &= \log\sqrt{(e^{4\pi}-1)t^2 + 1}\\ g(t) &= \begin{cases} 1, & t \le 1\\ 1 - e^{-1/(e^{t-1} -1)}, & t > 1\end{cases} \end{align} $$ and plotted over the region $|x|,|y| \le 2$. This $\rho(t)$ satisfies the first three requirement in the beginning of this answer with $r_c = 1$. For $r \gg r_c $, $\rho(t)$ tends to $2\pi$ exponentially as $t \to \infty$.
Back to the case of higher dimensions. Let's say we do have a $\rho$ that satisfies all four conditions. Consider following function:
$$V(x_1,x_2,\ldots,x_n) \stackrel{def}{=} U\left( x_1, \sqrt{x_2^2 + x_3^2 + \cdots + x_n^2} - (r_m+1)\right)$$
Since $B(0,r_m) \supset B(0,r_c)$, the map $\frac{\vec{\nabla}V}{|\vec{\nabla}V|}$ will be surjective over the hyper-torus $$ \bigg\{ (x_1,\ldots,x_n) : x_1^2 + \left(\sqrt{x_2^2 + x_3^2 + \cdots + x_n^2} - (r_m+1)\right)^2 \le r_m^2 \bigg\}$$ Since $U(x,y) = x$ whenever $x^2 + y^2 \ge r_m^2$, $V(x_1,\ldots,x_n)$ equals to $x_1$ over the complement of the hyper-torus. As that contains the $x_1$-axis, the square root in the definition of $V$ won't cause any problem. This makes $V$ smooth everywhere. Finally, it is easy to check $\nabla V$ never vanishes.
Update
About the question how this particular form of $\rho(t)$ is chosen.
First, $\rho(r)$ cannot change too rapidly or $\nabla U$ may become zero somewhere.
Second, we need $\rho(r)$ span across around a large enough interval to make $\hat{u}(r)$ surjective.
This suggest us to choose a $\rho(r)$ with $\rho'(r)$ as close to $\frac{1}{r}$ as possible. This means $$\rho(r) \sim \log r + constant.$$
In order for $U(x,y)$ to be smooth at the origin, I replace $\log r$ by $\log\sqrt{(e^{4\pi}-1)r^2 + 1}$. The constants are chosen so that when $r$ changes from $0$ to $1$, $\rho(r)$ pick up a change of $2\pi$ which we need.
Finally, we want $\rho(r)$ decay back to a constant multiple of $2\pi$ for large $r$ while keeping smooth all the way. We need this in order to extend the result to higher dimensions. This can be done using a smooth Partition of unity. The function $g(r)$ in the formula doesn't really have any specific meaning. It is simply something from trial and error which works.