Unexpected Behavior from Integrate with PrincipalValue -> True

Question

Bug introduced in 9.0 or earlier and persisting through 12.1 or later. Reported to Wolfram, Inc. as CASE:4544446.

The following code, taken from my answer to 97024,

Integrate[1/((1 - 4 ωp^2 - 4 I ωp γ) (1 - ωp^2 - 2 I ωp γ)^2 (ωp - ω)), 
    {ωp, -Infinity, Infinity}, Assumptions -> γ > 0 && ω > 0, PrincipalValue -> True]

runs for 41 minutes on Mathematica

$Version
(* 10.2.0 for Microsoft Windows (64-bit) (July 7, 2015) *)

before returning unevaluated, in the process twice giving the two error messages,

Limit::cas: Warning: contradictory assumption(s) False&&ωp>4096 encountered. >>
Integrate::idiv: Integral of ... does not converge on does not converge on {-∞,∞}. >>

neither of which is valid. Omitting the term, (1 - ωp^2 - 2 I ωp γ)^2 or setting ω -> 0 produces the same outcome, only faster. In fact, the only variant I have been able to make work is to set γ to a numerical value, here 3.

Integrate[1/((1 - 4 ωp^2 - 4 I ωp 3) (1 - ωp^2 - 2 I ωp 3)^2 (ωp - ω)), 
    {ωp, -Infinity, Infinity}, Assumptions -> ω > 0, PrincipalValue -> True]
(* -((I π)/((-1 + 4 ω (3 I + ω)) (-1 + ω (6 I + ω))^2)) *)

Is there some work-around for this problem, not involving evaluation of residues?

Issues persist in 10.3

With

$Version
(* 10.3.0 for Microsoft Windows (64-bit) (October 9, 2015) *)

a simpler version of the code above,

Integrate[1/((1 - 4 ωp^2 - 4 I ωp γ) (ωp - ω)), {ωp, -Infinity, Infinity}, 
    Assumptions -> γ > 0 && ω > 0, PrincipalValue -> True]

produces the same spurious error messages as before, although it takes twice as long to do so (about eight minutes). Incidentally, a still faster way to reproduce this problem is

Integrate[1/(-(ωp - 1/2 (-I γ - Sqrt[1 - γ^2])) (ωp - ω)), {ωp, -∞, ∞}, 
    Assumptions -> γ > 0 && ω ∈ Reals, PrincipalValue -> True]

which takes only seconds. Note that I also replaced ω > 0 by ω ∈ Reals to see whether it mattered. It does not.

Issues persist for 12.1

Mathematica 12.1 now can solve the first integral given under "Issues persist in 10.3" (in just over two minutes), but the answer is wrong.

(* ConditionalExpression[(Pi*(I*γ - Sqrt[1 - γ^2] + 2*ω))/(Sqrt[1 - γ^2]*
   (I + 4*(γ - I*ω)*ω)), γ < 1] *)

as can be seen from summing the residues,

