(EGJK) for compact objects of arbitrary shape and separation and applied to predict the force between a cylinder and a plate [10], and between two spheres [9]. EGJK’s agorithm is efficient to predict the Casimir energy between compact 3D objects of spheroidal or nearly- spheroidal shape. The algorithm leads to analytically tractable series solution. But it is not of practical use for general geometries due to the poor convergence rate for objects with corners. This limitation occurs since they used spherical basis functions to expand the current distribution. It could be overcome by expanding the currents with basis functions defined on a pair of adjacent triangles, which is known as the RWG basis [11, 12]. The RWG basis is vastly popular in the surface integral method in computational electromagnetics. After some manipulation, the evaluation of Casimir energy depends only on solving for all the eigenvalues of the impedance matrix obtained from the integral equation by the method of moments. 

The theoretical foundation for the Maxwell stress tensor approach was built by Lifshitz and his coworkers [13]: the net Casimir force on a body can be expressed as an integral over any closed surface around the body of the mean electromagnetic stress tensor, inte- grated over all frequencies. And the fluctuation-dissipation theorem states that the mean fluctuating field, which forms the Maxwell stress tensor, is related to the imaginary part of the dyadic Green’s function of the same geometry. The dyadic Green’s function can be evaluated by existing numerical methods in computational electromagnetics (electrodynam- ics) [14]. Rodriguez has demonstrated the applicability of this method by using a simple finite-difference frequency-domain (FDFD) method to calculate the dyadic Green’s func- tion in z -invariant structures [15, 16]. This approach relates the quantum electrodynamic phenomena with classical electromagnetics which motivates our piece of work: to evaluate the Casimir force among arbitrary objects using integral equation method. With surface integral equation method, the number of unknowns can be greatly reduced. Moreover, fast algorithms [17] are available to further improve the efficiency of the method. The compu- tational complexity of both the stress tensor approach and the path integral approach can be reduced to O ( N log N ) times the number of iterations at best, N being the number of surface unknowns. Thus their efficiency is comparable. The main difference between them is the physical quantity being calculated directly: path integral method calculates the energy while stress tensor approach calculates the force distribution (pressure). 

In a paper published earlier [20], we have demonstrated that surface integral equation 

