Why is gradient the direction of steepest ascent?

Each component of the gradient tells you how fast your function is changing with respect to the standard basis. It's not too far-fetched then to wonder, how fast the function might be changing with respect to some arbitrary direction? Letting $\vec v$ denote a unit vector, we can project along this direction in the natural way, namely via the dot product $\text{grad}( f(a))\cdot \vec v$. This is a fairly common definition of the directional derivative.

We can then ask in what direction is this quantity maximal? You'll recall that $$\text{grad}( f(a))\cdot \vec v = |\text{grad}( f(a))|| \vec v|\text{cos}(\theta)$$

Since $\vec v$ is unit, we have $|\text{grad}( f)|\text{cos}(\theta)$, which is maximal when $\cos(\theta)=1$, in particular when $\vec v$ points in the same direction as $\text{grad}(f(a))$.

Other answers are correct in using the directional derivative to show that the gradient is the direction of steepest ascent/descent. However, I think it is instructive to look at the definition of the directional derivative from first principles to understand why this is so (it is not arbitrarily defined to be the dot product of the gradient and the directional vector).

Let $f(\mathbf{x}):\mathbb{R}^n \rightarrow \mathbb{R}$. The partial derivatives of $f$ are the rates of change along the basis vectors of $\mathbf{x}$:

$\textrm{rate of change along }\mathbf{e}_i = \lim_{h\rightarrow 0} \frac{f(\mathbf{x} + h\mathbf{e}_i)- f(\mathbf{x})}{h} = \frac{\partial f}{\partial x_i}$

Each partial derivative is a scalar. It is simply a rate of change.

The gradient of $f$ is then defined as the vector:

$\nabla f = \sum_{i} \frac{\partial f}{\partial x_i} \mathbf{e}_i$

We can naturally extend the concept of the rate of change along a basis vector to a (unit) vector pointing in an arbitrary direction. Let $\mathbf{v}$ be such a vector, i.e., $\mathbf{v} = \sum_{i} \alpha_i \mathbf{e}_i$ where $\sum_{i} \alpha_i^2 = 1$. Then:

$\textrm{rate of change along }\mathbf{v} = \lim_{h\rightarrow 0} \frac{f(\mathbf{x} + h\mathbf{v}) - f(\mathbf{x})}{h}$

Again, this quantity is a scalar.

Now, it can be proven that if $f$ is differentiable at $\mathbf{x}$, the limit above evaluates to: $(\nabla f) \cdot \mathbf{v}$. This is a dot product of two vectors, which returns a scalar.

We know from linear algebra that the dot product is maximized when the two vectors point in the same direction. This means that the rate of change along an arbitrary vector $\mathbf{v}$ is maximized when $\mathbf{v}$ points in the same direction as the gradient. In other words, the gradient corresponds to the rate of steepest ascent/descent.

Consider a Taylor expansion of this function, $$f({\bf r}+{\bf\delta r})=f({\bf r})+(\nabla f)\cdot{\bf\delta r}+\ldots$$ The linear correction term $(\nabla f)\cdot{\bf\delta r}$ is maximized when ${\bf\delta r}$ is in the direction of $\nabla f$.