I strongly conjecture that the maximum absolute determinant of a Latin square can be attained by a circulant matrix. For example,
$$\pmatrix {5&4&2&3&1 \\ 1&5&4&2&3 \\ 3&1&5&4&2 \\ 2&3&1&5&4 \\ 4&2&3&1&5}$$
has determinant $2325$, which is indeed the maximum absolute value of the determinant of a Latin square of size $5 \times 5$. The sign of the determinant is not important because changing two rows always gives a Latin square with positive determinant, if the given has a negative determinant.
Is this conjecture true ?
If the conjecture is true, there would be a suitable way to find the maximum absolute value of the determinant of a Latin square of size $n \times n$. I tried to find out how an arbitrary Latin square can be transformed in the circulant form without changing its determinant, but without success.
If the conjecture is true, the maximal values are (on the left side, the top row of the matrix is given, the sign is not considered):
$$ 1\ 2\ \ \ \ \ \ \ 3$$
$$ 3\ 1\ 2\ \ \ \ \ \ 18$$
$$ 4\ 1\ 2\ 3\ \ \ \ \ 160$$
$$ 5\ 4\ 2\ 3\ 1\ \ \ \ \ \ 2 \ 325$$
$$ 6\ 5\ 3\ 2\ 4\ 1\ \ \ \ \ 41 \ 895$$
The next $4$ values are $961\ 772,\ 26\ 978\ 400,\ 929\ 587\ 995\ and \ 36\ 843\ 728\ 625$. So, an upper bound for the determinant of a $9 \times 9$-Sudoku-matrix would be $929\ 587\ 995$.
Best Answer
Let $A=[C_1,\cdots,C_n]\in M_n$ be a latin square, $M=A^TA$, $a_n=1\times n+2\times (n-1)+\cdots+n\times 1$, $b_n=\sum_{i=1}^ni^2=n(n+1)(2n+1)/6$ and $u=[1,\cdots,1]^T$. Note that $m_{i,j}=<C_i,C_j>\in[a_n,b_n]$ and $m_{i,i}=b_n$.
Moreover, for every $i$, $\sum_{j}m_{i,j}=n^2(n+1)^2/4$ and $c_n=\sum_{j|j\not= i}m_{i,j}=(n-1)n^2(n+1)(3n+2)/12$, which leads to the mean of the off-diagonal entries being $\mu_n=n(n+1)(3n+2)/12$.
Put $M=b_n.I+S$ where the diagonal of the symmetric matrix $S$ is $0$, $s_{i,j}\in [a_n,b_n]$ and $g_i(S)=\sum_{j|j\not= i}s_{i,j}-c_n=0$; we seek the max of $f:S\rightarrow \det(M)$ under the constraints $g_i(S)=0$. According to Lagrange, there is $(\lambda_i)_i\in\mathbb{R}^n$ s.t., for every $i\not= j$, $tr(adj(M)E_{i,j})+\lambda_i=0$, that is $z_{j,i}+\lambda_i=0=z_{i,j}+\lambda_i$ where $z_{i,j}$ is the $i,j$ entry of $adj(M)$; then, for every $i$, the $(z_{i,j})_{j|j\not= i}$ are equal and, since $adj(M)$ is symmetric, the $(z_{i,j})_{i\not= j}$ are equal. Finally, $M^{-1}$ is in the form $D+\alpha U$ where $D=diag(d_i)$ is diagonal, $\alpha\in\mathbb{R}$ and $U=[u_{i,j}]$ is defined by $u_{i,i}=0$ and , if $i\not= j$, then $u_{i,j}=1$. Therefore $MM^{-1}=(b_nI+S)(D+\alpha U)=b_nD+SD+\alpha b_n U+\alpha SU=I$. We consider the diagonals of RHS and LHS: for every $i$, $b_nd_i+\alpha c_n=1$; consequently, the $(d_i)_i$ are equal and $M^{-1}$ is in the form $M^{-1}=\beta I+\alpha U$. Thus $S$ is a polynomial in $U$ and all the non-zero entries of $S$ are equal to $\mu_n$.
Conclusion. The (ideal) maximum of $|\det(A)|$ is reached when the $(m_{i,j})_{i\not= j}$ are $\mu_n$. Note that, in general, $\mu_n$ is not an integer and the considered ideal maximum is not reached; yet, the idea is to find a matrix $A$ s.t. $A^TA$ is close to $b_nI+\mu_n U$.
The case $n=5$. Then $\mu_5=42.5$ and $\sqrt{\det(b_5I+\mu_5 U)}\approx 2343.75$. The real maximum $2325$ is reached (according to Peter) with a circulant matrix s.t. the $(m_{i,j})_{i\not= j}$ are $42$ or $43$.
The case $n=6$. Then $\mu_6=70$ but $\sqrt{\det(b_6I+\mu_6 U)}\approx 42439.23$ is not an integer. A candidate for the real maximum is $41895$; it is reached with a circulant matrix s.t. the $(m_{i,j})_{i\not= j}$ are $69,70,71$.
The case $n=9$. Then $\mu_9=217.5$ and $\sqrt{\det(b_9I+\mu_9 U)}\approx 934173632.8$. A candidate for the real maximum is $929587995$; it is reached with a circulant matrix and a sudoku; the last one is s.t. the $(m_{i,j})_{i\not= j}$ are $216,217,218,219$.
Remark. The above reasoning works, in the same way and with the same conclusion, when we replace $M$ with $M'=AA^T$. Thus, for such an ideal solution, $A^TA=AA^T$, that is, $A$ is a normal matrix. Thus, good candidates for the Graal are the circulant matrices -that are normal-.
EDIT. We consider the Cholesky decomposition $b_nI+\mu_n.U=C^TC$; then $A^TA=b_n+\mu_n.U$ iff $A=PC$ where $P\in O(n)$. Finally, the idea is to rotate the known triangular matrix $C$ to make it close to an integer matrix (of course, it is not easy). For example, when $n=9$, a candidate solution $B$ has $[1,2,4,3,8,9,5,6,7]$ as first row; clearly, $B=RC$, where $R$ is close to an orthogonal matrix $T$ (consider the polar decomposition of $R$); if we replace $R$ with $T$, then we obtain $[1.02,2.07,4.08,2.89,7.99,9.02,4.95,5.92,7.07]$ as new first row; of course, this is this last form that must be discovered. We may also require that $A$ is normal, that is equivalent to the linear equation in the unknown $P$: $MP-PCC^T=0$. Unfortunately, the associated kernel is big; for example,when $n=9$, its dimension is $65$.