Residue[1/((1 - 4 ωp^2 - 4 I ωp γ) (ωp - ω)), {ωp, #}] & /@ 
    (ωp /. Solve[((1 - 4 ωp^2 - 4 I ωp γ) (ωp - ω)) == 0, ωp])
(* {-(1/(2 Sqrt[1 - γ^2] (I γ + Sqrt[1 - γ^2] + 2 ω))),
    -(1/(2 Sqrt[1 - γ^2] (-I γ + Sqrt[1 - γ^2] - 2 ω))), 
    I/(I + 4 γ ω - 4 I ω^2)} *)

Simplify[2 π I (Total[%] - Last[%]/2)]
(* π/(I + 4 γ ω - 4 I ω^2) *)

The other integrals given in this question either return unevaluated with error messages or run indefinitely.

Answer 1

Interesting exercise. It can be done with some preparatory steps in a few seconds of computing time.

The OP asks for the integral

q := Integrate[h, {ωp, -∞, ∞}, 
  Assumptions -> {γ > 0, ω > 0}, PrincipalValue -> True]

where the integrand is

h = 1/((1 - 4 ωp^2 - 
     4 I ωp γ) (1 - ωp^2 - 
      2 I ωp γ)^2 (ωp - ω));

Compact answer

Letting

g = 1/((x - a) (x - b) (x - c) (x - d)^2 (x - f)^2);
rep = Thread[{b, c, d, d, f, f, a} -> (ωp /. Solve[0 == 1/h, ωp])];

we have

Simplify[h == (-1/4) g /. rep /. x -> ωp]

(*
Out[145]= True
*)

and the integral becomes

q = Timing[(-1/4) Simplify[
    Plus @@ Limit[
       Integrate[List @@ Apart[g], {x, -z, z}, PrincipalValue -> True, 
        Assumptions -> {z > 0, a ∈ Reals, Im[b] != 0, Im[c] != 0, 
          Im[d] != 0, Im[f] != 0}], z -> ∞] /. rep]]
{14.2116911`, -((
  2 π)/((I + 4 γ ω - 4 I ω^2) (2 γ ω - 
     I (-1 + ω^2))^2))}

and it is done in less that 15 seconds (Version 10.1. on a normal Windows Laptop).

Detailed derivation

We show here that the integral can in fact be calculated refraining from the use of evaluating residues. The main idea is to understand the integral as a limit z->∞ of the finite integral between -z and z. Also, an explicit partial fraction decomposition using Apart[] helps reduce the computing time.

Determining explcitly the roots of the denominator

sol = Solve[0 == 1/h, ωp]

(*
Out[70]= {{ωp -> 1/2 (-I γ - Sqrt[1 - γ^2])}, {ωp -> 
   1/2 (-I γ + Sqrt[1 - γ^2])}, {ωp -> -I γ - 
    I Sqrt[-1 + γ^2]}, {ωp -> -I γ - 
    I Sqrt[-1 + γ^2]}, {ωp -> -I γ + 
    I Sqrt[-1 + γ^2]}, {ωp -> -I γ + 
    I Sqrt[-1 + γ^2]}, {ωp -> ω}}
*)

Letting

p = ωp /. sol

(*
Out[88]= {1/2 (-I γ - Sqrt[1 - γ^2]), 
 1/2 (-I γ + Sqrt[1 - γ^2]), -I γ - 
  I Sqrt[-1 + γ^2], -I γ - 
  I Sqrt[-1 + γ^2], -I γ + 
  I Sqrt[-1 + γ^2], -I γ + I Sqrt[-1 + γ^2], ω}
*)

We find

Simplify[4 Times @@ (ωp - p) == -1/h]

(*
Out[89]= True
*)

Assuming general parameters defined in this replacement

rep = Thread[{b, c, d, d, f, f, a} -> p]

(*
Out[90]= {b -> 1/2 (-I γ - Sqrt[1 - γ^2]), 
 c -> 1/2 (-I γ + Sqrt[1 - γ^2]), 
 d -> -I γ - I Sqrt[-1 + γ^2], 
 d -> -I γ - I Sqrt[-1 + γ^2], 
 f -> -I γ + I Sqrt[-1 + γ^2], 
 f -> -I γ + I Sqrt[-1 + γ^2], a -> ω}
*)

We shall consider the integral of the more general expression for the integrand (for ease of reading we let ωp -> x)

g = 1/((x - a) (x - b) (x - c) (x - d)^2 (x - f)^2);

g /. rep /. x -> ωp // Simplify

(*
Out[97]= -(4/((ω - ωp) (-1 + 
    2 I γ ωp + ωp^2)^2 (-1 + 4 I γ ωp + 
    4 ωp^2)))
*)

The integral over g is then -4 times the integral in question, i.e. q = -g/4.

Before taking the integral It helps decomposing g into partial fractions, and splitting it into a list

ga = List @@ Apart[g]

(*
Out[107]= {1/((a - b) (a - c) (a - d)^2 (a - f)^2 (-a + x)), -(
  1/((a - b) (b - c) (b - d)^2 (b - f)^2 (-b + x))), -(
  1/((a - c) (-b + c) (c - d)^2 (c - f)^2 (-c + x))), -(
  1/((a - d) (-b + d) (-c + d) (d - f)^2 (-d + x)^2)), (
 2 a b c - 3 a b d - 3 a c d - 3 b c d + 4 a d^2 + 4 b d^2 + 4 c d^2 - 5 d^3 +
   a b f + a c f + b c f - 2 a d f - 2 b d f - 2 c d f + 
  3 d^2 f)/((a - d)^2 (-b + d)^2 (-c + d)^2 (d - f)^3 (-d + x)), -(
  1/((b - f) (-a + f) (-c + f) (-d + f)^2 (-f + x)^2)), (
 2 a b c + a b d + a c d + b c d - 3 a b f - 3 a c f - 3 b c f - 2 a d f - 
  2 b d f - 2 c d f + 4 a f^2 + 4 b f^2 + 4 c f^2 + 3 d f^2 - 
  5 f^3)/((a - f)^2 (-b + f)^2 (-c + f)^2 (-d + f)^3 (-f + x))}
*)

Now the integral over the list ga is quickly calculated

gai = Timing[
  Integrate[ga, {x, -z, z}, PrincipalValue -> True, 
   Assumptions -> {z > 0, a ∈ Reals, Im[b] != 0, Im[c] != 0, 
     Im[d] != 0, Im[f] != 0}]]

(*
Out[110]= {15.8029, {(-Log[-a - z] + 
   Log[-a + z])/((a - b) (a - c) (a - d)^2 (a - f)^2), (-Log[b - z] + 
   Log[b + z])/((a - b) (b - c) (b - d)^2 (b - f)^2), (-Log[c - z] + 
   Log[c + z])/((a - c) (-b + c) (c - d)^2 (c - f)^2), -((
   2 z)/((a - d) (-b + d) (-c + d) (d - f)^2 (d^2 - 
      z^2))), ((d (-3 b c + 4 (b + c) d - 5 d^2) + (b c - 2 (b + c) d + 
        3 d^2) f + 
     a (2 b c - 3 b d - 3 c d + 4 d^2 + (b + c - 2 d) f)) (Log[d - z] - 
     Log[d + z]))/((a - d)^2 (b - d)^2 (c - d)^2 (d - f)^3), -((
   2 z)/((b - f) (d - f)^2 (-a + f) (-c + f) (f^2 - z^2))), ((2 a b c + 
     a b d + a c d + 
     b c d - (3 b c + 3 a (b + c) + 2 (a + b + c) d) f + (4 (a + b + c) + 
        3 d) f^2 - 5 f^3) (Log[f - z] - Log[f + z]))/((a - f)^2 (b - f)^2 (c -
      f)^2 (-d + f)^3)}}
*)

Now taking the limit z->∞

gail = Timing[Limit[gi[[2]], z -> ∞]]

(*
Out[111]= {0.265202, {-((I π)/((a - b) (a - c) (a - d)^2 (a - f)^2)), -((
   I π)/((a - b) (b - c) (b - d)^2 (b - f)^2)), -((
   I π)/((a - c) (-b + c) (c - d)^2 (c - f)^2)), 0, (
  I (d (-3 b c + 4 (b + c) d - 5 d^2) + (b c - 2 (b + c) d + 3 d^2) f + 
     a (2 b c - 3 b d - 3 c d + 4 d^2 + (b + c - 2 d) f)) π)/((a - 
     d)^2 (b - d)^2 (c - d)^2 (d - f)^3), 0, (
  I (2 a b c + a b d + a c d + 
     b c d - (3 b c + 3 a (b + c) + 2 (a + b + c) d) f + (4 (a + b + c) + 
        3 d) f^2 - 5 f^3) π)/((a - f)^2 (b - f)^2 (c - f)^2 (-d + f)^3)}}
*)

and summing up the terms

gs = Plus @@ gail[[2]]

(*
Out[117]= -((I π)/((a - b) (a - c) (a - d)^2 (a - f)^2)) - (
 I π)/((a - b) (b - c) (b - d)^2 (b - f)^2) - (
 I π)/((a - c) (-b + c) (c - d)^2 (c - f)^2) + (
 I (2 a b c + a b d + a c d + 
    b c d - (3 b c + 3 a (b + c) + 2 (a + b + c) d) f + (4 (a + b + c) + 
       3 d) f^2 - 5 f^3) π)/((a - f)^2 (b - f)^2 (c - f)^2 (-d + f)^3) + (
 I (d (-3 b c + 4 (b + c) d - 5 d^2) + (b c - 2 (b + c) d + 3 d^2) f + 
    a (2 b c - 3 b d - 3 c d + 4 d^2 + (b + c - 2 d) f)) π)/((a - 
    d)^2 (b - d)^2 (c - d)^2 (d - f)^3)
*)

Replacing the parameters

gsrep = gs /. rep // Simplify

(*
Out[120]= (8 π)/((I + 4 γ ω - 4 I ω^2) (2 γ ω - 
   I (-1 + ω^2))^2)
*)

Hence the original integral is simply

q = -1/4 gsrep

$$q = -\frac{2 \pi }{\left(4 \gamma \omega -4 i \omega ^2+i\right) \left(2 \gamma \omega -i \left(\omega ^2-1\right)\right)^2}$$

Done.