Let $\{a_{n,k}\}_{\mathbb{N}\times\mathbb{N}}$ a double sequence of positive real numbers. Is the following equality $$\sum_{(n, k)\in\mathbb{N} \times\mathbb{N}}a_{n,k}=\sum_{n\in\mathbb{N}}\sum_{k \in \mathbb{ N }}a_{n,k}$$ a definition or a theorem?
Note. I wrote that the $a_{n,k}$ are nonnegative. Perhaps the condition of absolute convergence only serves to reverse the order of the series.
You can set things up either way. The problem with taking this as a definition is that it breaks the symmetry between the first and second copy of $\mathbb{N}$: in general, if $a_{n, k}$ is a double sequence of real numbers (not necessarily positive), it can happen that the sums $\sum_n \sum_k a_{n, k}$ and $\sum_k \sum_n a_{n, k}$ both converge but are not equal. So picking one of these as the definition of the double sum is, to my mind, conceptually unsatisfying. What would be more conceptually satisfying, to my mind, is a definition of the sum
$$\sum_{x \in X} a_x$$
where $a_x$ is indexed by elements of a set $X$, which does not depend on any auxiliary choice of something like an ordering on $X$, and so which manifestly preserves all symmetries (meaning that the sum can be reindexed arbitrarily and the value of the sum doesn't change). For a sum of non-negative real numbers the following definition can be used:
$$\boxed{ \sum_{x \in X} a_x = \sup_{S \subseteq X, |S| < \infty} \sum_{x \in S} a_x }.$$
In other words, we take the sum to be the supremum of the sum over all finite subsets. This is also the definition you get from specializing the definition of Lebesgue integration to the case of integrating a function with respect to counting measure.
This definition cannot safely be applied to a sum of possibly negative real numbers. Instead we can follow the usual Lebesgue integration trick of requiring everything to converge absolutely: we separate out the positive and negative terms, take their sums separately (in particular we require that both sums converge separately), then add them.
With these definitions the desired result is a theorem, Fubini's theorem, and furthermore because we have preserved all symmetries, if Fubini's theorem is true for $\sum_n \sum_k$ it must also be true for $\sum_k \sum_n$ and so the two sums must be equal, because we can just reindex $a_{n, k}$ and we must get the same sum over $\mathbb{N} \times \mathbb{N}$. (If the $a_{n, k}$ are allowed to be negative this is still true but it requires that the one-variable sum also be defined via "Lebesgue summation" above which means convergence must be absolute.)
It is a theorem.
Take a bijection $f:\mathbb N\to\mathbb N^2$, and suppose that $$\sum_m a_{f(m)}\tag{1}$$ exists and finite. Since the terms $a_{i,j}$ are all positive, (1) converges absolutely.
By the Dirichlet Rearrangement Theorem, for every bijection $g:\mathbb N\to \mathbb N$, $$\sum_m a_{g(f(m))}=\sum_m a_{f(m)}\tag{2}.$$
It follows that for every bijection $h:\mathbb N\to \mathbb N^2$, we may write $h=f\circ g$ for some bijection $g:\mathbb N\to \mathbb N$. Hence $$\sum_m a_{f(m)}=\sum_m a_{h(m)},\qquad \text{for every bijection } h:\mathbb N\to \mathbb N^2.\tag{3}$$ Therefore it makes sense to define $$\sum_{(n,k)\in \mathbb N^2} a_{n,k} = \sum_m a_{f(m)}$$.
When the terms are nonnegative, very simple arguments using basic properties of monotone sequences suffice to prove that (1) the double series converges if and only if (either) iterated series converges, and (2) the sums are equal.
The results (1) and (2) are also both true if the terms are not nonnegative but the convergence is absolute. The proof in that case is somewhat more difficult.
Suppose the double series $\sum_{(j, k)\in\mathbb{N} \times\mathbb{N}}a_{j,k}$ converges to $S$. By definition, for any $\epsilon > 0$ there exists a positive integer $M(\epsilon)$ such that if $m,n > M(\epsilon)$ then
$$S- \epsilon < \sum_{j=1}^m \sum_{k=1}^n a_{j,k} \leqslant S$$
For each $j$ we have $\sum_{k=1}^n a_{j,k} \leqslant S$ and since that partial sum is increasing, $\lim_{n\to \infty} \sum_{k=1}^n a_{j,k}$ exists.
Thus, for all $m > M(\epsilon)$ it follows that
$$S- \epsilon \leqslant \lim_{n\to\infty} \sum_{j=1}^m \sum_{k=1}^n a_{j,k}= \sum_{j=1}^m\sum_{k=1}^\infty a_{j,k}\leqslant S,$$
which proves that
$$\sum_{j=1}^\infty \sum_{k=1}^\infty a_{j,k}=S = \sum_{(j, k)\in\mathbb{N} \times\mathbb{N}}a_{j,k}$$
On the other hand, suppose that $\sum_{j=1}^\infty \sum_{k=1}^\infty a_{j,k}=S$. Assuming WLOG that $m \geqslant n$, we have
$$\tag{*}S_{n,n} :=\sum_{j=1}^n\sum_{k=1}^n a_{j,k} \leqslant \sum_{j=1}^m\sum_{k=1}^n a_{j,k}\leqslant \lim_{m \to \infty} \lim _{n \to \infty}\sum_{j=1}^m\sum_{k=1}^n a_{j,k} = S$$
Again since the terms $a_{j,k}$ are nonnegative, the sequence of square sums $S_{n,n}$ is increasing and bounded by (*) -- and, therefore, convergent to a limit $S' \leqslant S$. This means that for any $\epsilon >0$ there exists a positive integer $M(\epsilon)$ such that if $m,n > M(\epsilon)$, then $|S_{n,n} -S'| < \epsilon$ and $|S_{m,m} -S'| < \epsilon$. Furthermore, with $m,n > M(\epsilon)$ we have
$$\min(S_{n,n},S_{m,m}) \leqslant \sum_{j=1}^m \sum_{k=1}^n a_{j,k}\leqslant max(S_{n,n},S_{m,m})$$
and it follows that $\left|\sum_{j=1}^m \sum_{k=1}^n a_{j,k}-S'\right| < \epsilon$. Whence, the double series $ \sum_{(j, k)\in\mathbb{N} \times\mathbb{N}}a_{j,k}$ converges to $S'$. Using the previous argument, the double series and the iterated series must have the same sum $S' = S$.
