In a one-electron system such as the H atom, we are able to say that the orbitals are valid wavefunctions. But for complex systems with N-electrons, I have heard people say that they are just mathematical tools to help us build the wavefunction.
Thus, I am asking whether this complexity is a result of the electron-electron interaction in the N-electron system?
For a 1-electron system, wavefunctions and orbitals are the same thing. For a 2-electron system, or any N-electron system with N>1, we can use "orbitals" to build the wavefunction, which you described as an "added complexity".
The added "complexity" of using "orbitals" to build a wavefunction is introduced because we cannot analytically solve the Schroedinger equation for 2-electron systems. Your question therefore seems to implicitly be: Is the reason why we can't solve the SE analytically for 2-electron systems, because of the electron-electron interaction?
The Hartree-Fock approximate solution to a 2-electron Schroedinger equation is basically the solution that ignores electron-electron correlation, and yet we still use orbitals to build the wavefunction. Therefore the extra "complexity" is still there when we ignore electron-electron correlation. However, electron-electron correlation is more complex than just any electron-electron interaction, so perhaps you wonder, if we were to remove all electron-electron interactions, we wouldn't need to use orbitals anymore.
Let's look at another very simple system for which the Schroedinger equation cannot be solved analytically: $\ce{H2+}$. This is a 1-electron system, so the "electron-electron interaction" is zero, and yet we still cannot solve the Schroedinger equation analytically. However, despite needing to solve the Schroedinger equation for this 1-electron system using approximations, the added "complexity" of building a wavefunctions from orbitals is probably not appropriate because the system only has one electron anyway! Here I suppose the exact wavefunctions can be called "molecular orbitals" with a meaning different from the LCAO definition.
To simplify things, let's stick to having just one nucleus and two electrons (He-like systems). The Hamiltonian would be (if you will allow me to use atomic units, which I don't usually do):
$$\tag{1} -\frac{1}{2}\nabla^2_1 - \frac{1}{2}\nabla^2_2 - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{r_{12}}. $$
The electron-electron interaction is the $\frac{1}{r_{12}}$ term. Can we solve the Schroedinger equation exactly, without using orbitals, if we remove that term? Let's see:
$$\tag{1} -\frac{1}{2}\nabla^2_1 - \frac{1}{2}\nabla^2_2 - \frac{2}{r_1} - \frac{2}{r_2}, $$
This looks like it can be solved using separation of variables, but it's actually the same as the H2+ Hamiltonian but with one more kinetic energy term, making the solution even more complicated than what we'd get for H2+.
The electron-electron interaction term seems to be blamed here (although without justification/explanation):
"Unfortunately, the Coulomb repulsion terms make it impossible to find an exact solution to the Schrödinger equation for many-electron atoms and molecules even if there are only two electrons."
Also here, it credits the lack of the electron-electron repulsion term for being able to obtain an eigenvalue (not necessarily eigenfunction) of the H2+ Schroedinger equation:
"The energy E is the eigenvalue of the Schrödinger equation for the single electron. The equation can be solved in a relatively straightforward way due to the lack of electron–electron repulsion (electron correlation).
But then it goes on to say in the next sentence that the wavefunction would still be quite complicated:
"The wave equation (a partial differential equation) separates into two coupled ordinary differential equations when using prolate spheroidal coordinates instead of cartesian coordinates. The analytical solution of the equation, the wave function, is therefore proportional to a product of two infinite power series.[7]"
I suppose if you're okay with just having that product of two infinite power series (and numerically evaluating them on your computer) as your wavefunction, then you don't need to use approximate ansatze, but if you want something simpler, then the added complexity of using an anstaz would be valuable.
