¹

1 Introduction↩︎

Open-domain wave propagation is a fundamental challenge in computational science and engineering because the physical domain is unbounded while the numerical computation must be performed on a finite region. Schrödinger equations, nonlinear Schrödinger equations, acoustic waves, electromagnetic waves, elastic waves, and anisotropic hyperbolic systems are naturally posed on \(\mathbb{R}^d\) or on exterior domains. In computations, however, one must replace the unbounded domain by a finite computational region. This truncation introduces an artificial boundary, and the boundary treatment can strongly affect the solution. Outgoing waves may be reflected back into the physical region, periodically wrapped around the computational box, or over-damped by an absorbing layer. These artifacts are especially damaging for dispersive and anisotropic wave equations, where the outgoing character of a wave is not a purely local property of its trace on the boundary.

Classical open-boundary techniques include absorbing boundary conditions, nonreflecting boundary conditions, perfectly matched layers (PMLs), damping or sponge layers, and rational approximations of exact radiation operators. For homogeneous isotropic wave equations, these methods can be highly effective. Absorbing and radiation boundary conditions go back to the work of Engquist–Majda, Bayliss–Turkel, and related developments [1]–[6]. The PML, introduced by Bérenger for Maxwell’s equations, is among the most widely used artificial-boundary techniques [7], [8]. It may be interpreted as a complex coordinate stretching or, more generally, as a transformation of the dispersion relation. Complex scaling ideas also have a long history in scattering theory.

However, these methods are not uniformly reliable for anisotropic or multi-branch wave systems. In particular, for certain anisotropic problems the PML can be exponentially unstable at the PDE level. This phenomenon was observed for the linearized Euler equations about a mean flow and related systems [9]–[12]. The underlying difficulty is geometric: a local absorbing layer may shift parts of the spectrum into an unstable region, and a scalar local absorber does not necessarily distinguish different wave families and their group velocities. Similarly, high-order or optimized radiation conditions can be very effective when the radiation operator has a tractable scalar or isotropic structure, but their direct extension to general anisotropic or multi-branch systems is limited. Damping layers remain useful in practice, but they are not selective: they attenuate incoming, outgoing, glancing, acoustic, vortical, or dispersive content according to location rather than according to the physical propagation direction of each wave packet.

The time-dependent phase-space filter (TDPSF) provides a different viewpoint on open boundaries [13]–[16]. Instead of imposing a local boundary condition at the edge of the computational box, TDPSF identifies portions of the numerical solution that are localized near the boundary and whose group velocities point outward, and removes only those components before they can interact with the artificial boundary. The central criterion is microlocal rather than local: \[\text{remove a packet only if it is near the boundary and } v_g(k)\cdot n>0.\] For systems, this criterion must be applied branch by branch: \[v_{g,\ell}(k)\cdot n>0, \qquad v_{g,\ell}(k)=\nabla_k\omega_\ell(k).\] This branchwise group-velocity criterion is precisely what makes phase-space filtering suitable for anisotropic wave systems where PMLs, sponge layers, or scalar radiation conditions may fail or require problem-dependent stabilization.

Physics-informed neural networks (PINNs) provide a flexible residual-minimization framework for approximating solutions of partial differential equations from equation residuals, initial data, boundary data, and possibly observations [17]. Boundary and initial conditions in PINNs are commonly imposed either weakly through penalty terms or strongly through constrained trial spaces; see, for example, [18]–[20]. For wave problems on truncated domains, absorbing-boundary and PML-type treatments have also been incorporated into PINN losses or residual formulations [21]–[23]. These approaches remain local or layer-based: they act through boundary traces, boundary residuals, or damped/PML residuals, rather than through a microlocal classification of outgoing wave packets. For wave propagation, however, this local boundary-loss strategy inherits a basic limitation: radiation is microlocal, depending on both position and frequency, whereas a local boundary penalty acts only on traces or local residuals. A boundary loss does not directly distinguish incoming, outgoing, glancing, and branch-dependent modes.

This distinction is especially important for anisotropic or multi-branch systems. In such problems, phase velocity and group velocity need not be aligned, and different dispersion branches may propagate in different directions at the same frequency. These geometric effects are also responsible for known failures of local absorbing mechanisms; for example, PMLs can be exponentially unstable for some anisotropic systems and for the linearized Euler equations about a mean flow [9]–[12]. Thus a boundary treatment based only on local traces, or on a naive scalar sign condition involving \(k\cdot n\), can remove or retain the wrong components.

The purpose of this paper is to incorporate the microlocal time-dependent phase-space-filter mechanism [13]–[16] into a PINN framework. Let \[\Omega_L=[-L,L]^d\] denote the physical computational box and let \[\Omega_{L,w}=[-L-w,L+w]^d, \qquad \mathcal{B}_w=\Omega_{L,w}\setminus\Omega_L,\] denote the extended box and its buffer layer. On each time slab \[I_m=[t_m,t_{m+1}],\] a neural network approximates the interior evolution on \(I_m\times\Omega_{L,w}\). At the end of the slab, the network output is evaluated on a grid, and a phase-space filter removes the outgoing components located in the buffer. The filtered state is then used as the initial condition for the next slab: \[u_m^{\rm filt} \xrightarrow{\text{PINN on } I_m} u_{m+1}^{\rm raw} \xrightarrow{\text{phase-space filter}} u_{m+1}^{\rm filt}.\] The method can therefore be summarized as \[\boxed{ \text{slabwise neural residual solver} + \text{time-dependent phase-space open-boundary filter}. }\]

The phase-space filtering component of the present work follows the TDPSF philosophy. Outgoing radiation is not characterized by a local boundary trace alone. Instead, the solution is decomposed into localized wave packets, coherent states, or localized Fourier components. A packet near the artificial boundary is removed only if its group velocity points outward. In the scalar case, the outgoing criterion near the side \((j,s)\) is \[v_g(k)\cdot n_j^s>0.\] For systems with several dispersion branches \(\omega_\ell(k)\), the criterion becomes branchwise: \[v_{g,\ell}(k)\cdot n_j^s>0, \qquad v_{g,\ell}(k)=\nabla_k\omega_\ell(k).\] Thus the filter uses the physically relevant propagation direction of each packet and each branch, rather than imposing a universal scalar boundary condition.

We develop the framework for three classes of equations. First, for the linear Schrödinger equation with potential, \[\mathrm i\partial_t\psi = -\frac{1}{2}\Delta\psi+V(x)\psi+f(t,x),\] we assume that the potential is constant, negligible, or slowly varying in the buffer. The outgoing classifier is then constructed from the far-field Hamiltonian. In the isotropic case, \[\omega(k)=\frac{1}{2}|k|^2+V_\infty, \qquad v_g(k)=k.\] For an anisotropic effective mass tensor \(A=A^T>0\), \[\omega(k)=\frac{1}{2} k^TAk+V_\infty, \qquad v_g(k)=Ak.\] The latter case illustrates the basic anisotropic issue: the physical wave-packet velocity is \(Ak\), not necessarily \(k\).

Second, we extend the construction to nonlinear Schrödinger equations \[\mathrm i\partial_t\psi = -\frac{1}{2}\Delta\psi + V(x)\psi + \beta|\psi|^{2\sigma}\psi + f(t,x).\] The nonlinear dynamics is approximated by the PINN residual on each time slab, while the phase-space filter is applied only in the buffer, where the outgoing field is expected to be radiative. This separation is essential. The solution may be nonlinear in the physical region, but the open-boundary classification is reliable only when the boundary field is well approximated by outgoing linear wave packets. We therefore monitor nonlinear-buffer diagnostics such as \[N_{\rm buf}(t) = \int_{\mathcal{B}_w}|\psi(t,x)|^{2\sigma+2}\,\,\mathrm dx.\] Large values of this quantity indicate that a coherent nonlinear structure, such as a soliton, is entering the buffer and that the linear outgoing-packet classifier should not be trusted without warning. This gives the method a fail-gracefully mechanism: it reports when the filtering assumptions are questionable, rather than silently producing a spurious open-boundary solution.

Third, we consider anisotropic first-order wave systems whose far-field limit has the constant-coefficient form \[\partial_t U=\mathcal{H}U, \qquad \mathcal{H}=\sum_{r=1}^d A_r\partial_{x_r}+B.\] The far-field symbol defines a Hermitian dispersion matrix \(M(k)\). Its eigenvalues \(\omega_\ell(k)\) determine the dispersion branches, and the corresponding spectral projectors \(\Pi_\ell(k)\) determine the branch decomposition. The outgoing filter at side \((j,s)\) is built from the projectors \[P_j^s(k) = \sum_\ell S_\alpha\!\left(v_{g,\ell}(k)\cdot n_j^s-\gamma\right) \Pi_\ell(k),\] where \(S_\alpha\) is a smooth approximation of the Heaviside function and \(\gamma>0\) is a glancing buffer. This branchwise construction is essential in multi-branch anisotropic systems. For example, in linearized Euler equations about a uniform mean flow, acoustic and vortical branches propagate with different group velocities, and a vortical mode may be physically transported downstream even when its wave vector points upstream.

The main contribution of this paper is not a new neural-network architecture. Rather, the contribution is the coupling of a residual-based neural interior solver with a microlocal open-boundary mechanism. This distinction is important. When a forward problem has known constant coefficients and a spectral formulation is available, a classical FFT, split-step, or spectral solver combined with a time-dependent phase-space filter remains a natural and often superior choice. The proposed Microlocal PINNFilter framework is intended for settings in which the flexibility of a neural residual formulation is useful: nonlinear equations, sparse or off-grid observations, data assimilation, variable or unknown interior coefficients, and media or geometries for which the interior evolution is not conveniently diagonalized by FFT. In this framework, the phase-space filter supplies the radiation mechanism, while the PINN supplies the flexible interior approximation.

The numerical experiments are designed to separate these roles. The linear Schrödinger benchmark shows that phase-space-filtered PINNs avoid the over-damping and wraparound produced by standard boundary penalties. The potential-scattering benchmark shows that the filter removes waves that are outgoing at the artificial boundary while preserving physical reflection generated by the potential. The anisotropic Schrödinger benchmark demonstrates that the correct classifier is based on group velocity rather than a naive \(k\cdot n\) criterion. The cubic nonlinear Schrödinger benchmarks show that the method remains effective for both defocusing and focusing outgoing wave packets and reduces persistent nonlinear activity in the buffer. The soliton stress test demonstrates the diagnostic role of the nonlinear buffer energy: weak solitons are filtered reliably, while strongly nonlinear solitons trigger a warning. The linearized Euler incoming gust benchmark shows the importance of branchwise classification: a physically incoming vortical gust is preserved at the inflow and removed only after it becomes outgoing at the outflow. Finally, the localized acoustic-defect benchmark shows the structural advantage of the PINN step in a variable-coefficient inverse setting: the constant-coefficient FFT–TDPSF solver [13]–[16] cannot represent the defect, whereas Microlocal PINNFilter learns a localized sound-speed perturbation from sparse pressure-and-velocity sensors and improves the held-out forecast.

The paper is organized as follows. Section 2 introduces the computational geometry, coherent-state notation, and localized FFT-window filters. Section 3 develops the linear Schrödinger algorithm with potential. Section 4 extends the framework to nonlinear Schrödinger equations and discusses coherent-state deletion, FFT-window filtering, and nonlinear buffer diagnostics. Section 5 treats anisotropic multi-branch wave systems and defines the branchwise group-velocity filter. The numerical benchmarks are reported in Sections 6.1–6.7.

2 Computational geometry and phase-space notation↩︎

Let \(d\geq 1\). We use the nested boxes \[\Omega_L=[-L,L]^d, \qquad \Omega_{L,w}=[-L-w,L+w]^d, \qquad \mathcal{B}_w=\Omega_{L,w}\setminus\Omega_L,\] where \(w>0\) is the buffer width. The solution of physical interest is the restriction to \(\Omega_L\), or to a smaller observation region \[\Omega_{\rm obs}\Subset \Omega_L.\] The enlarged box \(\Omega_{L,w}\) is used only for the implementation of the phase-space filter.

For each coordinate direction \(j\in\{1,\ldots,d\}\) and each sign \(s\in\{+,-\}\), we define the side buffer \[\mathcal{B}_{w,j}^{s} = \begin{cases} \{x\in\Omega_{L,w}: x_j\in[L,L+w]\}, & s=+,\\[1mm] \{x\in\Omega_{L,w}: x_j\in[-L-w,-L]\}, & s=-. \end{cases}\] The corresponding outward unit normal is \[n_j^{s} = \begin{cases} e_j, & s=+,\\ -e_j, & s=-. \end{cases}\] Thus \(j\) labels the coordinate direction of the side, while \(s\) labels whether the side is the positive or negative face of the box.

2.1 Coherent-state and FFT-window phase-space filters↩︎

We first recall the transforms used in the phase-space filter. Let \(\mathcal{F}, \mathcal{F}^{-1}\) denote the spatial Fourier transform and its inverse \[(\mathcal{F} u)(k) = \widehat u(k) = \int_{\mathbb{R}^d} e^{-i x\cdot k}u(x)\,\mathrm dx,\;\;\; (\mathcal{F}^{-1}\widehat u)(x) = \frac{1}{(2\pi)^d} \int_{\mathbb{R}^d} e^{i x\cdot k}\widehat u(k)\,\mathrm dk.\] On a finite computational box these operators are implemented by FFT and inverse FFT.

Let \(g_\sigma\) be the normalized Gaussian window \[g_\sigma(x) = (\pi\sigma^2)^{-d/4} \exp\!\left(-\frac{|x|^2}{2\sigma^2}\right).\] The associated windowed Fourier transform, or coherent-state transform, is \[(\mathcal{G}_\sigma u)(q,p) = \int_{\mathbb{R}^d} u(x)\overline{g_\sigma(x-q)}e^{-ip\cdot x}\,\mathrm dx.\] Here \(q\in\mathbb{R}^d\) is the packet center in physical space and \(p\in\mathbb{R}^d\) is the packet frequency, or momentum. The adjoint operator \(\mathcal{G}_\sigma^*\) reconstructs a function from phase-space coefficients: \[(\mathcal{G}_\sigma^* a)(x) = \int_{\mathbb{R}^d}\int_{\mathbb{R}^d} a(q,p)g_\sigma(x-q)e^{ip\cdot x} \,\mathrm dq\,\mathrm dp,\] up to the normalization constant determined by the chosen Fourier convention.

