Distribution of log-likelihood gradient My question is simple: Is there any results regarding the distribution of log-likelihood function gradient? 
It may be asymptotic results as well.
 A: We have to view the log-likelihood not as a function of the parameters given the sample, but as a function of the random variables whose realization create the sample. In such a case, and to examine a concrete example, we look at the normal case.  
The log-likelihood of the sample of $n$ i.i.d variables each following a normal $N(\mu, \sigma^2)$ distribution, is
$$\ln L(\mu,\sigma^2,\mathbf X) = C -\frac n2\ln\sigma^2 -\frac 1{2\sigma^2}\sum_{i=1}^n(X_i-\mu)^2$$
The derivatives that form the gradient are
$$\frac{\partial \ln L}{\partial \mu} = \frac 1{\sigma^2}\sum_{i=1}^n(X_i-\mu) = \frac 1{\sigma}\sum_{i=1}^n\left(\frac{X_i-\mu}{\sigma}\right)  \sim N\left(0,\frac n{\sigma^2}\right) $$
We note that $$\frac {1}{\sqrt n}\cdot\frac{\partial \ln L}{\partial \mu} \sim N\left(0,\frac 1{\sigma^2}\right) \tag{1}$$
Also (we take the derivative with respect to $\sigma^2$ and not $\sigma$)
$$\frac{\partial \ln L}{\partial \sigma^2} = -\frac n  {2\sigma^2} + \frac 1{2\sigma^4}\sum_{i=1}^n(X_i-\mu)^2 = -\frac n  {2\sigma^2}  +\frac 1{2\sigma^2} \sum_{i=1}^n\left(\frac{X_i-\mu}{\sigma}\right)^2 $$
The sum follows a chi-square with $n$ degrees of freedom, so scaled by $1/2\sigma^2$ we get
$$\frac 1{2\sigma^2} \sum_{i=1}^n\left(\frac{X_i-\mu}{\sigma}\right)^2 = Q_n \sim \Gamma_d\left(\frac n2, \frac 1{\sigma^2}\right)$$
where $\Gamma_d()$ is a Gamma distribution with the shape, scale parametrization. 
We have $E(Q_n) = \frac n  {2\sigma^2}$ so we get
$$\frac{\partial \ln L}{\partial \sigma^2} = Q_n-E(Q_n)$$
Moreover $\text{Var}(Q_n) = \frac n{2\sigma^4}$. The Central Limit Theorem applies and so
$$\frac {\sigma^2\sqrt2}{\sqrt n}\cdot \frac{\partial \ln L}{\partial \sigma^2} \xrightarrow{d}N(0,1) \Rightarrow \frac {1}{\sqrt n}\cdot \frac{\partial \ln L}{\partial \sigma^2} \xrightarrow{d}N\left(0,\frac 1{2\sigma^4}\right) \tag{2}$$
Equations $(1)$ and $(2)$ give us the asymptotic marginal distributions of the elements of the gradient of an i.i.d. normal log-likelihood (the first one is also exact as a finite-sample distribution). But for the joint distribution, we need to consider the possible covariance of the two elements of the gradient. Since both have zero expected value, we have that
$$\text{Cov}\left( \frac {1}{\sqrt n}\frac{\partial \ln L}{\partial \mu}, \frac {1}{\sqrt n}\frac{\partial \ln L}{\partial \sigma^2}\right) =  \frac {1}{ n}E\left[\left(\frac 1{\sigma}\sum_{i=1}^n\left(\frac{X_i-\mu}{\sigma}\right)\right)\cdot\left(\frac 1{2\sigma^2} \sum_{i=1}^n\left(\frac{X_i-\mu}{\sigma}\right)^2-\frac n  {2\sigma^2}\right)\right]$$
All cross products, and all products with the constant term will have expected value zero, since the variables are independent, So we are left with
$$\text{Cov}\left( \frac {1}{\sqrt n}\frac{\partial \ln L}{\partial \mu}, \frac {1}{\sqrt n}\frac{\partial \ln L}{\partial \sigma^2}\right) = \frac 1nE\left[\frac 1{2\sigma^3}\sum_{i=1}^n\left(\frac{X_i-\mu}{\sigma}\right)^3\right]$$
But the standardized variables are standard normals raised to the 3d power. Then, we are looking at the 3d raw moment of the standard normal, which is zero. So we conclude that
$$\text{Cov}\left( \frac {1}{\sqrt n}\frac{\partial \ln L}{\partial \mu}, \frac {1}{\sqrt n}\frac{\partial \ln L}{\partial \sigma^2}\right) = 0$$
...which implies also that the asymptotic joint distribution will be normal. Therefore we obtain
$$\frac 1{\sqrt n}\left[\begin{matrix}
\frac{\partial \ln L}{\partial \mu}\\
\frac{\partial \ln L}{\partial \sigma^2}\\
\end{matrix}\right] \xrightarrow{d} N\left(\left[\begin{matrix}
0\\
0\\
\end{matrix}\right] ,\left[\begin{matrix}
\frac 1{\sigma^2} & 0\\
0 & \frac 1{2\sigma^4}\\\end{matrix}\right] \right)$$
