## 2 Principles of Parametric Array Loudspeakers  
As shown in Figure 1, in a typical scenario for reproducing an audible sound at frequency \(f_{\text{d}}\) with a PAL, the emitted signal can be expressed as  
\[p(0,t)=\sum_{n=1,2}P_{n}\cos(2\pi f_{n}t), \tag{1}\]  
where \(n=1,2\) denotes the index of the frequency components and \(P_{n}\) and \(f_{n}\) denote the initial sound pressure and the frequency of the corresponding component, respectively. In addition, it should be noted that \(|f_{1}-f_{2}|=f_{\text{d}}\).  
Because the attenuation is exponential, the ultrasound propagating along the propagation axis can be given as  
\[p(r,t)=\sum_{n=1,2}P_{n}e^{-a_{n}r}\cos(2\pi f_{n}t-k_{n}r), \tag{2}\]  
where \(r\) denotes the distance to the PAL and \(\alpha_{n}\) and \(k_{n}\) denote the attenuation rate and wavenumber of the frequency component indexed with \(n\), respectively.  
Based on [PERSON]'s theory [24], a virtual source at frequency is generated as the product of nonlinear interactions [34], whose intensity can be given as  
\[q(r,t)=\frac{\beta}{\rho^{2}c^{4}}\frac{\partial}{\partial t}p^{2}(r,t), \tag{3}\]  
where \(\beta\) and \(\rho\) denote the nonlinear coefficient and intensity of the air, respectively, and \(c\) denotes the speed of sound in air. By substituting Equation (2) into Equation (3), we find that the virtual source holds multiple frequency components, including \(f_{1}+f_{2}\) and \(|f_{1}-f_{2}|\). Here, we focus on the difference frequency \(f_{\rm d}=|f_{1}-f_{2}|\) as it will become the sound reproduced by a PAL with self-demodulation. The intensity of the virtual source at the difference frequency can be given as  
\[q_{\rm d}(r,t)=\frac{2\pi f_{\rm d}\beta P_{1}P_{2}}{\rho^{2}c^{4}e^{(\alpha_{ 1}+\alpha_{2})r}}\sin(2\pi f_{\rm d}t-k_{\rm d}r). \tag{4}\]  
Under [PERSON]'s theory [24], the demodulated sound can be regarded as a volume integral of virtual sources at the difference frequency \(f_{\rm d}\). Considering a far-field observation point at distance \(R_{0}\), we can approximate the volume distribution of the virtual source further as a line distribution along the propagation axis. As shown in Figure 1, the collimated ultrasound can be approximated as a cylinder with a base area \(S\) and height \(L\). In this case, the demodulated sound can be given as  
\[p_{\rm d}(R_{0},t) = \frac{\rho}{8\pi}\iiint\frac{e^{-j\delta_{\rm d}R(x,y,z)}}{R(x,y, z)}\frac{\partial q_{\rm d}(r,t)}{\partial t}dxdydz+\text{c.c.} \tag{5}\] \[\approx \frac{\rho S}{8\pi}\int_{0}^{L}\frac{e^{-j\delta_{\rm d}R(r)}}{R( r)}\frac{\partial q_{\rm d}(r,t)}{\partial t}dr+\text{c.c.},\]  
where \(j\) denotes imaginary unit \(\sqrt{-1}\) and c.c. denotes a complex conjugate. Equation (5) shows that the observed demodulated sound can be regarded as an accumulation of the sounds generated from a huge number of continuously distributed virtual sources. Therefore, in the proposed method, we extend this consideration and approximate the observed demodulated sound to an accumulation of sounds generated from finite discrete virtual sources.  
Figure 1: Geometry of a parametric array loudspeaker, an on-axis virtual source, and a far-field point for observing demodulated sound.