Summary: We prove the existence and asymptotic expansion of a large class of solutions to nonlinear Helmholtz equations of the form \((\Delta - \lambda^2) u = N[u]\), where \(\Delta = -\sum_j \partial^2_j\) is the Laplacian on \(\mathbb{R}^n, \lambda\) is a positive real number, and \(N[u]\) is a nonlinear operator depending polynomially on \(u\) and its derivatives of order up to order two. Nonlinear Helmholtz eigenfunctions with \(N[u]= \pm |u|^{p-1} u\) were first considered by S. Gutiérrez [Math. Ann. 328, No. 1–2, 1–25 (2004; Zbl 1109.35045)]. We show that for suitable nonlinearities and for every \(f \in H^{k+4}(\mathbb{S}^{n-1})\) of sufficiently small norm, there is a nonlinear Helmholtz function taking the form \(u(r, \omega) = r^{-(n-1)/2} (e^{-i\lambda r} f(\omega) + e^{+i\lambda r} b(\omega) + O(r^{-\epsilon}))\), as \(r \to \infty\), \( \epsilon > 0\), for some \(b \in H^k(\mathbb{S}^{n-1})\). Moreover, we prove the result in the general setting of asymptotically conic manifolds. The proof uses an elaboration of anisotropic Sobolev spaces defined by A. Vasy [Lond. Math. Soc. Lect. Note Ser. 443, 219–374 (2018; Zbl 1416.83026)], between which the Helmholtz operator \(\Delta - \lambda^2\) acts invertibly. These spaces have a variable spatial weight \(\mathsf{l}_\pm \), varying in phase space and distinguishing between the two “radial sets” corresponding to incoming oscillations, \(e^{-i\lambda r} \), and outgoing oscillations, \(e^{+i\lambda r} \). Our spaces have, in addition, module regularity with respect to two different “test modules” and have algebra (or pointwise multiplication) properties that allow us to treat nonlinearities \(N[u]\) of the form specified above.