I'm not actually sure whether "gamma regression" is officially defined (it doesn't appear to have a wikipedia page for example), but if I were to define it (and some googling around suggests I'm not alone here), I would define it as setting up my regression problem so that for a given input vector $\underline{x}$, I predict a value y. When I do this, I meant that I think that the true value of the target will be gamma distributed with mean y.
How does this differ from least squares? One setup which leads to OLS is that you say that for given $\underline{x}$, the target variable will be normally distributed and your prediction $y(\underline{x})$ is the mean of that distribution. Of course a normal distribution is parametrised through its mean and variance, but it turns out that you don't need to know the variance in order to calculate the cost function you need to optimise, and thus this parameter doesn't need to be passed to xgboost.
For the gamma distribution however, this is different. Let's go through the maths. For the gamma distribution parametrised as $\frac{1}{\Gamma (k)\theta ^{k}}x^{k-1}e^{-\frac{x}{\theta}}$, the mean is given by $k\theta$ and the variance by $k \theta ^{2}$
Thus let's reparametrise in terms of $\mu$ and $\theta$ so that the distribution is given by $\frac{1}{\Gamma (\frac{\mu}{\theta})\theta ^{\frac{\mu}{\theta}}}x^{\frac{\mu}{\theta}-1}e^{-\frac{x}{\theta}}$
So for a given dataset, if you predict a bunch of $\hat{y}_{i}$ for target values $y_{i}$, the likelihood is given by
$\prod _{i=1}^{N} \frac{1}{\Gamma (\frac{\hat{y}_{i}}{\theta})\theta ^{\frac{\hat{y}_{i}}{\theta}}}y_{i}^{\frac{\hat{y}_{i}}{\theta} -1}e^{-\frac{y_{i}}{\theta}}$
and thus the negative (xgboost assumes a cost function you're trying to minimise) log-likelihood is
$\sum _{i=1}^{N} \ln \Gamma (\frac{\hat{y}_{i}}{\theta}) + \frac{\hat{y}_{i}}{\theta}\ln \theta - (\frac{\hat{y}_{i}}{\theta}-1) \ln y_{i} + \frac{y_{i}}{\theta}$
Compare this to the Gaussian Regression case where the negative log likelihood is given by
$\frac{1}{\sigma ^{2}}\sum _{i=1}^{N} \left(y_{i} - \hat{y}_{i}\right)^{2}$
The the latter case, the $\frac{1}{\sigma ^{2}}$ is a constant term out the front. If you were doing linear regression or even xgboost without regularisation, this would mean that no matter what value you changed $\sigma$ to, the linear regressor/xgboost you trained would turn out to be exactly the same, so "Gaussian regression with $\sigma = 10$ and Gaussian regression with $\sigma = 1$ lead to the same predictions". This is no longer true when you have a regulariser, but you can always suck the value of $\sigma$ into the definition of the regulariser to get around this, and this is why the OLS formula never includes a $\sigma$ in it.
In the gamma case however, because of the $\theta$ factor contained in the $\Gamma$ function and the $\ln \theta$, you can't just pull the factor of $\theta$ outside of the summation.
For xgboost, you now need to pass it the elementwise first and second derivatives of the cost function wrt $\hat{y}_{i}$. This is where basic calculus doesn't get you all the way, you'll likely need to look up that the derivative of the logarithm of the gamma function is given by the digramma function $\psi (z)$.
The (elementwise) first derivative of the loss will be given by (by the xgboost definition which is $G_{i}=\frac{\partial L}{\partial \hat{y}_{i}}$):
$G_{i} = \frac{1}{\theta}\psi (\frac{\hat{y}_{i}}{\theta}) + \frac{1}{\theta}\ln \theta - \frac{1}{\theta}\ln y_{i} $
The second derivative will require derivatives of the digamma function, I don't know much about this but some googling tells me you need the trigamma function $\psi _{1}(z)$ which is the derivative of the digamma function, thus
$H_{i}=\frac{1}{\theta ^{2}}\psi_{1}(\frac{\hat{y}_{i}}{\theta})$
Again, note that you still have to supply $\theta$ up front as a hyperparameter, pass this to xgboost and then train a new xgboost model every time you wish to investigate another $\theta$
Finally, it's worth noting that I just did this derivation myself today, I haven't lifted anything other than the definition of the gamma distribution from elsewhere, so there could easily be minor algebra error, I'd feel more comfortable if somebody else independently verified my workings.
Edit: Alternately, you could parametrise the other way around:
You could use k as your free parameter and $\theta = \frac{\mu}{k}$, thus your gamma distribution is $\frac{1}{\Gamma(k)(\frac{\mu}{k})^{k}}x^{k-1}e^{-\frac{xk}{\mu}}$
and thus your negative log-likelihood is given by
$L = \sum_{i=1}^{N}\left[\ln \Gamma (k) + k\ln \hat{y}_{i} - k\ln k - (k-1)\ln y_{i} + \frac{y_{i}k}{\hat{y}_{i}}\right]$
this parametrisation is easier to differentiate, you get
$G_{i}=\frac{k}{\hat{y}_{i}} - \frac{y_{i}k}{\hat{y}_{i}^{2}}$
and
$H_{i}=-\frac{k}{\hat{y}_{i}^{2}}+ 2\frac{y_{i}k}{\hat{y}_{i}^{3}}$
but you still need to pass k as a hyperparameter
