Here is a bit of a hacky answer to your question in the affirmative
You observe that $k = \mathrm{gcd}(a,b)^2$. After plodding through various resources, it is because one can Viete Jump:
$$ (a,b) \mapsto \big(a_1, b_1\big) \mapsto \dots \mapsto (k,0) $$
and the gcd is conserved. Once I agree with you, let $a = \text{gcd}(a,b) \, c$ and $b = \text{gcd}(a,b) \, d$ so that
\begin{eqnarray} ab+1 &=& \frac{a^2+b^2}{\mathrm{gcd}(a,b)^2}= c^2 + d^2 \\
&=& \,\text{gcd}(a,b)^2 \, cd + 1
\end{eqnarray}
In this way there are two conditions $c, d$ might solve (where the two $\square$ are different) and $\text{gcd}(c,d)=1$:
\begin{eqnarray*}
c^2 - \square \, cd + d^2 &=& 1 \\
c^2 + d^2 &=& \square
\end{eqnarray*}
Hopefully these two equations leads to a contradiction.
Added 11/15 The answer is definitely no. Let $k = \mathrm{gcd}(a,b)$ We are trying to solve in integers:
\begin{eqnarray} c^2 - k\; cd + d^2 &=& 1 \\
c^2 + d^2 = \square
\end{eqnarray}
As I learned, the first can be solve with $(c,d) = (1,0)$ or $(k,1)$ and there are an infinite family of solutions using consecutive terms of a recursive sequence [2, 3]
$$ x_{n+1} = k \, x_n + x_{n-1} $$
There are sometimes ways to link Pythagorean triples to Pell equations [1] (Modular Tree of Pythagoras)
$$ x_{n+1}^2 + \frac{1}{2} x_n^2 < \sqrt{x_{n+1}^2 + x_n^2 } = x_{n+1}^2 \sqrt{1 + (x_n/x_{n+1})^2} < x_{n+1}^2 + \frac{1}{2} x_n^2 + 1$$
This cannot be an integer. So any time we solve the Pell equation, we cannot also solve Pythagoras. $\quad\quad\square$
Old Answer
This is discussed on Wikipedia's article on Vieta Jumping:
Nobody of the six members of the Australian problem committee could solve it. Two of the members were George Szekeres and his wife, both famous problem solvers and problem creators. [...] The problem committee submitted it to the jury of the XXIX IMO marked with a double asterisk, which meant a superhard problem, possibly too hard to pose. After a long discussion, the jury finally had the courage to choose it as the last problem of the competition. Eleven students gave perfect solutions.
Among the eleven was Bau Chau NgĂ´ (Fields Medal 2010). His work on the Fundamental Lemma also has a jumpy flavor [1, 2, 3] but is quite advanced.
The discussion on YouTube is helpful as well. These videos give a thorough discussion of different ways to solve
These may not directly solve your problem but provide historical context and indicate possible strategies.
In the Wikipedia article, the example of Viete Jumping is IMO 1988/6 -- the same as asked in the question:
Let $a,b$ be positive integers such that $ab+1$ divides $a^2 + b^2$ show that $ \frac{a^2 + b^2}{ab+1}$ is a perfect square.
and the solution goes in three steps
#1 Let $a, b \geq 0$ be solutions to $\frac{a^2 + b^2}{ab+1} = k$ such that $k$ is not a perfect square: $k \neq \square$
#2 Starting from $(a,b)$ we can try to generate another solution $(x,b)$ which solves the quadratic equation:
$$ x^2 - kb\, x +(b^2 - k) = 0 $$
The map $(a,b) \mapsto (a_1,b)$ is our Vieta jumping Since both $a, a_1$ are acceptable solutions we have:
$$ (x-a)(x-a_1) = x^2 - (a + a_1) x + aa_1 = 0$$
By the Viete equations (comparing the coefficients. We find out two things:
#3 If $a \geq b$ we can deduce that $a_1 \geq 0$ (is positive) and additionally $ a > a_1 \ge 0 $
Summary We've show that given two positive numbers $a,b$ solving $\frac{a^2+b^2}{ab+1}=k$ with $k \neq \square$ we can always find another solution $(a_1,b)$ solving the same equation with $a > a_1$.
Then the Viete jump consists of a map:
$$ (a,b) \mapsto \left\{\begin{array}{rc}
(b,a) & \text{ if } a \leq b \\
(\frac{b^2-k}{a},b) & \text{ if } a \geq b
\end{array}\right.$$
While this does not solve your problem -- to show that $ab+1 \neq \square$ -- it does indicate possibly where to start and some possible resources.
A quick use of Bezout formula shows that $ab+1$ should also divide
$$ \big[(a^2 + b^2) + 2(ab+1)\big] +
\big[(a^2 + b^2) - 2(ab+1)\big] = (a+b)^2 +(a-b)^2 $$
and this could lead to your contradiction.
Best Answer
I'll summarise progress here to prevent the comments becoming hard to read. I have made this community wiki, so feel free to add any progress made.
Parametric solutions found (by Ivan Neretin): $$\begin{align}m = k(2k+1), &\quad n=k(2k-1) \\ m = k(k-1)(k^2+k-1), &\quad n=k(k-1)(k+1)\end{align}$$ Solutions found outside of the above patterns (covers all solutions for $n \leq m \leq 2300000$):
Necessary conditions (by Geoff Robinson): Either $m$ and $n$ share a common factor, or they are both odd. If they are both odd, then $m+n \equiv 0$ mod $4$.
Moreover, $m$ and $n$ have the same parity. Proof: The denominator $B=m^2+n^2+m+n$ is even. The numerator $T=m^3+n^3$ is a multiple of $B$, so is even, so $m^3$ and $n^3$ have the same parity, so $m$ and $n$ have the same parity.
This means that if $d$ is odd, $a$ and $b$ are odd. I (Rosie F) notice that in most cases $d$ is even. But it is quite rare that either $a$ or $b$ is even.
We also can remark that the only coprime solution is $(m,n)=(3,1)$. To show this assume $(m,n)=1$. Also note that $m^3+n^3=(m+n)(m^2+n^2-mn)$, so the original quantity is an integer iff $$m^2+n^2+m+n\mid (m+n)(m+n+mn)$$ Of course this is true iff $m^2+n^2+m+n\mid \gcd(m^2+n^2+m+n, m+n)(m+n+mn)$, which simplifies to $m^2+n^2+m+n\mid \gcd(m+n, 2mn)(m+n+mn)$. Much like Geoff noted, this can only be the case when $\gcd(m+n, 2mn)>1$ since $m^2+n^2+m+n>m+n+mn$. However under the assumption that $\gcd(m,n)=1$, this is only the case when $m$ and $n$ are both odd and $\gcd(m+n, 2mn)=2$. Hence our original quantity is an integer only if $$m^2+n^2+m+n\mid 2(m+n+nm)$$ and $m$ and $n$ are odd. Assume $m^2+n^2+m+n\neq 2(m+n+nm)$, then for some prime $p$, $m^2+n^2+m+n\mid 2(m+n+nm)/p$ and thus $$m^2+n^2+m+n\le 2(m+n+nm)/p\le m+n+nm$$ which is false. Hence $m^2+n^2+m+n=2(m+n+nm)$ and $(m-n)^2=n+m$. This gives $m-n\mid n+m$ so $$m-n=\gcd(m-n,n+m)=\gcd(m+n, 2m)=2$$ Thus $m=n+2$. Subbing this into $(m-n)^2=m+n$ gives $(m,n)=(3,1)$. Thus if $\gcd(m,n)=1$ we must have $(m,n)=(3,1)$, which if we test works.