In computations, \((q,p)\) is sampled on a lattice \(\Lambda_q\times\Lambda_p.\) The reconstruction is then replaced by a discrete frame synthesis formula, \[\mathcal{G}_\sigma^* a \approx \sum_{q\in\Lambda_q} \sum_{p\in\Lambda_p} w_{q,p}\, a(q,p)\, g_\sigma(\,\cdot-q\,)e^{ip\cdot(\cdot)},\] where \(w_{q,p}\) are quadrature or frame weights. Thus, in the discrete coherent-state implementation, filtering means modifying the frame coefficients and then reconstructing the filtered state by this weighted synthesis formula.

2.2 Outgoing masks and branchwise group velocity↩︎

Let \(\omega_\ell(k)\) denote the dispersion relation of branch \(\ell\). The associated group velocity is \(v_{g,\ell}(k)=\nabla_k\omega_\ell(k).\) For scalar equations there is only one branch, and we write simply \(\omega(k), \qquad v_g(k)=\nabla_k\omega(k).\)

The basic principle of the filter is that a wave packet should be removed only if it is both close to the artificial boundary and moving out of the computational domain. In the scalar case, this condition is encoded by an outgoing phase-space mask \[\chi_{\rm out}(q,k) = \eta_{\rm bdry}(q)\, S_\delta\!\left(v_g(k)\cdot n(q)\right).\] Here \(\eta_{\rm bdry}\) is a smooth spatial cutoff supported in the boundary buffer, \(n(q)\) is the outward normal at the boundary point nearest to \(q\), and \(S_\delta\) is a smooth approximation of the Heaviside function. A typical choice is \[S_\delta(r)=\frac{1}{1+e^{-r/\delta}}.\] Thus \(\chi_{\rm out}(q,k)\approx 1\) for wave packets near the artificial boundary whose group velocity points outward, and \(\chi_{\rm out}(q,k)\approx 0\) for packets that are either away from the boundary or not outgoing.

For systems, the filtering is performed branch by branch. Let \(\Pi_\ell(k)\) be the spectral projector associated with the branch \(\omega_\ell(k)\). For a simple branch, \[\Pi_\ell(k)=d_\ell(k)d_\ell(k)^*,\] where \(d_\ell(k)\) is the normalized eigenvector satisfying \(M(k)d_\ell(k)=\omega_\ell(k)d_\ell(k)\). Near repeated eigenvalues, \(\Pi_\ell(k)\) is replaced by the projector onto the corresponding spectral cluster.

For side \((j,s)\), branch \(\ell\) is outgoing when \(v_{g,\ell}(k)\cdot n_j^s>0\), and we denote the corresponding outgoing Fourier multiplier by \(P_{\ell,j}^{s}(k)\). In the scalar case, \(P_j^s(k)=S_\delta(v_g(k)\cdot n_j^s)\), while for a system a typical branchwise multiplier is \(P_{\ell,j}^{s}(k)=S_\delta(v_{g,\ell}(k)\cdot n_j^s)\Pi_\ell(k)\). Thus \(P_{\ell,j}^{s}(k)\) selects the branch-\(\ell\) component whose group velocity points outward through side \((j,s)\).

2.3 Coherent-state filter↩︎

With the notation above, the coherent-state version of the outgoing filter is \[\label{eq:gabor-filter} \mathsf F u = \mathcal{G}_\sigma^* \big[(1-\chi_{\rm out})\mathcal{G}_\sigma u\big].\tag{1}\] Equivalently, one first decomposes \(u\) into localized wave packets, then multiplies the phase-space coefficients by \(1-\chi_{\rm out}\), and finally reconstructs the solution. Wave packets classified as outgoing are deleted, while interior, incoming, and tangential packets are retained.

2.4 Localized FFT-window filter↩︎

For higher-dimensional computations and for anisotropic systems, we use a localized FFT-window version of the same principle. Let \(\eta_j^s\) be a smooth spatial cutoff supported in the side buffer \(\mathcal{B}_{w,j}^{s}\): \[0\leq \eta_j^s\leq 1, \qquad \operatorname{supp}\eta_j^s\subset \mathcal{B}_{w,j}^{s}.\] If branch-dependent buffers or branch-dependent cutoffs are needed, we write \(\eta_{\ell,j}^{s}\) and \(\mathcal{B}_{w,\ell,j}^{s}\). Otherwise the branch index is suppressed on the spatial cutoff.

For a scalar equation, the localized outgoing component through side \((j,s)\) is \[\eta_j^s(x)\mathcal{F}^{-1} \left[ P_j^s(k)\mathcal{F}(\eta_j^s u)(k) \right].\] The corresponding FFT-window filter is \[\label{eq:simplified-filter-template-scalar} u \longmapsto u- \eta_j^s(x)\mathcal{F}^{-1} \left[ P_j^s(k)\mathcal{F}(\eta_j^s u)(k) \right].\tag{2}\] For systems, the branchwise version is \[\label{eq:simplified-filter-template-system} u \longmapsto u- \eta_{\ell,j}^{s}(x)\mathcal{F}^{-1} \left[ P_{\ell,j}^{s}(k)\mathcal{F}(\eta_{\ell,j}^{s}u)(k) \right].\tag{3}\] This should be read as follows. First, multiplication by \(\eta_{\ell,j}^{s}\) localizes the solution to the buffer layer near the side \((j,s)\). Second, the Fourier transform decomposes this localized piece into frequencies. Third, \(P_{\ell,j}^{s}(k)\) selects the outgoing branch-\(\ell\) components whose group velocity points through that side. Fourth, the inverse Fourier transform reconstructs the outgoing component in physical space. Finally, this component is multiplied again by \(\eta_{\ell,j}^{s}\) and subtracted from \(u\).

Applying this operation over all sides and branches gives the schematic filter \[\label{eq:fft-filter-all-sides-branches} \mathsf F u = u- \sum_{j=1}^d \sum_{s\in\{+,-\}} \sum_{\ell} \eta_{\ell,j}^{s}(x)\mathcal{F}^{-1} \left[ P_{\ell,j}^{s}(k)\mathcal{F}(\eta_{\ell,j}^{s}u)(k) \right].\tag{4}\] If the spatial cutoffs are independent of the branch, one simply sets \(\eta_{\ell,j}^{s}=\eta_j^s\). The FFT-window formulation is easier to implement in high dimensions than the full coherent-state frame because it uses localized FFTs and branchwise Fourier multipliers rather than a full phase-space decomposition.

2.5 Construction of the side windows↩︎

We now describe a concrete construction of the cutoff \(\eta_j^s\). A convenient choice is obtained by mollifying a strip indicator. Let \[G_\sigma(x) = \frac{1}{(\sigma\sqrt{\pi})^d} \exp\!\left(-\frac{|x|^2}{\sigma^2}\right),\] and let \(I_j^s\) be the indicator of a strip inside \(\mathcal{B}_{w,j}^{s}\). For example, for the positive side normal to the \(x_j\)-direction, set \[I_j^+(x) = \mathbf{1}_{\mathcal{S}_j^+}(x),\] where \[\mathcal{S}_j^+ = \left\{ x: x_j\in [L+w/3,L+2w/3], \quad x_m\in[-L-2w/3,L+2w/3]\;\text{for }m\neq j \right\}.\] For the negative side, set \[I_j^-(x) = \mathbf{1}_{\mathcal{S}_j^-}(x),\] where \[\mathcal{S}_j^- = \left\{ x: x_j\in [-L-2w/3,-L-w/3], \quad x_m\in[-L-2w/3,L+2w/3]\;\text{for }m\neq j \right\}.\] We then define \[\eta_j^s=G_\sigma*I_j^s.\] If necessary, the resulting function is smoothly truncated near corners or renormalized so that \(0\leq \eta_j^s\leq 1\). For systems, one may use the same window for all branches, \(\eta_{\ell,j}^{s}=\eta_j^s\), or introduce branch-dependent windows \(\eta_{\ell,j}^{s}\) if different branches require different buffer widths.

The parameter \(\sigma\) controls the space–frequency uncertainty tradeoff. A small value of \(\sigma\) localizes sharply in physical space but spreads the window in frequency. A large value of \(\sigma\) reduces frequency spreading but allows the tails of \(\eta_j^s\) to leak farther into the interior computational region. In practice, one chooses \(\sigma\), a frequency buffer \(k_b\), and a tolerance \(\delta\) so that \[\label{eq:sigma-window-bounds} C_1 k_b^{-1}\sqrt{\log(\delta^{-1})} \leq \sigma \leq C_2 \frac{w}{\sqrt{\log(\delta^{-1})}}.\tag{5}\] The lower bound limits frequency leakage at scale \(k_b\), while the upper bound keeps the spatial tails of the window inside the side buffer. Consequently one obtains the expected uncertainty scaling \[\label{eq:w-scaling} w\gtrsim k_b^{-1}\log(\delta^{-1}).\tag{6}\] The constants in 5 depend on the precise definition of the window, the spatial dimension, and the Fourier transform convention.

The nonlinear Schrödinger section below uses both viewpoints. The coherent-state deletion formulation is closest to the original time-dependent phase-space-filter idea: the numerical solution is decomposed into localized coherent states, outgoing coherent states are deleted, and the remaining coefficients are synthesized back to physical space. The FFT-window formulation is more convenient for higher dimensional and anisotropic systems, where the outgoing classification is naturally expressed by the branchwise multipliers \(P_{\ell,j}^{s}(k)\).

3 Linear Schrödinger equation with potential↩︎

We first consider the linear Schrödinger equation \[\label{eq:schrodinger-full} \mathrm i\partial_t\psi(t,x) = -\frac{1}{2}\Delta\psi(t,x)+V(x)\psi(t,x)+f(t,x), \qquad x\in\mathbb{R}^d,\tag{7}\] with initial condition \[\label{eq:schrodinger-ic} \psi(0,x)=\psi_0(x).\tag{8}\] The potential \(V\) is assumed to be real-valued. The forcing term \(f\) may be set to zero.

Under Assumption 1, the far-field dispersion relation is \[\label{eq:schrodinger-dispersion} \omega(k)=\frac{1}{2}|k|^2+V_\infty, \qquad v_g(k)=\nabla_k\omega(k)=k.\tag{9}\] In some applications the far-field Schrödinger operator may have direction-dependent dispersion. To model this situation, one may replace the isotropic Laplacian by an anisotropic constant-coefficient operator \(-\frac{1}{2}\nabla\cdot(A\nabla),\) where \(A\in\mathbb{R}^{d\times d}\) is a symmetric positive definite matrix. The matrix \(A\) represents an effective inverse mass tensor, or more generally an anisotropic dispersion tensor, in the far-field medium. Thus the corresponding anisotropic far-field equation is \[\mathrm i\partial_t\psi = -\frac{1}{2}\nabla\cdot(A\nabla\psi)+V_\infty\psi .\] For the plane wave \(\psi(t,x)=a e^{\mathrm i(k\cdot x-\omega t)}\), we obtain \(\omega(k)=\frac{1}{2} k^T A k+V_\infty\). Consequently the group velocity is \[\label{eq:anisotropic-sch-dispersion} v_g(k) = \nabla_k\omega(k) = Ak.\tag{10}\] Hence, in the anisotropic scalar Schrödinger case, outgoing packets must be classified using \(v_g(k)\cdot n_j^s=Ak\cdot n_j^s\), not merely using \(k\cdot n_j^s\). This distinction is important because the direction of energy or wave-packet propagation need not be parallel to the wave vector \(k\).

3.1 Neural representation and PDE residual↩︎

We write the neural approximation in real and imaginary parts: \[\psi_\theta(t,x)=a_\theta(t,x)+\mathrm ib_\theta(t,x),\] where \(a_\theta,b_\theta:\mathbb{R}\times\mathbb{R}^d\to\mathbb{R}\) are the two outputs of the neural network. The complex residual associated with 7 is \[\label{eq:schrodinger-complex-residual} \mathcal{R}_\theta = \mathrm i\partial_t\psi_\theta + \frac{1}{2}\Delta\psi_\theta - V\psi_\theta - f.\tag{11}\] If \(f=f_R+\mathrm if_I\), then the real and imaginary residuals are \[\begin{align} \mathcal{R}_\theta^R &= -\partial_t b_\theta +\frac{1}{2}\Delta a_\theta -Va_\theta -f_R, \tag{12}\\ \mathcal{R}_\theta^I &= \partial_t a_\theta +\frac{1}{2}\Delta b_\theta -Vb_\theta -f_I. \tag{13} \end{align}\] On a time slab \(I_m=[t_m,t_{m+1}],\) the residual loss is approximated by collocation: \[\label{eq:sch-residual-loss} \mathcal{L}_{\rm PDE}^{(m)}(\theta) = \frac{1}{N_r} \sum_{r=1}^{N_r} \left( |\mathcal{R}_\theta^R(\tau_r,x_r)|^2 + |\mathcal{R}_\theta^I(\tau_r,x_r)|^2 \right),\tag{14}\] where \((\tau_r,x_r)\in I_m\times\Omega_{L,w}.\) The collocation points should be denser in the physical box \(\Omega_L\), near the support of the potential, and in the side buffers \(\mathcal{B}_{w,j}^{s}\), \(j=1,\ldots,d\), \(s\in\{+,-\}\).

3.2 Outgoing frequency masks for the linear Schrödinger equation↩︎

