The problem
I'd like to (numerically) compute nested integrals of the type $$I(a;k) \equiv \left( \prod_{i=1}^k \int_0^\infty \mathrm{d}x_i \right) \exp \left( - \frac{1}{4a} \sum_{i=1}^k \left( x_i - x_{i+1} \right)^2 - \frac{a}{2} \sum_{i=1}^k (x_i + x_{i+1}) \right)$$ where $a > 0$ is a parameter and $x_{k+1} \equiv x_1$.
Analytic result for $k=1,2$
For $k = 1$ and $k = 2$ the integrals can be done analytically with results $$I(a;1) = 1/a, ~~~ I(a;2) = \sqrt{\frac{\pi}{2a}} ~ e^{a^3/2} ~ \text{Erfc}\left( \frac{a^{3/2}}{\sqrt{2}} \right),$$ but as far as I'm aware no closed-form expression is known for $k > 2$. One suggestion discussed in the comments was to write the argument of the exponential as a quadratic form $\boldsymbol{x}^\top A \boldsymbol{x} + \boldsymbol{J}^\top \boldsymbol{x}$ with $A/\boldsymbol{J}$ some $a$-dependent symmetric matrix/vector and attempt to use the same logic as in the derivation of the standard multiple-Gaussian integral (that is, orthogonally diagonalizing $A = O^\top D O$ and changing variables $\boldsymbol{y} = O \boldsymbol{x}$). This is complicated however due to the range of integration being $[0,\infty)^k$ and not $(-\infty,\infty)^k$ -- these are not Gaussian integrals.
My code so far and why it's wrong
So far I have written the following NIntegrate command:
int[a_, k_] :=
NIntegrate[
Exp[-1/(4 a) (Sum[(x[i] - x[i + 1])^2, {i, 1, k - 1}] + (x[k] -
x[1])^2) -
a/2 (Sum[(x[i] + x[i + 1]), {i, 1, k - 1}] + (x[k] + x[1]))],
Sequence @@ Table[{x[i], 0, \[Infinity]}, {i, k}] // Evaluate,
AccuracyGoal -> 3];
For values $a$ of order unity the integral should be of order unity, hence the unambitious AccuracyGoal. An accuracy of $10^{-3}$ for $a = O(1)$ would suffice for my purposes though it would be nice if it could be systematically improved. It seems like the integrals become more difficult to do when $a$ is small -- let's say we'd like $10^{-1} < a < 1$ but again ways to systematically enlarge the range would be nice.
To see that the above code gives the wrong answer we can compare this numerics with the exact result for the $k=2$ case:
$k=2$." />
For $a \lesssim 0.3$ the numerics clearly start failing beyond the error bar allowed by the AccuracyGoal setting, though Mathematica produces no error message. For $k=3,4$ an error message (a NIntegrate::slwcon) is produced when $a$ drops below some threshold value, but because of what we saw above it is doubtful that the code computes the correct value for the integral within the error bar even when $a$ is slightly larger and no error message is produced. I've tried increasing the WorkingPrecision but that gives an NIntegrate::precw error and doesn't solve the $k=2$ discrepancy.
Check on the correct solution
The numerical solution for $k > 2$ should be checked by computing also $$\mathrm{d} I(a;k)/\mathrm{d}a = \left( \prod_{i=1}^k \int_0^\infty \mathrm{d}x_i \right) \left[ \frac{1}{4a^2} \sum_{i=1}^k \left( x_i - x_{i+1} \right)^2 - \frac{1}{2} \sum_{i=1}^k (x_i + x_{i+1}) \right] \exp \left( - \frac{1}{4a} \sum_{i=1}^k \left( x_i - x_{i+1} \right)^2 - \frac{a}{2} \sum_{i=1}^k (x_i + x_{i+1}) \right),$$ computing its (numerical) integral from some reference value $a_0$ to a variable upper limit $b$ and showing that it agrees with $I(b;k)-I(a_0;k)$ (for some range of $b$ and some moderate values of $k$).
