A minimalist proof of the linearity of Ethe expected value operator
Indeed, there are a lot of inexplicit steps in this calculation! Let's define the random variable $Z = X + Y$ and derive $\mathsf E[Z]$:
\begin{align} &\mathsf E[Z]\\ &=\int_\Bbb R z~p(z) ~ \mathsf d z && \text{by definition of the expected value}\\ &=\int_\Bbb R z ~ \Big( \int_\Bbb R p(z \cap y) ~ dy \Big) \mathsf d z && \text{by the law of total probability}\\ &=\int_\Bbb R z ~ \Big( \int_\Bbb R p(z | y) ~ p(y) ~ dy \Big) \mathsf d z && \text{by definition of the conditional probability}\\ &=\int_\Bbb R z ~ \Big( \int_\Bbb R p(z - y | y) ~ p(y) ~ dy \Big) \mathsf d z && \text{since p(Z=x+y|Y=y) = p(X=x|Y=y)}\\ &=\int_\Bbb R \Big( \int_\Bbb R z ~ p(z - y | y) ~ p(y) ~ dy \Big) \mathsf d z && \text{by linearity of the integral}\\ &=\int_\Bbb R \Big( \int_\Bbb R (x + y) ~ p(x | y) ~ p(y) ~ dy \Big) \mathsf d x && \text{by substituting $x = z - y$ (note that $dx = dz$)} \end{align}
By linearity of the integral, this last line is equal to $$\int_\Bbb R x \Big( \int_\Bbb R p(x | y) p(y) dy \Big) \mathsf d x + \int_\Bbb R \Big( \int_\Bbb R y ~ p(x | y) p(y) dy \Big) \mathsf d x$$$$\int_\Bbb R x ~ \Big( \int_\Bbb R p(x | y) ~ p(y) ~ dy \Big) \mathsf d x + \int_\Bbb R \Big( \int_\Bbb R y ~ p(x | y) ~ p(y) ~ dy \Big) \mathsf d x$$
Let's deal with the first term: \begin{align} \int_\Bbb R x \Big( \int_\Bbb R p(x | y) p(y) dy \Big) \mathsf d x &=\int_\Bbb R x \Big( \int_\Bbb R p(x \cap y) dy \Big) \mathsf d x && \text{conditional prob}\\ &=\int_\Bbb R x p(x) \mathsf d x && \text{total prob}\\ &=\mathsf E[X] \end{align}\begin{align} \int_\Bbb R x ~ \Big( \int_\Bbb R p(x | y) p(y) dy \Big) \mathsf d x &=\int_\Bbb R x \Big( \int_\Bbb R p(x \cap y) dy \Big) \mathsf d x && \text{conditional prob}\\ &=\int_\Bbb R x p(x) \mathsf d x && \text{total prob}\\ &=\mathsf E[X] \end{align}
Let's deal with the second term. We'll make the reasonable assumption that X and Y are measurable so as to use Fubini's theorem. This theorem allows us to change the order of integrals: \begin{align} \int_\Bbb R \Big( \int_\Bbb R y ~ p(x | y) ~ p(y) ~ dy \Big) \mathsf d x &=\int_\Bbb R y \Big( \int_\Bbb R p(y) ~ p(x | y) ~ dx \Big) \mathsf d y &&\text{with Fubini's theorem}\\ &=\int_\Bbb R y \Big( \int_\Bbb R p(x \cap y) ~ dx \Big) \mathsf d y &&\text{conditional prob}\\ &=\int_\Bbb R y ~ p(y) ~ \mathsf d y &&\text{total prob}\\ &=\mathsf E[Y] \end{align}
About $dz$
Because in my point of view I may write \begin{align}\mathsf E(X+Y)&=\int_\Bbb R (x+y)~f_{X+Y}(x+y)\mathsf d(x+y)\\&= \int_\Bbb R (x+y)~f_{X+Y}(x+y)\mathsf (dx+dy) \\&= \int_\Bbb R ~f_{X+Y}(x+y) (xdx+ydx+xdy+ydy)\\&\vdots \end{align}
Can someone explain this? Especially second step suddenly introduce one more variable becoming 2 dimensional integral and the third step replace ๐๐+๐,๐(๐ง,๐ฆ) to ๐๐,๐(๐งโ๐ฆ,๐ฆ) . I am quite confused of the change of dz -> dzdy -> dxdy, which seems like change freely...
Indeed, as you correctly stated, one cannot change the order of these notations freely, nor can they separate $d(x+y)$ into $dx + dy$.
In my opinion, you can think of $dz$ as the part of these notations that indicates what is the name of the variable that "scans" the space over which you are integrating. Since you are integrating once over $\mathbb{R}$, $d(x+y)$ is a confusing notation which gives you the impression that you could split it in two ($dx + dy$), when in fact it is not the case.
Besides, when substituting a variable for another, it's in fact a function of the variable that you are substituting. Check the wikipedia article about substitution for what I think is a clear formulation: https://en.wikipedia.org/wiki/Integration_by_substitution#Statement_for_definite_integrals
In our case, we are substituting $x$ to $\phi(z) = z - y$, and $dx$ is equal to $\frac{d\phi}{dz}(z) ~ dz = 1 \cdot dz$.