For the side \((j,s)\), where \(j\in\{1,\ldots,d\}\) and \(s\in\{+,-\}\), the outward unit normal is \(n_j^s\). The outgoing condition for a scalar Schrödinger packet is \[\label{eq:sch-outgoing-condition} v_g(k)\cdot n_j^s>0.\tag{15}\] For the standard isotropic Schrödinger equation, \(v_g(k)=k\), so 15 becomes \(k\cdot n_j^s>0\), or equivalently \(s k_j>0\), where \(s=+1\) on the positive \(x_j\)-face and \(s=-1\) on the negative \(x_j\)-face.

To avoid aggressively filtering packets close to the glancing set \(v_g(k)\cdot n_j^s=0,\) we introduce a glancing buffer \(\gamma>0\) and define the smooth outgoing Fourier multiplier \[\label{eq:sch-frequency-mask} P_j^s(k) = S_\alpha\!\left(v_g(k)\cdot n_j^s-\gamma\right),\tag{16}\] where \(S_\alpha(r)=\frac{1}{1+\exp(-r/\alpha)}\) is a smoothed Heaviside function. The parameter \(\alpha>0\) controls the sharpness of the transition, while \(\gamma>0\) prevents packets near the glancing region from being removed too strongly.

For the anisotropic Schrödinger operator \(-\frac{1}{2}\nabla\cdot(A\nabla)\), the multiplier is still given by 16 , but with \(v_g(k)=Ak\). Thus the correct outgoing condition is \(Ak\cdot n_j^s>0\), not simply \(k\cdot n_j^s>0\).

3.3 Linear Schrödinger outgoing filter↩︎

Let \(u:\Omega_{L,w}\to\mathbb{C}\) be a grid function. For the side \((j,s)\), let \(\eta_j^s\) be the side window supported in \(\mathcal{B}_{w,j}^{s}\), as constructed in Section 2.5. The outgoing component through side \((j,s)\) is defined by \[\label{eq:sch-outgoing-component} \mathcal{O}_j^s u = \eta_j^s(x) \mathcal{F}^{-1} \left[ P_j^s(k)\, \mathcal{F}(\eta_j^s u)(k) \right](x).\tag{17}\] The corresponding one-side filter is \[\label{eq:sch-filter-one-side} \mathsf F_j^s u = u-\mathcal{O}_j^s u.\tag{18}\] Applying the filter over all sides gives the complete linear Schrödinger filter \[\label{eq:sch-complete-filter} \mathsf F_{\rm Sch}u = \left( \prod_{j=1}^{d} \mathsf F_j^{-}\mathsf F_j^{+} \right)u.\tag{19}\] The order of the product is not important at the formal level. Numerically, the side filters are applied sequentially. In corner regions, one may either allow the sequential filters to act or use a partition of unity to avoid double counting.

3.4 Slabwise PINN–TDPSF algorithm for the linear Schrödinger equation↩︎

Let \(0=t_0<t_1<\cdots<t_M=T, \qquad t_{m+1}-t_m\leq \Delta T_{\rm filt}.\) The filtered state at time \(t_m\) is denoted by \(\psi_m^{\rm filt}.\) Initially, \(\psi_0^{\rm filt} = \psi_0|_{\Omega_{L,w}}.\) On each slab \(I_m=[t_m,t_{m+1}]\), we either train a new neural network or warm-start from the previous slab.

The initial matching loss is \[\label{eq:sch-initial-loss} \mathcal{L}_{\rm init}^{(m)}(\theta) = \frac{1}{N_i} \sum_{r=1}^{N_i} \left| \psi_\theta(t_m,y_r)-\psi_m^{\rm filt}(y_r) \right|^2.\tag{20}\] The points \(y_r\) may be FFT grid points, random points, or interpolation points from the filtered spectral grid.

The basic slab loss is \[\label{eq:sch-slab-loss} \mathcal{L}_{\rm Sch}^{(m)}(\theta) = \mathcal{L}_{\rm PDE}^{(m)}(\theta) + \lambda_i\mathcal{L}_{\rm init}^{(m)}(\theta) + \lambda_{\rm per}\mathcal{L}_{\rm per}^{(m)}(\theta) + \lambda_{\rm reg}\mathcal{L}_{\rm reg}^{(m)}(\theta).\tag{21}\] The periodic term \(\mathcal{L}_{\rm per}^{(m)}\) is optional and is imposed only on the outer boundary of \(\Omega_{L,w}\) to stabilize the FFT-based implementation. It should not be interpreted as a physical periodic boundary condition on \(\Omega_L\). The role of the filter is to remove outgoing waves in the buffer before they reach the outer boundary of the enlarged computational box.

4 Nonlinear Schrödinger equations and the TDPSF mechanism↩︎

We now incorporate nonlinear Schrödinger equations. The model problem is \[\label{eq:nls-general} \mathrm i\partial_t\psi = -\frac{1}{2}\Delta\psi +V(x)\psi +\beta |\psi|^{2\sigma}\psi +f(t,x), \qquad x\in\mathbb{R}^d,\tag{22}\] where \(\beta\in\mathbb{R}\), \(\sigma>0\), and \(V\) is real-valued. The cubic NLS corresponds to \(\sigma=1\). More generally, the nonlinear term may be replaced by \[\label{eq:nls-general-g} g(t,x,\psi)\psi,\tag{23}\] provided that the far-field behavior in the buffer is asymptotically linear.

The TDPSF idea is especially natural for NLS. The solution may be strongly nonlinear in the physical region \(\Omega_L\), while the outgoing radiation near the artificial boundary is often well described by linear Schrödinger wave packets. Thus the nonlinear equation is solved in the extended box \(\Omega_{L,w}\), but the open-boundary treatment is performed by a phase-space criterion based on outgoing linear radiation in the buffer \(\mathcal{B}_w=\Omega_{L,w}\setminus\Omega_L.\)

4.1 PINN residual for NLS↩︎

We write the neural approximation as \[\psi_\theta=a_\theta+\mathrm ib_\theta, \qquad \rho_\theta=a_\theta^2+b_\theta^2.\] For 22 , define the complex residual \[\label{eq:nls-complex-residual} \mathcal{R}_\theta^{\rm NLS} = \mathrm i\partial_t\psi_\theta +\frac{1}{2}\Delta\psi_\theta -V\psi_\theta -\beta |\psi_\theta|^{2\sigma}\psi_\theta -f.\tag{25}\] If \(f=f_R+\mathrm if_I\), then the real and imaginary residuals are \[\begin{align} \mathcal{R}_{\theta,{\rm NLS}}^R &= -\partial_t b_\theta +\frac{1}{2}\Delta a_\theta -Va_\theta -\beta \rho_\theta^{\sigma}a_\theta -f_R, \tag{26}\\ \mathcal{R}_{\theta,{\rm NLS}}^I &= \partial_t a_\theta +\frac{1}{2}\Delta b_\theta -Vb_\theta -\beta \rho_\theta^{\sigma}b_\theta -f_I. \tag{27} \end{align}\] On the slab \(I_m=[t_m,t_{m+1}],\) the residual loss is approximated by collocation: \[\label{eq:nls-residual-loss} \mathcal{L}_{\rm PDE,NLS}^{(m)}(\theta) = \frac{1}{N_r} \sum_{r=1}^{N_r} \left( |\mathcal{R}_{\theta,{\rm NLS}}^R(\tau_r,x_r)|^2 + |\mathcal{R}_{\theta,{\rm NLS}}^I(\tau_r,x_r)|^2 \right),\tag{28}\] where \((\tau_r,x_r)\in I_m\times\Omega_{L,w}.\) The collocation points should be concentrated in the physical box \(\Omega_L\), near the support of the potential \(V\), near regions where \(|\psi_\theta|\) is large, and in the side buffers \(\mathcal{B}_{w,j}^{s}, j=1,\ldots,d,\quad s\in\{+,-\}.\)

The initial matching loss is \[\label{eq:nls-initial-loss} \mathcal{L}_{\rm init,NLS}^{(m)}(\theta) = \frac{1}{N_i} \sum_{r=1}^{N_i} \left| \psi_\theta(t_m,y_r)-\psi_m^{\rm filt}(y_r) \right|^2.\tag{29}\] Here \(\psi_m^{\rm filt}\) is the filtered state inherited from the previous slab, and \(y_r\in\Omega_{L,w}\) are grid points, random points, or interpolation points from the filtered spectral grid.

The basic NLS slab loss is \[\label{eq:nls-slab-loss} \mathcal{L}_{\rm NLS}^{(m)}(\theta) = \mathcal{L}_{\rm PDE,NLS}^{(m)}(\theta) + \lambda_i\mathcal{L}_{\rm init,NLS}^{(m)}(\theta) + \lambda_{\rm per}\mathcal{L}_{\rm per}^{(m)}(\theta) + \lambda_{\rm reg}\mathcal{L}_{\rm reg}^{(m)}(\theta).\tag{30}\] The term \(\mathcal{L}_{\rm per}^{(m)}\) is optional and is used only to make the extended box compatible with FFT-based filtering. It should be weighted mildly, because the physical boundary condition is not periodicity on \(\Omega_L\). The role of the phase-space filter is to remove outgoing radiation in the buffer before it reaches the outer boundary of \(\Omega_{L,w}\).

4.2 Coherent-state outgoing classification↩︎

Let \(\Lambda_q\subset\Omega_{L,w}\) be a position lattice and let \(\Lambda_p\subset[-K_{\max},K_{\max}]^d\) be a frequency lattice. At a filtering time, given a raw endpoint state \(u\), usually \(u=\psi_{m+1}^{\rm raw},\) we compute the coherent-state coefficients \[\label{eq:nls-coefficients} c_{q,p} = (\mathcal{G}_\sigma u)(q,p), \qquad (q,p)\in\Lambda_q\times\Lambda_p.\tag{31}\] Equivalently, if \(\phi_{q,p}^{\sigma}\) denotes the coherent state with center \(q\), momentum \(p\), and width \(\sigma\), then \(c_{q,p}=\langle u,\phi_{q,p}^{\sigma}\rangle\). To classify \(c_{q,p}\), we first determine whether \(q\) lies in the buffer \(\mathcal{B}_w\); if \(q\in\mathcal{B}_{w,j}^{s}\), we use the outward normal \(n(q)=n_j^s\). For the standard far-field Schrödinger operator, \(v_g(p)=p\), while for an anisotropic effective mass tensor \(A\), one uses \(v_g(p)=Ap\).

Given a glancing buffer \(\gamma>0\), the coefficient is classified as \[\begin{array}{lll} \text{outgoing} &\text{if}& q\in\mathcal{B}_w\;\text{and}\;v_g(p)\cdot n(q)\geq \gamma, \\[1mm] \text{incoming} &\text{if}& q\in\mathcal{B}_w\;\text{and}\;v_g(p)\cdot n(q)\leq -\gamma, \\[1mm] \text{glancing or ambiguous} &\text{if}& q\in\mathcal{B}_w\;\text{and}\;|v_g(p)\cdot n(q)|<\gamma, \\[1mm] \text{interior} &\text{if}& q\notin\mathcal{B}_w. \end{array}\] Only outgoing coefficients are deleted. Incoming, glancing, ambiguous, and interior coefficients are retained. The glancing or ambiguous coefficients are also recorded as a diagnostic.

A smooth version replaces the hard outgoing indicator by \[\label{eq:nls-smooth-outgoing-mask} \chi_{\rm out}(q,p) = \eta_{\rm bdry}(q) S_\alpha\!\left(v_g(p)\cdot n(q)-\gamma\right),\tag{32}\] where \(\eta_{\rm bdry}\) is a smooth cutoff supported in \(\mathcal{B}_w\), and \(S_\alpha(r)=\frac{1}{1+\exp(-r/\alpha)}\) is a smoothed Heaviside function. The filtered coherent-state coefficients are \[\label{eq:nls-filtered-coefficients} c_{q,p}^{\rm filt} = \big(1-\chi_{\rm out}(q,p)\big)c_{q,p}.\tag{33}\] If the coherent states form a tight frame with frame constant \(A_0\), the reconstructed filtered state is \[\label{eq:nls-coherent-reconstruction} u^{\rm filt}(x) = A_0^{-1} \sum_{(q,p)\in\Lambda_q\times\Lambda_p} c_{q,p}^{\rm filt}\phi_{q,p}^{\sigma}(x).\tag{34}\] For a non-tight frame, one uses the dual frame or solves the associated least-squares reconstruction problem.

4.3 Original-style NLS TDPSF step↩︎

The following algorithm describes the NLS phase-space-filter step in the spirit of the time-dependent phase-space filter. The interior propagator may be a spectral method, a finite difference method, or the slabwise PINN described below.

4.4 FFT-window version of the NLS filter↩︎

For high-dimensional implementations, one may replace the coherent-state deletion step by the localized FFT-window filter from Section 2.4. For the scalar NLS, define \[P_j^s(k) = S_\alpha\!\left(v_g(k)\cdot n_j^s-\gamma\right), \qquad v_g(k)=k\] for the standard far-field Schrödinger operator. For anisotropic mass tensor \(A\), use \(v_g(k)=Ak\). The outgoing component through side \((j,s)\) is \(\mathcal{O}_j^s u=\eta_j^s(x)\mathcal{F}^{-1} [P_j^s(k)\mathcal{F}(\eta_j^s u)(k)](x)\), and the corresponding one-side filter is \(\mathsf F_j^s u=u-\mathcal{O}_j^s u\). The complete NLS FFT-window filter is \[\label{eq:nls-fft-complete-filter} \mathsf F_{\rm NLS}u = \left( \prod_{j=1}^{d} \mathsf F_j^{-}\mathsf F_j^{+} \right)u.\tag{35}\] Equivalently, in additive schematic form, \[\mathsf F_{\rm NLS}u \approx u- \sum_{j=1}^d \sum_{s\in\{+,-\}} \eta_j^s(x)\mathcal{F}^{-1} \left[ P_j^s(k)\mathcal{F}(\eta_j^s u)(k) \right].\] As in the linear case, the product form means that the side filters are applied sequentially. In corner regions, one may use a partition of unity to avoid double counting.

