Multi-dimensional integral in the complex plane with poles and essential singularity

Question

I've spent the last week searching for a way to numerically calculate the pole contributions to this multi-dimensional complex contour integral while avoiding the singularity at z=0: $$ \oint_{C_1(0)}\oint_{C_p}\frac{z_{1}z_{2}}{z_{2}^{-1}p+(1-p)z_{1}-1} \, e^{\frac{1}{z_{1}} + z_{1} + \frac{1}{z_{2}} + z_{2}} \quad \text d z_1 \, \text d z_2 $$ where $p \in (0,1)$, $C_1(0)$ is a circle of radius 1 centered at 0 and $C_p$ is a contour surrounding the poles only.

I used Finding residues of multi-dimensional complex functions to compute this integral. When I give the whole integral to Mathematica, I don't have the contribution of the poles because the contour cannot be deformed due to the singularity. So I'm trying to find a way to exclude the singularity of the contour, maybe with a small contour around $z = 0$ using Piecewise or something similar to get only the contribution of the poles.

Code

NContourIntegrate[f_, par : (z_ -> g_), {t_, a_, b_}] := 
  NIntegrate[Evaluate[D[g, t] (f /. par) /. t -> t1], {t1, a, b}]  

Clear[Pinfz1];  
Pinfz1[p_?NumericQ, z2_?NumericQ] := 
  NContourIntegrate[1/(2*I*Pi)^2*z1*z2/(p + (1 - p)*z1*z2-z2)*E^(1/z1 + z1 + 1/z2 + z2), 
   z1 -> Exp[I t], 
   {t, 0, 2*Pi}]  

Clear[Pinfz];  
Pinfz[p_?NumericQ] := NContourIntegrate[Pinfz1[p, z2], z2 -> Exp[I t], {t, 0, 2*Pi }]  

Answer 1

I'm not familiar with the physics so I cannot say whether integrating over this set of 2 real dimension, in C^2, is what is wanted. I think the code below will cover the product space of the contours that are requested.

ii = z1*z2/(p/z2 + (1 - p) z1 - 1)*Exp[1/z1 + z1 + 1/z2 + z2];
i1 = (ii /. {z1 -> Exp[I*t1], z2 -> Exp[I*t2]})*I*Exp[I*t1]*I*
   Exp[I*t2];
f[pval_] := 
 Quiet[NIntegrate[
   Evaluate[i1 /. p -> pval], {t1, 0, 2*Pi}, {t2, 0, 2*Pi}]]
f[1/2]
f[1/3]

(* Out[114]= 89.3519154317 - 5.40968238685*10^-8 I

Out[115]= 62.3532282953 + 0.000124769719229 I *)

Also:

f[.1]

(* Out[116]= 34.4365223939 + 0.000414467508089 I *)

f[.9]

(* Out[117]= 170.924622109 + 0.00264992865898 I *)

Offhand I cannot say much about how accurate these might be. When I did them at higher working precision I got consistent results, but that does not prove anything.