It is relatively simple to write the log-density of the negative binomial distribution in terms of its mean, and then use this to get a log-likelihood expression for the negative binomial GLM. For all values $y=0,1,2,...$ the log-density of the negative binomial distribution is:
$$\log f(y|r,\theta) = \log {y+r-1 \choose y} + r \log (1-\theta) + y \log (\theta).$$
If we let $\mu \equiv \mathbb{E}(Y)$ denote the mean of the distribution then we have the relationships:
$$\mu = r \cdot \frac{\theta}{1-\theta}
\quad \quad \quad \theta = \frac{\mu}{r+\mu}.$$
Thus, we can re-parameterise the log-density in terms of its mean parameter to obtain the log-likelihood for a single data point:
$$\ell_y(r,\mu) \equiv \log f(y|r,\mu) = \log {y+r-1 \choose y} + r \log (r) + y \log (\mu) - (r+y) \log (r+\mu).$$
The derivatives with respect to $\mu$ are:
$$\begin{aligned}
\frac{d\ell_y}{d \mu} (r,\mu)
&= \frac{y}{\mu} - \frac{r+y}{r+\mu}, \\[6pt]
\frac{d^2 \ell_y}{d \mu^2} (r,\mu)
&= -\frac{y}{\mu^2} + \frac{r+y}{(r+\mu)^2}. \\[6pt]
\end{aligned}$$
Letting $\hat{\mu}$ be a critical point of the log-density, we obtain:
$$\frac{d^2 \ell_y}{d \mu^2} (r,\hat{\mu})
= -\frac{y}{\hat{\mu}} \Bigg[ \frac{1}{\hat{\mu}} + \frac{1}{r+\hat{\mu}} \Bigg] <0.$$
This shows that all the log-likelihood for a single data point has a single maximising value for its mean parameter, and so it is quasi-concave with respect to the mean parameter. Unfortunately, this simple result gets more muddied when we transition to the log-likelihood for multiple data points in the negative binomial GLM.
Negative binomial GLM: In the negative binomial GLM we have data $Y_i \sim \text{NegBin}(r, \mu_i)$ where the value $\mu_i$ is determined from a corresponding explanatory variable $\mathbf{x}_i$ using a link function. Written in general terms, using the mean vector $\boldsymbol{\mu}$ the log-likelihood of the data $\mathbf{y}=(y_1,...,y_n)$ is:
$$\ell_\mathbf{y} (r, \boldsymbol{\mu}) = \sum_i \log {y_i+r-1 \choose y_i} + n r \log (r) + \sum_i y_i \log (\mu_i) - \sum_i (r+y_i) \log (r+\mu_i).$$
In the standard negative binomial GLM we use the link function $\log \mu_i(\boldsymbol{\beta}) = \sum_k \beta_k \cdot x_{i,k}$, which gives the likelihood function:
$$\begin{aligned}
\ell_\mathbf{y} (r, \boldsymbol{\beta})
&= \sum_i \log {y_i+r-1 \choose y_i} + n r \log (r) \\
&\quad + \sum_i y_i \sum_k \beta_k \cdot x_{i,k} - \sum_i (r+y_i) \log (r+\mu_i(\boldsymbol{\beta})).
\end{aligned}$$
The derivatives of this function with respect to the model coefficients are:
$$\begin{aligned}
\frac{\partial \ell_\mathbf{y}}{\partial \beta_k} (r, \boldsymbol{\beta})
&= \sum_i y_i \cdot x_{i,k} - \sum_i (r+y_i) \cdot x_{i,k} \cdot \frac{\mu_i(\boldsymbol{\beta})}{r+\mu_i(\boldsymbol{\beta})}, \\[6pt]
\frac{\partial^2 \ell_\mathbf{y}}{\partial \beta_k \partial \beta_l} (r, \boldsymbol{\beta})
&= - \sum_i (r+y_i) \cdot \frac{\mu_i(\boldsymbol{\beta})}{r+\mu_i(\boldsymbol{\beta})} \Bigg( 1- \frac{1}{r+\mu_i(\boldsymbol{\beta})} \Bigg) \cdot x_{i,k} \cdot x_{i,l}. \\[6pt]
\end{aligned}$$
As you can see, these derivatives are more complicated than in the case of a single data point. This is owing to the fact that we now have $n$ data points, and we are looking at derivatives with respect to the model coefficients, rather than the mean value.
