  We analyze the theoretical properties of an adaptive Legendre-Galerkin method
in the multidimensional case. After the recent investigations for
Fourier-Galerkin methods in a periodic box and for Legendre-Galerkin methods in
the one dimensional setting, the present study represents a further step
towards a mathematically rigorous understanding of adaptive spectral/$hp$
discretizations of elliptic boundary-value problems. The main contribution of
the paper is a careful construction of a multidimensional Riesz basis in $H^1$,
based on a quasi-orthonormalization procedure. This allows us to design an
adaptive algorithm, to prove its convergence by a contraction argument, and to
discuss its optimality properties (in the sense of non-linear approximation
theory) in certain sparsity classes of Gevrey type.
