Update History.
Edit #1. The original answer did not take the condition $c \leq n$ into consideration. A new answer fixed this.
Edit #2. A minor error found in the formula for $A_n$. The net effect was that the previous answer estimated $A_{n-1}$ instead of $A_n$. This is now corrected.
Note that the $n$th term $A_n$ is given by
\begin{align*}
A_n
&= \sum_{1 \leq a < b < c \leq n} \mathbf{1}_{\{ c^2 < a^2 + b^2\}} \\
&= \sum_{1 \leq a < b \leq n} \sum_{c \in \mathbb{Z}} \mathbf{1}_{\{b+1 \leq c < \lceil \sqrt{a^2 + b^2} \rceil \text{ and } b \leq n\}} \\
&= \sum_{1 \leq a < b \leq n} \left( \left\lceil {\textstyle\sqrt{a^2 + b^2}} \right\rceil \wedge (n+1) - b - 1 \right),
\end{align*}
where $x \wedge y = \min\{x, y\}$. The presence of the term $\lceil\sqrt{a^2 + b^2} \rceil$ term suggests that lower-order terms in the asymptotic expansion of $A_n$ will be susceptible to the number-theoretic properties of the lattice points in circles, which is known to be a very hard problem. (See Gauss's circle problem.) So, I guess that there is no simple and exact formula for $A_n$.
1. Leading Term
First, we study the leading term in the asymptotic expansion of $A_n$. As $n \to \infty$,
\begin{align*}
\frac{A_n}{n^3}
&= \sum_{1 \leq a < b \leq n} \left( \biggl\lceil \sqrt{\Bigl(\frac{a}{n}\Bigr)^2 + \Bigl(\frac{b}{n}\Bigr)^2} \biggr\rceil \wedge \frac{n+1}{n} - \frac{b+1}{n} \right) \frac{1}{n^2} \\
&\to \int_{0}^{1} \int_{x}^{1} \bigl( {\textstyle \sqrt{x^2 + y^2}} \wedge 1 - y \bigr) \, \mathrm{d}y\mathrm{d}x \\
&= \frac{1}{6} - \frac{\pi}{24}.
\end{align*}
This shows that
$$ \bbox[color:navy;padding:5px;border:1px navy dotted;]{ A_n = \left( \frac{1}{6} - \frac{\pi}{24} \right) n^3 + o(n^3). } $$
2. Second Term
We now turn to extracting the second-order term. To analyze the behavior of $A_n$, we decompose it as
\begin{align*}
A_n
&= \sum_{1 \leq a < b \leq n} \left\lceil {\textstyle\sqrt{a^2 + b^2}} \right\rceil \wedge (n + 1)
- \sum_{1 \leq a < b \leq n} (b+1)
\\
&= \frac{1}{2} \sum_{1 \leq a, b \leq n} \left\lceil {\textstyle\sqrt{a^2 + b^2}} \right\rceil \wedge (n + 1)
- \frac{1}{2} \sum_{1 \leq a \leq n}\left\lceil \sqrt{2} a \right\rceil \wedge (n+1)
- \sum_{1 \leq a < b \leq n} (b+1)
\\
&= \frac{1}{8} \sum_{|a|, |b| \leq n} \left\lceil {\textstyle\sqrt{a^2 + b^2}} \right\rceil \wedge (n + 1)
- \frac{1}{2} \sum_{1 \leq a \leq n} a \\
&\qquad - \frac{1}{2} \sum_{1 \leq a \leq n} \left\lceil \sqrt{2} a \right\rceil \wedge (n + 1)
- \sum_{1 \leq a < b \leq n} (b+1). \tag{1}
\end{align*}
We first estimate easier parts. It is straightforward to check that
\begin{align*}
- \frac{1}{2} \sum_{1 \leq a \leq n} a
- \sum_{1 \leq a \leq b \leq n} (b+1)
&= -\frac{n^3}{3} + -\frac{3n^2}{4} + \frac{7n}{12}. \tag{2}
\end{align*}
Also,
\begin{align*}
- \frac{1}{2} \sum_{1 \leq a \leq n} \left\lceil \sqrt{2} a \right\rceil \wedge (n+1)
&= -\frac{n^2}{2} \left( \int_{0}^{1} ( \sqrt{2} x ) \wedge 1 \, \mathrm{d}x \right) + o(n^2) \\
&= \left( \frac{1}{4 \sqrt{2}} - \frac{1}{2} \right) n^2 + o(n^2). \tag{3}
\end{align*}
Now we turn to estimating the hardest part. Write $N(r) = \# \{z \in \mathbb{Z}^2 : \|z\| \leq r\}$ for the number of lattice points within a distance of $n$ from the origin. Then it is known that
$$ N(r) = \pi r^2 + o(r). \tag{4} $$
Using this, we decompose the first sum in $\text{(1)}$ as
\begin{align*}
&\sum_{|a|, |b| \leq n} \left\lceil {\textstyle\sqrt{a^2 + b^2}} \right\rceil \wedge (n+1) \\
&= \sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} \|z\|
+ \sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} (\left\lceil \|z\| \right\rceil - \|z\|)
+ \sum_{|a|, |b| \leq n} (n + 1) \cdot \mathbf{1}_{\{ a^2 + b^2 > n^2 \}}. \tag{5}
\end{align*}
Then the first sum in the right-hand side of $\text{(5)}$ is
\begin{align*}
\sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} \|z\|
&= \sum_{z \in \mathbb{Z}^2} \left( \int_{0}^{n} \mathbf{1}_{\{r < \|z\| \}} \, \mathrm{d}r \right) \mathbf{1}_{\{\|z\| \leq n\}}
= \int_{0}^{n} \left( \sum_{z \in \mathbb{Z}^2} \mathbf{1}_{\{r < \|z\| \leq n \}} \right) \, \mathrm{d}r \\
&= \int_{0}^{n} \left( N(n) - N(r) \right) \, \mathrm{d}r
\stackrel{\text{(4)}}= \frac{2\pi}{3} n^3 + o(n^2).
\end{align*}
Arguing similarly, we can show that $\{ \left\lceil \|z\| \right\rceil - \|z\| : z \in \mathbb{Z}^2 \text{ and } \|z\| \leq n\}$ is equidistributed on $[0, 1]$, hence the second sum in $\text{(5)}$ is
\begin{align*}
\sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} (\left\lceil \|z\| \right\rceil - \|z\|)
&= \frac{\pi}{2} n^2 + o(n^2).
\end{align*}
Finally, the last sum in $\text{(5)}$ reduces to
\begin{align*}
\sum_{|a|, |b| \leq n} (n+1) \cdot \mathbf{1}_{\{ a^2 + b^2 > n^2 \}}
&= (n+1) \left[ (2n+1)^2 - N(n) \right] \\
&= (4-\pi) n^3 + (8 -\pi) n^2 + o(n^2)
\end{align*}
Plugging these three estimates to $\text{(5)}$, we get
$$
\sum_{|a|, |b| \leq n} \left\lceil {\textstyle\sqrt{a^2 + b^2}} \right\rceil \wedge (n+1)
= \left(4 - \frac{\pi}{3} \right) n^3 - \left( 8 - \frac{\pi}{2}\right) n^2 + o(n^2). \tag{6}
$$
Then combining $\text{(2)}$, $\text{(3)}$, and $\text{(6)}$, it follows that
\begin{align*}
\bbox[color:navy;padding:5px;border:1px navy dotted;]{ A_n
= \left(\frac{1}{6} - \frac{\pi}{24} \right) n^3 + \left( - \frac{\pi}{16} + \frac{1}{4 \sqrt{2}} - \frac{1}{4}\right) n^2 + o(n^2). }
\end{align*}
3. Numerical Experiments
I tested the following ansatz for the asymptotic formula of $A_n$:
$$ A_n = c_3 n^3 + c_2 n^2 + c_{3/2} n^{3/2} + c_1 n + c_{1/2} \sqrt{n} + c_0 + E_n, \tag{7} $$
where the true values of $c_3$ and $c_2$ are
\begin{align*}
c_3 &= \frac{1}{6} - \frac{\pi}{24}
&& \approx 0.035766972767091948397\ldots, \\
c_2 &= - \frac{\pi}{16} + \frac{1}{4 \sqrt{2}} - \frac{1}{4}
&& \approx -0.26957284555272519630\ldots
\end{align*}
The graph below shows estimated valued of $c_3, c_2, c_{3/2}, c_1, c_{1/2}, c_0$ using the data $(A_1, \ldots, A_n)$ for $n$ between $100$ and $800$, using least squares fitting. In the first row, the dashed lines represent the exact value of $c_3$ and $c_2$.
This suggests that the third term in the asymptotic expansion of $A_n$ is indeed of the form $c_{3/2} x^{3/2}$ for some constant $c_{3/2}$. On the other hand, further lower-order terms seem suffering from a kind of pseudo-randomness.
Addendum: Proof of Equidistribution Theorem.
Lemma. For any Riemann integrable $f$, we have
$$ \lim_{n \to \infty} \frac{1}{\pi n^2} \sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} f(\left\lceil \|z\| \right\rceil - \|z\|)
= \int_{0}^{1} f(x) \, \mathrm{d}x. $$
Proof. We will follow the standard approach. If the claim can be proved for $f \in C(\mathbb{R}/\mathbb{Z})$, then the general claim for Riemann integrable $f$ will follow by approximating $f$ from above and below by continuous functions. Then by noting that the set of trigonometric polynomials are dense in $C(\mathbb{R}/\mathbb{Z})$ equipped with the uniform norm, the claim boils down to proving the Weyl's criterion
$$ \lim_{n \to \infty} \frac{1}{\pi n^2} \sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} e^{2\pi i \ell \|z\|}
= \mathbf{1}_{\{ \ell = 0\}} \qquad \forall \ell \in \mathbb{Z}. $$
Indeed, the sum in the limit is recast as
\begin{align*}
\sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} e^{2\pi i \ell \|z\|}
&= \sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| \leq n}} \left( 1 + \int_{0}^{\|z\|} 2\pi i \ell e^{2\pi i \ell r} \, \mathrm{d}r \right) \\
&= N(n) + \int_{0}^{n} \Biggl( \sum_{\substack{z \in \mathbb{Z}^2 \\ \|z\| < n}} \mathbf{1}_{\{ r < \|z\| \leq n\}} \Biggr) 2\pi i \ell e^{2\pi i \ell r} \, \mathrm{d}r \\
&= N(n) + \int_{0}^{n} (\pi (n^2 - r^2) + o(n) ) 2\pi i \ell e^{2\pi i \ell r} \, \mathrm{d}r \\
&= \begin{cases}
\pi n^2 + o(n^2), & \ell = 0, \\
o(n^2), & \ell \neq 0,
\end{cases}
\end{align*}
where we utilized $\text{(4)}$ in the penultimate step. Therefore the desired claim follows. $\square$
