Are there orthogonal Cauchy-like matrices with rational entries for any given size? This is inspired by a recent question about the existence of orthogonal Cauchy-like matrices. It is proved that there are indeed such matrices, i.e. there are vectors  $x,y,r,s\in\mathbb R^n$ such that for the matrix $C$ defined by
$$ 
C_{ij} = \frac{r_i s_j}{ x_i - y_j},
$$ we have $C C^T = I$.
In the pdf with the solution, the following fact is used: If $C$ is a $n \times n$ real Cauchy matrix $C$, i.e. there are vectors  $x,y\in\mathbb R^n$ such that $$C_{ij} = \frac{1}{ x_i - y_j},$$ then its inverse admits the factorization
$$C^{-1} = D_aC^TD_b,$$
where $D_a,D_b$ denote the diagonal matrices with the entries of $a, b\in\mathbb R^n$, which can be computed by explicit formulae obtained by Lagrangian interpolation. BTW, it is not clear to me why the RHS features $C^T=(( \frac{1}{ x_j - y_i}))_{i,j=1}^n$, as the given formula $(7)$ in Theorem 1 of the original article (accessible here with different notation, as it calls $a,b$ what we call $x,y$ here) seems to contain rather $((  { x_j - y_i}))_{i,j=1}^n$ than the reciprocals.
In any case, my question is

whether the vectors $x,y,r,s\in\mathbb R^n$ defining a Cauchy-like matrix can all have rational entries, i.e. whether there is a Cauchy matrix $C$ such that each entry of the corresponding vectors $a,b$ in the factorization of $C^{-1}$ quoted above is a (rational) square?

Note that the given construction shows that there are Cauchy matrices where $a,b$ have only positive entries and uses their square roots at some point. But for rational square roots, I wouldn't know how to go about that even for $n=2$.
 A: Let us start with four rational numbers $x_1, x_2; y_1,y_2$ so that using the cross ratio
$$
r=(x_1,x_2,y_1,y_2)
$$
the numbers $-r$ and $1-r$ are rational squares.
(Exchanging the mid components $x_2, y_1$ bring $r$ into $1-r$. Exchanging the first and/or last two components we obtain the multiplicative inverse, etc. - so we want in the following to set in evidence squares times cross ratio values known to be squares.)

For instance for $0,1;-1,8/17$ we are producing $r=-9/16$. It may be simpler to follow the construction based on this example.
Let $C$ be the associated Cauchy matrix. In the example:
$$
C=\begin{bmatrix}1 & -17/8\\1/2 & 17/9\end{bmatrix}\ .
$$
Let $a_1,a_2;b_1,b_2$ be the rational squares:
$$
\begin{aligned}
a_1 &= (x_1-y_1)^2 \; (x_1, y_1, y_2, x_2) \ ,\\
a_2 &= (x_1-y_2)^2 \; (x_1, y_2, y_1, x_2) \ ,\\[2mm]
b_1 &= 1\ ,\\
b_2 &= (-1)\; \left(\frac{x_2-y_1}{x_1-y_1}\right)^2 \;(x_1, x_2, y_1, y_2) \ .
\end{aligned}
$$
Then we have $C^{-1}=D_a\; C^T\; D_b$.
Let now $L$ be a the Cartan-like matrix of the shape $L=D_r\; C\; D_s$. Then:
$$
\begin{aligned}
L &=D_r\; C\; D_s\ ,\\
L^T &=D_s\; C^T\; D_r\ ,\\
L^{-1} 
&=D_s^{-1}\; C^{-1}\; D_r^{-1}\\
&=D_s^{-1}\; D_a C^T D_b\; D_r^{-1}
\ ,\\[3mm]
&\qquad\text{ and we want $L^T=L^{-1}$, i.e.}\\[3mm]
D_s\; C^T\; D_r &= D_s^{-1}\; D_a C^T D_b\; D_r^{-1}\text{ i.e.}\\
C^T &= D_s^{-2}\; D_a C^T D_b\; D_r^{-2}\ .
\end{aligned}
$$
Recall that $-r$ and $1-r$ both squares
implies $a,b$ squares, as wanted in the OP, so we can arrange for $s,t$ with rational entries.

To make the above easy to test, here is some sage code making the computations.
var('x1,x2,x3,y1,y2,y3');

def r(s,t,u,v):
    return (s-u)/(s-v)/(t-u)*(t-v)

a1 = (x1-y1)^2 * r(x1, y1, y2, x2)
a2 = (x1-y2)^2 * r(x1, y2, y1, x2)
b1 = 1
b2 = (-1) * (x2-y1)^2 / (x1-y1)^2 * r(x1, x2, y1, y2)

C = matrix([ [1/(x1-y1), 1/(x1-y2)] , [1/(x2-y1), 1/(x2-y2)] ])
Da = diagonal_matrix( [a1, a2] )
Db = diagonal_matrix( [b1, b2] )

print("Is C^-1 = Da * C^T * Db? %s"
      % bool(C^-1 == Da * C.transpose() *Db))
 

And we get:
Is C^-1 = Da * C^T * Db? True


We use now instead of general variables the special values:
x1, x2, y1, y2 = 0, 1, -1, 8/17

(just replace the first var line with the above, keep the next lines of the used code) and ask for the values of $a$, $b$:
sage: a1, a2, b1, b2
(16/25, 576/7225, 1, 9/4)

Now consider the matrix $L$ given by
sage: L = diagonal_matrix([sqrt(b1), sqrt(b2)]) * C * diagonal_matrix([sqrt(a1), sqrt(a2)])
sage: L
[ 4/5 -3/5]
[ 3/5  4/5]

Which is an orthogonal matrix. (I found this problem while searching for the tag elliptic-curves, but the above solution is maybe closer to K-theory.)
