only on bisecting y=0 plane


Integrate[1/Sqrt[(x-x1)^2+y1^2+z^2]^3,{y1,-a,a},Assumptions->Element[{x,x1,z,a},Reals]&&a>0]
%[[1]]
z Integrate[%,{x1,-a,a},Assumptions->Element[{x,z,a},Reals]&&a>0]

                            a (a - x)
Out[4]= 2 (ArcTan[------------------------------] + 
                            2            2    2
                  z Sqrt[2 a  - 2 a x + x  + z ]
 
                        a (a + x)
>      ArcTan[------------------------------])
                        2            2    2
              z Sqrt[2 a  + 2 a x + x  + z ]

omega=%[[1]] /. a->1
Plot[Evaluate[omega /. x->0],{z,-2,2}]
Plot[Evaluate[omega /. x->2],{z,-2,2}]
GraphicsArray[{%,%%}]
Show[%]

ContourPlot[omega, {x, -2, 2}, {z, .001, 2},PlotPoints->100,ContourShading ->False,Contours -> Table[i/2,{i,1,12}],PlotRangePadding ->None ,AspectRatio->Automatic] 

b={-D[omega,x],-D[omega,z]}
Show[VectorPlot[b,{x,-2,2},{z,.001,2},VectorPoints->{12,5},VectorScale->{.5,.1,Automatic},RegionFunction -> Function[{x, z, vx, vy, n}, Min[(x-1)^2,(x+1)^2]+z^2>.2],PlotRange->{{-2,2},{0,2}},AspectRatio->Automatic]]
Show[%,%%%]
GraphicsGrid[{{%},{%7}}]

along diagonal

Integrate[1/Sqrt[(x-x1)^2+y1^2+z^2]^3,{y1,-(d-x1),(d-x1)},Assumptions->Element[{x,x1,z,d},Reals]&&d>0]
%[[1]]
z Integrate[%,{x1,0,d},Assumptions->Element[{x,z,d},Reals]&&d>0]
omega1=%[[1]] /. d->Sqrt[2]
Integrate[1/Sqrt[(x-x1)^2+y1^2+z^2]^3,{y1,-(d+x1),(d+x1)},Assumptions->Element[{x,x1,z,d},Reals]&&d>0]
%[[1]]
z Integrate[%,{x1,-d,0},Assumptions->Element[{x,z,d},Reals]&&d>0]
omega2=%[[1]] /. d->Sqrt[2]
omega0=Mod[Re[omega1+omega2], 4 Pi]
I*(Log[(-2*z*(z - I*Sqrt[2 - 2*Sqrt[2]*x + x^2 + z^2]))/
    (-Sqrt[2] + x + I*z)^2] - 
  Log[(-2*z*(z + I*Sqrt[2 - 2*Sqrt[2]*x + x^2 + z^2]))/
    (Sqrt[2] - x + I*z)^2] - 
     Log[(I*z*(2 + 2*Sqrt[2]*x + x^2 + z^2 - I*z*Sqrt[2 + 2*Sqrt[2]*x + x^2 + z^2]))/
(Sqrt[2] + x - I*z)^2] + 
  Log[((-I)*z*(2 + 2*Sqrt[2]*x + x^2 + z^2 +  I*z*Sqrt[2 + 2*Sqrt[2]*x + x^2 + z^2]))/
(Sqrt[2] + x + I*z)^2] - 
       Log[(z*(2*I + (x + I*z)*(z - I*Sqrt[2 + x^2 + z^2]) -  I*Sqrt[2]*(-x - I*z + Sqrt[2 + x^2 + z^2])))/
    ((x + I*z)*(Sqrt[2] + x + I*z)^2)] - 
Log[((-2*I)*z*(2 + Sqrt[2]*(-x - I*z + Sqrt[2 + x^2 + z^2]) -  (x + I*z)*(I*z + Sqrt[2 + x^2 + z^2])))/
    ((x + I*z)*(-Sqrt[2] + x + I*z)^2)] + 
       Log[(z*(-2*I + (x - I*z)*(z + I*Sqrt[2 + x^2 + z^2]) +  I*Sqrt[2]*(-x + I*z + Sqrt[2 + x^2 + z^2])))/
    ((x - I*z)*(Sqrt[2] + x - I*z)^2)] + 
  Log[((2*I)*z*(2 - (x - I*z)*((-I)*z + Sqrt[2 + x^2 + z^2]) +  Sqrt[2]*(-x + I*z + Sqrt[2 + x^2 + z^2])))/
    ((x - I*z)*(Sqrt[2] - x + I*z)^2)])

ContourPlot[omega0, {x, -2, 2}, {z, .001, 2},PlotPoints->100,ContourShading ->False,Contours -> Table[i/2,{i,1,12}],PlotRangePadding ->None ,AspectRatio->Automatic] 
