Broadly speaking, the information you gave in your question is up-to-date: ever since Gauss, number theorists have believed there should be infinitely many real quadratic fields of class number one, but we are not any closer to being able to prove this than Gauss was, to the best of my knowledge. Moreover you have identified the key problem: the fact that a real quadratic field has an infinite unit group means that there is an additional quantity in the analytic class number formula, the regulator, which is not present in the imaginary quadratic case, and it is hard to separate out the respective contributions of the class number and regulator.
In recent years a trend of research has been towards making increasingly precise quantitative conjectures about the expected behavior of the class number -- or class group -- of a real quadratic field of discriminant $D$. A lot of this comes under the heading of Cohen-Lenstra heuristics.
Since you asked for some pointers to recent work, here are two interesting papers that I found (from a google search, I fully admit):
I. Stephens and Williams, Computation of real quadratic fields with class number one.
II. Ono, Indivisibility of class numbers of real quadratic fields. (Wayback Machine), doi:10.1023/A:1001533613223
There is certainly plenty more literature available. A MathSciNet search for "class number" AND "real quadratic field" calls up precisely [well, not for long!] 500 papers...
Hint:
Have you learnt about Pell's equation?
Try adding a multiple of $y^2$ (on both sides) to complete the square on the left.
Best Answer
EDIT: there have been many comments on my answers asking about the use of the word automorph. This is a real thing! I did not just make up a word. For a fixed quadratic form, you get a group of integral automorphs. In just two variables, there is a good recipe for finding all. In three or more variables, it is a mess. Sometimes this is called the orthogonal group of the form, or the isometry group. If you think about finding all real solutions of the basic matrix equation, $A^T F A = F,$ where $F$ is a symmetric matrix associated with a quadratic form, the group part may be clearer, especially when $F=I.$ If $F$ has all rational entries, it is reasonable to solve this with rational or $p$-adic entries in $A.$ Finally, when $F,$ or at least $2 F,$ has integer entries, it is reasonable to ask about $A$ with integer entries.
I would like people to know more about this. In dimension 2 this has considerable overlap with algebraic number theory for quadratic fields. In dimension 3 or more, the theory of quadratic forms, in this case indefinite forms, separates from algebraic number theory to a considerable extent.
ORIGINAL:Given nontrivial $\tau^2 - n \sigma^2 = 1,$ we get
$$ \left( \begin{array}{cc} \tau & \sigma \\ n \sigma & \tau \end{array} \right) \left( \begin{array}{cc} 1 & 0 \\ 0 & -n \end{array} \right) \left( \begin{array}{cc} \tau & n \sigma \\ \sigma & \tau \end{array} \right) = \left( \begin{array}{cc} 1 & 0 \\ 0 & -n \end{array} \right). $$
As a result, if $x^2 - n y^2 = k,$ then we get the same $k$ for $$ \left( \begin{array}{cc} \tau & n \sigma \\ \sigma & \tau \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} \tau x + n \sigma y \\ \sigma x + \tau y \end{array} \right). $$
This 2 by 2 matrix is called an automorph of the quadratic form.
Every indefinite form $f(x,y) = a x^2 + b x y + c y^2$ where $\Delta = b^2 - 4 a c$ is positive but not a square, has such an automorph, leading to infinitely many solutions. Indeed, given $\tau^2 - \Delta \sigma^2 = 4,$ we get
$$ \left( \begin{array}{cc} \frac{\tau - b \sigma}{2} & a \sigma \\ -c \sigma & \frac{\tau + b \sigma}{2} \end{array} \right) \left( \begin{array}{cc} a & \frac{b}{2} \\ \frac{b}{2} & c \end{array} \right) \left( \begin{array}{cc} \frac{\tau - b \sigma}{2} & -c \sigma \\ a \sigma & \frac{\tau + b \sigma}{2} \end{array} \right) = \left( \begin{array}{cc} a & \frac{b}{2} \\ \frac{b}{2} & c \end{array} \right). $$ Therefore, if we have $ a x^2 + b x y + c y^2 = k,$ we have another with $$ \left( \begin{array}{cc} \frac{\tau - b \sigma}{2} & -c \sigma \\ a \sigma & \frac{\tau + b \sigma}{2} \end{array} \right) \left( \begin{array}{c} x \\ y \end{array} \right) = \left( \begin{array}{c} \frac{\tau - b \sigma}{2} x - c \sigma y \\ a \sigma x + \frac{\tau + b \sigma}{2} y \end{array} \right). $$
For previous answers in which I show how to use an automorph, see Solve the Diophantine equation $ 3x^2 - 2y^2 =1 $
How to find solutions of $x^2-3y^2=-2$?
Books: H.E.Rose, A Course in Number Theory, chapter 9, section 3, especially pages 162-164 in the first edition.
Thomas W. Cusick and Mary E. Flahive, The Markoff and Lagrange Spectra appendix 3 on pages 91-92.
Duncan A. Buell, Binary Quadratic Forms chapter 3, section 2, pages 31-34.
William J. LeVeque, Topics in Number Theory, volume 2, pages 24-29. The two volumes are available as a one volume paperback, LeVeque Book
Leonard Eugene Dickson, Introduction to the Theory of Numbers, especially pages 111-112.