Ygthkb

In Atkins Physical Chemistry 10th edition in chapter number 10B it is given that within the Born-Oppenheimer approximation, the Schrödinger equation can be solved for the molecular species $ce{H2+}$ but to simplify things and to obtain a method more readily generalizable to more complex species, we use the linear combination of atomic orbitals concept.

But my question is: how do we know that the wave function obtained from the linear combination obtained (say of two 1s orbitals of each atom) satisfies the Hamiltonian of this new system involving two nuclei? (The modified Hamiltonian is also given on the same page in the book).

Or do we need to make some iterative calculation as in the Hartree-Fock self consistent field method.

PS. I am a beginner so please try to keep the answer technically simple if possible.

As you requested a technically simple answer, I'll describe every step thoroughly to be safe. I'm also assuming we're dealing with the time-independent Schrödinger equation, and will refer to it below just as the Schrödinger equation (I'm happy to explain why in the comments if you're interested).

By solving the Schrödinger equation for your system, you find the eigenstates/eigenfunctions of your hamiltonian, as well as the energy values of each of those eigenstates. Because your state function (the specific wave function you're interested in -- the one that describes your system in its current state) obeys the Hamiltonian of the system, then you are guaranteed that the state function can be exactly written (with no loss of accuracy or approximation) as a linear combination of the eigenfunctions you found by solving the Schrödinger equation.

For any system containing two or more electrons, the electron-electron interaction terms render it impossible to separate the Schrödinger equation into separable variables, meaning you'll need to use further approximations to solve the problem. However, within the Born-Oppenheimer approximation, the Schrödinger equation can be solved exactly when there's a single electron. That being said, the dihydrogen cation is a great example of where that doesn't mean it can be done easily.

A linear combination of atomic orbitals is a workaround to having to solve the Schrödinger equation for your system. Using them is similar in mindset to using the eigenfunctions to compose your state function. However, this time, there is no guarantee that your wave function of interest will exactly be a mixture of the eigenfunctions. This is because here, your eigenfunctions are for the hydrogen atoms, which is not exactly the system we are looking at! As such, you often need to use numerous basis functions to get closer to agreement with the exact solution. The usage of this substitution will certainly come with an optimization through a self consistent field or other theoretical method, in order to find the best possible combination of the eigenfunctions to reflect your state function (you wouldn't expect a 2p and a 1s atomic basis function to equally contribute to the dihydrogen ground state wave function).

TLDR: The atomic orbitals are not eigenfunctions of the dihydrogen cation, and so substituting the actual state function with a linear combination of atomic orbitals will be approximate. However, if you use the complete basis set (CBS) limit and use a sufficient level of theory, you can recover excellent agreement with the exact solution.

I am paraphrasing from an article in Wikipedia nnd applying it to your specific case.

To simplify things and to obtain a method more readily generalizable to more complex species, we use the linear combination of atomic orbitals concept.

The special thing about the ground state solution of the Schrödinger equation is that it has the lowest energy. The variational method uses the idea that any other wave function will have a higher energy.

If you have a bunch of wave functions and want to know which one is closest to the real one for the ground state, you calculate the energy each one of these would have (using the Hamiltonian for $ce{H2+}$ in your case). Instead of taking unrelated functions, you can also build trial wave functions from linear combinations of basis functions, in your case just the 1s functions of the two hydrogen atoms, each multiplied by a coefficient. This is called a linear combination. Then, you optimize those two coefficients until you get the lowest energy.

How do we know that the wave function obtained from the linear combination obtained (say of two 1s orbitals of each atom) satisfies the Hamiltonian of this new system involving two nuclei?

They don't satisfy the Hamiltonian but they come pretty close. We know that by comparing the calculated energy to experimental data or to more sophisticated computational results.

In your case, you get a pretty good solution (with two positive coefficients) and a pretty awful solution (with coefficients of opposite signs), the latter of which is called the anti-bonding case (i.e. an excited state that does not hold the atoms together).

Do we need to make some iterative calculation as in the Hartree-Fock self consistent field method.

No, you don't need to iterate because there is only one electron. With only two coefficients, you will be able to find the global minimum in a straightforward manner.

If you want better solutions, you add more basis functions to your trial function (like 2s, 2p etc, or something completely different). The calculation will get more complicated, but the solutions will be closer to the actual ground state wave function. Now it will be more difficult to find the set of coefficients that give you the wave function with the lowest energy, and there might be some iteration to search for the global minimum.

There are some tricks to cut down on computational time (or to do some of this qualitatively in your head or on paper). One central idea is to sort basis functions by symmetry, and keep different symmetries separate from each other (certain terms will be zero when the two basis functions appearing in the are zero). Another trick is to pre-compute parts of some terms, and then just vary the coefficients.

In your case no tricks are necessary because you are just combining two basis functions.