# Binaural Auralization of Room Acoustics with a Highly Scalable Wave-Based Acoustics Simulation
## 2 Theory
### Massively Parallel Room Acoustics Solver  
#### 2.1.1 Dispersion-Optimized Explicit TD-FEM  
Room-acoustics simulation in a time domain solves the following second-order scalar wave equation related to sound pressure \(p\) as  
\[\frac{\partial^{2}p(\mathbf{r},t)}{\partial t^{2}}-c_{0}^{2}\
abla^{2}p(\mathbf{r},t)= \rho_{0}c_{0}^{2}\frac{\partial q(\mathbf{r},t)}{\partial t}\delta(\mathbf{r}-\mathbf{r}_{ \text{s}}), \tag{1}\]  
with the following three boundary conditions (BCs): The rigid BC \(\Gamma_{0}\), the vibration BC \(\Gamma_{\text{V}}\), and the impedance BC \(\Gamma_{Z}\). \[\frac{\partial p(\mathbf{r},t)}{\partial n}= \begin{cases}0&\text{on }\Gamma_{0}\\ -\rho_{0}\bar{v}_{\text{n}}(\mathbf{r},t)&\text{on }\Gamma_{\text{V}}\\ -\frac{1}{c_{0}}\int_{-\infty}^{t}\ddot{y}(\mathbf{r},t-\tau)\dot{p}(\mathbf{r},\tau) \mathrm{d}\tau&\text{on }\Gamma_{Z}\end{cases}  \tag{2}\]  
Here, \(q\) and \(v_{\text{n}}\) respectively denote the volume velocity per unit volume and the vibration velocity in a normal direction on the boundary surface. \(\
abla^{2}\) and \(\delta\) respectively express the Laplacian operator and the Dirac delta function. In addition, \(c_{0}\) and \(\rho_{0}\) respectively represent the speed of sound and the air density. In addition, \(\mathbf{r}\) and \(\mathbf{r}_{\text{s}}\) respectively stand for the arbitrary coordinate vector in a computational domain and the coordinate vector of a point source. Symbol \(\cdot\) is the first-order time derivative of \(\partial/\partial t\). \(\ddot{y}\) stands for the specific acoustic admittance ratio of materials in the time domain. Impedance BC involves the convolution to address the frequency dependency of \(\ddot{y}\). Applying [PERSON] FE approximation into Equation (1) yields semi-discretized FE equation as  
\[\mathbf{M}\ddot{\mathbf{p}}+c_{0}^{2}\mathbf{K}\mathbf{p}+c_{0}\mathbf{C}(\ddot{y}*\dot{\mathbf{p}})= \mathbf{f}. \tag{3}\]  
Therein, \(\mathbf{M},\mathbf{K}\), and \(\mathbf{C}\) respectively represent the global consistent mass matrix, the global stiffness matrix, and the global dissipation matrix composed with each element matrix \(\mathbf{M}_{\text{e}}\), \(\mathbf{K}_{\text{e}}\), and \(\mathbf{C}_{\text{e}}\). \(\mathbf{p}\) and \(\mathbf{f}\) respectively represent the vectors for sound pressure and external forces. Symbols \(\cdot\) and \(*\) respectively symbolize the second-order time derivative \(\partial^{2}/\partial t^{2}\) and convolution operator. Dispersion-reduced TD-FEM [23; 24] solves Equation (3) directly to compute RIRs, but the resulting time-marching scheme shows an implicit algorithm, which requires the solution of a linear system of equations at each time step. The dispersion-optimized explicit TD-FEM [9; 10] solves the following simultaneous first-order ODEs, which is equivalent to Equation (3), with approximation of convolution relevant to the impedance BC by the auxiliary differential equations (ADE) method [25]:  
\[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\] \[\mathbf{D}\dot{\mathbf{v}}=\mathbf{f}-c_{0}^{2}\mathbf{K}\mathbf{p}-c_{0}\mathbf{C}(y_{ \infty}\dot{\mathbf{p}}+\sum_{i=1}^{n_{\text{cp}}}A_{i}\mathbf{\phi}_{i}+2\sum_{i=1}^{n _{\text{cp}}}(B_{i}\mathbf{\psi}_{i}^{(1)}+C_{i}\mathbf{\psi}_{i}^{(2)})),\] (5) \[\dot{\mathbf{\phi}}_{i}+\lambda_{i}\mathbf{\phi}_{i}=\dot{\mathbf{p}},\] (6) \[\dot{\mathbf{\psi}}_{i}^{(1)}+\alpha_{i}\mathbf{\psi}_{i}^{(1)}+\beta_{i} \mathbf{\psi}_{i}^{(2)}=\dot{\mathbf{p}},\] (7) \[\dot{\mathbf{\psi}}_{i}^{(2)}+\alpha_{i}\mathbf{\psi}_{i}^{(2)}-\beta_{i} \mathbf{\psi}_{i}^{(1)}=\mathbf{0}. \tag{8}\]  
In the equations mentioned above, \(\mathbf{D}\) represents the lumped mass matrix converted from \(\mathbf{M}\) using the row-sum method [26]. Additionally, \(\mathbf{v}\) denotes the auxiliary vector equivalent to the first derivative of \(\mathbf{p}\) with respect to time. The ADE method efficiently computes the convolution integral using an approximated acoustic admittance function form \(\hat{y}\) at angular frequency \(\omega\) in the frequency domain.