Return to Answer

Here is a partial answer, which I probably will edit later to fill in some steps. Let me know if there is something unclear. I will first try to resemble the core ideas of DFT in a few paragraphs, which hopefully also answers some of your questions.

1. The Hohenberg-Kohn Theorems:

The external potential operator for your system of $N$ (identical) fermions (e.g. electrons) is of the form $V=\sum\limits_{i=1}^N v(x_i)$. The function $v:\mathbb R^3\longrightarrow \mathbb R$ is the external potential, which every electron feels (e.g. the presence of the fixed ions of a lattice in some solid material). The effective potential also takes this form $V_\mathrm{eff} =\sum\limits_{i=1}^N v_{\mathrm {eff}}(x_i)$.

So to get the ground state density of interest from your non-interacting system, you have to "choose" $v_{\mathrm ext}$ appropriately. Note that $v$ and $v_\mathrm{eff}$ are functions of a single variable (or three if you want), irrespective of the number of particles $N$ you consider.

If you knew the effective potential, you could build the Hamiltonian of the corresponding non-interacting system, which by definition is $$H_{\mathrm{eff}}=T+V_{\mathrm{eff}}= \sum\limits_{j=1}^N h_i \quad ,$$ where $h_i =- \frac{\hbar^2}{2m}\nabla_i^2 + v_{\mathrm {eff}}(x_i)$. Without going into too much detail, we can show (with some assumptions, e.g. non-degeneracy) that the ground state is of the form of a Slater-determinant: If $$h\varphi_j=\epsilon_j \varphi_j \quad ,$$ then an eigenstate of $H_{\mathrm{eff}}$ is given by $\psi\sim \varphi_{j_1} \wedge \varphi_{j_2} \wedge \ldots \varphi_{j_N}$ with energy $E=\sum\limits_{k=1}^N \epsilon_{j_k}$. This also shows you that the ground state of this system is the one such that its corresponding energy is build from the $N$ lowest $\epsilon_j$'s.

Consequently, in order to solve the non-interacting Schrödinger equation with the Hamiltonian $H_\mathrm{eff}$, you can solve the single-particle Schrödinger equation(s) from which you take the solutions $\varphi_j$ with the $N$ lowest energies $\epsilon_j$ to build your ground state Slater determinant, from which you get the ground state density.

  But $v_\mathrm{eff}$ depends on the ground state density $n_0$, which we have to determine. This makes the single-particle Schrödinger equation(s) a non-linear problem and hence we have to solve the Kohn-Sham equations self-consistently.

1. The Hohenberg-Kohn Theorems:

The external potential operator for your system of $N$ (identical) fermions (e.g. electrons) is of the form $V=\sum\limits_{i=1}^N v(x_i)$. The function $v:\mathbb R^3\longrightarrow \mathbb R$ is the external potential, which every electron feels (e.g. the presence of the fixed ions of a lattice in some solid material). The effective potential also takes this form $V_\mathrm{eff} =\sum\limits_{i=1}^N v_{\mathrm {eff}}(x_i)$. The Hamiltonian of the non-interacting system is by definition $$H_{\mathrm{eff}}=T+V_{\mathrm{eff}}= \sum\limits_{i=1}^N h_i \quad ,$$ where $h_i =- \frac{\hbar^2}{2m}\nabla_i^2 + v_{\mathrm {eff}}(x_i)$. Without going into too much detail, we can show (with some assumptions, e.g. non-degeneracy) that the ground state is of the form of a Slater-determinant: If $$h\varphi_j=\epsilon_j \varphi_j \quad ,$$ then an eigenstate of $H_{\mathrm{eff}}$ is given by $\psi\sim \varphi_{j_1} \wedge \varphi_{j_2} \wedge \ldots \wedge \varphi_{j_N}$ with energy $E=\sum\limits_{k=1}^N \epsilon_{j_k}$. This also shows you that the ground state of this system is the one such that its corresponding energy is build from the $N$ lowest $\epsilon_j$'s.

Consequently, in order to solve the non-interacting Schrödinger equation with the Hamiltonian $H_\mathrm{eff}$, you can solve the single-particle Schrödinger equation from which you take the solutions $\varphi_j$ with the $N$ lowest energies $\epsilon_j$ to build your ground state Slater determinant, from which in turn you get the ground state density. But $v_\mathrm{eff}$ depends on the ground state density $n_0$, which we have to determine. This makes the corresponding single-particle Schrödinger equation a non-linear problem and hence we have to solve the Kohn-Sham equations self-consistently.

However, the conditional statement is important. We know that there exists some $V_\mathrm{eff}$ with the desired properties, but we don't know a priori how it looks like. There are ways however, explained in basically any reference I gave above, to construct this potential. What you will find is that $V_\mathrm{eff}$ depends on the functional derivative of the energy functional.

Answer to questions in the comments:

The external potential operator for your system of $N$ (identical) fermions (e.g. electrons) is of the form $\hat V=\sum\limits_{i=1}^N v(x_i)$. The function $v:\mathbb R^3\longrightarrow \mathbb R$ is the external potential, which every electron feels (e.g. the presence of the fixed ions of a lattice in some solid material). The effective potential also takes this form $V_\mathrm{eff} =\sum\limits_{i=1}^N v_{\mathrm {eff}}(x_i)$.

