<p style='text-indent:20px;'>We expose here a novel application of the so-called coupled complex boundary method – first put forward by Cheng et al. (2014) to deal with inverse source problems – in the framework of shape optimization for solving the exterior Bernoulli problem, a prototypical model of stationary free boundary problems. The idea of the method is to transform the overdetermined problem to a complex boundary value problem with a complex Robin boundary condition coupling the Dirichlet and Neumann boundary conditions on the free boundary. Then, we optimize the cost function constructed by the imaginary part of the solution in the whole domain in order to identify the free boundary. We also prove the existence of the shape derivative of the complex state with respect to the domain. Afterwards, we compute the shape gradient of the cost functional, and characterize its shape Hessian at the optimal domain under a strong, and then a mild regularity assumption on the domain. We then prove the ill-posedness of the proposed shape problem by showing that the latter expression is compact. Also, we devise an iterative algorithm based on a Sobolev gradient scheme via finite element method to solve the minimization problem. Finally, we illustrate the applicability of the method through several numerical examples, both in two and three spatial dimensions.</p>