4.5 NLS TDPSF–PINN algorithm↩︎

We now replace the generic interior propagator in Algorithm 4 by a slabwise PINN. The filtered state at the beginning of the slab \(I_m=[t_m,t_{m+1}]\) is \(\psi_m^{\rm filt}\). The network is trained by minimizing 30 . Then its endpoint is filtered by either the coherent-state NLS filter or the FFT-window NLS filter.

4.6 Practical NLS failure indicators↩︎

The NLS TDPSF–PINN algorithm should report the following warnings.

1. Large nonlinear buffer energy. If \(\int_{\mathcal{B}_w}|\psi|^{2\sigma+2}\,\mathrm dx\) is large, the far-field linear outgoing classifier is questionable.

2. Large glancing energy. If many coherent states satisfy \(|v_g(p)\cdot n(q)|<\gamma\), then the algorithm cannot confidently decide whether they are outgoing or incoming.

3. Large high-frequency tail. Since \(v_g(k)=k\) for the standard Schrödinger equation, high frequencies may cross the buffer rapidly; hence the time-step restriction must be based on the resolved maximum frequency \(K_{\max}\), or, in the anisotropic case, on \(v_{\max}=\sup_{|k|\leq K_{\max}}|v_g(k)|\).

4. Large coefficient mass near corners. Near corners, a packet may be outgoing relative to more than one side \((j,s)\). A partition of unity or a corner-specific normal should be used to avoid double counting.

5. Large mismatch after filtering. If \(\|\psi_{m+1}^{\rm raw}-\psi_{m+1}^{\rm filt}\|_{L^2(\Omega_L)}\) is not small, then the filter is affecting the physical box, indicating that the buffer is too narrow, the side window is too wide, or the slab length is too large.

These diagnostics are part of the value of the method: rather than silently producing a reflected or wrapped solution, the algorithm identifies when the phase-space classification is unreliable.

5 Anisotropic wave equations on unbounded domains↩︎

We next consider linear anisotropic wave systems. The far-field model is a constant-coefficient first-order system \[\label{eq:anisotropic-system} \partial_t U(t,x)=\mathcal{H} U(t,x), \qquad \mathcal{H}=\sum_{r=1}^d A_r\partial_{x_r}+B,\tag{36}\] where \(U:\mathbb{R}^d\to\mathbb{C}^q.\) We assume that the matrices \(A_r\in\mathbb{C}^{q\times q}\) are Hermitian and that \(B\in\mathbb{C}^{q\times q}\) is skew-Hermitian. More generally, a symmetrizable hyperbolic system \[T\partial_t U = \sum_{r=1}^d A_r\partial_{x_r}U+BU, \qquad T=T^*>0,\] can be reduced to the form 36 by the change of variables \(U\mapsto T^{1/2}U\).

The actual PDE may contain scatterers, lower-order terms, nonlinearities, or variable coefficients inside a compact region \(B_{\rm int}\Subset\Omega_L\). The filter, however, is constructed from the constant-coefficient far-field operator \(\mathcal{H}\) in the buffer \(\mathcal{B}_w=\Omega_{L,w}\setminus\Omega_L\). Thus the open-boundary mechanism is determined by the far-field dispersion relation and by the group velocities of the outgoing branches.

5.1 Dispersion branches and group velocities↩︎

