We study the optimization of Steklov eigenvalues with respect to a boundary density function \(\rho \) on a bounded Lipschitz domain \(\Omega \subset \mathbb {R}^N\) . We investigate the minimization and maximization of \(\lambda _k(\rho )\) , the kth Steklov eigenvalue, over admissible densities satisfying pointwise bounds and a fixed integral constraint. Our analysis covers both first and higher-order eigenvalues and applies to domains \(\Omega \) with general geometry and topology. We establish the existence of optimal solutions and provide structural characterizations: minimizers are bang--bang functions and may have disconnected support, while maximizers are not necessarily bang--bang. On circular domains, the minimization problem admits infinitely many minimizers generated by rotational symmetry, while the maximization problem has infinitely many distinct maximizers that are not symmetry-induced. We also show that the maps \(\rho \mapsto \lambda _k(\rho )\) and \(\rho \mapsto 1/\lambda _k(\rho )\) are generally neither convex nor concave, limiting the use of classical convex optimization tools. To address these challenges, we analyze the objective functional and introduce a Fréchet differentiable surrogate that enables the derivation of optimality conditions. We further design an efficient numerical algorithm, with experiments illustrating the difficulty of recovering optimal densities when they lack smoothness or exhibit oscillations.
