AbstractStability of the eigenfunctions of nonnegative selfadjoint second‐order linear elliptic operators subject to homogeneous Dirichlet boundary data under domain perturbation is investigated. Let Ω, \documentclass{article}\usepackage{amssymb}\begin{document}\pagestyle{empty}$\Omega ^{\prime } \subset \mathbb {R}^n$\end{document} be bounded open sets. The main result gives estimates for the variation of the eigenfunctions under perturbations Ω′ of Ω such that \documentclass{article}\usepackage{amssymb}\begin{document}\pagestyle{empty}$\Omega _{\varepsilon } = \lbrace x \in \Omega : {\rm dist}(x,\, \mathbb {R}^n \!\setminus \! \Omega ) > \varepsilon \rbrace \subset \Omega ^{\prime } \subset \overline{\Omega ^{\prime }} \subset \Omega$\end{document} in terms of powers of ε, where the parameter ε > 0 is sufficiently small. The estimates obtained here hold under some regularity assumptions on Ω, Ω′. They are obtained by using the notion of a gap between linear operators, which has been recently extended by the authors to differential operators defined on different open sets.