Closed form of ∫₀^∞ x^(2n) arctan(x)/(e^(aπx)+1) dx via negative-order polygamma

The integral in question is Ω_n(a) = ∫_0^∞ x^(2n) arctan(x) / (e^(aπx) + 1) dx, for a natural number n and a real a > 0. It looks like a standard Fermi–Dirac-weighted moment, except the arctangent spoils every textbook substitution. No residue argument gives a clean answer, and a computer algebra system returns closed forms for small n without exposing the pattern behind them.

The problem was posed on Mathematics Stack Exchange. The asker had reduced it to a single integral of a digamma function minus its own asymptotic expansion, and then stalled. The interesting part is how that last step was finished. The trick is an explicit antiderivative built from negative-order polygamma functions, and it turns a tail integral into a limit that Stirling's series can evaluate.

Turning the arctangent into a parameter integral

The first move is to write arctan(x) = ∫_0^1 x / (x²y² + 1) dy and swap the order of integration. Rescaling x → x/y in the inner integral gives Ω_n = ∫_0^1 y^(-2n-2) ∫_0^∞ x^(2n+1) / ((x² + 1)(e^(aπx/y) + 1)) dx dy. The inner integrand is now a rational function times a Fermi factor, which is much better behaved.

Polynomial division handles the rational part. For natural n, t^n/(t+1) = (-1)^n/(t+1) + Σ_{k=0}^{n-1} (-1)^(n-k-1) t^k. With t = x², the inner integral splits into one term containing x/(x²+1) and a finite sum of pure power moments ∫ x^(2k-1)/(e^(bx)+1) dx.

Both pieces are known. The power moments give (2k-1)! (1 - 2^(1-2k)) ζ(2k) / b^(2k) with b = aπ/y. The x/(x²+1) term evaluates to ½(ψ(a/(2y) + ½) - ln(a/(2y))), where ψ is the digamma function.

The reduced integral and why it resists

After substituting y → 1/y, rewriting ζ(2k) through Bernoulli numbers B_(2k), and setting x = ay/2, everything collapses to one line:

Ω_n = (-1)^n · 2^(2n) / a^(2n+1) · ∫_{a/2}^∞ x^(2n) [ ψ(x + ½) - ln x - ½ Σ_{k=1}^{n} (1 - 2^(1-2k))/k · B_(2k) x^(-2k) ] dx

The bracket is ψ(x + ½) - ln x with exactly the first n terms of its asymptotic expansion subtracted. That subtraction is what makes the integral converge against the growing x^(2n). It also means a naive term-by-term approach cancels to nothing useful.

The asker also tried the repeated integration-by-parts identity ∫ x^n f(x) dx = Σ_k (-1)^k k! C(n,k) x^(n-k) f^(-k-1)(x). The obstacle is the limit at infinity: each iterated antiderivative of ψ grows, and the divergent pieces have to cancel exactly against the subtracted Bernoulli terms.

The antiderivative: P_m and the duplication formula

The accepted answer packages that identity into a named function. For m ≥ 0, define P_m(z) = Σ_{j=0}^{m} (-1)^j j! C(m,j) z^(m-j) ψ^(-j-1)(z), with ψ^(-1)(z) = log Γ(z) and each ψ^(-j-2) an antiderivative of ψ^(-j-1). Differentiating term by term telescopes to P_m'(z) = z^m ψ(z).

The shifted argument x + ½ is removed with the duplication formula ψ(x + ½) = 2ψ(2x) - ψ(x) - 2 log 2. That lets ψ(x+½) be written using only P_(2n) at x and 2x. The answer then defines

A_n(x) = 2^(-2n) P_(2n)(2x) - P_(2n)(x) - x^(2n+1)(log x + 2 log 2)/(2n+1) + x^(2n+1)/(2n+1)² - Σ_{k=1}^{n} (1 - 2^(1-2k)) B_(2k) / (2k(2n-2k+1)) · x^(2n-2k+1)

and checks that A_n'(x) is exactly the integrand above. The integral is therefore lim A_n - A_n(a/2). Stirling's expansion shows A_n(x) → C_n, with C_0 = -log 2 / 2 and, for n ≥ 1, C_n = (-1)^n (1 - 2^(-2n)) (2n)! ζ(2n+1) / (2^(2n+1) π^(2n)).

The closed form

Putting it together gives Ω_n(a) = (-1)^n 2^(2n) / a^(2n+1) · (C_n - A_n(a/2)). The asker expanded this into an explicit sum:

Ω_n(a) = Σ_{j=0}^{2n} (-1)^(n+j+1) j! C(2n,j) a^(-j-1) [ ψ^(-j-1)(a) - 2^j ψ^(-j-1)(a/2) ]
       + (1 - 2^(-2n)) (2n)! ζ(2n+1) / (2 a^(2n+1) π^(2n))
       + (-1)^n ln(2a) / (2(2n+1))
       + (-1)^n Σ_{k=1}^{n} (2^(2k-1) - 1) B_(2k) / (2k (2n-2k+1) a^(2k))
       - (-1)^n / (2(2n+1)²)

This matches the structure of the Mathematica results for n = 1, 2, 3. Each has paired ψ^(-j-1)(a) and ψ^(-j-1)(a/2) terms, one odd zeta value ζ(2n+1), a ln(2a) term with coefficient ±1/(2(2n+1)), and a constant ∓1/(2(2n+1)²). Examples are 1/18 for n = 1 and -1/50 for n = 2. The n = 0 case is simpler: Ω_0 = -ln Γ((a+1)/2)/a + ln 2/(2a) + ln π/(2a) + ½ ln(a/2) - ½.

The series route and why it stalls

A second answer tried the obvious alternative. It expands 1/(e^(aπx)+1) = Σ_k (-1)^k e^(-aπ(k+1)x) and integrates term by term. After one integration by parts, each term becomes an integral of x^(2n+1)/(x²+1) against a generalized exponential integral E_(-2n).

Mathematica evaluates those terms as Meijer G-functions G^{4,1}_{2,4} with argument π²a²(k+1)²/4 and a clear parameter pattern in n. The author called it 'probably not very useful', and that is fair. You end up with an alternating infinite series of Meijer G values instead of a finite expression.

Where this approach does not fit, and what to watch

• Integer n only. The polynomial division step needs n ∈ ℕ. Non-integer powers of x break the finite sum of Bernoulli terms.
• a > 0 is essential. Both the Fermi factor expansion and the convergence at infinity depend on it.
• Normalization of ψ^(-k). The relation dψ^(-j-2)/dz = ψ^(-j-1) fixes each negative-order polygamma only up to integration constants. Those constants feed into P_m and therefore into C_n. Confirm that your system's convention, such as Mathematica's PolyGamma[-k, z], is the one the constants were derived under. Then check one value of a numerically against direct quadrature before you trust the formula.
• Not fully simplified. The asker suspected that ψ^(-j-1)(a) - 2^j ψ^(-j-1)(a/2) collapses to terms in ψ^(-j-1)((a+1)/2), which the Ω_0 form hints at. Treat that as an open simplification, not an established identity.

For small a, the expression involves large cancellations between terms scaling like a^(-2n-1). Evaluate it in extended precision, not double.

This article draws on How to find $\int_0^\infty \frac{x^{2n} \arctan(x)}{e^{a\pi x}+1} \, dx $? and answers to it (1, 2), on the Stack Exchange network, licensed under CC BY-SA 4.0.