Ndtyjky

For a system of two identical particles, where $r_1$ is the position vector of particle 1 and $r_2$ is the position vec. of particle 2, the wave function should be one of the plus or minus states:
begin{equation}
psi _ pm (r_1,r_2) = A [psi _a(r_1) psi_b (r_2) pm psi_b(r_1) psi_a(r_2)] ,
end{equation}
where $psi_a$ and $psi_b$ are the wave functions of particle 1 and 2 respectively [equation 5.10 of Griffith's Intro to Quantum Mechanics 2nd Ed.].

I see that this equation makes the wave function $psi_pm$ treat the two particles identically, but I don't know of any proof that it is actually the only way of writing this wave equation to treat them identically. For example, why not a wave function like:
begin{equation}
psi _ pm (r_1,r_2) = A sqrt{psi^2 _a(r_1) psi_b^2 (r_2) pm psi^2_b(r_1) psi^2_a(r_2)}~?
end{equation}

The requirement is
$$
psi(x_1,x_2) =
begin{cases}
psi(x_2,x_1) & text{for bosons} \
-psi(x_2,x_1) & text{for fermions}.
end{cases}
tag{1}
$$
This property is required by the spin-statistics theorem in relativistic quantum field theory. Since non-relativistic quantum mechanics is supposed to be an approximation to relativistic quantum field theory, we also enforce it in non-relativistic QM.
A special case of equation (1) is
$$
psi(x_1,x_2) approx
begin{cases}
f(x_1)g(x_2)+f(x_2)g(x_1) & text{for bosons} \
f(x_1)g(x_2)-f(x_2)g(x_1) & text{for fermions},
end{cases}
tag{2}
$$
but like Lewis Miller's answer said, this is only a special case. The general requirement is equation (1).

The square-root example written in the question does not satisfy the requirement (1).

This product form of two-particle wave function is only correct if the particles are not interacting. Nevertheless, it is often used as a first approximation, and if you take the expectation value of the true Hamiltonian and minimize it (by taking a variational derivative) you get the two-body Hartree-Fock equation which is often used to approximate the ground state energy and wave function for many-body systems of Fermions. This approximation is often called the mean field approximation.

You’re forgetting that the wave function must also satisfy the TISE. With this condition the combine wavefunctions must be the sum of a permutation of products.

If we have two particles, one in state $psi_a$ and the other in state $psi_b$, then the state vector would be $|psi_arangle|psi_brangle$.

However, if the particles are indistinguishable, then it is equally likely to have the opposite be true (i.e. the "first" particle in state $psi_b$ and the "second" particle in state $psi_a$). Therefore, we would want the entire state to be a linear combination of these two states, each with equal weight. Therefore we end up with

$$|Psirangle=|psi_arangle|psi_branglepm|psi_brangle|psi_arangle$$

If we choose to work in the $|r_1rangle|r_2rangle$ basis, then we end up with the expression you state.

I think the issue with your state is that it is not a "nice" linear combination of the states where one particle is in state $psi_a$ and the other in state $psi_b$. We need this if we want the postulate to hold that when $|psirangle=sum c_i|psi_nrangle$, we know that there is a probability of $|c_i|^2$ to measure the system to be in state $|psi_irangle$