I have been recently following works studying axion field. In many works studying axions, the authors assume the axion field is a bunch of plane wave, for example in this paper, equation 7 (though this paper is about vector particle, but it's similar for scalar particle like axion).
Typically, the authors will say the total axion field is something like the following $$a_{total}(t,\vec{x})=\sum_{i}a_i=\sum_{i}A_i\sin(\omega_i t-\vec{k}_i\cdot\vec{x}+\phi_i)$$ Where $A_{i,0}$ is the amplitude of the plane wave while $\phi_i$ is some random phase. $\omega_i=\sqrt{k_i^2+m_a^2}$
Typically, the authors claim $a_i$ to be the wave function of the $i^{th}$ axion particle (as claimed in the above cited paper), i.e. $A_i\sin(\omega_i t-\vec{k}_i\cdot\vec{x}+\phi_i)$. But it confuses me what does it mean to add together the wave functions of a bunch of independent particles. According to quantum mechanics, if these are the wave functions of different particles, shouldn't the total wave function be the direct product of them (not addition)? If on the other hand, $A_i\sin(\omega_i t-\vec{k}_i\cdot\vec{x}+\phi_i)$ is not the wave function of a single axion particle, what each plane wave physically represent?
Edit for the bounty: We know that there is coherent state $| C\rangle$ for the field, that gives a plane wave field expectation value (For construction of the state in field theory, see this paper), i.e. $\langle C_i|a(x,t)|C_i\rangle=A_i\sin(\omega_i t-\vec{k}_i\cdot\vec{x}+\phi_i)$. Then does this mean the field is in a mixed state described described by the density matrix $\rho=\sum_{i}|C_i\rangle \langle C_i|$?