# 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
The function form is expressed with the real poles \(\lambda_{i}\), the pairs of complex conjugate poles \(\alpha_{i}\pm\mathrm{j}\beta_{i}\) with the imaginary unit j, and real-valued parameters \(y_{\infty}\), \(A_{i}\), \(B_{i}\) and \(C_{i}\) as follows:  
\[\hat{y}(\omega)\approx y_{\infty}+\sum_{i=1}^{n_{\text{cp}}}\frac{A_{i}}{ \lambda_{i}+\mathrm{j}\omega}+\sum_{i=1}^{n_{\text{cp}}}(\frac{B_{i}-\mathrm{ j}C_{i}}{\alpha_{i}-\mathrm{j}\beta_{i}+j\omega}+\frac{B_{i}+\mathrm{j}C_{i}}{ \alpha_{i}+\mathrm{j}\beta_{i}+\mathrm{j}\omega}). \tag{9}\]  
Therein, \(n_{\text{rp}}\) and \(n_{\text{cp}}\) are the numbers of real poles and complex conjugate pole pairs. It is crucially important that the approximated admittance with multi-pole model be causal and passive for reasons of stability. The dispersion-optimized explicit TD-FEM minimizes the dispersion error at a specific frequency for sound waves propagating in the axial and diagonal directions of the Cartesian system with special integration points for calculating \(M_{\text{e}}\) and \(K_{\text{e}}\). The time-marching scheme is formulated using a dissipation-free three-step time integration method [10] to Equation (4), the first-order backward difference to \(\hat{v}\) and the first-order central difference to \(\hat{p}\) in Equation (5), and the Crank-Nicolson method to Equations (6)-(8). The complete expression of the time-marching scheme is given as  
\[\mathbf{p}^{n}=(2\mathbf{p}^{n-1}-2\mathbf{p}^{n-2}+\mathbf{p}^{n-3})\] \[+\Delta t\mathbf{D}^{-1}(w_{\text{opt}}\mathbf{M}\mathbf{v}^{n-1}+(1-2w_{\text {opt}})\mathbf{M}\mathbf{v}^{n-2}+w_{\text{opt}}\mathbf{M}\mathbf{v}^{n-3}), \tag{10}\] \[\left(\mathbf{I}+\frac{w_{\text{opt}}}{2}c_{0}\Delta t(y_{\infty}+ \sum_{i=1}^{n_{\text{rp}}}\frac{A_{i}\Delta t}{2+\Delta t\lambda_{i}}+2\sum_{ i=1}^{n_{\text{cp}}}\frac{B_{i}(2\Delta t+\alpha_{i}\Delta t^{2})+C_{i}\beta_{i} \Delta t^{2})}{4+4\alpha_{i}\Delta t+(\alpha_{i}^{2}+\beta_{i}^{2})\Delta t^{ 2}}\right)\mathbf{D}^{-1}\mathbf{C}\mathbf{D}^{-1}\mathbf{M}\right)\mathbf{v}^{n}\] \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]0 (11) \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]1 (12) \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]2 (13) \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]3 (14) \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]4\[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]5 \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]6 (16) \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]7 (17) \[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]8  
Therein, \(\Delta t\) denotes the time interval; in addition, \(n\) stands for the time step. \(w_{\rm opt}\) is the optimization parameter for minimizing temporal dispersion error at a specific frequency. In later numerical experiments, dispersion error optimizations are conducted following earlier reports of the relevant literature [9; 10]. All equations are explicit in the time-marching scheme presented above, except for Equation (11). Equation (11) must solve implicitly for stability reasons [27]. Fortunately, Equation (11) can be solved efficiently because we only solve the equation for impedance BC's related nodes. Therefore, most of the systems are solvable explicitly as  
\[\mathbf{D}\dot{\mathbf{p}}=\mathbf{M}\mathbf{v}, \tag{4}\]9  
An iterative solver, the conjugate residual method with a convergence tolerance of \(10^{-4}\), is an effective means of solving Equation (11).