Uplink grant-free non-orthogonal multiple access (NOMA) is a promising technology for massive connectivity with low latency and high energy efficiency. In code-domain NOMA schemes, the requirements boil down to the design of codebooks that contain a large number of spreading sequences with low peak-to-average power ratio (PAPR) while maintaining low coherence. When employing binary Golay sequences with guaranteed low PAPR in the design, the fundamental problem is to construct a large set of $n$-variable quadratic bent or near-bent functions in a particular form such that the difference of any two is bent for even $n$ or near-bent for odd $n$ to achieve optimally low coherence. In this work, we propose a theoretical construction of NOMA codebooks by applying a recursive approach to those particular quadratic bent functions in smaller dimensions. The proposed construction yields desired NOMA codebooks that contain $6\cdot N$ Golay sequences of length $N=2^{4m}$ for any positive integer $m$ and have the lowest possible coherence $1/\sqrt{N}$.