I will first talk about the way I understand the DFT.
Let $u:[0,2\pi)\to \mathbb{R}$ and consider the equally spaced sample of $u$ (sample size $N$): $\{u[0],u[1], \cdots u[N-1]\}$ where $u[n] = u(n\Delta x)$ for $n=0,1,\cdots N-1$ and $\Delta x = \frac{2\pi}{N}$. The Fourier coefficient of $u$ is defined by
$\displaystyle \hat{u}[k] = \frac{1}{2\pi}\int_{0}^{2\pi}u(x)e^{-ikx}dx\tag*{}, k\in \mathbb{Z}$
We can approximate the Fourier coefficient of $u$ by discretizing the integral:
$\displaystyle \hat{u}[k] \sim \frac{1}{2\pi} \sum_{n=0}^{N-1} u(n\Delta x)e^{-ikn\Delta x}\Delta x\tag*{}$
Substituting the $u[n]=u(n\Delta x), \Delta x = 2\pi/N$ would get us
$\displaystyle \hat{u}[k] \sim \frac{1}{N} \sum_{n=0}^{N-1} u[n]e^{-i\frac{2\pi kn}{N} }\tag*{}$
The right-hand side is called the discrete Fourier transform of $u$. Assuming that the $\hat{u}[k]$ vanishes for large enough $k$, we can assume that the range of the index $k$ is $k=-N/2, \cdots N/2-1$ for large enough even $N$.
Main Question.
Now let's say I want to approximate the Fourier coefficient of $u^2$ from the DFT sequence $\{\hat{u}[-N/2],\cdots,\hat{u}[N/2-1]\}$. Strictly, it is the discrete convolution:
$\displaystyle \hat{u^2}[k] = \sum_{j=-\infty}^{\infty}\hat{u}[j]\hat{u}[k-j] \tag*{}$
Using the assumption that $\hat{u}[k]$ vanishes outside of $-N/2, \cdots N/2-1$, $j$ should satisfy $-N/2\leq j\leq N/2-1$ and $-N/2\leq k-j \leq N/2-1$( or $k+1-N/2\leq j\leq k+N/2$). So
$\displaystyle \hat{u^2}[k] = \sum_{j=\textrm{max}\{-N/2,\textrm{ }k+1-N/2\}}^{\textrm{min}\{N/2-1,\textrm{ }k+N/2\}}\hat{u}[j]\hat{u}[k-j] \tag*{}$
which looks too complicated. I would imagine that this would be equal to the circular convolution, but I am not sure how to prove it.