This is meant to encourage Kevin Costello and Roland Bacher to do more numerical simulations to bound a(3).
I was impressed that Kevin found a maximum as large as 9/50, and that with 5 ones. I then computed the four essentially different forms that the determinant (or its negative) could take coming from a 3x3 matrix with 5 ones. Except for one of them, I found I could simplify the expressions by using the substitutions z_i = a_i - 1. I think when these forms are investigated, an explanation for the non-1 entries will be forthcoming. I also wonder if perhaps the non-one entries can be replaced with 4/9, -1/9, -5/9, and -1, which if it works would be explained by the heuristic of spacing out the non-1 values. A similar analysis might be done using 4 ones instead of 5.
ADDED 2011.12.23 Thanks to Will Orrick for catching the erroneous "four" above and reminding me of a fifth form. While the substitutions might help in analyzing the now 120 possible expressions for a determinant containing 5 ones, I now have a couple of nice observations that might help estimate a(3) and possibly a(n). I have convinced myself that a(n) is decreasing as n > 1 increases, but only have a suggestion as to how to show it.
For a(3), one of the configurations involving 4 ones has them arranged in a 2x2 square, leading to a determinant of (a-b)(c-d) and a tentative upper bound for this case of 1/4. While it is possible that 1/4 could be beaten by multisets which have only 3 or 2 ones in them, I am confident that the other posters are close to determining a(3) and that it will be less than 1/4.
For general $n$, one can assume (as Will Orrick noted) 1 is in the multiset of $n^2$ values and sort them and then place the sorted set in the array in various ways while preserving relations like $a_{i,j} > a_{i,j+1}$. Now one can subtract multiples of one row from another to get rows of small norm (when consider as vectors in $R^n$). It may then be feasible to show that the given matrix is equivalent via such row operations to one in which the product of the norms of the rows is at most 1, and possibly smaller. An alternative picture is to say that the rows belong to the region of the positive unit orthant which satisfies $ x_1 \leq x_2 \leq \ldots \leq x_n$, and that the angle between a group of, oh, say, $n/2$ of them is less than $\epsilon(n)$ which would give a similar bound on the determinant. END ADDED 2011.12.23
Gerhard "Ask Me About System Design" Paseman, 2011.12.06