The Hamiltonian of the non-interacting system by definition is $H_{\mathrm{eff}}=T+V_{\mathrm{eff}}= \sum\limits_{j=1}^N h_i$, with $h_i =- \frac{\hbar^2}{2m}\nabla_i^2 + v_{\mathrm {eff}}(x_i)$. Without going into too much detail, we can show (with some assumptions, e.g. non-degeneracy) that the ground state is of the form of a Slater-determinant: If $h\varphi_j=\epsilon_j \varphi_j$, then an eigenstate of $H_{\mathrm{eff}}$ is given by $\psi\sim \varphi_{j_1} \wedge \varphi_{j_2} \wedge \ldots \varphi_{j_N}$ with energy $E=\sum\limits_{k=1}^N \epsilon_{j_k}$. This also shows you that the ground state of this system is the one such that its corresponding energy is build from the $N$ lowest $\epsilon_j$'s.

Summarizing: To solve the non-interacting Schrödinger equation with the Hamiltonian $H_\mathrm{eff}$, you can solve the single-particle Schrödinger equation $h\varphi=\epsilon \varphi$, from which you take the solutions $\varphi_j$ with the $N$ lowest energies $\epsilon_j$ to build your Slater determinant/the ground state/ground state density.

But $v_\mathrm{eff}$ depends on the ground state density $n_0$, which we have to determine. Thus, we have to solve the Kohn-Sham equations self-consistently.

However, the conditional statement is important. We know that there exists some $V_\mathrm{eff}$ with the desired properties, but we don't know a priori how it looks like. There are ways however, explained in basically any reference I gave above, to construct this potential. What you will find is that $V_\mathrm{eff}$ depends on the functional derivative of the energy functional and in particular it depends on the ground state density.

Answer to questions in the comments:

The external potential operator for your system of $N$ (identical) fermions (e.g. electrons) is of the form $V=\sum\limits_{i=1}^N v(x_i)$. The function $v:\mathbb R^3\longrightarrow \mathbb R$ is the external potential, which every electron feels (e.g. the presence of the fixed ions of a lattice in some solid material). The effective potential also takes this form $V_\mathrm{eff} =\sum\limits_{i=1}^N v_{\mathrm {eff}}(x_i)$.

If you knew the effective potential, you could build the Hamiltonian of the corresponding non-interacting system, which by definition is $$H_{\mathrm{eff}}=T+V_{\mathrm{eff}}= \sum\limits_{j=1}^N h_i \quad ,$$ where $h_i =- \frac{\hbar^2}{2m}\nabla_i^2 + v_{\mathrm {eff}}(x_i)$. Without going into too much detail, we can show (with some assumptions, e.g. non-degeneracy) that the ground state is of the form of a Slater-determinant: If $$h\varphi_j=\epsilon_j \varphi_j \quad ,$$ then an eigenstate of $H_{\mathrm{eff}}$ is given by $\psi\sim \varphi_{j_1} \wedge \varphi_{j_2} \wedge \ldots \varphi_{j_N}$ with energy $E=\sum\limits_{k=1}^N \epsilon_{j_k}$. This also shows you that the ground state of this system is the one such that its corresponding energy is build from the $N$ lowest $\epsilon_j$'s.

Consequently, in order to solve the non-interacting Schrödinger equation with the Hamiltonian $H_\mathrm{eff}$, you can solve the single-particle Schrödinger equation(s) from which you take the solutions $\varphi_j$ with the $N$ lowest energies $\epsilon_j$ to build your ground state Slater determinant, from which you get the ground state density.

But $v_\mathrm{eff}$ depends on the ground state density $n_0$, which we have to determine. This makes the single-particle Schrödinger equation(s) a non-linear problem and hence we have to solve the Kohn-Sham equations self-consistently.

This then leads to the known fact that you have to solve the Kohn-Sham equations self-consistently. In other words, given some approximation for the energy functional, you may start with an initial guess for the ground state density $\tilde n_0$, solve the Kohn-Sham equations (basically the non-interacting Schrödinger equation rewritten) and get a new ground state density. Repeat the cycle until convergence (which might be a topic on its own).

So to get the ground state density of interest from your non-interacting system, you have to "choose" $v_{\mathrm ext}$ appropriately. Note that $v$ and $v_\mathrm{eff}$ are functions of a single variable (or three if you want), irrespective of the number of particles $N$ you consider.

The Kohn-Sham equations are nothing but the non-interacting Schrödinger equation in disguise:

Summarizing: To solve the non-interacting Schrödinger equation with the Hamiltonian $H_\mathrm{eff}$, you can solve the single-particle Schrödinger equation $h\varphi=\epsilon \varphi$, from which you take the solutions $\varphi_j$ with the $N$ lowest energies $\epsilon_j$ to build your Slater determinant/the ground state/ground state density.

This then leads to the known fact that you have to solve the Kohn-Sham equations self-consistently. In other words, given some approximation for the energy functional, you may start with an initial guess for the ground state density $\tilde n_0$, solve the Kohn-Sham equations and get a new ground state density. Repeat the cycle until convergence (which might be a topic on its own).

So to get the ground state density of interest from your non-interacting system, you have to "choose" $v_{\mathrm ext}$ appropriately. Note that $v$ and $v_\mathrm{eff}$ are functions of a single variable (or three if you want), irrespective of the number of particles $N$ you consider.

Summarizing: To solve the non-interacting Schrödinger equation with the Hamiltonian $H_\mathrm{eff}$, you can solve the single-particle Schrödinger equation $h\varphi=\epsilon \varphi$, from which you take the solutions $\varphi_j$ with the $N$ lowest energies $\epsilon_j$ to build your Slater determinant/the ground state/ground state density.

But $v_\mathrm{eff}$ depends on the ground state density $n_0$, which we have to determine. Thus, we have to solve the Kohn-Sham equations self-consistently.