First: $\tau$ is not a state variable but a forcing term (torque in the equation), so the perturbation $\tau = \tau_0 + \delta\tau$ is unnecessary.
Second: In principle, you need not necessarily choose a fixed point equilibrium to linearize about (though it certainly helps for any practical implementation of a controller).
First let's derive the general solution, then we'll see how it works for the specific problem.
Let $x$ be an arbitrary state vector in the Hilbert space $X$ and suppose the time evolution of $x$ is given via $\dot{x} = f(x,t)$ for a smooth mapping $f:X\rightarrow X$. As obareey states, $f$ enjoys a Taylor expansion about the function $y$ given by
$$
f(x) = f(y) + \partial_xf(y)(x-y) + \frac{1}{2}\partial^2_xf(y)(x-y)(x-y)^T+\cdots,
$$
where
$$
\partial_xf(y) = \frac{\partial f}{\partial x}|_{x = y},\ \ \partial^2_xf(y) = \frac{\partial^2 f}{\partial x\partial x^T}|_{x = y},
$$
etc. The linearization is of $f$ is defined to be the truncation of this series to first (that is, linear) order. This produces the linearized equation
$$
\dot{x} \approx f(y) + \partial_xf(y)(x-y).
$$
Since this approximation is based on the assumption $x = y + \delta$ where $|\delta|/|y| \ll 1$, it follows that $\dot{x} = \dot{y} + \dot{\delta} = f(y) +\dot{\delta}$, and so we have
$$
\dot{\delta} = \partial_xf(y)\delta.
$$
Since $y$ is presumed known then $x = y + \delta$ is solved when $\delta$ is known and $\partial_xf(y)$ is just a (possibly time-dependent) matrix $A(t)$. By this method we have reduced the nonlinear equation to the equation
$$
\dot{\delta} = A(t)\delta.
$$
Now on to your problem:
The state vector is $x = [\theta, \dot{\theta}]^T$ and the nonlinear model is
$$
\frac{d\theta}{dt} = f_\theta = \dot{\theta},\ \ \frac{d\dot{\theta}}{dt} = f_{\dot{\theta}} = \frac{\tau}{m\ell^2}- \frac{g}{\ell}\cos\theta.
$$
Choosing the constant reference $y = [\theta_0,\dot{\theta}_0]^T$, the linearization is equivalent to finding the partials
$$\begin{align}
\frac{\partial f_\theta}{\partial \theta} = 0, &\ \frac{\partial f_\theta}{\partial \dot{\theta}} = 1\\
\frac{\partial f_\dot{\theta}}{\partial \theta} = \frac{g}{\ell}\sin\theta_0, &\ \frac{\partial f_\dot{\theta}}{\partial \dot{\theta}}=0\\
\end{align}$$
and therefore the linearized system is of the form
$$
\dot{\delta} = \begin{bmatrix}
0 & 1\\
\frac{g}{\ell}\sin\theta_0 & 0\\
\end{bmatrix}\delta
$$
where $\delta = [\theta-\theta_0,\dot{\theta}-\dot{\theta}_0]^T$. Choosing the fixed point $\theta_0 = \cos^{-1}(\tau/mg\ell)$ and $\dot{\theta}_0 = 0$ we have
$$
\begin{align}
\frac{d\theta}{dt} &= \dot{\theta}\\
\frac{d\dot{\theta}}{dt} &= \frac{g}{\ell}\sin(\cos^{-1}(\tau/mg\ell))\theta - \frac{g}{\ell}\sin(\cos^{-1}(\tau/mg\ell))\cos^{-1}(\tau/mg\ell)
\end{align}
$$
