Trying to understand the introduction of the quantum formalism I have come across this question myself. I paste here the answer that I have given to myself and that convinces me the most. Please, if there are errors, let me know so I can continue learning:
We start by assuming a classical setup:
- we have pure and mixed states (which are probability distributions of pure states).
- the observables are real-valued functions, which can be seen like part of a c-star algebra $\mathcal{A}$. We can reverse our point of view and think of all this like if states were functionals on $\mathcal{A}$, specifically they are normalized, positive, linear functionals. That is, a mixed state $\omega$ is a probability distribution and an observable $f$ is a random variable, so we consider the pairing $\omega(f)=E_{\omega}(f)$ (the expected value). The pure states are those functionals which are also multiplicative (degenerated distributions, with all the probability concentrated on a single point). We decide to stick to this point of view, since it is equivalent to the original one (this is proven by the Gelfand-Naimark theorem and Riesz-Markov representation theorem), and because at the end of the day observables are what we manage in real life. If you can't measure the difference between two systems, you have no right to treat them as different. By the way, in mathematics there are lots of examples where a space is recovered from its algebra of functions, or at least from its sheaves. This is the spirit of Algebraic Geometry.
Now, some experimental facts, like for example the discrete emission spectrums, the Stern-Gerlach experiment, and so on, led people like Heisenberg to conclude that there were observables $A,B$ such that the mean square deviations $\Delta_{\omega}(A),\Delta_{\omega}(B)$ cannot be simultaneously reduced. For example, it was checked that for the position $q$ and momentum $p$ of electrons in any state $\omega$ it is satisfied the relation $$ \Delta_{\omega} (q)\Delta_{\omega} (p)\geq h/4\pi \equiv\hbar/2, \tag{1} $$ named Heisenberg uncertainty relation (think that when we watch an electron we modify its momentum). But this doesn't fit the description above, since in theory we could refine the apparatus and obtain a pure state $\omega$ such that $\Delta_{\omega} (q)=\Delta_{\omega} (p)=0$.
On the other hand, without having anything to do with this, observe the following algebraic development. Let $A=A^*$ and $B=B^*$ be two self-adjoint elements of an arbitrary c-star algebra such that, without lost of generality, satisfy $\omega(A)=\omega(B)=0$, being $\omega$ a pure state (normalized, positive, linear functional). For $\lambda \in \mathbb R$ we have $(A-i \lambda B)(A+i \lambda B)$ is a positive observable and positivity of $\omega$ implies $$ \omega\left(A^2\right)+|\lambda|^2 \omega\left(B^2\right)+i \lambda \omega([A, B]) \geq 0, $$ where we define, right now, $[A,B]:=AB-BA$. But if we think of this expression as a polynomial in $\lambda$ with real coefficients then we have: $$ 4 \omega\left(A^2\right) \omega\left(B^2\right) \geq|\omega(i[A, B])|^2, $$ and then $$ \Delta_\omega(A) \Delta_\omega(B) \geq \frac{1}{2}|\omega([A, B])|. \tag{2} $$
So a possible explanation for the Heisenberg relation (1) could be that the observables do not belong to a commutative c-star algebra but to a non commutative one!!!! We can introduce the canonical commutation relations in the algebra to explain that empirical fact.
What is happening here? As mathematicians (or as human beings), we elaborate theories assuming some things (and deducing others). At the beginning (classical physics), we chose modelling "experimental data outputs", or observables, as a commutative algebra: the functions defined on the phase space. This was the assumption. Since observables are tightly related to measurement apparatuses we thought that the operations of the algebra $\mathcal A$ should be given by the corresponding real number operations with the measured values.
But new experimental facts revealed that reality would be better modelled if we drop the "commutative product" hypothesis. We have that the observables "xposition" and "yposition" can be simultaneously measured, in real life, with precision enough, so we can let them behave "classically" from the operational point of view. But since "xposition" and "xmomentum" cannot be simultaneously measured with precision (as pointed out by Heisenberg (1)) then we should define multiplication in $\mathcal A$ in such a way that they do not commute, and then relation (2) would explain (1). In other words, the product of observables cannot be the pointwise product of the obtained values, but some other operation which only when restricted to the particular case of compatible observables is the pointwise product. Does this product have an interpretation in terms of the output of the measurement apparatuses? I am not aware of any.
The sum remains being commutative because we don't have the necessity of remove this hypothesis, and we can still stick to the idea of observables living inside a big c-star algebra. In case of compatible observables like "xposition" and "yposition", the sum should be what you naturally would think: a new observable whose meaning is self evident. All the observables compatible with, for example, "xposition" constitutes a subalgebra (let's denote it by $\mathcal{A}_{xposition}$) where we can perform sums and products and still remain inside. Now consider "xposition" a "xmomentum", for example. The sum is well defined even as a new observable, since the sum of self adjoint elements is self adjoint. But, what is its meaning?
To my understanding, the sum and product of observables belonging to "different subalgebras" is not naturally defined. It depends on the "big algebra" $(\mathcal A,+,\cdot)$ we choose. But keep an eye, this choice has to reconcile what we know from experiments:
- The sum and product on $\mathcal A$ are such that when restricted to any subalgebra $\mathcal A_B=C_{\mathcal A} (\{B\})$, where $B\in \mathcal A$ and $C$ denotes the centralizer, the sum and product of observables belonging to $\mathcal A_B$ corresponds to the sum and product of the outputs when they are interpreted as apparatuses.
- Some relations for elements belonging to different subalgebras, like $q\cdot p-p\cdot q=i\hbar$, have to be introduced "by hand", in order to algebraically justify Heisenberg uncertainty relations.
So the answer to "Why is addition of observables in quantum mechanics commutative?" is:
Because so far there has been no need to eliminate that hypothesis. But we cannot even be sure that we have more than several commutative algebras $\mathcal A_B$ tied together with some commutation relations. That is, not only can we not be sure that the sum of two arbitrary elements is commutative, but we cannot even be sure that such a sum exists. But as I understand it, at the moment there is no empirical contradiction in assuming it.
The following is a link to the same text but with hyperlinks that maybe help to follow the reasoning.