Here is a case where it is non-Abelian. I use $K$ of class number 3. If I use the Gross curve, it is Abelian. If I twist in $Q(\sqrt{-15})$, it is Abelian for every one I tried, maybe because it is one class per genus. My comments are not from an expert.
> K<s>:=QuadraticField(-23);
> jinv:=jInvariant((1+Sqrt(RealField(200)!-23))/2);
> jrel:=PowerRelation(jinv,3 : Al:="LLL");
> Kj<j>:=ext<K|jrel>;
> E:=EllipticCurve([-3*j/(j-1728),-2*j/(j-1728)]);
> HasComplexMultiplication(E);
true -23
> c4, c6 := Explode(cInvariants(E)); // random twist with this j
> f:=Polynomial([-c6/864,-c4/48,0,1]);
> poly:=DivisionPolynomial(E,3); // Linear x Linear x Quadratic
> R:=Roots(poly);
> Kj2:=ext<Kj|Polynomial([-Evaluate(f,R[1][1]),0,1])>;
> KK:=ext<Kj2|Polynomial([-Evaluate(f,R[2][1]),0,1])>;
> assert #DivisionPoints(ChangeRing(E,KK)!0,3) eq 3^2; // all E[3] here
> f:=Factorization(ChangeRing(DefiningPolynomial(AbsoluteField(KK)),K))[1][1];
> GaloisGroup(f); /* not immediate to compute */
Permutation group acting on a set of cardinality 12
Order = 48 = 2^4 * 3
> IsAbelian($1);
false
This group has $A_4$ and $Z_2^4$ as normal subgroups, but I don't know it's name if any.
PS. 5-torsion is too long to compute most often.
One of the references at your oeis sequence is
D. A. Buell, Binary Quadratic Forms. Springer-Verlag, NY, 1989, pp. 224-241. On page 103, he proves that the narrow class group you ask about is isomorphic to the class group of binary quadratic forms. On page 82, Buell points out that computations show about 80 percent of positive prime discriminants have class number one. Note these primes are $1 \pmod 4.$ And it is certainly conjectured that the list is infinite.
INSERTED::: I found the quote I really wanted, Buell top of page 79:
Aside from the exceptional
discriminants of very small magnitude,
an odd class number can only occur for
an odd prime discriminant.
The Cohen-Lenstra heuristics are discussed in http://en.wikipedia.org/wiki/Class_number_problem including this quote
For real fields they predict that
about 75.446% will have class number
1, a result that agrees with
computations
P. S. Buell let some errors into the tables. He sent me corrections for Table 2A, odd positive fundamental discriminants.
EDIT: found it! Buell page 101, Theorem 6.19, part (b): If the discriminant $\Delta$ of the quadratic field $\mathbf Q(\sqrt \Delta)$ is positive, and a solution exists to the equation $x^2 - \Delta y^2 = -4,$ then the class group and the narrow class group are isomorphic.
Then , on page 162, Theorem 9.3, If $ \Delta = p,$ with $p$ a prime congruent to 1 modulo 4, then equation (9.1) is solvable, and the class number is odd.
Meanwhile the mentioned equation is
(9.1) $ X^2 - \Delta Y^2 = -4 $
So for positive prime discriminants $p \equiv 1 \pmod 4,$ everything agrees, and we are discussing the narrow class group in this case.
There is a simple proof that (9.1) is solvable in the case mentioned, it should be in Mordell's book on Diophantine Equations, the proof should actually be for $-1$ rather than $-4.$
EDIT TOOO: I thought I might put in the example I mentioned, discriminant 9973 (prime), I have programs that take a given indefinite binary form to "reduced" status and then show the entire cycle of reduced forms, so here is what happens for $x^2 + x y - 2493 y^2,$
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle
Input three coefficients a b c for indef f(x,y)= a x^2 + b x y + c y^2
1 1 -2493
0 form 1 99 -43 delta -2
1 form -43 73 27 delta 3
2 form 27 89 -19 delta -4
3 form -19 63 79 delta 1
4 form 79 95 -3 delta -32
5 form -3 97 47 delta 2
6 form 47 91 -9 delta -10
7 form -9 89 57 delta 1
8 form 57 25 -41 delta -1
9 form -41 57 41 delta 1
10 form 41 25 -57 delta -1
11 form -57 89 9 delta 10
12 form 9 91 -47 delta -2
13 form -47 97 3 delta 32
14 form 3 95 -79 delta -1
15 form -79 63 19 delta 4
16 form 19 89 -27 delta -3
17 form -27 73 43 delta 2
18 form 43 99 -1 delta -99
19 form -1 99 43 delta 2
20 form 43 73 -27 delta -3
21 form -27 89 19 delta 4
22 form 19 63 -79 delta -1
23 form -79 95 3 delta 32
24 form 3 97 -47 delta -2
25 form -47 91 9 delta 10
26 form 9 89 -57 delta -1
27 form -57 25 41 delta 1
28 form 41 57 -41 delta -1
29 form -41 25 57 delta 1
30 form 57 89 -9 delta -10
31 form -9 91 47 delta 2
32 form 47 97 -3 delta -32
33 form -3 95 79 delta 1
34 form 79 63 -19 delta -4
35 form -19 89 27 delta 3
36 form 27 73 -43 delta -2
37 form -43 99 1 delta 99
38 form 1 99 -43
minimum was 1rep 1 0 disc 9973 dSqrt 99.864908752 M_Ratio 9973
Automorph, written on right of Gram matrix:
1938468307 1282756488
-569466216 1395887771
Trace: -960611218 gcd(a21, a22 - a11, a12) : 8
=========================================
where one can see that halfway through the cycle (step numbered 18) we have represented $-1$ as a first or third coefficient of an equivalent form. This is Buell's description of "reduced" and adjacent forms. The middle coefficient is always positive and tends to be comparatively large.
THere is an extreme behavior when an odd prime is a square plus 4, as $293 = 17^2 + 4,$
=========================================
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle
Input three coefficients a b c for indef f(x,y)= a x^2 + b x y + c y^2
1 1 -73
0 form 1 17 -1 delta -17
1 form -1 17 1 delta 17
2 form 1 17 -1
minimum was 1rep 1 0 disc 293 dSqrt 17.117242769 M_Ratio 293
Automorph, written on right of Gram matrix:
-1 -17
-17 -290
Trace: -291 gcd(a21, a22 - a11, a12) : 17
=========================================
There may be only a finite number of odd primes $p$ with $p = u^2 + 4,$ unknown.
As GH mentioned discriminant $p = 4 u^2 + 1,$ here is one with class number one:
=========================================
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle
Input three coefficients a b c for indef f(x,y)= a x^2 + b x y + c y^2
1 1 -169
0 form 1 25 -13 delta -1
1 form -13 1 13 delta 1
2 form 13 25 -1 delta -25
3 form -1 25 13 delta 1
4 form 13 1 -13 delta -1
5 form -13 25 1 delta 25
6 form 1 25 -13
minimum was 1rep 1 0 disc 677 dSqrt 26.019223663 M_Ratio 677
Automorph, written on right of Gram matrix:
-53 -1352
-104 -2653
Trace: -2706 gcd(a21, a22 - a11, a12) : 104
=========================================
Finally, the first odd number where we are a little surprised at the impossibility of the negative Pell equation is 205, in that $205 = 5 \cdot 41$ is not divisible by any prime $q \equiv 3 \pmod 4$ but $x^2 - 205 y^2 = -1$ and $x^2 - 205 y^2 = -4$ are impossible,
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle
Input three coefficients a b c for indef f(x,y)= a x^2 + b x y + c y^2
1 1 -51
0 form 1 13 -9 delta -1
1 form -9 5 5 delta 1
2 form 5 5 -9 delta -1
3 form -9 13 1 delta 13
4 form 1 13 -9
minimum was 1rep 1 0 disc 205 dSqrt 14.317821063 M_Ratio 205
Automorph, written on right of Gram matrix:
2 27
3 41
Trace: 43 gcd(a21, a22 - a11, a12) : 3
=========================================
Here we are using a theorem of Lagrange, that given an indefinite form with positive nonsquare discriminant $\Delta,$ then any nonzero integer $n$ primitively represented by the form, also satisfying $ |n| < \frac{1}{2} \sqrt \Delta,$ occurs as the first coefficient of at least one form in the cycle of adjacent equivalent reduced forms. This is Theorem 85 on page 111 of Introduction to the Theory of Numbers by Leonard Eugene Dickson (1929). I do not see that Buell proves (or even quotes) the full Lagrange result, but his Theorem 3.18 on page 42 does the result for target $-4$ that we use here. An unusual, and quite pretty, description of the numbers represented by an indefinite binary is in pages 18-23 of The Sensual Quadratic Form by John Horton Conway. The presentation includes a method that quickly finds all represented numbers up to some bound in absolute value, this is the Climbing Lemma on page 11.
Oh, $x^2 + x y - k y^2$ represents a superset of the numbers represented by $x^2 - (4k+1) y^2.$
Best Answer
I would say there are two issues: first of all, the good thing making the above machinery work is that the moduli problem you look at, represented by $Y_\text{nonsplit}(n)$, works for all quadratic fields at once. In other words, you are interested in looking at all possible curves with CM by some imaginary quadratic field and you end up looking at points always in the same moduli space. This is not the case for Hilbert-Blumenthal surfaces with real multiplication by a real quadratic $F$, who do possess a moduli space (stack, indeed) $\mathfrak{M}_F$ which heavily depends upon $F$.
But suppose this is only a technical issue, either by constructing some monstrous $\mathfrak{M}$ parametrizing abelian surfaces with al multiplication by some quadratic field not specified in the moduli problem; or by hoping, in your "pipe dream" that one proves that there are infinitely many $F$ with at least one element in $\mathfrak{M}_F(\mathbb{Q})$. Then a major (and, to my knowledge and understanding, not a merely technical one) problem is the following: the crucial step in the CM procedure, as you said, is the implication (now $K$ is imaginary quadratic) $$ h(K)=1\Leftrightarrow j(E_K)\in\mathbb{Q} $$ (it is indeed integral, but that won't matter, here) where $E_K$ is some/any elliptic curve with CM by $K$. This is false for real multiplication: to convince yourself of the failure of one implication, observe that there are elliptic curves with real multiplication by $\mathbb{Q}$ (what else?) which are not defined over the rationals although $h(\mathbb{Q})=1$. Likewise, the field of definition of the abelian surface tells you almost nohing on the arithmetic of the field of real multiplication, unlike the CM case, preventing you from reducing the class number problem to counting varieties (with bonus structure) defined over $\mathbb{Q}$ - i.e. rational points on some moduli space.