square that can be written as mean of as many pairs of squares as possible If $b^2=\displaystyle\frac{a_1^2+c_1^2}2$ and $b^2=\displaystyle\frac{a_2^2+c_2^2}2$... and $b^2=\displaystyle\frac{a_n^2+c_n^2}2$, what is the largest possible value of $n$, or can $n$ be arbitrarily large?
e.g for $n=7$, $$325^2=\frac{49^2+457^2}2=\frac{65^2+455^2}2=\frac{115^2+445^2}2=\frac{175^2+425^2}2$$
$$=\frac{221^2+403^2}2=\frac{235^2+395^2}2=\frac{287^2+359^2}2,$$
and
$$425^2=\frac{7^2+601^2}2=\frac{85^2+595^2}2=\frac{175^2+575^2}2=\frac{205^2+565^2}2$$
$$=\frac{289^2+527^2}2=\frac{329^2+503^2}2=\frac{355^2+485^2}2.$$
 A: You are asking for an integer ($2b^2)$ that can be written as a sum of two squares in as many ways as possible, one of which is a sum of two equal squares. By the result that I cited in this answer, there is no bound to this number of ways. Concretely, for instance $2\times25^n$ can be written as a sum of two squares in $4(2n+1)$ ways; these include $4$ ways of the form $(\pm5^n)^2+(\pm5^n)^2$ that you don't want to count and the remaining $8n$ possibilities come in symmetry classes of $8$ (by signs and order of terms), for a total of $n$ classes.
In terms of your question, the square $(5^n)^2$ can be written as the average of two squares in $n$ non-equivalent ways. To find these expressions, take the Gaussian integers $1+\mathbf i$ and $2n$ copies of  $2+\mathbf i$, conjugate $i$ of the latter for $0\leq i<n$ and multiply everything together; the resulting $n$ Gaussian integers are all of norm-squared equal to $2\times25^n$, and their real and imaginary parts provide $n$ non-equivalent pairs of numbers, the averages of whose squares is $25^n=(5^n)^2$.
This is not the most economic way to get lots of expressions; it would be better to combine distinct prime numbers congruent to $1$ modulo $4$, rather than to take powers of one of them namely $5$. This explains your examples $325=5^2\times13$ and $425=5^2\times17$. 
Added: The general formula for the number of solutions for writing $N^2$ as the average of a set of two squares of distinct positive numbers, in terms of the prime factorization of $N$, is as follows: only the primes congruent to $1$ modulo $4$ contribute; multiply together for every nonzero multiplicity $m$ of such a prime the numbers $2m+1$, subtract $1$ from the the (odd) product so obtained, and divide by $2$ (which accounts for the ignored order). So for  $N=325=5^2\times13$ and $N=425=5^2\times17$ one gets $\frac{5\times3-1}2=7$ solutions, as indicated. Another value with many solutions is $N=5\times13\times17=1105$, namely $\frac{3\times3\times3-1}2=13$ solutions.
If one wants to disqualify solutions with a nontrivial common factor, as Jyrki Lahtonen suggests (I wouldn't know why), then each appropriate prime only contributes a factor $2$ independently of its multiplicity rather than $2m+1$ (but subtracting $1$ is omitted). This is beacuse mixing a Gaussian integer in a product with its complex conjugate introduces an integer factor, which will be common to the real and imaginary parts. In this variant one retains only $\frac{2\times2}2=2$ solutions for $N\in\{325,425\}$, and only $\frac{2\times2\times2}2=4$ solutions for $N=1105$ (namely $(73,1561)$, $(367,1519)$, $(809,1337)$, $(1057,1151)$). Even with this restriction the number will still be unbounded as $N$ acquires more and more useful prime factors.
A: From a beautiful theorem of Jacobi, the number of ways of writing an integer as a sum of two squares is equal to four times its number of divisors of the form $4n+1$ minus four times its number of divisors of the form $4n+3$. You can find an elementary proof of this result in Hardy & Wright, or at the end of my answer to this question if you are familiar with the basics of $L$-functions.
Hence if we choose $b$ to be, say, only divisible by primes of the form $4n+1$ then we can find arbitrarily many solutions to $2b^2=a^2+b^2$ with a fixed value of $b$.
Note: It can be easily shown that every solution of $a^2+c^2=2b^2$ in integers can be written as
$$(n^2-2mn-m^2)^2+(n^2+2mn-m^2)^2 = 2(m^2+n^2)^2$$
where $m,n$ are any integers. In particular, $b$ is itself always a sum of two squares, and we see that the number of solutions of $a^2+c^2=2b^2$ for a fixed value of $b$ is at most equal to the number of ways of writing $b$ as a sum of two squares.
