The procedure of dimensional regularization for UV-divergent integrals is generally described as first evaluating the integral in dimensions low enough for it to converge, then "analytically continuing" this result in the number of dimensions $d$. I don't understand how this could possibly work conceptually, because a d-dimensional integral $I_d$ is only defined when $d$ is an integer greater than or equal to 1, so the domain of $I_d$ is discrete, and there's no way to analytically continue a function defined on a discrete set.
For example, in Srednicki's QFT book, the key equation from which all the dim reg results come is (pg. 101) "… the area $\Omega_d$ of the unit sphere in $d$ dimensions … is $\Omega_d = \frac{2 \pi ^{d/2}}{\Gamma \left( \frac{d}{2} \right) };$ (14.23)". (Note: see edit below.) But this is highly misleading at best. The area of the unit sphere in $d$ dimensions is $\frac{2 \pi^{d/2}}{\left( \frac{d}{2} – 1 \right) !}$ if $d$ is even and $\geq 2$, it is $\frac{2^d \pi^\frac{d-1}{2} \left( \frac{d-1}{2} \right)! }{(d-1)!}$ if $d$ is odd and $\geq 1$, and it is nothing at all if $d$ is not a positive integer. These formulas agree with Srednicki's when $d$ is a positive integer, but they avoid giving the misleading impression that there is a natural value to assign to $\Omega_d$ when it isn't.
Beyond purely mathematical objections, there's a practical ambiguity in this framework – how do you interpolate the factorial function to the complex plane? Srednicki chooses to do so via the Euler gamma function without any explanation. But there are other possible interpolations which seem equally natural – for example, the Hadamard gamma function or Luschny's factorial function. (See http://www.luschny.de/math/factorial/hadamard/HadamardsGammaFunction.html for more examples.) Why not use those?
In fact, these two alternative functions are both analytic everywhere, so you can't use them to extract the integral's pole structure, which you need in order to cancel the UV infinities. To me, this suggests that the final results of dim reg might be highly dependent on your choice of interpolation scheme, therefore requiring a justification for using the Euler gamma function. Could we prove to a dim reg skeptic that all results for physical observables are independent of the interpolation scheme? (Note that this is a stronger requirement than showing they are independent of the fictitious mass parameter $\tilde{\mu}$.)
(I know that the Bohr-Mollerup theorem shows that the Euler gamma function uniquely has certain "nice" properties, but I don't see why those properties are helpful for doing dim reg.)
I'm not looking for a hyper-technical treatment of dim reg, just a conceptual picture of what it even means to analytically continue a function from the discrete set of positive integers.
Edit: It appears that the details of exactly which field-theory results do and do not depend on the choice of regularization scheme are not well-understood; see this paper for one discussion.
Edit: kaylimekay correctly points out below that the relevant $\Gamma(d)$ term is actually the one that comes from the radial integral, not the one that comes from the angular integral. But I don't think that this really helps solve the problem at all.
The issue basically boils down to ambiguous notation; two qualitatively different "exponentiation" functions use the same notation. To be explicit, I'll give them two different names.
The first is the function $\mathrm{expNat}: \mathbb{R} \times \mathbb{N} \to \mathbb{R}$ defined by the usual "repeated multiplication" that we learn in Algebra 1 (or whenever):
$$\mathrm{expNat}(x, n) := \underbrace{x \times x \times \dots \times x}_{n \text{ times}}.$$
The second is the more complicated function $\mathrm{expReal}: \mathbb{R^+} \times \mathbb{R} \to \mathbb{R^+}$ (where the exponent can be an arbitrary real number) that we learn in Calculus 1 (or whenever). There are several equivalent ways to define this function, but for concreteness we'll take
$$\mathrm{expReal}(x, y) := \exp(y \ln(x)),$$
with $\exp(z)$ defined as the unique solution to the initial value problem $\exp'(z) \equiv \exp(z),\ \exp(0) = 1$, and $\ln(x)$ defined as its inverse.
The parameter that kaylimekay calls $\alpha$ comes from using Feynman's formula to rewrite the product of $\alpha$ different terms, and it is always a natural number. I'll rename it $n$ to make this clear. Rewritten in more explicit notation, the relevant radial integral is
$$
\int_0^\infty dk\, \frac{\mathrm{expNat}(k, d-1)}{\mathrm{expNat}(k^2 + \Delta^2, n)}.
$$
It's not clear to me why we can extend $\mathrm{expNat}$ to $\mathrm{expReal}$; as explained above, this isn't a legitimate analytic continuation because it starts from a discrete set. And there are other possible extensions as well; for example, it would be equally correct to write kaylimekay's expression as
$$
\int_0^\infty dk\, \frac{k^{(d-1) \cos(2 \pi d)}}{(k^2 + \Delta^2)^n},
$$
but writing it this way would suggest a different extension to the real domain for $d$, which (as far as I can tell) could give you different final answers after analytic continuation.
Best Answer
Every regularization scheme is somewhat arbitrary. There are three popular regularization schemes when it comes to path integrals and their associated perturbative divergent integrals: time slicing, mode regularization, and dimensional regularization.
Time slicing is the usual procedure used to derive the path integral, and it is the discretization of time into finite time intervals.
Mode regularization is essentially an UV cut off, i.e. the truncation of the high-energy modes in the Fourier expansion of the path.
Dimensional regularization is performed as you described exploiting (one) generalization of the factorial to complex numbers.
In any case, the regularization is a limiting procedure, a finite (or different from zero) parameter is introduced such that all the integrals become finite, then they're manipulated in a way that the result in the limit when the parameter goes to infinity (zero) remains finite. In principle, and in practice, the final result may be dependent on the scheme chosen. Therefore it is necessary to introduce counterterms such that all the results agree with each other. This is done somewhat ad hoc, but luckily the counterterms are fixed at a low (second) order in the perturbation expansion in many situations.
The procedures chosen are all in some sense arbitrary, for (at least for the moment) there is not a satisfactory and unambiguous mathematical definition of the involved path integrals/QFT perturbative expansions. The dimensional regularization is often preferred for essentially one reason (as far as I know), it is the easiest to deal with: the resulting counterterm in fact is relativistically covariant (and that is important in relativistic theories/in the presence of a curved background) and the additional vertices coming form the counterterm at higher loops are easy to compute.
Now my guess is that it could be possible to regularize also using one of the other "complex extensions" of the factorial you mentioned, but in all likelihood the resulting counterterms would be different and maybe not covariant.
For a more detailed discussion on regularization schemes I suggest to read this book of Bastianelli and van Nieuwnehuizen.