how to integrate $\int\underbrace{x^{x^{\cdot^{\cdot^x}}}}_ndx$

If you replace $x$ by $x+1$ then you have $ \int (x+1)^{(x+1)^{...^{(x+1)}}} dx $ and the exponential-tower has an interesting power series, whose coefficients at the leading terms become constant in spite of increasing height $h$ (that $n$ in your formula) . Then you can integrate termwise to have a power series for the integral.
For instance, the tower of iteration height $h=6$ has the power series $$ 1 + x + x^2 + 3/2 \cdot x^3 + 7/3 \cdot x^4 + 4 \cdot x^5 + 283/40 \cdot x^6 + 4321/360 \cdot x^7 + O(x^8)$$ where the first terms up to $4 \cdot x^5$ stay constant for all higher iterations/exponential towers. (I do not know the range of convergence at the moment, maybe it is $ \small \eta-1 \approx 0.4446... $ where $ \small \eta = \exp(\exp(-1))$ because of the range of convergence for the exponential-tower of infinite height).

The termwise integration gives $$ x + 1/2 \cdot x^2 + 1/3 \cdot x^3 + 3/8 \cdot x^4 + 7/15 \cdot x^5 + 2/3 \cdot x^6 + 283/280 \cdot x^7 + O(x^8) $$ for the indefinite integral of the height $h=6$ exponential tower.

I get, using Pari/GP by the "explicite" integration (and substitution $x+1$ for $x$)

intnum(x=-1,0,(1+x)^(1+x)^(1+x)^(1+x)^(1+x)^(1+x))
%379 = 0.710658941398

which should be the correct value wrt to truncation to the shown digits. The use of the power series gives the value $ 0.710452400137$ which is inaccurate from the fourth digit (although I also applied Euler-summation for the diverging terms which have alternating signs), so the power series should be in principle usable also for higher exponential towers for small integration bounds and a more improved summation-procedure.

Unfortunately, there's a replacement $x \to (x+1)$ inside the integral and I do not know whether this is legitimate (I'm nearly illiterate with integration)

MathWorld has already done for this.

Please see the formula (10) in http://mathworld.wolfram.com/PowerTower.html.

I have found this solution: $$\mathrm{\int{x^{x^{x^…}}dx=\sum_{n=1}^\infty (-1)^nn^{n-2} Q\left(n,-ln(x)\right)+C}}$$

This works as a definite integral for $$\mathrm{e^{-\frac1e}\le x\le e^\frac 1e}$$

Here is my

work