We deal with a regularity theorem for the free boundary of solutions to shape optimization problems with integral functionals, for which the energy of a domain $\Omega$ is obtained as the integral of a cost function $j(u,x)$ depending on the solution $u$ of a certain PDE problem on $\Omega$. The main feature of these functionals is that (in general) the minimality of a domain $\Omega$ cannot be translated into a variational problem for a single (real or vector valued) state function, which is the starting point in all the previous regularity results for optimal shapes.

In this paper we focus on the case of affine cost functions $j(u,x)=-g(x)u+Q(x)$, where $u$ is the solution of the PDE $-\Delta u=f$ with Dirichlet boundary conditions. We obtain the Lipschitz continuity and the non-degeneracy of the optimal $u$ from the inwards/outwards optimality of $\Omega$ and then we use the stability of $\Omega$ with respect to variations with smooth vector fields in order to study the blow-up limits of the state function $u$.

By performing a triple consecutive blow-up, we prove the existence of blow-up sequences converging to homogeneous stable solution of the one-phase Bernoulli problem and according to the blow-up limits, we decompose $\partial\Omega$ into a singular and a regular part. In order to estimate the Hausdorff dimension of the singular set of $\partial\Omega$ we give a new formulation of the notion of stability for the one-phase problem, which is preserved under blow-up limits and allows to develop a dimension reduction principle. Finally, by combining a higher order Boundary Harnack principle and a viscosity approach, we prove $C^\infty$ regularity of the regular part of the free boundary (if the data are smooth).