A brute-force search shows that the problem has no solution with at a least one tuple (xi,yi,zi) from the set {-N,...,N}, N = 40.
It shoud not be too hard to prove this in general (mathematically).

It's only a display decision when interface(typesetting=extended) .
For interface(typesetting=standard)  ==>  a^0.5

`@W`  is a name (symbol)  enclosed in back quotes. It could be used as a generic procedure, or maybe a procedure `@W` is defined somewhere. Such a name can be used in 2022 too!

See here the correct formula and its proof (it's very simple):

Proof of the Euler product formula for the Riemann zeta function - Wikipedia

It cannot be proved in Maple for obvious reasons.

It is easy to formulate the problem in Maple:

f:=x^(1/(x + 1/x)) + y^(1/(y + 1/y)) - exp(1);
F:=piecewise(f>0, 1, 0);
int(F, x=0..infinity, y=0..infinity); # No chance!

Even the approximation of the integral is not trivial, but it works.

Digits:=15:
evalf(Int(F, [x=2..7, y=2..7], method = _CubaCuhre, epsilon=1e-6));
#                        13.7592321343801

For a symbolic result we may try to exploit the symmetry, or to find a parametrization of the curve, but this is maths... 

eval([Pi/2<=argument(2*I*z), argument(2*I*z)<=Pi, 0<abs(4*I*z), abs(4*I*z)<=1], z=x+I*y); 
plots:-inequal(evalc(%), x = -1 .. 1, y = -1 .. 1, optionsopen = [linestyle = dash, color = red, thickness = 2]);

Remark.  Here argument(2*I*z)<=Pi  is superfluous. Note however that some authors/books use argument(z) in [0,2*Pi)  instead of (-Pi,Pi]. In such situations we have to adapt the maple code!

This is actually a visual representation of a Riemann surface for Ln (= multi-valued log)

col:=[green, cyan, red, maroon, grey]:
plots:-display( seq( plot3d([ r*cos(u), r*sin(u), u],
u=-Pi+2*k*Pi .. Pi+2*k*Pi, r=0..1,  labels=[Rez,Imz, 'Im(Ln(z))'], grid=[120,30], color=col[k+3], lightmodel=none), k=-2..2));

A few considerations.

The equality 
arctan(z/b)/b = int(1/(b^2+u^2), u = 0 .. z)     (1)

is well known when z, b are real and b<>0 (note that I changed the dummy variable, it's not a good idea to use the same letter here).
Your substitution does not make much sense.

Now about the complex case in (1), we must be careful with complex integration. Maple is used for formal (symbolic) computations but it is supposed that the user knows the maths behind these computations.

A complex integral such as (1) is usually a "complex path integral". The integrand is analytic (holomorphic} in C \ {+-I*b}. If the path is missing (as here), it usually means that the integral is independent of the path (i.e. depends only on the endpoints 0 and z). This is true (for our analytic function) if it is considered in a simply connected domain containing 0 and z but not containing the poles +-I*b. So, we must restrict the computations to such domains; for example, C \ {+-I*b} is not simply connected. Ignoring this aspect leads to wrong results or nonsenses.
The conclusion is that you must have knowledge in complex analysis to address such problems.

restart;
Vol:=(u,v,w,U,V,W)->sqrt(-u^2*U^4 + ((-W^2 + u^2 + v^2)*V^2 + (u^2 + w^2)*W^2 - u^4 + (v^2 + w^2)*u^2 - v^2*w^2)*U^2 - v^2*V^4 + ((v^2 + w^2)*W^2 + (v^2 - w^2)*u^2 - v^4 + v^2*w^2)*V^2 - W^2*(W^2*w^2 + (v^2 - w^2)*u^2 - v^2*w^2 + w^4))/12:
m:=20; nr:=0: nrmax:=50;
for U to m do
for V from U to m do
for W from V to U+V-1  do
for u to m do
for v from abs(W-u)+1 to W+u-1  do
for w from max(abs(W-u),abs(U-v))+1 to min(W+u,U+v)-1  do
  vv:=Vol(u,v,w,U,V,W);
  if vv::posint then nr++; lprint(nr, [u,v,w],[U,V,W],vol=vv)  fi;
  if nr>=nrmax then break 6 fi
od:od:od:od:od:od:

# Your example appears at the position 441  (for m=20)

                            m := 20

                          nrmax := 50

1, [8, 8, 7], [2, 3, 4], vol = 6
2, [17, 16, 16], [2, 3, 4], vol = 15
3, [6, 7, 8], [2, 4, 4], vol = 6
4, [6, 8, 7], [2, 4, 4], vol = 6
5, [7, 5, 6], [2, 4, 4], vol = 6
6, [7, 6, 5], [2, 4, 4], vol = 6
7, [4, 8, 7], [2, 5, 6], vol = 6
8, [7, 4, 4], [2, 5, 6], vol = 6
9, [16, 14, 13], [2, 5, 6], vol = 12
10, [4, 6, 5], [2, 7, 8], vol = 6
11, [6, 4, 4], [2, 7, 8], vol = 6
12, [8, 4, 3], [2, 7, 8], vol = 6
13, [12, 14, 13], [2, 7, 8], vol = 24
14, [10, 11, 12], [2, 8, 8], vol = 21
15, [10, 12, 11], [2, 8, 8], vol = 21
16, [11, 9, 10], [2, 8, 8], vol = 21
17, [11, 10, 9], [2, 8, 8], vol = 21
18, [17, 11, 12], [2, 8, 8], vol = 21
19, [17, 12, 11], [2, 8, 8], vol = 21
20, [8, 12, 11], [2, 9, 10], vol = 21
21, [11, 8, 8], [2, 9, 10], vol = 21
22, [8, 10, 9], [2, 11, 12], vol = 21
23, [10, 8, 8], [2, 11, 12], vol = 21
24, [17, 8, 8], [2, 11, 12], vol = 21
25, [12, 8, 7], [2, 13, 14], vol = 24
26, [16, 6, 5], [2, 13, 14], vol = 12
27, [16, 22, 21], [2, 13, 14], vol = 60
28, [17, 3, 4], [2, 16, 16], vol = 15
29, [17, 4, 3], [2, 16, 16], vol = 15
30, [12, 19, 19], [2, 19, 19], vol = 72
31, [8, 24, 23], [2, 19, 20], vol = 42
32, [17, 16, 17], [3, 5, 7], vol = 30
33, [4, 7, 8], [3, 6, 6], vol = 9
34, [4, 8, 7], [3, 6, 6], vol = 9
35, [14, 16, 17], [3, 6, 6], vol = 36
36, [14, 17, 16], [3, 6, 6], vol = 36
37, [4, 6, 6], [3, 7, 8], vol = 9
38, [7, 11, 10], [3, 7, 8], vol = 24
39, [8, 4, 2], [3, 7, 8], vol = 6
40, [11, 7, 8], [3, 7, 8], vol = 24
41, [16, 16, 17], [3, 7, 8], vol = 48
42, [16, 8, 10], [3, 8, 10], vol = 21
43, [11, 15, 16], [3, 9, 9], vol = 42
44, [11, 16, 15], [3, 9, 9], vol = 42
45, [7, 8, 7], [3, 10, 11], vol = 24
46, [9, 15, 14], [3, 11, 12], vol = 48
47, [8, 14, 13], [3, 13, 14], vol = 48
48, [9, 12, 11], [3, 14, 15], vol = 48
49, [9, 17, 15], [3, 15, 16], vol = 42
50, [11, 9, 9], [3, 15, 16], vol = 42

Most "complex plot" commands are very simple, but are convenient. Basically they just take the real and imaginary parts and call a regular plot.
For example, for f:=sin(x+I):

complexplot(f,x=-Pi..Pi);

is equivqlent to:

plot([Re(f),Im(f),x=-Pi..Pi]);

restart;
i2 := (x,y) -> -(1/2)*I*(exp(I*x)*(sin(x)/x)-exp(I*y)*(sin(y)/y))/(x-y);
i3_r := -(1/2)*I*(i2(y,z)-i2(y,x))/(z-x);
normal(eval(series(eval(i3_r,[x=t*x,y=t*y,z=t*z]),t=0),t=1));

On my rather old computer, for any of your graphs (having 21 vertices and 76 edges), VertexConnectivity needs about 0.5 sec.
So, it is easy to estimate the total time for your almost 15000 graphs.

It should be easy to write a compiled version for VertexConnectivity. It will be much faster (and a skilled user can do it), but probably the designers did not anticipate that someone will need it for so many graphs.

`simplify/size` is responsible for this. It calls `simplify/size/size` to compute the "math complexity" (see ?simplify/size) and chooses the simpler one. 

ex:=[t^2*x^2+t^2*y^2, t^2*x^2-t^2*y^2]:
`simplify/size/size`~(%);
#         [56, 59]

simplify~(ex);
#         [t^2*(x^2+y^2), t^2*x^2-t^2*y^2]

`simplify/size/size`~(%);
#        [38, 59]

Note that simplify(..., size=false) keeps both expressions unchanged.

If you want the rational exponent of x^r, where x is a variable and r is rational (including r=1)
then r=1 must be treated separately.

Exponent := proc(x::{name, name ^ rational}) 
  if x::name then 1 else op(2, x) fi 
end proc:
Exponent(x^(-3/2));   # -3/2
Exponent(Y[1,3]);     # 1

Note that r=0 is not possible because x^0 is automatically simplified to 1.