  Superpositions of plane waves are known to approximate well the solutions of
the Helmholtz equation. Their use in discretizations is typical of Trefftz
methods for Helmholtz problems, aiming to achieve high accuracy with a small
number of degrees of freedom. However, Trefftz methods lead to ill-conditioned
linear systems, and it is often impossible to obtain the desired accuracy in
floating-point arithmetic. In this paper we show that a judicious choice of
plane waves can ensure high-accuracy solutions in a numerically stable way, in
spite of having to solve such ill-conditioned systems.
  Numerical accuracy of plane wave methods is linked not only to the
approximation space, but also to the size of the coefficients in the plane wave
expansion. We show that the use of plane waves can lead to exponentially large
coefficients, regardless of the orientations and the number of plane waves, and
this causes numerical instability. We prove that all Helmholtz fields are
continuous superposition of evanescent plane waves, i.e., plane waves with
complex propagation vectors associated with exponential decay, and show that
this leads to bounded representations. We provide a constructive scheme to
select a set of real and complex-valued propagation vectors numerically. This
results in an explicit selection of plane waves and an associated Trefftz
method that achieves accuracy and stability.
  The theoretical analysis is provided for a two-dimensional domain with
circular shape. However, the principles are general and we conclude the paper
with a numerical experiment demonstrating practical applicability also for
polygonal domains.
