Small inter nuclear separation limit for Diatomic molecule Let’s take the a simple $H_2^+$ molecule, where there is only electron which is $r_a$ away from the first proton and $r_b$ away from the other one.
Let’s call the separation between the two protons $R$.
As $R\rightarrow \infty$, the electron will stick to one of the two protons, so the wavefunction will be: $$ \phi =  N_{\pm}(1s_a \pm 1s_b).$$
I recognise the two solutions as the gerade and ungerade orbitals. The $1s$ means ground state around each proton.
I can work out the normalisation constant to be $$ N_{\pm} = \sqrt{\frac{1}{2(1\pm S)}}, \, S = \int 1s_a 1s_b \mathrm{d}^3 r $$
Now in the limit of $R \rightarrow 0$, the gerade solution becomes just $1s$ which makes sense, but the ungerade is not defined - what happens to it?
 A: Think physically about the shape of the odd state.  It has a nodal plane at the midpoint of the line connecting the two nuclei; this is the only node to the wave function.  As the nuclei become coincident, that nodal plane ends up passing through the combined nucleus; since this is the only node, the resulting wave function is the $2p_{z}$ state, if the axis of the molecule was originally oriented along the $z$-axis.
You will not get this by taking a limit of $1s$ wave functions.  When the nuclei are very far apart, the overlap integral that determines the energy difference between the odd and even states is very small, and the wave function in the vicinity of each nucleus looks very close to a $1s$.  However, as you the distance between the nuclei gets comparable to the Bohr radius $a_{0}$ or smaller, the true wave function will need to be constructed out more than just $1s$ basis states.  When the separation is small comparable to $a_{0}$, you get a very complicated wave function, although it simplifies to just the $2p$ state once the separation gets very small compared to $a_{0}$.
A: Let $\Psi(\vec{r})$ be $1s_a$ wave function and $\Psi(\vec{R}+\vec{r})$ be $1s_b$ wave function. As $\vec{R} \to 0$ we have for ungerade state:
$$
\phi_{-}(\vec{r}) = N_{-}(\Psi(\vec{R}+\vec{r})-\Psi(\vec{r})) \approx N_{-} \nabla\Psi(\vec{r})\vec{R}
$$
Thus the limit (not normalized) is $\nabla\Psi(\vec{r})\vec{n}$, where $\vec{n} = \vec{R}/|\vec{R}|$. This function has nodal plane as Buzz said. Also it has discontinuity at $\vec{r} = 0$. But I think this function does not coincide with $2p$ state.
