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.$$

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.
