Consider 
$H\_\*(X\wedge Y;Z)$, where $X=Y=BZ/2$ for concreteness' sake. If we write $e\_i$ the generator of $H\_i(BZ/2;Z/2)$., we see that the $E\_2=E\_{\infty}$ term of the Bockstein spectral sequence
for $X\wedge Y$ is trivial, thus the permanent cycles are the image of $\beta$. As we have
$$\beta e\_{2j}=e\_{2j-1}$$ one can choose a basis of the set of permanent cycles as 
$$\{e\_{2i-1}\otimes e\_{2j-1}\}\cup \{e\_{2i}\otimes e\_{2j-1}+e\_{2i-1}\otimes e\_{2j}\}.$$
If we still denote by $e\_{2i-1}$ the lift of $e\_{2i-1}\in H\_{2i-1}(BZ/2;Z/2)$ to $H\_{2i-1}(BZ/2;Z)$, then we see that the elements of the first set lifts simply to the "products $e\_{2i-1}\otimes e\_{2j-1}$" of integral homology classes, whereas those in the second lift to "some sort of Massey products $\langle e\_{2i-1} , 2, e\_{2j-1}\rangle $".


Now, my questions are


1. Is there any reference for this kind of facts, that is the description of the homology of the product of spaces using "Massey products"?
2. Is there a setting in which one can "really" consider the obvious lifts of the elements $e\_{2i}\otimes e\_{2j-1}+e\_{2i-1}\otimes e\_{2j}$ as Massey product?



Homology of the product of spaces with integer coefficients and the Massey products
Here's a very general form. Suppose that we have six chain complexes $A\_0, A\_1, A\_2, A\_{01}, A\_{12}, A\_{012}$, with bilinear "multiplication" pairings of chain complexes:
$$
\begin{align\*}
A\_0 \otimes A\_1 &\to A\_{01}\\
A\_1 \otimes A\_2 &\to A\_{12}\\
A\_{01} \otimes A\_2 &\to A\_{012}\\
A\_0 \otimes A\_{12} &\to A\_{012}
\end{align\*}
$$
Let's assume that we also have a chain homotopy $H$ between the two composite maps $A\_0 \otimes A\_1 \otimes A\_2 \rightrightarrows A\_{012}$. Given cycles $x \in A\_0, y \in A\_1, z \in A\_2$ of degrees $p$, $q$, and $r$ respectively such that $[x \cdot y] = 0 \in H\_{p+q}(A\_{01})$ and $[y \cdot z] = 0 \in H\_{q+r}(A\_{12})$, then we can form a Massey product: choose elements $u$ and $v$ such that $\partial u = x \cdot y$ and $\partial v = y \cdot z$, and form
$$
\langle x,y,z\rangle = u \cdot z - (-1)^p x \cdot v + H(x\otimes y \otimes z)
$$
which represents an element in $H\_{p+q+r+1}(A\_{012})$. This depends in the usual way on the choices of $u$ and $v$, so there is indeterminacy. (If you use cohomological indexing then for elements in cohomological degrees $p$, $q$, and $r$ you get a Massey product in $H^{p+q+r-1}(A\_{012})$.)


The "usual" definition of a Massey product is when all six chain complexes are equal to a single complex $A$ equipped with a strictly associative multiplication (allowing us to choose $H = 0$) that makes it into a differential graded algebra.


In the case you describe, we can take


* $A\_0 = A\_{01} = C\_\*(X,pt)$ the relative singular chain complex of $X$,
* $A\_2 = A\_{12} = C\_\*(Y,pt)$ the relative singular chain complex of $Y$,
* $A\_{012} = C\_\*(X \wedge Y, pt)$ the relative singular chain complex of $X \wedge Y$, and
* $A\_1 = \Bbb Z$, the constant chain complex $\Bbb Z$ in degree 0.


Then two of the multiplications are just the unit isomorphism and the other two multiplications $C\_\*(X,\*) \otimes C\_\*(Y,\*)$ are induced by the Eilenberg-Zilber shuffle map. The unit isomorphism is associative and so we can take $H = 0$ in this case, recovering your formula.


If you like, the entire framework of these six objects, pairings, homotopy, and chosen elements asks for a DG-category (really an $A\_\infty$ DG-category) with four objects and chosen maps
$$
a \stackrel{z}{\to} b \stackrel{y}{\to} c \stackrel{x}{\to} d
$$
so that the double composites are null; in this description we can also think of this secondary operation $\langle x,y,z\rangle$ as related to Toda's "bracket" construction. The case of a differential graded algebra is the one-object DG-category case.