For a plane wave \((U(t,x)=d(k)\exp\!\big(\mathrm i(k\cdot x-\omega t)\big),\) the Fourier symbol of \(\mathcal{H}\) is \[\widehat{\mathcal{H}}(k) = \mathrm i\sum_{r=1}^d k_rA_r+B.\] Since \(\widehat{\mathcal{H}}(k)\) is skew-Hermitian, the matrix \[\label{eq:M-k-Hermitian} M(k)=-\mathrm i\widehat{\mathcal{H}}(k)\tag{37}\] is Hermitian. Let \[\label{eq:dispersion-eigenproblem} M(k)d_\ell(k)=\omega_\ell(k)d_\ell(k), \qquad \ell=1,\ldots,q,\tag{38}\] where the eigenvectors \(d_\ell(k)\) are chosen orthonormally whenever the eigenvalues are simple. The group velocity on branch \(\ell\) is \[\label{eq:anisotropic-group-velocity} v_{g,\ell}(k)=\nabla_k\omega_\ell(k).\tag{39}\] If the eigenvalue \(\omega_\ell(k)\) is simple, then the Hellmann–Feynman formula gives \[\label{eq:hellmann-feynman} \partial_{k_r}\omega_\ell(k) = d_\ell(k)^*\partial_{k_r}M(k)d_\ell(k), \qquad r=1,\ldots,d.\tag{40}\] At repeated eigenvalues, or near crossings, one should use spectral cluster projectors rather than individual eigenvectors.

This group-velocity criterion distinguishes isotropic from anisotropic propagation: in isotropic models, phase and group velocities are aligned, whereas in anisotropic models they need not be. Therefore filtering by the sign of \(k\cdot n_j^s\) is generally incorrect; one must filter by the sign of \(v_{g,\ell}(k)\cdot n_j^s\).

5.2 Branchwise outgoing projectors↩︎

Let \(\Pi_\ell(k)=d_\ell(k)d_\ell(k)^*\) be the spectral projector onto branch \(\ell\) when the eigenvalue is simple. Near multiple eigenvalues, \(\Pi_\ell(k)\) should be replaced by the projector onto the corresponding spectral cluster.

Choose a glancing buffer \(\gamma>0\) and a smoothing width \(\alpha>0\). For branch \(\ell\), side direction \(j\), and sign \(s\), define the scalar branch mask \[\label{eq:anisotropic-branch-mask} p_{\ell,j}^{s}(k) = S_\alpha\!\left(v_{g,\ell}(k)\cdot n_j^s-\gamma\right),\tag{42}\] where \(S_\alpha(r)=\frac{1}{1+\exp(-r/\alpha)}.\) Thus \(p_{\ell,j}^{s}(k)\approx 1\) for branch-\(\ell\) frequencies whose group velocity points outward through the side \((j,s)\), and \(p_{\ell,j}^{s}(k)\approx 0\) for non-outgoing or glancing frequencies.

The branchwise outgoing Fourier multiplier is \[\label{eq:anisotropic-branch-projector} P_{\ell,j}^{s}(k) = p_{\ell,j}^{s}(k)\Pi_\ell(k).\tag{43}\] The total outgoing projector for side \((j,s)\) is \[\label{eq:anisotropic-frequency-projector} P_j^s(k) = \sum_{\ell=1}^{q}P_{\ell,j}^{s}(k) = \sum_{\ell=1}^{q} p_{\ell,j}^{s}(k)\Pi_\ell(k).\tag{44}\] Equivalently, if \(D(k)\) is the unitary matrix whose columns are the eigenvectors \(d_\ell(k)\), then \[P_j^s(k) = D(k) \operatorname{diag} \big( p_{1,j}^{s}(k),\ldots,p_{q,j}^{s}(k) \big) D(k)^*.\] The precise placement of \(D(k)\) and \(D(k)^*\) depends only on whether the eigenvectors are stored as columns or rows; throughout the paper, \(P_j^s(k)\) denotes the projector acting on the physical vector \(\widehat U(k)\in\mathbb{C}^q\).

5.3 Outgoing wave filter for anisotropic systems↩︎

Let \(U:\Omega_{L,w}\to\mathbb{C}^q\) be a grid function. For side \((j,s)\), let \(\eta_j^s\) be the side window supported in \(\mathcal{B}_{w,j}^{s}\). If branch-dependent localization is needed, we write \(\eta_{\ell,j}^{s}\) and \(\mathcal{B}_{w,\ell,j}^{s}\). Otherwise the spatial window is shared by all branches and we set \(\eta_{\ell,j}^{s}=\eta_j^s.\)

Using the total outgoing projector \(P_j^s(k)\), the outgoing component through side \((j,s)\) is \[\label{eq:anisotropic-outgoing-component} \mathcal{O}_j^s U = \eta_j^s(x)\mathcal{F}^{-1} \left[ P_j^s(k)\mathcal{F}(\eta_j^s U)(k) \right](x).\tag{45}\] Here the Fourier transform is applied componentwise to the vector field \(U\), and \(P_j^s(k)\) acts as a \(q\times q\) matrix multiplier on the Fourier coefficient \(\widehat U(k)\).

The corresponding one-side filter is \[\label{eq:anisotropic-one-side-filter} \mathsf F_j^s U = U-\mathcal{O}_j^s U.\tag{46}\] The complete anisotropic phase-space filter is \[\label{eq:anisotropic-complete-filter} \mathsf F_{\rm aniso}U = \left( \prod_{j=1}^d \mathsf F_j^-\mathsf F_j^+ \right)U.\tag{47}\] The product means that the side filters are applied sequentially. At the formal level the order is immaterial; in computations, corner regions may be treated either by sequential filtering or by a partition of unity to avoid double counting.

Equivalently, in additive schematic form, \[\mathsf F_{\rm aniso}U \approx U- \sum_{j=1}^d \sum_{s\in\{+,-\}} \eta_j^s(x)\mathcal{F}^{-1} \left[ P_j^s(k)\mathcal{F}(\eta_j^s U)(k) \right].\] If branch-dependent windows are used, this may be written as \[\mathsf F_{\rm aniso}U \approx U- \sum_{j=1}^d \sum_{s\in\{+,-\}} \sum_{\ell=1}^q \eta_{\ell,j}^{s}(x)\mathcal{F}^{-1} \left[ P_{\ell,j}^{s}(k) \mathcal{F}(\eta_{\ell,j}^{s} U)(k) \right].\]

5.4 PINN residual for anisotropic systems↩︎

Let \(U_\theta:I_m\times\Omega_{L,w}\to\mathbb{C}^q\) be a neural approximation. For the far-field model 36 , the residual is \[\label{eq:aniso-residual-constant} \mathcal{R}_\theta = \partial_tU_\theta - \sum_{r=1}^d A_r\partial_{x_r}U_\theta - BU_\theta.\tag{48}\] For a more general interior model \[\label{eq:aniso-residual-variable-model} \partial_t U = \mathcal{H}_{\rm int}(x,\partial_x)U+F(t,x),\tag{49}\] we define \[\label{eq:aniso-residual-variable} \mathcal{R}_\theta = \partial_tU_\theta - \mathcal{H}_{\rm int}(x,\partial_x)U_\theta - F.\tag{50}\] The important requirement is that \(\mathcal{H}_{\rm int}\) coincides with the constant-coefficient far-field operator \(\mathcal{H}\) in the buffer where the filter is applied.

The residual loss on the slab \(I_m=[t_m,t_{m+1}]\) is \[\label{eq:aniso-residual-loss} \mathcal{L}_{\rm PDE}^{(m)}(\theta) = \frac{1}{N_r} \sum_{r=1}^{N_r} \left| \mathcal{R}_\theta(\tau_r,x_r) \right|_{\mathbb{C}^q}^2,\tag{51}\] where \((\tau_r,x_r)\in I_m\times\Omega_{L,w}.\) The initial matching loss is \[\label{eq:aniso-initial-loss} \mathcal{L}_{\rm init}^{(m)}(\theta) = \frac{1}{N_i} \sum_{r=1}^{N_i} \left| U_\theta(t_m,y_r)-U_m^{\rm filt}(y_r) \right|_{\mathbb{C}^q}^2.\tag{52}\] The slab loss is \[\label{eq:aniso-slab-loss} \mathcal{L}_{\rm aniso}^{(m)}(\theta) = \mathcal{L}_{\rm PDE}^{(m)}(\theta) + \lambda_i\mathcal{L}_{\rm init}^{(m)}(\theta) + \lambda_{\rm per}\mathcal{L}_{\rm per}^{(m)}(\theta) + \lambda_{\rm reg}\mathcal{L}_{\rm reg}^{(m)}(\theta).\tag{53}\]

6 Numerical results↩︎

6.1 Benchmark 1↩︎

We first test the method on the free one-dimensional Schrödinger equation \[\label{eq:benchmark-1-pde} \mathrm i\partial_t\psi=-\frac{1}{2}\partial_{xx}\psi, \qquad x\in\mathbb{R},\tag{54}\] with Gaussian wave-packet initial data \[\label{eq:benchmark-1-ic} \psi_0(x) = \exp\!\left(-\frac{(x-x_0)^2}{2\sigma^2}\right) \exp(\mathrm ik_0 x), \qquad x_0=-3,\quad k_0=3,\quad \sigma=1.\tag{55}\] The physical interval is \(\Omega_L=(-L,L)\), with \(L=4\), and the TDPSF–PINN computation is carried out on the extended interval \(\Omega_{L,w}=(-L-w,L+w)\), with \(w=5\). The final time is \(T=4\). The exact solution of 54 –55 is \[\label{eq:benchmark-1-exact} \psi(t,x) = \frac{\sigma}{\sqrt{\sigma^2+\mathrm it}} \exp\!\left( -\frac{(x-x_0-k_0 t)^2}{2(\sigma^2+\mathrm it)} \right) \exp\!\left(\mathrm ik_0 x-\frac{\mathrm i}{2}k_0^2 t\right),\tag{56}\] which provides a pointwise reference on \([0,T]\times\mathbb{R}\).

6.1.0.1 Methods compared.

We compare four boundary treatments.

1. Dirichlet PINN. A PINN is trained on \(\Omega_L\times[0,T]\) with a soft homogeneous Dirichlet penalty \[|\psi(t,-L)|^2+|\psi(t,L)|^2\] on the artificial boundary.

2. Periodic PINN. A PINN is trained on \(\Omega_L\times[0,T]\) with a periodic penalty matching both \(\psi\) and \(\partial_x\psi\) at \(x=-L\) and \(x=L\).

3. Absorbing-layer PINN. A PINN is trained on \(\Omega_L\times[0,T]\) using the modified residual corresponding to \[\mathrm i\partial_t\psi = -\frac{1}{2}\partial_{xx}\psi-\mathrm iV_{\rm abs}(x)\psi,\] where \(V_{\rm abs}\) is a smooth nonnegative bump of height \(4\) supported in the layer \(|x|\in(L-1,L).\)

4. TDPSF–PINN. The slabwise method of Algorithm 3 is applied on \(\Omega_{L,w}\times[0,T]\), using five slabs of length \(\Delta t=0.8\). At the endpoint of each slab we apply the FFT-window outgoing filter 17 , with glancing buffer \(\gamma=0.2\), frequency-mask smoothing \(\alpha=0.1\), and Gaussian side windows of width \(\sigma_{\rm win}=0.6\). The slab length satisfies \(\Delta t\leq w/(3v_{\max})\) for the resolved frequency content of the wave packet.

6.1.0.2 Common neural ansatz.

To isolate the effect of the boundary treatment, all four methods use the same neural architecture: a sinusoidal representation network with four hidden layers of width \(48\) and first-layer frequency \(\omega_0=4\). Because the exact solution contains the carrier oscillation \(\exp(\mathrm ik_0x-\frac{\mathrm i}{2}k_0^2t)\), we train the network on the slowly varying envelope \(\phi(t,x)=\psi(t,x)\exp(-\mathrm i(k_0x-\frac{1}{2}k_0^2t))\), and then recover \(\psi\) analytically. The envelope satisfies \[\label{eq:bench1-envelope-eq} \mathrm i\partial_t\phi = -\frac{1}{2}\partial_{xx}\phi -\mathrm ik_0\partial_x\phi.\tag{57}\] The real and imaginary parts of the residual are computed by automatic differentiation. Each network is trained using Adam with a cosine learning-rate schedule, followed by a local refinement stage. For the TDPSF–PINN, each slab is warm-started from the parameters obtained on the previous slab. Thus the four methods differ only in the boundary mechanism and in the slabwise filtering structure.

6.1.0.3 Diagnostics.

We report the interior \(L^2\) error on the observation interval \(\Omega_{\rm obs}=(-3,3)\): \[\label{eq:bench1-Eint} E_{\rm int}(t) = \|\psi_{\rm num}(t,\cdot)-\psi(t,\cdot)\|_{L^2(\Omega_{\rm obs})}.\tag{58}\] Here \(\psi_{\rm num}\) denotes the numerical approximation produced by one of the four methods.

We also report the in-box norm \[\label{eq:bench1-Ebox} E_{\rm box}(t) = \|\psi_{\rm num}(t,\cdot)\|_{L^2(\Omega_L)}.\tag{59}\] This diagnostic is not itself a reflection measure. Rather, it measures whether the method retains the correct amount of mass inside the physical interval. The corresponding exact value is \(E_{\rm box}^{\rm exact}(t) = \|\psi(t,\cdot)\|_{L^2(\Omega_L)}.\)

The pointwise phase error is measured at the probe point \(x_\ast=0\): \[\label{eq:bench1-Ephase} E_{\rm phase}(t) = \left| \arg\!\left( \frac{\psi_{\rm num}(t,x_\ast)}{\psi(t,x_\ast)} \right) \right|,\tag{60}\] and is reported only at times for which \(|\psi(t,x_\ast)|>0.05.\)

For Method D, we additionally monitor the mass-balance residual \[\label{eq:bench1-massbal} \Delta M(t) = M_{\rm box}(t)+M_{\rm removed,cum}(t)-M_{\rm box}(0),\tag{61}\] where \(M_{\rm box}(t) = \|\psi_{\rm num}(t,\cdot)\|_{L^2(\Omega_L)}^2.\) The term \(M_{\rm removed,cum}(t)\) is the cumulative squared \(\ell^2\)-norm of the outgoing components removed by the filter up to time \(t\). This quantity is not an exact conservation law, because the filter is applied on the extended box and the outgoing components are not perfectly orthogonal. It should therefore be interpreted as a diagnostic of filter consistency rather than an exact invariant.

6.1.0.4 Final-time summary.

Table 1 gives the diagnostics at \(t=T\). The exact value \(E_{\rm box}^{\rm exact}(T)\approx 0.27\) is nonzero because the dispersing Gaussian retains a small but non-negligible tail inside \(\Omega_L\) at the final time. A faithful open-boundary method should reproduce this value rather than drive the solution to zero.

6.1.0.5 Snapshots.

Figure 9 overlays \(|\psi(t,x)|^2\) for the four methods at \(t=0,1,2,3,4.\)

At \(t=1\), before the packet center reaches the right boundary, Methods A and C already underestimate the peak amplitude. This indicates that the artificial boundary treatment affects the global least-squares fit, not only the solution near the boundary. Method B preserves the peak amplitude more effectively at early times, but the periodic constraint produces a spurious component entering from the left once the right-going packet has exited. Method D remains close to the exact solution at all displayed times.

6.1.0.6 Space–time visualization.

Figure 10 shows \(|\psi(t,x)|^2\) in the \((t,x)\)-plane for the exact reference and the four numerical solutions. The exact solution appears as a right-moving dispersive packet. Method B exhibits the expected periodic wraparound: as the packet exits near \(x=L\), a spurious component enters near \(x=-L\). Methods A and C damp the packet prematurely as it approaches the artificial boundary. Method D reproduces the reference space–time profile most accurately.

Figure 11 shows the corresponding pointwise error \(|\psi_{\rm num}(t,x)-\psi(t,x)|^2\). For Method D the error remains uniformly small on the plotted scale. For Methods A and C the error is largest along the packet trajectory, confirming that the boundary penalty contaminates the interior approximation. For Method B the dominant error is the wraparound component.

6.1.0.7 Diagnostics over time.

Figure 12 shows the four diagnostics as functions of time. Panel (a) displays the interior error \(E_{\rm int}(t)\). Method D remains substantially more accurate than the three boundary-penalty methods throughout the time interval. Panel (b) reports the phase error at \(x_\ast=0\); Method D maintains phase coherence for a longer time, while Methods A–C eventually lose phase accuracy. Panel (c) compares \(E_{\rm box}(t)\) with the exact in-box norm \(E_{\rm box}^{\rm exact}(t)\). Method D tracks the exact curve within a few percent, whereas Methods A and C underpredict the in-box norm and Method B overpredicts it because of periodic wraparound.

6.1.0.8 Mass-balance diagnostic.

Panel (d) of Figure 12 reports the quantities in 61 for Method D. Up to the slab boundary near \(t=1.6\), the packet is essentially contained in \(\Omega_L\), so \(M_{\rm box}\) remains close to its initial value and \(M_{\rm removed,cum}\) is small. Between \(t=1.6\) and \(t=3.2\), the packet passes through the buffer. Over this interval, \(M_{\rm removed,cum}\) increases from approximately \(0.01\) to \(1.08\), while \(M_{\rm box}\) decreases from approximately \(1.66\) to \(0.24\). The quantity \[M_{\rm box}+M_{\rm removed,cum}\] tracks \(M_{\rm box}(0)\approx 1.63\) to within about \(20\%\).

The remaining discrepancy has two sources. First, the PINN endpoint approximation in the buffer is not exact. Second, a small portion of the spreading packet can move beyond the effective support of the window during a slab. Thus the mass-balance residual provides a useful diagnostic of the accuracy budget of the filtering step. This is consistent with the fail-gracefully philosophy of the method: the algorithm does not merely remove outgoing waves, but also reports the extent to which the filtering process is consistent with the expected open-domain mass balance.

6.1.0.9 Concluding observation.

Benchmark 1 shows that replacing local artificial-boundary penalties by a phase-space outgoing-wave deletion mechanism substantially improves the neural approximation of the open-domain Schrödinger flow. The Dirichlet and absorbing-layer PINNs overdamp the wave packet, while the periodic PINN introduces wraparound. In contrast, the TDPSF–PINN produces the smallest interior error and preserves the correct in-box norm. This indicates that, for this benchmark, the dominant error in the boundary-penalty PINNs is not the expressive capacity of the neural network, but the incompatibility between the artificial boundary condition and the outgoing dispersive dynamics.

None

None

None

None

6.2 Benchmark 2↩︎

The second benchmark adds a localized potential. In this setting, the artificial boundary must not corrupt the physical reflected and transmitted waves generated by the potential, while outgoing components should leave the computational box without producing spurious echoes. We consider \[\label{eq:bench2-pde} \mathrm i\partial_t\psi = -\frac{1}{2}\partial_{xx}\psi+V(x)\psi, \qquad x\in\mathbb{R},\tag{62}\] with Gaussian potential \[V(x)=V_0e^{-x^2/a^2}, \qquad V_0=4, \qquad a=1,\] and incoming Gaussian wave packet \[\label{eq:bench2-ic} \psi_0(x) = \exp\!\left(-\frac{(x-x_0)^2}{2\sigma^2}\right)e^{\mathrm ik_0x}, \qquad x_0=-5,\quad k_0=3,\quad \sigma=1.\tag{63}\] The physical interval is \(\Omega_L=(-L,L)\), with \(L=6\), and the TDPSF–PINN computation is performed on the extended interval \(\Omega_{L,w}=(-L-w,L+w)\), with \(w=5\). The final time is \(T=5\). Since the kinetic energy \(k_0^2/2=4.5\) is comparable to the barrier height \(V_0=4\), the packet undergoes strong above-barrier scattering. A high-resolution reference calculation on \((-30,30)\), described below, shows that by \(t=T\) the solution contains a transmitted component, a reflected component, and a slowly decaying remnant near the potential. Thus this benchmark tests whether the boundary treatment removes only waves that are outgoing at the artificial boundary, while preserving the physical scattering dynamics in the interior.

6.2.0.1 Reference solution.

No closed-form solution is available. We compute the reference solution by a second-order Strang split-step Fourier method on the large interval \((-30,30)\), using \(N_x=8192\) grid points and time step \(\Delta t_{\rm ref}=2\times 10^{-3}\). The reference domain is large enough that no wave reaches its boundary before \(t=T\), and the scheme conserves \(\|\psi(t,\cdot)\|_{L^2(\mathbb{R})}^2\) to five significant figures over the time interval considered. We therefore use the restriction of this large-domain solution to the observation interval \(\Omega_{\rm obs}=(-4,4)\) as the reference solution.

6.2.0.2 Methods compared.

We compare the same four boundary treatments as in Section 6.1.

1. Dirichlet PINN. A PINN is trained on \(\Omega_L\times[0,T]\) with a soft homogeneous Dirichlet penalty at \(x=\pm L\).

2. Periodic PINN. A PINN is trained on \(\Omega_L\times[0,T]\) with a periodic penalty matching both \(\psi\) and \(\partial_x\psi\) at \(x=-L\) and \(x=L\).

3. Absorbing-layer PINN. A PINN is trained on \(\Omega_L\times[0,T]\) using the residual of the complex-absorbing-potential equation \[\mathrm i\partial_t\psi = -\frac{1}{2}\partial_{xx}\psi+\bigl(V(x)-\mathrm iW(x)\bigr)\psi,\] where \(W\) is a smooth nonnegative absorber supported in \(|x|\in(L-1.2,L).\)

4. TDPSF–PINN. The slabwise method of Algorithm 3 is applied on \(\Omega_{L,w}\times[0,T]\), using eight slabs of length \(\Delta t=0.625\). At the end of each slab, we apply the FFT-window outgoing filter 17 .

All four methods use the same SIREN architecture, with four hidden layers of width \(56\) and first-layer frequency \(\omega_0=10\). The residual in each case includes the potential term \(V\psi\). Since the scattered solution contains both right-moving and left-moving components, we train the network directly on \(\psi\), rather than on a single carrier-envelope representation.

6.2.0.3 Stabilizing the optimization.

For this strongly scattering problem, the bare PINN objective can admit near-trivial approximations in which the wave amplitude decays spuriously in the interior. To improve conditioning, we add a mild mass-profile anchor at a small number of monitoring times. The target is chosen consistently with each method’s own governing equation and boundary mechanism. For the self-adjoint box methods, namely the Dirichlet and periodic PINNs, the target is the conserved initial box mass. For the absorbing-layer PINN, the target is the box-mass profile of the corresponding Schrödinger equation on \(\Omega_L\). For the TDPSF–PINN, the target on each slab is the slab-initial mass before the endpoint filter is applied. This anchor is used only to stabilize the optimization. It does not modify the PDE residual or the boundary mechanism being tested. In particular, we do not train any method using the open-domain reference mass profile.

6.2.0.4 Diagnostics.

We report the interior error \[\label{eq:bench2-Eint} E_{\rm int}(t) = \|\psi_{\rm num}(t,\cdot)-\psi_{\rm ref}(t,\cdot)\|_{L^2(\Omega_{\rm obs})}, \qquad \Omega_{\rm obs}=(-4,4),\tag{64}\] where \(\psi_{\rm num}\) denotes the numerical approximation and \(\psi_{\rm ref}\) denotes the large-domain reference solution.

We also report the in-box norm \[\label{eq:bench2-Ebox} E_{\rm box}(t) = \|\psi_{\rm num}(t,\cdot)\|_{L^2(\Omega_L)}.\tag{65}\] This is not a direct reflection measure. Rather, it measures whether the method keeps the correct amount of wave mass inside the physical box. The corresponding reference quantity is \[E_{\rm box}^{\rm ref}(t) = \|\psi_{\rm ref}(t,\cdot)\|_{L^2(\Omega_L)}.\]

The probe phase error is evaluated at \(x_\ast=2\), on the transmitted side, at times for which the reference amplitude at the probe is not too small: \[\label{eq:bench2-phase-error} E_{\rm phase}(t) = \left| \arg\!\left( \frac{\psi_{\rm num}(t,x_\ast)}{\psi_{\rm ref}(t,x_\ast)} \right) \right|.\tag{66}\]

Finally, we introduce the late-time artificial echo diagnostic \[\label{eq:bench2-echo} E_{\rm echo}(t) = \|\psi_{\rm num}(t,\cdot)\|_{L^2(\Omega_{\rm obs})},\tag{67}\] read against the reference value \(E_{\rm echo}^{\rm ref}(t) = \|\psi_{\rm ref}(t,\cdot)\|_{L^2(\Omega_{\rm obs})}.\) After the physical transmitted and reflected packets have largely left \(\Omega_{\rm obs}\), the excess \(E_{\rm echo}(t)-E_{\rm echo}^{\rm ref}(t)\) measures spurious signal re-entering the observation region from the artificial boundary.

6.2.0.5 Final-time summary.

Table 2 reports the diagnostics at \(t=T\). The TDPSF–PINN gives the smallest interior error, an in-box norm closest to the reference, and the smallest artificial echo. The Dirichlet and periodic PINNs produce the largest late-time echoes, due respectively to wall reflection and wraparound. The absorbing-layer PINN is better than the Dirichlet and periodic methods, but still leaves a visible residual echo.

6.2.0.6 Snapshots.

Figure 13 shows \(|\psi|^2\) at \(t=0,\;1.5,\;2.5,\;3.5,\;5,\) with the support of the potential shaded. During the interaction (\(t=1.5\) and \(t=2.5\)), the TDPSF–PINN tracks the reference as the packet splits into a transmitted lobe moving to the right and a reflected lobe moving to the left. By \(t=5\), the reference retains only a small remnant near the potential. The TDPSF–PINN reproduces this remnant, whereas the Dirichlet PINN displays a spurious peak near the right boundary and the periodic PINN shows a component re-entering from the left.

6.2.0.7 Space–time and error fields.

Figure 14 plots \(|\psi(t,x)|^2\) in the \((t,x)\)-plane. The reference solution shows the incoming packet, the scattering event near \((t,x)\approx(1.7,0)\), and two outgoing rays. The TDPSF–PINN reproduces this pattern, with both outgoing components leaving through \(x=\pm L\). The Dirichlet PINN reflects outgoing waves off the artificial walls, while the periodic PINN wraps the transmitted component back into the computational interval. The absorbing-layer PINN damps the outgoing field, but still leaves a visible discrepancy in the interior.

The pointwise error fields in Figure 15 show the same behavior quantitatively. The TDPSF–PINN error is localized mainly near the scattering event, while the other methods develop larger errors along reflected or wrapped boundary-generated rays.

6.2.0.8 Diagnostics over time.

Figure 16 shows the diagnostics as functions of time. Panel (a) shows that the TDPSF–PINN has the smallest post-scattering interior error. Panel (b) shows that its phase error at the probe point remains below approximately \(0.6\) radians, while the local-boundary methods lose phase coherence. Panel (c) compares the in-box norm \(E_{\rm box}(t)\) with the reference \(E_{\rm box}^{\rm ref}(t)\). The TDPSF–PINN follows the physical decay most closely; the Dirichlet and periodic PINNs remain artificially large because of boundary-generated echoes. Panel (d) shows the echo diagnostic: in the late-time window, where the physical wave has mostly left \(\Omega_{\rm obs}\), the TDPSF–PINN stays closest to the reference.

6.2.0.9 Mass balance for Method D.

Table 3 reports the TDPSF mass budget. For early times, the packet is mostly contained in the physical interval and the cumulative removed mass is negligible. As the reflected and transmitted components reach the left and right buffers, the cumulative removed mass increases while the in-box mass decreases. This indicates that the filter removes the physically outgoing components as they reach the artificial boundary.

6.2.0.10 Conclusion.

Benchmark 2 confirms the main advantage of the phase-space filter over local boundary treatments. The filter removes waves that are outgoing at the artificial boundary while preserving the physically reflected field generated by the potential. The TDPSF–PINN achieves the smallest interior error, the closest in-box norm, and the smallest late-time artificial echo among the four methods. In contrast, the Dirichlet and periodic PINNs fail through wall reflection and wraparound, respectively, while the absorbing-layer PINN still reflects a non-negligible portion of the outgoing field back into the interior.

None

None

None

None

6.3 Benchmark 3↩︎

The third benchmark tests the central anisotropic feature of the phase-space filter: outgoing wave packets must be classified by their group velocity, not merely by their wave vector. We consider the two-dimensional anisotropic Schrödinger equation \[\label{eq:bench3-pde} \mathrm i\partial_t\psi = -\frac{1}{2}\nabla\cdot(A\nabla\psi), \qquad A= \begin{pmatrix} a_1 & 0\\ 0 & a_2 \end{pmatrix}, \qquad a_1\neq a_2 .\tag{68}\] For this constant-coefficient operator, the dispersion relation and group velocity are \[\label{eq:bench3-disp} \omega(k)=\frac{1}{2} k^T A k, \qquad v_g(k)=\nabla_k\omega(k)=Ak .\tag{69}\] Thus, when \(a_1\neq a_2\), the physical propagation direction of a wave packet is generally not parallel to its wave vector \(k\). The microlocal outgoing condition at a boundary face with outward normal \(n\) is therefore \(v_g(k)\cdot n=(Ak)\cdot n>0\), not simply \(k\cdot n>0\). We take \(a_1=1\) and \(a_2=6\), and use the physical box \(\Omega_L=(-R,R)^2\) with \(R=6\), buffer width \(w=4\), observation box \(\Omega_{\rm obs}=(-4,4)^2\), and final time \(T=5\). The initial condition is the Gaussian wave packet \[\label{eq:bench3-ic} \psi_0(x) = \exp\!\left(-\frac{|x-x_0|^2}{2\sigma^2}\right) e^{\mathrm ik_0\cdot x}, \qquad x_0=(0,-2),\quad k_0=(0.25,0.5),\quad \sigma=1.5 .\tag{70}\] For these parameters, \(v_g(k_0)=Ak_0=(0.25,3.0).\) Thus the packet moves slowly in the \(x_1\)-direction and rapidly in the \(x_2\)-direction. Its center therefore travels almost vertically upward and reaches the top side \(x_2=R\) at approximately \(t_{\rm exit} = \frac{R-x_{0,2}}{(v_g(k_0))_2} = \frac{6-(-2)}{3} \approx 2.67.\)

6.3.0.1 Where the two filters differ.

For a diagonal positive definite matrix \(A\) on an axis-aligned square, the group-velocity and wave-vector tests have the same sign on each straight face. For example, on the top face \(n=(0,1)\), \[v_g(k)\cdot n=a_2k_2, \qquad k\cdot n=k_2,\] and these quantities have the same sign because \(a_2>0\). Therefore, in this benchmark, the difference between the correct group-velocity filter and the naive wave-vector filter is not the sign of the outgoing test, but the treatment of the glancing buffer.

The filter should retain modes that leave the domain too slowly, since such near-glancing components can remain dynamically relevant in the interior. The correct measure of the exit speed is \(v_g(k)\cdot n.\) A naive isotropic implementation instead measures this speed by \(k\cdot n\). On the top face, the correct group-velocity filter removes modes with \(a_2k_2>\gamma\), equivalently \(k_2>\gamma/a_2\), whereas the naive filter removes only modes with \(k_2>\gamma\). With \(\gamma=2\) and \(a_2=6\), the correct threshold is \(k_2>\gamma/a_2=1/3\), while the naive threshold is \(k_2>2\). The packet is centered at \(k_{0,2}=0.5\). It is therefore genuinely fast-outgoing in the physical sense, since \((v_g(k_0))_2=a_2k_{0,2}=3.0>\gamma\), but it is incorrectly classified as glancing by the naive rule, since \(k_{0,2}=0.5<\gamma\). Figure 18 illustrates this distinction in the \((k_1,k_2)\)-plane: the spectral mass of the packet lies in the group-velocity removal region, but below the naive wave-vector threshold.

6.3.0.2 Interior solver and reference.

This benchmark is designed to isolate the boundary-filter classification. The distinction being tested is a property of the phase-space filter and does not depend on the particular interior time-stepping method. Since the PINN interior solver was already tested in Benchmarks 1 and 2, here we use an exact spectral propagator for the interior evolution. This prevents errors from the interior neural approximation from obscuring the comparison between the two filters.

Specifically, on the extended grid \((-10,10)^2\) with \(256^2\) points, we advance the solution by the exact Fourier multiplier \(e^{-\mathrm i\omega(k)\Delta t}\), where \(\omega(k)=\frac{1}{2} k^T A k\), and apply the windowed FFT outgoing filter once per slab. We use \(\Delta t=0.625\) and eight slabs. The reference solution is available in closed form: since \(A\) is diagonal, the equation separates into two one-dimensional anisotropic Schrödinger evolutions, the packet center moves with velocity \(v_g(k_0)=Ak_0\), and the \(j\)-th coordinate spreads according to the complex width \(\sigma^2+\mathrm ia_jt\). This exact solution is used as ground truth on \(\Omega_{\rm obs}\).

6.3.0.3 Filters compared.

Both filters use the same buffer, spatial windows, filtering frequency, and threshold \(\gamma=2\). They differ only in the outgoing classifier: \[\label{eq:bench3-filters} \text{naive:}\quad k\cdot n(q)>\gamma, \qquad \text{group-velocity:}\quad v_g(k)\cdot n(q)=(Ak)\cdot n(q)>\gamma .\tag{71}\] On the left and right faces, \(a_1=1\), so the two filters coincide. On the top and bottom faces, \(a_2=6\), and the two filters differ by the factor \(a_2\) in the normal exit speed.

6.3.0.4 Results.

Table 4 reports the diagnostics. The group-velocity filter removes the outgoing packet as it reaches the top buffer and remains close to the exact open-domain reference. The naive filter removes almost no mass, because it classifies the dominant packet as glancing rather than outgoing. Consequently, spurious reflected mass returns to the physical box.

Figure 17 shows \(|\psi|^2\) at five times. In the reference solution, the anisotropic packet moves upward and exits through the top face. The group-velocity filter reproduces this behavior and leaves only the physically correct residual tail inside \(\Omega_L\). The naive filter fails to remove the packet at the top buffer; the unremoved outgoing wave reflects and a bright lobe re-enters the box, clearly visible at late times.

Figure 19 shows the time histories. The interior error of the naive filter grows to about \(1.2\), whereas the group-velocity filter remains substantially smaller. The phase error also separates the two filters: once the reflected wave returns, the naive phase error spikes above \(2\) radians, while the group-velocity filter stays below approximately \(0.7\) radians. The in-box norm shows the same behavior: the group-velocity filter follows the physical decay of the reference, whereas the naive curve plateaus as reflected mass accumulates. The mass-balance panel gives the clearest diagnostic: the naive filter removes almost no mass, while the group-velocity filter removes the outgoing packet as it reaches the buffer.

6.3.0.5 Conclusion.

Benchmark 3 verifies the microlocal principle underlying the anisotropic phase-space filter. In anisotropic media, outgoing behavior must be classified by the group velocity \(v_g(k)\), not by the wave vector \(k\). In this test, the naive wave-vector rule treats a fast-outgoing packet as glancing and fails to remove it. The group-velocity TDPSF filter correctly identifies the outgoing packet, suppresses artificial reflection, and produces a substantially more accurate interior solution.

None

None

None

6.4 Benchmark 4↩︎

The fourth benchmark considers the cubic nonlinear Schrödinger equation \[\label{eq:bench4-pde} \mathrm i\partial_t\psi = -\frac{1}{2}\partial_{xx}\psi+\beta|\psi|^2\psi, \qquad x\in\mathbb{R},\tag{72}\] in both the defocusing case \(\beta=+1\) and the focusing case \(\beta=-1\). The initial condition is a traveling Gaussian packet, \[\label{eq:bench4-ic} \psi_0(x) = A\exp\!\left(-\frac{(x-x_0)^2}{2\sigma^2}\right)e^{\mathrm ik_0x}, \qquad A=1.2,\quad x_0=-3,\quad k_0=3,\quad \sigma=1 .\tag{73}\] The physical interval is \(\Omega_L=(-4,4)\), the buffer width is \(w=5\), the observation region is \(\Omega_{\rm obs}=(-2.5,2.5)\), and the final time is \(T=3.5\).

This benchmark is the nonlinear analogue of Benchmark 1. Its purpose is to test whether the time-dependent phase-space filter can be coupled effectively with a nonlinear PINN residual solver. The amplitude \(A=1.2\) is large enough for the cubic term to influence the dynamics: the defocusing solution spreads mildly, while the focusing solution forms a sharper coherent packet before leaving the physical interval.

6.4.0.1 Reference and methods.

The reference solution is computed using a Strang split-step Fourier method for 72 on the large interval \((-40,40)\), with \(N_x=8192, \Delta t_{\rm ref}=10^{-3}.\) The domain is large enough that no wave reaches its boundary before \(T\). We compare three PINN-based methods.

1. Periodic PINN. A PINN is trained on \(\Omega_L\times[0,T]\) with periodic matching of \(\psi\) and \(\partial_x\psi\) at \(x=-4\) and \(x=4\).

2. Absorbing-layer PINN. A PINN is trained on \(\Omega_L\times[0,T]\) with a complex absorbing potential. The residual corresponds to \(\mathrm i\partial_t\psi=-\frac{1}{2}\partial_{xx}\psi+\beta|\psi|^2\psi-\mathrm iW(x)\psi\), where \(W\) is a quadratic absorber of strength \(W_0=3\), supported in the outer layer \(|x|\in(L-1.5,L)\).

3. TDPSF–PINN. The TDPSF–PINN is advanced in five slabs of length \(\Delta t=0.7.\) At the end of each slab, the FFT-window outgoing filter is applied. Each slab is trained with the temporal curriculum described in Section 6.2.

All three methods use the same SIREN architecture and the same carrier-envelope representation \(\psi(t,x)=\phi(t,x)e^{\mathrm i(k_0x-\frac{1}{2} k_0^2t)}\). Since \(|\psi|=|\phi|\), the cubic term in the envelope equation remains \(\beta|\phi|^2\phi\). As in the previous benchmarks, each PINN uses a method-consistent mass anchor. For the conservative periodic method the target is the constant initial mass. For the absorbing-layer method the target is the mass profile associated with the corresponding absorbing-layer NLS problem. For the TDPSF–PINN, the target on each slab is the slab initial mass before endpoint filtering. The open-domain reference is not used as a training target.

6.4.0.2 Diagnostics.

We report the interior error \(E_{\rm int}(t) = \|\psi_{\rm num}(t,\cdot)-\psi_{\rm ref}(t,\cdot)\|_{L^2(\Omega_{\rm obs})},\) the in-box norm \[\label{eq:bench4-Ebox} E_{\rm box}(t) = \|\psi_{\rm num}(t,\cdot)\|_{L^2(\Omega_L)},\tag{74}\] and the probe phase error at \(x_\ast=2\), whenever the reference amplitude at the probe is not too small. For the TDPSF–PINN, we also monitor the mass-balance residual \(M_{\rm box}(t)+M_{\rm removed,cum}(t)-M_{\rm box}(0)\), where \(M_{\rm box}(t)=\|\psi_{\rm num}(t,\cdot)\|_{L^2(\Omega_L)}^2\).

In addition, we monitor the nonlinear buffer quantity \[\label{eq:bench4-nbuf} N_{\rm buf}(t) = \int_{\mathcal{B}_w}|\psi_\theta(t,x)|^4\,\mathrm dx,\tag{75}\] where \(\mathcal{B}_w\) is the buffer region. This quantity measures the amount of nonlinear self-interaction present in the filtering layer. A large value of \(N_{\rm buf}\) does not automatically indicate failure: it may occur when a genuine nonlinear packet is physically crossing the buffer. However, if a local boundary treatment traps or wraps the packet, then \(N_{\rm buf}\) remains large after the outgoing component should have left. In contrast, the TDPSF mechanism should reduce \(N_{\rm buf}\) after the outgoing packet is filtered.

6.4.0.3 Results.

Table 5 reports the final-time diagnostics. In both the defocusing and focusing cases, the TDPSF–PINN gives the smallest interior error. The periodic method retains too much mass because the packet wraps around the computational interval. The absorbing-layer method reduces the in-box norm, but its interior error remains substantially larger than that of the TDPSF–PINN.

Figures 20 and 21 show \(|\psi|^2\) on the physical interval. In the defocusing case, the periodic PINN develops oscillatory structure re-entering from the left, whereas the PINN-TDPSF follows the reference spreading and rightward drift. In the focusing case, the reference packet self-focuses into a sharp coherent peak and exits the physical interval. The PINN-TDPSF reproduces this coherent peak and the subsequent decay of the in-box solution. The periodic PINN leaves a wrapped lobe, and the absorbing PINN leaves a distorted residual.

Figure 22 shows the time histories. Panel (a) reports the interior error. In the defocusing case, the TDPSF curve remains small throughout the run. In the focusing case, all methods experience increased error during the self-focusing event near \(t\approx1.6\), which is an interior effect common to all methods. After the packet reaches the boundary, the errors separate: the periodic error grows due to wraparound, the absorbing error plateaus, and the TDPSF error decreases as the outgoing packet is removed. Panel (b) shows that the TDPSF in-box norm tracks the reference decay, while the periodic norm remains artificially high. Panel (c) shows the nonlinear buffer quantity \(N_{\rm buf}\). The local treatments sustain large buffer activity because trapped or wrapped packets continue to self-interact near the boundary; the TDPSF curve follows the reference during physical exit and then drops once the outgoing packet is filtered.

6.4.0.4 A note on the comparison.

The PINN-TDPSF is advanced slabwise, whereas the periodic and absorbing PINNs are trained over the full time interval. Therefore, part of the early-time interior-error gap may reflect the better conditioning of a slabwise representation, not only the boundary mechanism. The boundary-specific evidence is nevertheless clear: the periodic wraparound, the post-interaction divergence in the focusing case, and the sustained nonlinear buffer activity of the local methods are all effects removed by the phase-space filter.

6.4.0.5 Conclusion.

Benchmark 4 shows that the time-dependent phase-space filter can be combined effectively with a nonlinear PINN residual solver. For both signs of the cubic nonlinearity, the PINN-TDPSF removes the outgoing wave packet, tracks the reference interior solution more accurately than the local boundary treatments, and suppresses the wrapped or reflected packets retained by the periodic and absorbing PINNs. The nonlinear buffer diagnostic further confirms that the phase-space filter prevents persistent nonlinear activity from accumulating in the boundary layer.

None

None

None

6.5 Benchmark 5↩︎

The fifth benchmark is a stress test for the nonlinear filtering mechanism. We consider the focusing cubic nonlinear Schrödinger equation \[\label{eq:bench5-pde} \mathrm i\partial_t\psi = -\frac{1}{2}\partial_{xx}\psi-|\psi|^2\psi, \qquad x\in\mathbb{R},\tag{76}\] with traveling soliton initial data (\(\eta\) is the focusing-soliton amplitude parameter) \[\label{eq:bench5-ic} \psi_0(x) = \eta\,\mathop{\mathrm{sech}}\!\big(\eta(x-x_0)\big)e^{\mathrm iv x}, \qquad x_0=-3, \qquad v=2.5.\tag{77}\] The physical interval is \(\Omega_L=(-4,4)\), the buffer width is \(w=5\), the observation region is \(\Omega_{\rm obs}=(-2.5,2.5)\), and the final time is \(T=3.5\). The exact traveling soliton is \[\label{eq:bench5-exact-soliton} \psi(t,x) = \eta\,\mathop{\mathrm{sech}}\!\big(\eta(x-x_0-vt)\big) \exp\!\left( \mathrm i\left( vx-\frac{1}{2}v^2t+\frac{1}{2}\eta^2t \right) \right).\tag{78}\] This exact solution is used as the reference.

A soliton is a coherent nonlinear object. The phase-space filter, however, is based on the far-field assumption that the boundary field is well described by a superposition of outgoing linear wave packets. If a strong soliton reaches the buffer while it is still highly nonlinear, this assumption is violated. The purpose of the present benchmark is therefore not only to test absorption, but also to test whether the method detects when its filtering assumptions are no longer reliable.

6.5.0.1 Reliability flag.

The soliton amplitude \(\eta\) controls the strength of the nonlinear interaction. For a sech soliton, \(\int_{\mathbb{R}}|\psi(t,x)|^4\,\mathrm dx = \frac{4}{3}\eta^3.\) Thus the nonlinear buffer diagnostic \[\label{eq:bench5-nbuf} N_{\rm buf}(t) = \int_{\mathcal{B}_w}|\psi_\theta(t,x)|^4\,\mathrm dx\tag{79}\] scales like \(\eta^3\) when the soliton is inside the buffer. We use the following reliability rule: \[N_{\rm buf}(t_m)>N_{\rm crit} \quad\Longrightarrow\quad \text{the filtering step at }t_m\text{ is flagged as unreliable}.\] In the experiments below we take \(N_{\rm crit}=1.\) This threshold is not used to tune the solution; it is a diagnostic for detecting whether the boundary field is too nonlinear for the linear outgoing-packet classifier.

We run two cases with the same geometry and velocity, \(\eta=0.6\) and \(\eta=2.2\). The first is a broad, weak soliton that is closer to a dispersive wave packet, while the second is a narrow, strongly bound soliton. The solver is the PINN-TDPSF used in the previous benchmarks, with five slabs of length \(\Delta t=0.7\), a carrier-envelope SIREN architecture, and causal slab training. The diagnostic \(N_{\rm buf}\) is evaluated from the network field at each slab endpoint before filtering.

6.5.0.2 Results.

Table 6 summarizes the two outcomes. For the weak soliton, the nonlinear buffer diagnostic remains below the threshold throughout the computation. The filter removes the outgoing packet with small error and produces an in-box norm close to the reference. In contrast, for the strong soliton, \(N_{\rm buf}\) exceeds the threshold as the soliton enters the buffer. The run is therefore flagged as unreliable, and the late-time solution confirms the warning: a substantial spurious residual remains in the box.

The weak soliton, shown in the top row of Figure 23, is broad and close to the radiative regime for which the filter is designed. The PINN-TDPSF follows the reference as the packet moves into the buffer, and the physical box empties in agreement with the exact solution.

The strong soliton, shown in the bottom row, is a sharply localized coherent structure. The network already has a larger interior error before the boundary interaction, and when the soliton reaches the buffer, the linear outgoing-packet assumption is no longer appropriate. The filter does not remove the structure cleanly and a reflected or residual component remains in the physical box. Crucially, this failure is detected by \(N_{\rm buf}\) without using the exact reference solution. Panel (c) of Figure 24 shows that the weak soliton remains below \(N_{\rm crit}\) at every filtering step, whereas the strong soliton crosses into the flagged region as it enters the buffer.

6.5.0.3 A note on the failure mode.

The strong case combines two difficulties. First, the PINN approximation itself has to represent a narrow coherent peak, which is a more demanding interior approximation problem. Second, the phase-space filter assumes that the boundary field is approximately linear and radiative, an assumption that is violated by a strongly bound soliton. Both effects reflect the same physical issue: the object reaching the buffer is too nonlinear and too coherent for a far-field linear outgoing-packet treatment. The nonlinear buffer diagnostic \(N_{\rm buf}\) is therefore a practical, reference-free indicator that the open-boundary filter should not be certified in this regime.

6.5.0.4 Conclusion.

Benchmark 5 demonstrates the intended two-sided behavior of the method. When the soliton is weak enough that it reaches the buffer in an approximately radiative state, the PINN-TDPSF removes it with small reflection and low interior error. When a strongly nonlinear soliton reaches the buffer, the diagnostic \(N_{\rm buf}\) exceeds its threshold and correctly warns that the outgoing classification is unreliable. Thus the method does not merely attempt to absorb every outgoing object; it also identifies when the assumptions behind the phase-space filter are violated.

None

None

6.6 Benchmark 6↩︎

The sixth benchmark tests the method on a genuine anisotropic, multi-branch wave system: the two-dimensional linearized Euler equations about a uniform mean flow. We consider the convected acoustic system \[\label{eq:bench6-pde} \partial_t U + A_1\partial_{x_1}U + A_2\partial_{x_2}U = 0, \qquad U=(p,u,v)^T,\tag{80}\] with unit background density and unit sound speed. The mean flow is \(M=(U_0,V_0).\) The far-field dispersion branches are \[\label{eq:bench6-branches} \omega_{\rm vort}(k)=M\cdot k, \qquad \omega_{\rm ac\pm}(k)=M\cdot k\pm |k|,\tag{81}\] with group velocities \[\label{eq:bench6-vg} v_{g,{\rm vort}}(k)=M, \qquad v_{g,{\rm ac}\pm}(k)=M\pm\frac{k}{|k|}, \qquad k\neq 0.\tag{82}\] The vorticity branch is the decisive case. It is convected with the mean flow \(M\), independently of the wave-vector direction. Thus a vortical mode whose wave vector points upstream can nevertheless be transported downstream. A boundary classifier based only on \(k\cdot n\) can therefore mislabel such modes. The correct criterion is the branchwise group-velocity condition \(v_{g,\ell}(k)\cdot n.\)

6.6.0.1 Incoming-gust test.

To stress the classifier, we use a configuration in which the artificial boundary must pass a physical incoming wave and remove it only later after it becomes outgoing. We take the mean flow to be \(M=(U_0,V_0)=(1.2,0)\). The initial data consist of an incoming vortical gust placed in the upstream buffer, centered at \((-7,0)\), together with an acoustic pressure pulse centered at the origin. The vortical gust is initialized as a pure vorticity mode, \(p=0\), \(u=-\partial_y\Psi\), and \(v=\partial_x\Psi\), where \(\Psi\) is a localized Gaussian stream function. The gust is convected through the physical box and exits through the downstream face. The physical box is \(\Omega_L=(-5,5)^2\), the buffer width is \(w=5\), the observation box is \(\Omega_{\rm obs}=(-3,3)^2\), and the final time is \(T=11\). At the upstream left face the vortical gust satisfies \(v_{g,{\rm vort}}\cdot n<0\), so it is physically incoming and must be retained. At the downstream right face the same branch satisfies \(v_{g,{\rm vort}}\cdot n>0\), so it is outgoing and should be removed. This is precisely the distinction that a branchwise group-velocity filter is designed to make.

6.6.0.2 Setup and interior solver.

As in Benchmark 3, the purpose of this experiment is to isolate the boundary-filter classification. The PINN interior solver was tested in Benchmarks 1, 2, and 4; here, to compare boundary treatments cleanly, we advance the interior with the exact spectral propagator associated with the constant-coefficient far-field system. The computation uses an extended grid with \(320^2\) points and fourteen slabs. The reference solution is computed on the larger domain \((-24,24)^2\), so that no wave reaches the reference boundary before \(T\). At each slab endpoint, the buffer field is localized by side windows, Fourier transformed, decomposed into branch components using the spectral projectors \(\Pi_\ell(k)\), with \(\ell\in\{\mathrm{vort},\mathrm{ac}+,\mathrm{ac}-\}\), and then filtered according to the branchwise outgoing criterion. For a side with outward normal \(n\), the idealized branchwise filter is \[\label{eq:bench6-filter} \widehat U(k) \longmapsto \sum_{\ell} \bigl(1-\chi_{{\rm out},\ell}(k;n)\bigr) \Pi_\ell(k)\widehat U(k), \qquad \chi_{{\rm out},\ell}(k;n) = \mathbf{1}_{\{v_{g,\ell}(k)\cdot n>\gamma\}},\tag{83}\] with glancing threshold \(\gamma=0.25.\) In the actual implementation, this operation is applied in each side buffer using the windowed FFT formulation described earlier.

6.6.0.3 Methods compared.

We compare four boundary treatments. The same interior spectral propagator is used for all four methods; only the boundary treatment is changed. This isolates the effect of the outgoing-wave classifier.

Method A: untreated/periodic boundary. The solution is advanced on the extended computational box without any absorbing or filtering mechanism. Since the spectral propagator is implemented on a finite box, this corresponds to periodic continuation of the field. Thus outgoing waves are not removed; after reaching the outer part of the computational domain, they re-enter from the opposite side. Method A is included as a baseline for the wraparound error that occurs when no effective open-boundary treatment is applied.

Method B: damped sponge layer. The equation is modified in the buffer by adding a scalar damping term. In the buffer, the system is evolved with \(\partial_t U + A_1\partial_{x_1}U + A_2\partial_{x_2}U = -\sigma_{\rm damp}(x)U,\) where \(\sigma_{\rm damp}\geq0\) is supported in \(\mathcal{B}_w=\Omega_{L,w}\setminus\Omega_L.\) For example, on each side one may take \(\sigma_{\rm damp}(x) = \sigma_0 \left( \frac{\operatorname{dist}(x,\Omega_L)}{w} \right)^2 \qquad \text{in } \mathcal{B}_w,\) with \(\sigma_{\rm damp}=0\) in \(\Omega_L\). This method damps all components entering the buffer, regardless of whether they are incoming, outgoing, acoustic, or vortical. It is a useful local absorbing baseline, but it is not selective: in the incoming-gust test it also attenuates physically incoming vortical content that should be allowed to enter the physical box.

Method C: branchwise group-velocity TDPSF. This is the proposed phase-space filter.

Method D: naive non-branchwise \(k\cdot n\) filter. This method uses the same windows, FFT grid, filtering times, and thresholds as Method C, but replaces the branchwise group-velocity criterion by a scalar wave-vector criterion. For side \((j,s)\), it uses \(P_{j,\rm naive}^s(k) = S_\alpha(k\cdot n_j^s-\gamma)I,\) where \(I\) is the identity matrix on the state vector \(U=(p,u,v)^T\). Thus it removes the same Fourier frequencies from all components and all branches, without projecting onto acoustic and vortical modes. This is intentionally incorrect for anisotropic or multi-branch systems. It is included to demonstrate that filtering by \(k\cdot n\) alone is not sufficient: for the vorticity branch, the physical transport direction is \(M\), independent of \(k\). Consequently, the naive filter may delete physically incoming vortical modes at the inflow and retain outgoing modes at the outflow, depending on the sign of their wave vector rather than their group velocity. Method D removes the components satisfying \(k\cdot n>\gamma\) from the full field, without decomposing into dispersion branches. It is included to isolate the effect of the branchwise group-velocity classification.

6.6.0.4 Results.

Table 7 and Figure 26 report the diagnostics. During the time window in which the gust occupies \(\Omega_{\rm obs}\), approximately \(3.3\leq t\leq 8.3,\) the branchwise TDPSF keeps the interior error small and its in-box energy tracks the reference. The vortical gust is convected through the box without being removed at the inflow and is removed only after it reaches the outflow.

The naive filter behaves differently. Because it uses \(k\cdot n\) rather than the vorticity branch velocity \(v_{g,{\rm vort}}=M\), it removes part of the physically incoming gust at the inflow. In this experiment it has already removed \(1.19\) units of incoming energy by \(t=0.8\). Consequently the gust enters the physical box corrupted, and the interior error rises. The damped layer also corrupts the incoming gust, because it attenuates fields in the buffer without distinguishing incoming from outgoing components. The untreated boundary passes the incoming gust, but leaves outgoing content to wrap around and contaminate the interior at later times.

Figure 25 shows the vorticity field. In the reference solution, the gust is convected through the box and exits cleanly. The branchwise TDPSF reproduces this behavior. The damped layer reduces the gust to a weak smear, while the naive filter distorts it by removing part of its spectral content at the inflow. The glancing energy for the branchwise TDPSF remains small throughout the run, as shown in Figure 26 (d), indicating that the near-tangential components retained by the filter do not accumulate in this test.

6.6.0.5 Conclusion.

Benchmark 6 confirms that the branchwise group-velocity filter can suppress artificial reflected or wrapped waves without corrupting a physical incoming wave. The incoming-gust configuration makes this requirement explicit: the open boundary must pass incoming content at the inflow and remove the same branch only when it becomes outgoing at the outflow. The branchwise TDPSF does so. A damping layer attenuates the incoming wave, and a non-branchwise \(k\cdot n\) filter removes part of it at the inflow.

None

None

6.7 Sparse-data recovery of a localized acoustic defect↩︎

This benchmark tests the structural role of the PINN step in Microlocal PINNFilter. The unknown is not a flexible time-dependent source, but a localized medium parameter, namely the sound speed \(c(x)\). Different sound-speed profiles generate different propagation and scattering patterns. Therefore, a constant-coefficient FFT propagator, while highly accurate for a homogeneous medium, has no mechanism for representing the localized defect. This setting is designed to highlight the advantage of combining a residual-based neural interior model with the TDPSF open-boundary filter.

6.7.0.1 Setup.

We consider the one-dimensional convected acoustic system \[\begin{align} \partial_t p+M\partial_x p+c(x)^2\partial_x u &= F_p^\dagger(t,x), \tag{84}\\ \partial_t u+M\partial_x u+\partial_x p &= 0, \tag{85} \end{align}\] with known mean flow \(M=0.4.\) The source is known and localized in space and time: \(F_p^\dagger(t,x)=\chi_{\rm src}(t) \exp(-(t-t_0)^2/\tau_0^2)\exp(-(x-\xi_0)^2/a^2)\), where \(t_0=1.0\), \(\tau_0=0.5\), \(\xi_0=-1\), and \(a=0.4\). Here \(\chi_{\rm src}\) is a smooth cutoff supported in \([0,T_{\rm src}]\), with \(T_{\rm src}=2.5\). The unknown coefficient is the localized sound-speed defect \(c^\dagger(x)=c_\infty+A_c\exp(-(x-x_c)^2/a_c^2)\), with \(c_\infty=1\), \(A_c=0.5\), \(x_c=0\), and \(a_c=0.6\). Thus the wave speed is increased by \(50\%\) at the center of the defect. The defect is localized away from the buffer, so the far-field medium remains constant-coefficient.

The learned coefficient is parameterized as \(c_\eta(x)=c_\infty+\chi_{\rm int}(x)\widetilde{c}_\eta(x)\), where \(\chi_{\rm int}\) is supported in the physical interval \(\Omega_L=(-44)\) and vanishes in the buffer \(\mathcal{B}_w\). Hence \(c_\eta(x)=c_\infty\) in \(\mathcal{B}_w\), which ensures that the TDPSF filter is constructed from the constant far-field dispersion. We sample \(80\) noisy pressure-and-velocity sensors randomly over \([0,T_{\rm train}]\times\Omega_{\rm obs} =[0,3.5]\times[-2.5,2.5]\). The data have the form \(y_r=(p^\dagger,u^\dagger)(t_r,x_r)+\varepsilon_r\), for \(r=1,\ldots,80\). Since \(T_{\rm train}=3.5>T_{\rm src}=2.5\), the training data include post-source free propagation through the medium. The forecast window is \((T_{\rm train},T]=(3.5,5]\).

6.7.0.2 Methods.

We compare three methods. The comparison is designed to separate three effects: the error caused by assuming a homogeneous medium, the benefit of learning the localized coefficient from sparse data, and the best-case forecast when the true coefficient is known.

• Constant-\(c\) FFT–TDPSF. This method assumes that the medium is homogeneous: \(c(x)\equiv c_\infty.\) It therefore solves the constant-coefficient convected acoustic system \[\partial_t p+M\partial_x p+c_\infty^2\partial_x u = F_p^\dagger(t,x), \qquad \partial_t u+M\partial_x u+\partial_x p=0.\] The field is advanced using a Fourier/spectral propagator, and the open boundary is treated using the TDPSF filter built from the same constant far-field speed \(c_\infty\). Thus this method has an accurate classical forward solver and an accurate open-boundary mechanism for a homogeneous medium. However, it has no degree of freedom with which to represent the localized defect \(c^\dagger(x)-c_\infty.\) It is therefore structurally misspecified for the present problem. The purpose of this baseline is to show what happens if one applies the classical FFT–TDPSF approach while assuming that the medium is uniform.

• Microlocal PINNFilter. This is the proposed method. It approximates both the acoustic field and the unknown coefficient: \[U(t,x)=(p(t,x),u(t,x))\approx U_\theta(t,x) = (p_\theta(t,x),u_\theta(t,x)), \qquad c^\dagger(x)\approx c_\eta(x).\] The learned coefficient is constrained to be constant in the buffer: \[c_\eta(x) = c_\infty+\chi_{\rm int}(x)\widetilde{c}_\eta(x), \qquad c_\eta(x)=c_\infty \quad\text{for }x\in\mathcal{B}_w.\] The PDE residual used in the training is the variable-coefficient residual \[\mathcal{R}_{p,\theta,\eta} = \partial_t p_\theta + M\partial_x p_\theta + c_\eta(x)^2\partial_x u_\theta - F_p^\dagger(t,x),\] \[\mathcal{R}_{u,\theta} = \partial_t u_\theta + M\partial_x u_\theta + \partial_x p_\theta.\] The loss includes the PDE residual, sparse sensor misfit, initial matching across slabs, and a mild regularization on \(c_\eta\), such as \(\|\partial_x c_\eta\|_{L^2(\Omega_L)}^2 + \|\partial_{xx}c_\eta\|_{L^2(\Omega_L)}^2.\) At the end of each slab, the raw neural endpoint is evaluated on the FFT grid and filtered using the TDPSF open-boundary operator constructed from the far-field speed \(c_\infty\). Thus the PINN part learns the variable interior medium, while the TDPSF part supplies the open-boundary radiation mechanism.

• Oracle variable-coefficient solver. This method uses the true coefficient \(c^\dagger(x)\) in the forward solver: \[\partial_t p+M\partial_x p+c^\dagger(x)^2\partial_x u = F_p^\dagger(t,x), \qquad \partial_t u+M\partial_x u+\partial_x p=0.\] It is not an inverse method and does not learn from sparse data. It is included only as a numerical floor for the best achievable prediction when the coefficient is known exactly. The oracle therefore separates coefficient-recovery error from boundary or time-propagation error.

6.7.0.3 Results.

The constant-\(c\) FFT solver is structurally unable to represent the defect: its propagator has no coefficient degree of freedom that could encode the localized scatterer. Microlocal PINNFilter learns \(c_\eta(x)\) from the sparse sensor data and recovers a defect with the correct location and qualitative shape. The recovered peak is approximately \(1.30\), compared with the true peak \(1.50\); see Figure 27 (a).

The consequences for the wave field are shown in Figure 27 (b) and Table 8. Microlocal PINNFilter tracks the reference field to within approximately \(3\%\)–\(4\%\) during the training interval and yields a held-out forecast error \(E_{\rm pred}=0.31\), compared with \(E_{\rm pred}=0.61\) for the constant-\(c\) baseline; the oracle method, with \(c^\dagger\) supplied, gives \(E_{\rm pred}\approx 0.01\), which represents the numerical reference floor.

6.7.0.4 Interpretation.

This benchmark exhibits the structural advantage of Microlocal PINNFilter most clearly. When the unknown is a flexible source, a classical inverse method may exploit source equivalence: different sources can radiate similar fields at the available sensor locations. Here the unknown is encoded in the medium dispersion itself. The constant-coefficient FFT propagator therefore has no equivalent source-like degree of freedom with which to hide its model error. Microlocal PINNFilter combines the PINN’s ability to represent variable interior physics with the TDPSF filter’s phase-space open-boundary mechanism, and hence addresses a regime that the constant-coefficient FFT solver cannot represent.

None

None

7 Conclusion↩︎

We have introduced a Microlocal PINNFilter framework for wave propagation on unbounded domains. The method combines a slabwise neural residual approximation for the interior evolution with a time-dependent phase-space filter in the buffer region. Its guiding principle is microlocal: a component is removed only when it is localized near the artificial boundary and its group velocity points outward. Thus the boundary treatment is not a local penalty imposed on traces, but a phase-space radiation mechanism.

The framework was developed for linear Schrödinger equations, nonlinear Schrödinger equations, and anisotropic wave systems. The numerical experiments support the main claims of the method. The phase-space-filtered PINN reduces over-damping, wall reflection, and periodic wraparound in the free Schrödinger benchmark; preserves physical reflection and transmission in potential scattering; correctly uses group velocity in anisotropic Schrödinger propagation; and handles both defocusing and focusing nonlinear Schrödinger packets. The soliton stress test shows that the nonlinear-buffer diagnostic can flag regimes in which the filtering assumptions are no longer reliable. For linearized Euler waves, the branchwise filter preserves a physically incoming vortical gust at the inflow and removes it only after it becomes outgoing. Finally, the sparse-data acoustic defect benchmark illustrates the role of the PINN step in a variable-coefficient inverse setting: a constant-coefficient FFT–TDPSF solver cannot represent the localized sound-speed defect, while Microlocal PINNFilter learns it from sparse data and improves the forecast.

The method is not intended to replace classical FFT, spectral, or split-step solvers [13]–[16] for known-coefficient forward problems. When the PDE is linear, constant-coefficient, and exactly diagonalizable by a Fourier method, classical propagation combined with TDPSF remains the natural choice. The value of Microlocal PINNFilter lies in embedding the TDPSF open-boundary mechanism into a residual-based neural framework, making it applicable to nonlinear dynamics, sparse observations, variable or unknown coefficients, and other settings where the interior physics is not available as an exact spectral propagator.

Overall, the results indicate that replacing local artificial-boundary penalties by a microlocal phase-space filtering mechanism can improve neural wave computations on unbounded domains, while preserving the flexibility of PINNs for interior forward, inverse, and data-driven problems.

References↩︎

---

1. N. G. Hien and M.-B. T are funded in part by NSF CAREER DMS-2303146, and NSF Grants DMS-2204795, DMS-2305523, DMS-2306379. A.S. is funded in part by NSF Grant DMS-2205931.↩︎